Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 649)
+++ /trunk/LMDZ.GENERIC/README	(revision 650)
@@ -661,3 +661,8 @@
    - rcp, cpp and mugaz can now be computed using gases.def in newstart
 - Correction of a bug arising in gcm.F when the solar days are long (thanks Melanie V.)
-
+- Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
+- Minor name changes in watercommon
+- Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk
+       (for radfixed=true)
+- Added consistency check in inifis
+- Moved 1d water initialization from physiqu to rcm1d
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 650)
@@ -172,7 +172,14 @@
 
 !     are particle radii fixed?
-      logical radfixed
+      logical, save ::  radfixed
       real maxrad, minrad
 
+!     Local water cloud optical properties
+      real, parameter ::  rad_froid=35.0e-6 
+      real, parameter ::  rad_chaud=10.0e-6 
+      real, parameter ::  coef_chaud=0.13
+      real, parameter ::  coef_froid=0.09
+      real zfice
+      
       real CBRT
       external CBRT
@@ -183,6 +190,4 @@
       REAL*8 muvarrad(L_LEVELS)
 
-      radfixed=.false.
-
 !=======================================================================
 !     Initialization on first call
@@ -195,7 +200,11 @@
       qsiaer(:,:,:)=0.0
       giaer(:,:,:) =0.0
+      radfixed=.false.
 
 
       if(firstcall) then
+         if(kastprof)then
+            radfixed=.true.
+         endif  
 
          call system('rm -f surf_vals_long.out')
@@ -305,7 +314,5 @@
       enddo
 
-      if(kastprof)then
-         radfixed=.true.
-      endif
+
 
       if(radfixed)then
@@ -320,5 +327,10 @@
                do ig=1,ngrid
                   !reffrad(ig,l,2) = 2.e-5 ! H2O ice
-                  reffrad(ig,l,2) = 5.e-6 ! H2O ice
+                   reffrad(ig,l,2) = rad_chaud ! H2O ice
+
+                   zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
+                   zfice = MIN(MAX(zfice,0.0),1.0)
+                   reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice
+                   nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice
                enddo
             enddo
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 650)
@@ -11,5 +11,4 @@
      &   , Nmix_co2, Nmix_h2o                                           &
      &   , dusttau                                                      &
-     &   , nonideal                                                     &
      &   , meanOLR                                                      &
      &   , specOLR                                                      &
Index: /trunk/LMDZ.GENERIC/libf/phystd/evap.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/evap.F	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/evap.F	(revision 650)
@@ -61,5 +61,5 @@
          ! ignoring qevap term creates huge difference when qevap large!!!
 
-         zdelta = MAX(0.,SIGN(1.,To-tevap(ig,l))) ! what is this?
+         zdelta = MAX(0.,SIGN(1.,T_h2O_ice_liq-tevap(ig,l))) ! what is this?
                                                   ! for division between water / liquid
 	 dqevap(ig,l) = MAX(0.0,qevap(ig,l,igcm_h2o_ice))/dtime
Index: /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 650)
@@ -3,5 +3,5 @@
      albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
 
-  use watercommon_h, only: To, RLFTT, rhowater, mx_eau_sol
+  use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
 
   implicit none
@@ -189,6 +189,6 @@
             hice(ig)    = 0.
 
-            if(twater(ig) .lt. To)then
-               E=min((To+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
+            if(twater(ig) .lt. T_h2O_ice_liq)then
+               E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
                hice(ig)        = E/(RLFTT*rhowater)
                hice(ig)        = max(hice(ig),0.0)
@@ -221,8 +221,8 @@
 
 !     melt the snow
-            if(ztsurf(ig).gt.To)then
+            if(ztsurf(ig).gt.T_h2O_ice_liq)then
                if(zqsurf(ig,iice).gt.1.0e-8)then
 
-                  a     = (ztsurf(ig)-To)*pcapcal(ig)/RLFTT
+                  a     = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
                   b     = zqsurf(ig,iice)
                   fsnoi = min(a,b)
@@ -240,5 +240,5 @@
                if(zqsurf(ig,iliq).gt.1.0e-8)then
 
-                  a     = -(ztsurf(ig)-To)*pcapcal(ig)/RLFTT
+                  a     = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
                   b     = zqsurf(ig,iliq)
                   
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 650)
@@ -219,12 +219,4 @@
          call getin("calladj",calladj)
          write(*,*) " calladj = ",calladj
-        
-
-! Test of incompatibility
-         if (enertest.and.nonideal) then
-            print*,'Energy conservation calculations currently
-     &           assume ideal gas!'
-            call abort
-         endif
 
          write(*,*) "call CO2 condensation ?"
@@ -232,5 +224,10 @@
          call getin("co2cond",co2cond)
          write(*,*) " co2cond = ",co2cond
-
+! Test of incompatibility
+         if (co2cond.and.(.not.tracer)) then
+            print*,'We need a CO2 ice tracer to condense CO2'
+            call abort
+         endif	 
+	 
          write(*,*) "CO2 supersaturation level ?"
          co2supsat=1.0 ! default value
Index: /trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90	(revision 650)
@@ -1,5 +1,5 @@
 subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb)
 
-  use watercommon_h, only: To, RLVTT, RCPD
+  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD
 
   implicit none
@@ -124,10 +124,4 @@
          DO i = 1, ngridmx
             v_cptt(i,k) = RCPD * local_t(i,k)
-         ENDDO
-      ENDDO
-
-      DO k = 1, nlayermx
-         DO i = 1, ngridmx
-            v_cptt(i,k) = RCPD * local_t(i,k)
             v_t = local_t(i,k)
             v_p = pplay(i,k)
@@ -142,5 +136,5 @@
 !         DO i = 1, ngridmx
 !            v_t = local_t(i,k)
-!            IF (v_t.LT.To) THEN
+!            IF (v_t.LT.T_h2O_ice_liq) THEN
 !               print*,'RHs=',q(i,k) / v_qs(i,k)
 !            ELSE
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 650)
@@ -8,5 +8,5 @@
  
       use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV
-      use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
+      use watercommon_h, only : RLVTT
       use gases_h
       use radcommon_h, only: sigma
@@ -370,6 +370,4 @@
       real vdifcncons(ngridmx)
       real cadjncons(ngridmx)
-      real qzero1D
-      save qzero1D
 
 !      double precision qsurf_hist(ngridmx,nqmx)
@@ -544,17 +542,6 @@
             ! initialise variables for water cycle
 
-            qzero1D=0.0
-
             if(enertest)then
                watertest = .true.
-            endif
-
-            if(ngrid.eq.1)then
-               qzero1D               = 0.0
-               qsurf(1,igcm_h2o_vap) = qzero1D
-               do l=1, nlayer
-                  pq(1,l,igcm_h2o_vap)=0.0
-                  pq(1,l,igcm_h2o_ice)=0.0
-               enddo
             endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/rain.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 650)
@@ -2,5 +2,5 @@
 
 
-  use watercommon_h, only: To, RLVTT, RCPD, RCPV, RV, RVTMP2
+  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2
 
   implicit none
@@ -82,5 +82,4 @@
 
       REAL zoliq(ngridmx)
-      REAL ztglace
       REAL zdz(ngridmx),zrho(ngridmx),ztot(ngridmx), zrhol(ngridmx)
       REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
@@ -140,7 +139,4 @@
       ENDDO
       ENDDO
-
-!     Determine the cold clouds by their temperature
-      ztglace = To - 15.0
 
 !     Initialise the outputs
@@ -219,5 +215,5 @@
             DO i = 1, ngridmx
                ttemp = zt(i,k)
-               IF (ttemp .ge. To) THEN
+               IF (ttemp .ge. T_h2O_ice_liq) THEN
                   lconvert=rainthreshold
                ELSEIF (ttemp .gt. t_crit) THEN
@@ -245,5 +241,5 @@
                   zrho(i)  = pplay(i,k) / ( zt(i,k) * R )
                   zdz(i)   = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
-                  zfrac(i) = (zt(i,k)-ztglace) / (To-ztglace)
+                  zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
                   zfrac(i) = MAX(zfrac(i), 0.0)
                   zfrac(i) = MIN(zfrac(i), 1.0)
@@ -291,5 +287,5 @@
             call abort
          endif
-         IF (t(i,1) .LT. To) THEN
+         IF (t(i,1) .LT. T_h2O_ice_liq) THEN
             dqssnow(i) = zrfl(i)
             dqsrain(i) = 0.0
Index: /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 650)
@@ -73,5 +73,6 @@
 !      REAL co2ice           ! co2ice layer (kg.m-2) !not used anymore
       integer :: i_co2_ice=0  ! tracer index of co2 ice
-      integer :: i_h2o_ice=0  ! tracer index of co2 ice
+      integer :: i_h2o_ice=0  ! tracer index of h2o ice
+      integer :: i_h2o_vap=0  ! tracer index of h2o vapor
       REAL emis             ! surface layer
       REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
@@ -111,4 +112,5 @@
       REAL cloudfrac(ngridmx,nlayermx)
       REAL hice(ngridmx),totcloudfrac(ngridmx)
+      REAL qzero1D   !initial water amount on the ground
 
 c=======================================================================
@@ -228,4 +230,5 @@
          i_co2_ice=0
          i_h2o_ice=0
+         i_h2o_vap=0
          do iq=1,nq
            if (tnom(iq)=="co2_ice") then
@@ -233,4 +236,6 @@
            elseif (tnom(iq)=="h2o_ice") then
              i_h2o_ice=iq
+           elseif (tnom(iq)=="h2o_vap") then
+             i_h2o_vap=iq
            endif
          enddo
@@ -275,5 +280,5 @@
       area(1)=1.E+0
       aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
-     !!! GEOPOTENTIAL. useless in 1D because control bu surface pressure
+     !!! GEOPOTENTIAL. useless in 1D because control by surface pressure
       phisfi(1)=0.E+0
      !!! LATITUDE. can be set. 
@@ -493,4 +498,9 @@
         qsurf(iq) = 0.
       ENDDO
+
+      if (water) then
+         qzero1D               = 0.0
+         qsurf(i_h2o_vap) = qzero1D
+      endif
 
 c   Initialisation pour prendre en compte les vents en 1-D
Index: /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 650)
@@ -7,5 +7,5 @@
           pdqdif,pdqsdif,lastcall)
 
-      use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol
+      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
       use radcommon_h, only : sigma
 
@@ -609,5 +609,5 @@
 	             else !dqsdif_total(ig).ge.0
                         !If water vapour is condensing, we must decide whether it forms ice or liquid.
-                        if(ztsrf(ig).gt.To)then
+                        if(ztsrf(ig).gt.T_h2O_ice_liq)then
                            pdqsdif(ig,iice_surf)=0.0
                            pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
Index: /trunk/LMDZ.GENERIC/libf/phystd/vdifc.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/vdifc.F	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/vdifc.F	(revision 650)
@@ -7,5 +7,5 @@
      &     pdqdif,pdqsdif,lastcall)
 
-      use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol
+      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
       use radcommon_h, only : sigma
 
@@ -644,5 +644,5 @@
 !     If water vapour is condensing, we must decide whether it forms ice or liquid.
                  if(dqsdif_total(ig).gt.0)then ! a bug was here!
-                    if(ztsrf2(ig).gt.To)then
+                    if(ztsrf2(ig).gt.T_h2O_ice_liq)then
                        pdqsdif(ig,iice)=0.0
                        pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
@@ -704,5 +704,5 @@
                Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice)
 
-               if(ztsrf2(ig).gt.To)then
+               if(ztsrf2(ig).gt.T_h2O_ice_liq)then
                   pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff
                else
Index: /trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90	(revision 650)
@@ -4,16 +4,18 @@
 
       real, parameter :: T_coup = 234.0
-      real, parameter :: To = 273.16
+      real, parameter :: T_h2O_ice_liq = 273.16
+      real, parameter :: T_h2O_ice_clouds = T_h2O_ice_liq-15. 
       real, parameter :: mH2O = 18.01528   
 
       ! benjamin additions
-      real, parameter :: RLVTT = 2.257E+6 ! Latent heat of sublimation (J kg-1) 
+      real, parameter :: RLVTT = 2.257E+6 ! Latent heat of vaporization (J kg-1) 
       real, parameter :: RLSTT = 2.257E+6 ! 2.591E+6 in reality ! Latent heat of sublimation (J kg-1)
       !real, parameter :: RLVTT = 0.0 
       !real, parameter :: RLSTT = 0.0
 
-      real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1) 
+      real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1) ! entails an energy sink but better description of albedo
       !real, parameter :: RLFTT = 0.0 
       real, parameter :: rhowater = 1.0E+3 ! mass of water (kg/m^3)
+      real, parameter :: capcal_h2o_liq = 4181.3 ! specific heat capacity of liquid water J/kg/K
       real, parameter :: mx_eau_sol = 150 ! mass of water (kg/m^2)
 
@@ -21,6 +23,5 @@
 !     lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
 
-      real epsi, RCPD, RCPV, RV, RVTMP2
-      save epsi, RCPD, RCPV, RV, RVTMP2
+      real, save :: epsi, RCPD, RCPV, RV, RVTMP2
 
       end module watercommon_h
Index: /trunk/LMDZ.GENERIC/libf/phystd/watersat.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/watersat.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/watersat.F90	(revision 650)
@@ -1,5 +1,5 @@
 subroutine watersat(T,p,qsat)
 
-  use watercommon_h, only: To, epsi
+  use watercommon_h, only: T_h2O_ice_liq, epsi
   implicit none
 
@@ -28,5 +28,5 @@
 ! x by epsi converts to mass mixing ratio
 
-  if (T.lt.To) then ! solid / vapour
+  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
      qsat = 100.0 * epsi * 10**(2.07023 - 0.00320991             &
           * T - 2484.896 / T + 3.56654 * alog10(T))
Index: /trunk/LMDZ.GENERIC/libf/phystd/watersat_grad.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/watersat_grad.F90	(revision 649)
+++ /trunk/LMDZ.GENERIC/libf/phystd/watersat_grad.F90	(revision 650)
@@ -1,5 +1,5 @@
 subroutine watersat_grad(T,qsat,dqsat)
 
-  use watercommon_h, only: T_coup, RLVTT, RCPD
+  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD,T_coup
   implicit none
 
@@ -7,6 +7,6 @@
 !     Purpose
 !     -------
-!     Compute the water mass mixing ratio at saturation (kg/kg)
-!     for a given pressure (Pa) and temperature (K)
+!     Compute the L/cp*d (q_sat)/d T
+!     for a given temperature (K)
 !
 !     Authors
@@ -22,5 +22,6 @@
   real dqsat
 
-  if (T.lt.T_coup) then ! solid / vapour
+!  if (T.lt.T_coup) then ! solid / vapour !why use T_coup?????????? JL12
+  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
      dqsat = RLVTT/RCPD*qsat*(3.56654/T             &
           +2484.896*LOG(10.)/T**2                   &
