Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 7199)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 7206)
@@ -152,6 +152,4 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: dz5                      !! Used for numerical calculation [-]
 !$OMP THREADPRIVATE(dz5)
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcs                      !! Saturation humidity [m3/m3]
-!$OMP THREADPRIVATE(mcs)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: QZ                       !! quartz content [-]
 !$OMP THREADPRIVATE(QZ)
@@ -231,5 +229,5 @@
   !! \n
   !_ ==============================================================================================================================
-  SUBROUTINE thermosoil_initialize (kjit,          kjpindex,   rest_id,                   &
+  SUBROUTINE thermosoil_initialize (kjit,          kjpindex,   rest_id,          mcs,     &
                                     temp_sol_new,  snow,       shumdiag_perma,            &
                                     soilcap,       soilflx,    stempdiag,                 &
@@ -244,4 +242,5 @@
     INTEGER(i_std), INTENT(in)                            :: kjpindex         !! Domain size (unitless)
     INTEGER(i_std),INTENT (in)                            :: rest_id          !! Restart file identifier (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcs              !! Saturated moisture content (m3/m3)
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: temp_sol_new     !! Surface temperature at the present time-step,
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: snow             !! Snow mass (kg)
@@ -370,7 +369,4 @@
     IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt','','')
 
-    ALLOCATE (mcs(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcs','','')
-
     ALLOCATE (QZ(nscm),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','')
@@ -379,5 +375,5 @@
     IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','')
     
-    ALLOCATE (so_capa_ice(nscm),stat=ier)
+    ALLOCATE (so_capa_ice(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_ice','','')
 
@@ -388,9 +384,7 @@
        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
        QZ(:) = QZ_usda(:)
        so_capa_dry_ns(:) = so_capa_dry_ns_usda(:)
-       mcs(:) = mcs_usda(:)
     CASE DEFAULT
        WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map'
@@ -504,7 +498,5 @@
 
     ! Calculate so_capa_ice
-    DO jsc = 1, nscm  
-       so_capa_ice(jsc) = so_capa_dry + mcs(jsc)*capa_ice*rho_ice 
-    END DO
+    so_capa_ice(:) = so_capa_dry + mcs(:)*capa_ice*rho_ice 
     IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice(:)=', so_capa_ice(:),' and capa_ice=',capa_ice
     
@@ -530,5 +522,5 @@
        CALL thermosoil_coef (&
             kjpindex,      temp_sol_new,    snow,           njsc, &
-            frac_snow_veg, frac_snow_nobio, totfrac_nobio,        &
+            mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio,        &
             snowdz,        snowrho,         snowtemp,       pb,   &
             ptn,                                                  &
@@ -582,5 +574,5 @@
 
   SUBROUTINE thermosoil_main (kjit, kjpindex, &
-       index, indexgrnd, &
+       index, indexgrnd, mcs, &
        temp_sol_new, snow, soilcap, soilflx, &
        shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
@@ -601,4 +593,5 @@
     INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd        !! Indeces of the points on the 3D map (vertical 
                                                                               !! dimension towards the ground) (unitless)
+    REAL(i_std),DIMENSION (kjpindex), INTENT (in)         :: mcs              !! Saturated moisture content (m3/m3)
     REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: temp_sol_new     !! Surface temperature at the present time-step,
                                                                               !! Ts @tex ($K$) @endtex
@@ -756,5 +749,5 @@
     CALL thermosoil_coef (&
          kjpindex,      temp_sol_new,    snow,         njsc, &
-         frac_snow_veg, frac_snow_nobio, totfrac_nobio,      &
+         mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio,      &
          snowdz,        snowrho,         snowtemp,     pb,   &
          ptn,                                                &
@@ -918,5 +911,5 @@
 
   SUBROUTINE thermosoil_coef (kjpindex,      temp_sol_new,    snow,           njsc, &
-                              frac_snow_veg, frac_snow_nobio, totfrac_nobio,        &
+                              mcs,           frac_snow_veg,   frac_snow_nobio,totfrac_nobio, &
                               snowdz,        snowrho,         snowtemp,       pb,   &
                               ptn,                                                  &
@@ -932,4 +925,5 @@
     REAL(r_std), DIMENSION (kjpindex), INTENT (in)         :: snow         !! snow mass @tex ($Kg$) @endtex
     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), INTENT (in)         :: mcs          !! Saturated moisture content (m3/m3)
     REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: frac_snow_veg   !! Snow cover fraction on vegeted area
     REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in)    :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
@@ -994,8 +988,8 @@
     ! Computation of the soil thermal properties; snow properties are also accounted for
     IF (ok_freeze_thermix) THEN
-       CALL thermosoil_getdiff( kjpindex, snow, ptn, njsc, snowrho, snowtemp, pb )
+       CALL thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb )
     ELSE
        ! Special case without soil freezing
-       CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )
+       CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb )
     ENDIF
 
@@ -1258,5 +1252,5 @@
 !================================================================================================================================
 
-  SUBROUTINE thermosoil_cond (kjpindex, njsc, smc, qz, sh2o, cnd)
+  SUBROUTINE thermosoil_cond (kjpindex, njsc, mcs, smc, qz, sh2o, cnd)
 
     !! 0. Variables and parameter declaration
@@ -1265,4 +1259,5 @@
     INTEGER(i_std), INTENT(in)                                 :: kjpindex      !! Domain size (unitless)
     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), INTENT(IN)              :: mcs           !! Saturated moisture content (m3/m3)
     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)
@@ -1297,5 +1292,5 @@
 
       !! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1)
-      gammd = (1. - mcs(jst))*2700.
+      gammd = (1. - mcs(ji))*2700.
       thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd)
 
@@ -1312,11 +1307,11 @@
       DO jg = 1,ngrnd      
         !! 1.4. saturation ratio
-        satratio(ji,jg) = smc(ji,jg) / mcs(jst)
+        satratio(ji,jg) = smc(ji,jg) / mcs(ji)
     
         !! 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) * 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))
+           xu(ji,jg) = xunfroz(ji,jg) * mcs(ji)  ! Unfrozen volume for saturation (porosity*xunfroz)
+           thksat(ji,jg) = thks ** (1. - mcs(ji))* THKICE ** (mcs(ji) - xu(ji,jg))* THKW ** (xu(ji,jg))
         ELSE
            ! this value will not be used since ake=0 for this case
@@ -1606,5 +1601,5 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, njsc, snowrho, snowtemp, pb )
+  SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb )
 
    !! 0. Variables and parameter declaration
@@ -1614,7 +1609,8 @@
     REAL(r_std),DIMENSION(kjpindex),INTENT (in)	        :: snow       !! Snow mass
     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), INTENT (in)       :: mcs        !! Saturated moisture content (m3/m3)
     REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho    !! Snow density
     REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp   !! Snow temperature (K)
-    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: pb         !! Surface presure (hPa)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: pb         !! Surface pressure (hPa)
     REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in)	:: ptn        !! Soil temperature profile
 
@@ -1640,5 +1636,5 @@
       DO jg = 1, ngrnd
          jst = njsc(ji)
-         pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
+         pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
          !
          ! 2. Calculate volumetric heat capacity with allowance for permafrost
@@ -1651,6 +1647,6 @@
 	    profil_froz(ji,jg) = 1.
  	    pcappa_supp(ji,jg)= 0.
-            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(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + so_capa_ice(ji) * tmc_layt(ji,jg) / mille / dlt(jg)
+            rho_tot = rho_soil * (1-mcs(ji)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg) 
             pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
 
@@ -1660,5 +1656,5 @@
 	    profil_froz(ji,jg) = 0.
  	    pcappa_supp(ji,jg)= 0.
-            rho_tot = rho_soil * (1-mcs(jst)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)
+            rho_tot = rho_soil * (1-mcs(ji)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg)
             pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
      	 ELSE
@@ -1668,20 +1664,20 @@
 
     	   IF (ok_freeze_thaw_latent_heat) THEN
-              pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + &
+              pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + &
                 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
-                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
+                so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + &
+                shum_ngrnd_perma(ji,jg)*mcs(ji)*lhf*rho_water/fr_dT
 	   ELSE
-              pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + &
+              pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + &
                 water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
-                so_capa_ice(jst) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
+                so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
            ENDIF
 
-           rho_tot =  rho_soil* (1-mcs(jst)) + &
+           rho_tot =  rho_soil* (1-mcs(ji)) + &
                 rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + &
                 rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx)
            pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot
 
- 	   pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*mcs(jst)*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg)
+ 	   pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*mcs(ji)*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg)
 
          ENDIF
@@ -1710,7 +1706,7 @@
     !
     IF (ok_freeze_thaw_latent_heat) THEN
-    	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa)
+    	CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa)
     ELSE
-    	CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa)
+    	CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa)
     ENDIF
 
@@ -1744,5 +1740,5 @@
 !_ ================================================================================================================================
 
-    SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, njsc, snowrho, snowtemp, pb )
+    SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb )
 
    !! 0. Variables and parameter declaration
@@ -1751,7 +1747,8 @@
       INTEGER(i_std), INTENT(in) :: kjpindex
       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), INTENT (in)       :: mcs      !! Saturated moisture content (m3/m3)
       REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho  !! Snow density
       REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
-      REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: pb       !! Surface presure (hPa)
+      REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: pb       !! Surface pressure (hPa)
 
 
@@ -1764,5 +1761,5 @@
          DO ji = 1,kjpindex
             jst = njsc(ji)
-            pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(jst)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
+            pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg)
             pcapa(ji,jg) = pcapa_tmp(ji, jg)
             pcapa_en(ji,jg) = pcapa_tmp(ji, jg)
@@ -1770,5 +1767,5 @@
       ENDDO
 
-      CALL thermosoil_cond (kjpindex, njsc, mc_layt, QZ, mcl_layt, pkappa)
+      CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa)
 
       IF (brk_flag == 1) THEN
