source: LMDZ4/trunk/libf/dyn3d/advy.F @ 988

Last change on this file since 988 was 524, checked in by lmdzadmin, 21 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 KB
RevLine 
[524]1!
2! $Header$
3!
4      SUBROUTINE advy(limit,dty,pbarv,sm,s0,sx,sy,sz)
5      IMPLICIT NONE
6
7CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
8C                                                                C
9C  first-order moments (SOM) advection of tracer in Y direction  C
10C                                                                C
11C  Source : Pascal Simon ( Meteo, CNRM )                         C
12C  Adaptation : A.A. (LGGE)                                      C
13C  Derniere Modif : 15/12/94 LAST
14C                                                                C
15C  sont les arguments d'entree pour le s-pg                      C
16C                                                                C
17C  argument de sortie du s-pg                                    C
18C                                                                C
19CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
20CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
21C
22C  Rem : Probleme aux poles il faut reecrire ce cas specifique
23C        Attention au sens de l'indexation
24C
25C  parametres principaux du modele
26C
27C
28#include "dimensions.h"
29#include "paramet.h"
30#include "comconst.h"
31#include "comvert.h"
32#include "comgeom2.h"
33 
34C  Arguments :
35C  ----------
36C  dty : frequence fictive d'appel du transport
37C  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
38
39      INTEGER lon,lat,niv
40      INTEGER i,j,jv,k,kp,l
41      INTEGER ntra
42      PARAMETER (ntra = 1)
43
44      REAL dty
45      REAL pbarv ( iip1,jjm, llm )
46
47C  moments: SM  total mass in each grid box
48C           S0  mass of tracer in each grid box
49C           Si  1rst order moment in i direction
50C
51      REAL SM(iip1,jjp1,llm)
52     +    ,S0(iip1,jjp1,llm,ntra)
53      REAL sx(iip1,jjp1,llm,ntra)
54     +    ,sy(iip1,jjp1,llm,ntra)
55     +    ,sz(iip1,jjp1,llm,ntra)
56
57
58C  Local :
59C  -------
60
61C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
62C  mass fluxes in kg
63C  declaration :
64
65      REAL VGRI(iip1,0:jjp1,llm)
66
67C  Rem : UGRI et WGRI ne sont pas utilises dans
68C  cette subroutine ( advection en y uniquement )
69C  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
70C
71C  the moments F are similarly defined and used as temporary
72C  storage for portions of the grid boxes in transit
73C
74      REAL F0(iim,0:jjp1,ntra),FM(iim,0:jjp1)
75      REAL FX(iim,jjm,ntra),FY(iim,jjm,ntra)
76      REAL FZ(iim,jjm,ntra)
77      REAL S00(ntra)
78      REAL SM0             ! Just temporal variable
79C
80C  work arrays
81C
82      REAL ALF(iim,0:jjp1),ALF1(iim,0:jjp1)
83      REAL ALFQ(iim,0:jjp1),ALF1Q(iim,0:jjp1)
84      REAL TEMPTM          ! Just temporal variable
85c
86C  Special pour poles
87c
88      REAL sbms,sfms,sfzs,sbmn,sfmn,sfzn
89      REAL sns0(ntra),snsz(ntra),snsm
90      REAL s1v(llm),slatv(llm)
91      REAL qy1(iim,llm,ntra),qylat(iim,llm,ntra)
92      REAL cx1(llm,ntra), cxLAT(llm,ntra)
93      REAL cy1(llm,ntra), cyLAT(llm,ntra)
94      REAL z1(iim), zcos(iim), zsin(iim)
95      real smpn,smps,s0pn,s0ps
96      REAL SSUM
97      EXTERNAL SSUM
98C
99      REAL sqi,sqf
100      LOGICAL LIMIT
101
102      lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
103      lat = jjp1        ! a cause des dim. differentes entre les
104      niv=llm
105
106C
107C  the moments Fi are used as temporary storage for
108C  portions of the grid boxes in transit at the current level
109C
110C  work arrays
111C
112
113      DO l = 1,llm
114         DO j = 1,jjm
115            DO i = 1,iip1 
116            vgri (i,j,llm+1-l)=-1.*pbarv(i,j,l) 
117            enddo
118         enddo
119         do i=1,iip1
120             vgri(i,0,l) = 0.
121             vgri(i,jjp1,l) = 0.
122         enddo
123      enddo
124
125      DO 1 L=1,NIV
126C
127C  place limits on appropriate moments before transport
128C      (if flux-limiting is to be applied)
129C
130      IF(.NOT.LIMIT) GO TO 11
131C
132      DO 10 JV=1,NTRA
133      DO 10 K=1,LAT
134      DO 100 I=1,LON
135         sy(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
136     +                           ABS(sy(I,K,L,JV))),sy(I,K,L,JV))
137 100  CONTINUE
138 10   CONTINUE
139C
140 11   CONTINUE
141C
142C  le flux a travers le pole Nord est traite separement
143C
144      SM0=0.
145      DO 20 JV=1,NTRA
146         S00(JV)=0.
147 20   CONTINUE
148C
149      DO 21 I=1,LON
150C
151         IF(VGRI(I,0,L).LE.0.) THEN
152           FM(I,0)=-VGRI(I,0,L)*DTY
153           ALF(I,0)=FM(I,0)/SM(I,1,L)
154           SM(I,1,L)=SM(I,1,L)-FM(I,0)
155           SM0=SM0+FM(I,0)
156         ENDIF
157C
158         ALFQ(I,0)=ALF(I,0)*ALF(I,0)
159         ALF1(I,0)=1.-ALF(I,0)
160         ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
161C
162 21   CONTINUE
163C
164      DO 22 JV=1,NTRA
165      DO 220 I=1,LON
166C
167         IF(VGRI(I,0,L).LE.0.) THEN
168C
169           F0(I,0,JV)=ALF(I,0)*
170     +               ( S0(I,1,L,JV)-ALF1(I,0)*sy(I,1,L,JV) )
171C
172           S00(JV)=S00(JV)+F0(I,0,JV)
173           S0(I,1,L,JV)=S0(I,1,L,JV)-F0(I,0,JV)
174           sy(I,1,L,JV)=ALF1Q(I,0)*sy(I,1,L,JV)
175           sx(I,1,L,JV)=ALF1 (I,0)*sx(I,1,L,JV)
176           sz(I,1,L,JV)=ALF1 (I,0)*sz(I,1,L,JV)
177C
178         ENDIF
179C
180 220  CONTINUE
181 22   CONTINUE
182C
183      DO 23 I=1,LON
184         IF(VGRI(I,0,L).GT.0.) THEN
185           FM(I,0)=VGRI(I,0,L)*DTY
186           ALF(I,0)=FM(I,0)/SM0
187         ENDIF
188 23   CONTINUE
189C
190      DO 24 JV=1,NTRA
191      DO 240 I=1,LON
192         IF(VGRI(I,0,L).GT.0.) THEN
193           F0(I,0,JV)=ALF(I,0)*S00(JV)
194         ENDIF
195 240  CONTINUE
196 24   CONTINUE
197C
198C  puts the temporary moments Fi into appropriate neighboring boxes
199C
200      DO 25 I=1,LON
201C
202         IF(VGRI(I,0,L).GT.0.) THEN
203           SM(I,1,L)=SM(I,1,L)+FM(I,0)
204           ALF(I,0)=FM(I,0)/SM(I,1,L)
205         ENDIF
206C
207         ALF1(I,0)=1.-ALF(I,0)
208C
209 25   CONTINUE
210C
211      DO 26 JV=1,NTRA
212      DO 260 I=1,LON
213C
214         IF(VGRI(I,0,L).GT.0.) THEN
215C
216         TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
217         S0(I,1,L,JV)=S0(I,1,L,JV)+F0(I,0,JV)
218         sy(I,1,L,JV)=ALF1(I,0)*sy(I,1,L,JV)+3.*TEMPTM
219C
220         ENDIF
221C
222 260  CONTINUE
223 26   CONTINUE
224C
225C  calculate flux and moments between adjacent boxes
226C  1- create temporary moments/masses for partial boxes in transit
227C  2- reajusts moments remaining in the box
228C
229C  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
230C
231      DO 30 K=1,LAT-1
232      KP=K+1
233      DO 300 I=1,LON
234C
235         IF(VGRI(I,K,L).LT.0.) THEN
236           FM(I,K)=-VGRI(I,K,L)*DTY
237           ALF(I,K)=FM(I,K)/SM(I,KP,L)
238           SM(I,KP,L)=SM(I,KP,L)-FM(I,K)
239         ELSE
240           FM(I,K)=VGRI(I,K,L)*DTY
241           ALF(I,K)=FM(I,K)/SM(I,K,L)
242           SM(I,K,L)=SM(I,K,L)-FM(I,K)
243         ENDIF
244C
245         ALFQ(I,K)=ALF(I,K)*ALF(I,K)
246         ALF1(I,K)=1.-ALF(I,K)
247         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
248C
249 300  CONTINUE
250 30   CONTINUE
251C
252      DO 31 JV=1,NTRA
253      DO 31 K=1,LAT-1
254      KP=K+1
255      DO 310 I=1,LON
256C
257         IF(VGRI(I,K,L).LT.0.) THEN
258C
259           F0(I,K,JV)=ALF (I,K)*
260     +                ( S0(I,KP,L,JV)-ALF1(I,K)*sy(I,KP,L,JV) )
261           FY(I,K,JV)=ALFQ(I,K)*sy(I,KP,L,JV)
262           FX(I,K,JV)=ALF (I,K)*sx(I,KP,L,JV)
263           FZ(I,K,JV)=ALF (I,K)*sz(I,KP,L,JV)
264C
265           S0(I,KP,L,JV)=S0(I,KP,L,JV)-F0(I,K,JV)
266           sy(I,KP,L,JV)=ALF1Q(I,K)*sy(I,KP,L,JV)
267           sx(I,KP,L,JV)=sx(I,KP,L,JV)-FX(I,K,JV)
268           sz(I,KP,L,JV)=sz(I,KP,L,JV)-FZ(I,K,JV)
269C
270         ELSE
271C
272           F0(I,K,JV)=ALF (I,K)*
273     +               ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
274           FY(I,K,JV)=ALFQ(I,K)*sy(I,K,L,JV)
275           FX(I,K,JV)=ALF(I,K)*sx(I,K,L,JV)
276           FZ(I,K,JV)=ALF(I,K)*sz(I,K,L,JV)
277C
278           S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,K,JV)
279           sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
280           sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,K,JV)
281           sz(I,K,L,JV)=sz(I,K,L,JV)-FZ(I,K,JV)
282C
283         ENDIF
284C
285 310  CONTINUE
286 31   CONTINUE
287C
288C  puts the temporary moments Fi into appropriate neighboring boxes
289C
290      DO 32 K=1,LAT-1
291      KP=K+1
292      DO 320 I=1,LON
293C
294         IF(VGRI(I,K,L).LT.0.) THEN
295           SM(I,K,L)=SM(I,K,L)+FM(I,K)
296           ALF(I,K)=FM(I,K)/SM(I,K,L)
297         ELSE
298           SM(I,KP,L)=SM(I,KP,L)+FM(I,K)
299           ALF(I,K)=FM(I,K)/SM(I,KP,L)
300         ENDIF
301C
302         ALF1(I,K)=1.-ALF(I,K)
303C
304 320  CONTINUE
305 32   CONTINUE
306C
307      DO 33 JV=1,NTRA
308      DO 33 K=1,LAT-1
309      KP=K+1
310      DO 330 I=1,LON
311C
312         IF(VGRI(I,K,L).LT.0.) THEN
313C
314         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
315         S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
316         sy(I,K,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,K,L,JV)
317     +               +3.*TEMPTM
318         sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,K,JV)
319         sz(I,K,L,JV)=sz(I,K,L,JV)+FZ(I,K,JV)
320C
321         ELSE
322C
323         TEMPTM=ALF(I,K)*S0(I,KP,L,JV)-ALF1(I,K)*F0(I,K,JV)
324         S0(I,KP,L,JV)=S0(I,KP,L,JV)+F0(I,K,JV)
325         sy(I,KP,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,KP,L,JV)
326     +                +3.*TEMPTM
327         sx(I,KP,L,JV)=sx(I,KP,L,JV)+FX(I,K,JV)
328         sz(I,KP,L,JV)=sz(I,KP,L,JV)+FZ(I,K,JV)
329C
330         ENDIF
331C
332 330  CONTINUE
333 33   CONTINUE
334C
335C  traitement special pour le pole Sud (idem pole Nord)
336C
337      K=LAT
338C
339      SM0=0.
340      DO 40 JV=1,NTRA
341         S00(JV)=0.
342 40   CONTINUE
343C
344      DO 41 I=1,LON
345C
346         IF(VGRI(I,K,L).GE.0.) THEN
347           FM(I,K)=VGRI(I,K,L)*DTY
348           ALF(I,K)=FM(I,K)/SM(I,K,L)
349           SM(I,K,L)=SM(I,K,L)-FM(I,K)
350           SM0=SM0+FM(I,K)
351         ENDIF
352C
353         ALFQ(I,K)=ALF(I,K)*ALF(I,K)
354         ALF1(I,K)=1.-ALF(I,K)
355         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
356C
357 41   CONTINUE
358C
359      DO 42 JV=1,NTRA
360      DO 420 I=1,LON
361C
362         IF(VGRI(I,K,L).GE.0.) THEN
363           F0 (I,K,JV)=ALF(I,K)*
364     +                ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
365           S00(JV)=S00(JV)+F0(I,K,JV)
366C
367           S0(I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
368           sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
369           sx(I,K,L,JV)=ALF1(I,K)*sx(I,K,L,JV)
370           sz(I,K,L,JV)=ALF1(I,K)*sz(I,K,L,JV)
371         ENDIF
372C
373 420  CONTINUE
374 42   CONTINUE
375C
376      DO 43 I=1,LON
377         IF(VGRI(I,K,L).LT.0.) THEN
378           FM(I,K)=-VGRI(I,K,L)*DTY
379           ALF(I,K)=FM(I,K)/SM0
380         ENDIF
381 43   CONTINUE
382C
383      DO 44 JV=1,NTRA
384      DO 440 I=1,LON
385         IF(VGRI(I,K,L).LT.0.) THEN
386           F0(I,K,JV)=ALF(I,K)*S00(JV)
387         ENDIF
388 440  CONTINUE
389 44   CONTINUE
390C
391C  puts the temporary moments Fi into appropriate neighboring boxes
392C
393      DO 45 I=1,LON
394C
395         IF(VGRI(I,K,L).LT.0.) THEN
396           SM(I,K,L)=SM(I,K,L)+FM(I,K)
397           ALF(I,K)=FM(I,K)/SM(I,K,L)
398         ENDIF
399C
400         ALF1(I,K)=1.-ALF(I,K)
401C
402 45   CONTINUE
403C
404      DO 46 JV=1,NTRA
405      DO 460 I=1,LON
406C
407         IF(VGRI(I,K,L).LT.0.) THEN
408C
409         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
410         S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
411         sy(I,K,L,JV)=ALF1(I,K)*sy(I,K,L,JV)+3.*TEMPTM
412C
413         ENDIF
414C
415 460  CONTINUE
416 46   CONTINUE
417C
418 1    CONTINUE
419C
420      RETURN
421      END
422
Note: See TracBrowser for help on using the repository browser.