My Project
 All Classes Files Functions Variables Macros
vlspltqs_p.F
Go to the documentation of this file.
1 c
2 c $Id: vlspltqs_p.F 1520 2011-05-23 11:37:09Z emillour $
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
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
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
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