source: LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F @ 2076

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