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

Last change on this file since 4625 was 2600, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.2 KB
Line 
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
33      include "dimensions.h"
34      include "paramet.h"
35      include "comgeom.h"
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
137      DO 500 l = 1,llm
138         DO 500 j = 1,jjp1
139            DO 500 i = 1,iip1 
140            wgri (i,j,llm+1-l) = w (i,j,l) 
141  500 CONTINUE
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
156      DO 1 K=1,LAT
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
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
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
182 100     CONTINUE
183 10   CONTINUE
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
195      DO 11 L=1,NIV-1
196      LP=L+1
197C
198      DO 110 I=1,LON
199C
200         IF(WGRI(I,K,L).LT.0.) THEN
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
217 110  CONTINUE
218C
219      DO 111 JV=1,NTRA
220      DO 1110 I=1,LON
221C
222         IF(WGRI(I,K,L).LT.0.) THEN
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
275 1110 CONTINUE
276 111  CONTINUE
277C
278 11   CONTINUE
279C
280C  puts the temporary moments Fi into appropriate neighboring boxes
281C
282      DO 12 L=1,NIV-1
283      LP=L+1
284C
285      DO 120 I=1,LON
286C
287         IF(WGRI(I,K,L).LT.0.) THEN
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
301 120  CONTINUE
302C
303      DO 121 JV=1,NTRA
304      DO 1210 I=1,LON
305C
306         IF(WGRI(I,K,L).LT.0.) THEN
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
344 1210 CONTINUE
345 121  CONTINUE
346C
347 12   CONTINUE
348C
349C  fin de la boucle principale sur les latitudes
350C
351 1    CONTINUE
352C
353      DO l = 1,llm
354      DO j = 1,jjp1
355          SM(iip1,j,l) = SM(1,j,l)
356          S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
357          SSX(iip1,j,l,ntra) = SSX(1,j,l,ntra)
358          SY(iip1,j,l,ntra) = SY(1,j,l,ntra)
359          SZ(iip1,j,l,ntra) = SZ(1,j,l,ntra)
360      ENDDO
361      ENDDO
362c                                                                                C-------------------------------------------------------------
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.