Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 874)
+++ trunk/LMDZ.GENERIC/README	(revision 875)
@@ -896,2 +896,9 @@
 - Cosmetic changes in optcv (mostly spaces and line breaks) 
     ... so that comparisons with optci are easy e.g. through vimdiff
+    
+== 11/02/2013 == JL
+- Updated moist convection scheme to handle situations with a large water vapor content 
+- Added a keyword to enable ocean runoff in callphys.def (activerunoff)
+
+ 
+
Index: trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 874)
+++ trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 875)
@@ -3,4 +3,5 @@
      albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
 
+  use ioipsl_getincom 
   use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
   USE surfdat_h
@@ -55,7 +56,7 @@
       parameter (oceantime=10*24*3600)
 
-      logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
-      logical activerunoff ! enable simple runoff scheme?
-      logical oceanalbvary ! ocean albedo varies with the diurnal cycle?
+      logical,save :: oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
+      logical,save :: activerunoff ! enable simple runoff scheme?
+      logical,save :: oceanalbvary ! ocean albedo varies with the diurnal cycle?
 
 !     Arguments
@@ -97,10 +98,16 @@
       data firstcall /.true./
 
-      oceanbulkavg=.false.
-      activerunoff=.false.
-      oceanalbvary=.false.
 
       if(firstcall)then
 
+         oceanbulkavg=.false.
+         oceanalbvary=.false.
+         write(*,*)"Activate runnoff into oceans?"
+         activerunoff=.false.
+         call getin("activerunoff",activerunoff)
+         write(*,*)" activerunoff = ",activerunoff
+	 
+	 
+	 
          ALLOCATE(runoff(ngrid))
 
Index: trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 874)
+++ trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 875)
@@ -50,11 +50,13 @@
       EXTERNAL CBRT
       INTEGER i, k , nn
-      INTEGER,PARAMETER :: nitermax=1000
-      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6
+      INTEGER,PARAMETER :: nitermax=2000
+      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6 
+      ! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by
+      !                   decreasing alpha and increasing nitermax accordingly
       REAL zt(ngrid), zq(ngrid)
       REAL zcond(ngrid),zcond_iter
       REAL zdelq(ngrid)
       REAL zqs(ngrid), zdqs(ngrid)
-      REAL psat_tmp
+      REAL psat_tmp,dlnpsat_tmp
       
 ! evaporation calculations
@@ -115,5 +117,5 @@
             Do nn=1,nitermax  
                call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
-	       call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
+	       call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
                zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i))	   
                   !zcond can be negative here
@@ -131,5 +133,5 @@
             zcond(i) = 0.0
             Do nn=1,nitermax  
-	       call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
+	       call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
                zcond_iter = MAX(0.0,alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i)))	   
                   !zcond always postive! cannot evaporate clouds!
Index: trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 874)
+++ trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 875)
@@ -1,5 +1,5 @@
 subroutine moistadj(ngrid, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
 
-  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
+  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, RCPV, Psat_water, Lcpdqsat_water
   USE tracer_h
 
@@ -59,6 +59,6 @@
       REAL cp_delta_t(nlayermx)
       REAL v_cptj(nlayermx), v_cptjk1, v_ssig
-      REAL v_cptt(ngrid,nlayermx), v_p, v_t, zpsat
-      REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx)
+      REAL v_cptt(ngrid,nlayermx), v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
+      REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx),zpsat(ngrid,nlayermx),zdlnpsat(ngrid,nlayermx)
       REAL zq1(ngrid), zq2(ngrid)
       REAL gamcpdz(ngrid,2:nlayermx)
@@ -122,6 +122,6 @@
             v_p = pplay(i,k)
 
-            call Psat_water(v_t,v_p,zpsat,zqs(i,k))
-	    call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
+            call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
+	    call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
          ENDDO
       ENDDO
@@ -132,8 +132,15 @@
             zdp = pplev(i,k)-pplev(i,k+1)
             zdpm = pplev(i,k-1)-pplev(i,k)
-            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
-                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
-	        * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
-                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+!            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
+!                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
+!	        * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
+!                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+! general case where water is not a trace gas (JL13)
+            v_zqs     = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp)
+	    v_cptt2   = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp)
+	    v_pratio  = ((1./(1.+zpsat(i,k-1)/pplay(i,k-1)))*zdpm + (1./(1.+zpsat(i,k)/pplay(i,k)))*zdp)/(zdpm+zdp)
+	    v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
+            gamcpdz(i,k) = v_pratio*( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
+                / ((1.- v_zqs) + v_zqs * RCPV/RCPD + v_zqs * v_pratio * v_dlnpsat)                    
          ENDDO
       ENDDO
@@ -200,6 +207,6 @@
          v_p = pplay(i,k)
          
-         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
-         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
+         call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
+         call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
 
 
@@ -226,8 +233,15 @@
          zdpm = pplev(i,k-1) - pplev(i,k)
          zdp = pplev(i,k) - pplev(i,k+1)
-         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)             &
-                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
-                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
-                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+!         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
+!             +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
+!	      * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
+!             / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+! general case where water is not a trace gas (JL13)
+         v_zqs     = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp)
+	 v_cptt2   = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp)
+	 v_pratio  = ((1./(1.+zpsat(i,k-1)/pplay(i,k-1)))*zdpm + (1./(1.+zpsat(i,k)/pplay(i,k)))*zdp)/(zdpm+zdp)
+	 v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
+         gamcpdz(i,k) = v_pratio*( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
+                / ((1.- v_zqs) + v_zqs * RCPV/RCPD + v_zqs * v_pratio * v_dlnpsat)                    
       ENDDO
 
Index: trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90	(revision 874)
+++ trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90	(revision 875)
@@ -1,6 +1,4 @@
       subroutine su_watercycle
 
-! to use  'getin'
-      use ioipsl_getincom 
       use watercommon_h
       implicit none
@@ -28,6 +26,5 @@
       RV = 1000.*8.314/mH2O ! caution! R is R/mugaz already!
 
-      RCPV   = 4. *RV	 ! could be more precise...
-      !RCPV   = 5/2 *8.31*1000/mH2O ! BC modif here: doesn't appear to be correct
+      RCPV   = 1.88e3 ! specific heat capacity of water vapor at 350K
 
       RVTMP2 = RCPV/RCPD-1. ! not currently used...
Index: trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 874)
+++ trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 875)
@@ -77,5 +77,5 @@
 
 !==================================================================
-      subroutine Lcpdqsat_water(T,p,psat,qsat,dqsat)
+      subroutine Lcpdqsat_water(T,p,psat,qsat,dqsat,dlnpsat)
 
          implicit none
@@ -84,5 +84,5 @@
 !     Purpose
 !     -------
-!     Compute L/cp*d (q_sat)/d T
+!     Compute dqsat=L/cp*d (q_sat)/d T and dlnpsat=L/cp d(ln Psat)/d T
 !     for a given temperature (K)! 
 !     Based on the Tetens formula from L.Li physical parametrization manual
@@ -98,5 +98,5 @@
   
 !        output
-         real dqsat
+         real dqsat,dlnpsat
 
 ! JL12 variables for tetens formula
@@ -125,4 +125,5 @@
 
          dqsat=RLVTT/RCPD*qsat*(p/(p-(1.-epsi)*psat))*dummy
+	 dlnpsat=RLVTT/RCPD*dummy
          return
       end subroutine Lcpdqsat_water
