source: LMDZ6/trunk/libf/dyn3d_common/advzp.F @ 5238

Last change on this file since 5238 was 5084, checked in by Laurent Fairhead, 3 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF
  • 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: 12.2 KB
RevLine 
[524]1!
2! $Header$
3!
4      SUBROUTINE ADVZP(LIMIT,DTZ,W,SM,S0,SSX,SY,SZ
5     .                 ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra )
6
7      IMPLICIT NONE
8
9CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
10C                                                                 C
11C  second-order moments (SOM) advection of tracer in Z direction  C
12C                                                                 C
13CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
14C                                                                 C
15C  Source : Pascal Simon ( Meteo, CNRM )                          C
16C  Adaptation : A.A. (LGGE)                                       C
17C  Derniere Modif : 19/11/95 LAST                                 C
18C                                                                 C
19C  sont les arguments d'entree pour le s-pg                       C
20C                                                                 C
21C  argument de sortie du s-pg                                     C
22C                                                                 C
23CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
24CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
25C
26C Rem : Probleme aux poles il faut reecrire ce cas specifique
27C        Attention au sens de l'indexation
28C
29
30C
31C  parametres principaux du modele
32C
[2600]33      include "dimensions.h"
34      include "paramet.h"
35      include "comgeom.h"
[524]36C
37C  Arguments :
38C  ----------
39C  dty : frequence fictive d'appel du transport
40C  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
41c
42        INTEGER lon,lat,niv
43        INTEGER i,j,jv,k,kp,l,lp
44        INTEGER ntra
45c        PARAMETER (ntra = 1)
46c
47        REAL dtz
48        REAL w ( iip1,jjp1,llm )
49c
50C  moments: SM  total mass in each grid box
51C           S0  mass of tracer in each grid box
52C           Si  1rst order moment in i direction
53C
54      REAL SM(iip1,jjp1,llm)
55     +    ,S0(iip1,jjp1,llm,ntra)
56      REAL SSX(iip1,jjp1,llm,ntra)
57     +    ,SY(iip1,jjp1,llm,ntra)
58     +    ,SZ(iip1,jjp1,llm,ntra)
59     +    ,SSXX(iip1,jjp1,llm,ntra)
60     +    ,SSXY(iip1,jjp1,llm,ntra)
61     +    ,SSXZ(iip1,jjp1,llm,ntra)
62     +    ,SYY(iip1,jjp1,llm,ntra)
63     +    ,SYZ(iip1,jjp1,llm,ntra)
64     +    ,SZZ(iip1,jjp1,llm,ntra)
65C
66C  Local :
67C  -------
68C
69C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
70C  mass fluxes in kg
71C  declaration :
72C
73      REAL WGRI(iip1,jjp1,0:llm)
74
75C Rem : UGRI et VGRI ne sont pas utilises dans
76C  cette subroutine ( advection en z uniquement )
77C  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
78C         attention a celui de WGRI
79C
80C  the moments F are similarly defined and used as temporary
81C  storage for portions of the grid boxes in transit
82C
83C  the moments Fij are used as temporary storage for
84C  portions of the grid boxes in transit at the current level
85C
86C  work arrays
87C
88C
89      REAL F0(iim,llm,ntra),FM(iim,llm)
90      REAL FX(iim,llm,ntra),FY(iim,llm,ntra)
91      REAL FZ(iim,llm,ntra)
92      REAL FXX(iim,llm,ntra),FXY(iim,llm,ntra)
93      REAL FXZ(iim,llm,ntra),FYY(iim,llm,ntra)
94      REAL FYZ(iim,llm,ntra),FZZ(iim,llm,ntra)
95      REAL S00(ntra)
96      REAL SM0             ! Just temporal variable
97C
98C  work arrays
99C
100      REAL ALF(iim),ALF1(iim)
101      REAL ALFQ(iim),ALF1Q(iim)
102      REAL ALF2(iim),ALF3(iim)
103      REAL ALF4(iim)
104      REAL TEMPTM          ! Just temporal variable
105      REAL SLPMAX,S1MAX,S1NEW,S2NEW
106c
107      REAL sqi,sqf
108      LOGICAL LIMIT
109
110      lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
111      lat = jjp1        ! a cause des dim. differentes entre les
112      niv = llm         !       tab. S et VGRI
113                   
114c-----------------------------------------------------------------
115C *** Test : diag de la qtite totale de traceur dans
116C            l'atmosphere avant l'advection en Y
117c
118      sqi = 0.
119      sqf = 0.
120c
121      DO l = 1,llm
122         DO j = 1,jjp1
123           DO i = 1,iim
124              sqi = sqi + S0(i,j,l,ntra)
125           END DO
126         END DO
127      END DO
128      PRINT*,'---------- DIAG DANS ADVZP - ENTREE --------'
129      PRINT*,'sqi=',sqi
130
131c-----------------------------------------------------------------
132C  Interface : adaptation nouveau modele
133C  -------------------------------------
134C
135C  Conversion des flux de masses en kg
136
[5084]137      DO 500 l = 1,llm
138         DO 500 j = 1,jjp1
139            DO 500 i = 1,iip1 
[524]140            wgri (i,j,llm+1-l) = w (i,j,l) 
[5084]141  500 CONTINUE
[524]142      do j=1,jjp1
143         do i=1,iip1
144            wgri(i,j,0)=0.
145         enddo
146      enddo
147c
148cAA rem : Je ne suis pas sur du signe 
149cAA       Je ne suis pas sur pour le 0:llm
150c
151c-----------------------------------------------------------------
152C---------------------- START HERE -------------------------------
153C
154C  boucle sur les latitudes
155C
[5084]156      DO 1 K=1,LAT
[524]157C
158C  place limits on appropriate moments before transport
159C      (if flux-limiting is to be applied)
160C
161      IF(.NOT.LIMIT) GO TO 101
162C
[5084]163      DO 10 JV=1,NTRA
164      DO 10 L=1,NIV
165         DO 100 I=1,LON
166            IF(S0(I,K,L,JV).GT.0.) THEN
[524]167              SLPMAX=S0(I,K,L,JV)
168              S1MAX =1.5*SLPMAX
169              S1NEW =AMIN1(S1MAX,AMAX1(-S1MAX,SZ(I,K,L,JV)))
170              S2NEW =AMIN1( 2.*SLPMAX-ABS(S1NEW)/3. ,
171     +                     AMAX1(ABS(S1NEW)-SLPMAX,SZZ(I,K,L,JV)) )
172              SZ (I,K,L,JV)=S1NEW
173              SZZ(I,K,L,JV)=S2NEW
174              SSXZ(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SSXZ(I,K,L,JV)))
175              SYZ(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SYZ(I,K,L,JV)))
176            ELSE
177              SZ (I,K,L,JV)=0.
178              SZZ(I,K,L,JV)=0.
179              SSXZ(I,K,L,JV)=0.
180              SYZ(I,K,L,JV)=0.
181            ENDIF
[5084]182 100     CONTINUE
183 10   CONTINUE
[524]184C
185 101  CONTINUE
186C
187C  boucle sur les niveaux intercouches de 1 a NIV-1
188C   (flux nul au sommet L=0 et a la base L=NIV)
189C
190C  calculate flux and moments between adjacent boxes
191C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
192C  1- create temporary moments/masses for partial boxes in transit
193C  2- reajusts moments remaining in the box
194C
[5084]195      DO 11 L=1,NIV-1
[524]196      LP=L+1
197C
[5084]198      DO 110 I=1,LON
[524]199C
[5084]200         IF(WGRI(I,K,L).LT.0.) THEN
[524]201           FM(I,L)=-WGRI(I,K,L)*DTZ
202           ALF(I)=FM(I,L)/SM(I,K,LP)
203           SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
204         ELSE
205           FM(I,L)=WGRI(I,K,L)*DTZ
206           ALF(I)=FM(I,L)/SM(I,K,L)
207           SM(I,K,L)=SM(I,K,L)-FM(I,L)
208         ENDIF
209C
210         ALFQ (I)=ALF(I)*ALF(I)
211         ALF1 (I)=1.-ALF(I)
212         ALF1Q(I)=ALF1(I)*ALF1(I)
213         ALF2 (I)=ALF1(I)-ALF(I)
214         ALF3 (I)=ALF(I)*ALFQ(I)
215         ALF4 (I)=ALF1(I)*ALF1Q(I)
216C
[5084]217 110  CONTINUE
[524]218C
[5084]219      DO 111 JV=1,NTRA
220      DO 1110 I=1,LON
[524]221C
[5084]222         IF(WGRI(I,K,L).LT.0.) THEN
[524]223C
224           F0 (I,L,JV)=ALF (I)* ( S0(I,K,LP,JV)-ALF1(I)*
225     +          ( SZ(I,K,LP,JV)-ALF2(I)*SZZ(I,K,LP,JV) ) )
226           FZ (I,L,JV)=ALFQ(I)*(SZ(I,K,LP,JV)-3.*ALF1(I)*SZZ(I,K,LP,JV))
227           FZZ(I,L,JV)=ALF3(I)*SZZ(I,K,LP,JV)
228           FXZ(I,L,JV)=ALFQ(I)*SSXZ(I,K,LP,JV)
229           FYZ(I,L,JV)=ALFQ(I)*SYZ(I,K,LP,JV)
230           FX (I,L,JV)=ALF (I)*(SSX(I,K,LP,JV)-ALF1(I)*SSXZ(I,K,LP,JV))
231           FY (I,L,JV)=ALF (I)*(SY(I,K,LP,JV)-ALF1(I)*SYZ(I,K,LP,JV))
232           FXX(I,L,JV)=ALF (I)*SSXX(I,K,LP,JV)
233           FXY(I,L,JV)=ALF (I)*SSXY(I,K,LP,JV)
234           FYY(I,L,JV)=ALF (I)*SYY(I,K,LP,JV)
235C
236           S0 (I,K,LP,JV)=S0 (I,K,LP,JV)-F0 (I,L,JV)
237           SZ (I,K,LP,JV)=ALF1Q(I)
238     +                   *(SZ(I,K,LP,JV)+3.*ALF(I)*SZZ(I,K,LP,JV))
239           SZZ(I,K,LP,JV)=ALF4 (I)*SZZ(I,K,LP,JV)
240           SSXZ(I,K,LP,JV)=ALF1Q(I)*SSXZ(I,K,LP,JV)
241           SYZ(I,K,LP,JV)=ALF1Q(I)*SYZ(I,K,LP,JV)
242           SSX (I,K,LP,JV)=SSX (I,K,LP,JV)-FX (I,L,JV)
243           SY (I,K,LP,JV)=SY (I,K,LP,JV)-FY (I,L,JV)
244           SSXX(I,K,LP,JV)=SSXX(I,K,LP,JV)-FXX(I,L,JV)
245           SSXY(I,K,LP,JV)=SSXY(I,K,LP,JV)-FXY(I,L,JV)
246           SYY(I,K,LP,JV)=SYY(I,K,LP,JV)-FYY(I,L,JV)
247C
248         ELSE
249C
250           F0 (I,L,JV)=ALF (I)*(S0(I,K,L,JV)
251     +           +ALF1(I) * (SZ(I,K,L,JV)+ALF2(I)*SZZ(I,K,L,JV)) )
252           FZ (I,L,JV)=ALFQ(I)*(SZ(I,K,L,JV)+3.*ALF1(I)*SZZ(I,K,L,JV))
253           FZZ(I,L,JV)=ALF3(I)*SZZ(I,K,L,JV)
254           FXZ(I,L,JV)=ALFQ(I)*SSXZ(I,K,L,JV)
255           FYZ(I,L,JV)=ALFQ(I)*SYZ(I,K,L,JV)
256           FX (I,L,JV)=ALF (I)*(SSX(I,K,L,JV)+ALF1(I)*SSXZ(I,K,L,JV))
257           FY (I,L,JV)=ALF (I)*(SY(I,K,L,JV)+ALF1(I)*SYZ(I,K,L,JV))
258           FXX(I,L,JV)=ALF (I)*SSXX(I,K,L,JV)
259           FXY(I,L,JV)=ALF (I)*SSXY(I,K,L,JV)
260           FYY(I,L,JV)=ALF (I)*SYY(I,K,L,JV)
261C
262           S0 (I,K,L,JV)=S0 (I,K,L,JV)-F0(I,L,JV)
263           SZ (I,K,L,JV)=ALF1Q(I)*(SZ(I,K,L,JV)-3.*ALF(I)*SZZ(I,K,L,JV))
264           SZZ(I,K,L,JV)=ALF4 (I)*SZZ(I,K,L,JV)
265           SSXZ(I,K,L,JV)=ALF1Q(I)*SSXZ(I,K,L,JV)
266           SYZ(I,K,L,JV)=ALF1Q(I)*SYZ(I,K,L,JV)
267           SSX (I,K,L,JV)=SSX (I,K,L,JV)-FX (I,L,JV)
268           SY (I,K,L,JV)=SY (I,K,L,JV)-FY (I,L,JV)
269           SSXX(I,K,L,JV)=SSXX(I,K,L,JV)-FXX(I,L,JV)
270           SSXY(I,K,L,JV)=SSXY(I,K,L,JV)-FXY(I,L,JV)
271           SYY(I,K,L,JV)=SYY(I,K,L,JV)-FYY(I,L,JV)
272C
273         ENDIF
274C
[5084]275 1110 CONTINUE
276 111  CONTINUE
[524]277C
[5084]278 11   CONTINUE
[524]279C
280C  puts the temporary moments Fi into appropriate neighboring boxes
281C
[5084]282      DO 12 L=1,NIV-1
[524]283      LP=L+1
284C
[5084]285      DO 120 I=1,LON
[524]286C
[5084]287         IF(WGRI(I,K,L).LT.0.) THEN
[524]288           SM(I,K,L)=SM(I,K,L)+FM(I,L)
289           ALF(I)=FM(I,L)/SM(I,K,L)
290         ELSE
291           SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
292           ALF(I)=FM(I,L)/SM(I,K,LP)
293         ENDIF
294C
295         ALF1(I)=1.-ALF(I)
296         ALFQ(I)=ALF(I)*ALF(I)
297         ALF1Q(I)=ALF1(I)*ALF1(I)
298         ALF2(I)=ALF(I)*ALF1(I)
299         ALF3(I)=ALF1(I)-ALF(I)
300C
[5084]301 120  CONTINUE
[524]302C
[5084]303      DO 121 JV=1,NTRA
304      DO 1210 I=1,LON
[524]305C
[5084]306         IF(WGRI(I,K,L).LT.0.) THEN
[524]307C
308           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
309           S0 (I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
310           SZZ(I,K,L,JV)=ALFQ(I)*FZZ(I,L,JV)+ALF1Q(I)*SZZ(I,K,L,JV)
311     +        +5.*( ALF2(I)*(FZ(I,L,JV)-SZ(I,K,L,JV))+ALF3(I)*TEMPTM )
312           SZ (I,K,L,JV)=ALF (I)*FZ (I,L,JV)+ALF1 (I)*SZ (I,K,L,JV)
313     +                  +3.*TEMPTM
314           SSXZ(I,K,L,JV)=ALF (I)*FXZ(I,L,JV)+ALF1 (I)*SSXZ(I,K,L,JV)
315     +              +3.*(ALF1(I)*FX (I,L,JV)-ALF  (I)*SSX (I,K,L,JV))
316           SYZ(I,K,L,JV)=ALF (I)*FYZ(I,L,JV)+ALF1 (I)*SYZ(I,K,L,JV)
317     +              +3.*(ALF1(I)*FY (I,L,JV)-ALF  (I)*SY (I,K,L,JV))
318           SSX (I,K,L,JV)=SSX (I,K,L,JV)+FX (I,L,JV)
319           SY (I,K,L,JV)=SY (I,K,L,JV)+FY (I,L,JV)
320           SSXX(I,K,L,JV)=SSXX(I,K,L,JV)+FXX(I,L,JV)
321           SSXY(I,K,L,JV)=SSXY(I,K,L,JV)+FXY(I,L,JV)
322           SYY(I,K,L,JV)=SYY(I,K,L,JV)+FYY(I,L,JV)
323C
324         ELSE
325C
326           TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)
327           S0 (I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
328           SZZ(I,K,LP,JV)=ALFQ(I)*FZZ(I,L,JV)+ALF1Q(I)*SZZ(I,K,LP,JV)
329     +        +5.*( ALF2(I)*(SZ(I,K,LP,JV)-FZ(I,L,JV))-ALF3(I)*TEMPTM )
330           SZ (I,K,LP,JV)=ALF (I)*FZ(I,L,JV)+ALF1(I)*SZ(I,K,LP,JV)
331     +                   +3.*TEMPTM
332           SSXZ(I,K,LP,JV)=ALF(I)*FXZ(I,L,JV)+ALF1(I)*SSXZ(I,K,LP,JV)
333     +                   +3.*(ALF(I)*SSX(I,K,LP,JV)-ALF1(I)*FX(I,L,JV))
334           SYZ(I,K,LP,JV)=ALF(I)*FYZ(I,L,JV)+ALF1(I)*SYZ(I,K,LP,JV)
335     +                   +3.*(ALF(I)*SY(I,K,LP,JV)-ALF1(I)*FY(I,L,JV))
336           SSX (I,K,LP,JV)=SSX (I,K,LP,JV)+FX (I,L,JV)
337           SY (I,K,LP,JV)=SY (I,K,LP,JV)+FY (I,L,JV)
338           SSXX(I,K,LP,JV)=SSXX(I,K,LP,JV)+FXX(I,L,JV)
339           SSXY(I,K,LP,JV)=SSXY(I,K,LP,JV)+FXY(I,L,JV)
340           SYY(I,K,LP,JV)=SYY(I,K,LP,JV)+FYY(I,L,JV)
341C
342         ENDIF
343C
[5084]344 1210 CONTINUE
345 121  CONTINUE
[524]346C
[5084]347 12   CONTINUE
[524]348C
349C  fin de la boucle principale sur les latitudes
350C
[5084]351 1    CONTINUE
[524]352C
353      DO l = 1,llm
354      DO j = 1,jjp1
355          SM(iip1,j,l) = SM(1,j,l)
[2600]356          S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
[524]357          SSX(iip1,j,l,ntra) = SSX(1,j,l,ntra)
[2600]358          SY(iip1,j,l,ntra) = SY(1,j,l,ntra)
[524]359          SZ(iip1,j,l,ntra) = SZ(1,j,l,ntra)
360      ENDDO
361      ENDDO
[2600]362c                                                                                C-------------------------------------------------------------
[524]363C *** Test : diag de la qqtite totale de tarceur
364C            dans l'atmosphere avant l'advection en z
365       DO l = 1,llm
366       DO j = 1,jjp1
367       DO i = 1,iim
368          sqf = sqf + S0(i,j,l,ntra)
369       ENDDO
370       ENDDO
371       ENDDO
372       PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'
373       PRINT*,'sqf=', sqf
374
375      RETURN
376      END
Note: See TracBrowser for help on using the repository browser.