My Project
 All Classes Files Functions Variables Macros
advz.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE advz(limit,dtz,w,sm,s0,sx,sy,sz)
5  IMPLICIT NONE
6 
7 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
8 C C
9 C first-order moments (FOM) advection of tracer in Z direction C
10 C C
11 C Source : Pascal Simon (Meteo,CNRM) C
12 C Adaptation : A.Armengaud (LGGE) juin 94 C
13 C C
14 C C
15 C sont des arguments d'entree pour le s-pg... C
16 C C
17 C dq est l'argument de sortie pour le s-pg C
18 C C
19 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
20 C
21 C parametres principaux du modele
22 C
23 #include "dimensions.h"
24 #include "paramet.h"
25 #include "comconst.h"
26 #include "comvert.h"
27 
28 C #include "traceur.h"
29 
30 C Arguments :
31 C -----------
32 C dtz : frequence fictive d'appel du transport
33 C w : flux de masse en z en Pa.m2.s-1
34 
35  INTEGER ntra
36  parameter(ntra = 1)
37 
38  REAL dtz
39  REAL w ( iip1,jjp1,llm )
40 
41 C moments: SM total mass in each grid box
42 C S0 mass of tracer in each grid box
43 C Si 1rst order moment in i direction
44 C
45  REAL sm(iip1,jjp1,llm)
46  + ,s0(iip1,jjp1,llm,ntra)
47  REAL sx(iip1,jjp1,llm,ntra)
48  + ,sy(iip1,jjp1,llm,ntra)
49  + ,sz(iip1,jjp1,llm,ntra)
50 
51 
52 C Local :
53 C -------
54 
55 C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
56 C mass fluxes in kg
57 C declaration :
58 
59  REAL wgri(iip1,jjp1,0:llm)
60 
61 C
62 C the moments F are used as temporary storage for
63 C portions of grid boxes in transit at the current latitude
64 C
65  REAL fm(iim,llm)
66  REAL f0(iim,llm,ntra),fx(iim,llm,ntra)
67  REAL fy(iim,llm,ntra),fz(iim,llm,ntra)
68 C
69 C work arrays
70 C
71  REAL alf(iim),alf1(iim),alfq(iim),alf1q(iim)
72  REAL temptm ! Just temporal variable
73  REAL sqi,sqf
74 C
75  LOGICAL limit
76  INTEGER lon,lat,niv
77  INTEGER i,j,jv,k,l,lp
78 
79  lon = iim
80  lat = jjp1
81  niv = llm
82 
83 C *** Test de passage d'arguments ******
84 
85 c DO 399 l = 1, llm
86 c DO 399 j = 1, jjp1
87 c DO 399 i = 1, iip1
88 c IF (S0(i,j,l,ntra) .lt. 0. ) THEN
89 c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
90 c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
91 c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
92 c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
93 c PRINT*, 'AIE !! debut ADVZ - pbl arg. passage dans ADVZ'
94 c STOP
95 c ENDIF
96  399 CONTINUE
97 
98 C-----------------------------------------------------------------
99 C *** Test : diag de la qqtite totale de traceur
100 C dans l'atmosphere avant l'advection en z
101  sqi = 0.
102  sqf = 0.
103 
104  DO l = 1,llm
105  DO j = 1,jjp1
106  DO i = 1,iim
107 cIM 240305 sqi = sqi + S0(i,j,l,9)
108  sqi = sqi + s0(i,j,l,ntra)
109  ENDDO
110  ENDDO
111  ENDDO
112  print*,'-------- DIAG DANS ADVZ - ENTREE ---------'
113  print*,'sqi=',sqi
114 
115 C-----------------------------------------------------------------
116 C Interface : adaptation nouveau modele
117 C -------------------------------------
118 C
119 C Conversion du flux de masse en kg.s-1
120 
121  DO 500 l = 1,llm
122  DO 500 j = 1,jjp1
123  DO 500 i = 1,iip1
124 c wgri (i,j,llm+1-l) = w (i,j,l) / g
125  wgri(i,j,llm+1-l) = w(i,j,l)
126 c wgri (i,j,0) = 0. ! a detruire ult.
127 c wgri (i,j,l) = 0.1 ! w (i,j,l)
128 c wgri (i,j,llm) = 0. ! a detruire ult.
129  500 CONTINUE
130  DO j = 1,jjp1
131  DO i = 1,iip1
132  wgri(i,j,0)=0.
133  enddo
134  enddo
135 
136 C-----------------------------------------------------------------
137 
138 C start here
139 C boucle sur les latitudes
140 C
141  DO 1 k=1,lat
142 C
143 C place limits on appropriate moments before transport
144 C (if flux-limiting is to be applied)
145 C
146  IF(.NOT.limit) go to 101
147 C
148  DO 10 jv=1,ntra
149  DO 10 l=1,niv
150  DO 100 i=1,lon
151  sz(i,k,l,jv)=sign(amin1(amax1(s0(i,k,l,jv),0.),
152  + abs(sz(i,k,l,jv))),sz(i,k,l,jv))
153  100 CONTINUE
154  10 CONTINUE
155 C
156  101 CONTINUE
157 C
158 C boucle sur les niveaux intercouches de 1 a NIV-1
159 C (flux nul au sommet L=0 et a la base L=NIV)
160 C
161 C calculate flux and moments between adjacent boxes
162 C (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
163 C 1- create temporary moments/masses for partial boxes in transit
164 C 2- reajusts moments remaining in the box
165 C
166  DO 11 l=1,niv-1
167  lp=l+1
168 C
169  DO 110 i=1,lon
170 C
171  IF(wgri(i,k,l).LT.0.) THEN
172  fm(i,l)=-wgri(i,k,l)*dtz
173  alf(i)=fm(i,l)/sm(i,k,lp)
174  sm(i,k,lp)=sm(i,k,lp)-fm(i,l)
175  ELSE
176  fm(i,l)=wgri(i,k,l)*dtz
177  alf(i)=fm(i,l)/sm(i,k,l)
178  sm(i,k,l)=sm(i,k,l)-fm(i,l)
179  ENDIF
180 C
181  alfq(i)=alf(i)*alf(i)
182  alf1(i)=1.-alf(i)
183  alf1q(i)=alf1(i)*alf1(i)
184 C
185  110 CONTINUE
186 C
187  DO 111 jv=1,ntra
188  DO 1110 i=1,lon
189 C
190  IF(wgri(i,k,l).LT.0.) THEN
191 C
192  f0(i,l,jv)=alf(i)*( s0(i,k,lp,jv)-alf1(i)*sz(i,k,lp,jv) )
193  fz(i,l,jv)=alfq(i)*sz(i,k,lp,jv)
194  fx(i,l,jv)=alf(i)*sx(i,k,lp,jv)
195  fy(i,l,jv)=alf(i)*sy(i,k,lp,jv)
196 C
197  s0(i,k,lp,jv)=s0(i,k,lp,jv)-f0(i,l,jv)
198  sz(i,k,lp,jv)=alf1q(i)*sz(i,k,lp,jv)
199  sx(i,k,lp,jv)=sx(i,k,lp,jv)-fx(i,l,jv)
200  sy(i,k,lp,jv)=sy(i,k,lp,jv)-fy(i,l,jv)
201 C
202  ELSE
203 C
204  f0(i,l,jv)=alf(i)*(s0(i,k,l,jv)+alf1(i)*sz(i,k,l,jv) )
205  fz(i,l,jv)=alfq(i)*sz(i,k,l,jv)
206  fx(i,l,jv)=alf(i)*sx(i,k,l,jv)
207  fy(i,l,jv)=alf(i)*sy(i,k,l,jv)
208 C
209  s0(i,k,l,jv)=s0(i,k,l,jv)-f0(i,l,jv)
210  sz(i,k,l,jv)=alf1q(i)*sz(i,k,l,jv)
211  sx(i,k,l,jv)=sx(i,k,l,jv)-fx(i,l,jv)
212  sy(i,k,l,jv)=sy(i,k,l,jv)-fy(i,l,jv)
213 C
214  ENDIF
215 C
216  1110 CONTINUE
217  111 CONTINUE
218 C
219  11 CONTINUE
220 C
221 C puts the temporary moments Fi into appropriate neighboring boxes
222 C
223  DO 12 l=1,niv-1
224  lp=l+1
225 C
226  DO 120 i=1,lon
227 C
228  IF(wgri(i,k,l).LT.0.) THEN
229  sm(i,k,l)=sm(i,k,l)+fm(i,l)
230  alf(i)=fm(i,l)/sm(i,k,l)
231  ELSE
232  sm(i,k,lp)=sm(i,k,lp)+fm(i,l)
233  alf(i)=fm(i,l)/sm(i,k,lp)
234  ENDIF
235 C
236  alf1(i)=1.-alf(i)
237  alfq(i)=alf(i)*alf(i)
238  alf1q(i)=alf1(i)*alf1(i)
239 C
240  120 CONTINUE
241 C
242  DO 121 jv=1,ntra
243  DO 1210 i=1,lon
244 C
245  IF(wgri(i,k,l).LT.0.) THEN
246 C
247  temptm=-alf(i)*s0(i,k,l,jv)+alf1(i)*f0(i,l,jv)
248  s0(i,k,l,jv)=s0(i,k,l,jv)+f0(i,l,jv)
249  sz(i,k,l,jv)=alf(i)*fz(i,l,jv)+alf1(i)*sz(i,k,l,jv)+3.*temptm
250  sx(i,k,l,jv)=sx(i,k,l,jv)+fx(i,l,jv)
251  sy(i,k,l,jv)=sy(i,k,l,jv)+fy(i,l,jv)
252 C
253  ELSE
254 C
255  temptm=alf(i)*s0(i,k,lp,jv)-alf1(i)*f0(i,l,jv)
256  s0(i,k,lp,jv)=s0(i,k,lp,jv)+f0(i,l,jv)
257  sz(i,k,lp,jv)=alf(i)*fz(i,l,jv)+alf1(i)*sz(i,k,lp,jv)
258  + +3.*temptm
259  sx(i,k,lp,jv)=sx(i,k,lp,jv)+fx(i,l,jv)
260  sy(i,k,lp,jv)=sy(i,k,lp,jv)+fy(i,l,jv)
261 C
262  ENDIF
263 C
264  1210 CONTINUE
265  121 CONTINUE
266 C
267  12 CONTINUE
268 C
269 C fin de la boucle principale sur les latitudes
270 C
271  1 CONTINUE
272 C
273 C-------------------------------------------------------------
274 C
275 C ----------- AA Test en fin de ADVX ------ Controle des S*
276 
277 c DO 9999 l = 1, llm
278 c DO 9999 j = 1, jjp1
279 c DO 9999 i = 1, iip1
280 c IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
281 c PRINT*, '-------------------'
282 c PRINT*, 'En fin de ADVZ'
283 c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
284 c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
285 c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
286 c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
287 c WRITE (*,*) 'On arrete !! - pbl en fin de ADVZ1'
288 c STOP
289 c ENDIF
290  9999 CONTINUE
291 
292 C *** ------------------- bouclage cyclique en X ------------
293 
294 c DO l = 1,llm
295 c DO j = 1,jjp1
296 c SM(iip1,j,l) = SM(1,j,l)
297 c S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
298 C sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
299 c sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
300 c sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
301 c ENDDO
302 c ENDDO
303 
304 C-------------------------------------------------------------
305 C *** Test : diag de la qqtite totale de traceur
306 C dans l'atmosphere avant l'advection en z
307  DO l = 1,llm
308  DO j = 1,jjp1
309  DO i = 1,iim
310 cIM 240305 sqf = sqf + S0(i,j,l,9)
311  sqf = sqf + s0(i,j,l,ntra)
312  ENDDO
313  ENDDO
314  ENDDO
315  print*,'-------- DIAG DANS ADVZ - SORTIE ---------'
316  print*,'sqf=', sqf
317 
318 C-------------------------------------------------------------
319  RETURN
320  END
321 C_______________________________________________________________
322 C_______________________________________________________________