source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advzp.F @ 5098

Last change on this file since 5098 was 5079, checked in by abarral, 2 months ago

[coherence with Fortran standards]
Replace obsolete DO with shared termination
(minor) replace obsolete bool operators

  • 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
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 l = 1,llm
138         DO j = 1,jjp1
139            DO i = 1,iip1 
140            wgri (i,j,llm+1-l) = w (i,j,l) 
141      END DO
142      END DO
143      END DO
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 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 JV=1,NTRA
166      DO L=1,NIV
167         DO I=1,LON
168            IF(S0(I,K,L,JV)>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      END DO
185      END DO
186      END DO
187C
188 101  CONTINUE
189C
190C  boucle sur les niveaux intercouches de 1 a NIV-1
191C   (flux nul au sommet L=0 et a la base L=NIV)
192C
193C  calculate flux and moments between adjacent boxes
194C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
195C  1- create temporary moments/masses for partial boxes in transit
196C  2- reajusts moments remaining in the box
197C
198      DO L=1,NIV-1
199      LP=L+1
200C
201      DO I=1,LON
202C
203         IF(WGRI(I,K,L)<0.) THEN
204           FM(I,L)=-WGRI(I,K,L)*DTZ
205           ALF(I)=FM(I,L)/SM(I,K,LP)
206           SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
207         ELSE
208           FM(I,L)=WGRI(I,K,L)*DTZ
209           ALF(I)=FM(I,L)/SM(I,K,L)
210           SM(I,K,L)=SM(I,K,L)-FM(I,L)
211         ENDIF
212C
213         ALFQ (I)=ALF(I)*ALF(I)
214         ALF1 (I)=1.-ALF(I)
215         ALF1Q(I)=ALF1(I)*ALF1(I)
216         ALF2 (I)=ALF1(I)-ALF(I)
217         ALF3 (I)=ALF(I)*ALFQ(I)
218         ALF4 (I)=ALF1(I)*ALF1Q(I)
219C
220      END DO
221C
222      DO JV=1,NTRA
223      DO I=1,LON
224C
225         IF(WGRI(I,K,L)<0.) THEN
226C
227           F0 (I,L,JV)=ALF (I)* ( S0(I,K,LP,JV)-ALF1(I)*
228     +          ( SZ(I,K,LP,JV)-ALF2(I)*SZZ(I,K,LP,JV) ) )
229           FZ (I,L,JV)=ALFQ(I)*(SZ(I,K,LP,JV)-3.*ALF1(I)*SZZ(I,K,LP,JV))
230           FZZ(I,L,JV)=ALF3(I)*SZZ(I,K,LP,JV)
231           FXZ(I,L,JV)=ALFQ(I)*SSXZ(I,K,LP,JV)
232           FYZ(I,L,JV)=ALFQ(I)*SYZ(I,K,LP,JV)
233           FX (I,L,JV)=ALF (I)*(SSX(I,K,LP,JV)-ALF1(I)*SSXZ(I,K,LP,JV))
234           FY (I,L,JV)=ALF (I)*(SY(I,K,LP,JV)-ALF1(I)*SYZ(I,K,LP,JV))
235           FXX(I,L,JV)=ALF (I)*SSXX(I,K,LP,JV)
236           FXY(I,L,JV)=ALF (I)*SSXY(I,K,LP,JV)
237           FYY(I,L,JV)=ALF (I)*SYY(I,K,LP,JV)
238C
239           S0 (I,K,LP,JV)=S0 (I,K,LP,JV)-F0 (I,L,JV)
240           SZ (I,K,LP,JV)=ALF1Q(I)
241     +                   *(SZ(I,K,LP,JV)+3.*ALF(I)*SZZ(I,K,LP,JV))
242           SZZ(I,K,LP,JV)=ALF4 (I)*SZZ(I,K,LP,JV)
243           SSXZ(I,K,LP,JV)=ALF1Q(I)*SSXZ(I,K,LP,JV)
244           SYZ(I,K,LP,JV)=ALF1Q(I)*SYZ(I,K,LP,JV)
245           SSX (I,K,LP,JV)=SSX (I,K,LP,JV)-FX (I,L,JV)
246           SY (I,K,LP,JV)=SY (I,K,LP,JV)-FY (I,L,JV)
247           SSXX(I,K,LP,JV)=SSXX(I,K,LP,JV)-FXX(I,L,JV)
248           SSXY(I,K,LP,JV)=SSXY(I,K,LP,JV)-FXY(I,L,JV)
249           SYY(I,K,LP,JV)=SYY(I,K,LP,JV)-FYY(I,L,JV)
250C
251         ELSE
252C
253           F0 (I,L,JV)=ALF (I)*(S0(I,K,L,JV)
254     +           +ALF1(I) * (SZ(I,K,L,JV)+ALF2(I)*SZZ(I,K,L,JV)) )
255           FZ (I,L,JV)=ALFQ(I)*(SZ(I,K,L,JV)+3.*ALF1(I)*SZZ(I,K,L,JV))
256           FZZ(I,L,JV)=ALF3(I)*SZZ(I,K,L,JV)
257           FXZ(I,L,JV)=ALFQ(I)*SSXZ(I,K,L,JV)
258           FYZ(I,L,JV)=ALFQ(I)*SYZ(I,K,L,JV)
259           FX (I,L,JV)=ALF (I)*(SSX(I,K,L,JV)+ALF1(I)*SSXZ(I,K,L,JV))
260           FY (I,L,JV)=ALF (I)*(SY(I,K,L,JV)+ALF1(I)*SYZ(I,K,L,JV))
261           FXX(I,L,JV)=ALF (I)*SSXX(I,K,L,JV)
262           FXY(I,L,JV)=ALF (I)*SSXY(I,K,L,JV)
263           FYY(I,L,JV)=ALF (I)*SYY(I,K,L,JV)
264C
265           S0 (I,K,L,JV)=S0 (I,K,L,JV)-F0(I,L,JV)
266           SZ (I,K,L,JV)=ALF1Q(I)*(SZ(I,K,L,JV)-3.*ALF(I)*SZZ(I,K,L,JV))
267           SZZ(I,K,L,JV)=ALF4 (I)*SZZ(I,K,L,JV)
268           SSXZ(I,K,L,JV)=ALF1Q(I)*SSXZ(I,K,L,JV)
269           SYZ(I,K,L,JV)=ALF1Q(I)*SYZ(I,K,L,JV)
270           SSX (I,K,L,JV)=SSX (I,K,L,JV)-FX (I,L,JV)
271           SY (I,K,L,JV)=SY (I,K,L,JV)-FY (I,L,JV)
272           SSXX(I,K,L,JV)=SSXX(I,K,L,JV)-FXX(I,L,JV)
273           SSXY(I,K,L,JV)=SSXY(I,K,L,JV)-FXY(I,L,JV)
274           SYY(I,K,L,JV)=SYY(I,K,L,JV)-FYY(I,L,JV)
275C
276         ENDIF
277C
278      END DO
279      END DO
280C
281      END DO
282C
283C  puts the temporary moments Fi into appropriate neighboring boxes
284C
285      DO L=1,NIV-1
286      LP=L+1
287C
288      DO I=1,LON
289C
290         IF(WGRI(I,K,L)<0.) THEN
291           SM(I,K,L)=SM(I,K,L)+FM(I,L)
292           ALF(I)=FM(I,L)/SM(I,K,L)
293         ELSE
294           SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
295           ALF(I)=FM(I,L)/SM(I,K,LP)
296         ENDIF
297C
298         ALF1(I)=1.-ALF(I)
299         ALFQ(I)=ALF(I)*ALF(I)
300         ALF1Q(I)=ALF1(I)*ALF1(I)
301         ALF2(I)=ALF(I)*ALF1(I)
302         ALF3(I)=ALF1(I)-ALF(I)
303C
304      END DO
305C
306      DO JV=1,NTRA
307      DO I=1,LON
308C
309         IF(WGRI(I,K,L)<0.) THEN
310C
311           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
312           S0 (I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
313           SZZ(I,K,L,JV)=ALFQ(I)*FZZ(I,L,JV)+ALF1Q(I)*SZZ(I,K,L,JV)
314     +        +5.*( ALF2(I)*(FZ(I,L,JV)-SZ(I,K,L,JV))+ALF3(I)*TEMPTM )
315           SZ (I,K,L,JV)=ALF (I)*FZ (I,L,JV)+ALF1 (I)*SZ (I,K,L,JV)
316     +                  +3.*TEMPTM
317           SSXZ(I,K,L,JV)=ALF (I)*FXZ(I,L,JV)+ALF1 (I)*SSXZ(I,K,L,JV)
318     +              +3.*(ALF1(I)*FX (I,L,JV)-ALF  (I)*SSX (I,K,L,JV))
319           SYZ(I,K,L,JV)=ALF (I)*FYZ(I,L,JV)+ALF1 (I)*SYZ(I,K,L,JV)
320     +              +3.*(ALF1(I)*FY (I,L,JV)-ALF  (I)*SY (I,K,L,JV))
321           SSX (I,K,L,JV)=SSX (I,K,L,JV)+FX (I,L,JV)
322           SY (I,K,L,JV)=SY (I,K,L,JV)+FY (I,L,JV)
323           SSXX(I,K,L,JV)=SSXX(I,K,L,JV)+FXX(I,L,JV)
324           SSXY(I,K,L,JV)=SSXY(I,K,L,JV)+FXY(I,L,JV)
325           SYY(I,K,L,JV)=SYY(I,K,L,JV)+FYY(I,L,JV)
326C
327         ELSE
328C
329           TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)
330           S0 (I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
331           SZZ(I,K,LP,JV)=ALFQ(I)*FZZ(I,L,JV)+ALF1Q(I)*SZZ(I,K,LP,JV)
332     +        +5.*( ALF2(I)*(SZ(I,K,LP,JV)-FZ(I,L,JV))-ALF3(I)*TEMPTM )
333           SZ (I,K,LP,JV)=ALF (I)*FZ(I,L,JV)+ALF1(I)*SZ(I,K,LP,JV)
334     +                   +3.*TEMPTM
335           SSXZ(I,K,LP,JV)=ALF(I)*FXZ(I,L,JV)+ALF1(I)*SSXZ(I,K,LP,JV)
336     +                   +3.*(ALF(I)*SSX(I,K,LP,JV)-ALF1(I)*FX(I,L,JV))
337           SYZ(I,K,LP,JV)=ALF(I)*FYZ(I,L,JV)+ALF1(I)*SYZ(I,K,LP,JV)
338     +                   +3.*(ALF(I)*SY(I,K,LP,JV)-ALF1(I)*FY(I,L,JV))
339           SSX (I,K,LP,JV)=SSX (I,K,LP,JV)+FX (I,L,JV)
340           SY (I,K,LP,JV)=SY (I,K,LP,JV)+FY (I,L,JV)
341           SSXX(I,K,LP,JV)=SSXX(I,K,LP,JV)+FXX(I,L,JV)
342           SSXY(I,K,LP,JV)=SSXY(I,K,LP,JV)+FXY(I,L,JV)
343           SYY(I,K,LP,JV)=SYY(I,K,LP,JV)+FYY(I,L,JV)
344C
345         ENDIF
346C
347      END DO
348      END DO
349C
350      END DO
351C
352C  fin de la boucle principale sur les latitudes
353C
354      END DO
355C
356      DO l = 1,llm
357      DO j = 1,jjp1
358          SM(iip1,j,l) = SM(1,j,l)
359          S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
360          SSX(iip1,j,l,ntra) = SSX(1,j,l,ntra)
361          SY(iip1,j,l,ntra) = SY(1,j,l,ntra)
362          SZ(iip1,j,l,ntra) = SZ(1,j,l,ntra)
363      ENDDO
364      ENDDO
365c                                                                                C-------------------------------------------------------------
366C *** Test : diag de la qqtite totale de tarceur
367C            dans l'atmosphere avant l'advection en z
368       DO l = 1,llm
369       DO j = 1,jjp1
370       DO i = 1,iim
371          sqf = sqf + S0(i,j,l,ntra)
372       ENDDO
373       ENDDO
374       ENDDO
375       PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'
376       PRINT*,'sqf=', sqf
377
378      RETURN
379      END
Note: See TracBrowser for help on using the repository browser.