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

Annotation of /trunk/dyn3d/advzp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 3 months ago) by guez
File size: 12460 byte(s)
Changed all ".f90" suffixes to ".f".
1 guez 3
2 guez 81 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advzp.F,v 1.1.1.1 2004/05/19
3     ! 12:53:06 lmdzadmin Exp $
4 guez 3
5 guez 81 SUBROUTINE advzp(limit, dtz, w, sm, s0, ssx, sy, sz, ssxx, ssxy, ssxz, syy, &
6     syz, szz, ntra)
7 guez 3
8 guez 81 USE dimens_m
9     USE paramet_m
10     USE comconst
11     USE disvert_m
12     USE comgeom
13     IMPLICIT NONE
14 guez 3
15 guez 81 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
16     ! C
17     ! second-order moments (SOM) advection of tracer in Z direction C
18     ! C
19     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
20     ! C
21     ! Source : Pascal Simon ( Meteo, CNRM ) C
22     ! Adaptation : A.A. (LGGE) C
23     ! Derniere Modif : 19/11/95 LAST C
24     ! C
25     ! sont les arguments d'entree pour le s-pg C
26     ! C
27     ! argument de sortie du s-pg C
28     ! C
29     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
30     ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
31 guez 3
32 guez 81 ! Rem : Probleme aux poles il faut reecrire ce cas specifique
33     ! Attention au sens de l'indexation
34    
35    
36    
37     ! parametres principaux du modele
38    
39    
40     ! Arguments :
41     ! ----------
42     ! dty : frequence fictive d'appel du transport
43     ! parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
44    
45     INTEGER lon, lat, niv
46     INTEGER i, j, jv, k, kp, l, lp
47     INTEGER ntra
48     ! PARAMETER (ntra = 1)
49    
50     REAL dtz
51     REAL w(iip1, jjp1, llm)
52    
53     ! moments: SM total mass in each grid box
54     ! S0 mass of tracer in each grid box
55     ! Si 1rst order moment in i direction
56    
57     REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
58     REAL ssx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra), &
59     sz(iip1, jjp1, llm, ntra), ssxx(iip1, jjp1, llm, ntra), &
60     ssxy(iip1, jjp1, llm, ntra), ssxz(iip1, jjp1, llm, ntra), &
61     syy(iip1, jjp1, llm, ntra), syz(iip1, jjp1, llm, ntra), &
62     szz(iip1, jjp1, llm, ntra)
63    
64     ! Local :
65     ! -------
66    
67     ! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
68     ! mass fluxes in kg
69     ! declaration :
70    
71     REAL wgri(iip1, jjp1, 0:llm)
72    
73     ! Rem : UGRI et VGRI ne sont pas utilises dans
74     ! cette subroutine ( advection en z uniquement )
75     ! Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
76     ! attention a celui de WGRI
77    
78     ! the moments F are similarly defined and used as temporary
79     ! storage for portions of the grid boxes in transit
80    
81     ! the moments Fij are used as temporary storage for
82     ! portions of the grid boxes in transit at the current level
83    
84     ! work arrays
85    
86    
87     REAL f0(iim, llm, ntra), fm(iim, llm)
88     REAL fx(iim, llm, ntra), fy(iim, llm, ntra)
89     REAL fz(iim, llm, ntra)
90     REAL fxx(iim, llm, ntra), fxy(iim, llm, ntra)
91     REAL fxz(iim, llm, ntra), fyy(iim, llm, ntra)
92     REAL fyz(iim, llm, ntra), fzz(iim, llm, ntra)
93     REAL s00(ntra)
94     REAL sm0 ! Just temporal variable
95    
96     ! work arrays
97    
98     REAL alf(iim), alf1(iim)
99     REAL alfq(iim), alf1q(iim)
100     REAL alf2(iim), alf3(iim)
101     REAL alf4(iim)
102     REAL temptm ! Just temporal variable
103     REAL slpmax, s1max, s1new, s2new
104    
105     REAL sqi, sqf
106     LOGICAL limit
107    
108     lon = iim ! rem : Il est possible qu'un pbl. arrive ici
109     lat = jjp1 ! a cause des dim. differentes entre les
110     niv = llm ! tab. S et VGRI
111    
112     ! -----------------------------------------------------------------
113     ! *** Test : diag de la qtite totale de traceur dans
114     ! l'atmosphere avant l'advection en Y
115    
116     sqi = 0.
117     sqf = 0.
118    
119     DO l = 1, llm
120     DO j = 1, jjp1
121     DO i = 1, iim
122     sqi = sqi + s0(i, j, l, ntra)
123 guez 3 END DO
124 guez 81 END DO
125     END DO
126     PRINT *, '---------- DIAG DANS ADVZP - ENTREE --------'
127     PRINT *, 'sqi=', sqi
128 guez 3
129 guez 81 ! -----------------------------------------------------------------
130     ! Interface : adaptation nouveau modele
131     ! -------------------------------------
132 guez 3
133 guez 81 ! Conversion des flux de masses en kg
134 guez 3
135 guez 81 DO l = 1, llm
136     DO j = 1, jjp1
137     DO i = 1, iip1
138     wgri(i, j, llm+1-l) = w(i, j, l)
139     END DO
140     END DO
141     END DO
142     DO j = 1, jjp1
143     DO i = 1, iip1
144     wgri(i, j, 0) = 0.
145     END DO
146     END DO
147    
148     ! AA rem : Je ne suis pas sur du signe
149     ! AA Je ne suis pas sur pour le 0:llm
150    
151     ! -----------------------------------------------------------------
152     ! ---------------------- START HERE -------------------------------
153    
154     ! boucle sur les latitudes
155    
156     DO k = 1, lat
157    
158     ! place limits on appropriate moments before transport
159     ! (if flux-limiting is to be applied)
160    
161     IF (.NOT. limit) GO TO 101
162    
163     DO jv = 1, ntra
164     DO l = 1, niv
165     DO i = 1, lon
166     IF (s0(i,k,l,jv)>0.) THEN
167     slpmax = s0(i, k, l, jv)
168     s1max = 1.5*slpmax
169     s1new = amin1(s1max, amax1(-s1max,sz(i,k,l,jv)))
170     s2new = amin1(2.*slpmax-abs(s1new)/3., amax1(abs( &
171     s1new)-slpmax,szz(i,k,l,jv)))
172     sz(i, k, l, jv) = s1new
173     szz(i, k, l, jv) = s2new
174     ssxz(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,ssxz(i,k,l,jv)))
175     syz(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,syz(i,k,l,jv)))
176     ELSE
177     sz(i, k, l, jv) = 0.
178     szz(i, k, l, jv) = 0.
179     ssxz(i, k, l, jv) = 0.
180     syz(i, k, l, jv) = 0.
181     END IF
182     END DO
183     END DO
184     END DO
185    
186     101 CONTINUE
187    
188     ! boucle sur les niveaux intercouches de 1 a NIV-1
189     ! (flux nul au sommet L=0 et a la base L=NIV)
190    
191     ! calculate flux and moments between adjacent boxes
192     ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
193     ! 1- create temporary moments/masses for partial boxes in transit
194     ! 2- reajusts moments remaining in the box
195    
196     DO l = 1, niv - 1
197     lp = l + 1
198    
199     DO i = 1, lon
200    
201     IF (wgri(i,k,l)<0.) THEN
202     fm(i, l) = -wgri(i, k, l)*dtz
203     alf(i) = fm(i, l)/sm(i, k, lp)
204     sm(i, k, lp) = sm(i, k, lp) - fm(i, l)
205     ELSE
206     fm(i, l) = wgri(i, k, l)*dtz
207     alf(i) = fm(i, l)/sm(i, k, l)
208     sm(i, k, l) = sm(i, k, l) - fm(i, l)
209     END IF
210    
211     alfq(i) = alf(i)*alf(i)
212     alf1(i) = 1. - alf(i)
213     alf1q(i) = alf1(i)*alf1(i)
214     alf2(i) = alf1(i) - alf(i)
215     alf3(i) = alf(i)*alfq(i)
216     alf4(i) = alf1(i)*alf1q(i)
217    
218     END DO
219    
220     DO jv = 1, ntra
221     DO i = 1, lon
222    
223     IF (wgri(i,k,l)<0.) THEN
224    
225     f0(i, l, jv) = alf(i)*(s0(i,k,lp,jv)-alf1(i)*(sz(i,k,lp, &
226     jv)-alf2(i)*szz(i,k,lp,jv)))
227     fz(i, l, jv) = alfq(i)*(sz(i,k,lp,jv)-3.*alf1(i)*szz(i,k,lp,jv))
228     fzz(i, l, jv) = alf3(i)*szz(i, k, lp, jv)
229     fxz(i, l, jv) = alfq(i)*ssxz(i, k, lp, jv)
230     fyz(i, l, jv) = alfq(i)*syz(i, k, lp, jv)
231     fx(i, l, jv) = alf(i)*(ssx(i,k,lp,jv)-alf1(i)*ssxz(i,k,lp,jv))
232     fy(i, l, jv) = alf(i)*(sy(i,k,lp,jv)-alf1(i)*syz(i,k,lp,jv))
233     fxx(i, l, jv) = alf(i)*ssxx(i, k, lp, jv)
234     fxy(i, l, jv) = alf(i)*ssxy(i, k, lp, jv)
235     fyy(i, l, jv) = alf(i)*syy(i, k, lp, jv)
236    
237     s0(i, k, lp, jv) = s0(i, k, lp, jv) - f0(i, l, jv)
238     sz(i, k, lp, jv) = alf1q(i)*(sz(i,k,lp,jv)+3.*alf(i)*szz(i,k,lp, &
239     jv))
240     szz(i, k, lp, jv) = alf4(i)*szz(i, k, lp, jv)
241     ssxz(i, k, lp, jv) = alf1q(i)*ssxz(i, k, lp, jv)
242     syz(i, k, lp, jv) = alf1q(i)*syz(i, k, lp, jv)
243     ssx(i, k, lp, jv) = ssx(i, k, lp, jv) - fx(i, l, jv)
244     sy(i, k, lp, jv) = sy(i, k, lp, jv) - fy(i, l, jv)
245     ssxx(i, k, lp, jv) = ssxx(i, k, lp, jv) - fxx(i, l, jv)
246     ssxy(i, k, lp, jv) = ssxy(i, k, lp, jv) - fxy(i, l, jv)
247     syy(i, k, lp, jv) = syy(i, k, lp, jv) - fyy(i, l, jv)
248    
249     ELSE
250    
251     f0(i, l, jv) = alf(i)*(s0(i,k,l,jv)+alf1(i)*(sz(i,k,l, &
252     jv)+alf2(i)*szz(i,k,l,jv)))
253     fz(i, l, jv) = alfq(i)*(sz(i,k,l,jv)+3.*alf1(i)*szz(i,k,l,jv))
254     fzz(i, l, jv) = alf3(i)*szz(i, k, l, jv)
255     fxz(i, l, jv) = alfq(i)*ssxz(i, k, l, jv)
256     fyz(i, l, jv) = alfq(i)*syz(i, k, l, jv)
257     fx(i, l, jv) = alf(i)*(ssx(i,k,l,jv)+alf1(i)*ssxz(i,k,l,jv))
258     fy(i, l, jv) = alf(i)*(sy(i,k,l,jv)+alf1(i)*syz(i,k,l,jv))
259     fxx(i, l, jv) = alf(i)*ssxx(i, k, l, jv)
260     fxy(i, l, jv) = alf(i)*ssxy(i, k, l, jv)
261     fyy(i, l, jv) = alf(i)*syy(i, k, l, jv)
262    
263     s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, l, jv)
264     sz(i, k, l, jv) = alf1q(i)*(sz(i,k,l,jv)-3.*alf(i)*szz(i,k,l,jv))
265     szz(i, k, l, jv) = alf4(i)*szz(i, k, l, jv)
266     ssxz(i, k, l, jv) = alf1q(i)*ssxz(i, k, l, jv)
267     syz(i, k, l, jv) = alf1q(i)*syz(i, k, l, jv)
268     ssx(i, k, l, jv) = ssx(i, k, l, jv) - fx(i, l, jv)
269     sy(i, k, l, jv) = sy(i, k, l, jv) - fy(i, l, jv)
270     ssxx(i, k, l, jv) = ssxx(i, k, l, jv) - fxx(i, l, jv)
271     ssxy(i, k, l, jv) = ssxy(i, k, l, jv) - fxy(i, l, jv)
272     syy(i, k, l, jv) = syy(i, k, l, jv) - fyy(i, l, jv)
273    
274     END IF
275    
276     END DO
277     END DO
278    
279     END DO
280    
281     ! puts the temporary moments Fi into appropriate neighboring boxes
282    
283     DO l = 1, niv - 1
284     lp = l + 1
285    
286     DO i = 1, lon
287    
288     IF (wgri(i,k,l)<0.) THEN
289     sm(i, k, l) = sm(i, k, l) + fm(i, l)
290     alf(i) = fm(i, l)/sm(i, k, l)
291     ELSE
292     sm(i, k, lp) = sm(i, k, lp) + fm(i, l)
293     alf(i) = fm(i, l)/sm(i, k, lp)
294     END IF
295    
296     alf1(i) = 1. - alf(i)
297     alfq(i) = alf(i)*alf(i)
298     alf1q(i) = alf1(i)*alf1(i)
299     alf2(i) = alf(i)*alf1(i)
300     alf3(i) = alf1(i) - alf(i)
301    
302     END DO
303    
304     DO jv = 1, ntra
305     DO i = 1, lon
306    
307     IF (wgri(i,k,l)<0.) THEN
308    
309     temptm = -alf(i)*s0(i, k, l, jv) + alf1(i)*f0(i, l, jv)
310     s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, l, jv)
311     szz(i, k, l, jv) = alfq(i)*fzz(i, l, jv) + &
312     alf1q(i)*szz(i, k, l, jv) + 5.*(alf2(i)*(fz(i,l,jv)-sz(i,k,l, &
313     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) + &
317     alf1(i)*ssxz(i, k, l, jv) + 3.*(alf1(i)*fx(i,l,jv)-alf(i)*ssx(i &
318     ,k,l,jv))
319     syz(i, k, l, jv) = alf(i)*fyz(i, l, jv) + &
320     alf1(i)*syz(i, k, l, jv) + 3.*(alf1(i)*fy(i,l,jv)-alf(i)*sy(i,k &
321     ,l,jv))
322     ssx(i, k, l, jv) = ssx(i, k, l, jv) + fx(i, l, jv)
323     sy(i, k, l, jv) = sy(i, k, l, jv) + fy(i, l, jv)
324     ssxx(i, k, l, jv) = ssxx(i, k, l, jv) + fxx(i, l, jv)
325     ssxy(i, k, l, jv) = ssxy(i, k, l, jv) + fxy(i, l, jv)
326     syy(i, k, l, jv) = syy(i, k, l, jv) + fyy(i, l, jv)
327    
328     ELSE
329    
330     temptm = alf(i)*s0(i, k, lp, jv) - alf1(i)*f0(i, l, jv)
331     s0(i, k, lp, jv) = s0(i, k, lp, jv) + f0(i, l, jv)
332     szz(i, k, lp, jv) = alfq(i)*fzz(i, l, jv) + &
333     alf1q(i)*szz(i, k, lp, jv) + 5.*(alf2(i)*(sz(i,k,lp,jv)-fz(i,l, &
334     jv))-alf3(i)*temptm)
335     sz(i, k, lp, jv) = alf(i)*fz(i, l, jv) + &
336     alf1(i)*sz(i, k, lp, jv) + 3.*temptm
337     ssxz(i, k, lp, jv) = alf(i)*fxz(i, l, jv) + &
338     alf1(i)*ssxz(i, k, lp, jv) + 3.*(alf(i)*ssx(i,k,lp,jv)-alf1(i)* &
339     fx(i,l,jv))
340     syz(i, k, lp, jv) = alf(i)*fyz(i, l, jv) + &
341     alf1(i)*syz(i, k, lp, jv) + 3.*(alf(i)*sy(i,k,lp,jv)-alf1(i)*fy &
342     (i,l,jv))
343     ssx(i, k, lp, jv) = ssx(i, k, lp, jv) + fx(i, l, jv)
344     sy(i, k, lp, jv) = sy(i, k, lp, jv) + fy(i, l, jv)
345     ssxx(i, k, lp, jv) = ssxx(i, k, lp, jv) + fxx(i, l, jv)
346     ssxy(i, k, lp, jv) = ssxy(i, k, lp, jv) + fxy(i, l, jv)
347     syy(i, k, lp, jv) = syy(i, k, lp, jv) + fyy(i, l, jv)
348    
349     END IF
350    
351     END DO
352     END DO
353    
354     END DO
355    
356     ! fin de la boucle principale sur les latitudes
357    
358     END DO
359    
360     DO l = 1, llm
361     DO j = 1, jjp1
362     sm(iip1, j, l) = sm(1, j, l)
363     s0(iip1, j, l, ntra) = s0(1, j, l, ntra)
364     ssx(iip1, j, l, ntra) = ssx(1, j, l, ntra)
365     sy(iip1, j, l, ntra) = sy(1, j, l, ntra)
366     sz(iip1, j, l, ntra) = sz(1, j, l, ntra)
367     END DO
368     END DO
369     ! C-------------------------------------------------------------
370     ! *** Test : diag de la qqtite totale de tarceur
371     ! dans l'atmosphere avant l'advection en z
372     DO l = 1, llm
373     DO j = 1, jjp1
374     DO i = 1, iim
375     sqf = sqf + s0(i, j, l, ntra)
376     END DO
377     END DO
378     END DO
379     PRINT *, '-------- DIAG DANS ADVZ - SORTIE ---------'
380     PRINT *, 'sqf=', sqf
381    
382     RETURN
383     END SUBROUTINE advzp

  ViewVC Help
Powered by ViewVC 1.1.21