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

Annotation of /trunk/Sources/dyn3d/advx.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 31 - (hide annotations)
Thu Apr 1 14:59:19 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/advx.f
File size: 12447 byte(s)
Split "vlsplt.f" in single-procedure files. Gathered the files in
directory "dyn3d/Vlsplt".

Defined "pbarum(:, 1, :)" and "pbarum(:, jjm + 1, :)" in procedure
"groupe".

1 guez 3 !
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 guez 31 REAL, intent(in):: pbaru ( iip1,jjp1,llm )
43 guez 3
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