/[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/dyn3d/advz.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/dyn3d/advz.f90 revision 81 by guez, Wed Mar 5 14:38:41 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 disvert_m    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 : diag de la qqtite totale de traceur    ! dans l'atmosphere avant l'advection en z
82  C            dans l'atmosphere avant l'advection en z    sqi = 0.
83        sqi = 0.    sqf = 0.
84        sqf = 0.  
85      DO l = 1, llm
86        DO l = 1,llm      DO j = 1, jjp1
87           DO j = 1,jjp1        DO i = 1, iim
88              DO i = 1,iim          ! IM 240305            sqi = sqi + S0(i,j,l,9)
89  cIM 240305            sqi = sqi + S0(i,j,l,9)          sqi = sqi + s0(i, j, l, ntra)
90                 sqi = sqi + S0(i,j,l,ntra)        END DO
91              ENDDO      END DO
92           ENDDO    END DO
93        ENDDO    PRINT *, '-------- DIAG DANS ADVZ - ENTREE ---------'
94        PRINT*,'-------- DIAG DANS ADVZ - ENTREE ---------'    PRINT *, 'sqi=', sqi
95        PRINT*,'sqi=',sqi  
96      ! -----------------------------------------------------------------
97  C-----------------------------------------------------------------    ! Interface : adaptation nouveau modele
98  C  Interface : adaptation nouveau modele    ! -------------------------------------
99  C  -------------------------------------  
100  C    ! Conversion du flux de masse en kg.s-1
101  C  Conversion du flux de masse en kg.s-1  
102      DO l = 1, llm
103        DO 500 l = 1,llm      DO j = 1, jjp1
104           DO 500 j = 1,jjp1        DO i = 1, iip1
105              DO 500 i = 1,iip1            ! wgri (i,j,llm+1-l) =  w (i,j,l) / g
106  c            wgri (i,j,llm+1-l) =  w (i,j,l) / g          wgri(i, j, llm+1-l) = w(i, j, l)
107                 wgri (i,j,llm+1-l) =  w (i,j,l)          ! wgri (i,j,0) = 0.                ! a detruire ult.
108  c             wgri (i,j,0) = 0.                ! a detruire ult.          ! wgri (i,j,l) = 0.1               !    w (i,j,l)
109  c             wgri (i,j,l) = 0.1               !    w (i,j,l)          ! wgri (i,j,llm) = 0.              ! a detruire ult.
110  c             wgri (i,j,llm) = 0.              ! a detruire ult.        END DO
111    500 CONTINUE      END DO
112           DO  j = 1,jjp1    END DO
113              DO i = 1,iip1      DO j = 1, jjp1
114                 wgri(i,j,0)=0.      DO i = 1, iip1
115              enddo        wgri(i, j, 0) = 0.
116           enddo      END DO
117      END DO
118  C-----------------------------------------------------------------  
119        ! -----------------------------------------------------------------
120  C  start here            
121  C  boucle sur les latitudes    ! start here
122  C    ! boucle sur les latitudes
123        DO 1 K=1,LAT  
124  C    DO k = 1, lat
125  C  place limits on appropriate moments before transport  
126  C      (if flux-limiting is to be applied)      ! place limits on appropriate moments before transport
127  C      ! (if flux-limiting is to be applied)
128        IF(.NOT.LIMIT) GO TO 101  
129  C      IF (.NOT. limit) GO TO 101
130        DO 10 JV=1,NTRA  
131        DO 10 L=1,NIV      DO jv = 1, ntra
132           DO 100 I=1,LON        DO l = 1, niv
133              sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),          DO i = 1, lon
134       +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))            sz(i, k, l, jv) = sign(amin1(amax1(s0(i,k,l,jv), &
135   100     CONTINUE              0.),abs(sz(i,k,l,jv))), sz(i,k,l,jv))
136   10   CONTINUE          END DO
137  C        END DO
138   101  CONTINUE      END DO
139  C  
140  C  boucle sur les niveaux intercouches de 1 a NIV-1  101 CONTINUE
141  C   (flux nul au sommet L=0 et a la base L=NIV)  
142  C      ! boucle sur les niveaux intercouches de 1 a NIV-1
143  C  calculate flux and moments between adjacent boxes      ! (flux nul au sommet L=0 et a la base L=NIV)
144  C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)  
145  C  1- create temporary moments/masses for partial boxes in transit      ! calculate flux and moments between adjacent boxes
146  C  2- reajusts moments remaining in the box      ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
147  C      ! 1- create temporary moments/masses for partial boxes in transit
148        DO 11 L=1,NIV-1      ! 2- reajusts moments remaining in the box
149        LP=L+1  
150  C      DO l = 1, niv - 1
151        DO 110 I=1,LON        lp = l + 1
152  C  
153           IF(WGRI(I,K,L).LT.0.) THEN        DO i = 1, lon
154             FM(I,L)=-WGRI(I,K,L)*DTZ  
155             ALF(I)=FM(I,L)/SM(I,K,LP)          IF (wgri(i,k,l)<0.) THEN
156             SM(I,K,LP)=SM(I,K,LP)-FM(I,L)            fm(i, l) = -wgri(i, k, l)*dtz
157           ELSE            alf(i) = fm(i, l)/sm(i, k, lp)
158             FM(I,L)=WGRI(I,K,L)*DTZ            sm(i, k, lp) = sm(i, k, lp) - fm(i, l)
159             ALF(I)=FM(I,L)/SM(I,K,L)          ELSE
160             SM(I,K,L)=SM(I,K,L)-FM(I,L)            fm(i, l) = wgri(i, k, l)*dtz
161           ENDIF            alf(i) = fm(i, l)/sm(i, k, l)
162  C            sm(i, k, l) = sm(i, k, l) - fm(i, l)
163           ALFQ (I)=ALF(I)*ALF(I)          END IF
164           ALF1 (I)=1.-ALF(I)  
165           ALF1Q(I)=ALF1(I)*ALF1(I)          alfq(i) = alf(i)*alf(i)
166  C          alf1(i) = 1. - alf(i)
167   110  CONTINUE          alf1q(i) = alf1(i)*alf1(i)
168  C  
169        DO 111 JV=1,NTRA        END DO
170        DO 1110 I=1,LON  
171  C        DO jv = 1, ntra
172           IF(WGRI(I,K,L).LT.0.) THEN          DO i = 1, lon
173  C  
174             F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )            IF (wgri(i,k,l)<0.) THEN
175             FZ(I,L,JV)=ALFQ(I)*sz(I,K,LP,JV)  
176             FX(I,L,JV)=ALF (I)*sx(I,K,LP,JV)              f0(i, l, jv) = alf(i)*(s0(i,k,lp,jv)-alf1(i)*sz(i,k,lp,jv))
177             FY(I,L,JV)=ALF (I)*sy(I,K,LP,JV)              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             S0(I,K,LP,JV)=S0(I,K,LP,JV)-F0(I,L,JV)              fy(i, l, jv) = alf(i)*sy(i, k, lp, jv)
180             sz(I,K,LP,JV)=ALF1Q(I)*sz(I,K,LP,JV)  
181             sx(I,K,LP,JV)=sx(I,K,LP,JV)-FX(I,L,JV)              s0(i, k, lp, jv) = s0(i, k, lp, jv) - f0(i, l, jv)
182             sy(I,K,LP,JV)=sy(I,K,LP,JV)-FY(I,L,JV)              sz(i, k, lp, jv) = alf1q(i)*sz(i, k, lp, jv)
183  C              sx(i, k, lp, jv) = sx(i, k, lp, jv) - fx(i, l, jv)
184           ELSE              sy(i, k, lp, jv) = sy(i, k, lp, jv) - fy(i, l, jv)
185  C  
186             F0(I,L,JV)=ALF (I)*(S0(I,K,L,JV)+ALF1(I)*sz(I,K,L,JV) )            ELSE
187             FZ(I,L,JV)=ALFQ(I)*sz(I,K,L,JV)  
188             FX(I,L,JV)=ALF (I)*sx(I,K,L,JV)              f0(i, l, jv) = alf(i)*(s0(i,k,l,jv)+alf1(i)*sz(i,k,l,jv))
189             FY(I,L,JV)=ALF (I)*sy(I,K,L,JV)              fz(i, l, jv) = alfq(i)*sz(i, k, l, jv)
190  C              fx(i, l, jv) = alf(i)*sx(i, k, l, jv)
191             S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,L,JV)              fy(i, l, jv) = alf(i)*sy(i, k, l, jv)
192             sz(I,K,L,JV)=ALF1Q(I)*sz(I,K,L,JV)  
193             sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,L,JV)              s0(i, k, l, jv) = s0(i, k, l, jv) - f0(i, l, jv)
194             sy(I,K,L,JV)=sy(I,K,L,JV)-FY(I,L,JV)              sz(i, k, l, jv) = alf1q(i)*sz(i, k, l, jv)
195  C              sx(i, k, l, jv) = sx(i, k, l, jv) - fx(i, l, jv)
196           ENDIF              sy(i, k, l, jv) = sy(i, k, l, jv) - fy(i, l, jv)
197  C  
198   1110 CONTINUE            END IF
199   111  CONTINUE  
200  C          END DO
201   11   CONTINUE        END DO
202  C  
203  C  puts the temporary moments Fi into appropriate neighboring boxes      END DO
204  C  
205        DO 12 L=1,NIV-1      ! puts the temporary moments Fi into appropriate neighboring boxes
206        LP=L+1  
207  C      DO l = 1, niv - 1
208        DO 120 I=1,LON        lp = l + 1
209  C  
210           IF(WGRI(I,K,L).LT.0.) THEN        DO i = 1, lon
211             SM(I,K,L)=SM(I,K,L)+FM(I,L)  
212             ALF(I)=FM(I,L)/SM(I,K,L)          IF (wgri(i,k,l)<0.) THEN
213           ELSE            sm(i, k, l) = sm(i, k, l) + fm(i, l)
214             SM(I,K,LP)=SM(I,K,LP)+FM(I,L)            alf(i) = fm(i, l)/sm(i, k, l)
215             ALF(I)=FM(I,L)/SM(I,K,LP)          ELSE
216           ENDIF            sm(i, k, lp) = sm(i, k, lp) + fm(i, l)
217  C            alf(i) = fm(i, l)/sm(i, k, lp)
218           ALF1(I)=1.-ALF(I)          END IF
219           ALFQ(I)=ALF(I)*ALF(I)  
220           ALF1Q(I)=ALF1(I)*ALF1(I)          alf1(i) = 1. - alf(i)
221  C          alfq(i) = alf(i)*alf(i)
222   120  CONTINUE          alf1q(i) = alf1(i)*alf1(i)
223  C  
224        DO 121 JV=1,NTRA        END DO
225        DO 1210 I=1,LON  
226  C        DO jv = 1, ntra
227           IF(WGRI(I,K,L).LT.0.) THEN          DO i = 1, lon
228  C  
229             TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)            IF (wgri(i,k,l)<0.) THEN
230             S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)  
231             sz(I,K,L,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,L,JV)+3.*TEMPTM              temptm = -alf(i)*s0(i, k, l, jv) + alf1(i)*f0(i, l, jv)
232             sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,L,JV)              s0(i, k, l, jv) = s0(i, k, l, jv) + f0(i, l, jv)
233             sy(I,K,L,JV)=sy(I,K,L,JV)+FY(I,L,JV)              sz(i, k, l, jv) = alf(i)*fz(i, l, jv) + alf1(i)*sz(i, k, l, jv) + &
234  C                3.*temptm
235           ELSE              sx(i, k, l, jv) = sx(i, k, l, jv) + fx(i, l, jv)
236  C              sy(i, k, l, jv) = sy(i, k, l, jv) + fy(i, l, jv)
237             TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)  
238             S0(I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)            ELSE
239             sz(I,K,LP,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,LP,JV)  
240       +                  +3.*TEMPTM              temptm = alf(i)*s0(i, k, lp, jv) - alf1(i)*f0(i, l, jv)
241             sx(I,K,LP,JV)=sx(I,K,LP,JV)+FX(I,L,JV)              s0(i, k, lp, jv) = s0(i, k, lp, jv) + f0(i, l, jv)
242             sy(I,K,LP,JV)=sy(I,K,LP,JV)+FY(I,L,JV)              sz(i, k, lp, jv) = alf(i)*fz(i, l, jv) + &
243  C                alf1(i)*sz(i, k, lp, jv) + 3.*temptm
244           ENDIF              sx(i, k, lp, jv) = sx(i, k, lp, jv) + fx(i, l, jv)
245  C              sy(i, k, lp, jv) = sy(i, k, lp, jv) + fy(i, l, jv)
246   1210 CONTINUE  
247   121  CONTINUE            END IF
248  C  
249   12   CONTINUE          END DO
250  C        END DO
251  C  fin de la boucle principale sur les latitudes  
252  C      END DO
253   1    CONTINUE  
254        ! fin de la boucle principale sur les latitudes
255  C *** ------------------- bouclage cyclique  en X ------------  
256            END DO
257  c      DO l = 1,llm  
258  c         DO j = 1,jjp1    ! *** ------------------- bouclage cyclique  en X ------------
259  c            SM(iip1,j,l) = SM(1,j,l)  
260  c            S0(iip1,j,l,ntra) = S0(1,j,l,ntra)    ! DO l = 1,llm
261  C            sx(iip1,j,l,ntra) = sx(1,j,l,ntra)    ! DO j = 1,jjp1
262  c            sy(iip1,j,l,ntra) = sy(1,j,l,ntra)    ! SM(iip1,j,l) = SM(1,j,l)
263  c            sz(iip1,j,l,ntra) = sz(1,j,l,ntra)    ! S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
264  c         ENDDO    ! sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
265  c      ENDDO    ! sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
266                ! sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
267  C-------------------------------------------------------------    ! ENDDO
268  C *** Test : diag de la qqtite totale de traceur    ! ENDDO
269  C            dans l'atmosphere avant l'advection en z  
270        DO l = 1,llm    ! -------------------------------------------------------------
271           DO j = 1,jjp1    ! *** Test : diag de la qqtite totale de traceur
272              DO i = 1,iim    ! dans l'atmosphere avant l'advection en z
273  cIM 240305            sqf = sqf + S0(i,j,l,9)    DO l = 1, llm
274                 sqf = sqf + S0(i,j,l,ntra)      DO j = 1, jjp1
275              ENDDO        DO i = 1, iim
276           ENDDO          ! IM 240305            sqf = sqf + S0(i,j,l,9)
277        ENDDO          sqf = sqf + s0(i, j, l, ntra)
278        PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'        END DO
279        PRINT*,'sqf=', sqf      END DO
280      END DO
281  C-------------------------------------------------------------    PRINT *, '-------- DIAG DANS ADVZ - SORTIE ---------'
282        RETURN    PRINT *, 'sqf=', sqf
283        END  
284  C_______________________________________________________________    ! -------------------------------------------------------------
285  C_______________________________________________________________    RETURN
286    END SUBROUTINE advz
287    ! _______________________________________________________________
288    ! _______________________________________________________________
Legend:
Removed from v.76  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21