source: LMDZ5/trunk/libf/dyn3dpar/vlsplt_p.F @ 5437

Last change on this file since 5437 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: 29.6 KB
RevLine 
[1520]1c
2c $Id: vlsplt_p.F 2603 2016-07-25 09:31:56Z evignon $
3c
4
[630]5      SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
6c
7c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
8c
9c    ********************************************************************
10c     Shema  d'advection " pseudo amont " .
11c    ********************************************************************
12c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
13c
14c   pente_max facteur de limitation des pentes: 2 en general
15c                                               0 pour un schema amont
16c   pbaru,pbarv,w flux de masse en u ,v ,w
17c   pdt pas de temps
18c
19c   --------------------------------------------------------------------
[1823]20      USE parallel_lmdz
[630]21      USE mod_hallo
22      USE Vampir
23      IMPLICIT NONE
24c
25#include "dimensions.h"
26#include "paramet.h"
27
28c
29c   Arguments:
30c   ----------
31      REAL masse(ip1jmp1,llm),pente_max
32c      REAL masse(iip1,jjp1,llm),pente_max
33      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
34      REAL q(ip1jmp1,llm)
35c      REAL q(iip1,jjp1,llm)
36      REAL w(ip1jmp1,llm),pdt
37c
38c      Local
39c   ---------
40c
41      INTEGER i,ij,l,j,ii
42      INTEGER ijlqmin,iqmin,jqmin,lqmin
43c
44      REAL zm(ip1jmp1,llm),newmasse
45      REAL mu(ip1jmp1,llm)
46      REAL mv(ip1jm,llm)
47      REAL mw(ip1jmp1,llm+1)
48      REAL zq(ip1jmp1,llm),zz
49      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
50      REAL second,temps0,temps1,temps2,temps3
51      REAL ztemps1,ztemps2,ztemps3
52      REAL zzpbar, zzw
53      LOGICAL testcpu
54      SAVE testcpu
55      SAVE temps1,temps2,temps3
56      INTEGER iminn,imaxx
57
58      REAL qmin,qmax
59      DATA qmin,qmax/0.,1.e33/
60      DATA testcpu/.false./
61      DATA temps1,temps2,temps3/0.,0.,0./
62      INTEGER ijb,ije
63      type(request) :: MyRequest1
64      type(request) :: MyRequest2
65
66      call SetTag(MyRequest1,100)
67      call SetTag(MyRequest2,101)
68     
69      zzpbar = 0.5 * pdt
70      zzw    = pdt
71     
72      ijb=ij_begin
73      ije=ij_end
74      if (pole_nord) ijb=ijb+iip1
75      if (pole_sud)  ije=ije-iip1
76     
77      DO l=1,llm
78        DO ij = ijb,ije
79            mu(ij,l)=pbaru(ij,l) * zzpbar
80        ENDDO
81      ENDDO
82     
83      ijb=ij_begin-iip1
84      ije=ij_end
85      if (pole_nord) ijb=ij_begin
86      if (pole_sud)  ije=ij_end-iip1
87     
88      DO l=1,llm
89        DO ij=ijb,ije
90           mv(ij,l)=pbarv(ij,l) * zzpbar
91        ENDDO
92      ENDDO
93     
94      ijb=ij_begin
95      ije=ij_end
96     
97      DO l=1,llm
98        DO ij=ijb,ije
99           mw(ij,l)=w(ij,l) * zzw
100        ENDDO
101      ENDDO
102
103      DO ij=ijb,ije
104         mw(ij,llm+1)=0.
105      ENDDO
106     
107c      CALL SCOPY(ijp1llm,q,1,zq,1)
108c      CALL SCOPY(ijp1llm,masse,1,zm,1)
109       
110       ijb=ij_begin
111       ije=ij_end
112       zq(ijb:ije,:)=q(ijb:ije,:)
113       zm(ijb:ije,:)=masse(ijb:ije,:)
114     
115     
116c       print*,'Entree vlx1'
117c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
118      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1)
119      call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end)
120      call VTb(VTHallo)
121      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
122      call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest1)
123      call SendRequest(MyRequest1)
124      call VTe(VTHallo)
125      call vlx_p(zq,pente_max,zm,mu,ij_begin+2*iip1,ij_end-2*iip1)
126c      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
127      call VTb(VTHallo)
128      call WaitRecvRequest(MyRequest1)
129      call VTe(VTHallo)
130
131     
132c       print*,'Sortie vlx1'
133c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
134
135c        print*,'Entree vly1'
136c      call exchange_hallo(zq,ip1jmp1,llm,2,2)
137c      call exchange_hallo(zm,ip1jmp1,llm,1,1)
138     
139      call vly_p(zq,pente_max,zm,mv)
140c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
141c       print*,'Sortie vly1'
142      call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
143      call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
144      call VTb(VTHallo)
145      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
146      call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest2)
147      call SendRequest(MyRequest2)
148      call VTe(VTHallo)
149      call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1)
150      call VTb(VTHallo)
151      call WaitRecvRequest(MyRequest2)
152           
153      call VTe(VTHallo)
154     
155c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
156
157
158     
159     
160      call vly_p(zq,pente_max,zm,mv)
161c       call minmaxq(zq,qmin,qmax,'apres vly     ')
162
163
164      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
165c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
166
167       
168      ijb=ij_begin
169      ije=ij_end
170       
171      DO l=1,llm
172         DO ij=ijb,ije
173           q(ij,l)=zq(ij,l)
174         ENDDO
175      ENDDO
176     
177     
178      DO l=1,llm
179         DO ij=ijb,ije-iip1+1,iip1
180            q(ij+iim,l)=q(ij,l)
181         ENDDO
182      ENDDO
183
184      call WaitSendRequest(MyRequest1)
185      call WaitSendRequest(MyRequest2)
186     
187      RETURN
188      END
189     
190     
191      SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x)
192
193c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
194c
195c    ********************************************************************
196c     Shema  d'advection " pseudo amont " .
197c    ********************************************************************
198c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
199c
200c
201c   --------------------------------------------------------------------
[1823]202      USE parallel_lmdz
[630]203      IMPLICIT NONE
204c
205#include "dimensions.h"
206#include "paramet.h"
207c
208c
209c   Arguments:
210c   ----------
211      REAL masse(ip1jmp1,llm),pente_max
212      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
213      REAL q(ip1jmp1,llm)
214      REAL w(ip1jmp1,llm)
215c
216c      Local
217c   ---------
218c
219      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
220      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
221c
222      REAL new_m,zu_m,zdum(ip1jmp1,llm)
223      REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
224      REAL zz(ip1jmp1)
225      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
226      REAL u_mq(ip1jmp1,llm)
227
[764]228      Logical extremum
[630]229
230      REAL      SSUM
231      EXTERNAL  SSUM
232
233      REAL z1,z2,z3
234
235      INTEGER ijb,ije,ijb_x,ije_x
236     
237c   calcul de la pente a droite et a gauche de la maille
238
239      ijb=ijb_x
240      ije=ije_x
241       
242      if (pole_nord.and.ijb==1) ijb=ijb+iip1
243      if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
244         
245      IF (pente_max.gt.-1.e-5) THEN
246c       IF (pente_max.gt.10) THEN
247
248c   calcul des pentes avec limitation, Van Leer scheme I:
249c   -----------------------------------------------------
250
251c   calcul de la pente aux points u
[764]252c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
[630]253         DO l = 1, llm
254           
255            DO ij=ijb,ije-1
256               dxqu(ij)=q(ij+1,l)-q(ij,l)
257c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
258c              sigu(ij)=u_m(ij,l)/masse(ij,l)
259            ENDDO
260            DO ij=ijb+iip1-1,ije,iip1
261               dxqu(ij)=dxqu(ij-iim)
262c              sigu(ij)=sigu(ij-iim)
263            ENDDO
264
265            DO ij=ijb,ije
266               adxqu(ij)=abs(dxqu(ij))
267            ENDDO
268
269c   calcul de la pente maximum dans la maille en valeur absolue
270
271            DO ij=ijb+1,ije
272               dxqmax(ij,l)=pente_max*
273     ,      min(adxqu(ij-1),adxqu(ij))
274c limitation subtile
275c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
276         
277
278            ENDDO
279
280            DO ij=ijb+iip1-1,ije,iip1
281               dxqmax(ij-iim,l)=dxqmax(ij,l)
282            ENDDO
283
284            DO ij=ijb+1,ije
285#ifdef CRAY
286               dxq(ij,l)=
287     ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
288#else
289               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
290                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
291               ELSE
292c   extremum local
293                  dxq(ij,l)=0.
294               ENDIF
295#endif
296               dxq(ij,l)=0.5*dxq(ij,l)
297               dxq(ij,l)=
298     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
299            ENDDO
300
301         ENDDO ! l=1,llm
[764]302c$OMP END DO NOWAIT
[630]303c       print*,'Ok calcul des pentes'
304
305      ELSE ! (pente_max.lt.-1.e-5)
306
307c   Pentes produits:
308c   ----------------
[764]309c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]310         DO l = 1, llm
311            DO ij=ijb,ije-1
312               dxqu(ij)=q(ij+1,l)-q(ij,l)
313            ENDDO
314            DO ij=ijb+iip1-1,ije,iip1
315               dxqu(ij)=dxqu(ij-iim)
316            ENDDO
317
318            DO ij=ijb+1,ije
319               zz(ij)=dxqu(ij-1)*dxqu(ij)
320               zz(ij)=zz(ij)+zz(ij)
321               IF(zz(ij).gt.0) THEN
322                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
323               ELSE
324c   extremum local
325                  dxq(ij,l)=0.
326               ENDIF
327            ENDDO
328
329         ENDDO
[764]330c$OMP END DO NOWAIT
[630]331      ENDIF ! (pente_max.lt.-1.e-5)
332
333c   bouclage de la pente en iip1:
334c   -----------------------------
[764]335c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]336      DO l=1,llm
337         DO ij=ijb+iip1-1,ije,iip1
338            dxq(ij-iim,l)=dxq(ij,l)
339         ENDDO
[764]340         DO ij=ijb,ije
[630]341            iadvplus(ij,l)=0
342         ENDDO
343
344      ENDDO
[764]345c$OMP END DO NOWAIT
[630]346c        print*,'Bouclage en iip1'
347
348c   calcul des flux a gauche et a droite
349
350#ifdef CRAY
[764]351c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]352      DO l=1,llm
353       DO ij=ijb,ije-1
354          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
355     ,                     1.+u_m(ij,l)/masse(ij+1,l),
356     ,                     u_m(ij,l))
357          zdum(ij,l)=0.5*zdum(ij,l)
358          u_mq(ij,l)=cvmgp(
359     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
360     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
361     ,                u_m(ij,l))
362          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
363       ENDDO
364      ENDDO
[764]365c$OMP END DO NOWAIT
[630]366#else
367c   on cumule le flux correspondant a toutes les mailles dont la masse
368c   au travers de la paroi pENDant le pas de temps.
369c       print*,'Cumule ....'
[764]370c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]371      DO l=1,llm
372       DO ij=ijb,ije-1
373c       print*,'masse(',ij,')=',masse(ij,l)
374          IF (u_m(ij,l).gt.0.) THEN
375             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
376             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
377          ELSE
378             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
379             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
380          ENDIF
381       ENDDO
382      ENDDO
[764]383c$OMP END DO NOWAIT
[630]384#endif
385c       stop
386
387c       go to 9999
388c   detection des points ou on advecte plus que la masse de la
389c   maille
[764]390c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]391      DO l=1,llm
392         DO ij=ijb,ije-1
393            IF(zdum(ij,l).lt.0) THEN
394               iadvplus(ij,l)=1
395               u_mq(ij,l)=0.
396            ENDIF
397         ENDDO
398      ENDDO
[764]399c$OMP END DO NOWAIT
[630]400c       print*,'Ok test 1'
[764]401c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]402      DO l=1,llm
403       DO ij=ijb+iip1-1,ije,iip1
404          iadvplus(ij,l)=iadvplus(ij-iim,l)
405       ENDDO
406      ENDDO
[764]407c$OMP END DO NOWAIT
[630]408c        print*,'Ok test 2'
409
410
411c   traitement special pour le cas ou on advecte en longitude plus que le
412c   contenu de la maille.
413c   cette partie est mal vectorisee.
414
415c  calcul du nombre de maille sur lequel on advecte plus que la maille.
416
417      n0=0
[764]418c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]419      DO l=1,llm
420         nl(l)=0
421         DO ij=ijb,ije
422            nl(l)=nl(l)+iadvplus(ij,l)
423         ENDDO
424         n0=n0+nl(l)
425      ENDDO
[764]426c$OMP END DO NOWAIT
[630]427cym      IF(n0.gt.1) THEN
[764]428cym      IF(n0.gt.0) THEN
[630]429
430c      PRINT*,'Nombre de points pour lesquels on advect plus que le'
431c     &       ,'contenu de la maille : ',n0
[764]432c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]433         DO l=1,llm
434            IF(nl(l).gt.0) THEN
435               iju=0
436c   indicage des mailles concernees par le traitement special
437               DO ij=ijb,ije
438                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
439                     iju=iju+1
440                     indu(iju)=ij
441                  ENDIF
442               ENDDO
443               niju=iju
444c              PRINT*,'niju,nl',niju,nl(l)
445
446c  traitement des mailles
447               DO iju=1,niju
448                  ij=indu(iju)
449                  j=(ij-1)/iip1+1
450                  zu_m=u_m(ij,l)
451                  u_mq(ij,l)=0.
452                  IF(zu_m.gt.0.) THEN
453                     ijq=ij
454                     i=ijq-(j-1)*iip1
455c   accumulation pour les mailles completements advectees
456                     do while(zu_m.gt.masse(ijq,l))
457                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
458                        zu_m=zu_m-masse(ijq,l)
459                        i=mod(i-2+iim,iim)+1
460                        ijq=(j-1)*iip1+i
461                     ENDDO
462c   ajout de la maille non completement advectee
463                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
464     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
465                  ELSE
466                     ijq=ij+1
467                     i=ijq-(j-1)*iip1
468c   accumulation pour les mailles completements advectees
469                     do while(-zu_m.gt.masse(ijq,l))
470                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
471                        zu_m=zu_m+masse(ijq,l)
472                        i=mod(i,iim)+1
473                        ijq=(j-1)*iip1+i
474                     ENDDO
475c   ajout de la maille non completement advectee
476                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
477     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
478                  ENDIF
479               ENDDO
480            ENDIF
481         ENDDO
[764]482c$OMP END DO NOWAIT
483cym      ENDIF  ! n0.gt.0
[630]4849999    continue
485
486
487c   bouclage en latitude
488c       print*,'Avant bouclage en latitude'
[764]489c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]490      DO l=1,llm
491        DO ij=ijb+iip1-1,ije,iip1
492           u_mq(ij,l)=u_mq(ij-iim,l)
493        ENDDO
494      ENDDO
[764]495c$OMP END DO NOWAIT
[630]496
497c   calcul des tENDances
[764]498c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]499      DO l=1,llm
500         DO ij=ijb+1,ije
501            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
502            q(ij,l)=(q(ij,l)*masse(ij,l)+
503     &      u_mq(ij-1,l)-u_mq(ij,l))
504     &      /new_m
505            masse(ij,l)=new_m
506         ENDDO
507c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
508         DO ij=ijb+iip1-1,ije,iip1
509            q(ij-iim,l)=q(ij,l)
510            masse(ij-iim,l)=masse(ij,l)
511         ENDDO
512      ENDDO
[764]513c$OMP END DO NOWAIT
[630]514c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
515c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
516
517
518      RETURN
519      END
520
521
522      SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v)
523c
524c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
525c
526c    ********************************************************************
527c     Shema  d'advection " pseudo amont " .
528c    ********************************************************************
529c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
530c     dq               sont des arguments de sortie pour le s-pg ....
531c
532c
533c   --------------------------------------------------------------------
[1823]534      USE parallel_lmdz
[2597]535      USE comconst_mod, ONLY: pi
[630]536      IMPLICIT NONE
537c
538#include "dimensions.h"
539#include "paramet.h"
540#include "comgeom.h"
541c
542c
543c   Arguments:
544c   ----------
545      REAL masse(ip1jmp1,llm),pente_max
546      REAL masse_adv_v( ip1jm,llm)
547      REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
548c
549c      Local
550c   ---------
551c
552      INTEGER i,ij,l
553c
554      REAL airej2,airejjm,airescb(iim),airesch(iim)
555      REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
556      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
557      REAL qbyv(ip1jm,llm)
558
[1520]559      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
[630]560c     REAL newq,oldmasse
561      Logical extremum,first,testcpu
562      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
563      SAVE temps0,temps1,temps2,temps3,temps4,temps5
[764]564c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
[630]565      SAVE first,testcpu
[764]566c$OMP THREADPRIVATE(first,testcpu)
[630]567
568      REAL convpn,convps,convmpn,convmps
569      real massepn,masseps,qpn,qps
570      REAL sinlon(iip1),sinlondlon(iip1)
571      REAL coslon(iip1),coslondlon(iip1)
572      SAVE sinlon,coslon,sinlondlon,coslondlon
[764]573c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
[630]574      SAVE airej2,airejjm
[764]575c$OMP THREADPRIVATE(airej2,airejjm)
[630]576c
577c
578      REAL      SSUM
579      EXTERNAL  SSUM
580
581      DATA first,testcpu/.true.,.false./
582      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
583      INTEGER ijb,ije
584
585      IF(first) THEN
586c         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
587         first=.false.
588         do i=2,iip1
589            coslon(i)=cos(rlonv(i))
590            sinlon(i)=sin(rlonv(i))
591            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
592            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
593         ENDDO
594         coslon(1)=coslon(iip1)
595         coslondlon(1)=coslondlon(iip1)
596         sinlon(1)=sinlon(iip1)
597         sinlondlon(1)=sinlondlon(iip1)
598         airej2 = SSUM( iim, aire(iip2), 1 )
599         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
600      ENDIF
601
602c
603c       PRINT*,'CALCUL EN LATITUDE'
604
[764]605c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]606      DO l = 1, llm
607c
608c   --------------------------------
609c      CALCUL EN LATITUDE
610c   --------------------------------
611
612c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
613c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
614c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
615     
616      if (pole_nord) then
617        DO i = 1, iim
618          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
619        ENDDO
620        qpns   = SSUM( iim,  airescb ,1 ) / airej2
621      endif
622     
623      if (pole_sud) then
624        DO i = 1, iim
625          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
626        ENDDO
627        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
628      endif
629     
630     
631
632c   calcul des pentes aux points v
633
634      ijb=ij_begin-2*iip1
635      ije=ij_end+iip1
636      if (pole_nord) ijb=ij_begin
637      if (pole_sud)  ije=ij_end-iip1
638     
639      DO ij=ijb,ije
640         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
641         adyqv(ij)=abs(dyqv(ij))
642      ENDDO
643
644c   calcul des pentes aux points scalaires
645      ijb=ij_begin-iip1
646      ije=ij_end+iip1
647      if (pole_nord) ijb=ij_begin+iip1
648      if (pole_sud)  ije=ij_end-iip1
649     
650      DO ij=ijb,ije
651         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
652         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
653         dyqmax(ij)=pente_max*dyqmax(ij)
654      ENDDO
655
656c   calcul des pentes aux poles
657      IF (pole_nord) THEN
658        DO ij=1,iip1
659           dyq(ij,l)=qpns-q(ij+iip1,l)
660        ENDDO
661       
662        dyn1=0.
663        dyn2=0.
664        DO ij=1,iim
665          dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
666          dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
667        ENDDO
668        DO ij=1,iip1
669          dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
670        ENDDO
671       
672        DO ij=1,iip1
673         dyq(ij,l)=0.
674        ENDDO
675c ym tout cela ne sert pas a grand chose
676      ENDIF
677     
678      IF (pole_sud) THEN
679
680        DO ij=1,iip1
681           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
682        ENDDO
683
684        dys1=0.
685        dys2=0.
686
687        DO ij=1,iim
688          dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
689          dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
690        ENDDO
691
692        DO ij=1,iip1
693          dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
694        ENDDO
695       
696        DO ij=1,iip1
697         dyq(ip1jm+ij,l)=0.
698        ENDDO
699c ym tout cela ne sert pas a grand chose
700      ENDIF
701
702c   filtrage de la derivee
703
704c   calcul des pentes limites aux poles
705c ym partie inutile
706c      goto 8888
707c      fn=1.
708c      fs=1.
709c      DO ij=1,iim
710c         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
711c            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
712c         ENDIF
713c      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
714c         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
715c         ENDIF
716c      ENDDO
717c      DO ij=1,iip1
718c         dyq(ij,l)=fn*dyq(ij,l)
719c         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
720c      ENDDO
721c 8888    continue
722
723
724CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
725C  En memoire de dIFferents tests sur la
726C  limitation des pentes aux poles.
727CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
728C     PRINT*,dyq(1)
729C     PRINT*,dyqv(iip1+1)
[1520]730C     appn=abs(dyq(1)/dyqv(iip1+1))
[630]731C     PRINT*,dyq(ip1jm+1)
732C     PRINT*,dyqv(ip1jm-iip1+1)
[1520]733C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[630]734C     DO ij=2,iim
[1520]735C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
736C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
[630]737C     ENDDO
[1520]738C     appn=min(pente_max/appn,1.)
739C     apps=min(pente_max/apps,1.)
[630]740C
741C
742C   cas ou on a un extremum au pole
743C
744C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
[1520]745C    &   appn=0.
[630]746C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
747C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
[1520]748C    &   apps=0.
[630]749C
750C   limitation des pentes aux poles
751C     DO ij=1,iip1
[1520]752C        dyq(ij)=appn*dyq(ij)
753C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
[630]754C     ENDDO
755C
756C   test
757C      DO ij=1,iip1
758C         dyq(iip1+ij)=0.
759C         dyq(ip1jm+ij-iip1)=0.
760C      ENDDO
761C      DO ij=1,ip1jmp1
762C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
763C      ENDDO
764C
765C changement 10 07 96
766C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
767C    &   THEN
768C        DO ij=1,iip1
769C           dyqmax(ij)=0.
770C        ENDDO
771C     ELSE
772C        DO ij=1,iip1
773C           dyqmax(ij)=pente_max*abs(dyqv(ij))
774C        ENDDO
775C     ENDIF
776C
777C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
778C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
779C    &THEN
780C        DO ij=ip1jm+1,ip1jmp1
781C           dyqmax(ij)=0.
782C        ENDDO
783C     ELSE
784C        DO ij=ip1jm+1,ip1jmp1
785C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
786C        ENDDO
787C     ENDIF
788C   fin changement 10 07 96
789CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
790
791c   calcul des pentes limitees
792      ijb=ij_begin-iip1
793      ije=ij_end+iip1
794      if (pole_nord) ijb=ij_begin+iip1
795      if (pole_sud)  ije=ij_end-iip1
796
797      DO ij=ijb,ije
798         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
799            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
800         ELSE
801            dyq(ij,l)=0.
802         ENDIF
803      ENDDO
804
805      ENDDO
[764]806c$OMP END DO NOWAIT
[630]807
808      ijb=ij_begin-iip1
809      ije=ij_end
810      if (pole_nord) ijb=ij_begin
811      if (pole_sud)  ije=ij_end-iip1
812
[764]813c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]814      DO l=1,llm
815       DO ij=ijb,ije
816          IF(masse_adv_v(ij,l).gt.0) THEN
817              qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
818     ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
819          ELSE
820              qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
821     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
822          ENDIF
823          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
824       ENDDO
825      ENDDO
[764]826c$OMP END DO NOWAIT
[630]827     
828      ijb=ij_begin
829      ije=ij_end
830      if (pole_nord) ijb=ij_begin+iip1
831      if (pole_sud)  ije=ij_end-iip1
832     
[764]833c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]834      DO l=1,llm
835         DO ij=ijb,ije
836            newmasse=masse(ij,l)
837     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
838     
839            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
840     &         /newmasse
841            masse(ij,l)=newmasse
842         ENDDO
843c.-. ancienne version
844c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
845c        convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
846         if (pole_nord) then
847           convpn=SSUM(iim,qbyv(1,l),1)
848           convmpn=ssum(iim,masse_adv_v(1,l),1)
849           massepn=ssum(iim,masse(1,l),1)
850           qpn=0.
851           do ij=1,iim
852              qpn=qpn+masse(ij,l)*q(ij,l)
853           enddo
854           qpn=(qpn+convpn)/(massepn+convmpn)
855           do ij=1,iip1
856              q(ij,l)=qpn
857           enddo
858         endif
859         
860c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
861c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
862         
863         if (pole_sud) then
864         
865           convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
866           convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
867           masseps=ssum(iim, masse(ip1jm+1,l),1)
868           qps=0.
869           do ij = ip1jm+1,ip1jmp1-1
870              qps=qps+masse(ij,l)*q(ij,l)
871           enddo
872           qps=(qps+convps)/(masseps+convmps)
873           do ij=ip1jm+1,ip1jmp1
874              q(ij,l)=qps
875           enddo
876         endif
877c.-. fin ancienne version
878
879c._. nouvelle version
880c        convpn=SSUM(iim,qbyv(1,l),1)
881c        convmpn=ssum(iim,masse_adv_v(1,l),1)
882c        oldmasse=ssum(iim,masse(1,l),1)
883c        newmasse=oldmasse+convmpn
884c        newq=(q(1,l)*oldmasse+convpn)/newmasse
885c        newmasse=newmasse/apoln
886c        DO ij = 1,iip1
887c           q(ij,l)=newq
888c           masse(ij,l)=newmasse*aire(ij)
889c        ENDDO
890c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
891c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
892c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
893c        newmasse=oldmasse+convmps
894c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
895c        newmasse=newmasse/apols
896c        DO ij = ip1jm+1,ip1jmp1
897c           q(ij,l)=newq
898c           masse(ij,l)=newmasse*aire(ij)
899c        ENDDO
900c._. fin nouvelle version
901      ENDDO
[764]902c$OMP END DO NOWAIT
[630]903
904      RETURN
905      END
906     
907     
908     
909      SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x)
910c
911c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
912c
913c    ********************************************************************
914c     Shema  d'advection " pseudo amont " .
915c    ********************************************************************
916c    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
917c     dq               sont des arguments de sortie pour le s-pg ....
918c
919c
920c   --------------------------------------------------------------------
[1823]921      USE parallel_lmdz
[630]922      IMPLICIT NONE
923c
924#include "dimensions.h"
925#include "paramet.h"
926c
927c
928c   Arguments:
929c   ----------
930      REAL masse(ip1jmp1,llm),pente_max
931      REAL q(ip1jmp1,llm)
932      REAL w(ip1jmp1,llm+1)
933c
934c      Local
935c   ---------
936c
937      INTEGER i,ij,l,j,ii
938c
[764]939      REAL,SAVE :: wq(ip1jmp1,llm+1)
940      REAL newmasse
[630]941
[764]942      REAL,SAVE :: dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm)
943      REAL dzqmax
[630]944      REAL sigw
945
946      LOGICAL testcpu
947      SAVE testcpu
[764]948c$OMP THREADPRIVATE(testcpu)
[630]949      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
950      SAVE temps0,temps1,temps2,temps3,temps4,temps5
[764]951c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
952
[630]953      REAL      SSUM
[764]954      EXTERNAL  SSUM
[630]955
956      DATA testcpu/.false./
957      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
958      INTEGER ijb,ije,ijb_x,ije_x
959c    On oriente tout dans le sens de la pression c'est a dire dans le
960c    sens de W
961
962#ifdef BIDON
963      IF(testcpu) THEN
964         temps0=second(0.)
965      ENDIF
966#endif
967
968      ijb=ijb_x
969      ije=ije_x
[764]970
971c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[630]972      DO l=2,llm
973         DO ij=ijb,ije
974            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
975            adzqw(ij,l)=abs(dzqw(ij,l))
976         ENDDO
977      ENDDO
[764]978c$OMP END DO
[630]979
[764]980c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]981      DO l=2,llm-1
982         DO ij=ijb,ije
983#ifdef CRAY
984            dzq(ij,l)=0.5*
985     ,      cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
986#else
987            IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
988                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
989            ELSE
990                dzq(ij,l)=0.
991            ENDIF
992#endif
993            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
994            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
995         ENDDO
996      ENDDO
[764]997c$OMP END DO NOWAIT
[630]998
[764]999c$OMP MASTER
[630]1000      DO ij=ijb,ije
1001         dzq(ij,1)=0.
1002         dzq(ij,llm)=0.
1003      ENDDO
[764]1004c$OMP END MASTER
1005c$OMP BARRIER
[630]1006#ifdef BIDON
1007      IF(testcpu) THEN
1008         temps1=temps1+second(0.)-temps0
1009      ENDIF
1010#endif
1011c ---------------------------------------------------------------
1012c   .... calcul des termes d'advection verticale  .......
1013c ---------------------------------------------------------------
1014
1015c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
1016
[764]1017c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]1018       DO l = 1,llm-1
1019         do  ij = ijb,ije
1020          IF(w(ij,l+1).gt.0.) THEN
1021             sigw=w(ij,l+1)/masse(ij,l+1)
1022             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
1023          ELSE
1024             sigw=w(ij,l+1)/masse(ij,l)
1025             wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
1026          ENDIF
1027         ENDDO
1028       ENDDO
[764]1029c$OMP END DO NOWAIT
[630]1030
[764]1031c$OMP MASTER
[630]1032       DO ij=ijb,ije
1033          wq(ij,llm+1)=0.
1034          wq(ij,1)=0.
1035       ENDDO
[764]1036c$OMP END MASTER
1037c$OMP BARRIER
[630]1038
[764]1039c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]1040      DO l=1,llm
1041         DO ij=ijb,ije
1042            newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
1043            q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
1044     &         /newmasse
1045            masse(ij,l)=newmasse
1046         ENDDO
1047      ENDDO
[764]1048c$OMP END DO NOWAIT
[630]1049
1050
1051      RETURN
1052      END
1053c      SUBROUTINE minmaxq(zq,qmin,qmax,comment)
1054c
1055c#include "dimensions.h"
1056c#include "paramet.h"
1057
1058c      CHARACTER*(*) comment
1059c      real qmin,qmax
1060c      real zq(ip1jmp1,llm)
1061
1062c      INTEGER jadrs(ip1jmp1), jbad, k, i
1063
1064
1065c      DO k = 1, llm
1066c         jbad = 0
1067c         DO i = 1, ip1jmp1
1068c         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
1069c            jbad = jbad + 1
1070c            jadrs(jbad) = i
1071c         ENDIF
1072c         ENDDO
1073c         IF (jbad.GT.0) THEN
1074c         PRINT*, comment
1075c         DO i = 1, jbad
1076cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
1077c         ENDDO
1078c         ENDIF
1079c      ENDDO
1080
1081c      return
1082c      end
1083
1084
1085      subroutine minmaxq_p(zq,qmin,qmax,comment)
1086
1087#include "dimensions.h"
1088#include "paramet.h"
1089
1090      character*20 comment
1091      real qmin,qmax
1092      real zq(ip1jmp1,llm)
1093      real zzq(iip1,jjp1,llm)
1094
1095      integer imin,jmin,lmin,ijlmin
1096      integer imax,jmax,lmax,ijlmax
1097
1098      integer ismin,ismax
1099
1100#ifdef isminismax
1101      call scopy (ip1jmp1*llm,zq,1,zzq,1)
1102
1103      ijlmin=ismin(ijp1llm,zq,1)
1104      lmin=(ijlmin-1)/ip1jmp1+1
1105      ijlmin=ijlmin-(lmin-1.)*ip1jmp1
1106      jmin=(ijlmin-1)/iip1+1
1107      imin=ijlmin-(jmin-1.)*iip1
1108      zqmin=zq(ijlmin,lmin)
1109
1110      ijlmax=ismax(ijp1llm,zq,1)
1111      lmax=(ijlmax-1)/ip1jmp1+1
1112      ijlmax=ijlmax-(lmax-1.)*ip1jmp1
1113      jmax=(ijlmax-1)/iip1+1
1114      imax=ijlmax-(jmax-1.)*iip1
1115      zqmax=zq(ijlmax,lmax)
1116
1117       if(zqmin.lt.qmin)
1118c     s     write(*,9999) comment,
1119     s     write(*,*) comment,
1120     s     imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
1121       if(zqmax.gt.qmax)
1122c     s     write(*,9999) comment,
1123     s     write(*,*) comment,
1124     s     imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
1125#endif
1126      return
11279999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
1128      end
1129
1130
1131
Note: See TracBrowser for help on using the repository browser.