source: LMDZ5/trunk/libf/dyn3d_common/advzp.F @ 1985

Last change on this file since 1985 was 1945, checked in by lguez, 11 years ago

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