trunk/dyn3d/advz.f

!
! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advz.F,v 1.2 2005/05/25 13:10:09 fairhead Exp $
!
      SUBROUTINE advz(limit,dtz,w,sm,s0,sx,sy,sz)
      use dimens_m
      use paramet_m
      use comconst
      use comvert
      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                C
C  first-order moments (FOM) advection of tracer in Z direction  C
C                                                                C
C  Source : Pascal Simon (Meteo,CNRM)                            C
C  Adaptation : A.Armengaud (LGGE) juin 94                       C
C                                                                C
C                                                                C
C  sont des arguments d'entree pour le s-pg...                   C
C                                                                C
C  dq est l'argument de sortie pour le s-pg                      C
C                                                                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  parametres principaux du modele
C

C  Arguments :
C  -----------
C  dtz : frequence fictive d'appel du transport 
C  w : flux de masse en z en Pa.m2.s-1

      INTEGER ntra
      PARAMETER (ntra = 1)

      REAL, intent(in):: dtz
      REAL w ( iip1,jjp1,llm )
    
C  moments: SM  total mass in each grid box
C           S0  mass of tracer in each grid box
C           Si  1rst order moment in i direction
C
      REAL SM(iip1,jjp1,llm)
     +    ,S0(iip1,jjp1,llm,ntra)
      REAL sx(iip1,jjp1,llm,ntra)
     +    ,sy(iip1,jjp1,llm,ntra)
     +    ,sz(iip1,jjp1,llm,ntra)


C  Local :
C  ------- 

C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
C  mass fluxes in kg
C  declaration :

      REAL WGRI(iip1,jjp1,0:llm)

C
C  the moments F are used as temporary  storage for 
C  portions of grid boxes in transit at the current latitude
C
      REAL FM(iim,llm)
      REAL F0(iim,llm,ntra),FX(iim,llm,ntra)
      REAL FY(iim,llm,ntra),FZ(iim,llm,ntra)
C
C  work arrays
C
      REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
      REAL TEMPTM            ! Just temporal variable
      REAL sqi,sqf
C
      LOGICAL LIMIT
      INTEGER lon,lat,niv
      INTEGER i,j,jv,k,l,lp

      lon = iim
      lat = jjp1
      niv = llm 

C *** Test : diag de la qqtite totale de traceur 
C            dans l'atmosphere avant l'advection en z
      sqi = 0.
      sqf = 0.

      DO l = 1,llm
         DO j = 1,jjp1
            DO i = 1,iim
cIM 240305            sqi = sqi + S0(i,j,l,9)
               sqi = sqi + S0(i,j,l,ntra)
            ENDDO
         ENDDO
      ENDDO
      PRINT*,'-------- DIAG DANS ADVZ - ENTREE ---------'
      PRINT*,'sqi=',sqi

C-----------------------------------------------------------------
C  Interface : adaptation nouveau modele
C  -------------------------------------
C
C  Conversion du flux de masse en kg.s-1

      DO 500 l = 1,llm
         DO 500 j = 1,jjp1
            DO 500 i = 1,iip1  
c            wgri (i,j,llm+1-l) =  w (i,j,l) / g 
               wgri (i,j,llm+1-l) =  w (i,j,l) 
c             wgri (i,j,0) = 0.                ! a detruire ult.
c             wgri (i,j,l) = 0.1               !    w (i,j,l) 
c             wgri (i,j,llm) = 0.              ! a detruire ult.
  500 CONTINUE
         DO  j = 1,jjp1
            DO i = 1,iip1  
               wgri(i,j,0)=0.
            enddo
         enddo

C-----------------------------------------------------------------
  
C  start here          
C  boucle sur les latitudes
C
      DO 1 K=1,LAT
C
C  place limits on appropriate moments before transport
C      (if flux-limiting is to be applied)
C
      IF(.NOT.LIMIT) GO TO 101
C
      DO 10 JV=1,NTRA
      DO 10 L=1,NIV
         DO 100 I=1,LON
            sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
     +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
 100     CONTINUE
 10   CONTINUE
C
 101  CONTINUE
C
C  boucle sur les niveaux intercouches de 1 a NIV-1
C   (flux nul au sommet L=0 et a la base L=NIV)
C
C  calculate flux and moments between adjacent boxes
C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
C  1- create temporary moments/masses for partial boxes in transit
C  2- reajusts moments remaining in the box
C
      DO 11 L=1,NIV-1
      LP=L+1
C
      DO 110 I=1,LON
C
         IF(WGRI(I,K,L).LT.0.) THEN
           FM(I,L)=-WGRI(I,K,L)*DTZ
           ALF(I)=FM(I,L)/SM(I,K,LP)
           SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
         ELSE
           FM(I,L)=WGRI(I,K,L)*DTZ
           ALF(I)=FM(I,L)/SM(I,K,L)
           SM(I,K,L)=SM(I,K,L)-FM(I,L)
         ENDIF
C
         ALFQ (I)=ALF(I)*ALF(I)
         ALF1 (I)=1.-ALF(I)
         ALF1Q(I)=ALF1(I)*ALF1(I)
C
 110  CONTINUE
C
      DO 111 JV=1,NTRA
      DO 1110 I=1,LON
C
         IF(WGRI(I,K,L).LT.0.) THEN
C
           F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )
           FZ(I,L,JV)=ALFQ(I)*sz(I,K,LP,JV)
           FX(I,L,JV)=ALF (I)*sx(I,K,LP,JV)
           FY(I,L,JV)=ALF (I)*sy(I,K,LP,JV)
C
           S0(I,K,LP,JV)=S0(I,K,LP,JV)-F0(I,L,JV)
           sz(I,K,LP,JV)=ALF1Q(I)*sz(I,K,LP,JV)
           sx(I,K,LP,JV)=sx(I,K,LP,JV)-FX(I,L,JV)
           sy(I,K,LP,JV)=sy(I,K,LP,JV)-FY(I,L,JV)
C
         ELSE
C
           F0(I,L,JV)=ALF (I)*(S0(I,K,L,JV)+ALF1(I)*sz(I,K,L,JV) )
           FZ(I,L,JV)=ALFQ(I)*sz(I,K,L,JV)
           FX(I,L,JV)=ALF (I)*sx(I,K,L,JV)
           FY(I,L,JV)=ALF (I)*sy(I,K,L,JV)
C
           S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,L,JV)
           sz(I,K,L,JV)=ALF1Q(I)*sz(I,K,L,JV)
           sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,L,JV)
           sy(I,K,L,JV)=sy(I,K,L,JV)-FY(I,L,JV)
C
         ENDIF
C
 1110 CONTINUE
 111  CONTINUE
C
 11   CONTINUE
C
C  puts the temporary moments Fi into appropriate neighboring boxes
C
      DO 12 L=1,NIV-1
      LP=L+1
C
      DO 120 I=1,LON
C
         IF(WGRI(I,K,L).LT.0.) THEN
           SM(I,K,L)=SM(I,K,L)+FM(I,L)
           ALF(I)=FM(I,L)/SM(I,K,L)
         ELSE
           SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
           ALF(I)=FM(I,L)/SM(I,K,LP)
         ENDIF
C
         ALF1(I)=1.-ALF(I)
         ALFQ(I)=ALF(I)*ALF(I)
         ALF1Q(I)=ALF1(I)*ALF1(I)
C
 120  CONTINUE
C
      DO 121 JV=1,NTRA
      DO 1210 I=1,LON
C
         IF(WGRI(I,K,L).LT.0.) THEN
C
           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
           S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
           sz(I,K,L,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,L,JV)+3.*TEMPTM
           sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,L,JV)
           sy(I,K,L,JV)=sy(I,K,L,JV)+FY(I,L,JV)
C
         ELSE
C
           TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)
           S0(I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
           sz(I,K,LP,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,LP,JV)
     +                  +3.*TEMPTM
           sx(I,K,LP,JV)=sx(I,K,LP,JV)+FX(I,L,JV)
           sy(I,K,LP,JV)=sy(I,K,LP,JV)+FY(I,L,JV)
C
         ENDIF
C
 1210 CONTINUE
 121  CONTINUE
C
 12   CONTINUE
C
C  fin de la boucle principale sur les latitudes
C
 1    CONTINUE

C *** ------------------- bouclage cyclique  en X ------------
      
c      DO l = 1,llm
c         DO j = 1,jjp1
c            SM(iip1,j,l) = SM(1,j,l)
c            S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
C            sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
c            sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
c            sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
c         ENDDO
c      ENDDO
           
C-------------------------------------------------------------
C *** Test : diag de la qqtite totale de traceur 
C            dans l'atmosphere avant l'advection en z
      DO l = 1,llm
         DO j = 1,jjp1
            DO i = 1,iim
cIM 240305            sqf = sqf + S0(i,j,l,9)
               sqf = sqf + S0(i,j,l,ntra)
            ENDDO
         ENDDO
      ENDDO
      PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'
      PRINT*,'sqf=', sqf

C-------------------------------------------------------------
      RETURN
      END
C_______________________________________________________________
C_______________________________________________________________
1	guez	3	!
2			! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advz.F,v 1.2 2005/05/25 13:10:09 fairhead Exp $
3			!
4			SUBROUTINE advz(limit,dtz,w,sm,s0,sx,sy,sz)
5			use dimens_m
6			use paramet_m
7			use comconst
8			use comvert
9			IMPLICIT NONE
10
11			CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
12			C C
13			C first-order moments (FOM) advection of tracer in Z direction C
14			C C
15			C Source : Pascal Simon (Meteo,CNRM) C
16			C Adaptation : A.Armengaud (LGGE) juin 94 C
17			C C
18			C C
19			C sont des arguments d'entree pour le s-pg... C
20			C C
21			C dq est l'argument de sortie pour le s-pg C
22			C C
23			CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
24			C
25			C parametres principaux du modele
26			C
27
28			C Arguments :
29			C -----------
30			C dtz : frequence fictive d'appel du transport
31			C w : flux de masse en z en Pa.m2.s-1
32
33			INTEGER ntra
34			PARAMETER (ntra = 1)
35
36	guez	28	REAL, intent(in):: dtz
37	guez	3	REAL w ( iip1,jjp1,llm )
38
39			C moments: SM total mass in each grid box
40			C S0 mass of tracer in each grid box
41			C Si 1rst order moment in i direction
42			C
43			REAL SM(iip1,jjp1,llm)
44			+ ,S0(iip1,jjp1,llm,ntra)
45			REAL sx(iip1,jjp1,llm,ntra)
46			+ ,sy(iip1,jjp1,llm,ntra)
47			+ ,sz(iip1,jjp1,llm,ntra)
48
49
50			C Local :
51			C -------
52
53			C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
54			C mass fluxes in kg
55			C declaration :
56
57			REAL WGRI(iip1,jjp1,0:llm)
58
59			C
60			C the moments F are used as temporary storage for
61			C portions of grid boxes in transit at the current latitude
62			C
63			REAL FM(iim,llm)
64			REAL F0(iim,llm,ntra),FX(iim,llm,ntra)
65			REAL FY(iim,llm,ntra),FZ(iim,llm,ntra)
66			C
67			C work arrays
68			C
69			REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
70			REAL TEMPTM ! Just temporal variable
71			REAL sqi,sqf
72			C
73			LOGICAL LIMIT
74			INTEGER lon,lat,niv
75			INTEGER i,j,jv,k,l,lp
76
77			lon = iim
78			lat = jjp1
79			niv = llm
80
81			C *** Test : diag de la qqtite totale de traceur
82			C dans l'atmosphere avant l'advection en z
83			sqi = 0.
84			sqf = 0.
85
86			DO l = 1,llm
87			DO j = 1,jjp1
88			DO i = 1,iim
89			cIM 240305 sqi = sqi + S0(i,j,l,9)
90			sqi = sqi + S0(i,j,l,ntra)
91			ENDDO
92			ENDDO
93			ENDDO
94			PRINT*,'-------- DIAG DANS ADVZ - ENTREE ---------'
95			PRINT*,'sqi=',sqi
96
97			C-----------------------------------------------------------------
98			C Interface : adaptation nouveau modele
99			C -------------------------------------
100			C
101			C Conversion du flux de masse en kg.s-1
102
103			DO 500 l = 1,llm
104			DO 500 j = 1,jjp1
105			DO 500 i = 1,iip1
106			c wgri (i,j,llm+1-l) = w (i,j,l) / g
107			wgri (i,j,llm+1-l) = w (i,j,l)
108			c wgri (i,j,0) = 0. ! a detruire ult.
109			c wgri (i,j,l) = 0.1 ! w (i,j,l)
110			c wgri (i,j,llm) = 0. ! a detruire ult.
111			500 CONTINUE
112			DO j = 1,jjp1
113			DO i = 1,iip1
114			wgri(i,j,0)=0.
115			enddo
116			enddo
117
118			C-----------------------------------------------------------------
119
120			C start here
121			C boucle sur les latitudes
122			C
123			DO 1 K=1,LAT
124			C
125			C place limits on appropriate moments before transport
126			C (if flux-limiting is to be applied)
127			C
128			IF(.NOT.LIMIT) GO TO 101
129			C
130			DO 10 JV=1,NTRA
131			DO 10 L=1,NIV
132			DO 100 I=1,LON
133			sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
134			+ ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
135			100 CONTINUE
136			10 CONTINUE
137			C
138			101 CONTINUE
139			C
140			C boucle sur les niveaux intercouches de 1 a NIV-1
141			C (flux nul au sommet L=0 et a la base L=NIV)
142			C
143			C calculate flux and moments between adjacent boxes
144			C (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
145			C 1- create temporary moments/masses for partial boxes in transit
146			C 2- reajusts moments remaining in the box
147			C
148			DO 11 L=1,NIV-1
149			LP=L+1
150			C
151			DO 110 I=1,LON
152			C
153			IF(WGRI(I,K,L).LT.0.) THEN
154			FM(I,L)=-WGRI(I,K,L)*DTZ
155			ALF(I)=FM(I,L)/SM(I,K,LP)
156			SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
157			ELSE
158			FM(I,L)=WGRI(I,K,L)*DTZ
159			ALF(I)=FM(I,L)/SM(I,K,L)
160			SM(I,K,L)=SM(I,K,L)-FM(I,L)
161			ENDIF
162			C
163			ALFQ (I)=ALF(I)*ALF(I)
164			ALF1 (I)=1.-ALF(I)
165			ALF1Q(I)=ALF1(I)*ALF1(I)
166			C
167			110 CONTINUE
168			C
169			DO 111 JV=1,NTRA
170			DO 1110 I=1,LON
171			C
172			IF(WGRI(I,K,L).LT.0.) THEN
173			C
174			F0(I,L,JV)=ALF (I)( S0(I,K,LP,JV)-ALF1(I)sz(I,K,LP,JV) )
175			FZ(I,L,JV)=ALFQ(I)*sz(I,K,LP,JV)
176			FX(I,L,JV)=ALF (I)*sx(I,K,LP,JV)
177			FY(I,L,JV)=ALF (I)*sy(I,K,LP,JV)
178			C
179			S0(I,K,LP,JV)=S0(I,K,LP,JV)-F0(I,L,JV)
180			sz(I,K,LP,JV)=ALF1Q(I)*sz(I,K,LP,JV)
181			sx(I,K,LP,JV)=sx(I,K,LP,JV)-FX(I,L,JV)
182			sy(I,K,LP,JV)=sy(I,K,LP,JV)-FY(I,L,JV)
183			C
184			ELSE
185			C
186			F0(I,L,JV)=ALF (I)(S0(I,K,L,JV)+ALF1(I)sz(I,K,L,JV) )
187			FZ(I,L,JV)=ALFQ(I)*sz(I,K,L,JV)
188			FX(I,L,JV)=ALF (I)*sx(I,K,L,JV)
189			FY(I,L,JV)=ALF (I)*sy(I,K,L,JV)
190			C
191			S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,L,JV)
192			sz(I,K,L,JV)=ALF1Q(I)*sz(I,K,L,JV)
193			sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,L,JV)
194			sy(I,K,L,JV)=sy(I,K,L,JV)-FY(I,L,JV)
195			C
196			ENDIF
197			C
198			1110 CONTINUE
199			111 CONTINUE
200			C
201			11 CONTINUE
202			C
203			C puts the temporary moments Fi into appropriate neighboring boxes
204			C
205			DO 12 L=1,NIV-1
206			LP=L+1
207			C
208			DO 120 I=1,LON
209			C
210			IF(WGRI(I,K,L).LT.0.) THEN
211			SM(I,K,L)=SM(I,K,L)+FM(I,L)
212			ALF(I)=FM(I,L)/SM(I,K,L)
213			ELSE
214			SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
215			ALF(I)=FM(I,L)/SM(I,K,LP)
216			ENDIF
217			C
218			ALF1(I)=1.-ALF(I)
219			ALFQ(I)=ALF(I)*ALF(I)
220			ALF1Q(I)=ALF1(I)*ALF1(I)
221			C
222			120 CONTINUE
223			C
224			DO 121 JV=1,NTRA
225			DO 1210 I=1,LON
226			C
227			IF(WGRI(I,K,L).LT.0.) THEN
228			C
229			TEMPTM=-ALF(I)S0(I,K,L,JV)+ALF1(I)F0(I,L,JV)
230			S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
231			sz(I,K,L,JV)=ALF(I)FZ(I,L,JV)+ALF1(I)sz(I,K,L,JV)+3.*TEMPTM
232			sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,L,JV)
233			sy(I,K,L,JV)=sy(I,K,L,JV)+FY(I,L,JV)
234			C
235			ELSE
236			C
237			TEMPTM=ALF(I)S0(I,K,LP,JV)-ALF1(I)F0(I,L,JV)
238			S0(I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
239			sz(I,K,LP,JV)=ALF(I)FZ(I,L,JV)+ALF1(I)sz(I,K,LP,JV)
240			+ +3.*TEMPTM
241			sx(I,K,LP,JV)=sx(I,K,LP,JV)+FX(I,L,JV)
242			sy(I,K,LP,JV)=sy(I,K,LP,JV)+FY(I,L,JV)
243			C
244			ENDIF
245			C
246			1210 CONTINUE
247			121 CONTINUE
248			C
249			12 CONTINUE
250			C
251			C fin de la boucle principale sur les latitudes
252			C
253			1 CONTINUE
254
255			C *** ------------------- bouclage cyclique en X ------------
256
257			c DO l = 1,llm
258			c DO j = 1,jjp1
259			c SM(iip1,j,l) = SM(1,j,l)
260			c S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
261			C sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
262			c sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
263			c sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
264			c ENDDO
265			c ENDDO
266
267			C-------------------------------------------------------------
268			C *** Test : diag de la qqtite totale de traceur
269			C dans l'atmosphere avant l'advection en z
270			DO l = 1,llm
271			DO j = 1,jjp1
272			DO i = 1,iim
273			cIM 240305 sqf = sqf + S0(i,j,l,9)
274			sqf = sqf + S0(i,j,l,ntra)
275			ENDDO
276			ENDDO
277			ENDDO
278			PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'
279			PRINT*,'sqf=', sqf
280
281			C-------------------------------------------------------------
282			RETURN
283			END
284			C_______________________________________________________________
285			C_______________________________________________________________