My Project
 All Classes Files Functions Variables Macros
advx.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE advx(limit,dtx,pbaru,sm,s0,
5  $ sx,sy,sz,lati,latf)
6  IMPLICIT NONE
7 
8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9 C C
10 C first-order moments (FOM) advection of tracer in X direction C
11 C C
12 C Source : Pascal Simon (Meteo,CNRM) C
13 C Adaptation : A.Armengaud (LGGE) juin 94 C
14 C C
15 C limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C
16 C sont des arguments d'entree pour le s-pg... C
17 C C
18 C sm,s0,sx,sy,sz C
19 C sont les arguments de sortie pour le s-pg C
20 C C
21 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22 C
23 C parametres principaux du modele
24 C
25 #include "dimensions.h"
26 #include "paramet.h"
27 #include "comconst.h"
28 #include "comvert.h"
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  REAL pbaru ( iip1,jjp1,llm )
43 
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  sqi = sqi + s0(i,j,l,ntra)
109  ENDDO
110  ENDDO
111  ENDDO
112  print*,'-------- DIAG DANS ADVX - ENTREE ---------'
113  print*,'sqi=',sqi
114 
115 
116 C Interface : adaptation nouveau modele
117 C -------------------------------------
118 C
119 C ---------------------------------------------------------
120 C Conversion des flux de masses en kg/s
121 C pbaru est en N/s d'ou :
122 C ugri est en kg/s
123 
124  DO 500 l = 1,llm
125  DO 500 j = 1,jjm+1
126  DO 500 i = 1,iip1
127 C ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
128  ugri(i,j,llm+1-l) = pbaru(i,j,l)
129  500 CONTINUE
130 
131 
132 C ---------------------------------------------------------
133 C ---------------------------------------------------------
134 C ---------------------------------------------------------
135 
136 C start here
137 C
138 C boucle principale sur les niveaux et les latitudes
139 C
140  DO 1 l=1,niv
141  DO 1 k=lati,latf
142 C
143 C initialisation
144 C
145 C program assumes periodic boundaries in X
146 C
147  DO 10 i=2,lon
148  smnew(i)=sm(i,k,l)+(ugri(i-1,k,l)-ugri(i,k,l))*dtx
149  10 CONTINUE
150  smnew(1)=sm(1,k,l)+(ugri(lon,k,l)-ugri(1,k,l))*dtx
151 C
152 C modifications for extended polar zones
153 C
154  numk=num(k)
155  lonk=lon/numk
156 C
157  IF(numk.GT.1) THEN
158 C
159  DO 111 i=1,lon
160  tm(i)=0.
161  111 CONTINUE
162  DO 112 jv=1,ntra
163  DO 1120 i=1,lon
164  t0(i,jv)=0.
165  tx(i,jv)=0.
166  ty(i,jv)=0.
167  tz(i,jv)=0.
168  1120 CONTINUE
169  112 CONTINUE
170 C
171  DO 11 i2=1,numk
172 C
173  DO 113 i=1,lonk
174  i3=(i-1)*numk+i2
175  tm(i)=tm(i)+sm(i3,k,l)
176  alf(i)=sm(i3,k,l)/tm(i)
177  alf1(i)=1.-alf(i)
178  113 CONTINUE
179 C
180  DO jv=1,ntra
181  DO i=1,lonk
182  i3=(i-1)*numk+i2
183  temptm=-alf(i)*t0(i,jv)+alf1(i)
184  $ *s0(i3,k,l,jv)
185  t0(i,jv)=t0(i,jv)+s0(i3,k,l,jv)
186  tx(i,jv)=alf(i) *sx(i3,k,l,jv)+
187  $ alf1(i)*tx(i,jv) +3.*temptm
188  ty(i,jv)=ty(i,jv)+sy(i3,k,l,jv)
189  tz(i,jv)=tz(i,jv)+sz(i3,k,l,jv)
190  ENDDO
191  ENDDO
192 C
193  11 CONTINUE
194 C
195  ELSE
196 C
197  DO 115 i=1,lon
198  tm(i)=sm(i,k,l)
199  115 CONTINUE
200  DO 116 jv=1,ntra
201  DO 1160 i=1,lon
202  t0(i,jv)=s0(i,k,l,jv)
203  tx(i,jv)=sx(i,k,l,jv)
204  ty(i,jv)=sy(i,k,l,jv)
205  tz(i,jv)=sz(i,k,l,jv)
206  1160 CONTINUE
207  116 CONTINUE
208 C
209  ENDIF
210 C
211  DO 117 i=1,lonk
212  uext(i)=ugri(i*numk,k,l)
213  117 CONTINUE
214 C
215 C place limits on appropriate moments before transport
216 C (if flux-limiting is to be applied)
217 C
218  IF(.NOT.limit) go to 13
219 C
220  DO 12 jv=1,ntra
221  DO 120 i=1,lonk
222  tx(i,jv)=sign(amin1(amax1(t0(i,jv),0.),abs(tx(i,jv))),tx(i,jv))
223  120 CONTINUE
224  12 CONTINUE
225 C
226  13 CONTINUE
227 C
228 C calculate flux and moments between adjacent boxes
229 C 1- create temporary moments/masses for partial boxes in transit
230 C 2- reajusts moments remaining in the box
231 C
232 C flux from IP to I if U(I).lt.0
233 C
234  DO 140 i=1,lonk-1
235  IF(uext(i).LT.0.) THEN
236  fm(i)=-uext(i)*dtx
237  alf(i)=fm(i)/tm(i+1)
238  tm(i+1)=tm(i+1)-fm(i)
239  ENDIF
240  140 CONTINUE
241 C
242  i=lonk
243  IF(uext(i).LT.0.) THEN
244  fm(i)=-uext(i)*dtx
245  alf(i)=fm(i)/tm(1)
246  tm(1)=tm(1)-fm(i)
247  ENDIF
248 C
249 C flux from I to IP if U(I).gt.0
250 C
251  DO 141 i=1,lonk
252  IF(uext(i).GE.0.) THEN
253  fm(i)=uext(i)*dtx
254  alf(i)=fm(i)/tm(i)
255  tm(i)=tm(i)-fm(i)
256  ENDIF
257  141 CONTINUE
258 C
259  DO 142 i=1,lonk
260  alfq(i)=alf(i)*alf(i)
261  alf1(i)=1.-alf(i)
262  alf1q(i)=alf1(i)*alf1(i)
263  142 CONTINUE
264 C
265  DO 150 jv=1,ntra
266  DO 1500 i=1,lonk-1
267 C
268  IF(uext(i).LT.0.) THEN
269 C
270  f0(i,jv)=alf(i)* ( t0(i+1,jv)-alf1(i)*tx(i+1,jv) )
271  fx(i,jv)=alfq(i)*tx(i+1,jv)
272  fy(i,jv)=alf(i)*ty(i+1,jv)
273  fz(i,jv)=alf(i)*tz(i+1,jv)
274 C
275  t0(i+1,jv)=t0(i+1,jv)-f0(i,jv)
276  tx(i+1,jv)=alf1q(i)*tx(i+1,jv)
277  ty(i+1,jv)=ty(i+1,jv)-fy(i,jv)
278  tz(i+1,jv)=tz(i+1,jv)-fz(i,jv)
279 C
280  ENDIF
281 C
282  1500 CONTINUE
283  150 CONTINUE
284 C
285  i=lonk
286  IF(uext(i).LT.0.) THEN
287 C
288  DO 151 jv=1,ntra
289 C
290  f0(i,jv)=alf(i)* ( t0(1,jv)-alf1(i)*tx(1,jv) )
291  fx(i,jv)=alfq(i)*tx(1,jv)
292  fy(i,jv)=alf(i)*ty(1,jv)
293  fz(i,jv)=alf(i)*tz(1,jv)
294 C
295  t0(1,jv)=t0(1,jv)-f0(i,jv)
296  tx(1,jv)=alf1q(i)*tx(1,jv)
297  ty(1,jv)=ty(1,jv)-fy(i,jv)
298  tz(1,jv)=tz(1,jv)-fz(i,jv)
299 C
300  151 CONTINUE
301 C
302  ENDIF
303 C
304  DO 152 jv=1,ntra
305  DO 1520 i=1,lonk
306 C
307  IF(uext(i).GE.0.) THEN
308 C
309  f0(i,jv)=alf(i)* ( t0(i,jv)+alf1(i)*tx(i,jv) )
310  fx(i,jv)=alfq(i)*tx(i,jv)
311  fy(i,jv)=alf(i)*ty(i,jv)
312  fz(i,jv)=alf(i)*tz(i,jv)
313 C
314  t0(i,jv)=t0(i,jv)-f0(i,jv)
315  tx(i,jv)=alf1q(i)*tx(i,jv)
316  ty(i,jv)=ty(i,jv)-fy(i,jv)
317  tz(i,jv)=tz(i,jv)-fz(i,jv)
318 C
319  ENDIF
320 C
321  1520 CONTINUE
322  152 CONTINUE
323 C
324 C puts the temporary moments Fi into appropriate neighboring boxes
325 C
326  DO 160 i=1,lonk
327  IF(uext(i).LT.0.) THEN
328  tm(i)=tm(i)+fm(i)
329  alf(i)=fm(i)/tm(i)
330  ENDIF
331  160 CONTINUE
332 C
333  DO 161 i=1,lonk-1
334  IF(uext(i).GE.0.) THEN
335  tm(i+1)=tm(i+1)+fm(i)
336  alf(i)=fm(i)/tm(i+1)
337  ENDIF
338  161 CONTINUE
339 C
340  i=lonk
341  IF(uext(i).GE.0.) THEN
342  tm(1)=tm(1)+fm(i)
343  alf(i)=fm(i)/tm(1)
344  ENDIF
345 C
346  DO 162 i=1,lonk
347  alf1(i)=1.-alf(i)
348  162 CONTINUE
349 C
350  DO 170 jv=1,ntra
351  DO 1700 i=1,lonk
352 C
353  IF(uext(i).LT.0.) THEN
354 C
355  temptm=-alf(i)*t0(i,jv)+alf1(i)*f0(i,jv)
356  t0(i,jv)=t0(i,jv)+f0(i,jv)
357  tx(i,jv)=alf(i)*fx(i,jv)+alf1(i)*tx(i,jv)+3.*temptm
358  ty(i,jv)=ty(i,jv)+fy(i,jv)
359  tz(i,jv)=tz(i,jv)+fz(i,jv)
360 C
361  ENDIF
362 C
363  1700 CONTINUE
364  170 CONTINUE
365 C
366  DO 171 jv=1,ntra
367  DO 1710 i=1,lonk-1
368 C
369  IF(uext(i).GE.0.) THEN
370 C
371  temptm=alf(i)*t0(i+1,jv)-alf1(i)*f0(i,jv)
372  t0(i+1,jv)=t0(i+1,jv)+f0(i,jv)
373  tx(i+1,jv)=alf(i)*fx(i,jv)+alf1(i)*tx(i+1,jv)+3.*temptm
374  ty(i+1,jv)=ty(i+1,jv)+fy(i,jv)
375  tz(i+1,jv)=tz(i+1,jv)+fz(i,jv)
376 C
377  ENDIF
378 C
379  1710 CONTINUE
380  171 CONTINUE
381 C
382  i=lonk
383  IF(uext(i).GE.0.) THEN
384  DO 172 jv=1,ntra
385  temptm=alf(i)*t0(1,jv)-alf1(i)*f0(i,jv)
386  t0(1,jv)=t0(1,jv)+f0(i,jv)
387  tx(1,jv)=alf(i)*fx(i,jv)+alf1(i)*tx(1,jv)+3.*temptm
388  ty(1,jv)=ty(1,jv)+fy(i,jv)
389  tz(1,jv)=tz(1,jv)+fz(i,jv)
390  172 CONTINUE
391  ENDIF
392 C
393 C retour aux mailles d'origine (passage des Tij aux Sij)
394 C
395  IF(numk.GT.1) THEN
396 C
397  DO 180 i2=1,numk
398 C
399  DO 180 i=1,lonk
400 C
401  i3=i2+(i-1)*numk
402  sm(i3,k,l)=smnew(i3)
403  alf(i)=smnew(i3)/tm(i)
404  tm(i)=tm(i)-smnew(i3)
405 C
406  alfq(i)=alf(i)*alf(i)
407  alf1(i)=1.-alf(i)
408  alf1q(i)=alf1(i)*alf1(i)
409 C
410  180 CONTINUE
411 C
412  DO jv=1,ntra
413  DO i=1,lonk
414 C
415  i3=i2+(i-1)*numk
416  s0(i3,k,l,jv)=alf(i)
417  $ * (t0(i,jv)-alf1(i)*tx(i,jv))
418  sx(i3,k,l,jv)=alfq(i)*tx(i,jv)
419  sy(i3,k,l,jv)=alf(i)*ty(i,jv)
420  sz(i3,k,l,jv)=alf(i)*tz(i,jv)
421 C
422 C reajusts moments remaining in the box
423 C
424  t0(i,jv)=t0(i,jv)-s0(i3,k,l,jv)
425  tx(i,jv)=alf1q(i)*tx(i,jv)
426  ty(i,jv)=ty(i,jv)-sy(i3,k,l,jv)
427  tz(i,jv)=tz(i,jv)-sz(i3,k,l,jv)
428  ENDDO
429  ENDDO
430 C
431 C
432  ELSE
433 C
434  DO 190 i=1,lon
435  sm(i,k,l)=tm(i)
436  190 CONTINUE
437  DO 191 jv=1,ntra
438  DO 1910 i=1,lon
439  s0(i,k,l,jv)=t0(i,jv)
440  sx(i,k,l,jv)=tx(i,jv)
441  sy(i,k,l,jv)=ty(i,jv)
442  sz(i,k,l,jv)=tz(i,jv)
443  1910 CONTINUE
444  191 CONTINUE
445 C
446  ENDIF
447 C
448  1 CONTINUE
449 C
450 C ----------- AA Test en fin de ADVX ------ Controle des S*
451 c OK
452 c DO 9998 l = 1, llm
453 c DO 9998 j = 1, jjp1
454 c DO 9998 i = 1, iip1
455 c IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
456 c PRINT*, '-------------------'
457 c PRINT*, 'En fin de ADVX'
458 c PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
459 c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
460 c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
461 c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
462 c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
463 c WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
464 cc STOP
465 c ENDIF
466 c 9998 CONTINUE
467 c
468 C ---------- bouclage cyclique
469  DO itrac=1,ntra
470  DO l = 1,llm
471  DO j = lati,latf
472  sm(iip1,j,l) = sm(1,j,l)
473  s0(iip1,j,l,itrac) = s0(1,j,l,itrac)
474  sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
475  sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
476  sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
477  END DO
478  END DO
479  ENDDO
480 
481 c ----------- qqtite totale de traceur dans tte l'atmosphere
482  DO l = 1, llm
483  DO j = 1, jjp1
484  DO i = 1, iim
485  sqf = sqf + s0(i,j,l,ntra)
486  END DO
487  END DO
488  END DO
489 c
490  print*,'------ DIAG DANS ADVX - SORTIE -----'
491  print*,'sqf=',sqf
492 c-------------
493 
494  RETURN
495  END
496 C_________________________________________________________________
497 C_________________________________________________________________