Index: trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90	(revision 785)
+++ trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90	(revision 786)
@@ -152,5 +152,5 @@
          band=1
 	 Do while(lam(1).lt. real(10000.0/BWNV(band+1)))
-	    if (band+1.ge.L_NSPECTV+1) exit
+	    if (band.gt.L_NSPECTV-1) exit
             band=band+1
 	 enddo
Index: trunk/LMDZ.GENERIC/libf/phystd/dimphys.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/dimphys.h	(revision 785)
+++ trunk/LMDZ.GENERIC/libf/phystd/dimphys.h	(revision 786)
@@ -9,4 +9,5 @@
 ! nsoilmx : number of subterranean layers
       !integer, parameter :: nsoilmx = 4 ! for a test 
-      integer, parameter :: nsoilmx = 18 ! for z1=0.02 cm, depth = 104.8 m 
+      integer, parameter :: nsoilmx = 18 ! for z1=0.0002 cm, depth = 18 m => mars case 
+      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 cm, depth = 104.8 m => earth case
 !-----------------------------------------------------------------------
Index: trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 785)
+++ trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 786)
@@ -51,5 +51,6 @@
       EXTERNAL CBRT
       INTEGER i, k , nn
-      INTEGER,PARAMETER :: niter=4
+      INTEGER,PARAMETER :: nitermax=1000
+      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6
       REAL zt(ngridmx), zq(ngridmx)
       REAL zcond(ngridmx),zcond_iter
@@ -92,43 +93,58 @@
 
          if(zt(i).le.15.) then
-	    print*,'in lsc',i,zt(i)
-	    zt(i)=15.   ! check too low temperatures
+	    print*,'in lsc',i,k,zt(i)
+!	    zt(i)=15.   ! check too low temperatures
          endif
-!         call watersat(zt(i),pplay(i,k),zqs(i))
          call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
  
-         zdelq(i) = ratqs * zq(i)
-	 if(zq(i)+zdelq(i).lt.0.999) then
-	    rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
+         zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1.e-12)
+	 rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
+!	 print*,zq(i),zdelq(i),zqs(i),rneb(i,k)
+	 if (rneb(i,k).lt.0.) then  !no clouds
+
+	    rneb(i,k)=0.
+	    zcond(i)=0.
+
+	 else if (rneb(i,k).gt.1.) then    !complete cloud cover, we start without evaporating
+
+	    rneb(i,k)=1.
+            zt(i)=pt(i,k)+pdt(i,k)*ptimestep
+	    zx_q(i) = pq(i,k,igcm_h2o_vap)+pdq(i,k,igcm_h2o_vap)*ptimestep
+	    dqevap(i,k)=0.
+!           iterative process to stabilize the scheme when large water amounts JL12
+            zcond(i) = 0.0
+            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))
+               zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i))	   
+                  !zcond can be negative here
+               zx_q(i) = zx_q(i) - zcond_iter
+	       zcond(i) = zcond(i) + zcond_iter
+	       zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
+	       if (ABS(zcond_iter/alpha).lt.qthreshold) exit
+	    End do ! niter
+	    zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_h2o_ice)+pdq(i,k,igcm_h2o_ice)*ptimestep))
+
+	 else   !standard case	    
+
 	    zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0 !water vapor in cloudy sky
-            if (rneb(i,k) .LE. 0.0) zx_q(i) = 0.0
-            if (rneb(i,k) .GE. 1.0) zx_q(i) = zq(i)
-            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
-	 else
-	    if(zq(i).gt.zqs(i)) then
-	       rneb(i,k)=1.
-	       zx_q(i)=zq(i)
-	    else
-	       rneb(i,k)=0.
-	       zx_q(i)=zqs(i)  !no condensation needed
-	    Endif
+!           iterative process to stabilize the scheme when large water amounts JL12
+            zcond(i) = 0.0
+            Do nn=1,nitermax  
+	       call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
+               zcond_iter = MAX(0.0,alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i)))	   
+                  !zcond always postive! cannot evaporate clouds!
+                  !this is why we must reevaporate before largescale
+               zx_q(i) = zx_q(i) - zcond_iter
+	       zcond(i) = zcond(i) + zcond_iter
+	       if (ABS(zcond_iter/alpha).lt.qthreshold) exit
+	       zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
+               call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
+	    End do ! niter
+
 	 Endif
 
-!        iterative process to stabilize the scheme when large water amounts JL12
-         zcond(i) = 0.0
-         Do nn=1,niter  
-!            call watersat_grad(zt(i),zqs(i),zdqs(i))
-	    call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
-            zcond_iter = MAX(0.0,(zx_q(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i)))	   
-               !zcond always postive! cannot evaporate clouds!
-               !this is why we must reevaporate before largescale
-            zx_q(i) = zx_q(i) - zcond_iter
-	    zcond(i) = zcond(i) + zcond_iter
-	    zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
-            call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
-	 End do ! niter
+         zcond(i) = zcond(i)*rneb(i,k)/ptimestep ! JL12
 
-
-         zcond(i) = zcond(i)/ptimestep ! added by RDW
       ENDDO
 
Index: trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 785)
+++ trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 786)
@@ -146,5 +146,5 @@
 	          zmassboil(ig)=pqs(ig,igcm_h2o_vap)/ptimestep
 	       endif
-	       zmassboil(ig)=zmassboil(ig)*0.3 !momentary, should be 1. JL12
+	       zmassboil(ig)=zmassboil(ig)*0.0 !momentary, should be 1. JL12
                pdqsmr(ig,igcm_h2o_vap)=-zmassboil(ig)
 	       pdtsrfmr(ig)=-zmassboil(ig)*RLVTT/pcapcal(ig)
Index: trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F	(revision 785)
+++ trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F	(revision 786)
@@ -133,4 +133,7 @@
         ! read <depth> coordinate
         if (interpol) then !put values in oldmlayer
+          if (.not.allocated(oldmlayer)) then
+             allocate(oldmlayer(dimlen),stat=ierr)
+	  endif
 #ifdef NC_DOUBLE
            ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldmlayer)
@@ -159,5 +162,6 @@
       ! default mlayer distribution, following a power law:
       !  mlayer(k)=lay1*alpha**(k-1/2)
-        lay1=2.e-4
+        lay1=2.e-4  !mars case (nsoilmx=18)
+!        lay1=3.e-2  !earth case (nsoilmx=13)
 	alpha=2
         do iloop=0,nsoil-1
Index: trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 785)
+++ trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 786)
@@ -51,4 +51,5 @@
          real,parameter :: Trefvaporization=35.86
          real,parameter :: Trefsublimation=7.66
+         real,parameter :: Tmin=8.
          real,parameter :: r3vaporization=17.269
          real,parameter :: r3sublimation=21.875
@@ -56,8 +57,11 @@
 ! checked vs. old watersat data 14/05/2012 by JL.
 
-         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
-            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation))
-         else                 ! liquid / vapour
-            psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization))
+         if (T.gt.T_h2O_ice_liq) then 
+            psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization)) ! liquid / vapour
+         else if (T.lt.Tmin) then
+	    print*, "careful, T<Tmin in psat water"
+            psat = Pref_solid_liquid*Exp(r3sublimation*(Tmin-T_h2O_ice_liq)/(Tmin-Trefsublimation)) ! min psat  
+         else                 
+            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation)) ! solid / vapour
          endif
          if(psat.gt.p) then
@@ -99,4 +103,5 @@
          real,parameter :: Pref_solid_liquid=611.14
          real,parameter :: Trefvaporization=35.86
+         real,parameter :: Tmin=8.
          real,parameter :: Trefsublimation=7.66
          real,parameter :: r3vaporization=17.269
@@ -110,11 +115,14 @@
 	 endif
 
-         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
-            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2
-         else                 ! liquid / vapour
-            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2
+         if (T.gt.T_h2O_ice_liq) then
+            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2  ! liquid / vapour
+         else if (T.lt.Tmin) then
+	    print*, "careful, T<Tmin in Lcp psat water"
+            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(Tmin-Trefsublimation)**2  ! solid / vapour
+         else               
+            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2  ! solid / vapour
          endif
 
-         dqsat=RLVTT/RCPD*qsat**2*(p/(epsi*psat))*dummy
+         dqsat=RLVTT/RCPD*qsat*(p/(p-(1.-epsi)*psat))*dummy
          return
       end subroutine Lcpdqsat_water
