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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 134 - (hide annotations)
Wed Apr 29 15:47:56 2015 UTC (9 years, 1 month ago) by guez
File size: 12403 byte(s)
Sources inside, compilation outside.
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 guez 105 INTEGER i, j, jv, k, l, lp
47 guez 81 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    
94     ! work arrays
95    
96     REAL alf(iim), alf1(iim)
97     REAL alfq(iim), alf1q(iim)
98     REAL alf2(iim), alf3(iim)
99     REAL alf4(iim)
100     REAL temptm ! Just temporal variable
101     REAL slpmax, s1max, s1new, s2new
102    
103     REAL sqi, sqf
104     LOGICAL limit
105    
106     lon = iim ! rem : Il est possible qu'un pbl. arrive ici
107     lat = jjp1 ! a cause des dim. differentes entre les
108     niv = llm ! tab. S et VGRI
109    
110     ! -----------------------------------------------------------------
111     ! *** Test : diag de la qtite totale de traceur dans
112     ! l'atmosphere avant l'advection en Y
113    
114     sqi = 0.
115     sqf = 0.
116    
117     DO l = 1, llm
118     DO j = 1, jjp1
119     DO i = 1, iim
120     sqi = sqi + s0(i, j, l, ntra)
121 guez 3 END DO
122 guez 81 END DO
123     END DO
124     PRINT *, '---------- DIAG DANS ADVZP - ENTREE --------'
125     PRINT *, 'sqi=', sqi
126 guez 3
127 guez 81 ! -----------------------------------------------------------------
128     ! Interface : adaptation nouveau modele
129     ! -------------------------------------
130 guez 3
131 guez 81 ! Conversion des flux de masses en kg
132 guez 3
133 guez 81 DO l = 1, llm
134     DO j = 1, jjp1
135     DO i = 1, iip1
136     wgri(i, j, llm+1-l) = w(i, j, l)
137     END DO
138     END DO
139     END DO
140     DO j = 1, jjp1
141     DO i = 1, iip1
142     wgri(i, j, 0) = 0.
143     END DO
144     END DO
145    
146     ! AA rem : Je ne suis pas sur du signe
147     ! AA Je ne suis pas sur pour le 0:llm
148    
149     ! -----------------------------------------------------------------
150     ! ---------------------- START HERE -------------------------------
151    
152     ! boucle sur les latitudes
153    
154     DO k = 1, lat
155    
156     ! place limits on appropriate moments before transport
157     ! (if flux-limiting is to be applied)
158    
159     IF (.NOT. limit) GO TO 101
160    
161     DO jv = 1, ntra
162     DO l = 1, niv
163     DO i = 1, lon
164     IF (s0(i,k,l,jv)>0.) THEN
165     slpmax = s0(i, k, l, jv)
166     s1max = 1.5*slpmax
167     s1new = amin1(s1max, amax1(-s1max,sz(i,k,l,jv)))
168     s2new = amin1(2.*slpmax-abs(s1new)/3., amax1(abs( &
169     s1new)-slpmax,szz(i,k,l,jv)))
170     sz(i, k, l, jv) = s1new
171     szz(i, k, l, jv) = s2new
172     ssxz(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,ssxz(i,k,l,jv)))
173     syz(i, k, l, jv) = amin1(slpmax, amax1(-slpmax,syz(i,k,l,jv)))
174     ELSE
175     sz(i, k, l, jv) = 0.
176     szz(i, k, l, jv) = 0.
177     ssxz(i, k, l, jv) = 0.
178     syz(i, k, l, jv) = 0.
179     END IF
180     END DO
181     END DO
182     END DO
183    
184     101 CONTINUE
185    
186     ! boucle sur les niveaux intercouches de 1 a NIV-1
187     ! (flux nul au sommet L=0 et a la base L=NIV)
188    
189     ! calculate flux and moments between adjacent boxes
190     ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
191     ! 1- create temporary moments/masses for partial boxes in transit
192     ! 2- reajusts moments remaining in the box
193    
194     DO l = 1, niv - 1
195     lp = l + 1
196    
197     DO i = 1, lon
198    
199     IF (wgri(i,k,l)<0.) THEN
200     fm(i, l) = -wgri(i, k, l)*dtz
201     alf(i) = fm(i, l)/sm(i, k, lp)
202     sm(i, k, lp) = sm(i, k, lp) - fm(i, l)
203     ELSE
204     fm(i, l) = wgri(i, k, l)*dtz
205     alf(i) = fm(i, l)/sm(i, k, l)
206     sm(i, k, l) = sm(i, k, l) - fm(i, l)
207     END IF
208    
209     alfq(i) = alf(i)*alf(i)
210     alf1(i) = 1. - alf(i)
211     alf1q(i) = alf1(i)*alf1(i)
212     alf2(i) = alf1(i) - alf(i)
213     alf3(i) = alf(i)*alfq(i)
214     alf4(i) = alf1(i)*alf1q(i)
215    
216     END DO
217    
218     DO jv = 1, ntra
219     DO i = 1, lon
220    
221     IF (wgri(i,k,l)<0.) THEN
222    
223     f0(i, l, jv) = alf(i)*(s0(i,k,lp,jv)-alf1(i)*(sz(i,k,lp, &
224     jv)-alf2(i)*szz(i,k,lp,jv)))
225     fz(i, l, jv) = alfq(i)*(sz(i,k,lp,jv)-3.*alf1(i)*szz(i,k,lp,jv))
226     fzz(i, l, jv) = alf3(i)*szz(i, k, lp, jv)
227     fxz(i, l, jv) = alfq(i)*ssxz(i, k, lp, jv)
228     fyz(i, l, jv) = alfq(i)*syz(i, k, lp, jv)
229     fx(i, l, jv) = alf(i)*(ssx(i,k,lp,jv)-alf1(i)*ssxz(i,k,lp,jv))
230     fy(i, l, jv) = alf(i)*(sy(i,k,lp,jv)-alf1(i)*syz(i,k,lp,jv))
231     fxx(i, l, jv) = alf(i)*ssxx(i, k, lp, jv)
232     fxy(i, l, jv) = alf(i)*ssxy(i, k, lp, jv)
233     fyy(i, l, jv) = alf(i)*syy(i, k, lp, jv)
234    
235     s0(i, k, lp, jv) = s0(i, k, lp, jv) - f0(i, l, jv)
236     sz(i, k, lp, jv) = alf1q(i)*(sz(i,k,lp,jv)+3.*alf(i)*szz(i,k,lp, &
237     jv))
238     szz(i, k, lp, jv) = alf4(i)*szz(i, k, lp, jv)
239     ssxz(i, k, lp, jv) = alf1q(i)*ssxz(i, k, lp, jv)
240     syz(i, k, lp, jv) = alf1q(i)*syz(i, k, lp, jv)
241     ssx(i, k, lp, jv) = ssx(i, k, lp, jv) - fx(i, l, jv)
242     sy(i, k, lp, jv) = sy(i, k, lp, jv) - fy(i, l, jv)
243     ssxx(i, k, lp, jv) = ssxx(i, k, lp, jv) - fxx(i, l, jv)
244     ssxy(i, k, lp, jv) = ssxy(i, k, lp, jv) - fxy(i, l, jv)
245     syy(i, k, lp, jv) = syy(i, k, lp, jv) - fyy(i, l, jv)
246    
247     ELSE
248    
249     f0(i, l, jv) = alf(i)*(s0(i,k,l,jv)+alf1(i)*(sz(i,k,l, &
250     jv)+alf2(i)*szz(i,k,l,jv)))
251     fz(i, l, jv) = alfq(i)*(sz(i,k,l,jv)+3.*alf1(i)*szz(i,k,l,jv))
252     fzz(i, l, jv) = alf3(i)*szz(i, k, l, jv)
253     fxz(i, l, jv) = alfq(i)*ssxz(i, k, l, jv)
254     fyz(i, l, jv) = alfq(i)*syz(i, k, l, jv)
255     fx(i, l, jv) = alf(i)*(ssx(i,k,l,jv)+alf1(i)*ssxz(i,k,l,jv))
256     fy(i, l, jv) = alf(i)*(sy(i,k,l,jv)+alf1(i)*syz(i,k,l,jv))
257     fxx(i, l, jv) = alf(i)*ssxx(i, k, l, jv)
258     fxy(i, l, jv) = alf(i)*ssxy(i, k, l, jv)
259     fyy(i, l, jv) = alf(i)*syy(i, k, l, jv)
260    
261     s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, l, jv)
262     sz(i, k, l, jv) = alf1q(i)*(sz(i,k,l,jv)-3.*alf(i)*szz(i,k,l,jv))
263     szz(i, k, l, jv) = alf4(i)*szz(i, k, l, jv)
264     ssxz(i, k, l, jv) = alf1q(i)*ssxz(i, k, l, jv)
265     syz(i, k, l, jv) = alf1q(i)*syz(i, k, l, jv)
266     ssx(i, k, l, jv) = ssx(i, k, l, jv) - fx(i, l, jv)
267     sy(i, k, l, jv) = sy(i, k, l, jv) - fy(i, l, jv)
268     ssxx(i, k, l, jv) = ssxx(i, k, l, jv) - fxx(i, l, jv)
269     ssxy(i, k, l, jv) = ssxy(i, k, l, jv) - fxy(i, l, jv)
270     syy(i, k, l, jv) = syy(i, k, l, jv) - fyy(i, l, jv)
271    
272     END IF
273    
274     END DO
275     END DO
276    
277     END DO
278    
279     ! puts the temporary moments Fi into appropriate neighboring boxes
280    
281     DO l = 1, niv - 1
282     lp = l + 1
283    
284     DO i = 1, lon
285    
286     IF (wgri(i,k,l)<0.) THEN
287     sm(i, k, l) = sm(i, k, l) + fm(i, l)
288     alf(i) = fm(i, l)/sm(i, k, l)
289     ELSE
290     sm(i, k, lp) = sm(i, k, lp) + fm(i, l)
291     alf(i) = fm(i, l)/sm(i, k, lp)
292     END IF
293    
294     alf1(i) = 1. - alf(i)
295     alfq(i) = alf(i)*alf(i)
296     alf1q(i) = alf1(i)*alf1(i)
297     alf2(i) = alf(i)*alf1(i)
298     alf3(i) = alf1(i) - alf(i)
299    
300     END DO
301    
302     DO jv = 1, ntra
303     DO i = 1, lon
304    
305     IF (wgri(i,k,l)<0.) THEN
306    
307     temptm = -alf(i)*s0(i, k, l, jv) + alf1(i)*f0(i, l, jv)
308     s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, l, jv)
309     szz(i, k, l, jv) = alfq(i)*fzz(i, l, jv) + &
310     alf1q(i)*szz(i, k, l, jv) + 5.*(alf2(i)*(fz(i,l,jv)-sz(i,k,l, &
311     jv))+alf3(i)*temptm)
312     sz(i, k, l, jv) = alf(i)*fz(i, l, jv) + alf1(i)*sz(i, k, l, jv) + &
313     3.*temptm
314     ssxz(i, k, l, jv) = alf(i)*fxz(i, l, jv) + &
315     alf1(i)*ssxz(i, k, l, jv) + 3.*(alf1(i)*fx(i,l,jv)-alf(i)*ssx(i &
316     ,k,l,jv))
317     syz(i, k, l, jv) = alf(i)*fyz(i, l, jv) + &
318     alf1(i)*syz(i, k, l, jv) + 3.*(alf1(i)*fy(i,l,jv)-alf(i)*sy(i,k &
319     ,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    
326     ELSE
327    
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) + &
331     alf1q(i)*szz(i, k, lp, jv) + 5.*(alf2(i)*(sz(i,k,lp,jv)-fz(i,l, &
332     jv))-alf3(i)*temptm)
333     sz(i, k, lp, jv) = alf(i)*fz(i, l, jv) + &
334     alf1(i)*sz(i, k, lp, jv) + 3.*temptm
335     ssxz(i, k, lp, jv) = alf(i)*fxz(i, l, jv) + &
336     alf1(i)*ssxz(i, k, lp, jv) + 3.*(alf(i)*ssx(i,k,lp,jv)-alf1(i)* &
337     fx(i,l,jv))
338     syz(i, k, lp, jv) = alf(i)*fyz(i, l, jv) + &
339     alf1(i)*syz(i, k, lp, jv) + 3.*(alf(i)*sy(i,k,lp,jv)-alf1(i)*fy &
340     (i,l,jv))
341     ssx(i, k, lp, jv) = ssx(i, k, lp, jv) + fx(i, l, jv)
342     sy(i, k, lp, jv) = sy(i, k, lp, jv) + fy(i, l, jv)
343     ssxx(i, k, lp, jv) = ssxx(i, k, lp, jv) + fxx(i, l, jv)
344     ssxy(i, k, lp, jv) = ssxy(i, k, lp, jv) + fxy(i, l, jv)
345     syy(i, k, lp, jv) = syy(i, k, lp, jv) + fyy(i, l, jv)
346    
347     END IF
348    
349     END DO
350     END DO
351    
352     END DO
353    
354     ! fin de la boucle principale sur les latitudes
355    
356     END DO
357    
358     DO l = 1, llm
359     DO j = 1, jjp1
360     sm(iip1, j, l) = sm(1, j, l)
361     s0(iip1, j, l, ntra) = s0(1, j, l, ntra)
362     ssx(iip1, j, l, ntra) = ssx(1, j, l, ntra)
363     sy(iip1, j, l, ntra) = sy(1, j, l, ntra)
364     sz(iip1, j, l, ntra) = sz(1, j, l, ntra)
365     END DO
366     END DO
367     ! C-------------------------------------------------------------
368     ! *** Test : diag de la qqtite totale de tarceur
369     ! dans l'atmosphere avant l'advection en z
370     DO l = 1, llm
371     DO j = 1, jjp1
372     DO i = 1, iim
373     sqf = sqf + s0(i, j, l, ntra)
374     END DO
375     END DO
376     END DO
377     PRINT *, '-------- DIAG DANS ADVZ - SORTIE ---------'
378     PRINT *, 'sqf=', sqf
379    
380     RETURN
381     END SUBROUTINE advzp

  ViewVC Help
Powered by ViewVC 1.1.21