source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn3d_common/advxp.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 23.0 KB
Line 
1!
2! $Header$
3!
4       SUBROUTINE ADVXP(LIMIT,DTX,PBARU,SM,S0,SSX,SY,SZ
5     .                ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra)
6       IMPLICIT NONE
7CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
8C                                                                 C
9C  second-order moments (SOM) advection of tracer in X direction  C
10C                                                                 C
11CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
12C
13C  parametres principaux du modele
14C
15!-----------------------------------------------------------------------
16!   INCLUDE 'dimensions.h'
17!
18!   dimensions.h contient les dimensions du modele
19!   ndm est tel que iim=2**ndm
20!-----------------------------------------------------------------------
21
22      INTEGER iim,jjm,llm,ndm
23
24      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
25
26!-----------------------------------------------------------------------
27!
28! $Header$
29!
30!
31!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
32!                 veillez  n'utiliser que des ! pour les commentaires
33!                 et  bien positionner les & des lignes de continuation
34!                 (les placer en colonne 6 et en colonne 73)
35!
36!
37!-----------------------------------------------------------------------
38!   INCLUDE 'paramet.h'
39
40      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
41      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
42      INTEGER  ijmllm,mvar
43      INTEGER jcfil,jcfllm
44
45      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
46     &    ,jjp1=jjm+1-1/jjm)
47      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
48      PARAMETER( kftd  = iim/2 -ndm )
49      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
50      PARAMETER( ip1jmi1= ip1jm - iip1 )
51      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
52      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
53      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
54
55!-----------------------------------------------------------------------
56!
57! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
58!
59!-----------------------------------------------------------------------
60! INCLUDE comconst.h
61
62      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
63     &                 iflag_top_bound,mode_top_bound
64      COMMON/comconstr/dtvr,daysec,                                     &
65     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
66     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
67     & ,dissip_pupstart  ,tau_top_bound,                                &
68     & daylen,molmass, ihf
69      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
70
71      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
72      REAL dtvr ! dynamical time step (in s)
73      REAL daysec !length (in s) of a standard day
74      REAL pi    ! something like 3.14159....
75      REAL dtphys ! (s) time step for the physics
76      REAL dtdiss ! (s) time step for the dissipation
77      REAL rad ! (m) radius of the planet
78      REAL r ! Reduced Gas constant r=R/mu
79             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
80      REAL cpp   ! Cp
81      REAL kappa ! kappa=R/Cp
82      REAL cotot
83      REAL unsim ! = 1./iim
84      REAL g ! (m/s2) gravity
85      REAL omeg ! (rad/s) rotation rate of the planet
86! Dissipation factors, for Earth model:
87      REAL dissip_factz,dissip_zref !dissip_deltaz
88! Dissipation factors, for other planets:
89      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
90      REAL dissip_pupstart
91      INTEGER iflag_top_bound,mode_top_bound
92      REAL tau_top_bound
93      REAL daylen ! length of solar day, in 'standard' day length
94      REAL molmass ! (g/mol) molar mass of the atmosphere
95
96      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
97      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
98
99
100!-----------------------------------------------------------------------
101!
102! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
103!
104!-----------------------------------------------------------------------
105!   INCLUDE 'comvert.h'
106
107      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
108     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
109     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
110
111      common/comverti/disvert_type, pressure_exner
112
113      real ap     ! hybrid pressure contribution at interlayers
114      real bp     ! hybrid sigma contribution at interlayer
115      real presnivs ! (reference) pressure at mid-layers
116      real dpres
117      real pa     ! reference pressure (Pa) at which hybrid coordinates
118                  ! become purely pressure
119      real preff  ! reference surface pressure (Pa)
120      real nivsigs
121      real nivsig
122      real aps    ! hybrid pressure contribution at mid-layers
123      real bps    ! hybrid sigma contribution at mid-layers
124      real scaleheight ! atmospheric (reference) scale height (km)
125      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
126                     ! preff and scaleheight
127
128      integer disvert_type ! type of vertical discretization:
129                           ! 1: Earth (default for planet_type==earth),
130                           !     automatic generation
131                           ! 2: Planets (default for planet_type!=earth),
132                           !     using 'z2sig.def' (or 'esasig.def) file
133
134      logical pressure_exner
135!     compute pressure inside layers using Exner function, else use mean
136!     of pressure values at interfaces
137
138 !-----------------------------------------------------------------------
139
140       INTEGER ntra
141c      PARAMETER (ntra = 1)
142C
143C  definition de la grille du modele
144C
145      REAL dtx
146      REAL pbaru ( iip1,jjp1,llm )
147C
148C  moments: SM  total mass in each grid box
149C           S0  mass of tracer in each grid box
150C           Si  1rst order moment in i direction
151C           Sij 2nd  order moment in i and j directions
152C
153      REAL SM(iip1,jjp1,llm)
154     +    ,S0(iip1,jjp1,llm,ntra)
155      REAL SSX(iip1,jjp1,llm,ntra)
156     +    ,SY(iip1,jjp1,llm,ntra)
157     +    ,SZ(iip1,jjp1,llm,ntra)
158      REAL SSXX(iip1,jjp1,llm,ntra)
159     +    ,SSXY(iip1,jjp1,llm,ntra)
160     +    ,SSXZ(iip1,jjp1,llm,ntra)
161     +    ,SYY(iip1,jjp1,llm,ntra)
162     +    ,SYZ(iip1,jjp1,llm,ntra)
163     +    ,SZZ(iip1,jjp1,llm,ntra)
164
165C  Local :
166C  -------
167
168C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
169C  mass fluxes in kg
170C  declaration :
171
172       REAL UGRI(iip1,jjp1,llm)
173
174C  Rem : VGRI et WGRI ne sont pas utilises dans
175C  cette subroutine ( advection en x uniquement )
176C
177C
178C  Tij are the moments for the current latitude and level
179C
180      REAL TM (iim)
181      REAL T0 (iim,NTRA),TX (iim,NTRA)
182      REAL TY (iim,NTRA),TZ (iim,NTRA)
183      REAL TXX(iim,NTRA),TXY(iim,NTRA)
184      REAL TXZ(iim,NTRA),TYY(iim,NTRA)
185      REAL TYZ(iim,NTRA),TZZ(iim,NTRA)
186C
187C  the moments F are similarly defined and used as temporary
188C  storage for portions of the grid boxes in transit
189C
190      REAL FM (iim)
191      REAL F0 (iim,NTRA),FX (iim,NTRA)
192      REAL FY (iim,NTRA),FZ (iim,NTRA)
193      REAL FXX(iim,NTRA),FXY(iim,NTRA)
194      REAL FXZ(iim,NTRA),FYY(iim,NTRA)
195      REAL FYZ(iim,NTRA),FZZ(iim,NTRA)
196C
197C  work arrays
198C
199      REAL ALF (iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
200      REAL ALF2(iim),ALF3(iim),ALF4(iim)
201C
202      REAL SMNEW(iim),UEXT(iim)
203      REAL sqi,sqf
204      REAL TEMPTM
205      REAL SLPMAX
206      REAL S1MAX,S1NEW,S2NEW
207
208      LOGICAL LIMIT
209      INTEGER NUM(jjp1),LONK,NUMK
210      INTEGER lon,lati,latf,niv
211      INTEGER i,i2,i3,j,jv,l,k,iter
212
213      lon = iim
214      lati=2
215      latf = jjm
216      niv = llm
217
218C *** Test de passage d'arguments ******
219
220c      DO 399 l = 1, llm
221c       DO 399 j = 1, jjp1
222c        DO 399 i = 1, iip1
223c         IF (S0(i,j,l,ntra) .lt. 0. ) THEN
224c         PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
225c            print*, 'SSX(',i,j,l,')=',SSX(i,j,l,ntra)
226c         print*, 'SY(',i,j,l,')=',SY(i,j,l,ntra)
227c         print*, 'SZ(',i,j,l,')=',SZ(i,j,l,ntra)
228c         PRINT*, 'AIE !! debut ADVXP - pbl arg. passage dans ADVXP'
229cc            STOP
230c         ENDIF
231c  399 CONTINUE
232
233C *** Test : diagnostique de la qtite totale de traceur
234C            dans l'atmosphere avant l'advection
235c
236      sqi =0.
237      sqf =0.
238c
239      DO l = 1, llm
240      DO j = 1, jjp1
241      DO i = 1, iim
242         sqi = sqi + S0(i,j,l,ntra)
243      END DO
244      END DO
245      END DO
246      PRINT*,'------ DIAG DANS ADVX2 - ENTREE -----'
247      PRINT*,'sqi=',sqi
248c test
249c  -------------------------------------
250        DO 300 j =1,jjp1
251         NUM(j) =1 
252 300  CONTINUE
253c       DO l=1,llm
254c      NUM(2,l)=6
255c      NUM(3,l)=6
256c      NUM(jjm-1,l)=6 
257c      NUM(jjm,l)=6
258c      ENDDO
259c        DO j=2,6
260c       NUM(j)=12
261c       ENDDO
262c       DO j=jjm-5,jjm-1 
263c       NUM(j)=12
264c       ENDDO
265
266C  Interface : adaptation nouveau modele
267C  -------------------------------------
268C
269C  ---------------------------------------------------------
270C  Conversion des flux de masses en kg/s
271C  pbaru est en N/s d'ou :
272C  ugri est en kg/s
273
274       DO 500 l = 1,llm
275       DO 500 j = 1,jjp1
276       DO 500 i = 1,iip1
277       ugri (i,j,llm+1-l) =pbaru (i,j,l)
278 500   CONTINUE
279
280C  ---------------------------------------------------------
281C  start here
282C
283C  boucle principale sur les niveaux et les latitudes
284C     
285      DO 1 L=1,NIV
286      DO 1 K=lati,latf
287
288C
289C  initialisation
290C
291C  program assumes periodic boundaries in X
292C
293      DO 10 I=2,LON
294         SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
295 10   CONTINUE
296      SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
297C
298C  modifications for extended polar zones
299C
300      NUMK=NUM(K)
301      LONK=LON/NUMK
302C
303      IF(NUMK.GT.1) THEN
304C
305      DO 111 I=1,LON
306         TM(I)=0.
307 111  CONTINUE
308      DO 112 JV=1,NTRA
309      DO 1120 I=1,LON
310         T0 (I,JV)=0.
311         TX (I,JV)=0.
312         TY (I,JV)=0.
313         TZ (I,JV)=0.
314         TXX(I,JV)=0.
315         TXY(I,JV)=0.
316         TXZ(I,JV)=0.
317         TYY(I,JV)=0.
318         TYZ(I,JV)=0.
319         TZZ(I,JV)=0.
320 1120 CONTINUE
321 112  CONTINUE
322C
323      DO 11 I2=1,NUMK
324C
325         DO 113 I=1,LONK
326            I3=(I-1)*NUMK+I2
327            TM(I)=TM(I)+SM(I3,K,L)
328            ALF(I)=SM(I3,K,L)/TM(I)
329            ALF1(I)=1.-ALF(I)
330            ALFQ(I)=ALF(I)*ALF(I)
331            ALF1Q(I)=ALF1(I)*ALF1(I)
332            ALF2(I)=ALF1(I)-ALF(I)
333            ALF3(I)=ALF(I)*ALF1(I)
334 113     CONTINUE
335C
336         DO 114 JV=1,NTRA
337         DO 1140 I=1,LONK
338            I3=(I-1)*NUMK+I2
339            TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*S0(I3,K,L,JV)
340            T0 (I,JV)=T0(I,JV)+S0(I3,K,L,JV)
341            TXX(I,JV)=ALFQ(I)*SSXX(I3,K,L,JV)+ALF1Q(I)*TXX(I,JV)
342     +        +5.*( ALF3(I)*(SSX(I3,K,L,JV)-TX(I,JV))+ALF2(I)*TEMPTM )
343            TX (I,JV)=ALF(I)*SSX(I3,K,L,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
344            TXY(I,JV)=ALF (I)*SSXY(I3,K,L,JV)+ALF1(I)*TXY(I,JV)
345     +           +3.*(ALF1(I)*SY (I3,K,L,JV)-ALF (I)*TY (I,JV))
346            TXZ(I,JV)=ALF (I)*SSXZ(I3,K,L,JV)+ALF1(I)*TXZ(I,JV)
347     +           +3.*(ALF1(I)*SZ (I3,K,L,JV)-ALF (I)*TZ (I,JV))
348            TY (I,JV)=TY (I,JV)+SY (I3,K,L,JV)
349            TZ (I,JV)=TZ (I,JV)+SZ (I3,K,L,JV)
350            TYY(I,JV)=TYY(I,JV)+SYY(I3,K,L,JV)
351            TYZ(I,JV)=TYZ(I,JV)+SYZ(I3,K,L,JV)
352            TZZ(I,JV)=TZZ(I,JV)+SZZ(I3,K,L,JV)
353 1140    CONTINUE
354 114     CONTINUE
355C
356 11   CONTINUE
357C
358      ELSE
359C
360      DO 115 I=1,LON
361         TM(I)=SM(I,K,L)
362 115  CONTINUE
363      DO 116 JV=1,NTRA
364      DO 1160 I=1,LON
365         T0 (I,JV)=S0 (I,K,L,JV)
366         TX (I,JV)=SSX (I,K,L,JV)
367         TY (I,JV)=SY (I,K,L,JV)
368         TZ (I,JV)=SZ (I,K,L,JV)
369         TXX(I,JV)=SSXX(I,K,L,JV)
370         TXY(I,JV)=SSXY(I,K,L,JV)
371         TXZ(I,JV)=SSXZ(I,K,L,JV)
372         TYY(I,JV)=SYY(I,K,L,JV)
373         TYZ(I,JV)=SYZ(I,K,L,JV)
374         TZZ(I,JV)=SZZ(I,K,L,JV)
375 1160 CONTINUE
376 116  CONTINUE
377C
378      ENDIF
379C
380      DO 117 I=1,LONK
381         UEXT(I)=UGRI(I*NUMK,K,L)
382 117  CONTINUE
383C
384C  place limits on appropriate moments before transport
385C      (if flux-limiting is to be applied)
386C
387      IF(.NOT.LIMIT) GO TO 13
388C
389      DO 12 JV=1,NTRA
390      DO 120 I=1,LONK
391        IF(T0(I,JV).GT.0.) THEN
392          SLPMAX=T0(I,JV)
393          S1MAX=1.5*SLPMAX
394          S1NEW=AMIN1(S1MAX,AMAX1(-S1MAX,TX(I,JV)))
395          S2NEW=AMIN1( 2.*SLPMAX-ABS(S1NEW)/3. ,
396     +                 AMAX1(ABS(S1NEW)-SLPMAX,TXX(I,JV)) )
397          TX (I,JV)=S1NEW
398          TXX(I,JV)=S2NEW
399          TXY(I,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,TXY(I,JV)))
400          TXZ(I,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,TXZ(I,JV)))
401        ELSE
402          TX (I,JV)=0.
403          TXX(I,JV)=0.
404          TXY(I,JV)=0.
405          TXZ(I,JV)=0.
406        ENDIF
407 120  CONTINUE
408 12   CONTINUE
409C
410 13   CONTINUE
411C
412C  calculate flux and moments between adjacent boxes
413C  1- create temporary moments/masses for partial boxes in transit
414C  2- reajusts moments remaining in the box
415C
416C  flux from IP to I if U(I).lt.0
417C
418      DO 140 I=1,LONK-1
419         IF(UEXT(I).LT.0.) THEN
420           FM(I)=-UEXT(I)*DTX
421           ALF(I)=FM(I)/TM(I+1)
422           TM(I+1)=TM(I+1)-FM(I)
423         ENDIF
424 140  CONTINUE
425C
426      I=LONK
427      IF(UEXT(I).LT.0.) THEN
428        FM(I)=-UEXT(I)*DTX
429        ALF(I)=FM(I)/TM(1)
430        TM(1)=TM(1)-FM(I)
431      ENDIF
432C
433C  flux from I to IP if U(I).gt.0
434C
435      DO 141 I=1,LONK
436         IF(UEXT(I).GE.0.) THEN
437           FM(I)=UEXT(I)*DTX
438           ALF(I)=FM(I)/TM(I)
439           TM(I)=TM(I)-FM(I)
440         ENDIF
441 141  CONTINUE
442C
443      DO 142 I=1,LONK
444         ALFQ(I)=ALF(I)*ALF(I)
445         ALF1(I)=1.-ALF(I)
446         ALF1Q(I)=ALF1(I)*ALF1(I)
447         ALF2(I)=ALF1(I)-ALF(I)
448         ALF3(I)=ALF(I)*ALFQ(I)
449         ALF4(I)=ALF1(I)*ALF1Q(I)
450 142  CONTINUE
451C
452      DO 150 JV=1,NTRA
453      DO 1500 I=1,LONK-1
454C
455         IF(UEXT(I).LT.0.) THEN
456C
457           F0 (I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*
458     +             ( TX(I+1,JV)-ALF2(I)*TXX(I+1,JV) ) )
459           FX (I,JV)=ALFQ(I)*(TX(I+1,JV)-3.*ALF1(I)*TXX(I+1,JV))
460           FXX(I,JV)=ALF3(I)*TXX(I+1,JV)
461           FY (I,JV)=ALF (I)*(TY(I+1,JV)-ALF1(I)*TXY(I+1,JV))
462           FZ (I,JV)=ALF (I)*(TZ(I+1,JV)-ALF1(I)*TXZ(I+1,JV))
463           FXY(I,JV)=ALFQ(I)*TXY(I+1,JV)
464           FXZ(I,JV)=ALFQ(I)*TXZ(I+1,JV)
465           FYY(I,JV)=ALF (I)*TYY(I+1,JV)
466           FYZ(I,JV)=ALF (I)*TYZ(I+1,JV)
467           FZZ(I,JV)=ALF (I)*TZZ(I+1,JV)
468C
469           T0 (I+1,JV)=T0(I+1,JV)-F0(I,JV)
470           TX (I+1,JV)=ALF1Q(I)*(TX(I+1,JV)+3.*ALF(I)*TXX(I+1,JV))
471           TXX(I+1,JV)=ALF4(I)*TXX(I+1,JV)
472           TY (I+1,JV)=TY (I+1,JV)-FY (I,JV)
473           TZ (I+1,JV)=TZ (I+1,JV)-FZ (I,JV)
474           TYY(I+1,JV)=TYY(I+1,JV)-FYY(I,JV)
475           TYZ(I+1,JV)=TYZ(I+1,JV)-FYZ(I,JV)
476           TZZ(I+1,JV)=TZZ(I+1,JV)-FZZ(I,JV)
477           TXY(I+1,JV)=ALF1Q(I)*TXY(I+1,JV)
478           TXZ(I+1,JV)=ALF1Q(I)*TXZ(I+1,JV)
479C
480         ENDIF
481C
482 1500 CONTINUE
483 150  CONTINUE
484C
485      I=LONK
486      IF(UEXT(I).LT.0.) THEN
487C
488        DO 151 JV=1,NTRA
489C
490           F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*
491     +             ( TX(1,JV)-ALF2(I)*TXX(1,JV) ) )
492           FX (I,JV)=ALFQ(I)*(TX(1,JV)-3.*ALF1(I)*TXX(1,JV))
493           FXX(I,JV)=ALF3(I)*TXX(1,JV)
494           FY (I,JV)=ALF (I)*(TY(1,JV)-ALF1(I)*TXY(1,JV))
495           FZ (I,JV)=ALF (I)*(TZ(1,JV)-ALF1(I)*TXZ(1,JV))
496           FXY(I,JV)=ALFQ(I)*TXY(1,JV)
497           FXZ(I,JV)=ALFQ(I)*TXZ(1,JV)
498           FYY(I,JV)=ALF (I)*TYY(1,JV)
499           FYZ(I,JV)=ALF (I)*TYZ(1,JV)
500           FZZ(I,JV)=ALF (I)*TZZ(1,JV)
501C
502           T0 (1,JV)=T0(1,JV)-F0(I,JV)
503           TX (1,JV)=ALF1Q(I)*(TX(1,JV)+3.*ALF(I)*TXX(1,JV))
504           TXX(1,JV)=ALF4(I)*TXX(1,JV)
505           TY (1,JV)=TY (1,JV)-FY (I,JV)
506           TZ (1,JV)=TZ (1,JV)-FZ (I,JV)
507           TYY(1,JV)=TYY(1,JV)-FYY(I,JV)
508           TYZ(1,JV)=TYZ(1,JV)-FYZ(I,JV)
509           TZZ(1,JV)=TZZ(1,JV)-FZZ(I,JV)
510           TXY(1,JV)=ALF1Q(I)*TXY(1,JV)
511           TXZ(1,JV)=ALF1Q(I)*TXZ(1,JV)
512C
513 151    CONTINUE
514C
515      ENDIF
516C
517      DO 152 JV=1,NTRA
518      DO 1520 I=1,LONK
519C
520         IF(UEXT(I).GE.0.) THEN
521C
522           F0 (I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*
523     +             ( TX(I,JV)+ALF2(I)*TXX(I,JV) ) )
524           FX (I,JV)=ALFQ(I)*(TX(I,JV)+3.*ALF1(I)*TXX(I,JV))
525           FXX(I,JV)=ALF3(I)*TXX(I,JV)
526           FY (I,JV)=ALF (I)*(TY(I,JV)+ALF1(I)*TXY(I,JV))
527           FZ (I,JV)=ALF (I)*(TZ(I,JV)+ALF1(I)*TXZ(I,JV))
528           FXY(I,JV)=ALFQ(I)*TXY(I,JV)
529           FXZ(I,JV)=ALFQ(I)*TXZ(I,JV)
530           FYY(I,JV)=ALF (I)*TYY(I,JV)
531           FYZ(I,JV)=ALF (I)*TYZ(I,JV)
532           FZZ(I,JV)=ALF (I)*TZZ(I,JV)
533C
534           T0 (I,JV)=T0(I,JV)-F0(I,JV)
535           TX (I,JV)=ALF1Q(I)*(TX(I,JV)-3.*ALF(I)*TXX(I,JV))
536           TXX(I,JV)=ALF4(I)*TXX(I,JV)
537           TY (I,JV)=TY (I,JV)-FY (I,JV)
538           TZ (I,JV)=TZ (I,JV)-FZ (I,JV)
539           TYY(I,JV)=TYY(I,JV)-FYY(I,JV)
540           TYZ(I,JV)=TYZ(I,JV)-FYZ(I,JV)
541           TZZ(I,JV)=TZZ(I,JV)-FZZ(I,JV)
542           TXY(I,JV)=ALF1Q(I)*TXY(I,JV)
543           TXZ(I,JV)=ALF1Q(I)*TXZ(I,JV)
544C
545         ENDIF
546C
547 1520 CONTINUE
548 152  CONTINUE
549C
550C  puts the temporary moments Fi into appropriate neighboring boxes
551C
552      DO 160 I=1,LONK
553         IF(UEXT(I).LT.0.) THEN
554           TM(I)=TM(I)+FM(I)
555           ALF(I)=FM(I)/TM(I)
556         ENDIF
557 160  CONTINUE
558C
559      DO 161 I=1,LONK-1
560         IF(UEXT(I).GE.0.) THEN
561           TM(I+1)=TM(I+1)+FM(I)
562           ALF(I)=FM(I)/TM(I+1)
563         ENDIF
564 161  CONTINUE
565C
566      I=LONK
567      IF(UEXT(I).GE.0.) THEN
568        TM(1)=TM(1)+FM(I)
569        ALF(I)=FM(I)/TM(1)
570      ENDIF
571C
572      DO 162 I=1,LONK
573         ALF1(I)=1.-ALF(I)
574         ALFQ(I)=ALF(I)*ALF(I)
575         ALF1Q(I)=ALF1(I)*ALF1(I)
576         ALF2(I)=ALF1(I)-ALF(I)
577         ALF3(I)=ALF(I)*ALF1(I)
578 162  CONTINUE
579C
580      DO 170 JV=1,NTRA
581      DO 1700 I=1,LONK
582C
583         IF(UEXT(I).LT.0.) THEN
584C
585           TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
586           T0 (I,JV)=T0(I,JV)+F0(I,JV)
587           TXX(I,JV)=ALFQ(I)*FXX(I,JV)+ALF1Q(I)*TXX(I,JV)
588     +          +5.*( ALF3(I)*(FX(I,JV)-TX(I,JV))+ALF2(I)*TEMPTM )
589           TX (I,JV)=ALF (I)*FX (I,JV)+ALF1(I)*TX (I,JV)+3.*TEMPTM
590           TXY(I,JV)=ALF (I)*FXY(I,JV)+ALF1(I)*TXY(I,JV)
591     +          +3.*(ALF1(I)*FY (I,JV)-ALF (I)*TY (I,JV))
592           TXZ(I,JV)=ALF (I)*FXZ(I,JV)+ALF1(I)*TXZ(I,JV)
593     +          +3.*(ALF1(I)*FZ (I,JV)-ALF (I)*TZ (I,JV))
594           TY (I,JV)=TY (I,JV)+FY (I,JV)
595           TZ (I,JV)=TZ (I,JV)+FZ (I,JV)
596           TYY(I,JV)=TYY(I,JV)+FYY(I,JV)
597           TYZ(I,JV)=TYZ(I,JV)+FYZ(I,JV)
598           TZZ(I,JV)=TZZ(I,JV)+FZZ(I,JV)
599C
600         ENDIF
601C
602 1700 CONTINUE
603 170  CONTINUE
604C
605      DO 171 JV=1,NTRA
606      DO 1710 I=1,LONK-1
607C
608         IF(UEXT(I).GE.0.) THEN
609C
610           TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
611           T0 (I+1,JV)=T0(I+1,JV)+F0(I,JV)
612           TXX(I+1,JV)=ALFQ(I)*FXX(I,JV)+ALF1Q(I)*TXX(I+1,JV)
613     +           +5.*( ALF3(I)*(TX(I+1,JV)-FX(I,JV))-ALF2(I)*TEMPTM )
614           TX (I+1,JV)=ALF(I)*FX (I  ,JV)+ALF1(I)*TX (I+1,JV)+3.*TEMPTM
615           TXY(I+1,JV)=ALF(I)*FXY(I  ,JV)+ALF1(I)*TXY(I+1,JV)
616     +            +3.*(ALF(I)*TY (I+1,JV)-ALF1(I)*FY (I  ,JV))
617           TXZ(I+1,JV)=ALF(I)*FXZ(I  ,JV)+ALF1(I)*TXZ(I+1,JV)
618     +            +3.*(ALF(I)*TZ (I+1,JV)-ALF1(I)*FZ (I  ,JV))
619           TY (I+1,JV)=TY (I+1,JV)+FY (I,JV)
620           TZ (I+1,JV)=TZ (I+1,JV)+FZ (I,JV)
621           TYY(I+1,JV)=TYY(I+1,JV)+FYY(I,JV)
622           TYZ(I+1,JV)=TYZ(I+1,JV)+FYZ(I,JV)
623           TZZ(I+1,JV)=TZZ(I+1,JV)+FZZ(I,JV)
624C
625         ENDIF
626C
627 1710 CONTINUE
628 171  CONTINUE
629C
630      I=LONK
631      IF(UEXT(I).GE.0.) THEN
632        DO 172 JV=1,NTRA
633           TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
634           T0 (1,JV)=T0(1,JV)+F0(I,JV)
635           TXX(1,JV)=ALFQ(I)*FXX(I,JV)+ALF1Q(I)*TXX(1,JV)
636     +         +5.*( ALF3(I)*(TX(1,JV)-FX(I,JV))-ALF2(I)*TEMPTM )
637           TX (1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
638           TXY(1,JV)=ALF(I)*FXY(I,JV)+ALF1(I)*TXY(1,JV)
639     +          +3.*(ALF(I)*TY (1,JV)-ALF1(I)*FY (I,JV))
640           TXZ(1,JV)=ALF(I)*FXZ(I,JV)+ALF1(I)*TXZ(1,JV)
641     +          +3.*(ALF(I)*TZ (1,JV)-ALF1(I)*FZ (I,JV))
642           TY (1,JV)=TY (1,JV)+FY (I,JV)
643           TZ (1,JV)=TZ (1,JV)+FZ (I,JV)
644           TYY(1,JV)=TYY(1,JV)+FYY(I,JV)
645           TYZ(1,JV)=TYZ(1,JV)+FYZ(I,JV)
646           TZZ(1,JV)=TZZ(1,JV)+FZZ(I,JV)
647 172    CONTINUE
648      ENDIF
649C
650C  retour aux mailles d'origine (passage des Tij aux Sij)
651C
652      IF(NUMK.GT.1) THEN
653C
654      DO 18 I2=1,NUMK
655C
656         DO 180 I=1,LONK
657C
658            I3=I2+(I-1)*NUMK
659            SM(I3,K,L)=SMNEW(I3)
660            ALF(I)=SMNEW(I3)/TM(I)
661            TM(I)=TM(I)-SMNEW(I3)
662C
663            ALFQ(I)=ALF(I)*ALF(I)
664            ALF1(I)=1.-ALF(I)
665            ALF1Q(I)=ALF1(I)*ALF1(I)
666            ALF2(I)=ALF1(I)-ALF(I)
667            ALF3(I)=ALF(I)*ALFQ(I)
668            ALF4(I)=ALF1(I)*ALF1Q(I)
669C
670 180     CONTINUE
671C
672         DO 181 JV=1,NTRA
673         DO 181 I=1,LONK
674C
675            I3=I2+(I-1)*NUMK
676            S0 (I3,K,L,JV)=ALF (I)* ( T0(I,JV)-ALF1(I)*
677     +              ( TX(I,JV)-ALF2(I)*TXX(I,JV) ) )
678            SSX (I3,K,L,JV)=ALFQ(I)*(TX(I,JV)-3.*ALF1(I)*TXX(I,JV))
679            SSXX(I3,K,L,JV)=ALF3(I)*TXX(I,JV)
680            SY (I3,K,L,JV)=ALF (I)*(TY(I,JV)-ALF1(I)*TXY(I,JV))
681            SZ (I3,K,L,JV)=ALF (I)*(TZ(I,JV)-ALF1(I)*TXZ(I,JV))
682            SSXY(I3,K,L,JV)=ALFQ(I)*TXY(I,JV)
683            SSXZ(I3,K,L,JV)=ALFQ(I)*TXZ(I,JV)
684            SYY(I3,K,L,JV)=ALF (I)*TYY(I,JV)
685            SYZ(I3,K,L,JV)=ALF (I)*TYZ(I,JV)
686            SZZ(I3,K,L,JV)=ALF (I)*TZZ(I,JV)
687C
688C   reajusts moments remaining in the box
689C
690            T0 (I,JV)=T0(I,JV)-S0(I3,K,L,JV)
691            TX (I,JV)=ALF1Q(I)*(TX(I,JV)+3.*ALF(I)*TXX(I,JV))
692            TXX(I,JV)=ALF4 (I)*TXX(I,JV)
693            TY (I,JV)=TY (I,JV)-SY (I3,K,L,JV)
694            TZ (I,JV)=TZ (I,JV)-SZ (I3,K,L,JV)
695            TYY(I,JV)=TYY(I,JV)-SYY(I3,K,L,JV)
696            TYZ(I,JV)=TYZ(I,JV)-SYZ(I3,K,L,JV)
697            TZZ(I,JV)=TZZ(I,JV)-SZZ(I3,K,L,JV)
698            TXY(I,JV)=ALF1Q(I)*TXY(I,JV)
699            TXZ(I,JV)=ALF1Q(I)*TXZ(I,JV)
700C
701 181     CONTINUE
702C
703 18   CONTINUE
704C
705      ELSE
706C
707      DO 190 I=1,LON
708         SM(I,K,L)=TM(I)
709 190  CONTINUE
710      DO 191 JV=1,NTRA
711      DO 1910 I=1,LON
712         S0 (I,K,L,JV)=T0 (I,JV)
713         SSX (I,K,L,JV)=TX (I,JV)
714         SY (I,K,L,JV)=TY (I,JV)
715         SZ (I,K,L,JV)=TZ (I,JV)
716         SSXX(I,K,L,JV)=TXX(I,JV)
717         SSXY(I,K,L,JV)=TXY(I,JV)
718         SSXZ(I,K,L,JV)=TXZ(I,JV)
719         SYY(I,K,L,JV)=TYY(I,JV)
720         SYZ(I,K,L,JV)=TYZ(I,JV)
721         SZZ(I,K,L,JV)=TZZ(I,JV)
722 1910 CONTINUE
723 191  CONTINUE
724C
725      ENDIF
726C
727 1    CONTINUE
728C
729C ----------- AA Test en fin de ADVX ------ Controle des S*
730
731c      DO 9999 l = 1, llm
732c      DO 9999 j = 1, jjp1
733c      DO 9999 i = 1, iip1
734c          IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
735c           PRINT*, '-------------------'
736c               PRINT*, 'En fin de ADVXP'
737c           PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
738c               print*, 'SSX(',i,j,l,')=',SSX(i,j,l,ntra)
739c           print*, 'SY(',i,j,l,')=',SY(i,j,l,ntra)
740c               print*, 'SZ(',i,j,l,')=',SZ(i,j,l,ntra)
741c            WRITE (*,*) 'On arrete !! - pbl en fin de ADVXP'
742c            STOP
743c           ENDIF
744c 9999 CONTINUE
745c ---------- bouclage cyclique
746
747      DO l = 1,llm
748      DO j = 1,jjp1
749         SM(iip1,j,l) = SM(1,j,l)
750         S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
751         SSX(iip1,j,l,ntra) = SSX(1,j,l,ntra)
752         SY(iip1,j,l,ntra) = SY(1,j,l,ntra)
753         SZ(iip1,j,l,ntra) = SZ(1,j,l,ntra)
754      END DO
755      END DO
756
757C ----------- qqtite totale de traceur dans tte l'atmosphere
758      DO l = 1, llm
759      DO j = 1, jjp1
760      DO i = 1, iim
761        sqf = sqf + S0(i,j,l,ntra)
762      END DO
763      END DO
764      END DO
765
766      PRINT*,'------ DIAG DANS ADVX2 - SORTIE -----'
767      PRINT*,'sqf=',sqf
768c-------------------------------------------------------------
769      RETURN
770      END
Note: See TracBrowser for help on using the repository browser.