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