My Project
 All Classes Files Functions Variables Macros
advzp.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE advzp(LIMIT,DTZ,W,SM,S0,SSX,SY,SZ
5  . ,ssxx,ssxy,ssxz,syy,syz,szz,ntra )
6 
7  IMPLICIT NONE
8 
9 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
10 C C
11 C second-order moments (SOM) advection of tracer in Z direction C
12 C C
13 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
14 C C
15 C Source : Pascal Simon ( Meteo, CNRM ) C
16 C Adaptation : A.A. (LGGE) C
17 C Derniere Modif : 19/11/95 LAST C
18 C C
19 C sont les arguments d'entree pour le s-pg C
20 C C
21 C argument de sortie du s-pg C
22 C C
23 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
24 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
25 C
26 C Rem : Probleme aux poles il faut reecrire ce cas specifique
27 C Attention au sens de l'indexation
28 C
29 
30 C
31 C parametres principaux du modele
32 C
33 #include "dimensions.h"
34 #include "paramet.h"
35 #include "comconst.h"
36 #include "comvert.h"
37 #include "comgeom.h"
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