My Project
 All Classes Files Functions Variables Macros
advyp.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE advyp(LIMIT,DTY,PBARV,SM,S0,SSX,SY,SZ
5  . ,ssxx,ssxy,ssxz,syy,syz,szz,ntra )
6  IMPLICIT NONE
7 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
8 C C
9 C second-order moments (SOM) advection of tracer in Y direction C
10 C C
11 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
12 C C
13 C Source : Pascal Simon ( Meteo, CNRM ) C
14 C Adaptation : A.A. (LGGE) C
15 C Derniere Modif : 19/10/95 LAST
16 C C
17 C sont les arguments d'entree pour le s-pg C
18 C C
19 C argument de sortie du s-pg C
20 C C
21 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
23 C
24 C Rem : Probleme aux poles il faut reecrire ce cas specifique
25 C Attention au sens de l'indexation
26 C
27 C parametres principaux du modele
28 C
29 C
30 #include "dimensions.h"
31 #include "paramet.h"
32 #include "comconst.h"
33 #include "comvert.h"
34 #include "comgeom.h"
35 
36 C Arguments :
37 C ----------
38 C dty : frequence fictive d'appel du transport
39 C parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
40 
41  INTEGER lon,lat,niv
42  INTEGER i,j,jv,k,kp,l
43  INTEGER ntra
44 C PARAMETER (ntra = 1)
45 
46  REAL dty
47  REAL pbarv ( iip1,jjm, llm )
48 
49 C moments: SM total mass in each grid box
50 C S0 mass of tracer in each grid box
51 C Si 1rst order moment in i direction
52 C
53  REAL sm(iip1,jjp1,llm)
54  + ,s0(iip1,jjp1,llm,ntra)
55  REAL ssx(iip1,jjp1,llm,ntra)
56  + ,sy(iip1,jjp1,llm,ntra)
57  + ,sz(iip1,jjp1,llm,ntra)
58  + ,ssxx(iip1,jjp1,llm,ntra)
59  + ,ssxy(iip1,jjp1,llm,ntra)
60  + ,ssxz(iip1,jjp1,llm,ntra)
61  + ,syy(iip1,jjp1,llm,ntra)
62  + ,syz(iip1,jjp1,llm,ntra)
63  + ,szz(iip1,jjp1,llm,ntra)
64 C
65 C Local :
66 C -------
67 
68 C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
69 C mass fluxes in kg
70 C declaration :
71 
72  REAL vgri(iip1,0:jjp1,llm)
73 
74 C Rem : UGRI et WGRI ne sont pas utilises dans
75 C cette subroutine ( advection en y uniquement )
76 C Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
77 C
78 C the moments F are similarly defined and used as temporary
79 C storage for portions of the grid boxes in transit
80 C
81 C the moments Fij are used as temporary storage for
82 C portions of the grid boxes in transit at the current level
83 C
84 C work arrays
85 C
86 C
87  REAL f0(iim,0:jjp1,ntra),fm(iim,0:jjp1)
88  REAL fx(iim,jjm,ntra),fy(iim,jjm,ntra)
89  REAL fz(iim,jjm,ntra)
90  REAL fxx(iim,jjm,ntra),fxy(iim,jjm,ntra)
91  REAL fxz(iim,jjm,ntra),fyy(iim,jjm,ntra)
92  REAL fyz(iim,jjm,ntra),fzz(iim,jjm,ntra)
93  REAL s00(ntra)
94  REAL sm0 ! Just temporal variable
95 C
96 C work arrays
97 C
98  REAL alf(iim,0:jjp1),alf1(iim,0:jjp1)
99  REAL alfq(iim,0:jjp1),alf1q(iim,0:jjp1)
100  REAL alf2(iim,0:jjp1),alf3(iim,0:jjp1)
101  REAL alf4(iim,0:jjp1)
102  REAL temptm ! Just temporal variable
103  REAL slpmax,s1max,s1new,s2new
104 c
105 C Special pour poles
106 c
107  REAL sbms,sfms,sfzs,sbmn,sfmn,sfzn
108  REAL sns0(ntra),snsz(ntra),snsm
109  REAL qy1(iim,llm,ntra),qylat(iim,llm,ntra)
110  REAL cx1(llm,ntra), cxlat(llm,ntra)
111  REAL cy1(llm,ntra), cylat(llm,ntra)
112  REAL z1(iim), zcos(iim), zsin(iim)
113  REAL ssum
114  EXTERNAL ssum
115 C
116  REAL sqi,sqf
117  LOGICAL limit
118 
119  lon = iim ! rem : Il est possible qu'un pbl. arrive ici
120  lat = jjp1 ! a cause des dim. differentes entre les
121  niv = llm ! tab. S et VGRI
122 
123 c-----------------------------------------------------------------
124 C initialisations
125 
126  sbms = 0.
127  sfms = 0.
128  sfzs = 0.
129  sbmn = 0.
130  sfmn = 0.
131  sfzn = 0.
132 
133 c-----------------------------------------------------------------
134 C *** Test : diag de la qtite totale de traceur dans
135 C l'atmosphere avant l'advection en Y
136 c
137  sqi = 0.
138  sqf = 0.
139 
140  DO l = 1,llm
141  DO j = 1,jjp1
142  DO i = 1,iim
143  sqi = sqi + s0(i,j,l,ntra)
144  END DO
145  END DO
146  END DO
147  print*,'---------- DIAG DANS ADVY - ENTREE --------'
148  print*,'sqi=',sqi
149 
150 c-----------------------------------------------------------------
151 C Interface : adaptation nouveau modele
152 C -------------------------------------
153 C
154 C Conversion des flux de masses en kg
155 C-AA 20/10/94 le signe -1 est necessaire car indexation opposee
156 
157  DO 500 l = 1,llm
158  DO 500 j = 1,jjm
159  DO 500 i = 1,iip1
160  vgri(i,j,llm+1-l)=-1.*pbarv(i,j,l)
161  500 CONTINUE
162 
163 CAA Initialisation de flux fictifs aux bords sup. des boites pol.
164 
165  DO l = 1,llm
166  DO i = 1,iip1
167  vgri(i,0,l) = 0.
168  vgri(i,jjp1,l) = 0.
169  ENDDO
170  ENDDO
171 c
172 c----------------- START HERE -----------------------
173 C boucle sur les niveaux
174 C
175  DO 1 l=1,niv
176 C
177 C place limits on appropriate moments before transport
178 C (if flux-limiting is to be applied)
179 C
180  IF(.NOT.limit) go to 11
181 C
182  DO 10 jv=1,ntra
183  DO 10 k=1,lat
184  DO 100 i=1,lon
185  IF(s0(i,k,l,jv).GT.0.) THEN
186  slpmax=amax1(s0(i,k,l,jv),0.)
187  s1max=1.5*slpmax
188  s1new=amin1(s1max,amax1(-s1max,sy(i,k,l,jv)))
189  s2new=amin1( 2.*slpmax-abs(s1new)/3. ,
190  + amax1(abs(s1new)-slpmax,syy(i,k,l,jv)) )
191  sy(i,k,l,jv)=s1new
192  syy(i,k,l,jv)=s2new
193  ssxy(i,k,l,jv)=amin1(slpmax,amax1(-slpmax,ssxy(i,k,l,jv)))
194  syz(i,k,l,jv)=amin1(slpmax,amax1(-slpmax,syz(i,k,l,jv)))
195  ELSE
196  sy(i,k,l,jv)=0.
197  syy(i,k,l,jv)=0.
198  ssxy(i,k,l,jv)=0.
199  syz(i,k,l,jv)=0.
200  ENDIF
201  100 CONTINUE
202  10 CONTINUE
203 C
204  11 CONTINUE
205 C
206 C le flux a travers le pole Nord est traite separement
207 C
208  sm0=0.
209  DO 20 jv=1,ntra
210  s00(jv)=0.
211  20 CONTINUE
212 C
213  DO 21 i=1,lon
214 C
215  IF(vgri(i,0,l).LE.0.) THEN
216  fm(i,0)=-vgri(i,0,l)*dty
217  alf(i,0)=fm(i,0)/sm(i,1,l)
218  sm(i,1,l)=sm(i,1,l)-fm(i,0)
219  sm0=sm0+fm(i,0)
220  ENDIF
221 C
222  alfq(i,0)=alf(i,0)*alf(i,0)
223  alf1(i,0)=1.-alf(i,0)
224  alf1q(i,0)=alf1(i,0)*alf1(i,0)
225  alf2(i,0)=alf1(i,0)-alf(i,0)
226  alf3(i,0)=alf(i,0)*alfq(i,0)
227  alf4(i,0)=alf1(i,0)*alf1q(i,0)
228 C
229  21 CONTINUE
230 c print*,'ADVYP 21'
231 C
232  DO 22 jv=1,ntra
233  DO 220 i=1,lon
234 C
235  IF(vgri(i,0,l).LE.0.) THEN
236 C
237  f0(i,0,jv)=alf(i,0)* ( s0(i,1,l,jv)-alf1(i,0)*
238  + ( sy(i,1,l,jv)-alf2(i,0)*syy(i,1,l,jv) ) )
239 C
240  s00(jv)=s00(jv)+f0(i,0,jv)
241  s0(i,1,l,jv)=s0(i,1,l,jv)-f0(i,0,jv)
242  sy(i,1,l,jv)=alf1q(i,0)*
243  + (sy(i,1,l,jv)+3.*alf(i,0)*syy(i,1,l,jv))
244  syy(i,1,l,jv)=alf4(i,0)*syy(i,1,l,jv)
245  ssx(i,1,l,jv)=alf1(i,0)*
246  + (ssx(i,1,l,jv)+alf(i,0)*ssxy(i,1,l,jv) )
247  sz(i,1,l,jv)=alf1(i,0)*
248  + (sz(i,1,l,jv)+alf(i,0)*ssxz(i,1,l,jv) )
249  ssxx(i,1,l,jv)=alf1(i,0)*ssxx(i,1,l,jv)
250  ssxz(i,1,l,jv)=alf1(i,0)*ssxz(i,1,l,jv)
251  szz(i,1,l,jv)=alf1(i,0)*szz(i,1,l,jv)
252  ssxy(i,1,l,jv)=alf1q(i,0)*ssxy(i,1,l,jv)
253  syz(i,1,l,jv)=alf1q(i,0)*syz(i,1,l,jv)
254 C
255  ENDIF
256 C
257  220 CONTINUE
258  22 CONTINUE
259 C
260  DO 23 i=1,lon
261  IF(vgri(i,0,l).GT.0.) THEN
262  fm(i,0)=vgri(i,0,l)*dty
263  alf(i,0)=fm(i,0)/sm0
264  ENDIF
265  23 CONTINUE
266 C
267  DO 24 jv=1,ntra
268  DO 240 i=1,lon
269  IF(vgri(i,0,l).GT.0.) THEN
270  f0(i,0,jv)=alf(i,0)*s00(jv)
271  ENDIF
272  240 CONTINUE
273  24 CONTINUE
274 C
275 C puts the temporary moments Fi into appropriate neighboring boxes
276 C
277 c print*,'av ADVYP 25'
278  DO 25 i=1,lon
279 C
280  IF(vgri(i,0,l).GT.0.) THEN
281  sm(i,1,l)=sm(i,1,l)+fm(i,0)
282  alf(i,0)=fm(i,0)/sm(i,1,l)
283  ENDIF
284 C
285  alfq(i,0)=alf(i,0)*alf(i,0)
286  alf1(i,0)=1.-alf(i,0)
287  alf1q(i,0)=alf1(i,0)*alf1(i,0)
288  alf2(i,0)=alf1(i,0)-alf(i,0)
289  alf3(i,0)=alf1(i,0)*alf(i,0)
290 C
291  25 CONTINUE
292 c print*,'av ADVYP 25'
293 C
294  DO 26 jv=1,ntra
295  DO 260 i=1,lon
296 C
297  IF(vgri(i,0,l).GT.0.) THEN
298 C
299  temptm=alf(i,0)*s0(i,1,l,jv)-alf1(i,0)*f0(i,0,jv)
300  s0(i,1,l,jv)=s0(i,1,l,jv)+f0(i,0,jv)
301  syy(i,1,l,jv)=alf1q(i,0)*syy(i,1,l,jv)
302  + +5.*( alf3(i,0)*sy(i,1,l,jv)-alf2(i,0)*temptm )
303  sy(i,1,l,jv)=alf1(i,0)*sy(i,1,l,jv)+3.*temptm
304  ssxy(i,1,l,jv)=alf1(i,0)*ssxy(i,1,l,jv)+3.*alf(i,0)*ssx(i,1,l,jv)
305  syz(i,1,l,jv)=alf1(i,0)*syz(i,1,l,jv)+3.*alf(i,0)*sz(i,1,l,jv)
306 C
307  ENDIF
308 C
309  260 CONTINUE
310  26 CONTINUE
311 C
312 C calculate flux and moments between adjacent boxes
313 C 1- create temporary moments/masses for partial boxes in transit
314 C 2- reajusts moments remaining in the box
315 C
316 C flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
317 C
318 c print*,'av ADVYP 30'
319  DO 30 k=1,lat-1
320  kp=k+1
321  DO 300 i=1,lon
322 C
323  IF(vgri(i,k,l).LT.0.) THEN
324  fm(i,k)=-vgri(i,k,l)*dty
325  alf(i,k)=fm(i,k)/sm(i,kp,l)
326  sm(i,kp,l)=sm(i,kp,l)-fm(i,k)
327  ELSE
328  fm(i,k)=vgri(i,k,l)*dty
329  alf(i,k)=fm(i,k)/sm(i,k,l)
330  sm(i,k,l)=sm(i,k,l)-fm(i,k)
331  ENDIF
332 C
333  alfq(i,k)=alf(i,k)*alf(i,k)
334  alf1(i,k)=1.-alf(i,k)
335  alf1q(i,k)=alf1(i,k)*alf1(i,k)
336  alf2(i,k)=alf1(i,k)-alf(i,k)
337  alf3(i,k)=alf(i,k)*alfq(i,k)
338  alf4(i,k)=alf1(i,k)*alf1q(i,k)
339 C
340  300 CONTINUE
341  30 CONTINUE
342 c print*,'ap ADVYP 30'
343 C
344  DO 31 jv=1,ntra
345  DO 31 k=1,lat-1
346  kp=k+1
347  DO 310 i=1,lon
348 C
349  IF(vgri(i,k,l).LT.0.) THEN
350 C
351  f0(i,k,jv)=alf(i,k)* ( s0(i,kp,l,jv)-alf1(i,k)*
352  + ( sy(i,kp,l,jv)-alf2(i,k)*syy(i,kp,l,jv) ) )
353  fy(i,k,jv)=alfq(i,k)*
354  + (sy(i,kp,l,jv)-3.*alf1(i,k)*syy(i,kp,l,jv))
355  fyy(i,k,jv)=alf3(i,k)*syy(i,kp,l,jv)
356  fx(i,k,jv)=alf(i,k)*
357  + (ssx(i,kp,l,jv)-alf1(i,k)*ssxy(i,kp,l,jv))
358  fz(i,k,jv)=alf(i,k)*
359  + (sz(i,kp,l,jv)-alf1(i,k)*syz(i,kp,l,jv))
360  fxy(i,k,jv)=alfq(i,k)*ssxy(i,kp,l,jv)
361  fyz(i,k,jv)=alfq(i,k)*syz(i,kp,l,jv)
362  fxx(i,k,jv)=alf(i,k)*ssxx(i,kp,l,jv)
363  fxz(i,k,jv)=alf(i,k)*ssxz(i,kp,l,jv)
364  fzz(i,k,jv)=alf(i,k)*szz(i,kp,l,jv)
365 C
366  s0(i,kp,l,jv)=s0(i,kp,l,jv)-f0(i,k,jv)
367  sy(i,kp,l,jv)=alf1q(i,k)*
368  + (sy(i,kp,l,jv)+3.*alf(i,k)*syy(i,kp,l,jv))
369  syy(i,kp,l,jv)=alf4(i,k)*syy(i,kp,l,jv)
370  ssx(i,kp,l,jv)=ssx(i,kp,l,jv)-fx(i,k,jv)
371  sz(i,kp,l,jv)=sz(i,kp,l,jv)-fz(i,k,jv)
372  ssxx(i,kp,l,jv)=ssxx(i,kp,l,jv)-fxx(i,k,jv)
373  ssxz(i,kp,l,jv)=ssxz(i,kp,l,jv)-fxz(i,k,jv)
374  szz(i,kp,l,jv)=szz(i,kp,l,jv)-fzz(i,k,jv)
375  ssxy(i,kp,l,jv)=alf1q(i,k)*ssxy(i,kp,l,jv)
376  syz(i,kp,l,jv)=alf1q(i,k)*syz(i,kp,l,jv)
377 C
378  ELSE
379 C
380  f0(i,k,jv)=alf(i,k)* ( s0(i,k,l,jv)+alf1(i,k)*
381  + ( sy(i,k,l,jv)+alf2(i,k)*syy(i,k,l,jv) ) )
382  fy(i,k,jv)=alfq(i,k)*
383  + (sy(i,k,l,jv)+3.*alf1(i,k)*syy(i,k,l,jv))
384  fyy(i,k,jv)=alf3(i,k)*syy(i,k,l,jv)
385  fx(i,k,jv)=alf(i,k)*(ssx(i,k,l,jv)+alf1(i,k)*ssxy(i,k,l,jv))
386  fz(i,k,jv)=alf(i,k)*(sz(i,k,l,jv)+alf1(i,k)*syz(i,k,l,jv))
387  fxy(i,k,jv)=alfq(i,k)*ssxy(i,k,l,jv)
388  fyz(i,k,jv)=alfq(i,k)*syz(i,k,l,jv)
389  fxx(i,k,jv)=alf(i,k)*ssxx(i,k,l,jv)
390  fxz(i,k,jv)=alf(i,k)*ssxz(i,k,l,jv)
391  fzz(i,k,jv)=alf(i,k)*szz(i,k,l,jv)
392 C
393  s0(i,k,l,jv)=s0(i,k,l,jv)-f0(i,k,jv)
394  sy(i,k,l,jv)=alf1q(i,k)*
395  + (sy(i,k,l,jv)-3.*alf(i,k)*syy(i,k,l,jv))
396  syy(i,k,l,jv)=alf4(i,k)*syy(i,k,l,jv)
397  ssx(i,k,l,jv)=ssx(i,k,l,jv)-fx(i,k,jv)
398  sz(i,k,l,jv)=sz(i,k,l,jv)-fz(i,k,jv)
399  ssxx(i,k,l,jv)=ssxx(i,k,l,jv)-fxx(i,k,jv)
400  ssxz(i,k,l,jv)=ssxz(i,k,l,jv)-fxz(i,k,jv)
401  szz(i,k,l,jv)=szz(i,k,l,jv)-fzz(i,k,jv)
402  ssxy(i,k,l,jv)=alf1q(i,k)*ssxy(i,k,l,jv)
403  syz(i,k,l,jv)=alf1q(i,k)*syz(i,k,l,jv)
404 C
405  ENDIF
406 C
407  310 CONTINUE
408  31 CONTINUE
409 c print*,'ap ADVYP 31'
410 C
411 C puts the temporary moments Fi into appropriate neighboring boxes
412 C
413  DO 32 k=1,lat-1
414  kp=k+1
415  DO 320 i=1,lon
416 C
417  IF(vgri(i,k,l).LT.0.) THEN
418  sm(i,k,l)=sm(i,k,l)+fm(i,k)
419  alf(i,k)=fm(i,k)/sm(i,k,l)
420  ELSE
421  sm(i,kp,l)=sm(i,kp,l)+fm(i,k)
422  alf(i,k)=fm(i,k)/sm(i,kp,l)
423  ENDIF
424 C
425  alfq(i,k)=alf(i,k)*alf(i,k)
426  alf1(i,k)=1.-alf(i,k)
427  alf1q(i,k)=alf1(i,k)*alf1(i,k)
428  alf2(i,k)=alf1(i,k)-alf(i,k)
429  alf3(i,k)=alf1(i,k)*alf(i,k)
430 C
431  320 CONTINUE
432  32 CONTINUE
433 c print*,'ap ADVYP 32'
434 C
435  DO 33 jv=1,ntra
436  DO 33 k=1,lat-1
437  kp=k+1
438  DO 330 i=1,lon
439 C
440  IF(vgri(i,k,l).LT.0.) THEN
441 C
442  temptm=-alf(i,k)*s0(i,k,l,jv)+alf1(i,k)*f0(i,k,jv)
443  s0(i,k,l,jv)=s0(i,k,l,jv)+f0(i,k,jv)
444  syy(i,k,l,jv)=alfq(i,k)*fyy(i,k,jv)+alf1q(i,k)*syy(i,k,l,jv)
445  + +5.*( alf3(i,k)*(fy(i,k,jv)-sy(i,k,l,jv))+alf2(i,k)*temptm )
446  sy(i,k,l,jv)=alf(i,k)*fy(i,k,jv)+alf1(i,k)*sy(i,k,l,jv)
447  + +3.*temptm
448  ssxy(i,k,l,jv)=alf(i,k)*fxy(i,k,jv)+alf1(i,k)*ssxy(i,k,l,jv)
449  + +3.*(alf1(i,k)*fx(i,k,jv)-alf(i,k)*ssx(i,k,l,jv))
450  syz(i,k,l,jv)=alf(i,k)*fyz(i,k,jv)+alf1(i,k)*syz(i,k,l,jv)
451  + +3.*(alf1(i,k)*fz(i,k,jv)-alf(i,k)*sz(i,k,l,jv))
452  ssx(i,k,l,jv)=ssx(i,k,l,jv)+fx(i,k,jv)
453  sz(i,k,l,jv)=sz(i,k,l,jv)+fz(i,k,jv)
454  ssxx(i,k,l,jv)=ssxx(i,k,l,jv)+fxx(i,k,jv)
455  ssxz(i,k,l,jv)=ssxz(i,k,l,jv)+fxz(i,k,jv)
456  szz(i,k,l,jv)=szz(i,k,l,jv)+fzz(i,k,jv)
457 C
458  ELSE
459 C
460  temptm=alf(i,k)*s0(i,kp,l,jv)-alf1(i,k)*f0(i,k,jv)
461  s0(i,kp,l,jv)=s0(i,kp,l,jv)+f0(i,k,jv)
462  syy(i,kp,l,jv)=alfq(i,k)*fyy(i,k,jv)+alf1q(i,k)*syy(i,kp,l,jv)
463  + +5.*( alf3(i,k)*(sy(i,kp,l,jv)-fy(i,k,jv))-alf2(i,k)*temptm )
464  sy(i,kp,l,jv)=alf(i,k)*fy(i,k,jv)+alf1(i,k)*sy(i,kp,l,jv)
465  + +3.*temptm
466  ssxy(i,kp,l,jv)=alf(i,k)*fxy(i,k,jv)+alf1(i,k)*ssxy(i,kp,l,jv)
467  + +3.*(alf(i,k)*ssx(i,kp,l,jv)-alf1(i,k)*fx(i,k,jv))
468  syz(i,kp,l,jv)=alf(i,k)*fyz(i,k,jv)+alf1(i,k)*syz(i,kp,l,jv)
469  + +3.*(alf(i,k)*sz(i,kp,l,jv)-alf1(i,k)*fz(i,k,jv))
470  ssx(i,kp,l,jv)=ssx(i,kp,l,jv)+fx(i,k,jv)
471  sz(i,kp,l,jv)=sz(i,kp,l,jv)+fz(i,k,jv)
472  ssxx(i,kp,l,jv)=ssxx(i,kp,l,jv)+fxx(i,k,jv)
473  ssxz(i,kp,l,jv)=ssxz(i,kp,l,jv)+fxz(i,k,jv)
474  szz(i,kp,l,jv)=szz(i,kp,l,jv)+fzz(i,k,jv)
475 C
476  ENDIF
477 C
478  330 CONTINUE
479  33 CONTINUE
480 c print*,'ap ADVYP 33'
481 C
482 C traitement special pour le pole Sud (idem pole Nord)
483 C
484  k=lat
485 C
486  sm0=0.
487  DO 40 jv=1,ntra
488  s00(jv)=0.
489  40 CONTINUE
490 C
491  DO 41 i=1,lon
492 C
493  IF(vgri(i,k,l).GE.0.) THEN
494  fm(i,k)=vgri(i,k,l)*dty
495  alf(i,k)=fm(i,k)/sm(i,k,l)
496  sm(i,k,l)=sm(i,k,l)-fm(i,k)
497  sm0=sm0+fm(i,k)
498  ENDIF
499 C
500  alfq(i,k)=alf(i,k)*alf(i,k)
501  alf1(i,k)=1.-alf(i,k)
502  alf1q(i,k)=alf1(i,k)*alf1(i,k)
503  alf2(i,k)=alf1(i,k)-alf(i,k)
504  alf3(i,k)=alf(i,k)*alfq(i,k)
505  alf4(i,k)=alf1(i,k)*alf1q(i,k)
506 C
507  41 CONTINUE
508 c print*,'ap ADVYP 41'
509 C
510  DO 42 jv=1,ntra
511  DO 420 i=1,lon
512 C
513  IF(vgri(i,k,l).GE.0.) THEN
514  f0(i,k,jv)=alf(i,k)* ( s0(i,k,l,jv)+alf1(i,k)*
515  + ( sy(i,k,l,jv)+alf2(i,k)*syy(i,k,l,jv) ) )
516  s00(jv)=s00(jv)+f0(i,k,jv)
517 C
518  s0(i,k,l,jv)=s0(i,k,l,jv)-f0(i,k,jv)
519  sy(i,k,l,jv)=alf1q(i,k)*
520  + (sy(i,k,l,jv)-3.*alf(i,k)*syy(i,k,l,jv))
521  syy(i,k,l,jv)=alf4(i,k)*syy(i,k,l,jv)
522  ssx(i,k,l,jv)=alf1(i,k)*(ssx(i,k,l,jv)-alf(i,k)*ssxy(i,k,l,jv))
523  sz(i,k,l,jv)=alf1(i,k)*(sz(i,k,l,jv)-alf(i,k)*syz(i,k,l,jv))
524  ssxx(i,k,l,jv)=alf1(i,k)*ssxx(i,k,l,jv)
525  ssxz(i,k,l,jv)=alf1(i,k)*ssxz(i,k,l,jv)
526  szz(i,k,l,jv)=alf1(i,k)*szz(i,k,l,jv)
527  ssxy(i,k,l,jv)=alf1q(i,k)*ssxy(i,k,l,jv)
528  syz(i,k,l,jv)=alf1q(i,k)*syz(i,k,l,jv)
529  ENDIF
530 C
531  420 CONTINUE
532  42 CONTINUE
533 c print*,'ap ADVYP 42'
534 C
535  DO 43 i=1,lon
536  IF(vgri(i,k,l).LT.0.) THEN
537  fm(i,k)=-vgri(i,k,l)*dty
538  alf(i,k)=fm(i,k)/sm0
539  ENDIF
540  43 CONTINUE
541 c print*,'ap ADVYP 43'
542 C
543  DO 44 jv=1,ntra
544  DO 440 i=1,lon
545  IF(vgri(i,k,l).LT.0.) THEN
546  f0(i,k,jv)=alf(i,k)*s00(jv)
547  ENDIF
548  440 CONTINUE
549  44 CONTINUE
550 C
551 C puts the temporary moments Fi into appropriate neighboring boxes
552 C
553  DO 45 i=1,lon
554 C
555  IF(vgri(i,k,l).LT.0.) THEN
556  sm(i,k,l)=sm(i,k,l)+fm(i,k)
557  alf(i,k)=fm(i,k)/sm(i,k,l)
558  ENDIF
559 C
560  alfq(i,k)=alf(i,k)*alf(i,k)
561  alf1(i,k)=1.-alf(i,k)
562  alf1q(i,k)=alf1(i,k)*alf1(i,k)
563  alf2(i,k)=alf1(i,k)-alf(i,k)
564  alf3(i,k)=alf1(i,k)*alf(i,k)
565 C
566  45 CONTINUE
567 c print*,'ap ADVYP 45'
568 C
569  DO 46 jv=1,ntra
570  DO 460 i=1,lon
571 C
572  IF(vgri(i,k,l).LT.0.) THEN
573 C
574  temptm=-alf(i,k)*s0(i,k,l,jv)+alf1(i,k)*f0(i,k,jv)
575  s0(i,k,l,jv)=s0(i,k,l,jv)+f0(i,k,jv)
576  syy(i,k,l,jv)=alf1q(i,k)*syy(i,k,l,jv)
577  + +5.*(-alf3(i,k)*sy(i,k,l,jv)+alf2(i,k)*temptm )
578  sy(i,k,l,jv)=alf1(i,k)*sy(i,k,l,jv)+3.*temptm
579  ssxy(i,k,l,jv)=alf1(i,k)*ssxy(i,k,l,jv)-3.*alf(i,k)*ssx(i,k,l,jv)
580  syz(i,k,l,jv)=alf1(i,k)*syz(i,k,l,jv)-3.*alf(i,k)*sz(i,k,l,jv)
581 C
582  ENDIF
583 C
584  460 CONTINUE
585  46 CONTINUE
586 c print*,'ap ADVYP 46'
587 C
588  1 CONTINUE
589 
590 c--------------------------------------------------
591 C bouclage cyclique horizontal .
592 
593  DO l = 1,llm
594  DO jv = 1,ntra
595  DO j = 1,jjp1
596  sm(iip1,j,l) = sm(1,j,l)
597  s0(iip1,j,l,jv) = s0(1,j,l,jv)
598  ssx(iip1,j,l,jv) = ssx(1,j,l,jv)
599  sy(iip1,j,l,jv) = sy(1,j,l,jv)
600  sz(iip1,j,l,jv) = sz(1,j,l,jv)
601  END DO
602  END DO
603  END DO
604 
605 c -------------------------------------------------------------------
606 C *** Test negativite:
607 
608 c DO jv = 1,ntra
609 c DO l = 1,llm
610 c DO j = 1,jjp1
611 c DO i = 1,iip1
612 c IF (s0( i,j,l,jv ).lt.0.) THEN
613 c PRINT*, '------ S0 < 0 en FIN ADVYP ---'
614 c PRINT*, 'S0(',i,j,l,jv,')=', S0(i,j,l,jv)
615 cc STOP
616 c ENDIF
617 c ENDDO
618 c ENDDO
619 c ENDDO
620 c ENDDO
621 
622 
623 c -------------------------------------------------------------------
624 C *** Test : diag de la qtite totale de traceur dans
625 C l'atmosphere avant l'advection en Y
626 
627  DO l = 1,llm
628  DO j = 1,jjp1
629  DO i = 1,iim
630  sqf = sqf + s0(i,j,l,ntra)
631  END DO
632  END DO
633  END DO
634  print*,'---------- DIAG DANS ADVY - SORTIE --------'
635  print*,'sqf=',sqf
636 c print*,'ap ADVYP fin'
637 
638 c-----------------------------------------------------------------
639 C
640  RETURN
641  END
642 
643 
644 
645 
646 
647 
648 
649 
650 
651 
652 
653