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

Last change on this file since 4735 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
EM

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