/[lmdze]/trunk/dyn3d/advz.f
ViewVC logotype

Diff of /trunk/dyn3d/advz.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/advz.f revision 28 by guez, Fri Mar 26 18:33:04 2010 UTC trunk/dyn3d/advz.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advz.F,v 1.2 2005/05/25 13:10:09 fairhead Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advz.F,v 1.2 2005/05/25 13:10:09
3  !  ! fairhead Exp $
4        SUBROUTINE advz(limit,dtz,w,sm,s0,sx,sy,sz)  
5        use dimens_m  SUBROUTINE advz(limit, dtz, w, sm, s0, sx, sy, sz)
6        use paramet_m    USE dimens_m
7        use comconst    USE paramet_m
8        use comvert    USE comconst
9        IMPLICIT NONE    USE disvert_m
10      IMPLICIT NONE
11  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  
12  C                                                                C    ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
13  C  first-order moments (FOM) advection of tracer in Z direction  C    ! C
14  C                                                                C    ! first-order moments (FOM) advection of tracer in Z direction  C
15  C  Source : Pascal Simon (Meteo,CNRM)                            C    ! C
16  C  Adaptation : A.Armengaud (LGGE) juin 94                       C    ! Source : Pascal Simon (Meteo,CNRM)                            C
17  C                                                                C    ! Adaptation : A.Armengaud (LGGE) juin 94                       C
18  C                                                                C    ! C
19  C  sont des arguments d'entree pour le s-pg...                   C    ! C
20  C                                                                C    ! sont des arguments d'entree pour le s-pg...                   C
21  C  dq est l'argument de sortie pour le s-pg                      C    ! C
22  C                                                                C    ! dq est l'argument de sortie pour le s-pg                      C
23  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC    ! C
24  C    ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
25  C  parametres principaux du modele  
26  C    ! parametres principaux du modele
27    
28  C  Arguments :  
29  C  -----------    ! Arguments :
30  C  dtz : frequence fictive d'appel du transport    ! -----------
31  C  w : flux de masse en z en Pa.m2.s-1    ! dtz : frequence fictive d'appel du transport
32      ! w : flux de masse en z en Pa.m2.s-1
33        INTEGER ntra  
34        PARAMETER (ntra = 1)    INTEGER ntra
35      PARAMETER (ntra=1)
36        REAL, intent(in):: dtz  
37        REAL w ( iip1,jjp1,llm )    REAL, INTENT (IN) :: dtz
38          REAL w(iip1, jjp1, llm)
39  C  moments: SM  total mass in each grid box  
40  C           S0  mass of tracer in each grid box    ! moments: SM  total mass in each grid box
41  C           Si  1rst order moment in i direction    ! S0  mass of tracer in each grid box
42  C    ! Si  1rst order moment in i direction
43        REAL SM(iip1,jjp1,llm)  
44       +    ,S0(iip1,jjp1,llm,ntra)    REAL sm(iip1, jjp1, llm), s0(iip1, jjp1, llm, ntra)
45        REAL sx(iip1,jjp1,llm,ntra)    REAL sx(iip1, jjp1, llm, ntra), sy(iip1, jjp1, llm, ntra), &
46       +    ,sy(iip1,jjp1,llm,ntra)      sz(iip1, jjp1, llm, ntra)
47       +    ,sz(iip1,jjp1,llm,ntra)  
48    
49      ! Local :
50  C  Local :    ! -------
51  C  -------  
52      ! mass fluxes across the boundaries (UGRI,VGRI,WGRI)
53  C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)    ! mass fluxes in kg
54  C  mass fluxes in kg    ! declaration :
55  C  declaration :  
56      REAL wgri(iip1, jjp1, 0:llm)
57        REAL WGRI(iip1,jjp1,0:llm)  
58    
59  C    ! the moments F are used as temporary  storage for
60  C  the moments F are used as temporary  storage for    ! portions of grid boxes in transit at the current latitude
61  C  portions of grid boxes in transit at the current latitude  
62  C    REAL fm(iim, llm)
63        REAL FM(iim,llm)    REAL f0(iim, llm, ntra), fx(iim, llm, ntra)
64        REAL F0(iim,llm,ntra),FX(iim,llm,ntra)    REAL fy(iim, llm, ntra), fz(iim, llm, ntra)
65        REAL FY(iim,llm,ntra),FZ(iim,llm,ntra)  
66  C    ! work arrays
67  C  work arrays  
68  C    REAL alf(iim), alf1(iim), alfq(iim), alf1q(iim)
69        REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)    REAL temptm ! Just temporal variable
70        REAL TEMPTM            ! Just temporal variable    REAL sqi, sqf
71        REAL sqi,sqf  
72  C    LOGICAL limit
73        LOGICAL LIMIT    INTEGER lon, lat, niv
74        INTEGER lon,lat,niv    INTEGER i, j, jv, k, l, lp
75        INTEGER i,j,jv,k,l,lp  
76      lon = iim
77        lon = iim    lat = jjp1
78        lat = jjp1    niv = llm
79        niv = llm  
80      ! *** Test : diag de la qqtite totale de traceur
81  C *** Test de passage d'arguments ******    ! dans l'atmosphere avant l'advection en z
82      sqi = 0.
83  c     DO 399 l = 1, llm    sqf = 0.
84  c     DO 399 j = 1, jjp1  
85  c     DO 399 i = 1, iip1    DO l = 1, llm
86  c        IF (S0(i,j,l,ntra) .lt. 0. ) THEN      DO j = 1, jjp1
87  c           PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)        DO i = 1, iim
88  c           print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)          ! IM 240305            sqi = sqi + S0(i,j,l,9)
89  c           print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)          sqi = sqi + s0(i, j, l, ntra)
90  c           print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)        END DO
91  c           PRINT*, 'AIE !! debut ADVZ - pbl arg. passage dans ADVZ'      END DO
92  c            STOP    END DO
93  c        ENDIF    PRINT *, '-------- DIAG DANS ADVZ - ENTREE ---------'
94    399 CONTINUE    PRINT *, 'sqi=', sqi
95    
96  C-----------------------------------------------------------------    ! -----------------------------------------------------------------
97  C *** Test : diag de la qqtite totale de traceur    ! Interface : adaptation nouveau modele
98  C            dans l'atmosphere avant l'advection en z    ! -------------------------------------
99        sqi = 0.  
100        sqf = 0.    ! Conversion du flux de masse en kg.s-1
101    
102        DO l = 1,llm    DO l = 1, llm
103           DO j = 1,jjp1      DO j = 1, jjp1
104              DO i = 1,iim        DO i = 1, iip1
105  cIM 240305            sqi = sqi + S0(i,j,l,9)          ! wgri (i,j,llm+1-l) =  w (i,j,l) / g
106                 sqi = sqi + S0(i,j,l,ntra)          wgri(i, j, llm+1-l) = w(i, j, l)
107              ENDDO          ! wgri (i,j,0) = 0.                ! a detruire ult.
108           ENDDO          ! wgri (i,j,l) = 0.1               !    w (i,j,l)
109        ENDDO          ! wgri (i,j,llm) = 0.              ! a detruire ult.
110        PRINT*,'-------- DIAG DANS ADVZ - ENTREE ---------'        END DO
111        PRINT*,'sqi=',sqi      END DO
112      END DO
113  C-----------------------------------------------------------------    DO j = 1, jjp1
114  C  Interface : adaptation nouveau modele      DO i = 1, iip1
115  C  -------------------------------------        wgri(i, j, 0) = 0.
116  C      END DO
117  C  Conversion du flux de masse en kg.s-1    END DO
118    
119        DO 500 l = 1,llm    ! -----------------------------------------------------------------
120           DO 500 j = 1,jjp1  
121              DO 500 i = 1,iip1      ! start here
122  c            wgri (i,j,llm+1-l) =  w (i,j,l) / g    ! boucle sur les latitudes
123                 wgri (i,j,llm+1-l) =  w (i,j,l)  
124  c             wgri (i,j,0) = 0.                ! a detruire ult.    DO k = 1, lat
125  c             wgri (i,j,l) = 0.1               !    w (i,j,l)  
126  c             wgri (i,j,llm) = 0.              ! a detruire ult.      ! place limits on appropriate moments before transport
127    500 CONTINUE      ! (if flux-limiting is to be applied)
128           DO  j = 1,jjp1  
129              DO i = 1,iip1        IF (.NOT. limit) GO TO 101
130                 wgri(i,j,0)=0.  
131              enddo      DO jv = 1, ntra
132           enddo        DO l = 1, niv
133            DO i = 1, lon
134  C-----------------------------------------------------------------            sz(i, k, l, jv) = sign(amin1(amax1(s0(i,k,l,jv), &
135                  0.),abs(sz(i,k,l,jv))), sz(i,k,l,jv))
136  C  start here                    END DO
137  C  boucle sur les latitudes        END DO
138  C      END DO
139        DO 1 K=1,LAT  
140  C  101 CONTINUE
141  C  place limits on appropriate moments before transport  
142  C      (if flux-limiting is to be applied)      ! boucle sur les niveaux intercouches de 1 a NIV-1
143  C      ! (flux nul au sommet L=0 et a la base L=NIV)
144        IF(.NOT.LIMIT) GO TO 101  
145  C      ! calculate flux and moments between adjacent boxes
146        DO 10 JV=1,NTRA      ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
147        DO 10 L=1,NIV      ! 1- create temporary moments/masses for partial boxes in transit
148           DO 100 I=1,LON      ! 2- reajusts moments remaining in the box
149              sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),  
150       +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))      DO l = 1, niv - 1
151   100     CONTINUE        lp = l + 1
152   10   CONTINUE  
153  C        DO i = 1, lon
154   101  CONTINUE  
155  C          IF (wgri(i,k,l)<0.) THEN
156  C  boucle sur les niveaux intercouches de 1 a NIV-1            fm(i, l) = -wgri(i, k, l)*dtz
157  C   (flux nul au sommet L=0 et a la base L=NIV)            alf(i) = fm(i, l)/sm(i, k, lp)
158  C            sm(i, k, lp) = sm(i, k, lp) - fm(i, l)
159  C  calculate flux and moments between adjacent boxes          ELSE
160  C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)            fm(i, l) = wgri(i, k, l)*dtz
161  C  1- create temporary moments/masses for partial boxes in transit            alf(i) = fm(i, l)/sm(i, k, l)
162  C  2- reajusts moments remaining in the box            sm(i, k, l) = sm(i, k, l) - fm(i, l)
163  C          END IF
164        DO 11 L=1,NIV-1  
165        LP=L+1          alfq(i) = alf(i)*alf(i)
166  C          alf1(i) = 1. - alf(i)
167        DO 110 I=1,LON          alf1q(i) = alf1(i)*alf1(i)
168  C  
169           IF(WGRI(I,K,L).LT.0.) THEN        END DO
170             FM(I,L)=-WGRI(I,K,L)*DTZ  
171             ALF(I)=FM(I,L)/SM(I,K,LP)        DO jv = 1, ntra
172             SM(I,K,LP)=SM(I,K,LP)-FM(I,L)          DO i = 1, lon
173           ELSE  
174             FM(I,L)=WGRI(I,K,L)*DTZ            IF (wgri(i,k,l)<0.) THEN
175             ALF(I)=FM(I,L)/SM(I,K,L)  
176             SM(I,K,L)=SM(I,K,L)-FM(I,L)              f0(i, l, jv) = alf(i)*(s0(i,k,lp,jv)-alf1(i)*sz(i,k,lp,jv))
177           ENDIF              fz(i, l, jv) = alfq(i)*sz(i, k, lp, jv)
178  C              fx(i, l, jv) = alf(i)*sx(i, k, lp, jv)
179           ALFQ (I)=ALF(I)*ALF(I)              fy(i, l, jv) = alf(i)*sy(i, k, lp, jv)
180           ALF1 (I)=1.-ALF(I)  
181           ALF1Q(I)=ALF1(I)*ALF1(I)              s0(i, k, lp, jv) = s0(i, k, lp, jv) - f0(i, l, jv)
182  C              sz(i, k, lp, jv) = alf1q(i)*sz(i, k, lp, jv)
183   110  CONTINUE              sx(i, k, lp, jv) = sx(i, k, lp, jv) - fx(i, l, jv)
184  C              sy(i, k, lp, jv) = sy(i, k, lp, jv) - fy(i, l, jv)
185        DO 111 JV=1,NTRA  
186        DO 1110 I=1,LON            ELSE
187  C  
188           IF(WGRI(I,K,L).LT.0.) THEN              f0(i, l, jv) = alf(i)*(s0(i,k,l,jv)+alf1(i)*sz(i,k,l,jv))
189  C              fz(i, l, jv) = alfq(i)*sz(i, k, l, jv)
190             F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )              fx(i, l, jv) = alf(i)*sx(i, k, l, jv)
191             FZ(I,L,JV)=ALFQ(I)*sz(I,K,LP,JV)              fy(i, l, jv) = alf(i)*sy(i, k, l, jv)
192             FX(I,L,JV)=ALF (I)*sx(I,K,LP,JV)  
193             FY(I,L,JV)=ALF (I)*sy(I,K,LP,JV)              s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, l, jv)
194  C              sz(i, k, l, jv) = alf1q(i)*sz(i, k, l, jv)
195             S0(I,K,LP,JV)=S0(I,K,LP,JV)-F0(I,L,JV)              sx(i, k, l, jv) = sx(i, k, l, jv) - fx(i, l, jv)
196             sz(I,K,LP,JV)=ALF1Q(I)*sz(I,K,LP,JV)              sy(i, k, l, jv) = sy(i, k, l, jv) - fy(i, l, jv)
197             sx(I,K,LP,JV)=sx(I,K,LP,JV)-FX(I,L,JV)  
198             sy(I,K,LP,JV)=sy(I,K,LP,JV)-FY(I,L,JV)            END IF
199  C  
200           ELSE          END DO
201  C        END DO
202             F0(I,L,JV)=ALF (I)*(S0(I,K,L,JV)+ALF1(I)*sz(I,K,L,JV) )  
203             FZ(I,L,JV)=ALFQ(I)*sz(I,K,L,JV)      END DO
204             FX(I,L,JV)=ALF (I)*sx(I,K,L,JV)  
205             FY(I,L,JV)=ALF (I)*sy(I,K,L,JV)      ! puts the temporary moments Fi into appropriate neighboring boxes
206  C  
207             S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,L,JV)      DO l = 1, niv - 1
208             sz(I,K,L,JV)=ALF1Q(I)*sz(I,K,L,JV)        lp = l + 1
209             sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,L,JV)  
210             sy(I,K,L,JV)=sy(I,K,L,JV)-FY(I,L,JV)        DO i = 1, lon
211  C  
212           ENDIF          IF (wgri(i,k,l)<0.) THEN
213  C            sm(i, k, l) = sm(i, k, l) + fm(i, l)
214   1110 CONTINUE            alf(i) = fm(i, l)/sm(i, k, l)
215   111  CONTINUE          ELSE
216  C            sm(i, k, lp) = sm(i, k, lp) + fm(i, l)
217   11   CONTINUE            alf(i) = fm(i, l)/sm(i, k, lp)
218  C          END IF
219  C  puts the temporary moments Fi into appropriate neighboring boxes  
220  C          alf1(i) = 1. - alf(i)
221        DO 12 L=1,NIV-1          alfq(i) = alf(i)*alf(i)
222        LP=L+1          alf1q(i) = alf1(i)*alf1(i)
223  C  
224        DO 120 I=1,LON        END DO
225  C  
226           IF(WGRI(I,K,L).LT.0.) THEN        DO jv = 1, ntra
227             SM(I,K,L)=SM(I,K,L)+FM(I,L)          DO i = 1, lon
228             ALF(I)=FM(I,L)/SM(I,K,L)  
229           ELSE            IF (wgri(i,k,l)<0.) THEN
230             SM(I,K,LP)=SM(I,K,LP)+FM(I,L)  
231             ALF(I)=FM(I,L)/SM(I,K,LP)              temptm = -alf(i)*s0(i, k, l, jv) + alf1(i)*f0(i, l, jv)
232           ENDIF              s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, l, jv)
233  C              sz(i, k, l, jv) = alf(i)*fz(i, l, jv) + alf1(i)*sz(i, k, l, jv) + &
234           ALF1(I)=1.-ALF(I)                3.*temptm
235           ALFQ(I)=ALF(I)*ALF(I)              sx(i, k, l, jv) = sx(i, k, l, jv) + fx(i, l, jv)
236           ALF1Q(I)=ALF1(I)*ALF1(I)              sy(i, k, l, jv) = sy(i, k, l, jv) + fy(i, l, jv)
237  C  
238   120  CONTINUE            ELSE
239  C  
240        DO 121 JV=1,NTRA              temptm = alf(i)*s0(i, k, lp, jv) - alf1(i)*f0(i, l, jv)
241        DO 1210 I=1,LON              s0(i, k, lp, jv) = s0(i, k, lp, jv) + f0(i, l, jv)
242  C              sz(i, k, lp, jv) = alf(i)*fz(i, l, jv) + &
243           IF(WGRI(I,K,L).LT.0.) THEN                alf1(i)*sz(i, k, lp, jv) + 3.*temptm
244  C              sx(i, k, lp, jv) = sx(i, k, lp, jv) + fx(i, l, jv)
245             TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)              sy(i, k, lp, jv) = sy(i, k, lp, jv) + fy(i, l, jv)
246             S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)  
247             sz(I,K,L,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,L,JV)+3.*TEMPTM            END IF
248             sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,L,JV)  
249             sy(I,K,L,JV)=sy(I,K,L,JV)+FY(I,L,JV)          END DO
250  C        END DO
251           ELSE  
252  C      END DO
253             TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)  
254             S0(I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)      ! fin de la boucle principale sur les latitudes
255             sz(I,K,LP,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,LP,JV)  
256       +                  +3.*TEMPTM    END DO
257             sx(I,K,LP,JV)=sx(I,K,LP,JV)+FX(I,L,JV)  
258             sy(I,K,LP,JV)=sy(I,K,LP,JV)+FY(I,L,JV)    ! *** ------------------- bouclage cyclique  en X ------------
259  C  
260           ENDIF    ! DO l = 1,llm
261  C    ! DO j = 1,jjp1
262   1210 CONTINUE    ! SM(iip1,j,l) = SM(1,j,l)
263   121  CONTINUE    ! S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
264  C    ! sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
265   12   CONTINUE    ! sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
266  C    ! sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
267  C  fin de la boucle principale sur les latitudes    ! ENDDO
268  C    ! ENDDO
269   1    CONTINUE  
270  C    ! -------------------------------------------------------------
271  C-------------------------------------------------------------    ! *** Test : diag de la qqtite totale de traceur
272  C    ! dans l'atmosphere avant l'advection en z
273  C ----------- AA Test en fin de ADVX ------ Controle des S*    DO l = 1, llm
274        DO j = 1, jjp1
275  c     DO 9999 l = 1, llm        DO i = 1, iim
276  c     DO 9999 j = 1, jjp1          ! IM 240305            sqf = sqf + S0(i,j,l,9)
277  c     DO 9999 i = 1, iip1          sqf = sqf + s0(i, j, l, ntra)
278  c        IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN        END DO
279  c           PRINT*, '-------------------'      END DO
280  c           PRINT*, 'En fin de ADVZ'    END DO
281  c           PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)    PRINT *, '-------- DIAG DANS ADVZ - SORTIE ---------'
282  c           print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)    PRINT *, 'sqf=', sqf
283  c           print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)  
284  c           print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)    ! -------------------------------------------------------------
285  c           WRITE (*,*) 'On arrete !! - pbl en fin de ADVZ1'    RETURN
286  c            STOP  END SUBROUTINE advz
287  c        ENDIF  ! _______________________________________________________________
288   9999 CONTINUE  ! _______________________________________________________________
   
 C *** ------------------- bouclage cyclique  en X ------------  
         
 c      DO l = 1,llm  
 c         DO j = 1,jjp1  
 c            SM(iip1,j,l) = SM(1,j,l)  
 c            S0(iip1,j,l,ntra) = S0(1,j,l,ntra)  
 C            sx(iip1,j,l,ntra) = sx(1,j,l,ntra)  
 c            sy(iip1,j,l,ntra) = sy(1,j,l,ntra)  
 c            sz(iip1,j,l,ntra) = sz(1,j,l,ntra)  
 c         ENDDO  
 c      ENDDO  
             
 C-------------------------------------------------------------  
 C *** Test : diag de la qqtite totale de traceur  
 C            dans l'atmosphere avant l'advection en z  
       DO l = 1,llm  
          DO j = 1,jjp1  
             DO i = 1,iim  
 cIM 240305            sqf = sqf + S0(i,j,l,9)  
                sqf = sqf + S0(i,j,l,ntra)  
             ENDDO  
          ENDDO  
       ENDDO  
       PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'  
       PRINT*,'sqf=', sqf  
   
 C-------------------------------------------------------------  
       RETURN  
       END  
 C_______________________________________________________________  
 C_______________________________________________________________  
Legend:
Removed from v.28  
changed lines
  Added in v.82
  ViewVC Help
Powered by ViewVC 1.1.21