source: LMDZ4/trunk/libf/dyn3dpar/advx.F @ 3551

Last change on this file since 3551 was 764, checked in by Laurent Fairhead, 18 years ago

Merge entre la version V3_conv et le HEAD
YM, JG, LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.0 KB
RevLine 
[630]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#include "comconst.h"
28#include "comvert.h"
29
30C  Arguments :
31C  -----------
32C  dtx : frequence fictive d'appel du transport
33C  pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
34
35       INTEGER ntra
36       PARAMETER (ntra = 1)
37
38C ATTENTION partout ou on trouve ntra, insertion de boucle
39C           possible dans l'avenir.
40
41      REAL dtx
42      REAL pbaru ( iip1,jjp1,llm )
43
44C  moments: SM  total mass in each grid box
45C           S0  mass of tracer in each grid box
46C           Si  1rst order moment in i direction
47C
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
53C  Local :
54C  -------
55
56C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
57C  mass fluxes in kg
58C  declaration :
59
60      REAL UGRI(iip1,jjp1,llm)
61
62C  Rem : VGRI et WGRI ne sont pas utilises dans
63C  cette subroutine ( advection en x uniquement )
64C
65C  Ti are the moments for the current latitude and level
66C
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
71C
72C  the moments F are similarly defined and used as temporary
73C  storage for portions of the grid boxes in transit
74C
75      REAL FM(iim)
76      REAL F0(iim,ntra),FX(iim,ntra)
77      REAL FY(iim,ntra),FZ(iim,ntra)
78C
79C  work arrays
80C
81      REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
82C
83      REAL SMNEW(iim),UEXT(iim)
84C
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
95C *** Test de passage d'arguments ******
96
97
98C  -------------------------------------
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
[764]108               sqi = sqi + S0(i,j,l,ntra)
[630]109            ENDDO
110         ENDDO
111      ENDDO
112      PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
113      PRINT*,'sqi=',sqi
114
115
116C  Interface : adaptation nouveau modele
117C  -------------------------------------
118C
119C  ---------------------------------------------------------
120C  Conversion des flux de masses en kg/s
121C  pbaru est en N/s d'ou :
122C  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 
127C            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
132C  ---------------------------------------------------------
133C  ---------------------------------------------------------
134C  ---------------------------------------------------------
135 
136C  start here         
137C
138C  boucle principale sur les niveaux et les latitudes
139C
140      DO 1 L=1,NIV
141      DO 1 K=lati,latf
142C
143C  initialisation
144C
145C  program assumes periodic boundaries in X
146C
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
151C
152C  modifications for extended polar zones
153C
154      NUMK=NUM(K)
155      LONK=LON/NUMK
156C
157      IF(NUMK.GT.1) THEN
158C
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
170C
171      DO 11 I2=1,NUMK
172C
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
179C
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
192C
193 11   CONTINUE
194C
195      ELSE
196C
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
208C
209      ENDIF
210C
211      DO 117 I=1,LONK
212         UEXT(I)=UGRI(I*NUMK,K,L)
213 117  CONTINUE
214C
215C  place limits on appropriate moments before transport
216C      (if flux-limiting is to be applied)
217C
218      IF(.NOT.LIMIT) GO TO 13
219C
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
225C
226 13   CONTINUE
227C
228C  calculate flux and moments between adjacent boxes
229C  1- create temporary moments/masses for partial boxes in transit
230C  2- reajusts moments remaining in the box
231C
232C  flux from IP to I if U(I).lt.0
233C
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
241C
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
248C
249C  flux from I to IP if U(I).gt.0
250C
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
258C
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
264C
265      DO 150 JV=1,NTRA
266      DO 1500 I=1,LONK-1
267C
268         IF(UEXT(I).LT.0.) THEN
269C
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)
274C
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)
279C
280         ENDIF
281C
282 1500 CONTINUE
283 150  CONTINUE
284C
285      I=LONK
286      IF(UEXT(I).LT.0.) THEN
287C
288        DO 151 JV=1,NTRA
289C
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)
294C
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)
299C
300 151    CONTINUE
301C
302      ENDIF
303C
304      DO 152 JV=1,NTRA
305      DO 1520 I=1,LONK
306C
307         IF(UEXT(I).GE.0.) THEN
308C
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)
313C
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)
318C
319         ENDIF
320C
321 1520 CONTINUE
322 152  CONTINUE
323C
324C  puts the temporary moments Fi into appropriate neighboring boxes
325C
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
332C
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
339C
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
345C
346      DO 162 I=1,LONK
347         ALF1(I)=1.-ALF(I)
348 162  CONTINUE
349C
350      DO 170 JV=1,NTRA
351      DO 1700 I=1,LONK
352C
353         IF(UEXT(I).LT.0.) THEN
354C
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)
360C
361         ENDIF
362C
363 1700 CONTINUE
364 170  CONTINUE
365C
366      DO 171 JV=1,NTRA
367      DO 1710 I=1,LONK-1
368C
369         IF(UEXT(I).GE.0.) THEN
370C
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)
376C
377         ENDIF
378C
379 1710 CONTINUE
380 171  CONTINUE
381C
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
392C
393C  retour aux mailles d'origine (passage des Tij aux Sij)
394C
395      IF(NUMK.GT.1) THEN
396C
397      DO 180 I2=1,NUMK
398C
399         DO 180 I=1,LONK
400C
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)
405C
406            ALFQ(I)=ALF(I)*ALF(I)
407            ALF1(I)=1.-ALF(I)
408            ALF1Q(I)=ALF1(I)*ALF1(I)
409C
410 180     CONTINUE
411C
412         DO  JV=1,NTRA
413         DO  I=1,LONK
414C
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)
421C
422C   reajusts moments remaining in the box
423C
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
430C
431C
432      ELSE
433C
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
445C
446      ENDIF
447C
448 1    CONTINUE
449C
450C ----------- AA Test en fin de ADVX ------ Controle des S*
451c OK
452c      DO 9998 l = 1, llm
453c      DO 9998 j = 1, jjp1
454c      DO 9998 i = 1, iip1
455c         IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
456c            PRINT*, '-------------------'
457c            PRINT*, 'En fin de ADVX'
458c            PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
459c            PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
460c            print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
461c            print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
462c            print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
463c            WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
464cc            STOP
465c         ENDIF
466c 9998 CONTINUE
467c
468C ---------- 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
481c ----------- qqtite totale de traceur dans tte l'atmosphere
482      DO l = 1, llm
483        DO j = 1, jjp1
484          DO i = 1, iim
[764]485             sqf = sqf + S0(i,j,l,ntra)
[630]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.