Index: LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F90	(revision 1999)
+++ LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F90	(revision 2056)
@@ -1,13 +1,15 @@
-SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
-    u, v, tra, h, lv, qnk, unk, vnk, hp, tv, tvp, ep, clw, sig, ment, qent, &
-    hent, uent, vent, nent, sigij, elij, supmax, ments, qents, traent)
-  ! **************************************************************
-  ! *
-  ! CV3P_MIXING : compute mixed draught properties and,         *
-  ! within a scaling factor, mixed draught        *
-  ! mass fluxes.                                  *
-  ! written by  : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
-  ! modified by :                                               *
-  ! **************************************************************
+SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
+                       ph, t, rr, rs, u, v, tra, h, lv, qnk, &
+                       unk, vnk, hp, tv, tvp, ep, clw, sig, &
+                       Ment, Qent, hent, uent, vent, nent, &
+                       Sigij, elij, supmax, Ments, Qents, traent)
+! **************************************************************
+! *
+! CV3P_MIXING : compute mixed draught properties and,         *
+! within a scaling factor, mixed draught        *
+! mass fluxes.                                  *
+! written by  : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
+! modified by :                                               *
+! **************************************************************
 
   IMPLICIT NONE
@@ -17,112 +19,114 @@
   include "YOMCST2.h"
 
-  ! inputs:
-  INTEGER ncum, nd, na, ntra, nloc
-  INTEGER icb(nloc), inb(nloc), nk(nloc)
-  REAL sig(nloc, nd)
-  REAL qnk(nloc), unk(nloc), vnk(nloc)
-  REAL ph(nloc, nd+1)
-  REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
-  REAL u(nloc, nd), v(nloc, nd)
-  REAL tra(nloc, nd, ntra) ! input of convect3
-  REAL lv(nloc, na)
-  REAL h(nloc, na) !liquid water static energy of environment
-  REAL hp(nloc, na) !liquid water static energy of air shed from adiab. asc.
-  REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na)
-
-  ! outputs:
-  REAL ment(nloc, na, na), qent(nloc, na, na)
-  REAL uent(nloc, na, na), vent(nloc, na, na)
-  REAL sigij(nloc, na, na), elij(nloc, na, na)
-  REAL supmax(nloc, na) ! Highest mixing fraction of mixed updraughts
-    ! with the sign of (h-hp)
-  REAL traent(nloc, nd, nd, ntra)
-  REAL ments(nloc, nd, nd), qents(nloc, nd, nd)
-  REAL hent(nloc, nd, nd)
-  INTEGER nent(nloc, nd)
-
-  ! local variables:
+!inputs:
+  INTEGER, INTENT (IN)                               :: ncum, nd, na
+  INTEGER, INTENT (IN)                               :: ntra, nloc
+  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
+  REAL, DIMENSION (nloc), INTENT (IN)                :: qnk, unk, vnk
+  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
+  REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environMent
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: tv, tvp
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, clw
+
+!outputs:
+  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Ment, Qent
+  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: uent, vent
+  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Sigij, elij
+  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: supmax(nloc, na) ! Highest mixing fraction of mixed
+                                                                         ! updraughts with the sign of (h-hp)
+  REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
+  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: Ments, Qents
+  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: hent
+  INTEGER, DIMENSION (nloc, nd), INTENT (OUT)        :: nent
+
+!local variables:
   INTEGER i, j, k, il, im, jm
   INTEGER num1, num2
-  REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
-  REAL alt, delp, delm
-  REAL qmixmax(nloc), rmixmax(nloc), sqmrmax(nloc)
-  REAL qmixmin(nloc), rmixmin(nloc), sqmrmin(nloc)
-  REAL signhpmh(nloc)
-  REAL sx(nloc), scrit2
-  REAL smid(nloc), sjmin(nloc), sjmax(nloc)
-  REAL sbef(nloc), sup(nloc), smin(nloc)
-  REAL asij(nloc), smax(nloc), scrit(nloc)
-  REAL sij(nloc, nd, nd)
-  REAL csum(nloc, nd)
-  REAL awat
-  LOGICAL lwork(nloc)
+  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
+  REAL                               :: alt, delp, delm
+  REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
+  REAL, DIMENSION (nloc)             :: Qmixmin, Rmixmin, sqmrmin
+  REAL, DIMENSION (nloc)             :: signhpmh
+  REAL, DIMENSION (nloc)             :: Sx
+  REAL                               :: Scrit2
+  REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
+  REAL, DIMENSION (nloc)             :: Sbef, sup, smin
+  REAL, DIMENSION (nloc)             :: ASij, smax, Scrit
+  REAL, DIMENSION (nloc, nd, nd)     :: Sij
+  REAL, DIMENSION (nloc, nd)         :: csum
+  REAL                               :: awat
+  LOGICAL, DIMENSION (nloc)          :: lwork
 
   REAL amxupcrit, df, ff
   INTEGER nstep
 
-  ! --   Mixing probability distribution functions
-
-  REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f
-
-  qcoef1(f) = tanh(f/gammas)
-  qcoef2(f) = (tanh(f/gammas)+gammas*log(cosh((1.-f)/gammas)/cosh(f/gammas)))
-  qff(f) = max(min(f,1.), 0.)
-  qfff(f) = min(qff(f), scut)
-  qmix1(f) = (tanh((qff(f)-fmax)/gammas)+qcoef1max)/qcoef2max
-  rmix1(f) = (gammas*log(cosh((qff(f)-fmax)/gammas))+qff(f)*qcoef1max)/ &
-    qcoef2max
-  qmix2(f) = -log(1.-qfff(f))/scut
-  rmix2(f) = (qfff(f)+(1.-qff(f))*log(1.-qfff(f)))/scut
-  qmix(f) = qqa1*qmix1(f) + qqa2*qmix2(f)
-  rmix(f) = qqa1*rmix1(f) + qqa2*rmix2(f)
+! --   Mixing probability distribution functions
+
+  REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
+
+  Qcoef1(F) = tanh(F/gammas)
+  Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
+  QFF(F) = max(min(F,1.), 0.)
+  QFFf(F) = min(QFF(F), scut)
+  Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
+  Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+QFF(F)*Qcoef1max)/Qcoef2max
+  Qmix2(F) = -log(1.-QFFf(F))/scut
+  Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
+  Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
+  Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
 
   INTEGER, SAVE :: ifrst
   DATA ifrst/0/
-  !$OMP THREADPRIVATE(ifrst)
-
-
-  ! =====================================================================
-  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
-  ! =====================================================================
-
-  ! -- Initialize mixing PDF coefficients
+!$OMP THREADPRIVATE(ifrst)
+
+
+! =====================================================================
+! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
+! =====================================================================
+
+! -- Initialize mixing PDF coefficients
   IF (ifrst==0) THEN
     ifrst = 1
-    qcoef1max = qcoef1(fmax)
-    qcoef2max = qcoef2(fmax)
+    Qcoef1max = Qcoef1(Fmax)
+    Qcoef2max = Qcoef2(Fmax)
 
   END IF
 
 
-  ! ori        do 360 i=1,ncum*nlp
+! ori        do 360 i=1,ncum*nlp
   DO j = 1, nl
     DO i = 1, ncum
       nent(i, j) = 0
-      ! in convect3, m is computed in cv3_closure
-      ! ori          m(i,1)=0.0
-    END DO
-  END DO
-
-  ! ori      do 400 k=1,nlp
-  ! ori       do 390 j=1,nlp
+! in convect3, m is computed in cv3_closure
+! ori          m(i,1)=0.0
+    END DO
+  END DO
+
+! ori      do 400 k=1,nlp
+! ori       do 390 j=1,nlp
   DO j = 1, nl
     DO k = 1, nl
       DO i = 1, ncum
-        qent(i, k, j) = rr(i, j)
+        Qent(i, k, j) = rr(i, j)
         uent(i, k, j) = u(i, j)
         vent(i, k, j) = v(i, j)
         elij(i, k, j) = 0.0
         hent(i, k, j) = 0.0
-        ! AC!            ment(i,k,j)=0.0
-        ! AC!            sij(i,k,j)=0.0
-      END DO
-    END DO
-  END DO
-
-  ! AC!
-  ment(1:ncum, 1:nd, 1:nd) = 0.0
-  sij(1:ncum, 1:nd, 1:nd) = 0.0
-  ! AC!
+!AC!            Ment(i,k,j)=0.0
+!AC!            Sij(i,k,j)=0.0
+      END DO
+    END DO
+  END DO
+
+!AC!
+  Ment(1:ncum, 1:nd, 1:nd) = 0.0
+  Sij(1:ncum, 1:nd, 1:nd) = 0.0
+!AC!
 
   DO k = 1, ntra
@@ -136,9 +140,9 @@
   END DO
 
-  ! =====================================================================
-  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
-  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
-  ! --- FRACTION (sij)
-  ! =====================================================================
+! =====================================================================
+! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
+! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
+! --- FRACTION (Sij)
+! =====================================================================
 
   DO i = minorig + 1, nl
@@ -146,6 +150,6 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
-            1)) .AND. (j<=inb(il))) THEN
+        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
+                         .AND. (j<=inb(il))) THEN
 
           rti = qnk(il) - ep(il, i)*clw(il, i)
@@ -155,22 +159,21 @@
           dei = denom
           IF (abs(dei)<0.01) dei = 0.01
-          sij(il, i, j) = anum/dei
-          sij(il, i, i) = 1.0
-          altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
+          Sij(il, i, j) = anum/dei
+          Sij(il, i, i) = 1.0
+          altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
           altem = altem/bf2
           cwat = clw(il, j)*(1.-ep(il,j))
-          stemp = sij(il, i, j)
+          stemp = Sij(il, i, j)
           IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
             anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
             denom = denom + lv(il, j)*(rr(il,i)-rti)
             IF (abs(denom)<0.01) denom = 0.01
-            sij(il, i, j) = anum/denom
-            altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - &
-              rs(il, j)
+            Sij(il, i, j) = anum/denom
+            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
             altem = altem - (bf2-1.)*cwat
           END IF
-          IF (sij(il,i,j)>0.0) THEN
-            ! cc                 ment(il,i,j)=m(il,i)
-            ment(il, i, j) = 1.
+          IF (Sij(il,i,j)>0.0) THEN
+!!!                 Ment(il,i,j)=m(il,i)
+            Ment(il, i, j) = 1.
             elij(il, i, j) = altem
             elij(il, i, j) = amax1(0.0, elij(il,i,j))
@@ -178,6 +181,6 @@
           END IF
 
-          sij(il, i, j) = amax1(0.0, sij(il,i,j))
-          sij(il, i, j) = amin1(1.0, sij(il,i,j))
+          Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
+          Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
         END IF ! new
       END DO
@@ -185,22 +188,20 @@
 
 
-    ! ***   if no air can entrain at level i assume that updraft detrains
-    ! ***
-    ! ***   at that level and calculate detrained air flux and properties
-    ! ***
-
-
-    ! @      do 170 i=icb(il),inb(il)
+! ***   if no air can entrain at level i assume that updraft detrains  ***
+! ***   at that level and calculate detrained air flux and properties  ***
+
+
+! @      do 170 i=icb(il),inb(il)
 
     DO il = 1, ncum
       IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
-        ! @      if(nent(il,i).eq.0)then
-        ! cc      ment(il,i,i)=m(il,i)
-        ment(il, i, i) = 1.
-        qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
+! @      if(nent(il,i).eq.0)then
+!!!       Ment(il,i,i)=m(il,i)
+        Ment(il, i, i) = 1.
+        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
         uent(il, i, i) = unk(il)
         vent(il, i, i) = vnk(il)
         elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
-        sij(il, i, i) = 0.0
+        Sij(il, i, i) = 0.0
       END IF
     END DO
@@ -220,19 +221,19 @@
     DO i = minorig, nl
       DO il = 1, ncum
-        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
-            inb(il))) THEN
-          sigij(il, i, j) = sij(il, i, j)
-        END IF
-      END DO
-    END DO
-  END DO
-  ! @      enddo
-
-  ! @170   continue
-
-  ! =====================================================================
-  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
-  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
-  ! =====================================================================
+        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
+            (i>=icb(il)) .AND. (i<=inb(il))) THEN
+          Sigij(il, i, j) = Sij(il, i, j)
+        END IF
+      END DO
+    END DO
+  END DO
+! @      enddo
+
+! @170   continue
+
+! =====================================================================
+! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
+! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
+! =====================================================================
 
   CALL zilch(csum, nloc*nd)
@@ -242,7 +243,7 @@
   END DO
 
-  ! ---------------------------------------------------------------
-  DO i = minorig + 1, nl !Loop on origin level "i"
-    ! ---------------------------------------------------------------
+! ---------------------------------------------------------------
+  DO i = minorig + 1, nl      !Loop on origin level "i"
+! ---------------------------------------------------------------
 
     num1 = 0
@@ -253,12 +254,12 @@
 
 
-    ! jyg1    Find maximum of SIJ for J>I, if any.
-
-    sx(:) = 0.
+!JYG1    Find maximum of SIJ for J>I, if any.
+
+    Sx(:) = 0.
 
     DO il = 1, ncum
       IF (i>=icb(il) .AND. i<=inb(il)) THEN
         signhpmh(il) = sign(1., hp(il,i)-h(il,i))
-        sbef(il) = max(0., signhpmh(il))
+        Sbef(il) = max(0., signhpmh(il))
       END IF
     END DO
@@ -267,8 +268,8 @@
       DO il = 1, ncum
         IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
-          IF (sbef(il)<sij(il,i,j)) THEN
-            sx(il) = max(sij(il,i,j), sx(il))
-          END IF
-          sbef(il) = sij(il, i, j)
+          IF (Sbef(il)<Sij(il,i,j)) THEN
+            Sx(il) = max(Sij(il,i,j), Sx(il))
+          END IF
+          Sbef(il) = Sij(il, i, j)
         END IF
       END DO
@@ -279,172 +280,169 @@
       IF (i>=icb(il) .AND. i<=inb(il)) THEN
         lwork(il) = (nent(il,i)/=0)
-        qp = qnk(il) - ep(il, i)*clw(il, i)
-        anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
-          (cpv-cpd)*t(il, i)*(qp-rr(il,i))
-        denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
-          (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
+        rti = qnk(il) - ep(il, i)*clw(il, i)
+        anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
+               (cpv-cpd)*t(il, i)*(rti-rr(il,i))
+        denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
+                (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
         IF (abs(denom)<0.01) denom = 0.01
-        scrit(il) = min(anum/denom, 1.)
-        alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
-
-        ! jyg1    Find new critical value Scrit2
-        ! such that : Sij > Scrit2  => mixed draught will detrain at J<I
-        ! Sij < Scrit2  => mixed draught will detrain at J>I
-
-        scrit2 = min(scrit(il), sx(il))*max(0., -signhpmh(il)) + &
-          scrit(il)*max(0., signhpmh(il))
-
-        scrit(il) = scrit2
-
-        ! jyg    Correction pour la nouvelle logique; la correction pour ALT
-        ! est un peu au hazard
-        IF (scrit(il)<=0.0) scrit(il) = 0.0
-        IF (alt<=0.0) scrit(il) = 1.0
+        Scrit(il) = min(anum/denom, 1.)
+        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
+
+!JYG1    Find new critical value Scrit2
+!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
+!                     Sij < Scrit2  => mixed draught will detrain at J>I
+
+        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
+                 Scrit(il)*max(0., signhpmh(il))
+
+        Scrit(il) = Scrit2
+
+!JYG    Correction pour la nouvelle logique; la correction pour ALT
+! est un peu au hazard
+        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
+        IF (alt<=0.0) Scrit(il) = 1.0
 
         smax(il) = 0.0
-        asij(il) = 0.0
-        sup(il) = 0. ! upper S-value reached by descending draughts
-      END IF
-    END DO
-
-    ! ---------------------------------------------------------------
-    DO j = minorig, nl !Loop on destination level "j"
-      ! ---------------------------------------------------------------
+        ASij(il) = 0.0
+        sup(il) = 0.      ! upper S-value reached by descending draughts
+      END IF
+    END DO
+
+! ---------------------------------------------------------------
+    DO j = minorig, nl         !Loop on destination level "j"
+! ---------------------------------------------------------------
 
       num2 = 0
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-          il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+            lwork(il)) num2 = num2 + 1
       END DO
       IF (num2<=0) GO TO 175
 
-      ! -----------------------------------------------
+! -----------------------------------------------
       IF (j>i) THEN
-        ! -----------------------------------------------
+! -----------------------------------------------
         DO il = 1, ncum
-          IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-              il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
-            IF (sij(il,i,j)>0.0) THEN
-              smid(il) = min(sij(il,i,j), scrit(il))
-              sjmax(il) = smid(il)
-              sjmin(il) = smid(il)
-              IF (smid(il)<smin(il) .AND. sij(il,i,j+1)<smid(il)) THEN
-                smin(il) = smid(il)
-                sjmax(il) = min((sij(il,i,j+1)+sij(il,i, &
-                  j))/2., sij(il,i,j), scrit(il))
-                sjmin(il) = max((sbef(il)+sij(il,i,j))/2., sij(il,i,j))
-                sjmin(il) = min(sjmin(il), scrit(il))
-                sbef(il) = sij(il, i, j)
+          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+              lwork(il)) THEN
+            IF (Sij(il,i,j)>0.0) THEN
+              Smid(il) = min(Sij(il,i,j), Scrit(il))
+              Sjmax(il) = Smid(il)
+              Sjmin(il) = Smid(il)
+              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
+                smin(il) = Smid(il)
+                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
+                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
+                Sjmin(il) = min(Sjmin(il), Scrit(il))
+                Sbef(il) = Sij(il, i, j)
               END IF
             END IF
           END IF
         END DO
-        ! -----------------------------------------------
+! -----------------------------------------------
       ELSE IF (j==i) THEN
-        ! -----------------------------------------------
+! -----------------------------------------------
         DO il = 1, ncum
-          IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-              il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
-            IF (sij(il,i,j)>0.0) THEN
-              smid(il) = 1.
-              sjmin(il) = max((sij(il,i,j-1)+smid(il))/2., scrit(il))*max(0., &
-                -signhpmh(il)) + min((sij(il,i,j+1)+smid(il))/2., scrit(il))* &
-                max(0., signhpmh(il))
-              sjmin(il) = max(sjmin(il), sup(il))
-              sjmax(il) = 1.
-
-              ! -           preparation des variables Scrit, Smin et Sbef
-              ! pour la partie j>i
-              scrit(il) = min(sjmin(il), sjmax(il), scrit(il))
+          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+              lwork(il)) THEN
+            IF (Sij(il,i,j)>0.0) THEN
+              Smid(il) = 1.
+              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
+                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
+              Sjmin(il) = max(Sjmin(il), sup(il))
+              Sjmax(il) = 1.
+
+! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
+              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
 
               smin(il) = 1.
-              sbef(il) = max(0., signhpmh(il))
-              supmax(il, i) = sign(scrit(il), -signhpmh(il))
-            END IF
-          END IF
-        END DO
-        ! -----------------------------------------------
+              Sbef(il) = max(0., signhpmh(il))
+              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
+            END IF
+          END IF
+        END DO
+! -----------------------------------------------
       ELSE IF (j<i) THEN
-        ! -----------------------------------------------
+! -----------------------------------------------
         DO il = 1, ncum
-          IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-              il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
-            IF (sij(il,i,j)>0.0) THEN
-              smid(il) = max(sij(il,i,j), scrit(il))
-              sjmax(il) = smid(il)
-              sjmin(il) = smid(il)
-              IF (smid(il)>smax(il) .AND. sij(il,i,j+1)>smid(il)) THEN
-                smax(il) = smid(il)
-                sjmax(il) = max((sij(il,i,j+1)+sij(il,i,j))/2., sij(il,i,j))
-                sjmax(il) = max(sjmax(il), scrit(il))
-                sjmin(il) = min((sbef(il)+sij(il,i,j))/2., sij(il,i,j))
-                sjmin(il) = max(sjmin(il), scrit(il))
-                sbef(il) = sij(il, i, j)
+          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+              lwork(il)) THEN
+            IF (Sij(il,i,j)>0.0) THEN
+              Smid(il) = max(Sij(il,i,j), Scrit(il))
+              Sjmax(il) = Smid(il)
+              Sjmin(il) = Smid(il)
+              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
+                smax(il) = Smid(il)
+                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
+                Sjmax(il) = max(Sjmax(il), Scrit(il))
+                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
+                Sjmin(il) = max(Sjmin(il), Scrit(il))
+                Sbef(il) = Sij(il, i, j)
               END IF
-              IF (abs(sjmin(il)-sjmax(il))>1.E-10) sup(il) = max(sjmin(il), &
-                sjmax(il), sup(il))
-            END IF
-          END IF
-        END DO
-        ! -----------------------------------------------
-      END IF
-      ! -----------------------------------------------
-
-
-      DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
-          IF (sij(il,i,j)>0.0) THEN
+              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
+                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
+            END IF
+          END IF
+        END DO
+! -----------------------------------------------
+      END IF
+! -----------------------------------------------
+
+
+      DO il = 1, ncum
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+            lwork(il)) THEN
+          IF (Sij(il,i,j)>0.0) THEN
             rti = qnk(il) - ep(il, i)*clw(il, i)
-            qmixmax(il) = qmix(sjmax(il))
-            qmixmin(il) = qmix(sjmin(il))
-            rmixmax(il) = rmix(sjmax(il))
-            rmixmin(il) = rmix(sjmin(il))
-            sqmrmax(il) = sjmax(il)*qmix(sjmax(il)) - rmix(sjmax(il))
-            sqmrmin(il) = sjmin(il)*qmix(sjmin(il)) - rmix(sjmin(il))
-
-            ment(il, i, j) = abs(qmixmax(il)-qmixmin(il))*ment(il, i, j)
-
-            ! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
-            IF (abs(qmixmax(il)-qmixmin(il))>1.E-10) THEN
-              sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/ &
-                (qmixmax(il)-qmixmin(il))
+            Qmixmax(il) = Qmix(Sjmax(il))
+            Qmixmin(il) = Qmix(Sjmin(il))
+            Rmixmax(il) = Rmix(Sjmax(il))
+            Rmixmin(il) = Rmix(Sjmin(il))
+            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
+            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
+
+            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
+
+! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
+            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
+              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
             ELSE
-              sigij(il, i, j) = 0.
-            END IF
-
-            ! --    Compute Qent, uent, vent according to the true mixing
-            ! fraction
-            qent(il, i, j) = (1.-sigij(il,i,j))*rti + &
-              sigij(il, i, j)*rr(il, i)
-            uent(il, i, j) = (1.-sigij(il,i,j))*unk(il) + &
-              sigij(il, i, j)*u(il, i)
-            vent(il, i, j) = (1.-sigij(il,i,j))*vnk(il) + &
-              sigij(il, i, j)*v(il, i)
-
-            ! --     Compute liquid water static energy of mixed draughts
-            ! IF (j .GT. i) THEN
-            ! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
-            ! awat=amax1(awat,0.0)
-            ! ELSE
-            ! awat = 0.
-            ! ENDIF
-            ! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
-            ! :         + Sigij(il,i,j)*H(il,i)
-            ! :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
-            ! IM 301008 beg
-            hent(il, i, j) = (1.-sigij(il,i,j))*hp(il, i) + &
-              sigij(il, i, j)*h(il, i)
-
-            elij(il, i, j) = qent(il, i, j) - rs(il, j)
-            elij(il, i, j) = elij(il, i, j) + ((h(il,j)-hent(il,i, &
-              j))*rs(il,j)*lv(il,j)/((cpd*(1.-qent(il,i,j))+ &
-              qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
-            elij(il, i, j) = elij(il, i, j)/(1.+lv(il,j)*lv(il,j)*rs(il,j)/(( &
-              cpd*(1.-qent(il,i,j))+qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
+              Sigij(il, i, j) = 0.
+            END IF
+
+! --    Compute Qent, uent, vent according to the true mixing fraction
+            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
+            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
+            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)
+
+! --     Compute liquid water static energy of mixed draughts
+!    IF (j .GT. i) THEN
+!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
+!      awat=amax1(awat,0.0)
+!    ELSE
+!      awat = 0.
+!    ENDIF
+!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
+!    :         + Sigij(il,i,j)*H(il,i)
+!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
+!IM 301008 beg
+            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
+
+            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
+            elij(il, i, j) = elij(il, i, j) + &
+                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
+                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
+            elij(il, i, j) = elij(il, i, j) / &
+                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
+                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
 
             elij(il, i, j) = max(elij(il,i,j), 0.)
 
-            elij(il, i, j) = min(elij(il,i,j), qent(il,i,j))
+            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
 
             IF (j>i) THEN
@@ -455,19 +453,16 @@
             END IF
 
-            ! print
-            ! *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
-            ! :         t(il,j))
-
-            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))* &
-              awat
-            ! IM 301008 end
-
-            ! print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ',
-            ! :               i,j,hent(il,i,j),sigij(il,i,j)
-
-            ! --      ASij is the integral of P(F) over the relevant F
-            ! interval
-            asij(il) = asij(il) + abs(qmixmax(il)*(1.-sjmax(il))+rmixmax(il)- &
-              qmixmin(il)*(1.-sjmin(il))-rmixmin(il))
+! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
+! :         t(il,j))
+
+            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
+!IM 301008 end
+
+! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
+! :               i,j,hent(il,i,j),Sigij(il,i,j)
+
+! --      ASij is the integral of P(F) over the relevant F interval
+            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
+                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
 
           END IF
@@ -476,33 +471,33 @@
       DO k = 1, ntra
         DO il = 1, ncum
-          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
-              1)) .AND. (j<=inb(il)) .AND. lwork(il)) THEN
-            IF (sij(il,i,j)>0.0) THEN
-              traent(il, i, j, k) = sigij(il, i, j)*tra(il, i, k) + &
-                (1.-sigij(il,i,j))*tra(il, nk(il), k)
-            END IF
-          END IF
-        END DO
-      END DO
-
-      ! --    If I=J (detrainement and entrainement at the same level), then
-      ! only the
-      ! --    adiabatic ascent part of the mixture is considered
+          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
+              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
+              lwork(il)) THEN
+            IF (Sij(il,i,j)>0.0) THEN
+              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
+                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
+            END IF
+          END IF
+        END DO
+      END DO
+
+! --    If I=J (detrainement and entrainement at the same level), then only the
+! --    adiabatic ascent part of the mixture is considered
       IF (i==j) THEN
         DO il = 1, ncum
-          IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-              il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
-            IF (sij(il,i,j)>0.0) THEN
+          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+              lwork(il)) THEN
+            IF (Sij(il,i,j)>0.0) THEN
               rti = qnk(il) - ep(il, i)*clw(il, i)
-              ! cc          Ment(il,i,i) =
-              ! m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
-              ment(il, i, i) = abs(qmixmax(il)*(1.-sjmax( &
-                il))+rmixmax(il)-qmixmin(il)*(1.-sjmin(il))-rmixmin(il))
-              qent(il, i, i) = rti
+!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
+              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
+                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
+              Qent(il, i, i) = rti
               uent(il, i, i) = unk(il)
               vent(il, i, i) = vnk(il)
               hent(il, i, i) = hp(il, i)
               elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
-              sigij(il, i, i) = 0.
+              Sigij(il, i, i) = 0.
             END IF
           END IF
@@ -510,7 +505,8 @@
         DO k = 1, ntra
           DO il = 1, ncum
-            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
-                1)) .AND. (j<=inb(il)) .AND. lwork(il)) THEN
-              IF (sij(il,i,j)>0.0) THEN
+            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
+                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
+                lwork(il)) THEN
+              IF (Sij(il,i,j)>0.0) THEN
                 traent(il, i, i, k) = tra(il, nk(il), k)
               END IF
@@ -521,10 +517,12 @@
       END IF
 
-175 END DO
+! ---------------------------------------------------------------
+175 END DO        ! End loop on destination level "j"
+! ---------------------------------------------------------------
 
     DO il = 1, ncum
       IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
-        asij(il) = amax1(1.0E-16, asij(il))
-        asij(il) = 1.0/asij(il)
+        ASij(il) = amax1(1.0E-16, ASij(il))
+        ASij(il) = 1.0/ASij(il)
         csum(il, i) = 0.0
       END IF
@@ -533,7 +531,7 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il)) THEN
-          ment(il, i, j) = ment(il, i, j)*asij(il)
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
+          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
         END IF
       END DO
@@ -542,7 +540,7 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il)) THEN
-          csum(il, i) = csum(il, i) + ment(il, i, j)
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
+          csum(il, i) = csum(il, i) + Ment(il, i, j)
         END IF
       END DO
@@ -550,15 +548,14 @@
 
     DO il = 1, ncum
-      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) &
-          THEN
-        ! cc     :     .and. csum(il,i).lt.m(il,i) ) then
+      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
+! cc     :     .and. csum(il,i).lt.m(il,i) ) then
         nent(il, i) = 0
-        ! cc        ment(il,i,i)=m(il,i)
-        ment(il, i, i) = 1.
-        qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
+! cc        Ment(il,i,i)=m(il,i)
+        Ment(il, i, i) = 1.
+        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
         uent(il, i, i) = unk(il)
         vent(il, i, i) = vnk(il)
         elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
-        sij(il, i, i) = 0.0
+        Sij(il, i, i) = 0.0
       END IF
     END DO ! il
@@ -566,7 +563,6 @@
     DO j = 1, ntra
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) &
-            THEN
-          ! cc     :     .and. csum(il,i).lt.m(il,i) ) then
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
+! cc     :     .and. csum(il,i).lt.m(il,i) ) then
           traent(il, i, i, j) = tra(il, nk(il), j)
         END IF
@@ -574,5 +570,8 @@
     END DO
 
-789 END DO
+! ---------------------------------------------------------------
+789 END DO              ! End loop on origin level "i"
+! ---------------------------------------------------------------
+
 
   RETURN
