Index: LMDZ6/branches/Amaury_dev/libf/phylmd/conflx.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/conflx.F90	(revision 5112)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/conflx.F90	(revision 5142)
@@ -1,8 +1,7 @@
-
 ! $Header$
 
 SUBROUTINE conflx(dtime, pres_h, pres_f, t, q, con_t, con_q, pqhfl, w, d_t, &
-    d_q, rain, snow, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
-    kdtop, pmflxr, pmflxs)
+        d_q, rain, snow, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
+        kdtop, pmflxr, pmflxs)
 
   USE dimphy
@@ -20,5 +19,5 @@
   ! Entree:
   REAL dtime ! pas d'integration (s)
-  REAL pres_h(klon, klev+1) ! pression half-level (Pa)
+  REAL pres_h(klon, klev + 1) ! pression half-level (Pa)
   REAL pres_f(klon, klev) ! pression full-level (Pa)
   REAL t(klon, klev) ! temperature (K)
@@ -39,6 +38,6 @@
   REAL rain(klon) ! pluie (mm/s)
   REAL snow(klon) ! neige (mm/s)
-  REAL pmflxr(klon, klev+1)
-  REAL pmflxs(klon, klev+1)
+  REAL pmflxr(klon, klev + 1)
+  REAL pmflxs(klon, klev + 1)
   INTEGER kcbot(klon) ! niveau du bas de la convection
   INTEGER kctop(klon) ! niveau du haut de la convection
@@ -53,5 +52,5 @@
   REAL d_t_bis(klon, klev)
   REAL d_q_bis(klon, klev)
-  REAL paprs(klon, klev+1)
+  REAL paprs(klon, klev + 1)
   REAL paprsf(klon, klev)
   REAL zgeom(klon, klev)
@@ -65,8 +64,7 @@
   REAL zde_u(klon, klev)
   REAL zde_d(klon, klev)
-  REAL zmflxr(klon, klev+1)
-  REAL zmflxs(klon, klev+1)
+  REAL zmflxr(klon, klev + 1)
+  REAL zmflxs(klon, klev + 1)
   ! AA
-
 
   INTEGER i, k
@@ -117,28 +115,28 @@
   DO k = 1, klev
     DO i = 1, klon
-      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)
-
-      zdelta = max(0., sign(1.,rtt-pt(i,k)))
-      zqsat = r2es*foeew(pt(i,k), zdelta)/paprsf(i, k)
+      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)
+
+      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)
+      zqsat = zqsat / (1. - retv * zqsat)
       pqs(i, k) = zqsat
     END DO
   END DO
   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))
+    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))
   END DO
   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))
+      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))
     END DO
   END DO
@@ -147,6 +145,6 @@
 
   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)
+          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)
 
   ! AA--------------------------------------------------------
@@ -158,6 +156,6 @@
   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)
+      d_q(i, klev + 1 - k) = dtime * d_q_bis(i, k)
+      d_t(i, klev + 1 - k) = dtime * d_t_bis(i, k)
     END DO
   END DO
@@ -172,6 +170,6 @@
   DO k = 2, klev
     DO i = 1, klon
-      pmfu(i, klev+2-k) = zmfu(i, k)
-      pmfd(i, klev+2-k) = zmfd(i, k)
+      pmfu(i, klev + 2 - k) = zmfu(i, k)
+      pmfd(i, klev + 2 - k) = zmfd(i, k)
     END DO
   END DO
@@ -179,6 +177,6 @@
   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)
+      pen_u(i, klev + 1 - k) = zen_u(i, k)
+      pde_u(i, klev + 1 - k) = zde_u(i, k)
     END DO
   END DO
@@ -186,6 +184,6 @@
   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)
+      pen_d(i, klev + 1 - k) = -zen_d(i, k + 1)
+      pde_d(i, klev + 1 - k) = -zde_d(i, k + 1)
     END DO
   END DO
@@ -193,22 +191,22 @@
   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)
-    END DO
-  END DO
-
+      pmflxr(i, klev + 2 - k) = zmflxr(i, k)
+      pmflxs(i, klev + 2 - k) = zmflxs(i, k)
+    END DO
+  END DO
 
 END SUBROUTINE conflx
 ! --------------------------------------------------------------------
 SUBROUTINE flxmain(pdtime, pten, pqen, pqsen, pqhfl, pap, paph, pgeo, ldland, &
-    ptte, pqte, pvervel, prsfc, pssfc, kcbot, kctop, kdtop, & ! *
-                                                              ! ldcum, ktype,
-    pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, dt_con, dq_con, pmflxr, pmflxs)
+        ptte, pqte, pvervel, prsfc, pssfc, kcbot, kctop, kdtop, & ! *
+        ! ldcum, ktype,
+        pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, dt_con, dq_con, pmflxr, pmflxs)
   USE dimphy
+  USE lmdz_YOECUMF
+
   IMPLICIT NONE
   ! ------------------------------------------------------------------
   include "YOMCST.h"
   include "YOETHF.h"
-  include "YOECUMF.h"
   ! ----------------------------------------------------------------
   REAL pten(klon, klev), pqen(klon, klev), pqsen(klon, klev)
@@ -216,5 +214,5 @@
   REAL pqte(klon, klev)
   REAL pvervel(klon, klev)
-  REAL pgeo(klon, klev), pap(klon, klev), paph(klon, klev+1)
+  REAL pgeo(klon, klev), pap(klon, klev), paph(klon, klev + 1)
   REAL pqhfl(klon)
 
@@ -234,6 +232,6 @@
   REAL zdqpbl(klon), zdqcv(klon), zdhpbl(klon)
   REAL zrfl(klon)
-  REAL pmflxr(klon, klev+1)
-  REAL pmflxs(klon, klev+1)
+  REAL pmflxr(klon, klev + 1)
+  REAL pmflxs(klon, klev + 1)
   INTEGER ilab(klon, klev), ictop0(klon)
   LOGICAL llo1
@@ -275,6 +273,6 @@
   ! ----------------------------------------------------------------------
   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)
+          ptu, pqu, ptd, pqd, pmfd, zmfds, zmfdq, zdmfdp, pmfu, zmfus, zmfuq, &
+          zdmfup, zdpmel, plu, plude, ilab, pen_u, pde_u, pen_d, pde_d)
   ! ---------------------------------------------------------------------
   ! determiner les valeurs au niveau de base de la tour convective
@@ -290,5 +288,5 @@
   k = 1
   DO i = 1, klon
-    zdqcv(i) = pqte(i, k)*(paph(i,k+1)-paph(i,k))
+    zdqcv(i) = pqte(i, k) * (paph(i, k + 1) - paph(i, k))
     zdhpbl(i) = 0.0
     zdqpbl(i) = 0.0
@@ -297,9 +295,9 @@
   DO k = 2, klev
     DO i = 1, klon
-      zdqcv(i) = zdqcv(i) + pqte(i, k)*(paph(i,k+1)-paph(i,k))
+      zdqcv(i) = zdqcv(i) + pqte(i, k) * (paph(i, k + 1) - paph(i, k))
       IF (k>=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))
+        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))
       END IF
     END DO
@@ -308,5 +306,5 @@
   DO i = 1, klon
     ktype(i) = 2
-    IF (zdqcv(i)>max(0.,-1.5*pqhfl(i)*rg)) ktype(i) = 1
+    IF (zdqcv(i)>max(0., -1.5 * pqhfl(i) * rg)) ktype(i) = 1
     ! cc         if (zdqcv(i).GT.MAX(0.,-1.1*pqhfl(i)*RG)) ktype(i) = 1
   END DO
@@ -319,7 +317,7 @@
     ikb = kcbot(i)
     zqumqe = pqu(i, ikb) + plu(i, ikb) - zqenh(i, ikb)
-    zdqmin = max(0.01*zqenh(i,ikb), 1.E-10)
+    zdqmin = max(0.01 * zqenh(i, ikb), 1.E-10)
     IF (zdqpbl(i)>0. .AND. zqumqe>zdqmin .AND. ldcum(i)) THEN
-      zmfub(i) = zdqpbl(i)/(rg*max(zqumqe,zdqmin))
+      zmfub(i) = zdqpbl(i) / (rg * max(zqumqe, zdqmin))
     ELSE
       zmfub(i) = 0.01
@@ -327,9 +325,9 @@
     END IF
     IF (ktype(i)==2) THEN
-      zdh = rcpd*(ptu(i,ikb)-ztenh(i,ikb)) + rlvtt*zqumqe
-      zdh = rg*max(zdh, 1.0E5*zdqmin)
-      IF (zdhpbl(i)>0. .AND. ldcum(i)) zmfub(i) = zdhpbl(i)/zdh
+      zdh = rcpd * (ptu(i, ikb) - ztenh(i, ikb)) + rlvtt * zqumqe
+      zdh = rg * max(zdh, 1.0E5 * zdqmin)
+      IF (zdhpbl(i)>0. .AND. ldcum(i)) zmfub(i) = zdhpbl(i) / zdh
     END IF
-    zmfmax = (paph(i,ikb)-paph(i,ikb-1))/(rg*pdtime)
+    zmfmax = (paph(i, ikb) - paph(i, ikb - 1)) / (rg * pdtime)
     zmfub(i) = min(zmfub(i), zmfmax)
     zentr(i) = entrscv
@@ -345,17 +343,17 @@
   DO i = 1, klon
     ikb = kcbot(i)
-    zhcbase(i) = rcpd*ptu(i, ikb) + zgeoh(i, ikb) + rlvtt*pqu(i, ikb)
+    zhcbase(i) = rcpd * ptu(i, ikb) + zgeoh(i, ikb) + rlvtt * pqu(i, ikb)
     ictop0(i) = kcbot(i) - 1
   END DO
 
-  zalvdcp = rlvtt/rcpd
+  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.)
+      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<ictop0(i) .AND. zhcbase(i)>zhhat) ictop0(i) = k
     END DO
@@ -365,7 +363,7 @@
 
   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)
+          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==0) GO TO 1000
 
@@ -395,9 +393,9 @@
     ! 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)
+            zmfub, zrfl, ptd, pqd, pmfd, zmfds, zmfdq, zdmfdp, kdtop, lddraf)
 
     ! calculer le panache descendant
     CALL flxddraf(ztenh, zqenh, zgeoh, paph, zrfl, ptd, pqd, pmfd, zmfds, &
-      zmfdq, zdmfdp, lddraf, pen_d, pde_d)
+            zmfdq, zdmfdp, lddraf, pen_d, pde_d)
 
     ! calculer de nouveau le flux de masse entrant a travers la base
@@ -410,22 +408,22 @@
         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)
+        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)>0. .AND. zqumqe>zdqmin .AND. ldcum(i) .AND. &
-            zmfub(i)<zmfmax) THEN
-          zmfub1(i) = zdqpbl(i)/(rg*max(zqumqe,zdqmin))
+                zmfub(i)<zmfmax) THEN
+          zmfub1(i) = zdqpbl(i) / (rg * max(zqumqe, zdqmin))
         ELSE
           zmfub1(i) = zmfub(i)
         END IF
         IF (ktype(i)==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)>0. .AND. ldcum(i)) zmfub1(i) = zdhpbl(i)/zdh
+          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)>0. .AND. ldcum(i)) zmfub1(i) = zdhpbl(i) / zdh
         END IF
-        IF (.NOT. ((ktype(i)==1 .OR. ktype(i)==2) .AND. abs(zmfub1(i)-zmfub(i &
-          ))<0.2*zmfub(i))) zmfub1(i) = zmfub(i)
+        IF (.NOT. ((ktype(i)==1 .OR. ktype(i)==2) .AND. abs(zmfub1(i) - zmfub(i &
+                ))<0.2 * zmfub(i))) zmfub1(i) = zmfub(i)
       END IF
     END DO
@@ -433,11 +431,11 @@
       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
+          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
         END IF
       END DO
@@ -453,7 +451,7 @@
   ! -----------------------------------------------------------------------
   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)
+          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)
 
   ! -----------------------------------------------------------------------
@@ -462,7 +460,7 @@
   ! -----------------------------------------------------------------------
   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)
+          kcbot, kctop, lddraf, kdtop, ktype, ldcum, pmfu, pmfd, zmfus, zmfds, &
+          zmfuq, zmfdq, zmful, plude, zdmfup, zdmfdp, pten, prsfc, pssfc, zdpmel, &
+          itopm2, pmflxr, pmflxs)
 
   ! ----------------------------------------------------------------------
@@ -470,12 +468,12 @@
   ! ----------------------------------------------------------------------
   CALL flxdtdq(pdtime, itopm2, paph, ldcum, pten, zmfus, zmfds, zmfuq, zmfdq, &
-    zmful, zdmfup, zdmfdp, zdpmel, dt_con, dq_con)
-
-1000 CONTINUE
+          zmful, zdmfup, zdmfdp, zdpmel, dt_con, dq_con)
+
+  1000 CONTINUE
 
 END SUBROUTINE flxmain
 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)
+        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 dimphy
   IMPLICIT NONE
@@ -493,5 +491,5 @@
   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 paph(klon, klev + 1) ! pression aux demi-niveaux
   REAL ptenh(klon, klev) ! temperature aux demi-niveaux
   REAL pqenh(klon, klev) ! humidite aux demi-niveaux
@@ -532,24 +530,24 @@
 
     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)
+      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.
     END DO
 
     iCALL = 0
-    CALL flxadjtq(paph(1,k), ptenh(1,k), pqsenh(1,k), llflag, icall)
-
-    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.)
-    END DO
-
-  END DO
-
-  DO i = 1, klon
-    ptenh(i, klev) = (rcpd*pten(i,klev)+pgeo(i,klev)-pgeoh(i,klev))/rcpd
+    CALL flxadjtq(paph(1, k), ptenh(1, k), pqsenh(1, k), llflag, icall)
+
+    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.)
+    END DO
+
+  END DO
+
+  DO i = 1, klon
+    ptenh(i, klev) = (rcpd * pten(i, klev) + pgeo(i, klev) - pgeoh(i, klev)) / rcpd
     pqenh(i, klev) = pqen(i, klev)
     ptenh(i, 1) = pten(i, 1)
@@ -560,6 +558,6 @@
   DO k = klev - 1, 2, -1
     DO 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
+      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
     END DO
   END DO
@@ -596,8 +594,7 @@
   END DO
 
-
 END SUBROUTINE flxini
 SUBROUTINE flxbase(ptenh, pqenh, pgeoh, paph, ptu, pqu, plu, ldcum, kcbot, &
-    klab)
+        klab)
   USE dimphy
   IMPLICIT NONE
@@ -617,5 +614,5 @@
   ! ----------------------------------------------------------------
   REAL ptenh(klon, klev), pqenh(klon, klev)
-  REAL pgeoh(klon, klev), paph(klon, klev+1)
+  REAL pgeoh(klon, klev), paph(klon, klev + 1)
 
   REAL ptu(klon, klev), pqu(klon, klev), plu(klon, klev)
@@ -643,7 +640,7 @@
     is = 0
     DO i = 1, klon
-      IF (klab(i,k+1)==1) is = is + 1
+      IF (klab(i, k + 1)==1) is = is + 1
       llflag(i) = .FALSE.
-      IF (klab(i,k+1)==1) llflag(i) = .TRUE.
+      IF (klab(i, k + 1)==1) llflag(i) = .TRUE.
     END DO
     IF (is==0) GO TO 290
@@ -651,8 +648,8 @@
     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
+        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>0.) klab(i, k) = 1
         zqold(i) = pqu(i, k)
@@ -661,12 +658,12 @@
 
     iCALL = 1
-    CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
-
-    DO i = 1, klon
-      IF (llflag(i) .AND. pqu(i,k)/=zqold(i)) THEN
+    CALL flxadjtq(paph(1, k), ptu(1, k), pqu(1, k), llflag, icall)
+
+    DO i = 1, klon
+      IF (llflag(i) .AND. pqu(i, k)/=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
+        zbuo = ptu(i, k) * (1. + retv * pqu(i, k)) - ptenh(i, k) * (1. + retv * pqenh(i, k) &
+                ) + 0.5
         IF (zbuo>0.) kcbot(i) = k
         IF (zbuo>0.) ldcum(i) = .TRUE.
@@ -674,13 +671,14 @@
     END DO
 
-290 END DO
-
+  290 END DO
 
 END SUBROUTINE flxbase
 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)
+        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 dimphy
+  USE lmdz_YOECUMF
+
   IMPLICIT NONE
   ! ----------------------------------------------------------------------
@@ -690,5 +688,4 @@
   include "YOMCST.h"
   include "YOETHF.h"
-  include "YOECUMF.h"
 
   REAL pdtime
@@ -696,5 +693,5 @@
   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 pap(klon, klev), paph(klon, klev + 1)
   REAL pqte(klon, klev)
   REAL pvervel(klon, klev) ! vitesse verticale en Pa/s
@@ -735,5 +732,5 @@
   DO k = klev, 3, -1
     DO i = 1, klon
-      IF (pvervel(i,k)<zwmax(i)) THEN
+      IF (pvervel(i, k)<zwmax(i)) THEN
         zwmax(i) = pvervel(i, k)
         klwmin(i) = k
@@ -758,5 +755,5 @@
       pdmfup(i, k) = 0.
       IF (.NOT. ldcum(i) .OR. ktype(i)==3) klab(i, k) = 0
-      IF (.NOT. ldcum(i) .AND. paph(i,k)<4.E4) kctop0(i) = k
+      IF (.NOT. ldcum(i) .AND. paph(i, k)<4.E4) kctop0(i) = k
     END DO
   END DO
@@ -766,5 +763,5 @@
       zdland(i) = 3.0E4
       zdphi = pgeoh(i, kctop0(i)) - pgeoh(i, kcbot(i))
-      IF (ptu(i,kctop0(i))>=ztglace) zdland(i) = zdphi
+      IF (ptu(i, kctop0(i))>=ztglace) zdland(i) = zdphi
       zdland(i) = max(3.0E4, zdland(i))
       zdland(i) = min(5.0E4, zdland(i))
@@ -782,6 +779,6 @@
     END IF
     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)
+    pmfus(i, klev) = pmfub(i) * (rcpd * ptu(i, klev) + pgeoh(i, klev))
+    pmfuq(i, klev) = pmfub(i) * pqu(i, klev)
   END DO
 
@@ -797,21 +794,21 @@
   DO k = klev - 1, 3, -1
 
-    IF (lmfmid .AND. k<klev-1) THEN
+    IF (lmfmid .AND. k<klev - 1) THEN
       DO i = 1, klon
-        IF (.NOT. ldcum(i) .AND. klab(i,k+1)==0 .AND. &
-            pqen(i,k)>0.9*pqsen(i,k) .AND. pap(i,k)/paph(i,klev+1)>0.4) 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)
+        IF (.NOT. ldcum(i) .AND. klab(i, k + 1)==0 .AND. &
+                pqen(i, k)>0.9 * pqsen(i, k) .AND. pap(i, k) / paph(i, klev + 1)>0.4) 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
+          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
+          klab(i, k + 1) = 1
           ktype(i) = 3
           pentr(i) = entrmid
@@ -822,8 +819,8 @@
     is = 0
     DO i = 1, klon
-      is = is + klab(i, k+1)
-      IF (klab(i,k+1)==0) klab(i, k) = 0
+      is = is + klab(i, k + 1)
+      IF (klab(i, k + 1)==0) klab(i, k) = 0
       llflag(i) = .FALSE.
-      IF (klab(i,k+1)>0) llflag(i) = .TRUE.
+      IF (klab(i, k + 1)>0) llflag(i) = .TRUE.
     END DO
     IF (is==0) GO TO 480
@@ -834,5 +831,5 @@
       pen_u(i, k) = 0.0
       pde_u(i, k) = 0.0
-      zrho(i) = paph(i, k+1)/(rd*ptenh(i,k+1))
+      zrho(i) = paph(i, k + 1) / (rd * ptenh(i, k + 1))
       zpbot(i) = paph(i, kcbot(i))
       zptop(i) = paph(i, kctop0(i))
@@ -841,24 +838,24 @@
     DO 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
+        zdprho = (paph(i, k + 1) - paph(i, k)) / (rg * zrho(i))
+        zentr = pentr(i) * pmfu(i, k + 1) * zdprho
         llo1 = k < kcbot(i)
         IF (llo1) pde_u(i, k) = zentr
-        zpmid = 0.5*(zpbot(i)+zptop(i))
-        llo2 = llo1 .AND. ktype(i) == 2 .AND. (zpbot(i)-paph(i,k)<0.2E5 .OR. &
-          paph(i,k)>zpmid)
+        zpmid = 0.5 * (zpbot(i) + zptop(i))
+        llo2 = llo1 .AND. ktype(i) == 2 .AND. (zpbot(i) - paph(i, k)<0.2E5 .OR. &
+                paph(i, k)>zpmid)
         IF (llo2) pen_u(i, k) = zentr
         llo2 = llo1 .AND. (ktype(i)==1 .OR. ktype(i)==3) .AND. &
-          (k>=max(klwmin(i),kctop0(i)+2) .OR. pap(i,k)>zpmid)
+                (k>=max(klwmin(i), kctop0(i) + 2) .OR. pap(i, k)>zpmid)
         IF (llo2) pen_u(i, k) = zentr
         llo1 = pen_u(i, k) > 0. .AND. (ktype(i)==1 .OR. ktype(i)==2)
         IF (llo1) THEN
-          fact = 1. + 3.*(1.-min(1.,(zpbot(i)-pap(i,k))/1.5E4))
-          zentr = zentr*fact
-          pen_u(i, k) = pen_u(i, k)*fact
-          pde_u(i, k) = pde_u(i, k)*fact
+          fact = 1. + 3. * (1. - min(1., (zpbot(i) - pap(i, k)) / 1.5E4))
+          zentr = zentr * fact
+          pen_u(i, k) = pen_u(i, k) * fact
+          pde_u(i, k) = pde_u(i, k) * fact
         END IF
-        IF (llo2 .AND. pqenh(i,k+1)>1.E-5) pen_u(i, k) = zentr + &
-          max(pqte(i,k), 0.)/pqenh(i, k+1)*zrho(i)*zdprho
+        IF (llo2 .AND. pqenh(i, k + 1)>1.E-5) pen_u(i, k) = zentr + &
+                max(pqte(i, k), 0.) / pqenh(i, k + 1) * zrho(i) * zdprho
       END IF
     END DO
@@ -871,26 +868,26 @@
       IF (llflag(i)) THEN
         IF (k<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)
+          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)
         END IF
-        pde_u(i, k) = min(pde_u(i,k), 0.75*pmfu(i,k+1))
+        pde_u(i, k) = min(pde_u(i, k), 0.75 * pmfu(i, k + 1))
         ! 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)
+        pmfu(i, k) = pmfu(i, k + 1) + pen_u(i, k) - pde_u(i, k)
         ! 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)))-pgeoh(i,k))/rcpd
-        ptu(i, k) = max(100., ptu(i,k))
-        ptu(i, k) = min(400., ptu(i,k))
+        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))) - 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
@@ -904,14 +901,14 @@
 
     iCALL = 1
-    CALL flxadjtq(paph(1,k), ptu(1,k), pqu(1,k), llflag, icall)
-
-    DO i = 1, klon
-      IF (llflag(i) .AND. pqu(i,k)/=zqold(i)) THEN
+    CALL flxadjtq(paph(1, k), ptu(1, k), pqu(1, k), llflag, icall)
+
+    DO i = 1, klon
+      IF (llflag(i) .AND. pqu(i, k)/=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)==1) zbuo = zbuo + 0.5
-        IF (zbuo>0. .AND. pmfu(i,k)>=0.1*pmfub(i)) THEN
+        zbuo = ptu(i, k) * (1. + retv * pqu(i, k)) - ptenh(i, k) * (1. + retv * pqenh(i, k) &
+                )
+        IF (klab(i, k + 1)==1) zbuo = zbuo + 0.5
+        IF (zbuo>0. .AND. pmfu(i, k)>=0.1 * pmfub(i)) THEN
           kctop(i) = k
           ldcum(i) = .TRUE.
@@ -919,7 +916,7 @@
           IF (ldland(i)) zdnoprc = zdland(i)
           zprcon = cprcon
-          IF ((zpbot(i)-paph(i,k))<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))
+          IF ((zpbot(i) - paph(i, k))<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
@@ -931,11 +928,11 @@
     DO 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)
-      END IF
-    END DO
-
-480 END DO
+        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)
+      END IF
+    END DO
+
+  480 END DO
   ! ----------------------------------------------------------------------
   ! DETERMINE CONVECTIVE FLUXES ABOVE NON-BUOYANCY LEVEL
@@ -945,5 +942,5 @@
   ! ----------------------------------------------------------------------
   DO i = 1, klon
-    IF (kctop(i)==klev-1) ldcum(i) = .FALSE.
+    IF (kctop(i)==klev - 1) ldcum(i) = .FALSE.
     kcbot(i) = max(kcbot(i), kctop(i))
   END DO
@@ -961,26 +958,28 @@
     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)
+      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))
+      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)
+      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)
     END IF
   END DO
 
-800 CONTINUE
+  800 CONTINUE
 
 END SUBROUTINE flxasc
 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)
+        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 dimphy
   USE lmdz_print_control, ONLY: prt_level
+  USE lmdz_YOECUMF
+
   IMPLICIT NONE
   ! ----------------------------------------------------------------------
@@ -990,5 +989,4 @@
   include "YOMCST.h"
   include "YOETHF.h"
-  include "YOECUMF.h"
 
   REAL cevapcu(klon, klev)
@@ -996,5 +994,5 @@
   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)
+  REAL paph(klon, klev + 1), pgeoh(klon, klev)
 
   REAL pap(klon, klev)
@@ -1011,5 +1009,5 @@
   REAL pdmfdp(klon, klev), maxpdmfdp(klon, klev)
   REAL prfl(klon), psfl(klon)
-  REAL pmflxr(klon, klev+1), pmflxs(klon, klev+1)
+  REAL pmflxr(klon, klev + 1), pmflxs(klon, klev + 1)
   INTEGER kcbot(klon), kctop(klon), ktype(klon)
   LOGICAL ldland(klon), ldcum(klon)
@@ -1028,6 +1026,6 @@
   DO k = 1, klev
     DO i = 1, klon
-      cevapcu(i, k) = 1.93E-6*261.*sqrt(1.E3/(38.3*0.293)*sqrt(0.5*(paph(i,k) &
-        +paph(i,k+1))/paph(i,klev+1)))*0.5/rg
+      cevapcu(i, k) = 1.93E-6 * 261. * sqrt(1.E3 / (38.3 * 0.293) * sqrt(0.5 * (paph(i, k) &
+              + paph(i, k + 1)) / paph(i, klev + 1))) * 0.5 / rg
     END DO
   END DO
@@ -1035,6 +1033,6 @@
   ! SPECIFY CONSTANTS
 
-  zcons1 = rcpd/(rlmlt*rg*pdtime)
-  zcons2 = 1./(rg*pdtime)
+  zcons1 = rcpd / (rlmlt * rg * pdtime)
+  zcons2 = 1. / (rg * pdtime)
   zcucov = 0.05
   ztmelp2 = rtt + 2.
@@ -1052,7 +1050,7 @@
   DO k = ktopm2, klev
     DO i = 1, klon
-      IF (ldcum(i) .AND. k>=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)
+      IF (ldcum(i) .AND. k>=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
@@ -1062,15 +1060,15 @@
         ! evaporee dans l'environnement)
 
-        IF (paph(i,kcbot(i))-paph(i,kctop(i))>=zdp .AND. pqen(i,k-1)>0.8* &
-          pqsen(i,k-1)) pdmfup(i, k-1) = pdmfup(i, k-1) + plude(i, k-1)
+        IF (paph(i, kcbot(i)) - paph(i, kctop(i))>=zdp .AND. pqen(i, k - 1)>0.8 * &
+                pqsen(i, k - 1)) pdmfup(i, k - 1) = pdmfup(i, k - 1) + plude(i, k - 1)
 
         IF (lddraf(i) .AND. k>=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)
+          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.
+          pdmfdp(i, k - 1) = 0.
         END IF
       ELSE
@@ -1079,10 +1077,10 @@
         pmfuq(i, k) = 0.
         pmful(i, k) = 0.
-        pdmfup(i, k-1) = 0.
-        plude(i, k-1) = 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.
+        pdmfdp(i, k - 1) = 0.
       END IF
     END DO
@@ -1093,10 +1091,10 @@
       IF (ldcum(i) .AND. k>kcbot(i)) THEN
         ikb = kcbot(i)
-        zzp = ((paph(i,klev+1)-paph(i,k))/(paph(i,klev+1)-paph(i,ikb)))
+        zzp = ((paph(i, klev + 1) - paph(i, k)) / (paph(i, klev + 1) - paph(i, ikb)))
         IF (ktype(i)==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
+        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
       END IF
     END DO
@@ -1116,29 +1114,29 @@
     DO i = 1, klon
       IF (ldcum(i)) THEN
-        IF (pmflxs(i,k)>0.0 .AND. pten(i,k)>ztmelp2) THEN
-          zfac = zcons1*(paph(i,k+1)-paph(i,k))
-          zsnmlt = min(pmflxs(i,k), zfac*(pten(i,k)-ztmelp2))
+        IF (pmflxs(i, k)>0.0 .AND. pten(i, k)>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)
+          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)
+          zqsat = zqsat / (1. - retv * zqsat)
           pqsen(i, k) = zqsat
         END IF
-        IF (pten(i,k)>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)
+        IF (pten(i, k)>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)
+          pmflxs(i, k + 1) = pmflxs(i, k) + pdmfup(i, k) + pdmfdp(i, k)
+          pmflxr(i, k + 1) = pmflxr(i, k)
         END IF
         ! si la precipitation est negative, on ajuste le plux du
         ! panache descendant pour eliminer la negativite
-        IF ((pmflxr(i,k+1)+pmflxs(i,k+1))<0.0) THEN
+        IF ((pmflxr(i, k + 1) + pmflxs(i, k + 1))<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
+          pmflxr(i, k + 1) = 0.0
+          pmflxs(i, k + 1) = 0.0
           pdpmel(i, k) = 0.0
         END IF
@@ -1174,12 +1172,12 @@
         zrfl = pmflxr(i, k) + pmflxs(i, k)
         IF (zrfl>1.0E-20) THEN
-          zrnew = (max(0.,sqrt(zrfl/zcucov)-cevapcu(i, &
-            k)*(paph(i,k+1)-paph(i,k))*max(0.,pqsen(i,k)-pqen(i,k))))**2* &
-            zcucov
-          zrmin = zrfl - zcucov*max(0., 0.8*pqsen(i,k)-pqen(i,k))*zcons2*( &
-            paph(i,k+1)-paph(i,k))
+          zrnew = (max(0., sqrt(zrfl / zcucov) - cevapcu(i, &
+                  k) * (paph(i, k + 1) - paph(i, k)) * max(0., pqsen(i, k) - pqen(i, k))))**2 * &
+                  zcucov
+          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)
+          zdrfl = min(0., zrfln - zrfl)
           ! jq At least the amount of precipiation needed to feed the
           ! downdraft
@@ -1187,10 +1185,10 @@
           ! 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))
+          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)>rtt) THEN
+          zdenom = 1.0 / max(1.0E-20, pmflxr(i, k) + pmflxs(i, k))
+          IF (pten(i, k)>rtt) THEN
             zpdr = pdmfdp(i, k)
             zpds = 0.0
@@ -1199,31 +1197,32 @@
             zpds = pdmfdp(i, k)
           END IF
-          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
+          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
+          pmflxr(i, k + 1) = 0.0
+          pmflxs(i, k + 1) = 0.0
           pdmfdp(i, k) = 0.0
           pdpmel(i, k) = 0.0
         END IF
-        IF (pmflxr(i,k)+pmflxs(i,k)<-1.E-26 .AND. prt_level>=1) WRITE (*, *) &
-          'precip. < 1e-16 ', pmflxr(i, k) + pmflxs(i, k)
-      END IF
-    END DO
-  END DO
-
-  DO i = 1, klon
-    prfl(i) = pmflxr(i, klev+1)
-    psfl(i) = pmflxs(i, klev+1)
-  END DO
-
+        IF (pmflxr(i, k) + pmflxs(i, k)<-1.E-26 .AND. prt_level>=1) WRITE (*, *) &
+                'precip. < 1e-16 ', pmflxr(i, k) + pmflxs(i, k)
+      END IF
+    END DO
+  END DO
+
+  DO i = 1, klon
+    prfl(i) = pmflxr(i, klev + 1)
+    psfl(i) = pmflxs(i, klev + 1)
+  END DO
 
 END SUBROUTINE flxflux
 SUBROUTINE flxdtdq(pdtime, ktopm2, paph, ldcum, pten, pmfus, pmfds, pmfuq, &
-    pmfdq, pmful, pdmfup, pdmfdp, pdpmel, dt_con, dq_con)
+        pmfdq, pmful, pdmfup, pdmfdp, pdpmel, dt_con, dq_con)
   USE dimphy
+  USE lmdz_YOECUMF
+
   IMPLICIT NONE
   ! ----------------------------------------------------------------------
@@ -1232,9 +1231,8 @@
   include "YOMCST.h"
   include "YOETHF.h"
-  include "YOECUMF.h"
   ! -----------------------------------------------------------------
   LOGICAL llo1
 
-  REAL pten(klon, klev), paph(klon, klev+1)
+  REAL pten(klon, klev), paph(klon, klev + 1)
   REAL pmfus(klon, klev), pmfuq(klon, klev), pmful(klon, klev)
   REAL pmfds(klon, klev), pmfdq(klon, klev)
@@ -1254,13 +1252,13 @@
     DO i = 1, klon
       IF (ldcum(i)) THEN
-        llo1 = (pten(i,k)-rtt) > 0.
+        llo1 = (pten(i, k) - rtt) > 0.
         zalv = rlstt
         IF (llo1) zalv = rlvtt
-        zdtdt = rg/(paph(i,k+1)-paph(i,k))/rcpd*(pmfus(i,k+1)-pmfus(i,k)+ &
-          pmfds(i,k+1)-pmfds(i,k)-rlmlt*pdpmel(i,k)-zalv*(pmful(i, &
-          k+1)-pmful(i,k)-pdmfup(i,k)-pdmfdp(i,k)))
+        zdtdt = rg / (paph(i, k + 1) - paph(i, k)) / rcpd * (pmfus(i, k + 1) - pmfus(i, k) + &
+                pmfds(i, k + 1) - pmfds(i, k) - rlmlt * pdpmel(i, k) - 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))
+        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
       END IF
@@ -1271,21 +1269,22 @@
   DO i = 1, klon
     IF (ldcum(i)) THEN
-      llo1 = (pten(i,k)-rtt) > 0.
+      llo1 = (pten(i, k) - rtt) > 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)))
+      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))
+      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
     END IF
   END DO
 
-
 END SUBROUTINE flxdtdq
 SUBROUTINE flxdlfs(ptenh, pqenh, pgeoh, paph, ptu, pqu, ldcum, kcbot, kctop, &
-    pmfub, prfl, ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)
+        pmfub, prfl, ptd, pqd, pmfd, pmfds, pmfdq, pdmfdp, kdtop, lddraf)
   USE dimphy
+  USE lmdz_YOECUMF
+
   IMPLICIT NONE
 
@@ -1307,9 +1306,8 @@
   include "YOMCST.h"
   include "YOETHF.h"
-  include "YOECUMF.h"
 
   REAL ptenh(klon, klev)
   REAL pqenh(klon, klev)
-  REAL pgeoh(klon, klev), paph(klon, klev+1)
+  REAL pgeoh(klon, klev), paph(klon, klev + 1)
   REAL ptu(klon, klev), pqu(klon, klev)
   REAL pmfub(klon)
@@ -1354,5 +1352,5 @@
       zqenwb(i, k) = pqenh(i, k)
       llo2(i) = ldcum(i) .AND. prfl(i) > 0. .AND. .NOT. lddraf(i) .AND. &
-        (k<kcbot(i) .AND. k>kctop(i))
+              (k<kcbot(i) .AND. k>kctop(i))
       IF (llo2(i)) is = is + 1
     END DO
@@ -1360,5 +1358,5 @@
 
     iCALL = 2
-    CALL flxadjtq(paph(1,k), ztenwb(1,k), zqenwb(1,k), llo2, icall)
+    CALL flxadjtq(paph(1, k), ztenwb(1, k), zqenwb(1, k), llo2, icall)
 
     ! ----------------------------------------------------------------------
@@ -1369,10 +1367,10 @@
     DO 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))
+        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<0. .AND. prfl(i)>10.*zmftop*zcond(i)) THEN
+        zmftop = -cmfdeps * pmfub(i)
+        IF (zbuo<0. .AND. prfl(i)>10. * zmftop * zcond(i)) THEN
           kdtop(i) = k
           lddraf(i) = .TRUE.
@@ -1380,19 +1378,20 @@
           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)
+          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)
         END IF
       END IF
     END DO
 
-290 END DO
-
+  290 END DO
 
 END SUBROUTINE flxdlfs
 SUBROUTINE flxddraf(ptenh, pqenh, pgeoh, paph, prfl, ptd, pqd, pmfd, pmfds, &
-    pmfdq, pdmfdp, lddraf, pen_d, pde_d)
+        pmfdq, pdmfdp, lddraf, pen_d, pde_d)
   USE dimphy
+  USE lmdz_YOECUMF
+
   IMPLICIT NONE
 
@@ -1414,8 +1413,7 @@
   include "YOMCST.h"
   include "YOETHF.h"
-  include "YOECUMF.h"
 
   REAL ptenh(klon, klev), pqenh(klon, klev)
-  REAL pgeoh(klon, klev), paph(klon, klev+1)
+  REAL pgeoh(klon, klev), paph(klon, klev + 1)
 
   REAL ptd(klon, klev), pqd(klon, klev)
@@ -1443,5 +1441,5 @@
     is = 0
     DO i = 1, klon
-      llo2(i) = lddraf(i) .AND. pmfd(i, k-1) < 0.
+      llo2(i) = lddraf(i) .AND. pmfd(i, k - 1) < 0.
       IF (llo2(i)) is = is + 1
     END DO
@@ -1450,6 +1448,6 @@
     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))
+        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
@@ -1462,6 +1460,6 @@
         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))
+          pde_d(i, k) = pmfd(i, itopde) * (paph(i, k) - paph(i, k - 1)) / &
+                  (paph(i, klev + 1) - paph(i, itopde))
         END IF
       END DO
@@ -1470,15 +1468,15 @@
     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))
+        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)
       END IF
@@ -1486,22 +1484,22 @@
 
     iCALL = 2
-    CALL flxadjtq(paph(1,k), ptd(1,k), pqd(1,k), llo2, icall)
+    CALL flxadjtq(paph(1, k), ptd(1, k), pqd(1, k), llo2, icall)
 
     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 < 0. .AND. (prfl(i)-pmfd(i,k)*zcond(i)>0.)
+        zbuo = ptd(i, k) * (1. + retv * pqd(i, k)) - ptenh(i, k) * (1. + retv * pqenh(i, k) &
+                )
+        llo1 = zbuo < 0. .AND. (prfl(i) - pmfd(i, k) * zcond(i)>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
+        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
       END IF
     END DO
 
-180 END DO
+  180 END DO
 
 END SUBROUTINE flxddraf
@@ -1530,9 +1528,8 @@
   include "FCTTRE.h"
 
-  z5alvcp = r5les*rlvtt/rcpd
-  z5alscp = r5ies*rlstt/rcpd
-  zalvdcp = rlvtt/rcpd
-  zalsdcp = rlstt/rcpd
-
+  z5alvcp = r5les * rlvtt / rcpd
+  z5alscp = r5ies * rlstt / rcpd
+  zalvdcp = rlvtt / rcpd
+  zalsdcp = rlstt / rcpd
 
   DO i = 1, klon
@@ -1542,15 +1539,15 @@
   DO 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)
+      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))
+      zcor = 1. / (1. - retv * zqsat)
+      zqsat = zqsat * zcor
+      zcond(i) = (pq(i) - zqsat) / (1. + foede(pt(i), zdelta, zcvm5, zqsat, zcor))
       IF (kcall==1) zcond(i) = max(zcond(i), 0.)
       IF (kcall==2) zcond(i) = min(zcond(i), 0.)
-      pt(i) = pt(i) + zldcp*zcond(i)
+      pt(i) = pt(i) + zldcp * zcond(i)
       pq(i) = pq(i) - zcond(i)
     END IF
@@ -1565,26 +1562,26 @@
   DO i = 1, klon
     IF (ldflag(i) .AND. zcond(i)/=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)
+      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
+      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
     END IF
   END DO
 
-230 CONTINUE
+  230 CONTINUE
 
 END SUBROUTINE flxadjtq
 SUBROUTINE flxsetup
+  USE lmdz_YOECUMF
+
   IMPLICIT NONE
 
   ! THIS ROUTINE DEFINES DISPOSABLE PARAMETERS FOR MASSFLUX SCHEME
-
-  include "YOECUMF.h"
 
   entrpen = 1.0E-4 ! ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
@@ -1605,4 +1602,3 @@
   lmfdudv = .TRUE.
 
-
 END SUBROUTINE flxsetup
