source: LMDZ.3.3/trunk/libf/dyn3d/vlsplt.F @ 346

Last change on this file since 346 was 2, checked in by lmdz, 25 years ago

Initial revision

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