libf/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	!
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	REAL, intent(in):: dtz
37	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_______________________________________________________________