Index: LMDZ5/trunk/libf/phylmd/cv3_vertmix.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cv3_vertmix.F90	(revision 2040)
+++ LMDZ5/trunk/libf/phylmd/cv3_vertmix.F90	(revision 2041)
@@ -1,4 +1,5 @@
-SUBROUTINE cv3_vertmix(len, nd, iflag, plim1, plim2, p, ph, t, q, u, v, w, &
-    wi, nk, tmix, thmix, qmix, qsmix, umix, vmix, plcl)
+SUBROUTINE cv3_vertmix(len, nd, iflag, plim1, plim2, p, ph, &
+                       t, q, u, v, w, &
+                       wi, nk, tmix, thmix, qmix, qsmix, umix, vmix, plcl)
   ! **************************************************************
   ! *
@@ -13,6 +14,9 @@
   ! ==============================================================
 
-  ! vertmix : determine theta et r du melange obtenu en brassant
-  ! adiabatiquement entre plim1 et plim2, avec une ponderation w.
+  ! vertmix : determines theta, t, q, qs, u and v of the mixture generated by
+  ! adiabatic mixing of air between plim1 and plim2 with weighting w.
+  ! If plim1 and plim2 fall within the same model layer, then theta, ... v
+  ! are those of that layer.
+  ! A minimum value (dpmin) is imposed upon plim1-plim2
 
   ! ===============================================================
@@ -22,32 +26,47 @@
   include "YOMCST.h"
   include "FCTTRE.h"
-  ! input :
-  INTEGER nd, len
-  INTEGER nk(len), iflag(len)
-  REAL t(len, nd), q(len, nd), w(nd)
-  REAL u(len, nd), v(len, nd)
-  REAL p(len, nd), ph(len, nd+1)
-  REAL plim1(len), plim2(len)
-  ! output :
-  REAL tmix(len), thmix(len), qmix(len), wi(len, nd)
-  REAL umix(len), vmix(len)
-  REAL qsmix(len)
-  REAL plcl(len)
-  ! internal variables :
-  INTEGER j1(len), j2(len), niflag7
-  REAL a, b
-  REAL ahm(len), dpw(len), coef(len)
-  REAL p1(len, nd), p2(len, nd)
-  REAL rdcp(len), a2(len), b2(len), pnk(len)
-  REAL rh(len), chi(len)
-  REAL cpn
-  REAL x, y, p0, p0m1, zdelta, zcor
-
+!inputs:
+  INTEGER, INTENT (IN)                      ::  nd, len
+  INTEGER, DIMENSION (len), INTENT (IN)     ::  nk
+  REAL, DIMENSION (nd), INTENT (IN)        ::  w
+  REAL, DIMENSION (len), INTENT (IN)        :: plim1, plim2
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: t, q
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: u, v
+  REAL, DIMENSION (len,nd), INTENT (IN)     :: p
+  REAL, DIMENSION (len,nd+1), INTENT (IN)   :: ph
+!input/output:
+  INTEGER, DIMENSION (len), INTENT (INOUT)  ::  iflag
+!outputs:
+  REAL, DIMENSION (len), INTENT (OUT)       :: tmix, thmix, qmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: umix, vmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: qsmix
+  REAL, DIMENSION (len), INTENT (OUT)       :: plcl
+  REAL, DIMENSION (len,nd), INTENT (OUT)    :: wi
+!internal variables :
   INTEGER i, j
-
+  INTEGER niflag7
+  INTEGER, DIMENSION(len)                   :: j1, j2
+  REAL                                      :: a, b
+  REAL                                      :: cpn
+  REAL                                      :: x, y, p0, p0m1, zdelta, zcor
+  REAL                                      :: dpmin=1.
+!$OMP THREADPRIVATE(dpmin)
+  REAL, DIMENSION(len)                      :: plim2p  ! = min(plim2(:),plim1(:)-dpmin)
+  REAL, DIMENSION(len)                      :: ahm, dpw, coef
+  REAL, DIMENSION(len)                      :: rdcp, a2, b2, pnk
+  REAL, DIMENSION(len)                      :: rh, chi
+  REAL, DIMENSION(len)                      :: eqwght
+  REAL, DIMENSION(len,nd)                   :: p1, p2
+
+
+!!  print *,' ->cv3_vertmix, plim1,plim2 ', plim1,plim2   !jyg
+  plim2p(:) = min(plim2(:),plim1(:)-dpmin)
+  j1(:)=nd
+  j2(:) = 0
   DO j = 1, nd
     DO i = 1, len
       IF (plim1(i)<=ph(i,j)) j1(i) = j
-      IF (plim2(i)>=ph(i,j+1) .AND. plim2(i)<ph(i,j)) j2(i) = j
+!!!      IF (plim2p(i)>=ph(i,j+1) .AND. plim2p(i)<ph(i,j)) j2(i) = j
+      IF (plim2p(i)< ph(i,j)) j2(i) = j
     END DO
   END DO
@@ -68,4 +87,5 @@
     pnk(i) = p(i, nk(i))
   END DO
+  eqwght(:) = 0.
 
   p0 = 1000.
@@ -73,6 +93,13 @@
 
   DO i = 1, len
-    coef(i) = 1./(plim1(i)-plim2(i))
-  END DO
+    IF (j2(i) < j1(i)) THEN
+      coef(i) = 1.
+      eqwght(i) = 1.
+    ELSE
+      coef(i) = 1./(plim1(i)-plim2p(i))
+    ENDIF
+  END DO
+
+!!  print *,'cv3_vertmix, j1,j2,coef ', j1,j2,coef  !jyg
 
   DO j = 1, nd
@@ -80,16 +107,20 @@
       IF (j>=j1(i) .AND. j<=j2(i)) THEN
         p1(i, j) = min(ph(i,j), plim1(i))
-        p2(i, j) = max(ph(i,j+1), plim2(i))
+        p2(i, j) = max(ph(i,j+1), plim2p(i))
         ! CRtest:couplage thermiques: deja normalise
         ! wi(i,j) = w(j)
         ! print*,'wi',wi(i,j)
-        wi(i, j) = w(j)*(p1(i,j)-p2(i,j))*coef(i)
+        wi(i, j) = w(j)*(p1(i,j)-p2(i,j))*coef(i)+eqwght(i)
         dpw(i) = dpw(i) + wi(i, j)
+
+!!  print *,'cv3_vertmix, j, wi(1,j),dpw ', j, wi(1,j),dpw  !jyg
+
       END IF
     END DO
   END DO
+
   ! CR:print
   ! do i=1,len
-  ! print*,'plim',plim1(i),plim2(i)
+  ! print*,'plim',plim1(i),plim2p(i)
   ! enddo
   DO j = 1, nd
@@ -108,4 +139,7 @@
     rdcp(i) = (rrd*(1.-qmix(i))+qmix(i)*rrv)/(cpd*(1.-qmix(i))+qmix(i)*cpv)
   END DO
+
+
+!!  print *,'cv3_vertmix, rdcp ', rdcp  !jyg
 
 
@@ -159,11 +193,11 @@
       rh(i) = max(rh(i), 0.)
       plcl(i) = pnk(i)*(rh(i)**chi(i))
-      IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag &
-        (i) = 8
+      IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) &
+          iflag(i) = 8
 
     ELSE
 
       niflag7 = niflag7 + 1
-      plcl(i) = plim2(i)
+      plcl(i) = plim2p(i)
 
     END IF ! iflag=7
@@ -172,4 +206,7 @@
 
   END DO
+
+!!  print *,' cv3_vertmix->'  !jyg
+
 
   RETURN