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

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

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

trunk/libf/phylmd/conflx.f revision 51 by guez, Tue Sep 20 09:14:34 2011 UTC trunk/libf/phylmd/Conflx/conflx.f90 revision 62 by guez, Thu Jul 26 14:37:37 2012 UTC
# Line 1  Line 1 
1  !  module conflx_m
2  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/conflx.F,v 1.1.1.1 2004/05/19 12:53:08 lmdzadmin Exp $  
3  !    IMPLICIT none
4        SUBROUTINE conflx (dtime,pres_h,pres_f,  
5       e                   t, q, con_t, con_q, pqhfl, w,  contains
6       s                   d_t, d_q, rain, snow,  
7       s                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,    SUBROUTINE conflx (dtime, pres_h, pres_f, t, q, con_t, con_q, pqhfl, w, &
8       s                   kcbot, kctop, kdtop, pmflxr, pmflxs)         d_t, d_q, rain, snow, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, &
9  c         kctop, kdtop, pmflxr, pmflxs)
10        use dimens_m  
11        use dimphy      ! From LMDZ4/libf/phylmd/conflx.F, version 1.1.1.1 2004/05/19 12:53:08
12        use SUPHEC_M  
13        use yoethf_m      ! Author: Z. X. Li (LMD/CNRS)
14        use fcttre      ! date: 1994/10/14
15        IMPLICIT none  
16  c======================================================================      ! Objet: schéma flux de masse pour la convection (schéma de
17  c Auteur(s): Z.X. Li (LMD/CNRS) date: 19941014      ! Tiedtke avec quelques modifications mineures)
18  c Objet: Schema flux de masse pour la convection  
19  c        (schema de Tiedtke avec qqs modifications mineures)      ! Décembre 1997 : prise en compte des modifications introduites
20  c Dec.97: Prise en compte des modifications introduites par      ! par Olivier Boucher et Alexandre Armengaud pour mélange et
21  c         Olivier Boucher et Alexandre Armengaud pour melange      ! lessivage des traceurs passifs.
22  c         et lessivage des traceurs passifs.  
23  c======================================================================      use flxmain_m, only: flxmain
24  c Entree:      USE dimphy, ONLY: klev, klon
25        REAL, intent(in):: dtime            ! pas d'integration (s)      USE suphec_m, ONLY: rd, retv, rtt
26        REAL, intent(in):: pres_h(klon,klev+1) ! pression half-level (Pa)      USE yoethf_m, ONLY: r2es
27        REAL, intent(in):: pres_f(klon,klev)! pression full-level (Pa)      USE fcttre, ONLY: foeew
28        REAL, intent(in):: t(klon,klev)     ! temperature (K)  
29        REAL q(klon,klev)     ! humidite specifique (g/g)      ! Entree:
30        REAL w(klon,klev)     ! vitesse verticale (Pa/s)      REAL, intent(in):: dtime            ! pas d'integration (s)
31        REAL con_t(klon,klev) ! convergence de temperature (K/s)      REAL, intent(in):: pres_h(klon, klev+1) ! pression half-level (Pa)
32        REAL con_q(klon,klev) ! convergence de l'eau vapeur (g/g/s)      REAL, intent(in):: pres_f(klon, klev)! pression full-level (Pa)
33        REAL pqhfl(klon)      ! evaporation (negative vers haut) mm/s      REAL, intent(in):: t(klon, klev)     ! temperature (K)
34  c Sortie:      REAL q(klon, klev)     ! humidite specifique (g/g)
35        REAL d_t(klon,klev)   ! incrementation de temperature      REAL w(klon, klev)     ! vitesse verticale (Pa/s)
36        REAL d_q(klon,klev)   ! incrementation d'humidite      REAL con_t(klon, klev) ! convergence de temperature (K/s)
37        REAL pmfu(klon,klev)  ! flux masse (kg/m2/s) panache ascendant      REAL con_q(klon, klev) ! convergence de l'eau vapeur (g/g/s)
38        REAL pmfd(klon,klev)  ! flux masse (kg/m2/s) panache descendant      REAL pqhfl(klon)      ! evaporation (negative vers haut) mm/s
39        REAL pen_u(klon,klev)  
40        REAL pen_d(klon,klev)      ! Sortie:
41        REAL pde_u(klon,klev)      REAL d_t(klon, klev)   ! incrementation de temperature
42        REAL pde_d(klon,klev)      REAL d_q(klon, klev)   ! incrementation d'humidite
43        REAL rain(klon)       ! pluie (mm/s)  
44        REAL snow(klon)       ! neige (mm/s)      REAL, intent(out):: pmfu(:, :) ! (klon, klev)
45        REAL pmflxr(klon,klev+1)      ! flux masse (kg/m2/s) panache ascendant
46        REAL pmflxs(klon,klev+1)      
47        INTEGER kcbot(klon)  ! niveau du bas de la convection      REAL, intent(out):: pmfd(:, :) ! (klon, klev)
48        INTEGER kctop(klon)  ! niveau du haut de la convection      ! flux masse (kg/m2/s) panache descendant
49        INTEGER kdtop(klon)  ! niveau du haut des downdrafts  
50  c Local:      REAL pen_u(klon, klev)
51        REAL pt(klon,klev)      REAL pen_d(klon, klev)
52        REAL pq(klon,klev)      REAL pde_u(klon, klev)
53        REAL pqs(klon,klev)      REAL pde_d(klon, klev)
54        REAL pvervel(klon,klev)      REAL rain(klon)       ! pluie (mm/s)
55        LOGICAL land(klon)      REAL snow(klon)       ! neige (mm/s)
56  c      REAL pmflxr(klon, klev+1)
57        REAL d_t_bis(klon,klev)      REAL pmflxs(klon, klev+1)
58        REAL d_q_bis(klon,klev)      INTEGER kcbot(klon)  ! niveau du bas de la convection
59        REAL paprs(klon,klev+1)      INTEGER kctop(klon)  ! niveau du haut de la convection
60        REAL paprsf(klon,klev)      INTEGER kdtop(klon)  ! niveau du haut des downdrafts
61        REAL zgeom(klon,klev)  
62        REAL zcvgq(klon,klev)      ! Local:
63        REAL zcvgt(klon,klev)  
64  cAA      REAL pt(klon, klev)
65        REAL zmfu(klon,klev)      REAL pq(klon, klev)
66        REAL zmfd(klon,klev)      REAL pqs(klon, klev)
67        REAL zen_u(klon,klev)      REAL pvervel(klon, klev)
68        REAL zen_d(klon,klev)      LOGICAL land(klon)
69        REAL zde_u(klon,klev)  
70        REAL zde_d(klon,klev)      REAL d_t_bis(klon, klev)
71        REAL zmflxr(klon,klev+1)      REAL d_q_bis(klon, klev)
72        REAL zmflxs(klon,klev+1)      REAL paprs(klon, klev+1)
73  cAA      REAL paprsf(klon, klev)
74        REAL zgeom(klon, klev)
75  c      REAL zcvgq(klon, klev)
76        INTEGER i, k      REAL zcvgt(klon, klev)
77        REAL zdelta, zqsat  
78  c      REAL zmfu(klon, klev)
79  c      REAL zmfd(klon, klev)
80  c initialiser les variables de sortie (pour securite)      REAL zen_u(klon, klev)
81        DO i = 1, klon      REAL zen_d(klon, klev)
82           rain(i) = 0.0      REAL zde_u(klon, klev)
83           snow(i) = 0.0      REAL zde_d(klon, klev)
84           kcbot(i) = 0      REAL zmflxr(klon, klev+1)
85           kctop(i) = 0      REAL zmflxs(klon, klev+1)
86           kdtop(i) = 0  
87        ENDDO      INTEGER i, k
88        DO k = 1, klev      REAL zdelta, zqsat
89        DO i = 1, klon  
90           d_t(i,k) = 0.0      !--------------------------------------------------------------------
91           d_q(i,k) = 0.0  
92           pmfu(i,k) = 0.0      ! initialiser les variables de sortie (pour securite)
93           pmfd(i,k) = 0.0      DO i = 1, klon
94           pen_u(i,k) = 0.0         rain(i) = 0.0
95           pde_u(i,k) = 0.0         snow(i) = 0.0
96           pen_d(i,k) = 0.0         kcbot(i) = 0
97           pde_d(i,k) = 0.0         kctop(i) = 0
98           zmfu(i,k) = 0.0         kdtop(i) = 0
99           zmfd(i,k) = 0.0      ENDDO
100           zen_u(i,k) = 0.0      DO k = 1, klev
101           zde_u(i,k) = 0.0         DO i = 1, klon
102           zen_d(i,k) = 0.0            d_t(i, k) = 0.0
103           zde_d(i,k) = 0.0            d_q(i, k) = 0.0
104        ENDDO            pmfu(i, k) = 0.0
105        ENDDO            pmfd(i, k) = 0.0
106        DO k = 1, klev+1            pen_u(i, k) = 0.0
107        DO i = 1, klon            pde_u(i, k) = 0.0
108           zmflxr(i,k) = 0.0            pen_d(i, k) = 0.0
109           zmflxs(i,k) = 0.0            pde_d(i, k) = 0.0
110        ENDDO            zmfu(i, k) = 0.0
111        ENDDO            zmfd(i, k) = 0.0
112  c            zen_u(i, k) = 0.0
113  c calculer la nature du sol (pour l'instant, ocean partout)            zde_u(i, k) = 0.0
114        DO i = 1, klon            zen_d(i, k) = 0.0
115           land(i) = .FALSE.            zde_d(i, k) = 0.0
116        ENDDO         ENDDO
117  c      ENDDO
118  c preparer les variables d'entree (attention: l'ordre des niveaux      DO k = 1, klev+1
119  c verticaux augmente du haut vers le bas)         DO i = 1, klon
120        DO k = 1, klev            zmflxr(i, k) = 0.0
121        DO i = 1, klon            zmflxs(i, k) = 0.0
          pt(i,k) = t(i,klev-k+1)  
          pq(i,k) = q(i,klev-k+1)  
          paprsf(i,k) = pres_f(i,klev-k+1)  
          paprs(i,k) = pres_h(i,klev+1-k+1)  
          pvervel(i,k) = w(i,klev+1-k)  
          zcvgt(i,k) = con_t(i,klev-k+1)  
          zcvgq(i,k) = con_q(i,klev-k+1)  
 c  
          zdelta=MAX(0.,SIGN(1.,RTT-pt(i,k)))  
          zqsat=R2ES*FOEEW ( pt(i,k), zdelta ) / paprsf(i,k)  
          zqsat=MIN(0.5,zqsat)  
          zqsat=zqsat/(1.-RETV  *zqsat)  
          pqs(i,k) = zqsat  
       ENDDO  
       ENDDO  
       DO i = 1, klon  
          paprs(i,klev+1) = pres_h(i,1)  
          zgeom(i,klev) = RD * pt(i,klev)  
      .                   / (0.5*(paprs(i,klev+1)+paprsf(i,klev)))  
      .                   * (paprs(i,klev+1)-paprsf(i,klev))  
       ENDDO  
       DO k = klev-1, 1, -1  
       DO i = 1, klon  
          zgeom(i,k) = zgeom(i,k+1)  
      .              + RD * 0.5*(pt(i,k+1)+pt(i,k)) / paprs(i,k+1)  
      .                   * (paprsf(i,k+1)-paprsf(i,k))  
       ENDDO  
       ENDDO  
 c  
 c appeler la routine principale  
 c  
       CALL flxmain(dtime, pt, pq, pqs, pqhfl,  
      .             paprsf, paprs, zgeom, land, zcvgt, zcvgq, pvervel,  
      .             rain, snow, kcbot, kctop, kdtop,  
      .             zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,  
      .             d_t_bis, d_q_bis, zmflxr, zmflxs)  
 C  
 cAA--------------------------------------------------------  
 cAA rem : De la meme facon que l'on effectue le reindicage  
 cAA       pour la temperature t et le champ q  
 cAA       on reindice les flux necessaires a la convection  
 cAA       des traceurs  
 cAA--------------------------------------------------------  
       DO k = 1, klev  
       DO i = 1, klon  
          d_q(i,klev+1-k) = dtime*d_q_bis(i,k)  
          d_t(i,klev+1-k) = dtime*d_t_bis(i,k)  
       ENDDO  
       ENDDO  
 c  
       DO i = 1, klon  
          pmfu(i,1)= 0.  
          pmfd(i,1)= 0.  
          pen_d(i,1)= 0.  
          pde_d(i,1)= 0.  
       ENDDO  
       
       DO k = 2, klev  
       DO i = 1, klon  
          pmfu(i,klev+2-k)= zmfu(i,k)  
          pmfd(i,klev+2-k)= zmfd(i,k)  
       ENDDO  
       ENDDO  
 c  
       DO k = 1, klev  
       DO i = 1, klon  
          pen_u(i,klev+1-k)=  zen_u(i,k)  
          pde_u(i,klev+1-k)=  zde_u(i,k)  
       ENDDO  
       ENDDO  
 c  
       DO k = 1, klev-1  
       DO i = 1, klon  
          pen_d(i,klev+1-k)= -zen_d(i,k+1)  
          pde_d(i,klev+1-k)= -zde_d(i,k+1)  
       ENDDO  
       ENDDO  
   
       DO k = 1, klev+1  
       DO i = 1, klon  
          pmflxr(i,klev+2-k)= zmflxr(i,k)  
          pmflxs(i,klev+2-k)= zmflxs(i,k)  
       ENDDO  
       ENDDO  
   
       RETURN  
       END  
 c--------------------------------------------------------------------  
       SUBROUTINE flxmain(pdtime, pten, pqen, pqsen, pqhfl, pap, paph,  
      .                   pgeo, ldland, ptte, pqte, pvervel,  
      .                   prsfc, pssfc, kcbot, kctop, kdtop,  
 c     *                   ldcum, ktype,  
      .                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,  
      .                   dt_con, dq_con, pmflxr, pmflxs)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
             use yoecumf  
       IMPLICIT none  
 C     ------------------------------------------------------------------  
 C     ----------------------------------------------------------------  
       REAL pten(klon,klev), pqen(klon,klev), pqsen(klon,klev)  
       REAL ptte(klon,klev)  
       REAL pqte(klon,klev)  
       REAL pvervel(klon,klev)  
       REAL pgeo(klon,klev), pap(klon,klev), paph(klon,klev+1)  
       REAL pqhfl(klon)  
 c  
       REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)  
       REAL plude(klon,klev)  
       REAL pmfu(klon,klev)  
       REAL prsfc(klon), pssfc(klon)  
       INTEGER  kcbot(klon), kctop(klon), ktype(klon)  
       LOGICAL  ldland(klon), ldcum(klon)  
 c  
       REAL ztenh(klon,klev), zqenh(klon,klev), zqsenh(klon,klev)  
       REAL zgeoh(klon,klev)  
       REAL zmfub(klon), zmfub1(klon)  
       REAL zmfus(klon,klev), zmfuq(klon,klev), zmful(klon,klev)  
       REAL zdmfup(klon,klev), zdpmel(klon,klev)  
       REAL zentr(klon), zhcbase(klon)  
       REAL zdqpbl(klon), zdqcv(klon), zdhpbl(klon)  
       REAL zrfl(klon)  
       REAL pmflxr(klon,klev+1)  
       REAL pmflxs(klon,klev+1)  
       INTEGER  ilab(klon,klev), ictop0(klon)  
       LOGICAL  llo1  
       REAL dt_con(klon,klev), dq_con(klon,klev)  
       REAL zmfmax, zdh  
       REAL, intent(in):: pdtime  
       real zqumqe, zdqmin, zalvdcp, zhsat, zzz  
       REAL zhhat, zpbmpt, zgam, zeps, zfac  
       INTEGER i, k, ikb, itopm2, kcum  
 c  
       REAL pen_u(klon,klev), pde_u(klon,klev)  
       REAL pen_d(klon,klev), pde_d(klon,klev)  
 c  
       REAL ptd(klon,klev), pqd(klon,klev), pmfd(klon,klev)  
       REAL zmfds(klon,klev), zmfdq(klon,klev), zdmfdp(klon,klev)  
       INTEGER kdtop(klon)  
       LOGICAL lddraf(klon)  
 C---------------------------------------------------------------------  
       LOGICAL firstcal  
       SAVE firstcal  
       DATA firstcal / .TRUE. /  
 C---------------------------------------------------------------------  
       IF (firstcal) THEN  
          CALL flxsetup  
          firstcal = .FALSE.  
       ENDIF  
 C---------------------------------------------------------------------  
       DO i = 1, klon  
          ldcum(i) = .FALSE.  
       ENDDO  
       DO k = 1, klev  
       DO i = 1, klon  
          dt_con(i,k) = 0.0  
          dq_con(i,k) = 0.0  
       ENDDO  
       ENDDO  
 c----------------------------------------------------------------------  
 c initialiser les variables et faire l'interpolation verticale  
 c----------------------------------------------------------------------  
       CALL flxini(pten, pqen, pqsen, pgeo,  
      .     paph, zgeoh, ztenh, zqenh, zqsenh,  
      .     ptu, pqu, ptd, pqd, pmfd, zmfds, zmfdq, zdmfdp,  
      .     pmfu, zmfus, zmfuq, zdmfup,  
      .     zdpmel, plu, plude, ilab, pen_u, pde_u, pen_d, pde_d)  
 c---------------------------------------------------------------------  
 c determiner les valeurs au niveau de base de la tour convective  
 c---------------------------------------------------------------------  
       CALL flxbase(ztenh, zqenh, zgeoh, paph,  
      *            ptu, pqu, plu, ldcum, kcbot, ilab)  
 c---------------------------------------------------------------------  
 c calculer la convergence totale de l'humidite et celle en provenance  
 c de la couche limite, plus precisement, la convergence integree entre  
 c le sol et la base de la convection. Cette derniere convergence est  
 c comparee avec l'evaporation obtenue dans la couche limite pour  
 c determiner le type de la convection  
 c---------------------------------------------------------------------  
       k=1  
       DO i = 1, klon  
          zdqcv(i) = pqte(i,k)*(paph(i,k+1)-paph(i,k))  
          zdhpbl(i) = 0.0  
          zdqpbl(i) = 0.0  
       ENDDO  
 c  
       DO k=2,klev  
       DO i = 1, klon  
           zdqcv(i)=zdqcv(i)+pqte(i,k)*(paph(i,k+1)-paph(i,k))  
           IF (k.GE.kcbot(i)) THEN  
              zdqpbl(i)=zdqpbl(i)+pqte(i,k)*(paph(i,k+1)-paph(i,k))  
              zdhpbl(i)=zdhpbl(i)+(RCPD*ptte(i,k)+RLVTT*pqte(i,k))  
      .                          *(paph(i,k+1)-paph(i,k))  
           ENDIF  
       ENDDO  
       ENDDO  
 c  
       DO i = 1, klon  
          ktype(i) = 2  
          if (zdqcv(i).GT.MAX(0.,-1.5*pqhfl(i)*RG)) ktype(i) = 1  
 ccc         if (zdqcv(i).GT.MAX(0.,-1.1*pqhfl(i)*RG)) ktype(i) = 1  
       ENDDO  
 c  
 c---------------------------------------------------------------------  
 c determiner le flux de masse entrant a travers la base.  
 c on ignore, pour l'instant, l'effet du panache descendant  
 c---------------------------------------------------------------------  
       DO i = 1, klon  
          ikb=kcbot(i)  
          zqumqe=pqu(i,ikb)+plu(i,ikb)-zqenh(i,ikb)  
          zdqmin=MAX(0.01*zqenh(i,ikb),1.E-10)  
          IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)) THEN  
             zmfub(i) = zdqpbl(i)/(RG*MAX(zqumqe,zdqmin))  
          ELSE  
             zmfub(i) = 0.01  
             ldcum(i)=.FALSE.  
          ENDIF  
          IF (ktype(i).EQ.2) THEN  
             zdh = RCPD*(ptu(i,ikb)-ztenh(i,ikb)) + RLVTT*zqumqe  
             zdh = RG * MAX(zdh,1.0E5*zdqmin)  
             IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub(i)=zdhpbl(i)/zdh  
          ENDIF  
          zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)  
          zmfub(i) = MIN(zmfub(i),zmfmax)  
          zentr(i) = ENTRSCV  
          IF (ktype(i).EQ.1) zentr(i) = ENTRPEN  
       ENDDO  
 C-----------------------------------------------------------------------  
 C DETERMINE CLOUD ASCENT FOR ENTRAINING PLUME  
 C-----------------------------------------------------------------------  
 c (A) calculer d'abord la hauteur "theorique" de la tour convective sans  
 c     considerer l'entrainement ni le detrainement du panache, sachant  
 c     ces derniers peuvent abaisser la hauteur theorique.  
 c  
       DO i = 1, klon  
          ikb=kcbot(i)  
          zhcbase(i)=RCPD*ptu(i,ikb)+zgeoh(i,ikb)+RLVTT*pqu(i,ikb)  
          ictop0(i)=kcbot(i)-1  
       ENDDO  
 c  
       zalvdcp=RLVTT/RCPD  
       DO k=klev-1,3,-1  
       DO i = 1, klon  
          zhsat=RCPD*ztenh(i,k)+zgeoh(i,k)+RLVTT*zqsenh(i,k)  
          zgam=R5LES*zalvdcp*zqsenh(i,k)/  
      .        ((1.-RETV  *zqsenh(i,k))*(ztenh(i,k)-R4LES)**2)  
          zzz=RCPD*ztenh(i,k)*0.608  
          zhhat=zhsat-(zzz+zgam*zzz)/(1.+zgam*zzz/RLVTT)*  
      .               MAX(zqsenh(i,k)-zqenh(i,k),0.)  
          IF(k.LT.ictop0(i).AND.zhcbase(i).GT.zhhat) ictop0(i)=k  
       ENDDO  
       ENDDO  
 c  
 c (B) calculer le panache ascendant  
 c  
       CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,  
      .     pgeo, zgeoh, pap, paph, pqte, pvervel,  
      .     ldland, ldcum, ktype, ilab,  
      .     ptu, pqu, plu, pmfu, zmfub, zentr,  
      .     zmfus, zmfuq, zmful, plude, zdmfup,  
      .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)  
       IF (kcum.EQ.0) GO TO 1000  
 C  
 C verifier l'epaisseur de la convection et changer eventuellement  
 c le taux d'entrainement/detrainement  
 C  
       DO i = 1, klon  
          zpbmpt=paph(i,kcbot(i))-paph(i,kctop(i))  
          IF(ldcum(i).AND.ktype(i).EQ.1.AND.zpbmpt.LT.2.E4)ktype(i)=2  
          IF(ldcum(i)) ictop0(i)=kctop(i)  
          IF(ktype(i).EQ.2) zentr(i)=ENTRSCV  
       ENDDO  
 c  
       IF (lmfdd) THEN  ! si l'on considere le panache descendant  
 c  
 c calculer la precipitation issue du panache ascendant pour  
 c determiner l'existence du panache descendant dans la convection  
       DO i = 1, klon  
          zrfl(i)=zdmfup(i,1)  
       ENDDO  
       DO k=2,klev  
       DO i = 1, klon  
          zrfl(i)=zrfl(i)+zdmfup(i,k)  
       ENDDO  
       ENDDO  
 c  
 c determiner le LFS (level of free sinking: niveau de plonge libre)  
       CALL flxdlfs(ztenh, zqenh, zgeoh, paph, ptu, pqu,  
      *     ldcum,    kcbot,    kctop,    zmfub,    zrfl,  
      *     ptd,      pqd,  
      *     pmfd,     zmfds,    zmfdq,    zdmfdp,  
      *     kdtop,    lddraf)  
 c  
 c calculer le panache descendant  
       CALL flxddraf(ztenh,    zqenh,  
      *     zgeoh,    paph,     zrfl,  
      *     ptd,      pqd,  
      *     pmfd,     zmfds,    zmfdq,    zdmfdp,  
      *     lddraf, pen_d, pde_d)  
 c  
 c calculer de nouveau le flux de masse entrant a travers la base  
 c de la convection, sachant qu'il a ete modifie par le panache  
 c descendant  
       DO i = 1, klon  
       IF (lddraf(i)) THEN  
          ikb = kcbot(i)  
          llo1 = PMFD(i,ikb).LT.0.  
          zeps = 0.  
          IF ( llo1 ) zeps = CMFDEPS  
          zqumqe = pqu(i,ikb)+plu(i,ikb)-  
      .            zeps*pqd(i,ikb)-(1.-zeps)*zqenh(i,ikb)  
          zdqmin = MAX(0.01*zqenh(i,ikb),1.E-10)  
          zmfmax = (paph(i,ikb)-paph(i,ikb-1)) / (RG*pdtime)  
          IF (zdqpbl(i).GT.0..AND.zqumqe.GT.zdqmin.AND.ldcum(i)  
      .       .AND.zmfub(i).LT.zmfmax) THEN  
             zmfub1(i) = zdqpbl(i) / (RG*MAX(zqumqe,zdqmin))  
          ELSE  
             zmfub1(i) = zmfub(i)  
          ENDIF  
          IF (ktype(i).EQ.2) THEN  
             zdh = RCPD*(ptu(i,ikb)-zeps*ptd(i,ikb)-  
      .            (1.-zeps)*ztenh(i,ikb))+RLVTT*zqumqe  
             zdh = RG * MAX(zdh,1.0E5*zdqmin)  
             IF (zdhpbl(i).GT.0..AND.ldcum(i))zmfub1(i)=zdhpbl(i)/zdh  
          ENDIF  
          IF ( .NOT.((ktype(i).EQ.1.OR.ktype(i).EQ.2).AND.  
      .              ABS(zmfub1(i)-zmfub(i)).LT.0.2*zmfub(i)) )  
      .      zmfub1(i) = zmfub(i)  
       ENDIF  
       ENDDO  
       DO k = 1, klev  
       DO i = 1, klon  
       IF (lddraf(i)) THEN  
          zfac = zmfub1(i)/MAX(zmfub(i),1.E-10)  
          pmfd(i,k) = pmfd(i,k)*zfac  
          zmfds(i,k) = zmfds(i,k)*zfac  
          zmfdq(i,k) = zmfdq(i,k)*zfac  
          zdmfdp(i,k) = zdmfdp(i,k)*zfac  
          pen_d(i,k) = pen_d(i,k)*zfac  
          pde_d(i,k) = pde_d(i,k)*zfac  
       ENDIF  
       ENDDO  
       ENDDO  
       DO i = 1, klon  
          IF (lddraf(i)) zmfub(i)=zmfub1(i)  
       ENDDO  
 c  
       ENDIF   ! fin de test sur lmfdd  
 c  
 c-----------------------------------------------------------------------  
 c calculer de nouveau le panache ascendant  
 c-----------------------------------------------------------------------  
       CALL flxasc(pdtime,ztenh, zqenh, pten, pqen, pqsen,  
      .     pgeo, zgeoh, pap, paph, pqte, pvervel,  
      .     ldland, ldcum, ktype, ilab,  
      .     ptu, pqu, plu, pmfu, zmfub, zentr,  
      .     zmfus, zmfuq, zmful, plude, zdmfup,  
      .     kcbot, kctop, ictop0, kcum, pen_u, pde_u)  
 c  
 c-----------------------------------------------------------------------  
 c determiner les flux convectifs en forme finale, ainsi que  
 c la quantite des precipitations  
 c-----------------------------------------------------------------------  
       CALL flxflux(pdtime, pqen, pqsen, ztenh, zqenh, pap, paph,  
      .     ldland, zgeoh, kcbot, kctop, lddraf, kdtop, ktype, ldcum,  
      .     pmfu, pmfd, zmfus, zmfds, zmfuq, zmfdq, zmful, plude,  
      .     zdmfup, zdmfdp, pten, prsfc, pssfc, zdpmel, itopm2,  
      .     pmflxr, pmflxs)  
 c  
 c----------------------------------------------------------------------  
 c calculer les tendances pour T et Q  
 c----------------------------------------------------------------------  
       CALL flxdtdq(itopm2, paph, ldcum, pten,  
      e     zmfus, zmfds, zmfuq, zmfdq, zmful, zdmfup, zdmfdp, zdpmel,  
      s     dt_con,dq_con)  
 c  
  1000 CONTINUE  
       RETURN  
       END  
       SUBROUTINE flxini(pten, pqen, pqsen, pgeo, paph, pgeoh, ptenh,  
      .           pqenh, pqsenh, ptu, pqu, ptd, pqd, pmfd, pmfds, pmfdq,  
      .           pdmfdp, pmfu, pmfus, pmfuq, pdmfup, pdpmel, plu, plude,  
      .           klab,pen_u, pde_u, pen_d, pde_d)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
       IMPLICIT none  
 C----------------------------------------------------------------------  
 C THIS ROUTINE INTERPOLATES LARGE-SCALE FIELDS OF T,Q ETC.  
 C TO HALF LEVELS (I.E. GRID FOR MASSFLUX SCHEME),  
 C AND INITIALIZES VALUES FOR UPDRAFTS  
 C----------------------------------------------------------------------  
 C  
       REAL pten(klon,klev)   ! temperature (environnement)  
       REAL pqen(klon,klev)   ! humidite (environnement)  
       REAL pqsen(klon,klev)  ! humidite saturante (environnement)  
       REAL pgeo(klon,klev)   ! geopotentiel (g * metre)  
       REAL pgeoh(klon,klev)  ! geopotentiel aux demi-niveaux  
       REAL paph(klon,klev+1) ! pression aux demi-niveaux  
       REAL ptenh(klon,klev)  ! temperature aux demi-niveaux  
       REAL pqenh(klon,klev)  ! humidite aux demi-niveaux  
       REAL pqsenh(klon,klev) ! humidite saturante aux demi-niveaux  
 C  
       REAL ptu(klon,klev)    ! temperature du panache ascendant (p-a)  
       REAL pqu(klon,klev)    ! humidite du p-a  
       REAL plu(klon,klev)    ! eau liquide du p-a  
       REAL pmfu(klon,klev)   ! flux de masse du p-a  
       REAL pmfus(klon,klev)  ! flux de l'energie seche dans le p-a  
       REAL pmfuq(klon,klev)  ! flux de l'humidite dans le p-a  
       REAL pdmfup(klon,klev) ! quantite de l'eau precipitee dans p-a  
       REAL plude(klon,klev)  ! quantite de l'eau liquide jetee du  
 c                              p-a a l'environnement  
       REAL pdpmel(klon,klev) ! quantite de neige fondue  
 c  
       REAL ptd(klon,klev)    ! temperature du panache descendant (p-d)  
       REAL pqd(klon,klev)    ! humidite du p-d  
       REAL pmfd(klon,klev)   ! flux de masse du p-d  
       REAL pmfds(klon,klev)  ! flux de l'energie seche dans le p-d  
       REAL pmfdq(klon,klev)  ! flux de l'humidite dans le p-d  
       REAL pdmfdp(klon,klev) ! quantite de precipitation dans p-d  
 c  
       REAL pen_u(klon,klev) ! quantite de masse entrainee pour p-a  
       REAL pde_u(klon,klev) ! quantite de masse detrainee pour p-a  
       REAL pen_d(klon,klev) ! quantite de masse entrainee pour p-d  
       REAL pde_d(klon,klev) ! quantite de masse detrainee pour p-d  
 C  
       INTEGER  klab(klon,klev)  
       LOGICAL  llflag(klon)  
       INTEGER k, i, icall  
       REAL zzs  
 C----------------------------------------------------------------------  
 C SPECIFY LARGE SCALE PARAMETERS AT HALF LEVELS  
 C ADJUST TEMPERATURE FIELDS IF STATICLY UNSTABLE  
 C----------------------------------------------------------------------  
       DO 130 k = 2, klev  
 c  
       DO i = 1, klon  
          pgeoh(i,k)=pgeo(i,k)+(pgeo(i,k-1)-pgeo(i,k))*0.5  
          ptenh(i,k)=(MAX(RCPD*pten(i,k-1)+pgeo(i,k-1),  
      .             RCPD*pten(i,k)+pgeo(i,k))-pgeoh(i,k))/RCPD  
          pqsenh(i,k)=pqsen(i,k-1)  
          llflag(i)=.TRUE.  
       ENDDO  
 c  
       icall=0  
       CALL flxadjtq(paph(1,k),ptenh(1,k),pqsenh(1,k),llflag,icall)  
 c  
       DO i = 1, klon  
          pqenh(i,k)=MIN(pqen(i,k-1),pqsen(i,k-1))  
      .               +(pqsenh(i,k)-pqsen(i,k-1))  
          pqenh(i,k)=MAX(pqenh(i,k),0.)  
       ENDDO  
 c  
   130 CONTINUE  
 C  
       DO 140 i = 1, klon  
          ptenh(i,klev)=(RCPD*pten(i,klev)+pgeo(i,klev)-  
      1                   pgeoh(i,klev))/RCPD  
          pqenh(i,klev)=pqen(i,klev)  
          ptenh(i,1)=pten(i,1)  
          pqenh(i,1)=pqen(i,1)  
          pgeoh(i,1)=pgeo(i,1)  
   140 CONTINUE  
 c  
       DO 160 k = klev-1, 2, -1  
       DO 150 i = 1, klon  
          zzs = MAX(RCPD*ptenh(i,k)+pgeoh(i,k),  
      .             RCPD*ptenh(i,k+1)+pgeoh(i,k+1))  
          ptenh(i,k) = (zzs-pgeoh(i,k))/RCPD  
   150 CONTINUE  
   160 CONTINUE  
 C  
 C-----------------------------------------------------------------------  
 C INITIALIZE VALUES FOR UPDRAFTS AND DOWNDRAFTS  
 C-----------------------------------------------------------------------  
       DO k = 1, klev  
       DO i = 1, klon  
          ptu(i,k) = ptenh(i,k)  
          pqu(i,k) = pqenh(i,k)  
          plu(i,k) = 0.  
          pmfu(i,k) = 0.  
          pmfus(i,k) = 0.  
          pmfuq(i,k) = 0.  
          pdmfup(i,k) = 0.  
          pdpmel(i,k) = 0.  
          plude(i,k) = 0.  
 c  
          klab(i,k) = 0  
 c  
          ptd(i,k) = ptenh(i,k)  
          pqd(i,k) = pqenh(i,k)  
          pmfd(i,k) = 0.0  
          pmfds(i,k) = 0.0  
          pmfdq(i,k) = 0.0  
          pdmfdp(i,k) = 0.0  
 c  
          pen_u(i,k) = 0.0  
          pde_u(i,k) = 0.0  
          pen_d(i,k) = 0.0  
          pde_d(i,k) = 0.0  
       ENDDO  
       ENDDO  
 C  
       RETURN  
       END  
       SUBROUTINE flxbase(ptenh, pqenh, pgeoh, paph,  
      *     ptu, pqu, plu, ldcum, kcbot, klab)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
       IMPLICIT none  
 C----------------------------------------------------------------------  
 C THIS ROUTINE CALCULATES CLOUD BASE VALUES (T AND Q)  
 C  
 C INPUT ARE ENVIRONM. VALUES OF T,Q,P,PHI AT HALF LEVELS.  
 C IT RETURNS CLOUD BASE VALUES AND FLAGS AS FOLLOWS;  
 C   klab=1 FOR SUBCLOUD LEVELS  
 C   klab=2 FOR CONDENSATION LEVEL  
 C  
 C LIFT SURFACE AIR DRY-ADIABATICALLY TO CLOUD BASE  
 C (NON ENTRAINING PLUME,I.E.CONSTANT MASSFLUX)  
 C----------------------------------------------------------------------  
 C       ----------------------------------------------------------------  
       REAL ptenh(klon,klev), pqenh(klon,klev)  
       REAL pgeoh(klon,klev), paph(klon,klev+1)  
 C  
       REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)  
       INTEGER  klab(klon,klev), kcbot(klon)  
 C  
       LOGICAL llflag(klon), ldcum(klon)  
       INTEGER i, k, icall, is  
       REAL zbuo, zqold(klon)  
 C----------------------------------------------------------------------  
 C INITIALIZE VALUES AT LIFTING LEVEL  
 C----------------------------------------------------------------------  
       DO i = 1, klon  
          klab(i,klev)=1  
          kcbot(i)=klev-1  
          ldcum(i)=.FALSE.  
       ENDDO  
 C----------------------------------------------------------------------  
 C DO ASCENT IN SUBCLOUD LAYER,  
 C CHECK FOR EXISTENCE OF CONDENSATION LEVEL,  
 C ADJUST T,Q AND L ACCORDINGLY  
 C CHECK FOR BUOYANCY AND SET FLAGS  
 C----------------------------------------------------------------------  
       DO 290 k = klev-1, 2, -1  
 c  
       is = 0  
       DO i = 1, klon  
          IF (klab(i,k+1).EQ.1) is = is + 1  
          llflag(i) = .FALSE.  
          IF (klab(i,k+1).EQ.1) llflag(i) = .TRUE.  
       ENDDO  
       IF (is.EQ.0) GOTO 290  
 c  
       DO i = 1, klon  
       IF(llflag(i)) THEN  
          pqu(i,k) = pqu(i,k+1)  
          ptu(i,k) = ptu(i,k+1)+(pgeoh(i,k+1)-pgeoh(i,k))/RCPD  
          zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-  
      .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5  
          IF (zbuo.GT.0.) klab(i,k) = 1  
          zqold(i) = pqu(i,k)  
       ENDIF  
       ENDDO  
 c  
       icall=1  
       CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)  
 c  
       DO i = 1, klon  
       IF (llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN  
          klab(i,k) = 2  
          plu(i,k) = plu(i,k) + zqold(i)-pqu(i,k)  
          zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-  
      .          ptenh(i,k)*(1.+RETV*pqenh(i,k))+0.5  
          IF (zbuo.GT.0.) kcbot(i) = k  
          IF (zbuo.GT.0.) ldcum(i) = .TRUE.  
       ENDIF  
       ENDDO  
 c  
   290 CONTINUE  
 c  
       RETURN  
       END  
       SUBROUTINE flxasc(pdtime, ptenh, pqenh, pten, pqen, pqsen,  
      .     pgeo, pgeoh, pap, paph, pqte, pvervel,  
      .     ldland, ldcum, ktype, klab, ptu, pqu, plu,  
      .     pmfu, pmfub, pentr, pmfus, pmfuq,  
      .     pmful, plude, pdmfup, kcbot, kctop, kctop0, kcum,  
      .     pen_u, pde_u)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
             use yoecumf  
       IMPLICIT none  
 C----------------------------------------------------------------------  
 C THIS ROUTINE DOES THE CALCULATIONS FOR CLOUD ASCENTS  
 C FOR CUMULUS PARAMETERIZATION  
 C----------------------------------------------------------------------  
 C  
       REAL, intent(in):: pdtime  
       REAL pten(klon,klev), ptenh(klon,klev)  
       REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)  
       REAL pgeo(klon,klev), pgeoh(klon,klev)  
       REAL pap(klon,klev), paph(klon,klev+1)  
       REAL pqte(klon,klev)  
       REAL pvervel(klon,klev) ! vitesse verticale en Pa/s  
 C  
       REAL pmfub(klon), pentr(klon)  
       REAL ptu(klon,klev), pqu(klon,klev), plu(klon,klev)  
       REAL plude(klon,klev)  
       REAL pmfu(klon,klev), pmfus(klon,klev)  
       REAL pmfuq(klon,klev), pmful(klon,klev)  
       REAL pdmfup(klon,klev)  
       INTEGER ktype(klon), klab(klon,klev), kcbot(klon), kctop(klon)  
       INTEGER kctop0(klon)  
       LOGICAL ldland(klon), ldcum(klon)  
 C  
       REAL pen_u(klon,klev), pde_u(klon,klev)  
       REAL zqold(klon)  
       REAL zdland(klon)  
       LOGICAL llflag(klon)  
       INTEGER k, i, is, icall, kcum  
       REAL ztglace, zdphi, zqeen, zseen, zscde, zqude  
       REAL zmfusk, zmfuqk, zmfulk, zbuo, zdnoprc, zprcon, zlnew  
 c  
       REAL zpbot(klon), zptop(klon), zrho(klon)  
       REAL zdprho, zentr, zpmid, zmftest, zmfmax  
       LOGICAL llo1, llo2  
 c  
       REAL zwmax(klon), zzzmb  
       INTEGER klwmin(klon) ! level of maximum vertical velocity  
 C----------------------------------------------------------------------  
       ztglace = RTT - 13.  
 c  
 c Chercher le niveau ou la vitesse verticale est maximale:  
       DO i = 1, klon  
          klwmin(i) = klev  
          zwmax(i) = 0.0  
       ENDDO  
       DO k = klev, 3, -1  
       DO i = 1, klon  
       IF (pvervel(i,k).LT.zwmax(i)) THEN  
          zwmax(i) = pvervel(i,k)  
          klwmin(i) = k  
       ENDIF  
       ENDDO  
       ENDDO  
 C----------------------------------------------------------------------  
 C SET DEFAULT VALUES  
 C----------------------------------------------------------------------  
       DO i = 1, klon  
          IF (.NOT.ldcum(i)) ktype(i)=0  
       ENDDO  
 c  
       DO k=1,klev  
       DO i = 1, klon  
          plu(i,k)=0.  
          pmfu(i,k)=0.  
          pmfus(i,k)=0.  
          pmfuq(i,k)=0.  
          pmful(i,k)=0.  
          plude(i,k)=0.  
          pdmfup(i,k)=0.  
          IF(.NOT.ldcum(i).OR.ktype(i).EQ.3) klab(i,k)=0  
          IF(.NOT.ldcum(i).AND.paph(i,k).LT.4.E4) kctop0(i)=k  
       ENDDO  
       ENDDO  
 c  
       DO i = 1, klon  
       IF (ldland(i)) THEN  
          zdland(i)=3.0E4  
          zdphi=pgeoh(i,kctop0(i))-pgeoh(i,kcbot(i))  
          IF (ptu(i,kctop0(i)).GE.ztglace) zdland(i)=zdphi  
          zdland(i)=MAX(3.0E4,zdland(i))  
          zdland(i)=MIN(5.0E4,zdland(i))  
       ENDIF  
       ENDDO  
 C  
 C Initialiser les valeurs au niveau d'ascendance  
 C  
       DO i = 1, klon  
          kctop(i) = klev-1  
          IF (.NOT.ldcum(i)) THEN  
             kcbot(i) = klev-1  
             pmfub(i) = 0.  
             pqu(i,klev) = 0.  
          ENDIF  
          pmfu(i,klev) = pmfub(i)  
          pmfus(i,klev) = pmfub(i)*(RCPD*ptu(i,klev)+pgeoh(i,klev))  
          pmfuq(i,klev) = pmfub(i)*pqu(i,klev)  
       ENDDO  
 c  
       DO i = 1, klon  
          ldcum(i) = .FALSE.  
       ENDDO  
 C----------------------------------------------------------------------  
 C  DO ASCENT: SUBCLOUD LAYER (klab=1) ,CLOUDS (klab=2)  
 C  BY DOING FIRST DRY-ADIABATIC ASCENT AND THEN  
 C  BY ADJUSTING T,Q AND L ACCORDINGLY IN *flxadjtq*,  
 C  THEN CHECK FOR BUOYANCY AND SET FLAGS ACCORDINGLY  
 C----------------------------------------------------------------------  
       DO 480 k = klev-1,3,-1  
 c  
       IF (LMFMID .AND. k.LT.klev-1 .AND. k.GT.klev/2) THEN  
          DO i = 1, klon  
          IF (.NOT.ldcum(i) .AND. klab(i,k+1).EQ.0 .AND.  
      .       pqen(i,k).GT.0.9*pqsen(i,k)) THEN  
             ptu(i,k+1) = pten(i,k) +(pgeo(i,k)-pgeoh(i,k+1))/RCPD  
             pqu(i,k+1) = pqen(i,k)  
             plu(i,k+1) = 0.0  
             zzzmb = MAX(CMFCMIN, -pvervel(i,k)/RG)  
             zmfmax = (paph(i,k)-paph(i,k-1))/(RG*pdtime)  
             pmfub(i) = MIN(zzzmb,zmfmax)  
             pmfu(i,k+1) = pmfub(i)  
             pmfus(i,k+1) = pmfub(i)*(RCPD*ptu(i,k+1)+pgeoh(i,k+1))  
             pmfuq(i,k+1) = pmfub(i)*pqu(i,k+1)  
             pmful(i,k+1) = 0.0  
             pdmfup(i,k+1) = 0.0  
             kcbot(i) = k  
             klab(i,k+1) = 1  
             ktype(i) = 3  
             pentr(i) = ENTRMID  
          ENDIF  
          ENDDO  
       ENDIF  
 c  
       is = 0  
       DO i = 1, klon  
          is = is + klab(i,k+1)  
          IF (klab(i,k+1) .EQ. 0) klab(i,k) = 0  
          llflag(i) = .FALSE.  
          IF (klab(i,k+1) .GT. 0) llflag(i) = .TRUE.  
       ENDDO  
       IF (is .EQ. 0) GOTO 480  
 c  
 c calculer le taux d'entrainement et de detrainement  
 c  
       DO i = 1, klon  
          pen_u(i,k) = 0.0  
          pde_u(i,k) = 0.0  
          zrho(i)=paph(i,k+1)/(RD*ptenh(i,k+1))  
          zpbot(i)=paph(i,kcbot(i))  
          zptop(i)=paph(i,kctop0(i))  
       ENDDO  
 c  
       DO 125 i = 1, klon  
       IF(ldcum(i)) THEN  
          zdprho=(paph(i,k+1)-paph(i,k))/(RG*zrho(i))  
          zentr=pentr(i)*pmfu(i,k+1)*zdprho  
          llo1=k.LT.kcbot(i)  
          IF(llo1) pde_u(i,k)=zentr  
          zpmid=0.5*(zpbot(i)+zptop(i))  
          llo2=llo1.AND.ktype(i).EQ.2.AND.  
      .        (zpbot(i)-paph(i,k).LT.0.2E5.OR.  
      .         paph(i,k).GT.zpmid)  
          IF(llo2) pen_u(i,k)=zentr  
          llo2=llo1.AND.(ktype(i).EQ.1.OR.ktype(i).EQ.3).AND.  
      .        (k.GE.MAX(klwmin(i),kctop0(i)+2).OR.pap(i,k).GT.zpmid)  
          IF(llo2) pen_u(i,k)=zentr  
          llo1=pen_u(i,k).GT.0..AND.(ktype(i).EQ.1.OR.ktype(i).EQ.2)  
          IF(llo1) THEN  
             zentr=zentr*(1.+3.*(1.-MIN(1.,(zpbot(i)-pap(i,k))/1.5E4)))  
             pen_u(i,k)=pen_u(i,k)*(1.+3.*(1.-MIN(1.,  
      .                 (zpbot(i)-pap(i,k))/1.5E4)))  
             pde_u(i,k)=pde_u(i,k)*(1.+3.*(1.-MIN(1.,  
      .                 (zpbot(i)-pap(i,k))/1.5E4)))  
          ENDIF  
          IF(llo2.AND.pqenh(i,k+1).GT.1.E-5)  
      .   pen_u(i,k)=zentr+MAX(pqte(i,k),0.)/pqenh(i,k+1)*  
      .              zrho(i)*zdprho  
       ENDIF  
   125 CONTINUE  
 c  
 C----------------------------------------------------------------------  
 c DO ADIABATIC ASCENT FOR ENTRAINING/DETRAINING PLUME  
 C----------------------------------------------------------------------  
 c  
       DO 420 i = 1, klon  
       IF (llflag(i)) THEN  
          IF (k.LT.kcbot(i)) THEN  
             zmftest = pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)  
             zmfmax = MIN(zmftest,(paph(i,k)-paph(i,k-1))/(RG*pdtime))  
             pen_u(i,k)=MAX(pen_u(i,k)-MAX(0.0,zmftest-zmfmax),0.0)  
          ENDIF  
          pde_u(i,k)=MIN(pde_u(i,k),0.75*pmfu(i,k+1))  
 c calculer le flux de masse du niveau k a partir de celui du k+1  
          pmfu(i,k)=pmfu(i,k+1)+pen_u(i,k)-pde_u(i,k)  
 c calculer les valeurs Su, Qu et l du niveau k dans le panache montant  
          zqeen=pqenh(i,k+1)*pen_u(i,k)  
          zseen=(RCPD*ptenh(i,k+1)+pgeoh(i,k+1))*pen_u(i,k)  
          zscde=(RCPD*ptu(i,k+1)+pgeoh(i,k+1))*pde_u(i,k)  
          zqude=pqu(i,k+1)*pde_u(i,k)  
          plude(i,k)=plu(i,k+1)*pde_u(i,k)  
          zmfusk=pmfus(i,k+1)+zseen-zscde  
          zmfuqk=pmfuq(i,k+1)+zqeen-zqude  
          zmfulk=pmful(i,k+1)    -plude(i,k)  
          plu(i,k)=zmfulk*(1./MAX(CMFCMIN,pmfu(i,k)))  
          pqu(i,k)=zmfuqk*(1./MAX(CMFCMIN,pmfu(i,k)))  
          ptu(i,k)=(zmfusk*(1./MAX(CMFCMIN,pmfu(i,k)))-  
      1               pgeoh(i,k))/RCPD  
          ptu(i,k)=MAX(100.,ptu(i,k))  
          ptu(i,k)=MIN(400.,ptu(i,k))  
          zqold(i)=pqu(i,k)  
       ELSE  
          zqold(i)=0.0  
       ENDIF  
   420 CONTINUE  
 c  
 C----------------------------------------------------------------------  
 c DO CORRECTIONS FOR MOIST ASCENT BY ADJUSTING T,Q AND L  
 C----------------------------------------------------------------------  
 c  
       icall = 1  
       CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)  
 C  
       DO 440 i = 1, klon  
       IF(llflag(i).AND.pqu(i,k).NE.zqold(i)) THEN  
          klab(i,k) = 2  
          plu(i,k) = plu(i,k)+zqold(i)-pqu(i,k)  
          zbuo = ptu(i,k)*(1.+RETV*pqu(i,k))-  
      .          ptenh(i,k)*(1.+RETV*pqenh(i,k))  
          IF (klab(i,k+1).EQ.1) zbuo=zbuo+0.5  
          IF (zbuo.GT.0..AND.pmfu(i,k).GE.0.1*pmfub(i)) THEN  
             kctop(i) = k  
             ldcum(i) = .TRUE.  
             zdnoprc = 1.5E4  
             IF (ldland(i)) zdnoprc = zdland(i)  
             zprcon = CPRCON  
             IF ((zpbot(i)-paph(i,k)).LT.zdnoprc) zprcon = 0.0  
             zlnew=plu(i,k)/(1.+zprcon*(pgeoh(i,k)-pgeoh(i,k+1)))  
             pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))  
             plu(i,k)=zlnew  
          ELSE  
             klab(i,k)=0  
             pmfu(i,k)=0.  
          ENDIF  
       ENDIF  
   440 CONTINUE  
       DO 455 i = 1, klon  
       IF (llflag(i)) THEN  
          pmful(i,k)=plu(i,k)*pmfu(i,k)  
          pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)  
          pmfuq(i,k)=pqu(i,k)*pmfu(i,k)  
       ENDIF  
   455 CONTINUE  
 C  
   480 CONTINUE  
 C----------------------------------------------------------------------  
 C DETERMINE CONVECTIVE FLUXES ABOVE NON-BUOYANCY LEVEL  
 C    (NOTE: CLOUD VARIABLES LIKE T,Q AND L ARE NOT  
 C           AFFECTED BY DETRAINMENT AND ARE ALREADY KNOWN  
 C           FROM PREVIOUS CALCULATIONS ABOVE)  
 C----------------------------------------------------------------------  
       DO i = 1, klon  
          IF (kctop(i).EQ.klev-1) ldcum(i) = .FALSE.  
          kcbot(i) = MAX(kcbot(i),kctop(i))  
       ENDDO  
 c  
       ldcum(1)=ldcum(1)  
 c  
       is = 0  
       DO i = 1, klon  
          if (ldcum(i)) is = is + 1  
       ENDDO  
       kcum = is  
       IF (is.EQ.0) GOTO 800  
 c  
       DO 530 i = 1, klon  
       IF (ldcum(i)) THEN  
          k=kctop(i)-1  
          pde_u(i,k)=(1.-CMFCTOP)*pmfu(i,k+1)  
          plude(i,k)=pde_u(i,k)*plu(i,k+1)  
          pmfu(i,k)=pmfu(i,k+1)-pde_u(i,k)  
          zlnew=plu(i,k)  
          pdmfup(i,k)=MAX(0.,(plu(i,k)-zlnew)*pmfu(i,k))  
          plu(i,k)=zlnew  
          pmfus(i,k)=(RCPD*ptu(i,k)+pgeoh(i,k))*pmfu(i,k)  
          pmfuq(i,k)=pqu(i,k)*pmfu(i,k)  
          pmful(i,k)=plu(i,k)*pmfu(i,k)  
          plude(i,k-1)=pmful(i,k)  
       ENDIF  
   530 CONTINUE  
 C  
   800 CONTINUE  
       RETURN  
       END  
       SUBROUTINE flxflux(pdtime, pqen, pqsen, ptenh, pqenh, pap  
      .  ,  paph, ldland, pgeoh, kcbot, kctop, lddraf, kdtop  
      .  ,  ktype, ldcum, pmfu, pmfd, pmfus, pmfds  
      .  ,  pmfuq, pmfdq, pmful, plude, pdmfup, pdmfdp  
      .  ,  pten, prfl, psfl, pdpmel, ktopm2  
      .  ,  pmflxr, pmflxs)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
       use fcttre  
             use yoecumf  
       IMPLICIT none  
 C----------------------------------------------------------------------  
 C THIS ROUTINE DOES THE FINAL CALCULATION OF CONVECTIVE  
 C FLUXES IN THE CLOUD LAYER AND IN THE SUBCLOUD LAYER  
 C----------------------------------------------------------------------  
 C  
       REAL cevapcu(klev)  
 C     -----------------------------------------------------------------  
       REAL pqen(klon,klev), pqenh(klon,klev), pqsen(klon,klev)  
       REAL pten(klon,klev), ptenh(klon,klev)  
       REAL paph(klon,klev+1), pgeoh(klon,klev)  
 c  
       REAL pap(klon,klev)  
       REAL ztmsmlt, zdelta, zqsat  
 C  
       REAL pmfu(klon,klev), pmfus(klon,klev)  
       REAL pmfd(klon,klev), pmfds(klon,klev)  
       REAL pmfuq(klon,klev), pmful(klon,klev)  
       REAL pmfdq(klon,klev)  
       REAL plude(klon,klev)  
       REAL pdmfup(klon,klev), pdpmel(klon,klev)  
 cjq The variable maxpdmfdp(klon) has been introduced by Olivier Boucher  
 cjq 14/11/00 to fix the problem with the negative precipitation.        
       REAL pdmfdp(klon,klev), maxpdmfdp(klon,klev)  
       REAL prfl(klon), psfl(klon)  
       REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)  
       INTEGER  kcbot(klon), kctop(klon), ktype(klon)  
       LOGICAL  ldland(klon), ldcum(klon)  
       INTEGER k, kp, i  
       REAL zcons1, zcons2, zcucov, ztmelp2  
       REAL, intent(in):: pdtime  
       real zdp, zzp, zfac, zsnmlt, zrfl, zrnew  
       REAL zrmin, zrfln, zdrfl  
       REAL zpds, zpdr, zdenom  
       INTEGER ktopm2, itop, ikb  
 c  
       LOGICAL lddraf(klon)  
       INTEGER kdtop(klon)  
 c  
 c  
       DO 101 k=1,klev  
       CEVAPCU(k)=1.93E-6*261.*SQRT(1.E3/(38.3*0.293)  
      1 *SQRT(0.5*(paph(1,k)+paph(1,k+1))/paph(1,klev+1)) ) * 0.5/RG  
  101  CONTINUE  
 c  
 c SPECIFY CONSTANTS  
 c  
       zcons1 = RCPD/(RLMLT*RG*pdtime)  
       zcons2 = 1./(RG*pdtime)  
       zcucov = 0.05  
       ztmelp2 = RTT + 2.  
 c  
 c DETERMINE FINAL CONVECTIVE FLUXES  
 c  
       itop=klev  
       DO 110 i = 1, klon  
          itop=MIN(itop,kctop(i))  
          IF (.NOT.ldcum(i) .OR. kdtop(i).LT.kctop(i)) lddraf(i)=.FALSE.  
          IF(.NOT.ldcum(i)) ktype(i)=0  
   110 CONTINUE  
 c  
       ktopm2=itop-2  
       DO 120 k=ktopm2,klev  
       DO 115 i = 1, klon  
       IF(ldcum(i).AND.k.GE.kctop(i)-1) THEN  
          pmfus(i,k)=pmfus(i,k)-pmfu(i,k)*  
      .                (RCPD*ptenh(i,k)+pgeoh(i,k))  
          pmfuq(i,k)=pmfuq(i,k)-pmfu(i,k)*pqenh(i,k)  
          zdp = 1.5E4  
          IF ( ldland(i) ) zdp = 3.E4  
 c  
 c        l'eau liquide detrainee est precipitee quand certaines  
 c        conditions sont reunies (sinon, elle est consideree  
 c        evaporee dans l'environnement)  
 c  
          IF(paph(i,kcbot(i))-paph(i,kctop(i)).GE.zdp.AND.  
      .      pqen(i,k-1).GT.0.8*pqsen(i,k-1))  
      .      pdmfup(i,k-1)=pdmfup(i,k-1)+plude(i,k-1)  
 c  
          IF(lddraf(i).AND.k.GE.kdtop(i)) THEN  
             pmfds(i,k)=pmfds(i,k)-pmfd(i,k)*  
      .                   (RCPD*ptenh(i,k)+pgeoh(i,k))  
             pmfdq(i,k)=pmfdq(i,k)-pmfd(i,k)*pqenh(i,k)  
          ELSE  
             pmfd(i,k)=0.  
             pmfds(i,k)=0.  
             pmfdq(i,k)=0.  
             pdmfdp(i,k-1)=0.  
          END IF  
       ELSE  
          pmfu(i,k)=0.  
          pmfus(i,k)=0.  
          pmfuq(i,k)=0.  
          pmful(i,k)=0.  
          pdmfup(i,k-1)=0.  
          plude(i,k-1)=0.  
          pmfd(i,k)=0.  
          pmfds(i,k)=0.  
          pmfdq(i,k)=0.  
          pdmfdp(i,k-1)=0.  
       ENDIF  
   115 CONTINUE  
   120 CONTINUE  
 c  
       DO 130 k=ktopm2,klev  
       DO 125 i = 1, klon  
       IF(ldcum(i).AND.k.GT.kcbot(i)) THEN  
          ikb=kcbot(i)  
          zzp=((paph(i,klev+1)-paph(i,k))/  
      .        (paph(i,klev+1)-paph(i,ikb)))  
          IF (ktype(i).EQ.3) zzp = zzp**2  
          pmfu(i,k)=pmfu(i,ikb)*zzp  
          pmfus(i,k)=pmfus(i,ikb)*zzp  
          pmfuq(i,k)=pmfuq(i,ikb)*zzp  
          pmful(i,k)=pmful(i,ikb)*zzp  
       ENDIF  
   125 CONTINUE  
   130 CONTINUE  
 c  
 c CALCULATE RAIN/SNOW FALL RATES  
 c CALCULATE MELTING OF SNOW  
 c CALCULATE EVAPORATION OF PRECIP  
 c  
       DO k = 1, klev+1  
       DO i = 1, klon  
          pmflxr(i,k) = 0.0  
          pmflxs(i,k) = 0.0  
       ENDDO  
       ENDDO  
       DO k = ktopm2, klev  
       DO i = 1, klon  
       IF (ldcum(i)) THEN  
          IF (pmflxs(i,k).GT.0.0 .AND. pten(i,k).GT.ztmelp2) THEN  
             zfac=zcons1*(paph(i,k+1)-paph(i,k))  
             zsnmlt=MIN(pmflxs(i,k),zfac*(pten(i,k)-ztmelp2))  
             pdpmel(i,k)=zsnmlt  
             ztmsmlt=pten(i,k)-zsnmlt/zfac  
             zdelta=MAX(0.,SIGN(1.,RTT-ztmsmlt))  
             zqsat=R2ES*FOEEW(ztmsmlt, zdelta) / pap(i,k)  
             zqsat=MIN(0.5,zqsat)  
             zqsat=zqsat/(1.-RETV  *zqsat)  
             pqsen(i,k) = zqsat  
          ENDIF  
          IF (pten(i,k).GT.RTT) THEN  
          pmflxr(i,k+1)=pmflxr(i,k)+pdmfup(i,k)+pdmfdp(i,k)+pdpmel(i,k)  
          pmflxs(i,k+1)=pmflxs(i,k)-pdpmel(i,k)  
          ELSE  
            pmflxs(i,k+1)=pmflxs(i,k)+pdmfup(i,k)+pdmfdp(i,k)  
            pmflxr(i,k+1)=pmflxr(i,k)  
          ENDIF  
 c        si la precipitation est negative, on ajuste le plux du  
 c        panache descendant pour eliminer la negativite  
          IF ((pmflxr(i,k+1)+pmflxs(i,k+1)).LT.0.0) THEN  
             pdmfdp(i,k) = -pmflxr(i,k)-pmflxs(i,k)-pdmfup(i,k)  
             pmflxr(i,k+1) = 0.0  
             pmflxs(i,k+1) = 0.0  
             pdpmel(i,k) = 0.0  
          ENDIF  
       ENDIF  
       ENDDO  
       ENDDO  
 c  
 cjq The new variable is initialized here.  
 cjq It contains the humidity which is fed to the downdraft  
 cjq by evaporation of precipitation in the column below the base  
 cjq of convection.  
 cjq  
 cjq In the former version, this term has been subtracted from precip  
 cjq as well as the evaporation.  
 cjq        
       DO k = 1, klev  
       DO i = 1, klon  
          maxpdmfdp(i,k)=0.0  
       ENDDO  
       ENDDO  
       DO k = 1, klev  
        DO kp = k, klev  
         DO i = 1, klon  
          maxpdmfdp(i,k)=maxpdmfdp(i,k)+pdmfdp(i,kp)  
         ENDDO  
122         ENDDO         ENDDO
123        ENDDO      ENDDO
124  cjq End of initialization  
125  c            ! calculer la nature du sol (pour l'instant, ocean partout)
126        DO k = ktopm2, klev      DO i = 1, klon
127        DO i = 1, klon         land(i) = .FALSE.
128        IF (ldcum(i) .AND. k.GE.kcbot(i)) THEN      ENDDO
129           zrfl = pmflxr(i,k) + pmflxs(i,k)  
130           IF (zrfl.GT.1.0E-20) THEN      ! preparer les variables d'entree (attention: l'ordre des niveaux
131              zrnew=(MAX(0.,SQRT(zrfl/zcucov)-      ! verticaux augmente du haut vers le bas)
132       .            CEVAPCU(k)*(paph(i,k+1)-paph(i,k))*      DO k = 1, klev
133       .            MAX(0.,pqsen(i,k)-pqen(i,k))))**2*zcucov         DO i = 1, klon
134              zrmin=zrfl-zcucov*MAX(0.,0.8*pqsen(i,k)-pqen(i,k))            pt(i, k) = t(i, klev-k+1)
135       .            *zcons2*(paph(i,k+1)-paph(i,k))            pq(i, k) = q(i, klev-k+1)
136              zrnew=MAX(zrnew,zrmin)            paprsf(i, k) = pres_f(i, klev-k+1)
137              zrfln=MAX(zrnew,0.)            paprs(i, k) = pres_h(i, klev+1-k+1)
138              zdrfl=MIN(0.,zrfln-zrfl)            pvervel(i, k) = w(i, klev+1-k)
139  cjq At least the amount of precipiation needed to feed the downdraft            zcvgt(i, k) = con_t(i, klev-k+1)
140  cjq with humidity below the base of convection has to be left and can't            zcvgq(i, k) = con_q(i, klev-k+1)
141  cjq be evaporated (surely the evaporation can't be positive):              
142              zdrfl=MAX(zdrfl,            zdelta=MAX(0., SIGN(1., RTT-pt(i, k)))
143       .            MIN(-pmflxr(i,k)-pmflxs(i,k)-maxpdmfdp(i,k),0.0))            zqsat=R2ES*FOEEW ( pt(i, k), zdelta ) / paprsf(i, k)
144  cjq End of insertion            zqsat=MIN(0.5, zqsat)
145  c                        zqsat=zqsat/(1.-RETV  *zqsat)
146              zdenom=1.0/MAX(1.0E-20,pmflxr(i,k)+pmflxs(i,k))            pqs(i, k) = zqsat
147              IF (pten(i,k).GT.RTT) THEN         ENDDO
148                 zpdr = pdmfdp(i,k)      ENDDO
149                 zpds = 0.0      DO i = 1, klon
150              ELSE         paprs(i, klev+1) = pres_h(i, 1)
151                 zpdr = 0.0         zgeom(i, klev) = RD * pt(i, klev) &
152                 zpds = pdmfdp(i,k)              / (0.5*(paprs(i, klev+1)+paprsf(i, klev))) &
153              ENDIF              * (paprs(i, klev+1)-paprsf(i, klev))
154              pmflxr(i,k+1) = pmflxr(i,k) + zpdr + pdpmel(i,k)      ENDDO
155       .                    + zdrfl*pmflxr(i,k)*zdenom      DO k = klev-1, 1, -1
156              pmflxs(i,k+1) = pmflxs(i,k) + zpds - pdpmel(i,k)         DO i = 1, klon
157       .                    + zdrfl*pmflxs(i,k)*zdenom            zgeom(i, k) = zgeom(i, k+1) &
158              pdmfup(i,k) = pdmfup(i,k) + zdrfl                 + RD * 0.5*(pt(i, k+1)+pt(i, k)) / paprs(i, k+1) &
159           ELSE                 * (paprsf(i, k+1)-paprsf(i, k))
160              pmflxr(i,k+1) = 0.0         ENDDO
161              pmflxs(i,k+1) = 0.0      ENDDO
162              pdmfdp(i,k) = 0.0  
163              pdpmel(i,k) = 0.0      ! appeler la routine principale
164           ENDIF          
165           if (pmflxr(i,k) + pmflxs(i,k).lt.-1.e-26)      CALL flxmain(dtime, pt, pq, pqs, pqhfl, paprsf, paprs, zgeom, land, &
166       .    write(*,*) 'precip. < 1e-16 ',pmflxr(i,k) + pmflxs(i,k)           zcvgt, zcvgq, pvervel, rain, snow, kcbot, kctop, kdtop, zmfu, zmfd, &
167        ENDIF           zen_u, zde_u, zen_d, zde_d, d_t_bis, d_q_bis, zmflxr, zmflxs)
168        ENDDO  
169        ENDDO      ! De la même façon que l'on effectue le réindiçage pour la
170  c      ! température t et le champ q, on réindice les flux nécessaires à
171        DO 210 i = 1, klon      ! la convection des traceurs.
172           prfl(i) = pmflxr(i,klev+1)      DO k = 1, klev
173           psfl(i) = pmflxs(i,klev+1)         DO i = 1, klon
174    210 CONTINUE            d_q(i, klev+1-k) = dtime*d_q_bis(i, k)
175  c            d_t(i, klev+1-k) = dtime*d_t_bis(i, k)
176        RETURN         ENDDO
177        END      ENDDO
178        SUBROUTINE flxdtdq(ktopm2, paph, ldcum, pten  
179       .  ,  pmfus, pmfds, pmfuq, pmfdq, pmful, pdmfup, pdmfdp      DO i = 1, klon
180       .  ,  pdpmel, dt_con, dq_con)         pmfu(i, 1)= 0.
181        use dimens_m         pmfd(i, 1)= 0.
182        use dimphy         pen_d(i, 1)= 0.
183        use SUPHEC_M         pde_d(i, 1)= 0.
184        use yoethf_m      ENDDO
185              use yoecumf  
186        IMPLICIT none      DO k = 2, klev
187  c----------------------------------------------------------------------         DO i = 1, klon
188  c calculer les tendances T et Q            pmfu(i, klev+2-k)= zmfu(i, k)
189  c----------------------------------------------------------------------            pmfd(i, klev+2-k)= zmfd(i, k)
190  C     -----------------------------------------------------------------         ENDDO
191        LOGICAL  llo1      ENDDO
192  C  
193        REAL pten(klon,klev), paph(klon,klev+1)      DO k = 1, klev
194        REAL pmfus(klon,klev), pmfuq(klon,klev), pmful(klon,klev)         DO i = 1, klon
195        REAL pmfds(klon,klev), pmfdq(klon,klev)            pen_u(i, klev+1-k)=  zen_u(i, k)
196        REAL pdmfup(klon,klev)            pde_u(i, klev+1-k)=  zde_u(i, k)
197        REAL pdmfdp(klon,klev)         ENDDO
198        REAL pdpmel(klon,klev)      ENDDO
199        LOGICAL ldcum(klon)  
200        REAL dt_con(klon,klev), dq_con(klon,klev)      DO k = 1, klev-1
201  c         DO i = 1, klon
202        INTEGER ktopm2            pen_d(i, klev+1-k)= -zen_d(i, k+1)
203  c            pde_d(i, klev+1-k)= -zde_d(i, k+1)
204        INTEGER i, k         ENDDO
205        REAL zalv, zdtdt, zdqdt      ENDDO
206  c  
207        DO 210 k=ktopm2,klev-1      DO k = 1, klev+1
208        DO 220 i = 1, klon         DO i = 1, klon
209        IF (ldcum(i)) THEN            pmflxr(i, klev+2-k)= zmflxr(i, k)
210           llo1 = (pten(i,k)-RTT).GT.0.            pmflxs(i, klev+2-k)= zmflxs(i, k)
211           zalv = RLSTT         ENDDO
212           IF (llo1) zalv = RLVTT      ENDDO
213           zdtdt=RG/(paph(i,k+1)-paph(i,k))/RCPD  
214       .        *(pmfus(i,k+1)-pmfus(i,k)    END SUBROUTINE conflx
215       .         +pmfds(i,k+1)-pmfds(i,k)  
216       .          -RLMLT*pdpmel(i,k)  end module conflx_m
      .          -zalv*(pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))  
      .         )  
          dt_con(i,k)=zdtdt  
          zdqdt=RG/(paph(i,k+1)-paph(i,k))  
      .        *(pmfuq(i,k+1)-pmfuq(i,k)  
      .         +pmfdq(i,k+1)-pmfdq(i,k)  
      .          +pmful(i,k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k))  
          dq_con(i,k)=zdqdt  
       ENDIF  
   220 CONTINUE  
   210 CONTINUE  
 C  
       k = klev  
       DO 230 i = 1, klon  
       IF (ldcum(i)) THEN  
          llo1 = (pten(i,k)-RTT).GT.0.  
          zalv = RLSTT  
          IF (llo1) zalv = RLVTT  
          zdtdt=-RG/(paph(i,k+1)-paph(i,k))/RCPD  
      .         *(pmfus(i,k)+pmfds(i,k)+RLMLT*pdpmel(i,k)  
      .           -zalv*(pmful(i,k)+pdmfup(i,k)+pdmfdp(i,k)))  
          dt_con(i,k)=zdtdt  
          zdqdt=-RG/(paph(i,k+1)-paph(i,k))  
      .            *(pmfuq(i,k)+pmfdq(i,k)+pmful(i,k)  
      .             +pdmfup(i,k)+pdmfdp(i,k))  
          dq_con(i,k)=zdqdt  
       ENDIF  
   230 CONTINUE  
 C  
       RETURN  
       END  
       SUBROUTINE flxdlfs(ptenh, pqenh, pgeoh, paph, ptu, pqu,  
      .     ldcum, kcbot, kctop, pmfub, prfl, ptd, pqd,  
      .     pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
             use yoecumf  
       IMPLICIT none  
 C  
 C----------------------------------------------------------------------  
 C THIS ROUTINE CALCULATES LEVEL OF FREE SINKING FOR  
 C CUMULUS DOWNDRAFTS AND SPECIFIES T,Q,U AND V VALUES  
 C  
 C TO PRODUCE LFS-VALUES FOR CUMULUS DOWNDRAFTS  
 C FOR MASSFLUX CUMULUS PARAMETERIZATION  
 C  
 C INPUT ARE ENVIRONMENTAL VALUES OF T,Q,U,V,P,PHI  
 C AND UPDRAFT VALUES T,Q,U AND V AND ALSO  
 C CLOUD BASE MASSFLUX AND CU-PRECIPITATION RATE.  
 C IT RETURNS T,Q,U AND V VALUES AND MASSFLUX AT LFS.  
 C  
 C CHECK FOR NEGATIVE BUOYANCY OF AIR OF EQUAL PARTS OF  
 C MOIST ENVIRONMENTAL AIR AND CLOUD AIR.  
 C----------------------------------------------------------------------  
 C  
       REAL ptenh(klon,klev)  
       REAL pqenh(klon,klev)  
       REAL pgeoh(klon,klev), paph(klon,klev+1)  
       REAL ptu(klon,klev), pqu(klon,klev)  
       REAL pmfub(klon)  
       REAL prfl(klon)  
 C  
       REAL ptd(klon,klev), pqd(klon,klev)  
       REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)  
       REAL pdmfdp(klon,klev)  
       INTEGER  kcbot(klon), kctop(klon), kdtop(klon)  
       LOGICAL  ldcum(klon), lddraf(klon)  
 C  
       REAL ztenwb(klon,klev), zqenwb(klon,klev), zcond(klon)  
       REAL zttest, zqtest, zbuo, zmftop  
       LOGICAL  llo2(klon)  
       INTEGER i, k, is, icall  
 C----------------------------------------------------------------------  
       DO i= 1, klon  
          lddraf(i)=.FALSE.  
          kdtop(i)=klev+1  
       ENDDO  
 C  
 C----------------------------------------------------------------------  
 C DETERMINE LEVEL OF FREE SINKING BY  
 C DOING A SCAN FROM TOP TO BASE OF CUMULUS CLOUDS  
 C  
 C FOR EVERY POINT AND PROCEED AS FOLLOWS:  
 C     (1) DETEMINE WET BULB ENVIRONMENTAL T AND Q  
 C     (2) DO MIXING WITH CUMULUS CLOUD AIR  
 C     (3) CHECK FOR NEGATIVE BUOYANCY  
 C  
 C THE ASSUMPTION IS THAT AIR OF DOWNDRAFTS IS MIXTURE  
 C OF 50% CLOUD AIR + 50% ENVIRONMENTAL AIR AT WET BULB  
 C TEMPERATURE (I.E. WHICH BECAME SATURATED DUE TO  
 C EVAPORATION OF RAIN AND CLOUD WATER)  
 C----------------------------------------------------------------------  
 C  
       DO 290 k = 3, klev-3  
 C  
       is=0  
       DO 212 i= 1, klon  
          ztenwb(i,k)=ptenh(i,k)  
          zqenwb(i,k)=pqenh(i,k)  
          llo2(i) = ldcum(i).AND.prfl(i).GT.0.  
      .             .AND..NOT.lddraf(i)  
      .             .AND.(k.LT.kcbot(i).AND.k.GT.kctop(i))  
          IF ( llo2(i) ) is = is + 1  
   212 CONTINUE  
       IF(is.EQ.0) GO TO 290  
 C  
       icall=2  
       CALL flxadjtq(paph(1,k), ztenwb(1,k), zqenwb(1,k), llo2, icall)  
 C  
 C----------------------------------------------------------------------  
 C DO MIXING OF CUMULUS AND ENVIRONMENTAL AIR  
 C AND CHECK FOR NEGATIVE BUOYANCY.  
 C THEN SET VALUES FOR DOWNDRAFT AT LFS.  
 C----------------------------------------------------------------------  
       DO 222 i= 1, klon  
       IF (llo2(i)) THEN  
          zttest=0.5*(ptu(i,k)+ztenwb(i,k))  
          zqtest=0.5*(pqu(i,k)+zqenwb(i,k))  
          zbuo=zttest*(1.+RETV*zqtest)-  
      .        ptenh(i,k)*(1.+RETV  *pqenh(i,k))  
          zcond(i)=pqenh(i,k)-zqenwb(i,k)  
          zmftop=-CMFDEPS*pmfub(i)  
          IF (zbuo.LT.0..AND.prfl(i).GT.10.*zmftop*zcond(i)) THEN  
             kdtop(i)=k  
             lddraf(i)=.TRUE.  
             ptd(i,k)=zttest  
             pqd(i,k)=zqtest  
             pmfd(i,k)=zmftop  
             pmfds(i,k)=pmfd(i,k)*(RCPD*ptd(i,k)+pgeoh(i,k))  
             pmfdq(i,k)=pmfd(i,k)*pqd(i,k)  
             pdmfdp(i,k-1)=-0.5*pmfd(i,k)*zcond(i)  
             prfl(i)=prfl(i)+pdmfdp(i,k-1)  
          ENDIF  
       ENDIF  
   222 CONTINUE  
 c  
   290 CONTINUE  
 C  
       RETURN  
       END  
       SUBROUTINE flxddraf(ptenh, pqenh, pgeoh, paph, prfl,  
      .           ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp,  
      .           lddraf, pen_d, pde_d)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
             use yoecumf  
       IMPLICIT none  
 C  
 C----------------------------------------------------------------------  
 C          THIS ROUTINE CALCULATES CUMULUS DOWNDRAFT DESCENT  
 C  
 C          TO PRODUCE THE VERTICAL PROFILES FOR CUMULUS DOWNDRAFTS  
 C          (I.E. T,Q,U AND V AND FLUXES)  
 C  
 C          INPUT IS T,Q,P,PHI,U,V AT HALF LEVELS.  
 C          IT RETURNS FLUXES OF S,Q AND EVAPORATION RATE  
 C          AND U,V AT LEVELS WHERE DOWNDRAFT OCCURS  
 C  
 C          CALCULATE MOIST DESCENT FOR ENTRAINING/DETRAINING PLUME BY  
 C          A) MOVING AIR DRY-ADIABATICALLY TO NEXT LEVEL BELOW AND  
 C          B) CORRECTING FOR EVAPORATION TO OBTAIN SATURATED STATE.  
 C  
 C----------------------------------------------------------------------  
 C  
       REAL ptenh(klon,klev), pqenh(klon,klev)  
       REAL pgeoh(klon,klev), paph(klon,klev+1)  
 C  
       REAL ptd(klon,klev), pqd(klon,klev)  
       REAL pmfd(klon,klev), pmfds(klon,klev), pmfdq(klon,klev)  
       REAL pdmfdp(klon,klev)  
       REAL prfl(klon)  
       LOGICAL lddraf(klon)  
 C  
       REAL pen_d(klon,klev), pde_d(klon,klev), zcond(klon)  
       LOGICAL llo2(klon), llo1  
       INTEGER i, k, is, icall, itopde  
       REAL zentr, zseen, zqeen, zsdde, zqdde, zmfdsk, zmfdqk, zdmfdp  
       REAL zbuo  
 C----------------------------------------------------------------------  
 C CALCULATE MOIST DESCENT FOR CUMULUS DOWNDRAFT BY  
 C       (A) CALCULATING ENTRAINMENT RATES, ASSUMING  
 C           LINEAR DECREASE OF MASSFLUX IN PBL  
 C       (B) DOING MOIST DESCENT - EVAPORATIVE COOLING  
 C           AND MOISTENING IS CALCULATED IN *flxadjtq*  
 C       (C) CHECKING FOR NEGATIVE BUOYANCY AND  
 C           SPECIFYING FINAL T,Q,U,V AND DOWNWARD FLUXES  
 C  
       DO 180 k = 3, klev  
 c  
       is = 0  
       DO i = 1, klon  
          llo2(i)=lddraf(i).AND.pmfd(i,k-1).LT.0.  
          IF (llo2(i)) is = is + 1  
       ENDDO  
       IF (is.EQ.0) GOTO 180  
 c  
       DO i = 1, klon  
       IF (llo2(i)) THEN  
          zentr = ENTRDD*pmfd(i,k-1)*RD*ptenh(i,k-1)/  
      .           (RG*paph(i,k-1))*(paph(i,k)-paph(i,k-1))  
          pen_d(i,k) = zentr  
          pde_d(i,k) = zentr  
       ENDIF  
       ENDDO  
 c  
       itopde = klev-2  
       IF (k.GT.itopde) THEN  
          DO i = 1, klon  
          IF (llo2(i)) THEN  
             pen_d(i,k)=0.  
             pde_d(i,k)=pmfd(i,itopde)*  
      .      (paph(i,k)-paph(i,k-1))/(paph(i,klev+1)-paph(i,itopde))  
          ENDIF  
          ENDDO  
       ENDIF  
 C  
       DO i = 1, klon  
       IF (llo2(i)) THEN  
          pmfd(i,k) = pmfd(i,k-1)+pen_d(i,k)-pde_d(i,k)  
          zseen = (RCPD*ptenh(i,k-1)+pgeoh(i,k-1))*pen_d(i,k)  
          zqeen = pqenh(i,k-1)*pen_d(i,k)  
          zsdde = (RCPD*ptd(i,k-1)+pgeoh(i,k-1))*pde_d(i,k)  
          zqdde = pqd(i,k-1)*pde_d(i,k)  
          zmfdsk = pmfds(i,k-1)+zseen-zsdde  
          zmfdqk = pmfdq(i,k-1)+zqeen-zqdde  
          pqd(i,k) = zmfdqk*(1./MIN(-CMFCMIN,pmfd(i,k)))  
          ptd(i,k) = (zmfdsk*(1./MIN(-CMFCMIN,pmfd(i,k)))-  
      .               pgeoh(i,k))/RCPD  
          ptd(i,k) = MIN(400.,ptd(i,k))  
          ptd(i,k) = MAX(100.,ptd(i,k))  
          zcond(i) = pqd(i,k)  
       ENDIF  
       ENDDO  
 C  
       icall = 2  
       CALL flxadjtq(paph(1,k), ptd(1,k), pqd(1,k), llo2, icall)  
 C  
       DO i = 1, klon  
       IF (llo2(i)) THEN  
          zcond(i) = zcond(i)-pqd(i,k)  
          zbuo = ptd(i,k)*(1.+RETV  *pqd(i,k))-  
      .          ptenh(i,k)*(1.+RETV  *pqenh(i,k))  
          llo1 = zbuo.LT.0..AND.(prfl(i)-pmfd(i,k)*zcond(i).GT.0.)  
          IF (.not.llo1) pmfd(i,k) = 0.0  
          pmfds(i,k) = (RCPD*ptd(i,k)+pgeoh(i,k))*pmfd(i,k)  
          pmfdq(i,k) = pqd(i,k)*pmfd(i,k)  
          zdmfdp = -pmfd(i,k)*zcond(i)  
          pdmfdp(i,k-1) = zdmfdp  
          prfl(i) = prfl(i)+zdmfdp  
       ENDIF  
       ENDDO  
 c  
   180 CONTINUE  
       RETURN  
       END  
       SUBROUTINE flxadjtq(pp, pt, pq, ldflag, kcall)  
       use dimens_m  
       use dimphy  
       use SUPHEC_M  
       use yoethf_m  
       use fcttre  
       IMPLICIT none  
 c======================================================================  
 c Objet: ajustement entre T et Q  
 c======================================================================  
 C NOTE: INPUT PARAMETER kcall DEFINES CALCULATION AS  
 C        kcall=0    ENV. T AND QS IN*CUINI*  
 C        kcall=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)  
 C        kcall=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)  
 C  
 C  
       REAL pt(klon), pq(klon), pp(klon)  
       LOGICAL ldflag(klon)  
       INTEGER kcall  
 c  
       REAL zcond(klon), zcond1  
       REAL Z5alvcp, z5alscp, zalvdcp, zalsdcp  
       REAL zdelta, zcvm5, zldcp, zqsat, zcor  
       INTEGER is, i  
 C  
       z5alvcp = r5les*RLVTT/RCPD  
       z5alscp = r5ies*RLSTT/RCPD  
       zalvdcp = rlvtt/RCPD  
       zalsdcp = rlstt/RCPD  
 C  
   
       DO i = 1, klon  
          zcond(i) = 0.0  
       ENDDO  
   
       DO 210 i =1, klon  
       IF (ldflag(i)) THEN  
          zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))  
          zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp  
          zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp  
          zqsat = R2ES*FOEEW(pt(i),zdelta) / pp(i)  
          zqsat = MIN(0.5,zqsat)  
          zcor = 1./(1.-RETV*zqsat)  
          zqsat = zqsat*zcor  
          zcond(i) = (pq(i)-zqsat)  
      .     / (1. + FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor))  
          IF (kcall.EQ.1) zcond(i) = MAX(zcond(i),0.)  
          IF (kcall.EQ.2) zcond(i) = MIN(zcond(i),0.)  
          pt(i) = pt(i) + zldcp*zcond(i)  
          pq(i) = pq(i) - zcond(i)  
       ENDIF  
   210 CONTINUE  
 C  
       is = 0  
       DO i =1, klon  
          IF (zcond(i).NE.0.) is = is + 1  
       ENDDO  
       IF (is.EQ.0) GOTO 230  
 C  
       DO 220 i = 1, klon  
       IF(ldflag(i).AND.zcond(i).NE.0.) THEN  
          zdelta = MAX(0.,SIGN(1.,RTT-pt(i)))  
          zcvm5 = z5alvcp*(1.-zdelta) + zdelta*z5alscp  
          zldcp = zalvdcp*(1.-zdelta) + zdelta*zalsdcp  
          zqsat = R2ES* FOEEW(pt(i),zdelta) / pp(i)  
          zqsat = MIN(0.5,zqsat)  
          zcor = 1./(1.-RETV*zqsat)  
          zqsat = zqsat*zcor  
          zcond1 = (pq(i)-zqsat)  
      .     / (1. + FOEDE(pt(i),zdelta,zcvm5,zqsat,zcor))  
          pt(i) = pt(i) + zldcp*zcond1  
          pq(i) = pq(i) - zcond1  
       ENDIF  
   220 CONTINUE  
 C  
   230 CONTINUE  
       RETURN  
       END  
       SUBROUTINE flxsetup  
             use yoecumf  
       IMPLICIT none  
 C  
 C     THIS ROUTINE DEFINES DISPOSABLE PARAMETERS FOR MASSFLUX SCHEME  
 C  
 C  
       ENTRPEN=1.0E-4  ! ENTRAINMENT RATE FOR PENETRATIVE CONVECTION  
       ENTRSCV=3.0E-4  ! ENTRAINMENT RATE FOR SHALLOW CONVECTION  
       ENTRMID=1.0E-4  ! ENTRAINMENT RATE FOR MIDLEVEL CONVECTION  
       ENTRDD =2.0E-4  ! ENTRAINMENT RATE FOR DOWNDRAFTS  
       CMFCTOP=0.33  ! RELATIVE CLOUD MASSFLUX AT LEVEL ABOVE NONBUO LEVEL  
       CMFCMAX=1.0  ! MAXIMUM MASSFLUX VALUE ALLOWED FOR UPDRAFTS ETC  
       CMFCMIN=1.E-10  ! MINIMUM MASSFLUX VALUE (FOR SAFETY)  
       CMFDEPS=0.3  ! FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS  
       CPRCON =2.0E-4  ! CONVERSION FROM CLOUD WATER TO RAIN  
       RHCDD=1.  ! RELATIVE SATURATION IN DOWNDRAFRS (NO LONGER USED)  
 c                 (FORMULATION IMPLIES SATURATION)  
       LMFPEN = .TRUE.  
       LMFSCV = .TRUE.  
       LMFMID = .TRUE.  
       LMFDD = .TRUE.  
       LMFDUDV = .TRUE.  
 c  
       RETURN  
       END  

Legend:
Removed from v.51  
changed lines
  Added in v.62

  ViewVC Help
Powered by ViewVC 1.1.21