source: trunk/libf/dyn3dpar/vlspltqs_p.F @ 8

Last change on this file since 8 was 8, checked in by emillour, 14 years ago

Debut de mise a jour de la dynamique parallele par rapport aux modifs dans la partie sequentielle.

Mais NON TESTE , car pas (encore) possibilite de compiler et faire tourner cas simple (type newtonien sans physique).

Voir commit_v8.log pour les details.

Ehouarn

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