Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 1006)
+++ trunk/LMDZ.GENERIC/README	(revision 1016)
@@ -940,4 +940,5 @@
 == 15/05/2013 == JL
 - correction in radiative scheme to enforce double precision
+- corrects calculation of ISR
 
 == 22/05/2013 == EM
@@ -991,2 +992,8 @@
   compute the sponge quenching analytically rather than via Forward Euler
   approximation.
+
+== 07/08/2013 == JL
+- Water cycle in double precision (largescale+moistadj)
+- Improved wate rayleigh.
+- First step for rayleigh with variable species. Now, just need to change optcv.
+- changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers
Index: trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 1016)
@@ -13,4 +13,5 @@
 !     ------- 
 !     Robin Wordsworth (2010)
+!     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
 !
 !     Called by
@@ -25,5 +26,5 @@
 
       use radinc_h, only: L_NSPECTV
-      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep
+      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
       use gases_h
 
@@ -39,5 +40,5 @@
 
       logical typeknown
-      real*8 tauvar,tausum,tauwei,bwidth,bstart
+      real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart
       double precision df
 
@@ -54,7 +55,7 @@
       do igas=1,ngasmx
          if(igas.eq.vgas)then
-            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
-            'as it is variable.'
-         elseif(gfrac(igas).lt.5.e-2)then
+            print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
+	 endif
+	 if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
             print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
             'as its mixing ratio is less than 0.05.' 
@@ -68,5 +69,6 @@
                tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
             elseif(igas.eq.igas_H2O)then
-               tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 
+!               tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 
+               tauconsti(igas) = 1.98017E-10/(g*mugaz*100.) 
             elseif(igas.eq.igas_H2)then
                tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
@@ -83,5 +85,5 @@
             endif
 
-            if(gfrac(igas).eq.1.0)then
+            if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then
                print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.'
                typeknown=.true.
@@ -101,4 +103,6 @@
          tausum = 0.0
          tauwei = 0.0
+         tausumvar = 0.0
+         tauweivar = 0.0
          bstart = 10000.0/BWNV(N+1)
          bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
@@ -107,6 +111,7 @@
 
             tauvar=0.0
+            tauvarvar=0.0
             do igas=1,ngasmx
-               if((igas.eq.vgas).or.(gfrac(igas).lt.5.e-2))then
+               if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
                   ! ignore variable gas in Rayleigh calculation
                   ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
@@ -118,5 +123,6 @@
                      tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
                   elseif(igas.eq.igas_H2O)then
-                     tauvari(igas) = 1.0/wl**4 ! to be improved...
+!                     tauvari(igas) = 1.0/wl**4 ! to be improved...
+                     tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4 
                   elseif(igas.eq.igas_H2)then
                      tauvari(igas) = 1.0/wl**4 
@@ -130,5 +136,10 @@
                endif
 
-               tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas)
+               if(igas.eq.vgas) then
+                  tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas)
+                  tauvar=tauvar+0.0*0.0*gfrac(igas)
+               else
+                  tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas)
+               endif
 
             enddo
@@ -138,7 +149,10 @@
             tauwei=tauwei+df
             tausum=tausum+tauvar*df
+            tauweivar=tauweivar+df
+            tausumvar=tausumvar+tauvarvar*df
          
          enddo
          TAURAY(N)=tausum/tauwei
+         TAURAYVAR(N)=tausumvar/tauweivar
          ! we multiply by scalep here (100) because plev, which is used in optcv,
          ! is in units of mBar, so we need to convert.
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1016)
@@ -460,4 +460,8 @@
       endif
 
+!!! JL13: in the following, I changed some indices in the interpolations so that the model rsults are less dependent on the number of layers
+!!!    the older verions are commented with the commetn !JL13index 
+
+
 !-----------------------------------------------------------------------
 !     Water vapour (to be generalised for other gases eventually)
@@ -468,6 +472,7 @@
          do l=1,nlayer
             qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
-            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2    
-            ! Average approximation as for temperature...
+            qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)    
+!JL13index            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2    
+!JL13index            ! Average approximation as for temperature...
          end do
          qvar(1)=qvar(2)
@@ -497,12 +502,10 @@
          if(RH.lt.0.0) RH=0.0
 
-         ptemp = pplev(ig,1)
-         Ttemp = tsurf(ig)
-         call watersat(Ttemp,ptemp,qsat)
-
-         !qvar(2*nlayermx+1)=qsat      ! fully saturated everywhere
+!         ptemp = pplev(ig,1)
+!         Ttemp = tsurf(ig)
+!         call watersat(Ttemp,ptemp,qsat)
+
          qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
-         !qvar=0.005                   ! completely constant profile (JL)
-
+ 
       else
          do k=1,L_LEVELS
@@ -528,50 +531,54 @@
          END DO
 
-         !do k=1,L_LEVELS
-         !   qvar(k) = 0.0
-         !end do
-         !print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
-      endif
-
-
-      if(kastprof.and.(ngasmx.gt.1))then
-
+         if(ngasmx.gt.1)then
+
+            DO l=1,nlayer
+               muvarrad(2*l)   = muvar(ig,nlayer+2-l)
+               muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
+                                muvar(ig,max(nlayer+1-l,1)))/2
+            END DO
+      
+            muvarrad(1) = muvarrad(2)
+            muvarrad(2*nlayermx+1)=muvar(ig,1)
+
+            print*,'Recalculating qvar with VARIABLE epsi for kastprof'
+            print*,'Assumes that the variable gas is H2O!!!'
+            print*,'Assumes that there is only one tracer'
+            !i_var=igcm_h2o_vap
+            i_var=1
+            if(nq.gt.1)then
+               print*,'Need 1 tracer only to run kcm1d.e' 
+               stop
+            endif
+            do l=1,nlayer
+               vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi)) 
+               !vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O !JL to be changed
+            end do
+
+            do l=1,nlayer
+               qvar(2*l)   = vtmp(nlayer+1-l)
+               qvar(2*l+1) = vtmp(nlayer+1-l)
+!               qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
+            end do
+            qvar(1)=qvar(2)
+
+            print*,'Warning: reducing qvar in callcorrk.'
+            print*,'Temperature profile no longer consistent ', &
+                            'with saturated H2O. qsat=',satval
+            do k=1,L_LEVELS
+               qvar(k) = qvar(k)*satval
+            end do
+
+         endif
+      else ! if kastprof
          DO l=1,nlayer
             muvarrad(2*l)   = muvar(ig,nlayer+2-l)
-            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
-                                muvar(ig,max(nlayer+1-l,1)))/2
+            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
          END DO
       
          muvarrad(1) = muvarrad(2)
-         muvarrad(2*nlayermx+1)=muvar(ig,1)
-
-         print*,'Recalculating qvar with VARIABLE epsi for kastprof'
-         print*,'Assumes that the variable gas is H2O!!!'
-         print*,'Assumes that there is only one tracer'
-         !i_var=igcm_h2o_vap
-         i_var=1
-         if(nq.gt.1)then
-            print*,'Need 1 tracer only to run kcm1d.e' 
-            stop
-         endif
-         do l=1,nlayer
-            vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
-         end do
-
-         do l=1,nlayer
-            qvar(2*l)   = vtmp(nlayer+1-l)
-            qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
-         end do
-         qvar(1)=qvar(2)
-
-         print*,'Warning: reducing qvar in callcorrk.'
-         print*,'Temperature profile no longer consistent ', &
-                            'with saturated H2O.'
-         do k=1,L_LEVELS
-            qvar(k) = qvar(k)*satval
-         end do
-
+         muvarrad(2*nlayermx+1)=muvar(ig,1)         
       endif
-
+      
       ! Keep values inside limits for which we have radiative transfer coefficients
       if(L_REFVAR.gt.1)then ! there was a bug here!
@@ -607,11 +614,17 @@
       
       DO l=1,L_NLAYRAD-1
-         tmid(2*l+1) = tlevrad(2*l+1)
-         tmid(2*l+2) = tlevrad(2*l+1)
-         pmid(2*l+1) = plevrad(2*l+1)
-         pmid(2*l+2) = plevrad(2*l+1)
+         tmid(2*l+1) = tlevrad(2*l)
+         tmid(2*l+2) = tlevrad(2*l)
+         pmid(2*l+1) = plevrad(2*l)
+         pmid(2*l+2) = plevrad(2*l)
+!JL13index         tmid(2*l+1) = tlevrad(2*l+1)
+!JL13index         tmid(2*l+2) = tlevrad(2*l+1)
+!JL13index         pmid(2*l+1) = plevrad(2*l+1)
+!JL13index         pmid(2*l+2) = plevrad(2*l+1)
       END DO
-      pmid(L_LEVELS) = plevrad(L_LEVELS)
-      tmid(L_LEVELS) = tlevrad(L_LEVELS)
+      pmid(L_LEVELS) = plevrad(L_LEVELS-1)
+      tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
+!JL13index      pmid(L_LEVELS) = plevrad(L_LEVELS)
+!JL13index      tmid(L_LEVELS) = tlevrad(L_LEVELS)
 
       ! test for out-of-bounds pressure
@@ -635,11 +648,30 @@
             !print*,'WARNING, OVERRIDING FOR TEST'
             call abort
+	    !tlevrad(k)=tgasmin
          elseif(tlevrad(k).gt.tgasmax)then
             print*,'Maximum temperature is outside the radiative'
             print*,'transfer kmatrix bounds, exiting.'
 	    print*,'level,grid,T',k,ig,tlevrad(k)
-            print*,'WARNING, OVERRIDING FOR TEST'
-            !call abort
-	    tlevrad(k)=tgasmax
+            !print*,'WARNING, OVERRIDING FOR TEST'
+            call abort
+	    !tlevrad(k)=tgasmax
+         endif
+      enddo
+      do k=1,L_NLAYRAD+1
+         if(tmid(k).lt.tgasmin)then
+            print*,'Minimum temperature is outside the radiative'
+            print*,'transfer kmatrix bounds, exiting.'
+            print*,"k=",k," tlevrad(k)=",tlevrad(k)
+            print*,"tgasmin=",tgasmin
+            !print*,'WARNING, OVERRIDING FOR TEST'
+            call abort
+	    !tmid(k)=tgasmin
+         elseif(tmid(k).gt.tgasmax)then
+            print*,'Maximum temperature is outside the radiative'
+            print*,'transfer kmatrix bounds, exiting.'
+	    print*,'level,grid,T',k,ig,tlevrad(k)
+            !print*,'WARNING, OVERRIDING FOR TEST'
+            call abort
+	    !tmid(k)=tgasmax
          endif
       enddo
Index: trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F	(revision 1016)
@@ -105,5 +105,4 @@
       ! -- and stabilizes integrations
       NT    = int(TLEV(2*L+1)*NTfac) - NTstar+1
-
       !! For deep, opaque, thick first layers (e.g. Saturn)
       !! what is below works much better, not unstable, ...
Index: trunk/LMDZ.GENERIC/libf/phystd/largescale.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/largescale.F90	(revision 1016)
@@ -2,6 +2,9 @@
                         pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
 
+
+!     to use  'getin'
+      use ioipsl_getincom 
       use watercommon_h, only : RLVTT, RCPD, RVTMP2,  &
-          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
+          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_waterDP,Lcpdqsat_waterDP
       USE tracer_h
       IMPLICIT none
@@ -33,5 +36,5 @@
       REAL pplay(ngrid,nlayermx)   ! pression au milieu de couche
       REAL pt(ngrid,nlayermx)      ! temperature (K)
-      real pq(ngrid,nlayermx,nq) ! tracer mixing ratio (kg/kg)
+      REAL pq(ngrid,nlayermx,nq) ! tracer mixing ratio (kg/kg)
       REAL pdt(ngrid,nlayermx)     ! physical temperature tenedency (K/s)
       REAL pdq(ngrid,nlayermx,nq)! physical tracer tenedency (K/s)
@@ -43,6 +46,5 @@
 
 !     Options du programme
-      REAL ratqs   ! determine largeur de la distribution de vapeur
-      PARAMETER (ratqs=0.2)
+      REAL, SAVE :: ratqs   ! determine largeur de la distribution de vapeur
 
 !     Variables locales
@@ -50,13 +52,13 @@
       EXTERNAL CBRT
       INTEGER i, k , nn
-      INTEGER,PARAMETER :: nitermax=2000
-      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6 
+      INTEGER,PARAMETER :: nitermax=5000
+      DOUBLE PRECISION,PARAMETER :: alpha=.1,qthreshold=1.d-8
       ! 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,dlnpsat_tmp
+      DOUBLE PRECISION zt(ngrid), zq(ngrid)
+      DOUBLE PRECISION zcond(ngrid),zcond_iter
+      DOUBLE PRECISION zdelq(ngrid)
+      DOUBLE PRECISION zqs(ngrid), zdqs(ngrid)
+      DOUBLE PRECISION local_p,psat_tmp,dlnpsat_tmp,Lcp
       
 ! evaporation calculations
@@ -65,7 +67,17 @@
       REAL tevap(ngrid,nlayermx)
 
-      REAL zcor(ngrid), zdelta(ngrid), zcvm5(ngrid)
-      REAL zx_q(ngrid)
-      REAL Nmix_local,zfice
+      DOUBLE PRECISION zx_q(ngrid)
+      LOGICAL,SAVE :: firstcall=.true.
+
+
+      IF (firstcall) THEN
+
+	 write(*,*) "value for ratqs? "
+         ratqs=0.2 ! default value
+         call getin("ratqs",ratqs)
+         write(*,*) " ratqs = ",ratqs
+
+         firstcall = .false.
+      ENDIF
 
 !     GCM -----> subroutine variables, initialisation of outputs
@@ -75,4 +87,5 @@
       pdqliqlsc(1:ngrid,1:nlayermx) = 0.0
       rneb(1:ngrid,1:nlayermx) = 0.0
+      Lcp=RLVTT/RCPD
 
 
@@ -93,13 +106,13 @@
       DO i = 1, ngrid
 
+         local_p=pplay(i,k)
          if(zt(i).le.15.) then
 	    print*,'in lsc',i,k,zt(i)
 !	    zt(i)=15.   ! check too low temperatures
          endif
-         call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
+         call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i))
  
-         zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1.e-12)
+         zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1.d-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
 
@@ -107,6 +120,5 @@
 	    zcond(i)=0.
 
-	 else if (rneb(i,k).gt.1.) then    !complete cloud cover, we start without evaporating
-
+	 else if ((rneb(i,k).gt.0.99).or.(ratqs.lt.1.e-6)) then    !complete cloud cover, we start without evaporating
 	    rneb(i,k)=1.
             zt(i)=pt(i,k)+pdt(i,k)*ptimestep
@@ -114,14 +126,16 @@
 	    dqevap(i,k)=0.
 !           iterative process to stabilize the scheme when large water amounts JL12
-            zcond(i) = 0.0
+            zcond(i) = 0.0d0
             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),dlnpsat_tmp)
-               zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i))	   
+               call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i))
+	       call Lcpdqsat_waterDP(zt(i),local_p,psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
+               zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+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
+	       zt(i) = zt(i) + zcond_iter*Lcp
+	       if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
+!	       if (ABS(zcond_iter/alpha).lt.qthreshold) exit
+	       if (nn.eq.nitermax) print*,'itermax in largescale'
 	    End do ! niter
 	    zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_h2o_ice)+pdq(i,k,igcm_h2o_ice)*ptimestep))
@@ -129,17 +143,19 @@
 	 else   !standard case	    
 
-	    zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0 !water vapor in cloudy sky
+	    zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0d0 !water vapor in cloudy sky
 !           iterative process to stabilize the scheme when large water amounts JL12
-            zcond(i) = 0.0
+            zcond(i) = 0.0d0
             Do nn=1,nitermax  
-	       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)))	   
+	       call Lcpdqsat_waterDP(zt(i),local_p,psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
+               zcond_iter = MAX(0.0d0,alpha*(zx_q(i)-zqs(i))/(1.d0+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*rneb(i,k)
-               call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
+	       if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
+!	       if (ABS(zcond_iter/alpha).lt.qthreshold) exit
+	       zt(i) = zt(i) + zcond_iter*Lcp*rneb(i,k)
+               call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i))
+	       if (nn.eq.nitermax) print*,'itermax in largescale'
 	    End do ! niter
 
@@ -153,15 +169,8 @@
          pdqvaplsc(1:ngrid,k)  = dqevap(1:ngrid,k) - zcond(1:ngrid)
          pdqliqlsc(1:ngrid,k) = - pdqvaplsc(1:ngrid,k)
-         pdtlsc(1:ngrid,k)  = pdqliqlsc(1:ngrid,k)*RLVTT/RCPD
+         pdtlsc(1:ngrid,k)  = pdqliqlsc(1:ngrid,k)*real(Lcp)
 
    Enddo ! k= nlayermx, 1, -1
-   
-      !print*,'qsat=',zqs
-      !print*,'q=',q
-      !print*,'dq=',pdqvaplsc*ptimestep
-      !print*,'dT in LS=',pdtlsc*ptimestep
-
-      !print*,'rice=',rice
-      !print*,'rneb=',rneb
+ 
 
       return
Index: trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 1016)
@@ -56,17 +56,17 @@
       LOGICAL itest(ngrid)
       REAL delta_q(ngrid, nlayermx)
-      REAL cp_new_t(nlayermx)
+      DOUBLE PRECISION :: cp_new_t(nlayermx), v_cptt(ngrid,nlayermx)
       REAL cp_delta_t(nlayermx)
-      REAL v_cptj(nlayermx), v_cptjk1, v_ssig
-      REAL v_cptt(ngrid,nlayermx), v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
+      DOUBLE PRECISION :: v_cptj(nlayermx), v_cptjk1, v_ssig
+      REAL 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)
-      REAL zdp, zdpm
+      DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayermx)
+      DOUBLE PRECISION :: zdp, zdpm
 
       REAL zsat ! super-saturation
       REAL zflo ! flotabilite
 
-      REAL local_q(ngrid,nlayermx),local_t(ngrid,nlayermx)
+      DOUBLE PRECISION :: local_q(ngrid,nlayermx),local_t(ngrid,nlayermx)
 
       REAL zdelta, zcor, zcvm5
@@ -139,8 +139,8 @@
             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_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (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)                    
+            gamcpdz(i,k) = ( (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_pratio + v_zqs  * v_dlnpsat)                
          ENDDO
       ENDDO
@@ -210,7 +210,4 @@
          call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
 
-
-
-!          print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k)
       ENDDO
     Enddo ! nn=1,niter
@@ -240,8 +237,8 @@
          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_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (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)                    
+         gamcpdz(i,k) = ( (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_pratio + v_zqs  * v_dlnpsat)                
       ENDDO
 
@@ -271,11 +268,4 @@
  9999 CONTINUE ! loop over all the points
 
-!      print*,'k1=',k1
-!      print*,'k2=',k2
-
-!      print*,'local_t=',local_t
-!      print*,'v_cptt=',v_cptt
-!      print*,'gamcpdz=',gamcpdz
-
 !-----------------------------------------------------------------------
 ! Determine the cloud fraction (hypothese: la nebulosite a lieu
@@ -317,6 +307,4 @@
       ENDDO
       ENDDO
-
-!      print*,'local_q BEFORE=',local_q
 
       DO k = 1, nlayermx
Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 1016)
@@ -121,8 +121,9 @@
      if(kastprof)then
         dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
-        U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
+        U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
      else
-        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
-        U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
+        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))*mugaz/muvar(k)
+        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
+	    !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
      endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 1016)
@@ -108,9 +108,10 @@
      ! if we have continuum opacities, we need dz
      if(kastprof)then
-        dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K))
-        U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k) 
+        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
+        U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
      else
-        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
-        U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
+        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))*mugaz/muvar(k)
+        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
+	    !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
      endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1016)
@@ -9,5 +9,5 @@
  
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
-      use watercommon_h, only : RLVTT, Psat_water
+      use watercommon_h, only : RLVTT, Psat_water,epsi
       use gases_h
       use radcommon_h, only: sigma
@@ -329,5 +329,5 @@
       real dEzRadsw(ngrid,nlayermx),dEzRadlw(ngrid,nlayermx),dEzdiff(ngrid,nlayermx)
       real dEdiffs(ngrid),dEdiff(ngrid)
-      real madjdE(ngrid), lscaledE(ngrid)
+      real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayermx), lscaledEz(ngrid,nlayermx)
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
       
@@ -711,5 +711,11 @@
                  call abort
               endif
-              muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
+	      if(water) then
+                  muvar(1:ngrid,1:nlayermx)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayermx,igcm_h2o_vap)) 
+		  muvar(1:ngrid,nlayermx+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayermx,igcm_h2o_vap)) 
+		  ! take into account water effect on mean molecular weight
+	      else
+                  muvar(1:ngrid,1:nlayermx+1)=mugaz
+	      endif 
       
 !             standard callcorrk
@@ -1059,5 +1065,5 @@
 	       dtmoist(1:ngrid,1:nlayermx)=0.
 
-               call moistadj(ngrid,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
+	       call moistadj(ngrid,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
 
                pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap)   &
@@ -1071,4 +1077,5 @@
                if(enertest)then
                   dEtot=cpp*SUM(massarea(:,:)*dtmoist(:,:))/totarea
+		  madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:)
                   do ig=1,ngrid
                      madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
@@ -1089,5 +1096,4 @@
 !        Large scale condensation/evaporation
 !        --------------------------------
-
                call largescale(ngrid,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
 
@@ -1099,4 +1105,5 @@
                ! test energy conservation
                if(enertest)then
+                  lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:)
                   do ig=1,ngrid
                      lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
@@ -1534,7 +1541,5 @@
          do l = 1, nlayer
             do ig=1,ngrid
-!               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
                call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
-
                RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
             enddo
@@ -1698,4 +1703,10 @@
            call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
            call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
+!           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
+!           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
+!           call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
+!           call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
+!           call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
+!           call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
            call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
@@ -1715,8 +1726,9 @@
 	  endif 
           if(watercond) then
+!	     call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 
+!	     call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
 	     call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
 	     call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
-	     call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
-!	     call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
+	     call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
 !	     call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
           endif 
@@ -1759,4 +1771,5 @@
              call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
              call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
+	     call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
           endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 1016)
@@ -36,4 +36,6 @@
 !     TAURAY     - Array (NSPECTV elements) of the pressure-independent
 !                  part of Rayleigh scattering optical depth.
+!     TAURAYVAR  - Array (NSPECTV elements) of the pressure-independent
+!                  part of Rayleigh scattering optical depth for the variable gas.
 !     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
 !                  each temperature, pressure, and spectral interval
@@ -61,5 +63,5 @@
       REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI)
       REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV)
-      REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV)
+      REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV)
 
       REAL*8 blami(L_NSPECTI+1)
Index: trunk/LMDZ.GENERIC/libf/phystd/rain.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 1016)
@@ -66,5 +66,5 @@
       INTEGER,PARAMETER :: ninter=5
 
-      logical,parameter :: evap_prec=.true. ! Does the rain evaporate?
+      logical,save :: evap_prec ! Does the rain evaporate?
 
 !     for simple scheme
@@ -142,4 +142,9 @@
 
 	 endif
+
+         write(*,*) "re-evaporate precipitations?"
+         evap_prec=.true. ! default value
+         call getin("evap_prec",evap_prec)
+         write(*,*) " evap_prec = ",evap_prec
 
          firstcall = .false.
Index: trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1016)
@@ -543,9 +543,9 @@
                   zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
 
-                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
+                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
 		      + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
                       + zfluxq(ig,1)*dryness(ig)*RLVTT*((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
 
-                  z2(ig) = pcapcal(ig) + cpp*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
+                  z2(ig) = pcapcal(ig) + cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
                       + zdplanck(ig)+zfluxq(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*(1.0-zdq(ig,1))
 
Index: trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 1016)
@@ -175,4 +175,111 @@
       end subroutine Tsat_water
 
+!==================================================================
+      subroutine Psat_waterDP(T,p,psat,qsat)
+
+         implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the saturation vapor pressure and mass mixing ratio at saturation (kg/kg)
+!     for a given pressure (Pa) and temperature (K)
+!     DOUBLE PRECISION
+!     Based on the Tetens formula from L.Li physical parametrization manual
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+
+!        input
+         double precision, intent(in) :: T, p
+  
+!        output
+         double precision psat,qsat
+
+! JL12 variables for tetens formula
+         double precision,parameter :: Pref_solid_liquid=611.14d0
+         double precision,parameter :: Trefvaporization=35.86D0
+         double precision,parameter :: Trefsublimation=7.66d0
+         double precision,parameter :: Tmin=8.d0
+         double precision,parameter :: r3vaporization=17.269d0
+         double precision,parameter :: r3sublimation=21.875d0
+
+! checked vs. old watersat data 14/05/2012 by JL.
+
+         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
+            qsat=1.d0
+         else
+            qsat=epsi*psat/(p-(1.d0-epsi)*psat)
+         endif
+         return
+      end subroutine Psat_waterDP
+
+
+
+
+!==================================================================
+      subroutine Lcpdqsat_waterDP(T,p,psat,qsat,dqsat,dlnpsat)
+
+         implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     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
+!
+!     Authors
+!     -------
+!     Jeremy Leconte (2012)
+!
+!==================================================================
+
+!        input
+         double precision T, p, psat, qsat
+  
+!        output
+         double precision dqsat,dlnpsat
+
+! JL12 variables for tetens formula
+         double precision,parameter :: Pref_solid_liquid=611.14d0
+         double precision,parameter :: Trefvaporization=35.86d0
+         double precision,parameter :: Tmin=8.d0
+         double precision,parameter :: Trefsublimation=7.66d0
+         double precision,parameter :: r3vaporization=17.269d0
+         double precision,parameter :: r3sublimation=21.875d0
+
+         double precision :: dummy
+
+         if (psat.gt.p) then
+	    dqsat=0.d0
+	    return
+	 endif
+
+         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*(p/(p-(1.d0-epsi)*psat))*dummy
+	 dlnpsat=RLVTT/RCPD*dummy
+         return
+      end subroutine Lcpdqsat_waterDP
+
 
 end module watercommon_h
