source: LMDZ4/trunk/libf/dyn3d/advx.F @ 5416

Last change on this file since 5416 was 644, checked in by Laurent Fairhead, 20 years ago

Synchronisation avec tous les diagnostiques de Ionela IM
Inclusion du slab ocean IM
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.1 KB
RevLine 
[524]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
[644]108cIM 240305            sqi = sqi + S0(i,j,l,9)
109               sqi = sqi + S0(i,j,l,ntra)
[524]110            ENDDO
111         ENDDO
112      ENDDO
113      PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
114      PRINT*,'sqi=',sqi
115
116
117C  Interface : adaptation nouveau modele
118C  -------------------------------------
119C
120C  ---------------------------------------------------------
121C  Conversion des flux de masses en kg/s
122C  pbaru est en N/s d'ou :
123C  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 
128C            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
133C  ---------------------------------------------------------
134C  ---------------------------------------------------------
135C  ---------------------------------------------------------
136 
137C  start here         
138C
139C  boucle principale sur les niveaux et les latitudes
140C
141      DO 1 L=1,NIV
142      DO 1 K=lati,latf
143C
144C  initialisation
145C
146C  program assumes periodic boundaries in X
147C
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
152C
153C  modifications for extended polar zones
154C
155      NUMK=NUM(K)
156      LONK=LON/NUMK
157C
158      IF(NUMK.GT.1) THEN
159C
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
171C
172      DO 11 I2=1,NUMK
173C
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
180C
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
193C
194 11   CONTINUE
195C
196      ELSE
197C
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
209C
210      ENDIF
211C
212      DO 117 I=1,LONK
213         UEXT(I)=UGRI(I*NUMK,K,L)
214 117  CONTINUE
215C
216C  place limits on appropriate moments before transport
217C      (if flux-limiting is to be applied)
218C
219      IF(.NOT.LIMIT) GO TO 13
220C
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
226C
227 13   CONTINUE
228C
229C  calculate flux and moments between adjacent boxes
230C  1- create temporary moments/masses for partial boxes in transit
231C  2- reajusts moments remaining in the box
232C
233C  flux from IP to I if U(I).lt.0
234C
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
242C
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
249C
250C  flux from I to IP if U(I).gt.0
251C
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
259C
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
265C
266      DO 150 JV=1,NTRA
267      DO 1500 I=1,LONK-1
268C
269         IF(UEXT(I).LT.0.) THEN
270C
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)
275C
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)
280C
281         ENDIF
282C
283 1500 CONTINUE
284 150  CONTINUE
285C
286      I=LONK
287      IF(UEXT(I).LT.0.) THEN
288C
289        DO 151 JV=1,NTRA
290C
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)
295C
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)
300C
301 151    CONTINUE
302C
303      ENDIF
304C
305      DO 152 JV=1,NTRA
306      DO 1520 I=1,LONK
307C
308         IF(UEXT(I).GE.0.) THEN
309C
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)
314C
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)
319C
320         ENDIF
321C
322 1520 CONTINUE
323 152  CONTINUE
324C
325C  puts the temporary moments Fi into appropriate neighboring boxes
326C
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
333C
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
340C
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
346C
347      DO 162 I=1,LONK
348         ALF1(I)=1.-ALF(I)
349 162  CONTINUE
350C
351      DO 170 JV=1,NTRA
352      DO 1700 I=1,LONK
353C
354         IF(UEXT(I).LT.0.) THEN
355C
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)
361C
362         ENDIF
363C
364 1700 CONTINUE
365 170  CONTINUE
366C
367      DO 171 JV=1,NTRA
368      DO 1710 I=1,LONK-1
369C
370         IF(UEXT(I).GE.0.) THEN
371C
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)
377C
378         ENDIF
379C
380 1710 CONTINUE
381 171  CONTINUE
382C
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
393C
394C  retour aux mailles d'origine (passage des Tij aux Sij)
395C
396      IF(NUMK.GT.1) THEN
397C
398      DO 180 I2=1,NUMK
399C
400         DO 180 I=1,LONK
401C
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)
406C
407            ALFQ(I)=ALF(I)*ALF(I)
408            ALF1(I)=1.-ALF(I)
409            ALF1Q(I)=ALF1(I)*ALF1(I)
410C
411 180     CONTINUE
412C
413         DO  JV=1,NTRA
414         DO  I=1,LONK
415C
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)
422C
423C   reajusts moments remaining in the box
424C
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
431C
432C
433      ELSE
434C
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
446C
447      ENDIF
448C
449 1    CONTINUE
450C
451C ----------- AA Test en fin de ADVX ------ Controle des S*
452c OK
453c      DO 9998 l = 1, llm
454c      DO 9998 j = 1, jjp1
455c      DO 9998 i = 1, iip1
456c         IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
457c            PRINT*, '-------------------'
458c            PRINT*, 'En fin de ADVX'
459c            PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
460c            PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
461c            print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
462c            print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
463c            print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
464c            WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
465cc            STOP
466c         ENDIF
467c 9998 CONTINUE
468c
469C ---------- 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
482c ----------- qqtite totale de traceur dans tte l'atmosphere
483      DO l = 1, llm
484        DO j = 1, jjp1
485          DO i = 1, iim
[644]486cIM 240405          sqf = sqf + S0(i,j,l,9)
487             sqf = sqf + S0(i,j,l,ntra)
[524]488          END DO 
489        END DO
490      END DO
491c
492      PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
493      PRINT*,'sqf=',sqf
494c-------------
495
496      RETURN
497      END
498C_________________________________________________________________
499C_________________________________________________________________
Note: See TracBrowser for help on using the repository browser.