source: LMDZ5/branches/testing/libf/dyn3d_common/advx.F @ 2791

Last change on this file since 2791 was 2641, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2593:2640 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.1 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE  advx(limit,dtx,pbaru,sm,s0,
5     $     sx,sy,sz,lati,latf)
6      IMPLICIT NONE
7
8CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9C                                                                C
10C  first-order moments (FOM) advection of tracer in X direction  C
11C                                                                C
12C  Source : Pascal Simon (Meteo,CNRM)                            C
13C  Adaptation : A.Armengaud (LGGE) juin 94                       C
14C                                                                C
15C  limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz                       C
16C  sont des arguments d'entree pour le s-pg...                   C
17C                                                                C
18C  sm,s0,sx,sy,sz                                                C
19C  sont les arguments de sortie pour le s-pg                     C
20C                                                                C
21CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22C
23C  parametres principaux du modele
24C
25      include "dimensions.h"
26      include "paramet.h"
27
28C  Arguments :
29C  -----------
30C  dtx : frequence fictive d'appel du transport
31C  pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
32
33       INTEGER ntra
34       PARAMETER (ntra = 1)
35
36C ATTENTION partout ou on trouve ntra, insertion de boucle
37C           possible dans l'avenir.
38
39      REAL dtx
40      REAL pbaru ( iip1,jjp1,llm )
41
42C  moments: SM  total mass in each grid box
43C           S0  mass of tracer in each grid box
44C           Si  1rst order moment in i direction
45C
46      REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
47      REAL sx(iip1,jjp1,llm,ntra)
48     $    ,sy(iip1,jjp1,llm,ntra)
49      REAL sz(iip1,jjp1,llm,ntra)
50
51C  Local :
52C  -------
53
54C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
55C  mass fluxes in kg
56C  declaration :
57
58      REAL UGRI(iip1,jjp1,llm)
59
60C  Rem : VGRI et WGRI ne sont pas utilises dans
61C  cette subroutine ( advection en x uniquement )
62C
63C  Ti are the moments for the current latitude and level
64C
65      REAL TM(iim)
66      REAL T0(iim,ntra),TX(iim,ntra)
67      REAL TY(iim,ntra),TZ(iim,ntra)
68      REAL TEMPTM                ! just a temporary variable
69C
70C  the moments F are similarly defined and used as temporary
71C  storage for portions of the grid boxes in transit
72C
73      REAL FM(iim)
74      REAL F0(iim,ntra),FX(iim,ntra)
75      REAL FY(iim,ntra),FZ(iim,ntra)
76C
77C  work arrays
78C
79      REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
80C
81      REAL SMNEW(iim),UEXT(iim)
82C
83      REAL sqi,sqf
84
85      LOGICAL LIMIT
86      INTEGER NUM(jjp1),LONK,NUMK
87      INTEGER lon,lati,latf,niv
88      INTEGER i,i2,i3,j,jv,l,k,itrac
89
90      lon = iim
91      niv = llm
92
93C *** Test de passage d'arguments ******
94
95
96C  -------------------------------------
97      DO 300 j = 1,jjp1
98         NUM(j) = 1
99  300 CONTINUE
100      sqi = 0.
101      sqf = 0.
102
103      DO l = 1,llm
104         DO j = 1,jjp1
105            DO i = 1,iim
106cIM 240305            sqi = sqi + S0(i,j,l,9)
107               sqi = sqi + S0(i,j,l,ntra)
108            ENDDO
109         ENDDO
110      ENDDO
111      PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
112      PRINT*,'sqi=',sqi
113
114
115C  Interface : adaptation nouveau modele
116C  -------------------------------------
117C
118C  ---------------------------------------------------------
119C  Conversion des flux de masses en kg/s
120C  pbaru est en N/s d'ou :
121C  ugri est en kg/s
122
123      DO 500 l = 1,llm
124         DO 500 j = 1,jjm+1
125            DO 500 i = 1,iip1 
126C            ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
127             ugri (i,j,llm+1-l) = pbaru (i,j,l)
128  500 CONTINUE
129
130
131C  ---------------------------------------------------------
132C  ---------------------------------------------------------
133C  ---------------------------------------------------------
134 
135C  start here         
136C
137C  boucle principale sur les niveaux et les latitudes
138C
139      DO 1 L=1,NIV
140      DO 1 K=lati,latf
141C
142C  initialisation
143C
144C  program assumes periodic boundaries in X
145C
146      DO 10 I=2,LON
147         SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
148 10   CONTINUE
149      SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
150C
151C  modifications for extended polar zones
152C
153      NUMK=NUM(K)
154      LONK=LON/NUMK
155C
156      IF(NUMK.GT.1) THEN
157C
158      DO 111 I=1,LON
159         TM(I)=0.
160 111  CONTINUE
161      DO 112 JV=1,NTRA
162      DO 1120 I=1,LON
163         T0(I,JV)=0.
164         TX(I,JV)=0.
165         TY(I,JV)=0.
166         TZ(I,JV)=0.
167 1120 CONTINUE
168 112  CONTINUE
169C
170      DO 11 I2=1,NUMK
171C
172         DO 113 I=1,LONK
173            I3=(I-1)*NUMK+I2
174            TM(I)=TM(I)+SM(I3,K,L)
175            ALF(I)=SM(I3,K,L)/TM(I)
176            ALF1(I)=1.-ALF(I)
177 113     CONTINUE
178C
179         DO  JV=1,NTRA
180         DO  I=1,LONK
181            I3=(I-1)*NUMK+I2
182            TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)
183     $          *S0(I3,K,L,JV)
184            T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
185            TX(I,JV)=ALF(I)  *sx(I3,K,L,JV)+
186     $       ALF1(I)*TX(I,JV) +3.*TEMPTM
187            TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
188            TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
189         ENDDO
190         ENDDO
191C
192 11   CONTINUE
193C
194      ELSE
195C
196      DO 115 I=1,LON
197         TM(I)=SM(I,K,L)
198 115  CONTINUE
199      DO 116 JV=1,NTRA
200      DO 1160 I=1,LON
201         T0(I,JV)=S0(I,K,L,JV)
202         TX(I,JV)=sx(I,K,L,JV)
203         TY(I,JV)=sy(I,K,L,JV)
204         TZ(I,JV)=sz(I,K,L,JV)
205 1160 CONTINUE
206 116  CONTINUE
207C
208      ENDIF
209C
210      DO 117 I=1,LONK
211         UEXT(I)=UGRI(I*NUMK,K,L)
212 117  CONTINUE
213C
214C  place limits on appropriate moments before transport
215C      (if flux-limiting is to be applied)
216C
217      IF(.NOT.LIMIT) GO TO 13
218C
219      DO 12 JV=1,NTRA
220      DO 120 I=1,LONK
221        TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
222 120  CONTINUE
223 12   CONTINUE
224C
225 13   CONTINUE
226C
227C  calculate flux and moments between adjacent boxes
228C  1- create temporary moments/masses for partial boxes in transit
229C  2- reajusts moments remaining in the box
230C
231C  flux from IP to I if U(I).lt.0
232C
233      DO 140 I=1,LONK-1
234         IF(UEXT(I).LT.0.) THEN
235           FM(I)=-UEXT(I)*DTX
236           ALF(I)=FM(I)/TM(I+1)
237           TM(I+1)=TM(I+1)-FM(I)
238         ENDIF
239 140  CONTINUE
240C
241      I=LONK
242      IF(UEXT(I).LT.0.) THEN
243        FM(I)=-UEXT(I)*DTX
244        ALF(I)=FM(I)/TM(1)
245        TM(1)=TM(1)-FM(I)
246      ENDIF
247C
248C  flux from I to IP if U(I).gt.0
249C
250      DO 141 I=1,LONK
251         IF(UEXT(I).GE.0.) THEN
252           FM(I)=UEXT(I)*DTX
253           ALF(I)=FM(I)/TM(I)
254           TM(I)=TM(I)-FM(I)
255         ENDIF
256 141  CONTINUE
257C
258      DO 142 I=1,LONK
259         ALFQ(I)=ALF(I)*ALF(I)
260         ALF1(I)=1.-ALF(I)
261         ALF1Q(I)=ALF1(I)*ALF1(I)
262 142  CONTINUE
263C
264      DO 150 JV=1,NTRA
265      DO 1500 I=1,LONK-1
266C
267         IF(UEXT(I).LT.0.) THEN
268C
269           F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
270           FX(I,JV)=ALFQ(I)*TX(I+1,JV)
271           FY(I,JV)=ALF (I)*TY(I+1,JV)
272           FZ(I,JV)=ALF (I)*TZ(I+1,JV)
273C
274           T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
275           TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
276           TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
277           TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
278C
279         ENDIF
280C
281 1500 CONTINUE
282 150  CONTINUE
283C
284      I=LONK
285      IF(UEXT(I).LT.0.) THEN
286C
287        DO 151 JV=1,NTRA
288C
289           F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
290           FX (I,JV)=ALFQ(I)*TX(1,JV)
291           FY (I,JV)=ALF (I)*TY(1,JV)
292           FZ (I,JV)=ALF (I)*TZ(1,JV)
293C
294           T0(1,JV)=T0(1,JV)-F0(I,JV)
295           TX(1,JV)=ALF1Q(I)*TX(1,JV)
296           TY(1,JV)=TY(1,JV)-FY(I,JV)
297           TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
298C
299 151    CONTINUE
300C
301      ENDIF
302C
303      DO 152 JV=1,NTRA
304      DO 1520 I=1,LONK
305C
306         IF(UEXT(I).GE.0.) THEN
307C
308           F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
309           FX(I,JV)=ALFQ(I)*TX(I,JV)
310           FY(I,JV)=ALF (I)*TY(I,JV)
311           FZ(I,JV)=ALF (I)*TZ(I,JV)
312C
313           T0(I,JV)=T0(I,JV)-F0(I,JV)
314           TX(I,JV)=ALF1Q(I)*TX(I,JV)
315           TY(I,JV)=TY(I,JV)-FY(I,JV)
316           TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
317C
318         ENDIF
319C
320 1520 CONTINUE
321 152  CONTINUE
322C
323C  puts the temporary moments Fi into appropriate neighboring boxes
324C
325      DO 160 I=1,LONK
326         IF(UEXT(I).LT.0.) THEN
327           TM(I)=TM(I)+FM(I)
328           ALF(I)=FM(I)/TM(I)
329         ENDIF
330 160  CONTINUE
331C
332      DO 161 I=1,LONK-1
333         IF(UEXT(I).GE.0.) THEN
334           TM(I+1)=TM(I+1)+FM(I)
335           ALF(I)=FM(I)/TM(I+1)
336         ENDIF
337 161  CONTINUE
338C
339      I=LONK
340      IF(UEXT(I).GE.0.) THEN
341        TM(1)=TM(1)+FM(I)
342        ALF(I)=FM(I)/TM(1)
343      ENDIF
344C
345      DO 162 I=1,LONK
346         ALF1(I)=1.-ALF(I)
347 162  CONTINUE
348C
349      DO 170 JV=1,NTRA
350      DO 1700 I=1,LONK
351C
352         IF(UEXT(I).LT.0.) THEN
353C
354           TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
355           T0(I,JV)=T0(I,JV)+F0(I,JV)
356           TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
357           TY(I,JV)=TY(I,JV)+FY(I,JV)
358           TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
359C
360         ENDIF
361C
362 1700 CONTINUE
363 170  CONTINUE
364C
365      DO 171 JV=1,NTRA
366      DO 1710 I=1,LONK-1
367C
368         IF(UEXT(I).GE.0.) THEN
369C
370           TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
371           T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
372           TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
373           TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
374           TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
375C
376         ENDIF
377C
378 1710 CONTINUE
379 171  CONTINUE
380C
381      I=LONK
382      IF(UEXT(I).GE.0.) THEN
383        DO 172 JV=1,NTRA
384           TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
385           T0(1,JV)=T0(1,JV)+F0(I,JV)
386           TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
387           TY(1,JV)=TY(1,JV)+FY(I,JV)
388           TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
389 172    CONTINUE
390      ENDIF
391C
392C  retour aux mailles d'origine (passage des Tij aux Sij)
393C
394      IF(NUMK.GT.1) THEN
395C
396      DO 180 I2=1,NUMK
397C
398         DO 180 I=1,LONK
399C
400            I3=I2+(I-1)*NUMK
401            SM(I3,K,L)=SMNEW(I3)
402            ALF(I)=SMNEW(I3)/TM(I)
403            TM(I)=TM(I)-SMNEW(I3)
404C
405            ALFQ(I)=ALF(I)*ALF(I)
406            ALF1(I)=1.-ALF(I)
407            ALF1Q(I)=ALF1(I)*ALF1(I)
408C
409 180     CONTINUE
410C
411         DO  JV=1,NTRA
412         DO  I=1,LONK
413C
414            I3=I2+(I-1)*NUMK
415            S0(I3,K,L,JV)=ALF (I)
416     $       * (T0(I,JV)-ALF1(I)*TX(I,JV))
417            sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
418            sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
419            sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
420C
421C   reajusts moments remaining in the box
422C
423            T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
424            TX(I,JV)=ALF1Q(I)*TX(I,JV)
425            TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
426            TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
427          ENDDO
428          ENDDO
429C
430C
431      ELSE
432C
433      DO 190 I=1,LON
434         SM(I,K,L)=TM(I)
435 190  CONTINUE
436      DO 191 JV=1,NTRA
437      DO 1910 I=1,LON
438         S0(I,K,L,JV)=T0(I,JV)
439         sx(I,K,L,JV)=TX(I,JV)
440         sy(I,K,L,JV)=TY(I,JV)
441         sz(I,K,L,JV)=TZ(I,JV)
442 1910 CONTINUE
443 191  CONTINUE
444C
445      ENDIF
446C
447 1    CONTINUE
448C
449C ----------- AA Test en fin de ADVX ------ Controle des S*
450c OK
451c      DO 9998 l = 1, llm
452c      DO 9998 j = 1, jjp1
453c      DO 9998 i = 1, iip1
454c         IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
455c            PRINT*, '-------------------'
456c            PRINT*, 'En fin de ADVX'
457c            PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
458c            PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
459c            print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
460c            print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
461c            print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
462c            WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
463cc            STOP
464c         ENDIF
465c 9998 CONTINUE
466c
467C ---------- bouclage cyclique
468      DO itrac=1,ntra
469      DO l = 1,llm
470        DO j = lati,latf
471           SM(iip1,j,l) = SM(1,j,l)
472           S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
473           sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
474           sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
475           sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
476        END DO
477      END DO
478      ENDDO
479
480c ----------- qqtite totale de traceur dans tte l'atmosphere
481      DO l = 1, llm
482        DO j = 1, jjp1
483          DO i = 1, iim
484cIM 240405          sqf = sqf + S0(i,j,l,9)
485             sqf = sqf + S0(i,j,l,ntra)
486          END DO 
487        END DO
488      END DO
489c
490      PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
491      PRINT*,'sqf=',sqf
492c-------------
493
494      RETURN
495      END
496C_________________________________________________________________
497C_________________________________________________________________
Note: See TracBrowser for help on using the repository browser.