Index: /branches/ORCHIDEE_2_2/ORCHIDEE/orchidee.default
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/orchidee.default	(revision 7198)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/orchidee.default	(revision 7199)
@@ -1406,11 +1406,5 @@
 
 # DRY_SOIL_HEAT_COND ([W.m^{-2}.K^{-1}] ) :  Dry soil Thermal Conductivity of soils 	{OK_SECHIBA}
-DRY_SOIL_HEAT_COND =  0.40  
-
-# WET_SOIL_HEAT_CAPACITY ([J.m^{-3}.K^{-1}]) :  Wet soil Heat capacity of soils  	{OK_SECHIBA}
-WET_SOIL_HEAT_CAPACITY =  3.03e+6 
-
-# WET_SOIL_HEAT_COND ([W.m^{-2}.K^{-1}]) :  Wet soil Thermal Conductivity of soils 	{OK_SECHIBA }
-WET_SOIL_HEAT_COND =  1.89  
+DRY_SOIL_HEAT_COND =  0.40 
 
 # SNOW_HEAT_COND ([W.m^{-2}.K^{-1}]) :  Thermal Conductivity of snow 	{OK_SECHIBA  }
@@ -1441,7 +1435,4 @@
 OK_FREEZE_THAW_LATENT_HEAT =  FALSE  
 
-# POROS ([-] ) :  Soil porosity  	{OK_SECHIBA}
-POROS =  0.41 
-
 # fr_dT ([K] ) :  Freezing window     	{OK_SECHIBA}
 fr_dT =  2.0 
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil.f90	(revision 7198)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil.f90	(revision 7199)
@@ -105,38 +105,4 @@
 
 
-       !Config Key   = WET_SOIL_HEAT_CAPACITY
-       !Config Desc  = Wet soil Heat capacity of soils 
-       !Config If    = OK_SECHIBA
-       !Config Def   = 3.03e+6
-       !Config Help  = 
-       !Config Units = [J.m^{-3}.K^{-1}]
-       CALL getin_p("WET_SOIL_HEAT_CAPACITY",so_capa_wet)
-
-       !! Check parameter value (correct range)
-       IF ( so_capa_wet <= zero ) THEN
-          CALL ipslerr_p(error_level, "config_soil_parameters.", &
-               &     "Wrong parameter value for WET_SOIL_HEAT_CAPACITY.", &
-               &     "This parameter should be positive. ", &
-               &     "Please, check parameter value in run.def. ")
-       END IF
-
-
-       !Config Key   = WET_SOIL_HEAT_COND
-       !Config Desc  = Wet soil Thermal Conductivity of soils
-       !Config If    = OK_SECHIBA 
-       !Config Def   = 1.89 
-       !Config Help  = 
-       !Config Units = [W.m^{-2}.K^{-1}]
-       CALL getin_p("WET_SOIL_HEAT_COND",so_cond_wet)
-
-       !! Check parameter value (correct range)
-       IF ( so_cond_wet <= zero ) THEN
-          CALL ipslerr_p(error_level, "config_soil_parameters.", &
-               &     "Wrong parameter value for WET_SOIL_HEAT_COND.", &
-               &     "This parameter should be positive. ", &
-               &     "Please, check parameter value in run.def. ")
-       END IF
-
-
        !Config Key   = SNOW_HEAT_COND
        !Config Desc  = Thermal Conductivity of snow
@@ -290,14 +256,4 @@
     ok_freeze_thaw_latent_heat = .FALSE.
     CALL getin_p ('OK_FREEZE_THAW_LATENT_HEAT',ok_freeze_thaw_latent_heat)
-
-
-    !Config Key = POROS
-    !Config Desc = Soil porosity 
-    !Config If = OK_SECHIBA
-    !Config Def = 0.41
-    !Config Help = From USDA classification, mean value
-    !Config Units = [-] 
-    poros=0.41
-    CALL getin_p('POROS',poros)
 
 
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil_var.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil_var.f90	(revision 7198)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil_var.f90	(revision 7199)
@@ -80,10 +80,4 @@
                                                         !! @tex $(W.m^{-2}.K^{-1})$ @endtex
 !$OMP THREADPRIVATE(so_cond_dry)
-  REAL(r_std), SAVE :: so_capa_wet = 3.03e+6            !! Wet soil Heat capacity of soils 
-                                                        !! @tex $(J.m^{-3}.K^{-1})$ @endtex
-!$OMP THREADPRIVATE(so_capa_wet)
-  REAL(r_std), SAVE :: so_cond_wet = 1.89               !! Wet soil Thermal Conductivity of soils 
-                                                        !! @tex $(W.m^{-2}.K^{-1})$ @endtex 
-!$OMP THREADPRIVATE(so_cond_wet)
   REAL(r_std), SAVE :: sn_cond = 0.3                    !! Thermal Conductivity of snow 
                                                         !! @tex $(W.m^{-2}.K^{-1})$ @endtex  
@@ -163,7 +157,4 @@
  & (/ 0.1_r_std, 0.1_r_std, 0.1_r_std /)                               !!  @tex $(m^{3} m^{-3})$ @endtex
 
-  REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: SMCMAX_fao = &          !! porosity
- & (/ 0.41_r_std, 0.43_r_std, 0.41_r_std /)                            !! & (/ 0.434_r_std, 0.439_r_std, 0.465_r_std /) !!noah lsm
-
   REAL(r_std),PARAMETER,DIMENSION(nscm_fao) :: QZ_fao = &              !! QUARTZ CONTENT (SOIL TYPE DEPENDENT)
  & (/ 0.60_r_std, 0.40_r_std, 0.35_r_std /)                            !! Peters et al [1998]
@@ -247,10 +238,4 @@
  &    0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std, &
  &    0.1_r_std /) ! oxisols
-
-  REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: SMCMAX_usda = &          !! porosity for thermosoil
- & (/ 0.43_r_std, 0.41_r_std, 0.41_r_std, 0.45_r_std, &
- &    0.46_r_std, 0.43_r_std, 0.39_r_std, 0.43_r_std, &
- &    0.41_r_std, 0.38_r_std, 0.36_r_std, 0.38_r_std, &
- &    0.503_r_std  /) ! oxisols
  
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: QZ_usda = &              !! QUARTZ CONTENT (SOIL TYPE DEPENDENT)
@@ -283,6 +268,4 @@
   LOGICAL, SAVE        :: read_reftemp                !! Flag to initialize soil temperature using climatological temperature
 !$OMP THREADPRIVATE(read_reftemp)
-  REAL(r_std), SAVE    :: poros                       !! Soil porosity (from USDA classification, mean value)(-)
-!$OMP THREADPRIVATE(poros)
   REAL(r_std), SAVE    :: fr_dT                       !! Freezing window (K)
 !$OMP THREADPRIVATE(fr_dT)
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_var.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_var.f90	(revision 7198)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_var.f90	(revision 7199)
@@ -380,6 +380,4 @@
   !! Heat capacity
   REAL(r_std), PARAMETER :: capa_ice = 2.228*1.E3       !! Heat capacity of ice (J/kg/K)
-  REAL(r_std), SAVE      :: so_capa_ice                 !! Heat capacity of saturated frozen soil (J/K/m3)
-!$OMP THREADPRIVATE(so_capa_ice)
   REAL(r_std), PARAMETER :: rho_water = 1000.           !! Density of water (kg/m3)
   REAL(r_std), PARAMETER :: rho_ice = 920.              !! Density of ice (kg/m3)
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 7198)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 7199)
@@ -145,7 +145,8 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz              !! Frozen fraction of the soil on hydrological levels (-)
 !$OMP THREADPRIVATE(profil_froz)
-  REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat                  !! Accumulated latent heat for the whole soil (J)
+  REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat                  !! Latent heat released or consumed in the freezing/thawing processes summed vertically
+                                                                              !! for the whole soil (J/m2) and on the whole simulation to check/correct energy conservation
 !$OMP THREADPRIVATE(e_soil_lat)
-  REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp               !! Additional heat capacity due to soil freezing for each soil layer (J/K)
+  REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp               !! Additional surfacic heat capacity due to soil freezing for each soil layer (J/K/m2)
 !$OMP THREADPRIVATE(pcappa_supp)    
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: dz5                      !! Used for numerical calculation [-]
@@ -153,10 +154,10 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcs                      !! Saturation humidity [m3/m3]
 !$OMP THREADPRIVATE(mcs)
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: SMCMAX                   !! Soil porosity [m3/m3]
-!$OMP THREADPRIVATE(SMCMAX)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: QZ                       !! quartz content [-]
 !$OMP THREADPRIVATE(QZ)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: so_capa_dry_ns           !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1} 
 !$OMP THREADPRIVATE(so_capa_dry_ns)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: so_capa_ice              !! Heat capacity of saturated frozen soil (J/K/m3)
+!$OMP THREADPRIVATE(so_capa_ice)   
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layt                  !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface
 !$OMP THREADPRIVATE(mc_layt)
@@ -271,5 +272,5 @@
 
     !! 0.4 Local variables
-    INTEGER(i_std)                                        :: ier, i, jg
+    INTEGER(i_std)                                        :: ier, i, jg, jsc
     LOGICAL                                               :: calculate_coef   !! Local flag to initialize variables by call to thermosoil_coef   
 !_ ================================================================================================================================
@@ -372,7 +373,4 @@
     IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcs','','')
 
-    ALLOCATE (SMCMAX(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of SMCMAX','','')
-
     ALLOCATE (QZ(nscm),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','')
@@ -380,4 +378,7 @@
     ALLOCATE (so_capa_dry_ns(nscm),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','')
+    
+    ALLOCATE (so_capa_ice(nscm),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_ice','','')
 
 
@@ -385,11 +386,9 @@
     SELECTCASE (nscm)
     CASE (3)
-       SMCMAX(:) = SMCMAX_fao(:)
        QZ(:) = QZ_fao(:)
        so_capa_dry_ns(:) = so_capa_dry_ns_fao(:)
        mcs(:) = mcs_fao(:)
     CASE (13) !Salma changed from 12 to 13 for the new class Oxisols
-       SMCMAX(:) = SMCMAX_usda(:)
-      QZ(:) = QZ_usda(:)
+       QZ(:) = QZ_usda(:)
        so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
        mcs(:) = mcs_usda(:)
@@ -505,6 +504,8 @@
 
     ! Calculate so_capa_ice
-    so_capa_ice = so_capa_dry + poros*capa_ice*rho_ice
-    IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice=', so_capa_ice,' using poros=',poros,' and capa_ice=',capa_ice
+    DO jsc = 1, nscm  
+       so_capa_ice(jsc) = so_capa_dry + mcs(jsc)*capa_ice*rho_ice 
+    END DO
+    IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice(:)=', so_capa_ice(:),' and capa_ice=',capa_ice
     
     ! Computing some usefull constants for the numerical scheme
@@ -1257,5 +1258,5 @@
 !================================================================================================================================
 
-  SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, smcmax, sh2o, cnd)
+  SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, sh2o, cnd)
 
     !! 0. Variables and parameter declaration
@@ -1265,6 +1266,5 @@
     INTEGER(i_std), DIMENSION (kjpindex), INTENT (in)          :: njsc          !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
     REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN)        :: smc           !! Volumetric Soil Moisture Content (m3/m3)
-    REAL(r_std), DIMENSION (nscm), INTENT(IN)                  :: qz            !! Quartz Content (Soil Type Dependent) (0-1)               
-    REAL(r_std), DIMENSION (nscm), INTENT(IN)                  :: smcmax        !! Soil Porosity (0-1)
+    REAL(r_std), DIMENSION (nscm), INTENT(IN)                  :: qz            !! Quartz Content (Soil Type Dependent) (0-1)
     REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN)        :: sh2o          !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o
     
@@ -1297,5 +1297,5 @@
 
       !! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
-      gammd = (1. - smcmax(jst))*2700.
+      gammd = (1. - mcs(jst))*2700.
       thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd)
 
@@ -1312,11 +1312,11 @@
       DO jg = 1,ngrnd      
         !! 1.4. saturation ratio
-        satratio(ji,jg) = smc(ji,jg) / smcmax(jst)
+        satratio(ji,jg) = smc(ji,jg) / mcs(jst)
     
         !! 1.5. Saturated Thermal Conductivity (thksat)
         IF ( smc(ji,jg) > min_sechiba ) THEN
            xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg)   ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen))
-           xu(ji,jg) = xunfroz(ji,jg) * smcmax(jst)  ! Unfrozen volume for saturation (porosity*xunfroz)
-           thksat(ji,jg) = thks ** (1. - smcmax(jst))* THKICE ** (smcmax(jst) - xu(ji,jg))* THKW ** (xu(ji,jg))
+           xu(ji,jg) = xunfroz(ji,jg) * mcs(jst)  ! Unfrozen volume for saturation (porosity*xunfroz)
+           thksat(ji,jg) = thks ** (1. - mcs(jst))* THKICE ** (mcs(jst) - xu(ji,jg))* THKW ** (xu(ji,jg))
         ELSE
            ! this value will not be used since ake=0 for this case
@@ -1622,6 +1622,5 @@
     REAL						:: xx         !! Unfrozen fraction of the soil
     REAL(r_std), DIMENSION(kjpindex)             	:: snow_h
-    REAL(r_std), DIMENSION(kjpindex,ngrnd) 		:: zx1, zx2    
-    REAL			   			:: cap_iw     !! Heat capacity of ice/water mixture 
+    REAL(r_std), DIMENSION(kjpindex,ngrnd) 		:: zx1, zx2
     INTEGER						:: ji,jg
     INTEGER                                             :: jst
@@ -1652,5 +1651,5 @@
 	    profil_froz(ji,jg) = 1.
  	    pcappa_supp(ji,jg)= 0.
-            pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + so_capa_ice * tmc_layt(ji,jg) / mille / dlt(jg)
+            pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + so_capa_ice(jst) * tmc_layt(ji,jg) / mille / dlt(jg)
             rho_tot = rho_soil * (1-mcs(jst)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg) 
             pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
@@ -1668,15 +1667,13 @@
            profil_froz(ji,jg) = (1. - xx)
 
-	   ! net heat capacity of the ice/water mixture
-     	   cap_iw = xx * so_capa_wet + (1.-xx) * so_capa_ice
     	   IF (ok_freeze_thaw_latent_heat) THEN
               pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + &
                 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
-                so_capa_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &
-                shum_ngrnd_perma(ji,jg)*poros*lhf*rho_water/fr_dT
+                so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &
+                shum_ngrnd_perma(ji,jg)*mcs(jst)*lhf*rho_water/fr_dT
 	   ELSE
               pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + &
                 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
-                so_capa_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
+                so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
            ENDIF
 
@@ -1686,5 +1683,5 @@
            pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
 
- 	   pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*poros*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg)
+ 	   pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*mcs(jst)*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg)
 
          ENDIF
@@ -1713,7 +1710,7 @@
     !
     IF (ok_freeze_thaw_latent_heat) THEN
-    	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX, mcl_layt*(1-profil_froz), pkappa)
+    	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa)
     ELSE
-    	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX, mcl_layt, pkappa)
+    	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa)
     ENDIF
 
@@ -1773,5 +1770,5 @@
       ENDDO
 
-      CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, SMCMAX, mcl_layt, pkappa)
+      CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa)
 
       IF (brk_flag == 1) THEN
