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

Contents of /trunk/libf/dyn3d/advx.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: 12433 byte(s)
Initial import
1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advx.F,v 1.2 2005/05/25 13:10:09 fairhead Exp $
3 !
4 SUBROUTINE advx(limit,dtx,pbaru,sm,s0,
5 $ sx,sy,sz,lati,latf)
6 use dimens_m
7 use paramet_m
8 use comconst
9 use comvert
10 IMPLICIT NONE
11
12 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
13 C C
14 C first-order moments (FOM) advection of tracer in X direction C
15 C C
16 C Source : Pascal Simon (Meteo,CNRM) C
17 C Adaptation : A.Armengaud (LGGE) juin 94 C
18 C C
19 C limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C
20 C sont des arguments d'entree pour le s-pg... C
21 C C
22 C sm,s0,sx,sy,sz C
23 C sont les arguments de sortie pour le s-pg C
24 C C
25 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
26 C
27 C parametres principaux du modele
28 C
29
30 C Arguments :
31 C -----------
32 C dtx : frequence fictive d'appel du transport
33 C pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
34
35 INTEGER ntra
36 PARAMETER (ntra = 1)
37
38 C ATTENTION partout ou on trouve ntra, insertion de boucle
39 C possible dans l'avenir.
40
41 REAL dtx
42 REAL pbaru ( iip1,jjp1,llm )
43
44 C moments: SM total mass in each grid box
45 C S0 mass of tracer in each grid box
46 C Si 1rst order moment in i direction
47 C
48 REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
49 REAL sx(iip1,jjp1,llm,ntra)
50 $ ,sy(iip1,jjp1,llm,ntra)
51 REAL sz(iip1,jjp1,llm,ntra)
52
53 C Local :
54 C -------
55
56 C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
57 C mass fluxes in kg
58 C declaration :
59
60 REAL UGRI(iip1,jjp1,llm)
61
62 C Rem : VGRI et WGRI ne sont pas utilises dans
63 C cette subroutine ( advection en x uniquement )
64 C
65 C Ti are the moments for the current latitude and level
66 C
67 REAL TM(iim)
68 REAL T0(iim,ntra),TX(iim,ntra)
69 REAL TY(iim,ntra),TZ(iim,ntra)
70 REAL TEMPTM ! just a temporary variable
71 C
72 C the moments F are similarly defined and used as temporary
73 C storage for portions of the grid boxes in transit
74 C
75 REAL FM(iim)
76 REAL F0(iim,ntra),FX(iim,ntra)
77 REAL FY(iim,ntra),FZ(iim,ntra)
78 C
79 C work arrays
80 C
81 REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
82 C
83 REAL SMNEW(iim),UEXT(iim)
84 C
85 REAL sqi,sqf
86
87 LOGICAL LIMIT
88 INTEGER NUM(jjp1),LONK,NUMK
89 INTEGER lon,lati,latf,niv
90 INTEGER i,i2,i3,j,jv,l,k,itrac
91
92 lon = iim
93 niv = llm
94
95 C *** Test de passage d'arguments ******
96
97
98 C -------------------------------------
99 DO 300 j = 1,jjp1
100 NUM(j) = 1
101 300 CONTINUE
102 sqi = 0.
103 sqf = 0.
104
105 DO l = 1,llm
106 DO j = 1,jjp1
107 DO i = 1,iim
108 cIM 240305 sqi = sqi + S0(i,j,l,9)
109 sqi = sqi + S0(i,j,l,ntra)
110 ENDDO
111 ENDDO
112 ENDDO
113 PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
114 PRINT*,'sqi=',sqi
115
116
117 C Interface : adaptation nouveau modele
118 C -------------------------------------
119 C
120 C ---------------------------------------------------------
121 C Conversion des flux de masses en kg/s
122 C pbaru est en N/s d'ou :
123 C ugri est en kg/s
124
125 DO 500 l = 1,llm
126 DO 500 j = 1,jjm+1
127 DO 500 i = 1,iip1
128 C ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
129 ugri (i,j,llm+1-l) = pbaru (i,j,l)
130 500 CONTINUE
131
132
133 C ---------------------------------------------------------
134 C ---------------------------------------------------------
135 C ---------------------------------------------------------
136
137 C start here
138 C
139 C boucle principale sur les niveaux et les latitudes
140 C
141 DO 1 L=1,NIV
142 DO 1 K=lati,latf
143 C
144 C initialisation
145 C
146 C program assumes periodic boundaries in X
147 C
148 DO 10 I=2,LON
149 SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
150 10 CONTINUE
151 SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
152 C
153 C modifications for extended polar zones
154 C
155 NUMK=NUM(K)
156 LONK=LON/NUMK
157 C
158 IF(NUMK.GT.1) THEN
159 C
160 DO 111 I=1,LON
161 TM(I)=0.
162 111 CONTINUE
163 DO 112 JV=1,NTRA
164 DO 1120 I=1,LON
165 T0(I,JV)=0.
166 TX(I,JV)=0.
167 TY(I,JV)=0.
168 TZ(I,JV)=0.
169 1120 CONTINUE
170 112 CONTINUE
171 C
172 DO 11 I2=1,NUMK
173 C
174 DO 113 I=1,LONK
175 I3=(I-1)*NUMK+I2
176 TM(I)=TM(I)+SM(I3,K,L)
177 ALF(I)=SM(I3,K,L)/TM(I)
178 ALF1(I)=1.-ALF(I)
179 113 CONTINUE
180 C
181 DO JV=1,NTRA
182 DO I=1,LONK
183 I3=(I-1)*NUMK+I2
184 TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)
185 $ *S0(I3,K,L,JV)
186 T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
187 TX(I,JV)=ALF(I) *sx(I3,K,L,JV)+
188 $ ALF1(I)*TX(I,JV) +3.*TEMPTM
189 TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
190 TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
191 ENDDO
192 ENDDO
193 C
194 11 CONTINUE
195 C
196 ELSE
197 C
198 DO 115 I=1,LON
199 TM(I)=SM(I,K,L)
200 115 CONTINUE
201 DO 116 JV=1,NTRA
202 DO 1160 I=1,LON
203 T0(I,JV)=S0(I,K,L,JV)
204 TX(I,JV)=sx(I,K,L,JV)
205 TY(I,JV)=sy(I,K,L,JV)
206 TZ(I,JV)=sz(I,K,L,JV)
207 1160 CONTINUE
208 116 CONTINUE
209 C
210 ENDIF
211 C
212 DO 117 I=1,LONK
213 UEXT(I)=UGRI(I*NUMK,K,L)
214 117 CONTINUE
215 C
216 C place limits on appropriate moments before transport
217 C (if flux-limiting is to be applied)
218 C
219 IF(.NOT.LIMIT) GO TO 13
220 C
221 DO 12 JV=1,NTRA
222 DO 120 I=1,LONK
223 TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
224 120 CONTINUE
225 12 CONTINUE
226 C
227 13 CONTINUE
228 C
229 C calculate flux and moments between adjacent boxes
230 C 1- create temporary moments/masses for partial boxes in transit
231 C 2- reajusts moments remaining in the box
232 C
233 C flux from IP to I if U(I).lt.0
234 C
235 DO 140 I=1,LONK-1
236 IF(UEXT(I).LT.0.) THEN
237 FM(I)=-UEXT(I)*DTX
238 ALF(I)=FM(I)/TM(I+1)
239 TM(I+1)=TM(I+1)-FM(I)
240 ENDIF
241 140 CONTINUE
242 C
243 I=LONK
244 IF(UEXT(I).LT.0.) THEN
245 FM(I)=-UEXT(I)*DTX
246 ALF(I)=FM(I)/TM(1)
247 TM(1)=TM(1)-FM(I)
248 ENDIF
249 C
250 C flux from I to IP if U(I).gt.0
251 C
252 DO 141 I=1,LONK
253 IF(UEXT(I).GE.0.) THEN
254 FM(I)=UEXT(I)*DTX
255 ALF(I)=FM(I)/TM(I)
256 TM(I)=TM(I)-FM(I)
257 ENDIF
258 141 CONTINUE
259 C
260 DO 142 I=1,LONK
261 ALFQ(I)=ALF(I)*ALF(I)
262 ALF1(I)=1.-ALF(I)
263 ALF1Q(I)=ALF1(I)*ALF1(I)
264 142 CONTINUE
265 C
266 DO 150 JV=1,NTRA
267 DO 1500 I=1,LONK-1
268 C
269 IF(UEXT(I).LT.0.) THEN
270 C
271 F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
272 FX(I,JV)=ALFQ(I)*TX(I+1,JV)
273 FY(I,JV)=ALF (I)*TY(I+1,JV)
274 FZ(I,JV)=ALF (I)*TZ(I+1,JV)
275 C
276 T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
277 TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
278 TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
279 TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
280 C
281 ENDIF
282 C
283 1500 CONTINUE
284 150 CONTINUE
285 C
286 I=LONK
287 IF(UEXT(I).LT.0.) THEN
288 C
289 DO 151 JV=1,NTRA
290 C
291 F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
292 FX (I,JV)=ALFQ(I)*TX(1,JV)
293 FY (I,JV)=ALF (I)*TY(1,JV)
294 FZ (I,JV)=ALF (I)*TZ(1,JV)
295 C
296 T0(1,JV)=T0(1,JV)-F0(I,JV)
297 TX(1,JV)=ALF1Q(I)*TX(1,JV)
298 TY(1,JV)=TY(1,JV)-FY(I,JV)
299 TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
300 C
301 151 CONTINUE
302 C
303 ENDIF
304 C
305 DO 152 JV=1,NTRA
306 DO 1520 I=1,LONK
307 C
308 IF(UEXT(I).GE.0.) THEN
309 C
310 F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
311 FX(I,JV)=ALFQ(I)*TX(I,JV)
312 FY(I,JV)=ALF (I)*TY(I,JV)
313 FZ(I,JV)=ALF (I)*TZ(I,JV)
314 C
315 T0(I,JV)=T0(I,JV)-F0(I,JV)
316 TX(I,JV)=ALF1Q(I)*TX(I,JV)
317 TY(I,JV)=TY(I,JV)-FY(I,JV)
318 TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
319 C
320 ENDIF
321 C
322 1520 CONTINUE
323 152 CONTINUE
324 C
325 C puts the temporary moments Fi into appropriate neighboring boxes
326 C
327 DO 160 I=1,LONK
328 IF(UEXT(I).LT.0.) THEN
329 TM(I)=TM(I)+FM(I)
330 ALF(I)=FM(I)/TM(I)
331 ENDIF
332 160 CONTINUE
333 C
334 DO 161 I=1,LONK-1
335 IF(UEXT(I).GE.0.) THEN
336 TM(I+1)=TM(I+1)+FM(I)
337 ALF(I)=FM(I)/TM(I+1)
338 ENDIF
339 161 CONTINUE
340 C
341 I=LONK
342 IF(UEXT(I).GE.0.) THEN
343 TM(1)=TM(1)+FM(I)
344 ALF(I)=FM(I)/TM(1)
345 ENDIF
346 C
347 DO 162 I=1,LONK
348 ALF1(I)=1.-ALF(I)
349 162 CONTINUE
350 C
351 DO 170 JV=1,NTRA
352 DO 1700 I=1,LONK
353 C
354 IF(UEXT(I).LT.0.) THEN
355 C
356 TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
357 T0(I,JV)=T0(I,JV)+F0(I,JV)
358 TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
359 TY(I,JV)=TY(I,JV)+FY(I,JV)
360 TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
361 C
362 ENDIF
363 C
364 1700 CONTINUE
365 170 CONTINUE
366 C
367 DO 171 JV=1,NTRA
368 DO 1710 I=1,LONK-1
369 C
370 IF(UEXT(I).GE.0.) THEN
371 C
372 TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
373 T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
374 TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
375 TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
376 TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
377 C
378 ENDIF
379 C
380 1710 CONTINUE
381 171 CONTINUE
382 C
383 I=LONK
384 IF(UEXT(I).GE.0.) THEN
385 DO 172 JV=1,NTRA
386 TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
387 T0(1,JV)=T0(1,JV)+F0(I,JV)
388 TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
389 TY(1,JV)=TY(1,JV)+FY(I,JV)
390 TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
391 172 CONTINUE
392 ENDIF
393 C
394 C retour aux mailles d'origine (passage des Tij aux Sij)
395 C
396 IF(NUMK.GT.1) THEN
397 C
398 DO 180 I2=1,NUMK
399 C
400 DO 180 I=1,LONK
401 C
402 I3=I2+(I-1)*NUMK
403 SM(I3,K,L)=SMNEW(I3)
404 ALF(I)=SMNEW(I3)/TM(I)
405 TM(I)=TM(I)-SMNEW(I3)
406 C
407 ALFQ(I)=ALF(I)*ALF(I)
408 ALF1(I)=1.-ALF(I)
409 ALF1Q(I)=ALF1(I)*ALF1(I)
410 C
411 180 CONTINUE
412 C
413 DO JV=1,NTRA
414 DO I=1,LONK
415 C
416 I3=I2+(I-1)*NUMK
417 S0(I3,K,L,JV)=ALF (I)
418 $ * (T0(I,JV)-ALF1(I)*TX(I,JV))
419 sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
420 sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
421 sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
422 C
423 C reajusts moments remaining in the box
424 C
425 T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
426 TX(I,JV)=ALF1Q(I)*TX(I,JV)
427 TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
428 TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
429 ENDDO
430 ENDDO
431 C
432 C
433 ELSE
434 C
435 DO 190 I=1,LON
436 SM(I,K,L)=TM(I)
437 190 CONTINUE
438 DO 191 JV=1,NTRA
439 DO 1910 I=1,LON
440 S0(I,K,L,JV)=T0(I,JV)
441 sx(I,K,L,JV)=TX(I,JV)
442 sy(I,K,L,JV)=TY(I,JV)
443 sz(I,K,L,JV)=TZ(I,JV)
444 1910 CONTINUE
445 191 CONTINUE
446 C
447 ENDIF
448 C
449 1 CONTINUE
450 C
451 C ----------- AA Test en fin de ADVX ------ Controle des S*
452 c OK
453 c DO 9998 l = 1, llm
454 c DO 9998 j = 1, jjp1
455 c DO 9998 i = 1, iip1
456 c IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
457 c PRINT*, '-------------------'
458 c PRINT*, 'En fin de ADVX'
459 c PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
460 c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
461 c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
462 c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
463 c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
464 c WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
465 cc STOP
466 c ENDIF
467 c 9998 CONTINUE
468 c
469 C ---------- bouclage cyclique
470 DO itrac=1,ntra
471 DO l = 1,llm
472 DO j = lati,latf
473 SM(iip1,j,l) = SM(1,j,l)
474 S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
475 sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
476 sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
477 sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
478 END DO
479 END DO
480 ENDDO
481
482 c ----------- qqtite totale de traceur dans tte l'atmosphere
483 DO l = 1, llm
484 DO j = 1, jjp1
485 DO i = 1, iim
486 cIM 240405 sqf = sqf + S0(i,j,l,9)
487 sqf = sqf + S0(i,j,l,ntra)
488 END DO
489 END DO
490 END DO
491 c
492 PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
493 PRINT*,'sqf=',sqf
494 c-------------
495
496 RETURN
497 END
498 C_________________________________________________________________
499 C_________________________________________________________________

  ViewVC Help
Powered by ViewVC 1.1.21