LMDZ
vlspltqs_p.F
Go to the documentation of this file.
1 c
2 c $Id: vlspltqs_p.F 1907 2013-11-26 13:10:46Z lguez $
3 c
4  SUBROUTINE vlspltqs_p ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5  , p,pk,teta )
6 c
7 c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
8 c
9 c ********************************************************************
10 c Shema d'advection " pseudo amont " .
11 c + test sur humidite specifique: Q advecte< Qsat aval
12 c (F. Codron, 10/99)
13 c ********************************************************************
14 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
15 c
16 c pente_max facteur de limitation des pentes: 2 en general
17 c 0 pour un schema amont
18 c pbaru,pbarv,w flux de masse en u ,v ,w
19 c pdt pas de temps
20 c
21 c teta temperature potentielle, p pression aux interfaces,
22 c pk exner au milieu des couches necessaire pour calculer Qsat
23 c --------------------------------------------------------------------
24  USE parallel_lmdz
25  USE mod_hallo
26  USE vampir
27  IMPLICIT NONE
28 
29 c
30 #include "dimensions.h"
31 #include "paramet.h"
32 #include "logic.h"
33 #include "comvert.h"
34 #include "comconst.h"
35 
36 c
37 c Arguments:
38 c ----------
39  REAL masse(ip1jmp1,llm),pente_max
40  REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
41  REAL q(ip1jmp1,llm)
42  REAL w(ip1jmp1,llm),pdt
43  REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
44 c
45 c Local
46 c ---------
47 c
48  INTEGER i,ij,l,j,ii
49 c
50  REAL qsat(ip1jmp1,llm)
51  REAL zm(ip1jmp1,llm)
52  REAL mu(ip1jmp1,llm)
53  REAL mv(ip1jm,llm)
54  REAL mw(ip1jmp1,llm+1)
55  REAL zq(ip1jmp1,llm)
56  REAL temps1,temps2,temps3
57  REAL zzpbar, zzw
58  LOGICAL testcpu
59  SAVE testcpu
60  SAVE temps1,temps2,temps3
61 
62  REAL qmin,qmax
63  DATA qmin,qmax/0.,1.e33/
64  DATA testcpu/.false./
65  DATA temps1,temps2,temps3/0.,0.,0./
66 
67 c--pour rapport de melange saturant--
68 
69  REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
70  REAL ptarg,pdelarg,foeew,zdelta
71  REAL tempe(ip1jmp1)
72  INTEGER ijb,ije
73  type(request) :: MyRequest1
74  type(request) :: MyRequest2
75 
76 c fonction psat(T)
77 
78  foeew( ptarg,pdelarg ) = exp(
79  * (r3les*(1.-pdelarg)+r3ies*pdelarg) * (ptarg-rtt)
80  * / (ptarg-(r4les*(1.-pdelarg)+r4ies*pdelarg)) )
81 
82  r2es = 380.11733
83  r3les = 17.269
84  r3ies = 21.875
85  r4les = 35.86
86  r4ies = 7.66
87  retv = 0.6077667
88  rtt = 273.16
89 
90 c-- Calcul de Qsat en chaque point
91 c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
92 c pour eviter une exponentielle.
93 
94  call settag(myrequest1,100)
95  call settag(myrequest2,101)
96 
97  ijb=ij_begin-iip1
98  ije=ij_end+iip1
99  if (pole_nord) ijb=ij_begin
100  if (pole_sud) ije=ij_end
101 
102 
103  DO l = 1, llm
104  DO ij = ijb, ije
105  tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
106  ENDDO
107  DO ij = ijb, ije
108  zdelta = max( 0., sign(1., rtt - tempe(ij)) )
109  play = 0.5*(p(ij,l)+p(ij,l+1))
110  qsat(ij,l) = min(0.5, r2es* foeew(tempe(ij),zdelta) / play )
111  qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
112  ENDDO
113  ENDDO
114 
115 c PRINT*,'Debut vlsplt version debug sans vlyqs'
116 
117  zzpbar = 0.5 * pdt
118  zzw = pdt
119 
120  ijb=ij_begin
121  ije=ij_end
122  if (pole_nord) ijb=ijb+iip1
123  if (pole_sud) ije=ije-iip1
124 
125 
126  DO l=1,llm
127  DO ij = ijb,ije
128  mu(ij,l)=pbaru(ij,l) * zzpbar
129  ENDDO
130  ENDDO
131 
132  ijb=ij_begin-iip1
133  ije=ij_end
134  if (pole_nord) ijb=ij_begin
135  if (pole_sud) ije=ij_end-iip1
136 
137  DO l=1,llm
138  DO ij=ijb,ije
139  mv(ij,l)=pbarv(ij,l) * zzpbar
140  ENDDO
141  ENDDO
142 
143  ijb=ij_begin
144  ije=ij_end
145 
146  DO l=1,llm
147  DO ij=ijb,ije
148  mw(ij,l)=w(ij,l) * zzw
149  ENDDO
150  ENDDO
151 
152  DO ij=ijb,ije
153  mw(ij,llm+1)=0.
154  ENDDO
155 
156 c CALL SCOPY(ijp1llm,q,1,zq,1)
157 c CALL SCOPY(ijp1llm,masse,1,zm,1)
158 
159  ijb=ij_begin
160  ije=ij_end
161  zq(ijb:ije,1:llm)=q(ijb:ije,1:llm)
162  zm(ijb:ije,1:llm)=masse(ijb:ije,1:llm)
163 
164 
165  call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_begin,ij_begin+2*iip1-1)
166  call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_end-2*iip1+1,ij_end)
167 
168  call vtb(vthallo)
169  call register_hallo(zq,ip1jmp1,llm,2,2,2,2,myrequest1)
170  call register_hallo(zm,ip1jmp1,llm,1,1,1,1,myrequest1)
171  call sendrequest(myrequest1)
172  call vte(vthallo)
173 
174  call vlxqs_p(zq,pente_max,zm,mu,qsat,
175  . ij_begin+2*iip1,ij_end-2*iip1)
176 
177  call vtb(vthallo)
178  call waitrecvrequest(myrequest1)
179  call vte(vthallo)
180 
181  call vlyqs_p(zq,pente_max,zm,mv,qsat)
182 
183  call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
184  call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
185 
186  call vtb(vthallo)
187  call register_hallo(zq,ip1jmp1,llm,2,2,2,2,myrequest2)
188  call register_hallo(zm,ip1jmp1,llm,1,1,1,1,myrequest2)
189  call sendrequest(myrequest2)
190  call vte(vthallo)
191 
192  call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1)
193 
194  call vtb(vthallo)
195  call waitrecvrequest(myrequest2)
196  call vte(vthallo)
197 
198  call vlyqs_p(zq,pente_max,zm,mv,qsat)
199 
200 
201  call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_begin,ij_end)
202 
203 
204  ijb=ij_begin
205  ije=ij_end
206 
207  DO l=1,llm
208  DO ij=ijb,ije
209  q(ij,l)=zq(ij,l)
210  ENDDO
211  ENDDO
212 
213  DO l=1,llm
214  DO ij=ijb,ije-iip1+1,iip1
215  q(ij+iim,l)=q(ij,l)
216  ENDDO
217  ENDDO
218 
219  call waitsendrequest(myrequest1)
220  call waitsendrequest(myrequest2)
221 
222  RETURN
223  END
224  SUBROUTINE vlxqs_p(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
225 c
226 c Auteurs: P.Le Van, F.Hourdin, F.Forget
227 c
228 c ********************************************************************
229 c Shema d'advection " pseudo amont " .
230 c ********************************************************************
231 c
232 c --------------------------------------------------------------------
233  USE parallel_lmdz
234  IMPLICIT NONE
235 c
236 #include "dimensions.h"
237 #include "paramet.h"
238 #include "logic.h"
239 #include "comvert.h"
240 #include "comconst.h"
241 c
242 c
243 c Arguments:
244 c ----------
245  REAL masse(ip1jmp1,llm),pente_max
246  REAL u_m( ip1jmp1,llm )
247  REAL q(ip1jmp1,llm)
248  REAL qsat(ip1jmp1,llm)
249 c
250 c Local
251 c ---------
252 c
253  INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
254  INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
255 c
256  REAL new_m,zu_m,zdum(ip1jmp1,llm)
257  REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
258  REAL zz(ip1jmp1)
259  REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
260  REAL u_mq(ip1jmp1,llm)
261 
262  REAL SSUM
263 
264 
265  INTEGER ijb,ije,ijb_x,ije_x
266 
267 
268 c calcul de la pente a droite et a gauche de la maille
269 
270 c ijb=ij_begin
271 c ije=ij_end
272 
273  ijb=ijb_x
274  ije=ije_x
275 
276  if (pole_nord.and.ijb==1) ijb=ijb+iip1
277  if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
278 
279  IF (pente_max.gt.-1.e-5) THEN
280 c IF (pente_max.gt.10) THEN
281 
282 c calcul des pentes avec limitation, Van Leer scheme I:
283 c -----------------------------------------------------
284 
285 c calcul de la pente aux points u
286 
287 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
288  DO l = 1, llm
289  DO ij=ijb,ije-1
290  dxqu(ij)=q(ij+1,l)-q(ij,l)
291 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
292 c sigu(ij)=u_m(ij,l)/masse(ij,l)
293  ENDDO
294  DO ij=ijb+iip1-1,ije,iip1
295  dxqu(ij)=dxqu(ij-iim)
296 c sigu(ij)=sigu(ij-iim)
297  ENDDO
298 
299  DO ij=ijb,ije
300  adxqu(ij)=abs(dxqu(ij))
301  ENDDO
302 
303 c calcul de la pente maximum dans la maille en valeur absolue
304 
305  DO ij=ijb+1,ije
306  dxqmax(ij,l)=pente_max*
307  , min(adxqu(ij-1),adxqu(ij))
308 c limitation subtile
309 c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
310 
311 
312  ENDDO
313 
314  DO ij=ijb+iip1-1,ije,iip1
315  dxqmax(ij-iim,l)=dxqmax(ij,l)
316  ENDDO
317 
318  DO ij=ijb+1,ije
319 #ifdef CRAY
320  dxq(ij,l)=
321  , cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
322 #else
323  IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
324  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
325  ELSE
326 c extremum local
327  dxq(ij,l)=0.
328  ENDIF
329 #endif
330  dxq(ij,l)=0.5*dxq(ij,l)
331  dxq(ij,l)=
332  , sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
333  ENDDO
334 
335  ENDDO ! l=1,llm
336 c$OMP END DO NOWAIT
337 
338  ELSE ! (pente_max.lt.-1.e-5)
339 
340 c Pentes produits:
341 c ----------------
342 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
343  DO l = 1, llm
344  DO ij=ijb,ije-1
345  dxqu(ij)=q(ij+1,l)-q(ij,l)
346  ENDDO
347  DO ij=ijb+iip1-1,ije,iip1
348  dxqu(ij)=dxqu(ij-iim)
349  ENDDO
350 
351  DO ij=ijb+1,ije
352  zz(ij)=dxqu(ij-1)*dxqu(ij)
353  zz(ij)=zz(ij)+zz(ij)
354  IF(zz(ij).gt.0) THEN
355  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
356  ELSE
357 c extremum local
358  dxq(ij,l)=0.
359  ENDIF
360  ENDDO
361 
362  ENDDO
363 c$OMP END DO NOWAIT
364  ENDIF ! (pente_max.lt.-1.e-5)
365 
366 c bouclage de la pente en iip1:
367 c -----------------------------
368 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
369  DO l=1,llm
370  DO ij=ijb+iip1-1,ije,iip1
371  dxq(ij-iim,l)=dxq(ij,l)
372  ENDDO
373 
374  DO ij=ijb,ije
375  iadvplus(ij,l)=0
376  ENDDO
377 
378  ENDDO
379 c$OMP END DO NOWAIT
380 
381  if (pole_nord) THEN
382 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
383  DO l=1,llm
384  iadvplus(1:iip1,l)=0
385  ENDDO
386 c$OMP END DO NOWAIT
387  endif
388 
389  if (pole_sud) THEN
390 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
391  DO l=1,llm
392  iadvplus(ip1jm+1:ip1jmp1,l)=0
393  ENDDO
394 c$OMP END DO NOWAIT
395  endif
396 
397 c calcul des flux a gauche et a droite
398 
399 #ifdef CRAY
400 c--pas encore modification sur Qsat
401 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
402  DO l=1,llm
403  DO ij=ijb,ije-1
404  zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
405  , 1.+u_m(ij,l)/masse(ij+1,l),
406  , u_m(ij,l))
407  zdum(ij,l)=0.5*zdum(ij,l)
408  u_mq(ij,l)=cvmgp(
409  , q(ij,l)+zdum(ij,l)*dxq(ij,l),
410  , q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
411  , u_m(ij,l))
412  u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
413  ENDDO
414  ENDDO
415 c$OMP END DO NOWAIT
416 
417 #else
418 c on cumule le flux correspondant a toutes les mailles dont la masse
419 c au travers de la paroi pENDant le pas de temps.
420 c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
421 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
422  DO l=1,llm
423  DO ij=ijb,ije-1
424  IF (u_m(ij,l).gt.0.) THEN
425  zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
426  u_mq(ij,l)=u_m(ij,l)*
427  $ min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
428  ELSE
429  zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
430  u_mq(ij,l)=u_m(ij,l)*
431  $ min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
432  ENDIF
433  ENDDO
434  ENDDO
435 c$OMP END DO NOWAIT
436 #endif
437 
438 
439 c detection des points ou on advecte plus que la masse de la
440 c maille
441 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
442  DO l=1,llm
443  DO ij=ijb,ije-1
444  IF(zdum(ij,l).lt.0) THEN
445  iadvplus(ij,l)=1
446  u_mq(ij,l)=0.
447  ENDIF
448  ENDDO
449  ENDDO
450 c$OMP END DO NOWAIT
451 
452 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
453  DO l=1,llm
454  DO ij=ijb+iip1-1,ije,iip1
455  iadvplus(ij,l)=iadvplus(ij-iim,l)
456  ENDDO
457  ENDDO
458 c$OMP END DO NOWAIT
459 
460 
461 
462 c traitement special pour le cas ou on advecte en longitude plus que le
463 c contenu de la maille.
464 c cette partie est mal vectorisee.
465 
466 c pas d'influence de la pression saturante (pour l'instant)
467 
468 c calcul du nombre de maille sur lequel on advecte plus que la maille.
469 
470  n0=0
471 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
472  DO l=1,llm
473  nl(l)=0
474  DO ij=ijb,ije
475  nl(l)=nl(l)+iadvplus(ij,l)
476  ENDDO
477  n0=n0+nl(l)
478  ENDDO
479 c$OMP END DO NOWAIT
480 
481 cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
482 cym IF(n0.gt.1) THEN
483 cym IF(n0.gt.0) THEN
484 ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
485 ccc & ,'contenu de la maille : ',n0
486 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
487  DO l=1,llm
488  IF(nl(l).gt.0) THEN
489  iju=0
490 c indicage des mailles concernees par le traitement special
491  DO ij=ijb,ije
492  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
493  iju=iju+1
494  indu(iju)=ij
495  ENDIF
496  ENDDO
497  niju=iju
498 c PRINT*,'niju,nl',niju,nl(l)
499 
500 c traitement des mailles
501  DO iju=1,niju
502  ij=indu(iju)
503  j=(ij-1)/iip1+1
504  zu_m=u_m(ij,l)
505  u_mq(ij,l)=0.
506  IF(zu_m.gt.0.) THEN
507  ijq=ij
508  i=ijq-(j-1)*iip1
509 c accumulation pour les mailles completements advectees
510  do while(zu_m.gt.masse(ijq,l))
511  u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
512  zu_m=zu_m-masse(ijq,l)
513  i=mod(i-2+iim,iim)+1
514  ijq=(j-1)*iip1+i
515  ENDDO
516 c ajout de la maille non completement advectee
517  u_mq(ij,l)=u_mq(ij,l)+zu_m*
518  & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
519  ELSE
520  ijq=ij+1
521  i=ijq-(j-1)*iip1
522 c accumulation pour les mailles completements advectees
523  do while(-zu_m.gt.masse(ijq,l))
524  u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
525  zu_m=zu_m+masse(ijq,l)
526  i=mod(i,iim)+1
527  ijq=(j-1)*iip1+i
528  ENDDO
529 c ajout de la maille non completement advectee
530  u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
531  & 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
532  ENDIF
533  ENDDO
534  ENDIF
535  ENDDO
536 c$OMP END DO NOWAIT
537 cym ENDIF ! n0.gt.0
538 
539 
540 
541 c bouclage en latitude
542 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
543  DO l=1,llm
544  DO ij=ijb+iip1-1,ije,iip1
545  u_mq(ij,l)=u_mq(ij-iim,l)
546  ENDDO
547  ENDDO
548 c$OMP END DO NOWAIT
549 
550 c calcul des tendances
551 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
552  DO l=1,llm
553  DO ij=ijb+1,ije
554  new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
555  q(ij,l)=(q(ij,l)*masse(ij,l)+
556  & u_mq(ij-1,l)-u_mq(ij,l))
557  & /new_m
558  masse(ij,l)=new_m
559  ENDDO
560 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
561  DO ij=ijb+iip1-1,ije,iip1
562  q(ij-iim,l)=q(ij,l)
563  masse(ij-iim,l)=masse(ij,l)
564  ENDDO
565  ENDDO
566 c$OMP END DO NOWAIT
567 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
568 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
569 
570 
571  RETURN
572  END
573  SUBROUTINE vlyqs_p(q,pente_max,masse,masse_adv_v,qsat)
574 c
575 c Auteurs: P.Le Van, F.Hourdin, F.Forget
576 c
577 c ********************************************************************
578 c Shema d'advection " pseudo amont " .
579 c ********************************************************************
580 c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
581 c qsat est un argument de sortie pour le s-pg ....
582 c
583 c
584 c --------------------------------------------------------------------
585  USE parallel_lmdz
586  IMPLICIT NONE
587 c
588 #include "dimensions.h"
589 #include "paramet.h"
590 #include "logic.h"
591 #include "comvert.h"
592 #include "comconst.h"
593 #include "comgeom.h"
594 c
595 c
596 c Arguments:
597 c ----------
598  REAL masse(ip1jmp1,llm),pente_max
599  REAL masse_adv_v( ip1jm,llm)
600  REAL q(ip1jmp1,llm)
601  REAL qsat(ip1jmp1,llm)
602 c
603 c Local
604 c ---------
605 c
606  INTEGER i,ij,l
607 c
608  REAL airej2,airejjm,airescb(iim),airesch(iim)
609  REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
610  REAL adyqv(ip1jm),dyqmax(ip1jmp1)
611  REAL qbyv(ip1jm,llm)
612 
613  REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
614 c REAL newq,oldmasse
615  Logical first
616  SAVE first
617 c$OMP THREADPRIVATE(first)
618  REAL convpn,convps,convmpn,convmps
619  REAL sinlon(iip1),sinlondlon(iip1)
620  REAL coslon(iip1),coslondlon(iip1)
621  SAVE sinlon,coslon,sinlondlon,coslondlon
622  SAVE airej2,airejjm
623 c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
624 c$OMP THREADPRIVATE(airej2,airejjm)
625 c
626 c
627  REAL SSUM
628 
629  DATA first/.true./
630  INTEGER ijb,ije
631 
632  IF(first) THEN
633  print*,'Shema Amont nouveau appele dans Vanleer '
634  first=.false.
635  do i=2,iip1
636  coslon(i)=cos(rlonv(i))
637  sinlon(i)=sin(rlonv(i))
638  coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
639  sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
640  ENDDO
641  coslon(1)=coslon(iip1)
642  coslondlon(1)=coslondlon(iip1)
643  sinlon(1)=sinlon(iip1)
644  sinlondlon(1)=sinlondlon(iip1)
645  airej2 = ssum( iim, aire(iip2), 1 )
646  airejjm= ssum( iim, aire(ip1jm -iim), 1 )
647  ENDIF
648 
649 c
650 
651 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
652  DO l = 1, llm
653 c
654 c --------------------------------
655 c CALCUL EN LATITUDE
656 c --------------------------------
657 
658 c On commence par calculer la valeur du traceur moyenne sur le premier cercle
659 c de latitude autour du pole (qpns pour le pole nord et qpsn pour
660 c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
661 
662  if (pole_nord) then
663  DO i = 1, iim
664  airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
665  ENDDO
666  qpns = ssum( iim, airescb ,1 ) / airej2
667  endif
668 
669  if (pole_sud) then
670  DO i = 1, iim
671  airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
672  ENDDO
673  qpsn = ssum( iim, airesch ,1 ) / airejjm
674  endif
675 
676 
677 c calcul des pentes aux points v
678 
679  ijb=ij_begin-2*iip1
680  ije=ij_end+iip1
681  if (pole_nord) ijb=ij_begin
682  if (pole_sud) ije=ij_end-iip1
683 
684  DO ij=ijb,ije
685  dyqv(ij)=q(ij,l)-q(ij+iip1,l)
686  adyqv(ij)=abs(dyqv(ij))
687  ENDDO
688 
689 
690 c calcul des pentes aux points scalaires
691 
692  ijb=ij_begin-iip1
693  ije=ij_end+iip1
694  if (pole_nord) ijb=ij_begin+iip1
695  if (pole_sud) ije=ij_end-iip1
696 
697  DO ij=ijb,ije
698  dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
699  dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
700  dyqmax(ij)=pente_max*dyqmax(ij)
701  ENDDO
702 
703  IF (pole_nord) THEN
704 
705 c calcul des pentes aux poles
706  DO ij=1,iip1
707  dyq(ij,l)=qpns-q(ij+iip1,l)
708  ENDDO
709 
710 c filtrage de la derivee
711  dyn1=0.
712  dyn2=0.
713  DO ij=1,iim
714  dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
715  dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
716  ENDDO
717  DO ij=1,iip1
718  dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
719  ENDDO
720 
721 c calcul des pentes limites aux poles
722  fn=1.
723  DO ij=1,iim
724  IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
725  fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
726  ENDIF
727  ENDDO
728 
729  DO ij=1,iip1
730  dyq(ij,l)=fn*dyq(ij,l)
731  ENDDO
732 
733  ENDIF
734 
735  IF (pole_sud) THEN
736 
737  DO ij=1,iip1
738  dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
739  ENDDO
740 
741  dys1=0.
742  dys2=0.
743 
744  DO ij=1,iim
745  dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
746  dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
747  ENDDO
748 
749  DO ij=1,iip1
750  dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
751  ENDDO
752 
753 c calcul des pentes limites aux poles
754  fs=1.
755  DO ij=1,iim
756  IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
757  fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
758  ENDIF
759  ENDDO
760 
761  DO ij=1,iip1
762  dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
763  ENDDO
764 
765  ENDIF
766 
767 
768 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
769 C En memoire de dIFferents tests sur la
770 C limitation des pentes aux poles.
771 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
772 C PRINT*,dyq(1)
773 C PRINT*,dyqv(iip1+1)
774 C appn=abs(dyq(1)/dyqv(iip1+1))
775 C PRINT*,dyq(ip1jm+1)
776 C PRINT*,dyqv(ip1jm-iip1+1)
777 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
778 C DO ij=2,iim
779 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
780 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
781 C ENDDO
782 C appn=min(pente_max/appn,1.)
783 C apps=min(pente_max/apps,1.)
784 C
785 C
786 C cas ou on a un extremum au pole
787 C
788 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
789 C & appn=0.
790 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
791 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
792 C & apps=0.
793 C
794 C limitation des pentes aux poles
795 C DO ij=1,iip1
796 C dyq(ij)=appn*dyq(ij)
797 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
798 C ENDDO
799 C
800 C test
801 C DO ij=1,iip1
802 C dyq(iip1+ij)=0.
803 C dyq(ip1jm+ij-iip1)=0.
804 C ENDDO
805 C DO ij=1,ip1jmp1
806 C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
807 C ENDDO
808 C
809 C changement 10 07 96
810 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
811 C & THEN
812 C DO ij=1,iip1
813 C dyqmax(ij)=0.
814 C ENDDO
815 C ELSE
816 C DO ij=1,iip1
817 C dyqmax(ij)=pente_max*abs(dyqv(ij))
818 C ENDDO
819 C ENDIF
820 C
821 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
822 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
823 C &THEN
824 C DO ij=ip1jm+1,ip1jmp1
825 C dyqmax(ij)=0.
826 C ENDDO
827 C ELSE
828 C DO ij=ip1jm+1,ip1jmp1
829 C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
830 C ENDDO
831 C ENDIF
832 C fin changement 10 07 96
833 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
834 
835 c calcul des pentes limitees
836  ijb=ij_begin-iip1
837  ije=ij_end+iip1
838  if (pole_nord) ijb=ij_begin+iip1
839  if (pole_sud) ije=ij_end-iip1
840 
841  DO ij=ijb,ije
842  IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
843  dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
844  ELSE
845  dyq(ij,l)=0.
846  ENDIF
847  ENDDO
848 
849  ENDDO
850 c$OMP END DO NOWAIT
851 
852  ijb=ij_begin-iip1
853  ije=ij_end
854  if (pole_nord) ijb=ij_begin
855  if (pole_sud) ije=ij_end-iip1
856 
857 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
858  DO l=1,llm
859  DO ij=ijb,ije
860  IF( masse_adv_v(ij,l).GT.0. ) THEN
861  qbyv(ij,l)= min( qsat(ij+iip1,l), q(ij+iip1,l ) +
862  , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
863  ELSE
864  qbyv(ij,l)= min( qsat(ij,l), q(ij,l) - dyq(ij,l) *
865  , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
866  ENDIF
867  qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
868  ENDDO
869  ENDDO
870 c$OMP END DO NOWAIT
871 
872  ijb=ij_begin
873  ije=ij_end
874  if (pole_nord) ijb=ij_begin+iip1
875  if (pole_sud) ije=ij_end-iip1
876 
877 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
878  DO l=1,llm
879  DO ij=ijb,ije
880  newmasse=masse(ij,l)
881  & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
882  q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
883  & /newmasse
884  masse(ij,l)=newmasse
885  ENDDO
886 c.-. ancienne version
887 
888  IF (pole_nord) THEN
889 
890  convpn=ssum(iim,qbyv(1,l),1)/apoln
891  convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
892  DO ij = 1,iip1
893  newmasse=masse(ij,l)+convmpn*aire(ij)
894  q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
895  & newmasse
896  masse(ij,l)=newmasse
897  ENDDO
898 
899  ENDIF
900 
901  IF (pole_sud) THEN
902 
903  convps = -ssum(iim,qbyv(ip1jm-iim,l),1)/apols
904  convmps = -ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
905  DO ij = ip1jm+1,ip1jmp1
906  newmasse=masse(ij,l)+convmps*aire(ij)
907  q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
908  & newmasse
909  masse(ij,l)=newmasse
910  ENDDO
911 
912  ENDIF
913 c.-. fin ancienne version
914 
915 c._. nouvelle version
916 c convpn=SSUM(iim,qbyv(1,l),1)
917 c convmpn=ssum(iim,masse_adv_v(1,l),1)
918 c oldmasse=ssum(iim,masse(1,l),1)
919 c newmasse=oldmasse+convmpn
920 c newq=(q(1,l)*oldmasse+convpn)/newmasse
921 c newmasse=newmasse/apoln
922 c DO ij = 1,iip1
923 c q(ij,l)=newq
924 c masse(ij,l)=newmasse*aire(ij)
925 c ENDDO
926 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
927 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
928 c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
929 c newmasse=oldmasse+convmps
930 c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
931 c newmasse=newmasse/apols
932 c DO ij = ip1jm+1,ip1jmp1
933 c q(ij,l)=newq
934 c masse(ij,l)=newmasse*aire(ij)
935 c ENDDO
936 c._. fin nouvelle version
937  ENDDO
938 c$OMP END DO NOWAIT
939  RETURN
940  END
iip2
!$Header iip2
Definition: paramet.h:14
ip1jmp1
!$Header llmm1 INTEGER ip1jmp1
Definition: paramet.h:14
apols
!$Header!CDK comgeom COMMON comgeom apols
Definition: comgeom.h:8
mod_hallo::register_hallo
subroutine register_hallo(Field, ij, ll, RUp, Rdown, SUp, SDown, a_request)
Definition: mod_hallo.F90:875
llmp1
!$Header llmp1
Definition: paramet.h:14
vampir
Definition: vampir.F90:1
parallel_lmdz::ij_end
integer, save ij_end
Definition: parallel_lmdz.F90:23
parallel_lmdz::pole_sud
logical, save pole_sud
Definition: parallel_lmdz.F90:25
vlspltqs_p
subroutine vlspltqs_p(q, pente_max, masse, w, pbaru, pbarv, pdt, p, pk, teta)
Definition: vlspltqs_p.F:6
vampir::vtb
subroutine vtb(number)
Definition: vampir.F90:52
llm
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm day day1 day day1 *dt_dice endif!time annee_ref dt_dice swup_dice vg_dice omega_dice tg_prof vg_profd w_profd omega_profd!do llm!print llm l llm
Definition: 1D_interp_cases.h:103
pi
!$Id mode_top_bound COMMON comconstr && pi
Definition: comconst.h:7
aire
!$Header!CDK comgeom COMMON comgeom aire
Definition: comgeom.h:25
ip1jm
!$Header llmm1 INTEGER ip1jm
Definition: paramet.h:14
vlxqs_p
subroutine vlxqs_p(q, pente_max, masse, u_m, qsat, ijb_x, ije_x)
Definition: vlspltqs_p.F:225
apoln
!$Header!CDK comgeom COMMON comgeom apoln
Definition: comgeom.h:8
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
parallel_lmdz::pole_nord
logical, save pole_nord
Definition: parallel_lmdz.F90:24
parallel_lmdz
Definition: parallel_lmdz.F90:4
cpp
!$Id mode_top_bound COMMON comconstr cpp
Definition: comconst.h:7
vampir::vthallo
integer, parameter vthallo
Definition: vampir.F90:7
vlyqs_p
subroutine vlyqs_p(q, pente_max, masse, masse_adv_v, qsat)
Definition: vlspltqs_p.F:574
rlonu
!$Header!CDK comgeom COMMON comgeom rlonu
Definition: comgeom.h:25
mod_hallo::sendrequest
subroutine sendrequest(a_Request)
Definition: mod_hallo.F90:1072
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
parallel_lmdz::ij_begin
integer, save ij_begin
Definition: parallel_lmdz.F90:22
mod_hallo::waitrecvrequest
subroutine waitrecvrequest(a_Request)
Definition: mod_hallo.F90:1346
vlz_p
subroutine vlz_p(q, pente_max, masse, w, ijb_x, ije_x)
Definition: vlsplt_p.F:918
vampir::vte
subroutine vte(number)
Definition: vampir.F90:69
mod_hallo::waitsendrequest
subroutine waitsendrequest(a_Request)
Definition: mod_hallo.F90:1290
iim
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
mod_hallo::settag
subroutine settag(a_request, tag)
Definition: mod_hallo.F90:180
mod_hallo::request
Definition: mod_hallo.F90:38
mod_hallo
Definition: mod_hallo.F90:1
rlonv
!$Header!CDK comgeom COMMON comgeom rlonv
Definition: comgeom.h:25