/[lmdze]/trunk/phylmd/Conflx/flxflux.f
ViewVC logotype

Diff of /trunk/phylmd/Conflx/flxflux.f

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

revision 52 by guez, Fri Sep 23 12:28:01 2011 UTC revision 70 by guez, Mon Jun 24 15:39:52 2013 UTC
# Line 1  Line 1 
1        SUBROUTINE flxflux(pdtime, pqen, pqsen, ptenh, pqenh, pap &  module flxflux_m
2          ,  paph, ldland, pgeoh, kcbot, kctop, lddraf, kdtop &  
3          ,  ktype, ldcum, pmfu, pmfd, pmfus, pmfds &    IMPLICIT none
4          ,  pmfuq, pmfdq, pmful, plude, pdmfup, pdmfdp &  
5          ,  pten, prfl, psfl, pdpmel, ktopm2 &  contains
6          ,  pmflxr, pmflxs)  
7        use dimens_m    SUBROUTINE flxflux(pdtime, pqen, pqsen, ptenh, pqenh, pap, paph, ldland, &
8        use dimphy         pgeoh, kcbot, kctop, lddraf, kdtop, ktype, ldcum, pmfu, pmfd, pmfus, &
9        use SUPHEC_M         pmfds, pmfuq, pmfdq, pmful, plude, pdmfup, pdmfdp, pten, prfl, psfl, &
10        use yoethf_m         pdpmel, ktopm2, pmflxr, pmflxs)
11        use fcttre  
12              use yoecumf      ! This routine does the final calculation of convective fluxes in
13        IMPLICIT none      ! the cloud layer and in the subcloud layer.
14  !----------------------------------------------------------------------  
15  ! THIS ROUTINE DOES THE FINAL CALCULATION OF CONVECTIVE      USE dimphy, ONLY: klev, klon
16  ! FLUXES IN THE CLOUD LAYER AND IN THE SUBCLOUD LAYER      USE suphec_m, ONLY: rcpd, retv, rg, rlmlt, rtt
17  !----------------------------------------------------------------------      USE yoethf_m, ONLY: r2es
18  !      USE fcttre, ONLY: foeew
19        REAL cevapcu(klev)  
20  !     -----------------------------------------------------------------      REAL, intent(in):: pdtime
21        REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)      REAL, intent(in):: pqen(klon, klev)
22        REAL pten(klon,klev), ptenh(klon,klev)      real, intent(inout):: pqsen(klon, klev)
23        REAL paph(klon,klev+1), pgeoh(klon,klev)      REAL, intent(in):: ptenh(klon, klev)
24  !      REAL, intent(in):: pqenh(klon, klev)
25        REAL pap(klon,klev)      REAL pap(klon, klev)
26        REAL ztmsmlt, zdelta, zqsat      REAL paph(klon, klev+1)
27  !      LOGICAL ldland(klon)
28        REAL pmfu(klon,klev), pmfus(klon,klev)      real pgeoh(klon, klev)
29        REAL pmfd(klon,klev), pmfds(klon,klev)      INTEGER, intent(in):: kcbot(klon), kctop(klon)
30        REAL pmfuq(klon,klev), pmful(klon,klev)      LOGICAL lddraf(klon)
31        REAL pmfdq(klon,klev)      INTEGER kdtop(klon)
32        REAL plude(klon,klev)      INTEGER, intent(inout):: ktype(klon)
33        REAL pdmfup(klon,klev), pdpmel(klon,klev)      LOGICAL, intent(in):: ldcum(klon)
34  !jq The variable maxpdmfdp(klon) has been introduced by Olivier Boucher      REAL pmfu(klon, klev)
35  !jq 14/11/00 to fix the problem with the negative precipitation.      REAL pmfd(klon, klev)
36        REAL pdmfdp(klon,klev), maxpdmfdp(klon,klev)      real pmfus(klon, klev)
37        REAL prfl(klon), psfl(klon)      real pmfds(klon, klev)
38        REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)      REAL pmfuq(klon, klev)
39        INTEGER  kcbot(klon), kctop(klon), ktype(klon)      REAL pmfdq(klon, klev)
40        LOGICAL  ldland(klon), ldcum(klon)      real pmful(klon, klev)
41        INTEGER k, kp, i      REAL plude(klon, klev)
42        REAL zcons1, zcons2, zcucov, ztmelp2      REAL pdmfup(klon, klev)
43        REAL, intent(in):: pdtime      REAL pdmfdp(klon, klev)
44        real zdp, zzp, zfac, zsnmlt, zrfl, zrnew      REAL, intent(in):: pten(klon, klev)
45        REAL zrmin, zrfln, zdrfl      REAL prfl(klon), psfl(klon)
46        REAL zpds, zpdr, zdenom      real pdpmel(klon, klev)
47        INTEGER ktopm2, itop, ikb      INTEGER, intent(out):: ktopm2
48  !      REAL pmflxr(klon, klev+1), pmflxs(klon, klev+1)
49        LOGICAL lddraf(klon)  
50        INTEGER kdtop(klon)      ! Local:
51  !      REAL cevapcu(klev)
52  !      REAL ztmsmlt, zdelta, zqsat
53        DO 101 k=1,klev  
54        CEVAPCU(k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293) &      !jq The variable maxpdmfdp(klon) has been introduced by Olivier Boucher
55         *SQRT(0.5*(paph(1,k)+paph(1,k+1))/paph(1,klev+1)) ) * 0.5/RG      !jq 14/11/00 to fix the problem with the negative precipitation.
56   101  CONTINUE      real maxpdmfdp(klon, klev)
57  !  
58  ! SPECIFY CONSTANTS      INTEGER k, kp, i
59  !      REAL zcons1, zcons2, zcucov, ztmelp2
60        zcons1 = RCPD/(RLMLT*RG*pdtime)      real zdp, zzp, zfac, zsnmlt, zrfl, zrnew
61        zcons2 = 1./(RG*pdtime)      REAL zrmin, zrfln, zdrfl
62        zcucov = 0.05      REAL zpds, zpdr, zdenom
63        ztmelp2 = RTT + 2.      INTEGER ikb
64  !  
65  ! DETERMINE FINAL CONVECTIVE FLUXES      !---------------------------------------------------------------
66  !  
67        itop=klev      DO k=1, klev
68        DO 110 i = 1, klon         CEVAPCU(k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293) &
69           itop=MIN(itop,kctop(i))              *SQRT(0.5*(paph(1, k)+paph(1, k+1))/paph(1, klev+1))) * 0.5/RG
70           IF (.NOT.ldcum(i) .OR. kdtop(i).LT.kctop(i)) lddraf(i)=.FALSE.      end DO
71           IF(.NOT.ldcum(i)) ktype(i)=0  
72    110 CONTINUE      ! SPECIFY CONSTANTS
73  !  
74        ktopm2=itop-2      zcons1 = RCPD/(RLMLT*RG*pdtime)
75        DO 120 k=ktopm2,klev      zcons2 = 1./(RG*pdtime)
76        DO 115 i = 1, klon      zcucov = 0.05
77        IF(ldcum(i).AND.k.GE.kctop(i)-1) THEN      ztmelp2 = RTT + 2.
78           pmfus(i,k)=pmfus(i,k)-pmfu(i,k)* &  
79                        (RCPD*ptenh(i,k)+pgeoh(i,k))      ! DETERMINE FINAL CONVECTIVE FLUXES
80           pmfuq(i,k)=pmfuq(i,k)-pmfu(i,k)*pqenh(i,k)  
81           zdp = 1.5E4      DO i = 1, klon
82           IF ( ldland(i) ) zdp = 3.E4         IF (.NOT. ldcum(i) .OR. kdtop(i) < kctop(i)) lddraf(i)=.FALSE.
83  !         IF (.NOT. ldcum(i)) ktype(i) = 0
84  !        l'eau liquide detrainee est precipitee quand certaines      end DO
85  !        conditions sont reunies (sinon, elle est consideree  
86  !        evaporee dans l'environnement)      ktopm2 = min(klev, minval(kctop)) - 2
87  !  
88           IF(paph(i,kcbot(i))-paph(i,kctop(i)).GE.zdp.AND. &      DO k = ktopm2, klev
89              pqen(i,k-1).GT.0.8*pqsen(i,k-1)) &         DO i = 1, klon
90              pdmfup(i,k-1)=pdmfup(i,k-1)+plude(i,k-1)            IF (ldcum(i) .AND. k >= kctop(i) - 1) THEN
91  !               pmfus(i, k) = pmfus(i, k) &
92           IF(lddraf(i).AND.k.GE.kdtop(i)) THEN                    - pmfu(i, k) * (RCPD * ptenh(i, k) + pgeoh(i, k))
93              pmfds(i,k)=pmfds(i,k)-pmfd(i,k)* &               pmfuq(i, k)=pmfuq(i, k)-pmfu(i, k)*pqenh(i, k)
94                           (RCPD*ptenh(i,k)+pgeoh(i,k))               zdp = 1.5E4
95              pmfdq(i,k)=pmfdq(i,k)-pmfd(i,k)*pqenh(i,k)               IF (ldland(i)) zdp = 3.E4
96           ELSE  
97              pmfd(i,k)=0.               ! L'eau liquide détrainée est precipitée quand certaines
98              pmfds(i,k)=0.               ! conditions sont réunies (sinon, elle est considérée
99              pmfdq(i,k)=0.               ! évaporée dans l'environnement)
100              pdmfdp(i,k-1)=0.  
101           END IF               IF (paph(i, kcbot(i)) - paph(i, kctop(i)) >= zdp &
102        ELSE                    .AND. pqen(i, k - 1) > 0.8 * pqsen(i, k - 1)) &
103           pmfu(i,k)=0.                    pdmfup(i, k - 1) = pdmfup(i, k - 1) + plude(i, k - 1)
104           pmfus(i,k)=0.  
105           pmfuq(i,k)=0.               IF (lddraf(i).AND.k >= kdtop(i)) THEN
106           pmful(i,k)=0.                  pmfds(i, k)=pmfds(i, k)-pmfd(i, k)* &
107           pdmfup(i,k-1)=0.                       (RCPD*ptenh(i, k)+pgeoh(i, k))
108           plude(i,k-1)=0.                  pmfdq(i, k)=pmfdq(i, k)-pmfd(i, k)*pqenh(i, k)
109           pmfd(i,k)=0.               ELSE
110           pmfds(i,k)=0.                  pmfd(i, k)=0.
111           pmfdq(i,k)=0.                  pmfds(i, k)=0.
112           pdmfdp(i,k-1)=0.                  pmfdq(i, k)=0.
113        ENDIF                  pdmfdp(i, k-1)=0.
114    115 CONTINUE               END IF
115    120 CONTINUE            ELSE
116  !               pmfu(i, k)=0.
117        DO 130 k=ktopm2,klev               pmfus(i, k)=0.
118        DO 125 i = 1, klon               pmfuq(i, k)=0.
119        IF(ldcum(i).AND.k.GT.kcbot(i)) THEN               pmful(i, k)=0.
120           ikb=kcbot(i)               pdmfup(i, k-1)=0.
121           zzp=((paph(i,klev+1)-paph(i,k))/ &               plude(i, k-1)=0.
122                (paph(i,klev+1)-paph(i,ikb)))               pmfd(i, k)=0.
123           IF (ktype(i).EQ.3) zzp = zzp**2               pmfds(i, k)=0.
124           pmfu(i,k)=pmfu(i,ikb)*zzp               pmfdq(i, k)=0.
125           pmfus(i,k)=pmfus(i,ikb)*zzp               pdmfdp(i, k-1)=0.
126           pmfuq(i,k)=pmfuq(i,ikb)*zzp            ENDIF
127           pmful(i,k)=pmful(i,ikb)*zzp         end DO
128        ENDIF      end DO
129    125 CONTINUE  
130    130 CONTINUE      DO k=ktopm2, klev
131  !         DO i = 1, klon
132  ! CALCULATE RAIN/SNOW FALL RATES            IF (ldcum(i) .AND. k > kcbot(i)) THEN
133  ! CALCULATE MELTING OF SNOW               ikb=kcbot(i)
134  ! CALCULATE EVAPORATION OF PRECIP               zzp = ((paph(i, klev + 1) - paph(i, k)) &
135  !                    / (paph(i, klev + 1) - paph(i, ikb)))
136        DO k = 1, klev+1               IF (ktype(i) == 3) zzp = zzp**2
137        DO i = 1, klon               pmfu(i, k)=pmfu(i, ikb)*zzp
138           pmflxr(i,k) = 0.0               pmfus(i, k)=pmfus(i, ikb)*zzp
139           pmflxs(i,k) = 0.0               pmfuq(i, k)=pmfuq(i, ikb)*zzp
140        ENDDO               pmful(i, k)=pmful(i, ikb)*zzp
141        ENDDO            ENDIF
142        DO k = ktopm2, klev         end DO
143        DO i = 1, klon      end DO
144        IF (ldcum(i)) THEN  
145           IF (pmflxs(i,k).GT.0.0 .AND. pten(i,k).GT.ztmelp2) THEN      ! CALCULATE RAIN/SNOW FALL RATES
146              zfac=zcons1*(paph(i,k+1)-paph(i,k))      ! CALCULATE MELTING OF SNOW
147              zsnmlt=MIN(pmflxs(i,k),zfac*(pten(i,k)-ztmelp2))      ! CALCULATE EVAPORATION OF PRECIP
148              pdpmel(i,k)=zsnmlt  
149              ztmsmlt=pten(i,k)-zsnmlt/zfac      DO k = 1, klev+1
150              zdelta=MAX(0.,SIGN(1.,RTT-ztmsmlt))         DO i = 1, klon
151              zqsat=R2ES*FOEEW(ztmsmlt, zdelta) / pap(i,k)            pmflxr(i, k) = 0.0
152              zqsat=MIN(0.5,zqsat)            pmflxs(i, k) = 0.0
153              zqsat=zqsat/(1.-RETV  *zqsat)         ENDDO
154              pqsen(i,k) = zqsat      ENDDO
155           ENDIF      DO k = ktopm2, klev
156           IF (pten(i,k).GT.RTT) THEN         DO i = 1, klon
157           pmflxr(i,k+1)=pmflxr(i,k)+pdmfup(i,k)+pdmfdp(i,k)+pdpmel(i,k)            IF (ldcum(i)) THEN
158           pmflxs(i,k+1)=pmflxs(i,k)-pdpmel(i,k)               IF (pmflxs(i, k) > 0.0 .AND. pten(i, k) > ztmelp2) THEN
159           ELSE                  zfac=zcons1*(paph(i, k+1)-paph(i, k))
160             pmflxs(i,k+1)=pmflxs(i,k)+pdmfup(i,k)+pdmfdp(i,k)                  zsnmlt=MIN(pmflxs(i, k), zfac*(pten(i, k)-ztmelp2))
161             pmflxr(i,k+1)=pmflxr(i,k)                  pdpmel(i, k)=zsnmlt
162           ENDIF                  ztmsmlt=pten(i, k)-zsnmlt/zfac
163  !        si la precipitation est negative, on ajuste le plux du                  zdelta=MAX(0., SIGN(1., RTT-ztmsmlt))
164  !        panache descendant pour eliminer la negativite                  zqsat = R2ES * FOEEW(ztmsmlt, zdelta) / pap(i, k)
165           IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN                  zqsat = MIN(0.5, zqsat)
166              pdmfdp(i,k) = -pmflxr(i,k)-pmflxs(i,k)-pdmfup(i,k)                  zqsat = zqsat / (1. - RETV * zqsat)
167              pmflxr(i,k+1) = 0.0                  pqsen(i, k) = zqsat
168              pmflxs(i,k+1) = 0.0               ENDIF
169              pdpmel(i,k) = 0.0               IF (pten(i, k) > RTT) THEN
170           ENDIF                  pmflxr(i, k + 1) = pmflxr(i, k) + pdmfup(i, k) + pdmfdp(i, k) &
171        ENDIF                       + pdpmel(i, k)
172        ENDDO                  pmflxs(i, k+1)=pmflxs(i, k)-pdpmel(i, k)
173        ENDDO               ELSE
174  !                  pmflxs(i, k+1)=pmflxs(i, k)+pdmfup(i, k)+pdmfdp(i, k)
175  !jq The new variable is initialized here.                  pmflxr(i, k+1)=pmflxr(i, k)
176  !jq It contains the humidity which is fed to the downdraft               ENDIF
177  !jq by evaporation of precipitation in the column below the base               ! si la precipitation est negative, on ajuste le plux du
178  !jq of convection.               ! panache descendant pour eliminer la negativite
179  !jq               IF ((pmflxr(i, k+1)+pmflxs(i, k+1)) < 0.0) THEN
180  !jq In the former version, this term has been subtracted from precip                  pdmfdp(i, k) = -pmflxr(i, k)-pmflxs(i, k)-pdmfup(i, k)
181  !jq as well as the evaporation.                  pmflxr(i, k+1) = 0.0
182  !jq                  pmflxs(i, k+1) = 0.0
183        DO k = 1, klev                  pdpmel(i, k) = 0.0
184        DO i = 1, klon               ENDIF
185           maxpdmfdp(i,k)=0.0            ENDIF
186        ENDDO         ENDDO
187        ENDDO      ENDDO
188        DO k = 1, klev  
189        ! The new variable is initialized here.  It contains the
190        ! humidity which is fed to the downdraft by evaporation of
191        ! precipitation in the column below the base of convection.
192    
193        ! In the former version, this term has been subtracted from precip
194        ! as well as the evaporation.
195    
196        DO k = 1, klev
197           DO i = 1, klon
198              maxpdmfdp(i, k)=0.0
199           ENDDO
200        ENDDO
201        DO k = 1, klev
202         DO kp = k, klev         DO kp = k, klev
203          DO i = 1, klon            DO i = 1, klon
204           maxpdmfdp(i,k)=maxpdmfdp(i,k)+pdmfdp(i,kp)               maxpdmfdp(i, k)=maxpdmfdp(i, k)+pdmfdp(i, kp)
205          ENDDO            ENDDO
206           ENDDO
207        ENDDO
208        ! End of initialization
209    
210        DO k = ktopm2, klev
211           DO i = 1, klon
212              IF (ldcum(i) .AND. k >= kcbot(i)) THEN
213                 zrfl = pmflxr(i, k) + pmflxs(i, k)
214                 IF (zrfl > 1.0E-20) THEN
215                    zrnew = (MAX(0., SQRT(zrfl / zcucov) - CEVAPCU(k) &
216                         * (paph(i, k + 1) - paph(i, k)) &
217                         * MAX(0., pqsen(i, k) - pqen(i, k))))**2 * zcucov
218                    zrmin = zrfl - zcucov &
219                         * MAX(0., 0.8 * pqsen(i, k) - pqen(i, k)) &
220                         * zcons2 * (paph(i, k + 1) - paph(i, k))
221                    zrnew=MAX(zrnew, zrmin)
222                    zrfln=MAX(zrnew, 0.)
223                    zdrfl=MIN(0., zrfln-zrfl)
224                    ! At least the amount of precipiation needed to feed
225                    ! the downdraft with humidity below the base of
226                    ! convection has to be left and can't be evaporated
227                    ! (surely the evaporation can't be positive):
228                    zdrfl=MAX(zdrfl, &
229                         MIN(-pmflxr(i, k)-pmflxs(i, k)-maxpdmfdp(i, k), 0.0))
230    
231                    zdenom=1.0/MAX(1.0E-20, pmflxr(i, k)+pmflxs(i, k))
232                    IF (pten(i, k) > RTT) THEN
233                       zpdr = pdmfdp(i, k)
234                       zpds = 0.0
235                    ELSE
236                       zpdr = 0.0
237                       zpds = pdmfdp(i, k)
238                    ENDIF
239                    pmflxr(i, k+1) = pmflxr(i, k) + zpdr + pdpmel(i, k) &
240                         + zdrfl*pmflxr(i, k)*zdenom
241                    pmflxs(i, k+1) = pmflxs(i, k) + zpds - pdpmel(i, k) &
242                         + zdrfl*pmflxs(i, k)*zdenom
243                    pdmfup(i, k) = pdmfup(i, k) + zdrfl
244                 ELSE
245                    pmflxr(i, k+1) = 0.0
246                    pmflxs(i, k+1) = 0.0
247                    pdmfdp(i, k) = 0.0
248                    pdpmel(i, k) = 0.0
249                 ENDIF
250                 if (pmflxr(i, k) + pmflxs(i, k) < -1.e-26) &
251                      write(*, *) 'precip. < 1e-16 ', pmflxr(i, k) + pmflxs(i, k)
252              ENDIF
253         ENDDO         ENDDO
254        ENDDO      ENDDO
255  !jq End of initialization  
256  !      DO i = 1, klon
257        DO k = ktopm2, klev         prfl(i) = pmflxr(i, klev+1)
258        DO i = 1, klon         psfl(i) = pmflxs(i, klev+1)
259        IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN      end DO
260           zrfl = pmflxr(i,k) + pmflxs(i,k)  
261           IF (zrfl.GT.1.0E-20) THEN      RETURN
262              zrnew=(MAX(0.,SQRT(zrfl/zcucov)- &    END SUBROUTINE flxflux
263                    CEVAPCU(k)*(paph(i,k+1)-paph(i,k))* &  
264                    MAX(0.,pqsen(i,k)-pqen(i,k))))**2*zcucov  end module flxflux_m
             zrmin=zrfl-zcucov*MAX(0.,0.8*pqsen(i,k)-pqen(i,k)) &  
                   *zcons2*(paph(i,k+1)-paph(i,k))  
             zrnew=MAX(zrnew,zrmin)  
             zrfln=MAX(zrnew,0.)  
             zdrfl=MIN(0.,zrfln-zrfl)  
 !jq At least the amount of precipiation needed to feed the downdraft  
 !jq with humidity below the base of convection has to be left and can't  
 !jq be evaporated (surely the evaporation can't be positive):  
             zdrfl=MAX(zdrfl, &  
                   MIN(-pmflxr(i,k)-pmflxs(i,k)-maxpdmfdp(i,k),0.0))  
 !jq End of insertion  
 !  
             zdenom=1.0/MAX(1.0E-20,pmflxr(i,k)+pmflxs(i,k))  
             IF (pten(i,k).GT.RTT) THEN  
                zpdr = pdmfdp(i,k)  
                zpds = 0.0  
             ELSE  
                zpdr = 0.0  
                zpds = pdmfdp(i,k)  
             ENDIF  
             pmflxr(i,k+1) = pmflxr(i,k) + zpdr + pdpmel(i,k) &  
                           + zdrfl*pmflxr(i,k)*zdenom  
             pmflxs(i,k+1) = pmflxs(i,k) + zpds - pdpmel(i,k) &  
                           + zdrfl*pmflxs(i,k)*zdenom  
             pdmfup(i,k) = pdmfup(i,k) + zdrfl  
          ELSE  
             pmflxr(i,k+1) = 0.0  
             pmflxs(i,k+1) = 0.0  
             pdmfdp(i,k) = 0.0  
             pdpmel(i,k) = 0.0  
          ENDIF  
          if (pmflxr(i,k) + pmflxs(i,k).lt.-1.e-26)  &  
           write(*,*) 'precip. < 1e-16 ',pmflxr(i,k) + pmflxs(i,k)  
       ENDIF  
       ENDDO  
       ENDDO  
 !  
       DO 210 i = 1, klon  
          prfl(i) = pmflxr(i,klev+1)  
          psfl(i) = pmflxs(i,klev+1)  
   210 CONTINUE  
 !  
       RETURN  
       END  

Legend:
Removed from v.52  
changed lines
  Added in v.70

  ViewVC Help
Powered by ViewVC 1.1.21