/[lmdze]/trunk/libf/dyn3d/advzp.f
ViewVC logotype

Contents of /trunk/libf/dyn3d/advzp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (show annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 2 months ago) by guez
File size: 12472 byte(s)
Initial import
1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advzp.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $
3 !
4 SUBROUTINE ADVZP(LIMIT,DTZ,W,SM,S0,SSX,SY,SZ
5 . ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra )
6
7 use dimens_m
8 use paramet_m
9 use comconst
10 use comvert
11 use comgeom
12 IMPLICIT NONE
13
14 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
15 C C
16 C second-order moments (SOM) advection of tracer in Z direction C
17 C C
18 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
19 C C
20 C Source : Pascal Simon ( Meteo, CNRM ) C
21 C Adaptation : A.A. (LGGE) C
22 C Derniere Modif : 19/11/95 LAST C
23 C C
24 C sont les arguments d'entree pour le s-pg C
25 C C
26 C argument de sortie du s-pg C
27 C C
28 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
29 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
30 C
31 C Rem : Probleme aux poles il faut reecrire ce cas specifique
32 C Attention au sens de l'indexation
33 C
34
35 C
36 C parametres principaux du modele
37 C
38 C
39 C Arguments :
40 C ----------
41 C dty : frequence fictive d'appel du transport
42 C parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
43 c
44 INTEGER lon,lat,niv
45 INTEGER i,j,jv,k,kp,l,lp
46 INTEGER ntra
47 c PARAMETER (ntra = 1)
48 c
49 REAL dtz
50 REAL w ( iip1,jjp1,llm )
51 c
52 C moments: SM total mass in each grid box
53 C S0 mass of tracer in each grid box
54 C Si 1rst order moment in i direction
55 C
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)
67 C
68 C Local :
69 C -------
70 C
71 C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
72 C mass fluxes in kg
73 C declaration :
74 C
75 REAL WGRI(iip1,jjp1,0:llm)
76
77 C Rem : UGRI et VGRI ne sont pas utilises dans
78 C cette subroutine ( advection en z uniquement )
79 C Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
80 C attention a celui de WGRI
81 C
82 C the moments F are similarly defined and used as temporary
83 C storage for portions of the grid boxes in transit
84 C
85 C the moments Fij are used as temporary storage for
86 C portions of the grid boxes in transit at the current level
87 C
88 C work arrays
89 C
90 C
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
99 C
100 C work arrays
101 C
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
108 c
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
116 c-----------------------------------------------------------------
117 C *** Test : diag de la qtite totale de traceur dans
118 C l'atmosphere avant l'advection en Y
119 c
120 sqi = 0.
121 sqf = 0.
122 c
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
133 c-----------------------------------------------------------------
134 C Interface : adaptation nouveau modele
135 C -------------------------------------
136 C
137 C 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
149 c
150 cAA rem : Je ne suis pas sur du signe
151 cAA Je ne suis pas sur pour le 0:llm
152 c
153 c-----------------------------------------------------------------
154 C---------------------- START HERE -------------------------------
155 C
156 C boucle sur les latitudes
157 C
158 DO 1 K=1,LAT
159 C
160 C place limits on appropriate moments before transport
161 C (if flux-limiting is to be applied)
162 C
163 IF(.NOT.LIMIT) GO TO 101
164 C
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
186 C
187 101 CONTINUE
188 C
189 C boucle sur les niveaux intercouches de 1 a NIV-1
190 C (flux nul au sommet L=0 et a la base L=NIV)
191 C
192 C calculate flux and moments between adjacent boxes
193 C (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
194 C 1- create temporary moments/masses for partial boxes in transit
195 C 2- reajusts moments remaining in the box
196 C
197 DO 11 L=1,NIV-1
198 LP=L+1
199 C
200 DO 110 I=1,LON
201 C
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
211 C
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)
218 C
219 110 CONTINUE
220 C
221 DO 111 JV=1,NTRA
222 DO 1110 I=1,LON
223 C
224 IF(WGRI(I,K,L).LT.0.) THEN
225 C
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)
237 C
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)
249 C
250 ELSE
251 C
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)
263 C
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)
274 C
275 ENDIF
276 C
277 1110 CONTINUE
278 111 CONTINUE
279 C
280 11 CONTINUE
281 C
282 C puts the temporary moments Fi into appropriate neighboring boxes
283 C
284 DO 12 L=1,NIV-1
285 LP=L+1
286 C
287 DO 120 I=1,LON
288 C
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
296 C
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)
302 C
303 120 CONTINUE
304 C
305 DO 121 JV=1,NTRA
306 DO 1210 I=1,LON
307 C
308 IF(WGRI(I,K,L).LT.0.) THEN
309 C
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)
325 C
326 ELSE
327 C
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)
343 C
344 ENDIF
345 C
346 1210 CONTINUE
347 121 CONTINUE
348 C
349 12 CONTINUE
350 C
351 C fin de la boucle principale sur les latitudes
352 C
353 1 CONTINUE
354 C
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
364 c C-------------------------------------------------------------
365 C *** Test : diag de la qqtite totale de tarceur
366 C 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

  ViewVC Help
Powered by ViewVC 1.1.21