/[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

trunk/libf/phylmd/Conflx/flxflux.f90 revision 52 by guez, Fri Sep 23 12:28:01 2011 UTC trunk/phylmd/Conflx/flxflux.f revision 254 by guez, Mon Feb 5 10:39:38 2018 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, mfu, pmfd, mfus, &
9        use SUPHEC_M         pmfds, mfuq, pmfdq, mful, 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, intent(in):: pap(klon, klev)
26        REAL ztmsmlt, zdelta, zqsat      REAL, intent(in):: 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 mfu(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, intent(inout):: mfus(klon, klev)
37        REAL prfl(klon), psfl(klon)      real pmfds(klon, klev)
38        REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)      REAL mfuq(klon, klev)
39        INTEGER  kcbot(klon), kctop(klon), ktype(klon)      REAL pmfdq(klon, klev)
40        LOGICAL  ldland(klon), ldcum(klon)      real mful(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, 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  !               mfus(i, k) = mfus(i, k) &
92           IF(lddraf(i).AND.k.GE.kdtop(i)) THEN                    - mfu(i, k) * (RCPD * ptenh(i, k) + pgeoh(i, k))
93              pmfds(i,k)=pmfds(i,k)-pmfd(i,k)* &               mfuq(i, k)=mfuq(i, k)-mfu(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  !               mfu(i, k)=0.
117        DO 130 k=ktopm2,klev               mfus(i, k)=0.
118        DO 125 i = 1, klon               mfuq(i, k)=0.
119        IF(ldcum(i).AND.k.GT.kcbot(i)) THEN               mful(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               mfu(i, k)=mfu(i, ikb)*zzp
138           pmflxr(i,k) = 0.0               mfus(i, k)=mfus(i, ikb)*zzp
139           pmflxs(i,k) = 0.0               mfuq(i, k)=mfuq(i, ikb)*zzp
140        ENDDO               mful(i, k)=mful(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                  zqsat = R2ES * FOEEW(ztmsmlt, RTT >= ztmsmlt) / pap(i, k)
164  !        panache descendant pour eliminer la negativite                  zqsat = MIN(0.5, zqsat)
165           IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN                  zqsat = zqsat / (1. - RETV * zqsat)
166              pdmfdp(i,k) = -pmflxr(i,k)-pmflxs(i,k)-pdmfup(i,k)                  pqsen(i, k) = zqsat
167              pmflxr(i,k+1) = 0.0               ENDIF
168              pmflxs(i,k+1) = 0.0               IF (pten(i, k) > RTT) THEN
169              pdpmel(i,k) = 0.0                  pmflxr(i, k + 1) = pmflxr(i, k) + pdmfup(i, k) + pdmfdp(i, k) &
170           ENDIF                       + pdpmel(i, k)
171        ENDIF                  pmflxs(i, k + 1)=pmflxs(i, k)-pdpmel(i, k)
172        ENDDO               ELSE
173        ENDDO                  pmflxs(i, k + 1)=pmflxs(i, k) + pdmfup(i, k) + pdmfdp(i, k)
174  !                  pmflxr(i, k + 1)=pmflxr(i, k)
175  !jq The new variable is initialized here.               ENDIF
176  !jq It contains the humidity which is fed to the downdraft               ! si la precipitation est negative, on ajuste le plux du
177  !jq by evaporation of precipitation in the column below the base               ! panache descendant pour eliminer la negativite
178  !jq of convection.               IF ((pmflxr(i, k + 1) + pmflxs(i, k + 1)) < 0.0) THEN
179  !jq                  pdmfdp(i, k) = -pmflxr(i, k)-pmflxs(i, k)-pdmfup(i, k)
180  !jq In the former version, this term has been subtracted from precip                  pmflxr(i, k + 1) = 0.0
181  !jq as well as the evaporation.                  pmflxs(i, k + 1) = 0.0
182  !jq                  pdpmel(i, k) = 0.0
183        DO k = 1, klev               ENDIF
184        DO i = 1, klon            ENDIF
185           maxpdmfdp(i,k)=0.0         ENDDO
186        ENDDO      ENDDO
187        ENDDO  
188        DO k = 1, klev      ! The new variable is initialized here.  It contains the
189        ! humidity which is fed to the downdraft by evaporation of
190        ! precipitation in the column below the base of convection.
191    
192        ! In the former version, this term has been subtracted from precip
193        ! as well as the evaporation.
194    
195        DO k = 1, klev
196           DO i = 1, klon
197              maxpdmfdp(i, k)=0.0
198           ENDDO
199        ENDDO
200        DO k = 1, klev
201         DO kp = k, klev         DO kp = k, klev
202          DO i = 1, klon            DO i = 1, klon
203           maxpdmfdp(i,k)=maxpdmfdp(i,k)+pdmfdp(i,kp)               maxpdmfdp(i, k)=maxpdmfdp(i, k) + pdmfdp(i, kp)
204          ENDDO            ENDDO
205           ENDDO
206        ENDDO
207        ! End of initialization
208    
209        DO k = ktopm2, klev
210           DO i = 1, klon
211              IF (ldcum(i) .AND. k >= kcbot(i)) THEN
212                 zrfl = pmflxr(i, k) + pmflxs(i, k)
213                 IF (zrfl > 1.0E-20) THEN
214                    zrnew = (MAX(0., SQRT(zrfl / zcucov) - CEVAPCU(k) &
215                         * (paph(i, k + 1) - paph(i, k)) &
216                         * MAX(0., pqsen(i, k) - pqen(i, k))))**2 * zcucov
217                    zrmin = zrfl - zcucov &
218                         * MAX(0., 0.8 * pqsen(i, k) - pqen(i, k)) &
219                         * zcons2 * (paph(i, k + 1) - paph(i, k))
220                    zrnew=MAX(zrnew, zrmin)
221                    zrfln=MAX(zrnew, 0.)
222                    zdrfl=MIN(0., zrfln-zrfl)
223                    ! At least the amount of precipiation needed to feed
224                    ! the downdraft with humidity below the base of
225                    ! convection has to be left and can't be evaporated
226                    ! (surely the evaporation can't be positive):
227                    zdrfl=MAX(zdrfl, &
228                         MIN(-pmflxr(i, k)-pmflxs(i, k)-maxpdmfdp(i, k), 0.0))
229    
230                    zdenom=1.0/MAX(1.0E-20, pmflxr(i, k) + pmflxs(i, k))
231                    IF (pten(i, k) > RTT) THEN
232                       zpdr = pdmfdp(i, k)
233                       zpds = 0.0
234                    ELSE
235                       zpdr = 0.0
236                       zpds = pdmfdp(i, k)
237                    ENDIF
238                    pmflxr(i, k + 1) = pmflxr(i, k) + zpdr + pdpmel(i, k) &
239                         + zdrfl*pmflxr(i, k)*zdenom
240                    pmflxs(i, k + 1) = pmflxs(i, k) + zpds - pdpmel(i, k) &
241                         + zdrfl*pmflxs(i, k)*zdenom
242                    pdmfup(i, k) = pdmfup(i, k) + zdrfl
243                 ELSE
244                    pmflxr(i, k + 1) = 0.0
245                    pmflxs(i, k + 1) = 0.0
246                    pdmfdp(i, k) = 0.0
247                    pdpmel(i, k) = 0.0
248                 ENDIF
249                 if (pmflxr(i, k) + pmflxs(i, k) < -1.e-26) &
250                      write(*, *) 'precip. < 1e-16 ', pmflxr(i, k) + pmflxs(i, k)
251              ENDIF
252         ENDDO         ENDDO
253        ENDDO      ENDDO
254  !jq End of initialization  
255  !      DO i = 1, klon
256        DO k = ktopm2, klev         prfl(i) = pmflxr(i, klev + 1)
257        DO i = 1, klon         psfl(i) = pmflxs(i, klev + 1)
258        IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN      end DO
259           zrfl = pmflxr(i,k) + pmflxs(i,k)  
260           IF (zrfl.GT.1.0E-20) THEN      RETURN
261              zrnew=(MAX(0.,SQRT(zrfl/zcucov)- &    END SUBROUTINE flxflux
262                    CEVAPCU(k)*(paph(i,k+1)-paph(i,k))* &  
263                    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.254

  ViewVC Help
Powered by ViewVC 1.1.21