source: LMDZ5/trunk/libf/dyn3dpar/vlspltqs_p.F @ 1786

Last change on this file since 1786 was 1520, checked in by Ehouarn Millour, 14 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.6 KB
Line 
1c
2c $Id: vlspltqs_p.F 1520 2011-05-23 11:37:09Z idelkadi $
3c
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      INTEGER ijb,ije
73      type(request) :: MyRequest1
74      type(request) :: MyRequest2
75
76c    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
90c-- Calcul de Qsat en chaque point
91c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
92c   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
115c      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
156c      CALL SCOPY(ijp1llm,q,1,zq,1)
157c      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)
225c
226c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
227c
228c    ********************************************************************
229c     Shema  d'advection " pseudo amont " .
230c    ********************************************************************
231c
232c   --------------------------------------------------------------------
233      USE parallel
234      IMPLICIT NONE
235c
236#include "dimensions.h"
237#include "paramet.h"
238#include "logic.h"
239#include "comvert.h"
240#include "comconst.h"
241c
242c
243c   Arguments:
244c   ----------
245      REAL masse(ip1jmp1,llm),pente_max
246      REAL u_m( ip1jmp1,llm )
247      REAL q(ip1jmp1,llm)
248      REAL qsat(ip1jmp1,llm)
249c
250c      Local
251c   ---------
252c
253      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
254      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
255c
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
268c   calcul de la pente a droite et a gauche de la maille
269
270c      ijb=ij_begin
271c      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
280c     IF (pente_max.gt.10) THEN
281
282c   calcul des pentes avec limitation, Van Leer scheme I:
283c   -----------------------------------------------------
284
285c   calcul de la pente aux points u
286
287c$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)
291c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
292c              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)
296c              sigu(ij)=sigu(ij-iim)
297            ENDDO
298
299            DO ij=ijb,ije
300               adxqu(ij)=abs(dxqu(ij))
301            ENDDO
302
303c   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))
308c limitation subtile
309c    ,      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
326c   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
336c$OMP END DO NOWAIT
337
338      ELSE ! (pente_max.lt.-1.e-5)
339
340c   Pentes produits:
341c   ----------------
342c$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
357c   extremum local
358                  dxq(ij,l)=0.
359               ENDIF
360            ENDDO
361
362         ENDDO
363c$OMP END DO NOWAIT
364      ENDIF ! (pente_max.lt.-1.e-5)
365
366c   bouclage de la pente en iip1:
367c   -----------------------------
368c$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
379c$OMP END DO NOWAIT
380     
381      if (pole_nord) THEN
382c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
383        DO l=1,llm     
384          iadvplus(1:iip1,l)=0
385        ENDDO
386c$OMP END DO NOWAIT
387      endif
388     
389      if (pole_sud)  THEN
390c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
391        DO l=1,llm
392          iadvplus(ip1jm+1:ip1jmp1,l)=0
393        ENDDO
394c$OMP END DO NOWAIT
395      endif
396       
397c   calcul des flux a gauche et a droite
398
399#ifdef CRAY
400c--pas encore modification sur Qsat
401c$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
415c$OMP END DO NOWAIT
416
417#else
418c   on cumule le flux correspondant a toutes les mailles dont la masse
419c   au travers de la paroi pENDant le pas de temps.
420c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
421c$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
435c$OMP END DO NOWAIT
436#endif
437
438
439c   detection des points ou on advecte plus que la masse de la
440c   maille
441c$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
450c$OMP END DO NOWAIT
451
452c$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
458c$OMP END DO NOWAIT
459
460
461
462c   traitement special pour le cas ou on advecte en longitude plus que le
463c   contenu de la maille.
464c   cette partie est mal vectorisee.
465
466c   pas d'influence de la pression saturante (pour l'instant)
467
468c  calcul du nombre de maille sur lequel on advecte plus que la maille.
469
470      n0=0
471c$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
479c$OMP END DO NOWAIT
480
481cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
482cym      IF(n0.gt.1) THEN
483cym        IF(n0.gt.0) THEN
484ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
485ccc     &       ,'contenu de la maille : ',n0
486c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
487         DO l=1,llm
488            IF(nl(l).gt.0) THEN
489               iju=0
490c   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
498c              PRINT*,'niju,nl',niju,nl(l)
499
500c  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
509c   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
516c   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
522c   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
529c   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
536c$OMP END DO NOWAIT
537cym      ENDIF  ! n0.gt.0
538
539
540
541c   bouclage en latitude
542c$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
548c$OMP END DO NOWAIT
549
550c   calcul des tendances
551c$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
560c   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
566c$OMP END DO NOWAIT
567c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
568c     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)
574c
575c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
576c
577c    ********************************************************************
578c     Shema  d'advection " pseudo amont " .
579c    ********************************************************************
580c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
581c     qsat             est   un argument de sortie pour le s-pg ....
582c
583c
584c   --------------------------------------------------------------------
585      USE parallel
586      IMPLICIT NONE
587c
588#include "dimensions.h"
589#include "paramet.h"
590#include "logic.h"
591#include "comvert.h"
592#include "comconst.h"
593#include "comgeom.h"
594c
595c
596c   Arguments:
597c   ----------
598      REAL masse(ip1jmp1,llm),pente_max
599      REAL masse_adv_v( ip1jm,llm)
600      REAL q(ip1jmp1,llm)
601      REAL qsat(ip1jmp1,llm)
602c
603c      Local
604c   ---------
605c
606      INTEGER i,ij,l
607c
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
614c     REAL newq,oldmasse
615      Logical first
616      SAVE first
617c$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
623c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
624c$OMP THREADPRIVATE(airej2,airejjm)
625c
626c
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
649c
650
651c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
652      DO l = 1, llm
653c
654c   --------------------------------
655c      CALCUL EN LATITUDE
656c   --------------------------------
657
658c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
659c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
660c    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
677c   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
690c   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
705c   calcul des pentes aux poles
706        DO ij=1,iip1
707           dyq(ij,l)=qpns-q(ij+iip1,l)
708        ENDDO
709
710c   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
721c   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       
753c   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
768CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
769C  En memoire de dIFferents tests sur la
770C  limitation des pentes aux poles.
771CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
772C     PRINT*,dyq(1)
773C     PRINT*,dyqv(iip1+1)
774C     appn=abs(dyq(1)/dyqv(iip1+1))
775C     PRINT*,dyq(ip1jm+1)
776C     PRINT*,dyqv(ip1jm-iip1+1)
777C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
778C     DO ij=2,iim
779C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
780C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
781C     ENDDO
782C     appn=min(pente_max/appn,1.)
783C     apps=min(pente_max/apps,1.)
784C
785C
786C   cas ou on a un extremum au pole
787C
788C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
789C    &   appn=0.
790C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
791C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
792C    &   apps=0.
793C
794C   limitation des pentes aux poles
795C     DO ij=1,iip1
796C        dyq(ij)=appn*dyq(ij)
797C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
798C     ENDDO
799C
800C   test
801C      DO ij=1,iip1
802C         dyq(iip1+ij)=0.
803C         dyq(ip1jm+ij-iip1)=0.
804C      ENDDO
805C      DO ij=1,ip1jmp1
806C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
807C      ENDDO
808C
809C changement 10 07 96
810C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
811C    &   THEN
812C        DO ij=1,iip1
813C           dyqmax(ij)=0.
814C        ENDDO
815C     ELSE
816C        DO ij=1,iip1
817C           dyqmax(ij)=pente_max*abs(dyqv(ij))
818C        ENDDO
819C     ENDIF
820C
821C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
822C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
823C    &THEN
824C        DO ij=ip1jm+1,ip1jmp1
825C           dyqmax(ij)=0.
826C        ENDDO
827C     ELSE
828C        DO ij=ip1jm+1,ip1jmp1
829C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
830C        ENDDO
831C     ENDIF
832C   fin changement 10 07 96
833CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
834
835c   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
850c$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
857c$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
870c$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
877c$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
886c.-. 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
913c.-. fin ancienne version
914
915c._. nouvelle version
916c        convpn=SSUM(iim,qbyv(1,l),1)
917c        convmpn=ssum(iim,masse_adv_v(1,l),1)
918c        oldmasse=ssum(iim,masse(1,l),1)
919c        newmasse=oldmasse+convmpn
920c        newq=(q(1,l)*oldmasse+convpn)/newmasse
921c        newmasse=newmasse/apoln
922c        DO ij = 1,iip1
923c           q(ij,l)=newq
924c           masse(ij,l)=newmasse*aire(ij)
925c        ENDDO
926c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
927c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
928c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
929c        newmasse=oldmasse+convmps
930c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
931c        newmasse=newmasse/apols
932c        DO ij = ip1jm+1,ip1jmp1
933c           q(ij,l)=newq
934c           masse(ij,l)=newmasse*aire(ij)
935c        ENDDO
936c._. fin nouvelle version
937      ENDDO
938c$OMP END DO NOWAIT
939      RETURN
940      END
Note: See TracBrowser for help on using the repository browser.