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

Last change on this file since 1992 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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