diff --git a/src/SURFEX/av_pgd_param.F90 b/src/SURFEX/av_pgd_param.F90
index 6e7e253844610c73d9da3a439fadcd7511c7de4e..462ee4c915bd3cd265e044dd947fb6317f573945 100644
--- a/src/SURFEX/av_pgd_param.F90
+++ b/src/SURFEX/av_pgd_param.F90
@@ -108,7 +108,7 @@ REAL, DIMENSION(SIZE(PFIELD,1))   :: ZSUM_WEIGHT_PATCH
 REAL, DIMENSION(SIZE(PFIELD,1))   :: ZWORK
 REAL, DIMENSION(SIZE(PFIELD,1))   :: ZDZ
 !
-REAL, DIMENSION(31) :: ZCOUNT
+REAL, DIMENSION(0:31) :: ZCOUNT
 INTEGER, DIMENSION(SIZE(PFIELD,1))  :: NMASK
 INTEGER ::  PATCH_LIST(NVEGTYPE)
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -169,7 +169,7 @@ DO JV=1,NVEGTYPE
     ELSEIF (HSFTYPE=='TRE'.OR.HSFTYPE=='GRT') THEN
       IF (JV==NVT_TEBD.OR.JV==NVT_BONE.OR.JV==NVT_TRBE.OR.JV==NVT_TRBD.OR.&
           JV==NVT_TEBE.OR.JV==NVT_TENE.OR.JV==NVT_BOBD.OR.JV==NVT_BOND.OR.&
-          JV==NVT_SHRB) ZWEIGHT(JI,JV) = PVEGTYPE(JI,JV)
+          JV==NVT_SHRB) ZWEIGHT(JI,JV) = PVEGTYPE(IMASK,JV)
     ELSE
       CALL ABOR1_SFX('AV_PGD_PARAM_1D: WEIGHTING FUNCTION FOR VEGTYPE NOT ALLOWED')
     ENDIF
@@ -262,7 +262,7 @@ SELECT CASE (HATYPE)
       IF (ALL(ZCOUNT(:)==0.)) THEN
         ZWORK(JJ) = NUNDEF
       ELSE
-        ZWORK(JJ) = FLOAT(MAXLOC(ZCOUNT,1))
+        ZWORK(JJ) = FLOAT(MAXLOC(ZCOUNT,1)-1)
       ENDIF
     END DO
 !
diff --git a/src/SURFEX/diag_isba_initn.F90 b/src/SURFEX/diag_isba_initn.F90
index c4716f23d7738a8f7d8adeaf35c8010f24e2e0a9..50404eda5bc44c819ef61a126b87400167352a4c 100644
--- a/src/SURFEX/diag_isba_initn.F90
+++ b/src/SURFEX/diag_isba_initn.F90
@@ -124,7 +124,7 @@ TYPE(DIAG_EVAP_ISBA_t) :: YDE
 TYPE(DIAG_MISC_ISBA_t) :: YDM
 TYPE(ISBA_P_t), POINTER :: PK
 REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK
-LOGICAL :: GCUMUL, GDIM
+LOGICAL :: GCUMUL, GDIM, GDIM2
 INTEGER :: JP
 INTEGER           :: IVERSION, IBUG
 INTEGER           :: IRESP          ! IRESP  : return-code if a problem appears
@@ -448,7 +448,9 @@ IF ( GCUMUL ) THEN
     CALL READ_SURF(HPROGRAM,'VERSION',IVERSION,IRESP)
     CALL READ_SURF(HPROGRAM,'BUG ',IBUG,IRESP)
     !
-    GDIM = (IVERSION>8 .OR. IVERSION==8 .AND. IBUG>0)
+    GDIM2 = (IVERSION>8 .OR. IVERSION==8 .AND. IBUG>0)
+    GDIM = GDIM2
+    IF (GDIM2) CALL READ_SURF(HPROGRAM,'SPLIT_PATCH',GDIM,IRESP)    
     !
 #ifdef SFX_OL
     IF(DE%LWATER_BUDGET .OR. (LRESTART .AND. .NOT.DGO%LRESET_BUDGETC))THEN
diff --git a/src/SURFEX/e_budget_meb.F90 b/src/SURFEX/e_budget_meb.F90
index ffb942c2a7a178f6239b1652a877e589d94ac0c4..092705bc67ff92e59b850704b546e11aebd4d4e3 100644
--- a/src/SURFEX/e_budget_meb.F90
+++ b/src/SURFEX/e_budget_meb.F90
@@ -80,7 +80,7 @@
 !!    -------------
 !!      Original    22/01/11 
 !!                  10/10/14 (A. Boone) Removed understory vegetation
-!!
+!!                  13/09/18 (A. Boone) Add litter layer option to test-Tg computation
 !-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -92,14 +92,14 @@ USE MODD_DIAG_n, ONLY : DIAG_t
 USE MODD_DIAG_EVAP_ISBA_n, ONLY : DIAG_EVAP_ISBA_t
 USE MODD_DIAG_MISC_ISBA_n, ONLY : DIAG_MISC_ISBA_t
 !
-USE MODD_CSTS,                    ONLY : XLVTT, XLSTT, XTT, XCPD, XCPV, XCL, &
+USE MODD_CSTS,                    ONLY : XLVTT, XLSTT, XTT, XCPD, XCPV, XCL,  &
                                          XDAY, XPI, XLMTT, XRHOLW
 USE MODD_SURF_ATM,                ONLY : LCPL_ARP
 USE MODD_SURF_PAR,                ONLY : XUNDEF
 USE MODD_SNOW_METAMO,             ONLY : XSNOWDZMIN
-!
+
 USE MODE_THERMOS
-USE MODE_MEB,                     ONLY : SFC_HEATCAP_VEG
+USE MODE_MEB,                     ONLY : SFC_HEATCAP_VEG, MEBLITTER_THRM
 USE MODE_SNOW3L,                  ONLY : SNOW3LHOLD
 !
 USE MODI_TRIDIAG_GROUND_RM_COEFS
@@ -308,7 +308,7 @@ REAL, DIMENSION(SIZE(DMK%XSNOWDZ,1),SIZE(DMK%XSNOWDZ,2)):: ZSNOW_COEF_A, ZSNOW_C
 !
 REAL, DIMENSION(SIZE(DMK%XSNOWDZ,1),SIZE(DMK%XSNOWDZ,2)):: ZWHOLDMAX
 !
-REAL, DIMENSION(SIZE(PD_G,1),SIZE(PD_G,2)+SIZE(DMK%XSNOWDZ,2)) :: ZD, ZT, ZHCAPZ, ZCONDZ,         &
+REAL, DIMENSION(SIZE(PD_G,1),SIZE(PD_G,2)+SIZE(DMK%XSNOWDZ,2)+1) :: ZD, ZT, ZHCAPZ, ZCONDZ,       &
                                              ZCOEF_A, ZCOEF_B, ZSOURCE
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -343,8 +343,15 @@ ZSOIL_COEF_A(:,:) = 0.0
 ZSOIL_COEF_B(:,:) = 0.0
 ZSNOW_COEF_A(:,:) = 0.0
 ZSNOW_COEF_B(:,:) = 0.0
-
-!-------------------------------------------------------------------------------
+!
+ZD(:,:)           = 0.0
+ZT(:,:)           = 0.0
+ZHCAPZ(:,:)       = 0.0
+ZCONDZ(:,:)       = 0.0
+ZSOURCE(:,:)      = 0.0
+ZCOEF_A(:,:)      = 0.0
+ZCOEF_B(:,:)      = 0.0
+!
 !
 !*       1.     Some variables/coefficients needed for coupling or solution
 !               -----------------------------------------------------------
@@ -429,7 +436,7 @@ IF(IO%CISBA == 'DIF')THEN
 ! quite robust. Note, this only corresponds to the snow-covered part of grid box,
 ! so it is more accurate as the snow fraction approaches unity.
 ! Starting from snowpack surface downward to base of ground:
-!
+!   
    JL                     = 1
    ZD(:,JL)               = DMK%XSNOWDZ(:,1)
    ZT(:,JL)               = ZTNO(:,1)
@@ -446,6 +453,13 @@ IF(IO%CISBA == 'DIF')THEN
          ZSOURCE(JJ,JL)   = DEK%XSWNET_NS(JJ)*(PTAU_N(JJ,JK-1)-PTAU_N(JJ,JK))
       ENDDO
    ENDDO
+   IF(IO%LMEB_LITTER)THEN
+      JL                  = JL + 1
+      ZD(:,JL)            = PEK%XGNDLITTER(:)
+      ZT(:,JL)            = PEK%XTL(:)
+      CALL MEBLITTER_THRM(PEK%XWRL,PEK%XWRLI,PEK%XGNDLITTER,ZHCAPZ(:,JL),ZCONDZ(:,JL)) 
+      ZSOURCE(:,JL)       = 0.
+   ENDIF
    JL                     = JL + 1
    ZD(:,JL)               = PD_G(:,1)
    ZT(:,JL)               = ZTGO(:,1)
@@ -465,8 +479,9 @@ IF(IO%CISBA == 'DIF')THEN
 !
 ! Get coefficients from upward sweep (starting from soil base to snow surface):
 !
-   CALL TRIDIAG_GROUND_RM_COEFS(PTSTEP, ZD, ZT, ZHCAPZ, ZCONDZ, ZSOURCE, &
-                                PTDEEP_A, KK%XTDEEP, ZTCONDA_DELZ_N, ZCOEF_A, ZCOEF_B)
+   CALL TRIDIAG_GROUND_RM_COEFS(PTSTEP, ZD(:,1:JL), ZT(:,1:JL), ZHCAPZ(:,1:JL), ZCONDZ(:,1:JL), &
+                                ZSOURCE(:,1:JL),PTDEEP_A, KK%XTDEEP, ZTCONDA_DELZ_N,            &
+                                ZCOEF_A(:,1:JL), ZCOEF_B(:,1:JL))
 !
    ZSNOW_COEF_A(:,2)  = ZCOEF_A(:,2)
    ZSNOW_COEF_B(:,2)  = ZCOEF_B(:,2)
diff --git a/src/SURFEX/ini_var_from_data.F90 b/src/SURFEX/ini_var_from_data.F90
index 755b57888ba78c96d104fa660cd1b0c80967befa..c57f2c1a8c7f641d817d671137831236a1761aba 100644
--- a/src/SURFEX/ini_var_from_data.F90
+++ b/src/SURFEX/ini_var_from_data.F90
@@ -195,7 +195,8 @@ IF (HFTYP(1)=='DIRTYP') THEN
 
 ELSE
 
-  IF (.NOT.ALL(LEN_TRIM(HFNAM(:))/=0) .AND. .NOT.ALL(LEN_TRIM(HFNAM(2:))==0)) THEN
+        
+  IF (.NOT.ALL(LEN_TRIM(HFNAM(:))/=0) .AND. COUNT(LEN_TRIM(HFNAM(:))/=0)>1) THEN
     DO JV=1,SIZE(PFIELD,2)
       IF (LEN_TRIM(HFNAM(JV))==0) THEN
         DO JV2=JV-1,1,-1
diff --git a/src/SURFEX/isba_meb.F90 b/src/SURFEX/isba_meb.F90
index 40b95a3cae7d8f72d21f3bfa0447227573e318e2..511af70b79a50f336f1a5eae08476de1065f668f 100644
--- a/src/SURFEX/isba_meb.F90
+++ b/src/SURFEX/isba_meb.F90
@@ -627,13 +627,15 @@ WHERE(PSW_RAD(:) > ZSWRAD_MIN) ! Sun is up...approx
 !
 ELSEWHERE
 
-! Sun is down:
-   
-   DK%XALBT(:)                = 1.0
-   ZSWUP(:)                   = PSW_RAD(:)
+! Sun is down: (below threshold)
+! radiation amounts quite small, so make a simple approximation here:   
+
+   DK%XALBT(:)                = ZALBV(:)
+   ZSWUP(:)                   = DK%XALBT(:)*PSW_RAD(:)
+
    DEK%XSWDOWN_GN(:)          = 0.
    DEK%XSWNET_G(:)            = 0.
-   DEK%XSWNET_V(:)            = 0.
+   DEK%XSWNET_V(:)            = (1.-DK%XALBT(:))*PSW_RAD(:)
    DEK%XSWNET_N(:)            = 0.
    DEK%XSWNET_NS(:)           = 0.
    ZTAU_N(:,SIZE(PEK%TSNOW%WSNOW,2)) = 0.
@@ -918,7 +920,7 @@ ZVEGFACT(:) = ZSIGMA_F(:)*(1.0-PPALPHAN(:)*PEK%XPSN(:))
 ! snowpack and part falling onto snow-free understory.
 !
 !
- CALL HYDRO_VEG(IO%CRAIN, PTSTEP, KK%XMUF, ZRR, DEK%XLEV_CV, DEK%XLETR_CV,          &
+ CALL HYDRO_VEG(IO%CRAIN, PTSTEP, KK%XMUF, ZRR, DEK%XLEV_CV, DEK%XLETR_CV,        &
                 ZVEGFACT, ZPSNCV, PEK%XWR, ZWRMAX, ZRRSFC, DEK%XDRIP, DEK%XRRVEG, &
                 PK%XLVTT  )
 !
@@ -930,7 +932,7 @@ ZVEGFACT(:) = ZSIGMA_F(:)*(1.0-PPALPHAN(:)*PEK%XPSN(:))
  IF (PRESENT(PBLOWSNW_FLUX)) THEN
    CALL SNOW3L_ISBA(IO, G, PK, PEK, DK, DEK, DMK, OMEB, HIMPLICIT_WIND,       &
                   TPTIME, PTSTEP, PK%XVEGTYPE_PATCH,  ZTGL, ZCTSFC,         &
-                  ZSOILHCAPZ, ZSOILCONDZ(:,1), PPS, PTA, PSW_RAD, PQA,      &
+                  ZSOILHCAPZ, ZSOILCONDZ(:,1), PPS, PEK%XTC, DEK%XSWDOWN_GN, PEK%XQC,&
                   PVMOD, PLW_RAD, ZRRSFC, DEK%XSR_GN, PRHOA, ZUREF, PEXNS,  &
                   PEXNA, PDIRCOSZW, ZZREF, ZALBG, ZD_G, ZDZG, PPEW_A_COEF,  &
                   PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF,       &
@@ -941,7 +943,7 @@ ZVEGFACT(:) = ZSIGMA_F(:)*(1.0-PPALPHAN(:)*PEK%XPSN(:))
  ELSE
    CALL SNOW3L_ISBA(IO, G, PK, PEK, DK, DEK, DMK, OMEB, HIMPLICIT_WIND,       &
                   TPTIME, PTSTEP, PK%XVEGTYPE_PATCH,  ZTGL, ZCTSFC,         &
-                  ZSOILHCAPZ, ZSOILCONDZ(:,1), PPS, PTA, PSW_RAD, PQA,      &
+                  ZSOILHCAPZ, ZSOILCONDZ(:,1), PEK%XTC, DEK%XSWDOWN_GN, PEK%XQC, &
                   PVMOD, PLW_RAD, ZRRSFC, DEK%XSR_GN, PRHOA, ZUREF, PEXNS,  &
                   PEXNA, PDIRCOSZW, ZZREF, ZALBG, ZD_G, ZDZG, PPEW_A_COEF,  &
                   PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF,       &
@@ -1675,16 +1677,6 @@ INTEGER                               :: INJ, INL, JJ, JL
 !
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
-!*      0.3    declarations of local parameters
-!
-REAL, PARAMETER                       :: Z1 = 45.0       !litter bulk density (kg/m3)
-REAL, PARAMETER                       :: Z2 = 0.1        !coeff for litter conductivity (W/(mK))
-REAL, PARAMETER                       :: Z3 = 0.03       !coeff for litter conductivity
-REAL, PARAMETER                       :: Z4 = 0.95       !litter porosity       (m3/m3)
-REAL, PARAMETER                       :: Z5 = 0.12       !litter field capacity (m3/m3)
-!
-REAL, DIMENSION(SIZE(PEK%XWG,1))      :: ZWORK
-!
 !-------------------------------------------------------------------------------
 !
 IF (LHOOK) CALL DR_HOOK('ISBA_MEB:PREP_MEB_SOIL',0,ZHOOK_HANDLE)
@@ -1692,14 +1684,11 @@ IF (LHOOK) CALL DR_HOOK('ISBA_MEB:PREP_MEB_SOIL',0,ZHOOK_HANDLE)
 INJ  = SIZE(PK%XDG,1)
 INL  = SIZE(PK%XDG,2)
 !
-ZWORK(:) = 0.0
 IF(OMEB_LITTER)THEN
    PTGL(:,1)                  = PEK%XTL(:)
-   ZWORK(:)                   = PEK%XWRL(:)/(XRHOLW*PEK%XGNDLITTER(:))
-   PSOILHCAPL(:,1)            = XOMSPH*Z1 + (XCL*XRHOLW)*ZWORK(:) + (XCI*XRHOLI/XRHOLW)*PEK%XWRLI(:)/PEK%XGNDLITTER(:)
-   PSOILCONDL(:,1)            = Z2 + Z3 * ZWORK(:)
-   PWSATL(:,1)                = Z4
-   PWFCL(:,1)                 = Z5
+   CALL MEBLITTER_THRM(PEK%XWRL,PEK%XWRLI,PEK%XGNDLITTER,PSOILHCAPL(:,1),PSOILCONDL(:,1))
+   PWSATL(:,1)                = XLITTER_HYD_Z4
+   PWFCL(:,1)                 = XLITTER_HYD_Z5
    PD_GL(:,1)                 = PEK%XGNDLITTER(:)
    PDZGL(:,1)                 = PEK%XGNDLITTER(:)
    PCTSFC(:)                  = 1. / (PSOILHCAPL(:,1) * PEK%XGNDLITTER(:))
diff --git a/src/SURFEX/modd_meb_par.F90 b/src/SURFEX/modd_meb_par.F90
index 3a0909f7152da8b903b7a8202302002490bdd70d..17169f4c2afc11c73abb94da48ae185a344b3ac3 100644
--- a/src/SURFEX/modd_meb_par.F90
+++ b/src/SURFEX/modd_meb_par.F90
@@ -24,6 +24,8 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original       09/2013
+!!      13/09/18 (A. Boone) Added Litter thermal and hydrological parameters
+!!                          herein
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -56,6 +58,15 @@ REAL,    PARAMETER   :: XSW_WGHT_VIS = 0.48
 !
 REAL,    PARAMETER   :: XSW_WGHT_NIR = 0.52
 !
+! Litter thermal (THRM) and hydrological (HYD) properties
+! -------------------------------------------------------
+!
+REAL,    PARAMETER   :: XLITTER_THRM_Z1 = 45.00   !litter bulk density (kg/m3)
+REAL,    PARAMETER   :: XLITTER_THRM_Z2 =  0.10   !coeff for litter conductivity (W/(mK))
+REAL,    PARAMETER   :: XLITTER_THRM_Z3 =  0.03   !coeff for litter conductivity
+REAL,    PARAMETER   :: XLITTER_HYD_Z4  =  0.95   !litter porosity       (m3/m3)
+REAL,    PARAMETER   :: XLITTER_HYD_Z5  =  0.12   !litter field capacity (m3/m3)
+!
 !-------------------------------------------------------------------------------
 !
 END MODULE MODD_MEB_PAR
diff --git a/src/SURFEX/mode_meb.F90 b/src/SURFEX/mode_meb.F90
index 8665dd58adc64845058437da69ace5fc47b5db6f..68bb4fa2f29a69118ba992ad2c369f96e8ebba9e 100644
--- a/src/SURFEX/mode_meb.F90
+++ b/src/SURFEX/mode_meb.F90
@@ -32,6 +32,7 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original        18/01/11
+!!                      13/09/18 Added litter thermal computations here
 !-----------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -68,11 +69,224 @@ INTERFACE MEB_SHIELD_FACTOR
    MODULE PROCEDURE MEB_SHIELD_FACTOR_0D
 END INTERFACE
 !
+INTERFACE MEBLITTER_THRM
+   MODULE PROCEDURE MEBLITTER_THRM_2D
+   MODULE PROCEDURE MEBLITTER_THRM_1D
+   MODULE PROCEDURE MEBLITTER_THRM_0D
+END INTERFACE
+!
 !-------------------------------------------------------------------------------
 CONTAINS
 !
 !####################################################################
 !####################################################################
+!####################################################################
+      SUBROUTINE MEBLITTER_THRM_2D(PWRL,PWRLI,PGNDLITTER,PHCAP,PCOND) 
+!
+!!    PURPOSE
+!!    -------
+!     Calculation of litter thermal properties (based on Napoly et al., 2017)
+!
+!!    REFERENCE
+!!    ---------
+!!      Napoly et al., 2017: GMD
+!!
+!!    AUTHOR
+!!    ------
+!!!     A. Boone                 * Meteo-France/CNRS *       
+!!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original        13/09/18
+!
+!-----------------------------------------------------------------------------
+!
+USE MODD_MEB_PAR,                 ONLY : XLITTER_THRM_Z1, XLITTER_THRM_Z2,   &
+                                         XLITTER_THRM_Z3        
+USE MODD_ISBA_PAR,                ONLY : XOMSPH 
+USE MODD_CSTS,                    ONLY : XRHOLW, XRHOLI, XCI, XCL
+
+
+USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
+USE PARKIND1  ,ONLY : JPRB
+!
+IMPLICIT NONE
+!
+!*      0.1    declarations of arguments
+!
+REAL, DIMENSION(:,:), INTENT(IN)   :: PWRL, PWRLI, PGNDLITTER
+!                                     PWRL  = litter liquid water content (kg/m2)
+!                                     PWRLI = litter liquid water equivalent ice content (kg/m2)
+!                                     PGNDLITTER = litter layer thickness (m)
+!
+REAL, DIMENSION(:,:), INTENT(OUT)  :: PHCAP, PCOND
+!                                     PHCAP = litter heat capacity [J/(m3 K)]
+!                                     PCOND = litter thermal conductivity [W/(m K)]
+!
+!*      0.2    declaration of local variables
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+REAL, DIMENSION(SIZE(PGNDLITTER,1),SIZE(PGNDLITTER,2))  :: ZWORK
+!
+!-------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_MEB:MEBLITTER_THRM_2D',0,ZHOOK_HANDLE)
+!
+ZWORK(:,:)  = PWRL(:,:)/(XRHOLW*PGNDLITTER(:,:))
+!
+! Litter heat capacity:
+! 
+PHCAP(:,:)  = XOMSPH*XLITTER_THRM_Z1 + ( (XCL*XRHOLW)*ZWORK(:,:) +    &
+              (XCI*XRHOLI/XRHOLW)*PWRLI(:,:)/PGNDLITTER(:,:) )
+!
+! Litter thermal conductivity:
+!
+PCOND(:,:)  = XLITTER_THRM_Z2 + ( XLITTER_THRM_Z3 * ZWORK(:,:) )
+!
+IF (LHOOK) CALL DR_HOOK('MODE_MEB:MEBLITTER_THRM_2D',1,ZHOOK_HANDLE)
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MEBLITTER_THRM_2D
+!####################################################################
+!####################################################################
+!####################################################################
+      SUBROUTINE MEBLITTER_THRM_1D(PWRL,PWRLI,PGNDLITTER,PHCAP,PCOND) 
+!
+!!    PURPOSE
+!!    -------
+!     Calculation of litter thermal properties (based on Napoly et al., 2017)
+!
+!!    REFERENCE
+!!    ---------
+!!      Napoly et al., 2017: GMD
+!!
+!!    AUTHOR
+!!    ------
+!!!     A. Boone                 * Meteo-France/CNRS *       
+!!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original        13/09/18
+!
+!-----------------------------------------------------------------------------
+!
+USE MODD_MEB_PAR,                 ONLY : XLITTER_THRM_Z1, XLITTER_THRM_Z2,   &
+                                         XLITTER_THRM_Z3        
+USE MODD_ISBA_PAR,                ONLY : XOMSPH 
+USE MODD_CSTS,                    ONLY : XRHOLW, XRHOLI, XCI, XCL
+
+
+USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
+USE PARKIND1  ,ONLY : JPRB
+!
+IMPLICIT NONE
+!
+!*      0.1    declarations of arguments
+!
+REAL, DIMENSION(:), INTENT(IN)     :: PWRL, PWRLI, PGNDLITTER
+!                                     PWRL  = litter liquid water content (kg/m2)
+!                                     PWRLI = litter liquid water equivalent ice content (kg/m2)
+!                                     PGNDLITTER = litter layer thickness (m)
+!
+REAL, DIMENSION(:), INTENT(OUT)    :: PHCAP, PCOND
+!                                     PHCAP = litter heat capacity [J/(m3 K)]
+!                                     PCOND = litter thermal conductivity [W/(m K)]
+!
+!*      0.2    declaration of local variables
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+REAL, DIMENSION(SIZE(PGNDLITTER))  :: ZWORK
+!
+!-------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_MEB:MEBLITTER_THRM_1D',0,ZHOOK_HANDLE)
+!
+ZWORK(:)  = PWRL(:)/(XRHOLW*PGNDLITTER(:))
+!
+! Litter heat capacity:
+! 
+PHCAP(:)  = XOMSPH*XLITTER_THRM_Z1 + ( (XCL*XRHOLW)*ZWORK(:) +    &
+            (XCI*XRHOLI/XRHOLW)*PWRLI(:)/PGNDLITTER(:) )
+!
+! Litter thermal conductivity:
+!
+PCOND(:)  = XLITTER_THRM_Z2 + ( XLITTER_THRM_Z3 * ZWORK(:) )
+!
+IF (LHOOK) CALL DR_HOOK('MODE_MEB:MEBLITTER_THRM_1D',1,ZHOOK_HANDLE)
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MEBLITTER_THRM_1D
+!####################################################################
+!####################################################################
+!####################################################################
+      SUBROUTINE MEBLITTER_THRM_0D(PWRL,PWRLI,PGNDLITTER,PHCAP,PCOND) 
+!
+!!    PURPOSE
+!!    -------
+!     Calculation of litter thermal properties (based on Napoly et al., 2017)
+!
+!!    REFERENCE
+!!    ---------
+!!      Napoly et al., 2017: GMD
+!!
+!!    AUTHOR
+!!    ------
+!!!     A. Boone                 * Meteo-France/CNRS *       
+!!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original        13/09/18
+!
+!-----------------------------------------------------------------------------
+!
+USE MODD_MEB_PAR,                 ONLY : XLITTER_THRM_Z1, XLITTER_THRM_Z2,   &
+                                         XLITTER_THRM_Z3        
+USE MODD_ISBA_PAR,                ONLY : XOMSPH 
+USE MODD_CSTS,                    ONLY : XRHOLW, XRHOLI, XCI, XCL
+
+
+USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
+USE PARKIND1  ,ONLY : JPRB
+!
+IMPLICIT NONE
+!
+!*      0.1    declarations of arguments
+!
+REAL,               INTENT(IN)     :: PWRL, PWRLI, PGNDLITTER
+!                                     PWRL  = litter liquid water content (kg/m2)
+!                                     PWRLI = litter liquid water equivalent ice content (kg/m2)
+!                                     PGNDLITTER = litter layer thickness (m)
+!
+REAL,               INTENT(OUT)    :: PHCAP, PCOND
+!                                     PHCAP = litter heat capacity [J/(m3 K)]
+!                                     PCOND = litter thermal conductivity [W/(m K)]
+!
+!*      0.2    declaration of local variables
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+REAL                               :: ZWORK
+!
+!-------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_MEB:MEBLITTER_THRM_0D',0,ZHOOK_HANDLE)
+!
+ZWORK  = PWRL/(XRHOLW*PGNDLITTER)
+!
+! Litter heat capacity:
+! 
+PHCAP  = XOMSPH*XLITTER_THRM_Z1 +  (XCL*XRHOLW)*ZWORK +    &
+            (XCI*XRHOLI/XRHOLW)*PWRLI/PGNDLITTER 
+!
+! Litter thermal conductivity:
+!
+PCOND  = XLITTER_THRM_Z2 + XLITTER_THRM_Z3 * ZWORK 
+!
+IF (LHOOK) CALL DR_HOOK('MODE_MEB:MEBLITTER_THRM_0D',1,ZHOOK_HANDLE)
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE MEBLITTER_THRM_0D
+!####################################################################
+!####################################################################
 !####################################################################
       FUNCTION MEBPALPHAN_3D(PSNOWDEPTH,PH_VEG) RESULT(PPALPHAN)
 !
diff --git a/src/SURFEX/prep_snow_extern.F90 b/src/SURFEX/prep_snow_extern.F90
index 43f669649b1f9b1dc3daab7b746c7ca1e2183104..8b69ea59c204364ce632b6e9b90efe18dc85de10 100644
--- a/src/SURFEX/prep_snow_extern.F90
+++ b/src/SURFEX/prep_snow_extern.F90
@@ -151,6 +151,12 @@ CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'FULL  ')
 CALL READ_SURF(HFILEPGDTYPE,'VERSION',IVERSION_PGD,IRESP,HDIR='-')
 CALL READ_SURF(HFILEPGDTYPE,'BUG',IBUGFIX_PGD,IRESP,HDIR='-')
 GOLD_NAME=(IVERSION_PGD<7 .OR. (IVERSION_PGD==7 .AND. IBUGFIX_PGD<3))
+!
+!
+!-------------------------------------------------------------------------------------
+!
+!*      2.     Reading of grid
+!              ---------------
 !
  CALL PREP_GRID_EXTERN(GCP,HFILEPGDTYPE,KLUOUT,CINGRID_TYPE,CINTERP_TYPE,INI)
 !
@@ -160,6 +166,8 @@ CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
 !
 CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,YMASK)
 !
+IF (NRANK/=NPIO) INI = 0
+!
 ALLOCATE(ZMASK(INI))
 IF (IVERSION_PGD>=7) THEN 
   IF (YAREA(1:4)=='VEG ') THEN
@@ -199,13 +207,6 @@ END IF
 !
 !-------------------------------------------------------------------------------------
 !
-!*      2.     Reading of grid
-!              ---------------
-!
-IF (NRANK/=NPIO) INI = 0
-!
-!-------------------------------------------------------------------------------------
-!
 !*      4.     Reading of snow data
 !              ---------------------
 !
@@ -242,6 +243,7 @@ DO JP = 1,IPATCH
   IF (NRANK==NPIO) THEN
     !
     SELECT CASE (HSURF(1:3))
+
       CASE ('WWW')
         IF (OSNOW_IDEAL) THEN
           IF (JP<=1) ALLOCATE(PFIELD(INI,KLAYER,IPATCH))
@@ -260,125 +262,132 @@ DO JP = 1,IPATCH
   !*      6.     Albedo
   !              ------
   !
-    CASE ('ALB')
-      IF (JP<=1) ALLOCATE(PFIELD(INI,1,IPATCH))
-      PFIELD(:,1,JP) = TZSNOW%ALB(:)
+      CASE ('ALB')
+        IF (JP<=1) ALLOCATE(PFIELD(INI,1,IPATCH))
+        PFIELD(:,1,JP) = TZSNOW%ALB(:)
   !
   !-------------------------------------------------------------------------------------
   !
   !*      7.     Total depth to snow grid
   !              ------------------------
   !
-    CASE ('DEP')
-      IF (OSNOW_IDEAL) THEN    
-        IF (JP<=1) ALLOCATE(PFIELD(INI,KLAYER,IPATCH))  
-        PFIELD(:,:,JP) = TZSNOW%WSNOW(:,1:KLAYER)/TZSNOW%RHO(:,1:KLAYER)
-        WHERE(TZSNOW%WSNOW(:,1:KLAYER)==XUNDEF) PFIELD(:,:,JP)=XUNDEF
-      ELSE     
-        ALLOCATE(ZD(INI))
-        ZD(:) = 0.0
-        DO JL=1,TZSNOW%NLAYER
-          WHERE (TZSNOW%WSNOW(:,JL)/=XUNDEF)
-            ZD(:) = ZD(:) + TZSNOW%WSNOW(:,JL)/TZSNOW%RHO(:,JL)
-          ENDWHERE
-        END DO
-        IF (JP<=1) ALLOCATE(PFIELD(INI,1,IPATCH))
-        PFIELD(:,1,JP) = ZD(:)
-        DEALLOCATE(ZD)
-      ENDIF
+      CASE ('DEP')
+        IF (OSNOW_IDEAL) THEN    
+          IF (JP<=1) ALLOCATE(PFIELD(INI,KLAYER,IPATCH))  
+          PFIELD(:,:,JP) = TZSNOW%WSNOW(:,1:KLAYER)/TZSNOW%RHO(:,1:KLAYER)
+          WHERE(TZSNOW%WSNOW(:,1:KLAYER)==XUNDEF) PFIELD(:,:,JP)=XUNDEF
+        ELSE     
+          ALLOCATE(ZD(INI))
+          ZD(:) = 0.0
+          DO JL=1,TZSNOW%NLAYER
+            WHERE (TZSNOW%WSNOW(:,JL)/=XUNDEF)
+              ZD(:) = ZD(:) + TZSNOW%WSNOW(:,JL)/TZSNOW%RHO(:,JL)
+            ENDWHERE
+          END DO
+          IF (JP<=1) ALLOCATE(PFIELD(INI,1,IPATCH))
+          PFIELD(:,1,JP) = ZD(:)
+          DEALLOCATE(ZD)
+        ENDIF
   !
   !-------------------------------------------------------------------------------------
   !
   !*      8.     Density or heat profile
   !              -----------------------
   !
-    CASE ('RHO','HEA','SG1','SG2','HIS','AGE')
+      CASE ('RHO','HEA','SG1','SG2','HIS','AGE')
   !
-      SELECT CASE (TZSNOW%SCHEME)
-        CASE ('D95','1-L','EBA')
-          IF (JP<=1) ALLOCATE(PFIELD(INI,1,IPATCH))   
-          !* computes output physical variable
-          IF (HSURF(1:3)=='RHO') PFIELD(:,1,JP) = TZSNOW%RHO(:,1)
-          IF (HSURF(1:3)=='HEA') THEN
-            IF (TZSNOW%SCHEME=='D95'.OR.TZSNOW%SCHEME=='EBA') PFIELD(:,1,JP) = XTT-2.
-            IF (TZSNOW%SCHEME=='1-L') PFIELD(:,1,JP) = TZSNOW%T(:,1)
-          END IF
-          IF (HSURF(1:3)=='SG1') PFIELD(:,1,JP) = -20.0
-          IF (HSURF(1:3)=='SG2') PFIELD(:,1,JP) = 80.0
-          IF (HSURF(1:3)=='HIS') PFIELD(:,1,JP) = 0.0
-          IF (HSURF(1:3)=='AGE') PFIELD(:,1,JP) = 3.0
+        SELECT CASE (TZSNOW%SCHEME)
 
-        CASE ('3-L','CRO')
-          ALLOCATE(ZFIELD(INI,TZSNOW%NLAYER))
-          !* input physical variable
-          IF (HSURF(1:3)=='RHO') ZFIELD(:,:) = TZSNOW%RHO (:,1:TZSNOW%NLAYER)
-          IF (HSURF(1:3)=='AGE') ZFIELD(:,:) = TZSNOW%AGE(:,1:TZSNOW%NLAYER)
-          IF (TZSNOW%SCHEME=='CRO')THEN
-            IF (HSURF(1:3)=='SG1') ZFIELD(:,:) = TZSNOW%GRAN1(:,1:TZSNOW%NLAYER)
-            IF (HSURF(1:3)=='SG2') ZFIELD(:,:) = TZSNOW%GRAN2(:,1:TZSNOW%NLAYER)
-            IF (HSURF(1:3)=='HIS') ZFIELD(:,:) = TZSNOW%HIST(:,1:TZSNOW%NLAYER)
-          ELSE
-           IF (HSURF(1:3)=='SG1') ZFIELD(:,:) = -20.0
-           IF (HSURF(1:3)=='SG2') ZFIELD(:,:) = 80.0
-           IF (HSURF(1:3)=='HIS') ZFIELD(:,:) = 0.0                  
-          ENDIF    
-          !
-          IF ( HSURF(1:3)=='HEA') THEN
-            ALLOCATE(ZHEAT(INI,TZSNOW%NLAYER))
-            ZHEAT(:,:) = TZSNOW%HEAT(:,1:TZSNOW%NLAYER)
-            CALL SNOW_HEAT_TO_T_WLIQ(ZHEAT,TZSNOW%RHO,ZFIELD)
-            WHERE (ZFIELD>XTT.AND.ZFIELD/=XUNDEF) ZFIELD = XTT
-            DEALLOCATE(ZHEAT)
-          ENDIF
-          !
-          IF (OSNOW_IDEAL) THEN 
-            IF (JP<=1) ALLOCATE(PFIELD(INI,KLAYER,IPATCH))                  
-            PFIELD(:,:,JP) = ZFIELD(:,:)
-          ELSE
-            !
-            IF (JP<=1) ALLOCATE(PFIELD(INI,NGRID_LEVEL,IPATCH))
-            !* input snow layer thickness
-            ALLOCATE(ZDEPTH(INI,TZSNOW%NLAYER))
-            ZDEPTH(:,:) = TZSNOW%WSNOW(:,:)/TZSNOW%RHO(:,:)
+          CASE ('D95','1-L','EBA')
+            IF (JP<=1) ALLOCATE(PFIELD(INI,1,IPATCH))   
+            !* computes output physical variable
+            IF (HSURF(1:3)=='RHO') PFIELD(:,1,JP) = TZSNOW%RHO(:,1)
+            IF (HSURF(1:3)=='HEA') THEN
+              IF (TZSNOW%SCHEME=='D95'.OR.TZSNOW%SCHEME=='EBA') PFIELD(:,1,JP) = XTT-2.
+              IF (TZSNOW%SCHEME=='1-L') PFIELD(:,1,JP) = TZSNOW%T(:,1)
+            END IF
+            IF (HSURF(1:3)=='SG1') PFIELD(:,1,JP) = -20.0
+            IF (HSURF(1:3)=='SG2') PFIELD(:,1,JP) = 80.0
+            IF (HSURF(1:3)=='HIS') PFIELD(:,1,JP) = 0.0
+            IF (HSURF(1:3)=='AGE') PFIELD(:,1,JP) = 3.0
+
+          CASE ('3-L','CRO')
+            ALLOCATE(ZFIELD(INI,TZSNOW%NLAYER))
+            !* input physical variable
+            IF (HSURF(1:3)=='RHO') ZFIELD(:,:) = TZSNOW%RHO (:,1:TZSNOW%NLAYER)
+            IF (HSURF(1:3)=='AGE') ZFIELD(:,:) = TZSNOW%AGE(:,1:TZSNOW%NLAYER)
+            IF (TZSNOW%SCHEME=='CRO')THEN
+              IF (HSURF(1:3)=='SG1') ZFIELD(:,:) = TZSNOW%GRAN1(:,1:TZSNOW%NLAYER)
+              IF (HSURF(1:3)=='SG2') ZFIELD(:,:) = TZSNOW%GRAN2(:,1:TZSNOW%NLAYER)
+              IF (HSURF(1:3)=='HIS') ZFIELD(:,:) = TZSNOW%HIST(:,1:TZSNOW%NLAYER)
+            ELSE
+             IF (HSURF(1:3)=='SG1') ZFIELD(:,:) = -20.0
+             IF (HSURF(1:3)=='SG2') ZFIELD(:,:) = 80.0
+             IF (HSURF(1:3)=='HIS') ZFIELD(:,:) = 0.0                  
+            ENDIF    
             !
-            !* total depth
-            ALLOCATE(ZD(INI))
-            ZD(:) = 0.
-            DO JL=1,TZSNOW%NLAYER
-              ZD(:) = ZD(:) + ZDEPTH(:,JL)
-            ENDDO
+            IF ( HSURF(1:3)=='HEA') THEN
+              ALLOCATE(ZHEAT(INI,TZSNOW%NLAYER))
+              ZHEAT(:,:) = TZSNOW%HEAT(:,1:TZSNOW%NLAYER)
+              CALL SNOW_HEAT_TO_T_WLIQ(ZHEAT,TZSNOW%RHO,ZFIELD)
+              WHERE (ZFIELD>XTT.AND.ZFIELD/=XUNDEF) ZFIELD = XTT
+              DEALLOCATE(ZHEAT)
+            ENDIF
             !
-            !* input normalized grid
-            ALLOCATE(ZGRID(INI,TZSNOW%NLAYER))
-            DO JI=1,INI
-              IF(ZD(JI)==0.0)THEN
-                DO JL = 1,TZSNOW%NLAYER
-                  ZGRID(JI,JL)=REAL(JL)/REAL(TZSNOW%NLAYER)
-                ENDDO
-              ELSE
-                DO JL = 1,TZSNOW%NLAYER
-                  IF(JL==1)THEN
-                    ZGRID(JI,JL)=ZDEPTH(JI,JL)/ ZD(JI)
-                  ELSE
-                    ZGRID(JI,JL) = ZGRID(JI,JL-1) + ZDEPTH(JI,JL)/ZD(JI)
-                  ENDIF
-                ENDDO
-              ENDIF
-            ENDDO
-            DEALLOCATE(ZDEPTH)
-            DEALLOCATE(ZD)
-            !    
-            ! * interpolation of profile onto fine normalized snow grid
-            CALL INTERP_GRID_NAT(ZGRID(:,:),ZFIELD(:,:),    &
-                               XGRID_SNOW(:), PFIELD(:,:,JP))
-            DEALLOCATE(ZGRID)
-          ENDIF
-          DEALLOCATE(ZFIELD)
+            IF (OSNOW_IDEAL) THEN
+
+              IF (JP<=1) ALLOCATE(PFIELD(INI,KLAYER,IPATCH))                  
+              PFIELD(:,:,JP) = ZFIELD(:,:)
+
+            ELSE
+              !
+              IF (JP<=1) ALLOCATE(PFIELD(INI,NGRID_LEVEL,IPATCH))
+              !* input snow layer thickness
+              ALLOCATE(ZDEPTH(INI,TZSNOW%NLAYER))
+              ZDEPTH(:,:) = TZSNOW%WSNOW(:,:)/TZSNOW%RHO(:,:)
+              !
+              !* total depth
+              ALLOCATE(ZD(INI))
+              ZD(:) = 0.
+              DO JL=1,TZSNOW%NLAYER
+                ZD(:) = ZD(:) + ZDEPTH(:,JL)
+              ENDDO
+              !
+              !* input normalized grid
+              ALLOCATE(ZGRID(INI,TZSNOW%NLAYER))
+              DO JI=1,INI
+                IF(ZD(JI)==0.0)THEN
+                  DO JL = 1,TZSNOW%NLAYER
+                    ZGRID(JI,JL)=REAL(JL)/REAL(TZSNOW%NLAYER)
+                  ENDDO
+                ELSE
+                  DO JL = 1,TZSNOW%NLAYER
+                    IF(JL==1)THEN
+                      ZGRID(JI,JL)=ZDEPTH(JI,JL)/ ZD(JI)
+                    ELSE
+                      ZGRID(JI,JL) = ZGRID(JI,JL-1) + ZDEPTH(JI,JL)/ZD(JI)
+                    ENDIF
+                  ENDDO
+                ENDIF
+              ENDDO
+              DEALLOCATE(ZDEPTH)
+              DEALLOCATE(ZD)
+              !    
+              ! * interpolation of profile onto fine normalized snow grid
+              CALL INTERP_GRID_NAT(ZGRID(:,:),ZFIELD(:,:),XGRID_SNOW(:), PFIELD(:,:,JP))
+              DEALLOCATE(ZGRID)
+
+            ENDIF
+            DEALLOCATE(ZFIELD)
 
         END SELECT
-        !* put field form patch to all vegtypes    
+            !* put field form patch to all vegtypes    
     END SELECT
   !
+  ELSE
+    !
+    ALLOCATE(PFIELD(0,0,0))
+    !
   ENDIF
   !
   CALL DEALLOC_GR_SNOW(TZSNOW)
diff --git a/src/SURFEX/prep_teb_extern.F90 b/src/SURFEX/prep_teb_extern.F90
index 99da5d3b121432fde259094d8ef90ff349ad0d5a..bdcacd7a4eb6e38644e19fa4f4dcc748e420039f 100644
--- a/src/SURFEX/prep_teb_extern.F90
+++ b/src/SURFEX/prep_teb_extern.F90
@@ -128,8 +128,12 @@ IF (NRANK/=NPIO) INI = 0
  CALL TOWN_PRESENCE(HFILEPGDTYPE,GTEB,HDIR='-')
 !
 ALLOCATE(ZMASK(INI))
-IF (IVERSION_PGD>=7.AND.GTEB) THEN 
-  YRECFM='FRAC_TOWN'
+IF (IVERSION_PGD>=7) THEN
+  IF (GTEB) THEN 
+     YRECFM='FRAC_TOWN'
+  ELSE
+     YRECFM='FRAC_NATURE'
+  END IF  
   CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZMASK,IRESP,HDIR='A')
 ELSE
   ZMASK(:) = 1.
diff --git a/src/SURFEX/read_gr_snow.F90 b/src/SURFEX/read_gr_snow.F90
index c0d6330bc8e8ecffe1645ca6c257aeecc336a1bc..19e93845761813b7738bb2fe14425b0b11c8fc7c 100644
--- a/src/SURFEX/read_gr_snow.F90
+++ b/src/SURFEX/read_gr_snow.F90
@@ -135,7 +135,7 @@ GVERSION = (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3)
 !              -------------------
 !
 ISURFTYPE_LEN=LEN_TRIM(HSURFTYPE)
-!
+! 
 IF (KPATCH<=1) THEN
 
   IF (IVERSION <=2 .OR. (IVERSION==3 .AND. IBUGFIX<=4)) THEN
@@ -225,8 +225,8 @@ IF (IVERSION >= 7 .AND. HSURFTYPE=='VEG'.AND.KPATCH==1)  &
 !-------------------------------------------------------------------------------
 !
 !
-!*       5.    Snow reservoir
-!              --------------
+!*       5.    Snow reservoir and Snow density
+!              -------------------------------
 !
 ALLOCATE(ZWORK(KLU,INPATCH))
 !
@@ -243,7 +243,7 @@ IF (TPSNOW%SCHEME=='1-L' .OR. TPSNOW%SCHEME=='D95' .OR. TPSNOW%SCHEME=='EBA' &
   CALL READ_LAYERS(GVERSION,TPSNOW%NLAYER,YDIR,HPREFIX,YFMT,"WSNOW",HSURFTYPE,TPSNOW%WSNOW)
   CALL READ_LAYERS(GVERSION,TPSNOW%NLAYER,YDIR,HPREFIX,YFMT,"RSNOW",HSURFTYPE,TPSNOW%RHO)
   !
-  !*       7.    Snow temperature
+  !*       6.    Snow temperature
   !              ----------------
   !
   IF (TPSNOW%SCHEME=='1-L') THEN
@@ -252,19 +252,23 @@ IF (TPSNOW%SCHEME=='1-L' .OR. TPSNOW%SCHEME=='D95' .OR. TPSNOW%SCHEME=='EBA' &
     !
   ENDIF
   !
-  !*       8.    Heat content
+  !*       7.    Heat content
   !              ------------
   !
   IF (TPSNOW%SCHEME=='3-L' .OR. TPSNOW%SCHEME=='CRO') THEN
     !
     CALL READ_LAYERS(GVERSION,TPSNOW%NLAYER,YDIR,HPREFIX,YFMT,"HSNOW",HSURFTYPE,TPSNOW%HEAT)
     !
+    !
+    !*       8.    Historical parameter
+    !              -------------------
+
     IF (TPSNOW%SCHEME=='CRO') THEN
       !
       CALL READ_LAYERS(GVERSION,TPSNOW%NLAYER,YDIR,HPREFIX,YFMT,"SHIST",HSURFTYPE,TPSNOW%HIST)
       !
-      !*       9.    Snow Gran1
-      !              ------------
+      !*       9.    Snow Gran1 and  Snow Gran2 
+      !              ----------------------------
       !
       IF (GVERSION) THEN
         YFMT = "(A2,A1"//YFMT0//')'       
@@ -277,10 +281,10 @@ IF (TPSNOW%SCHEME=='1-L' .OR. TPSNOW%SCHEME=='D95' .OR. TPSNOW%SCHEME=='EBA' &
       !
     ENDIF
     !
+    !*       10.    Age parameter
+    !              -------------------
+    !
     IF ((TPSNOW%SCHEME=='3-L'.AND.IVERSION>=8) .OR. TPSNOW%SCHEME=='CRO') THEN
-      !*       12.    Age parameter
-      !              -------------------
-      !
       IF (GVERSION) THEN
         YFMT = "(A3"//YFMT0//')'         
       ELSE
@@ -303,8 +307,15 @@ IF (TPSNOW%SCHEME=='1-L' .OR. TPSNOW%SCHEME=='D95' .OR. TPSNOW%SCHEME=='EBA' &
     !
   ENDIF
   !
-  WRITE(YFMT,'(A5,I1,A2,I1,A1)') '(A4,A',ISURFTYPE_LEN,',A',IPAT_LEN,')'
-  WRITE(YRECFM,YFMT) 'ASN_',ADJUSTL(HSURFTYPE(:LEN_TRIM(HSURFTYPE))),ADJUSTL(YPAT)
+  !*       11.    Albedo
+  !              --------
+  IF (IVERSION<7 .OR. IVERSION==7 .AND. IBUGFIX<3) THEN
+    WRITE(YFMT,'(A5,I1,A2,I1,A1)')     '(A6,A',ISURFTYPE_LEN,',A',IPAT_LEN,')'
+    WRITE(YRECFM,YFMT) 'ASNOW_',ADJUSTL(HSURFTYPE(:LEN_TRIM(HSURFTYPE))),ADJUSTL(YPAT)
+  ELSE
+    WRITE(YFMT,'(A5,I1,A2,I1,A1)')     '(A4,A',ISURFTYPE_LEN,',A',IPAT_LEN,')'
+    WRITE(YRECFM,YFMT) 'ASN_',ADJUSTL(HSURFTYPE(:LEN_TRIM(HSURFTYPE))),ADJUSTL(YPAT)
+  ENDIF
   IF (GVERSION) YRECFM=ADJUSTL(HPREFIX//YRECFM)
   IF (GDIM2) THEN
     CALL READ_SURF(HPROGRAM,YRECFM,ZWORK(:,1),IRESP,HDIR=YDIR)
@@ -313,7 +324,6 @@ IF (TPSNOW%SCHEME=='1-L' .OR. TPSNOW%SCHEME=='D95' .OR. TPSNOW%SCHEME=='EBA' &
     CALL READ_SURF(HPROGRAM,YRECFM,ZWORK,IRESP,HDIR=YDIR)
     CALL PACK_SAME_RANK(KMASK_P,ZWORK(:,MAX(1,KPATCH)),TPSNOW%ALB(:))
   ENDIF
-!  IF (KLU>0) print*,YRECFM,minval(TPSNOW%ALB(:)),maxval(TPSNOW%ALB(:))
   !
 ENDIF
 !
diff --git a/src/SURFEX/read_isban.F90 b/src/SURFEX/read_isban.F90
index 8af417896027ac21554691310abaf94ea52071fe..a25ae69149f430380a601e404c5603dc7918f596 100644
--- a/src/SURFEX/read_isban.F90
+++ b/src/SURFEX/read_isban.F90
@@ -679,7 +679,7 @@ IF (IO%CRESPSL=='CNT') THEN
     YRECFM='LIGN_STR'//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
     CALL MAKE_CHOICE_ARRAY(HPROGRAM, IO%NPATCH, GDIM, YRECFM, ZWORK)
     DO JP = 1,IO%NPATCH
-      CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XSOILCARB(:,JNLITTLEVS))
+      CALL PACK_SAME_RANK(NP%AL(JP)%NR_P,ZWORK(:,JP),NPE%AL(JP)%XLIGNIN_STRUC(:,JNLITTLEVS))
     ENDDO     
   END DO
 !
diff --git a/src/SURFEX/read_nam_pgd_isba.F90 b/src/SURFEX/read_nam_pgd_isba.F90
index 85b23f284ab955d6b17bbfc4dabb3891e2cbd93d..9b90444d55c9adcb49c0f0fe30a26b954ff204a5 100644
--- a/src/SURFEX/read_nam_pgd_isba.F90
+++ b/src/SURFEX/read_nam_pgd_isba.F90
@@ -234,7 +234,7 @@ YPERMFILETYPE    = '      '
 YRUNOFFBFILETYPE = '      '
 YWDRAINFILETYPE  = '      ' 
 YPHFILETYPE      = '      '
-YPHFILETYPE      = '      '
+YFERTFILETYPE    = '      '
 !
 LIMP_CLAY        = .FALSE.
 LIMP_SAND        = .FALSE.
diff --git a/src/SURFEX/read_nam_pgd_isba_meb.F90 b/src/SURFEX/read_nam_pgd_isba_meb.F90
index 0cf00afaa9e7b9450dfb1548b237f8d0fdf297bc..7ccdfd3ea17cb8dc10c3d0300842d186d78f2e0a 100644
--- a/src/SURFEX/read_nam_pgd_isba_meb.F90
+++ b/src/SURFEX/read_nam_pgd_isba_meb.F90
@@ -43,7 +43,7 @@
 USE MODI_OPEN_NAMELIST
 USE MODI_CLOSE_NAMELIST
 !
-USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE_ECOSG
+USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE
 !
 USE MODE_POS_SURF
 !
@@ -75,7 +75,7 @@ LOGICAL                           :: GFOUND    ! flag when namelist is present
 !*    0.3    Declaration of namelists
 !            ------------------------
 !
-LOGICAL, DIMENSION(NVEGTYPE_ECOSG) :: LMEB_PATCH
+LOGICAL, DIMENSION(NVEGTYPE) :: LMEB_PATCH
 LOGICAL                :: LFORC_MEASURE
 LOGICAL                :: LMEB_LITTER
 LOGICAL                :: LMEB_GNDRES
diff --git a/src/SURFEX/read_pgd_isba_parn.F90 b/src/SURFEX/read_pgd_isba_parn.F90
index a54f1248d00e74112b5c94eaab5512336259f4ec..9926e7698acfecef17758b18dd9c14654fb37e2f 100644
--- a/src/SURFEX/read_pgd_isba_parn.F90
+++ b/src/SURFEX/read_pgd_isba_parn.F90
@@ -496,15 +496,17 @@ IF (.NOT.OLAND_USE) THEN
     ENDIF
   ENDIF
   GTIME(:) = .FALSE.
-  DO JT = 1,SIZE(GTIME)
-    IF (JT_BEG==JT_END .AND. TPDATE_END%YEAR==TPDATE_BEG%YEAR+1) THEN
-      GTIME(JT) = .TRUE.
-    ELSEIF (JT_BEG<=JT_END .AND. JT>=JT_BEG .AND. JT<=JT_END) THEN
-      GTIME(JT) = .TRUE.
-    ELSEIF (JT_BEG>JT_END .AND. (JT>=JT_BEG .OR. JT<=JT_END)) THEN
-      GTIME(JT) = .TRUE. 
-    ENDIF
-  ENDDO 
+  IF (TPDATE_END%YEAR>=TPDATE_BEG%YEAR+2) THEN ! if the run contains one whole year
+    GTIME(:) = .TRUE. ! all periods must be read
+  ELSEIF (TPDATE_END%YEAR==TPDATE_BEG%YEAR+1) THEN ! if the run passes one year
+    DO JT = 1,SIZE(GTIME)
+      IF (JT>=JT_BEG .OR. JT<=JT_END) GTIME(JT) = .TRUE. ! all periods after JT_BEG and before JT_END are read
+    ENDDO
+  ELSE ! if the run stays in one only year
+    DO JT = 1,SIZE(GTIME)
+      IF (JT>=JT_BEG .AND. JT<=JT_END) GTIME(JT) = .TRUE.  ! all periods between JT_BEG and JT_END are read
+    ENDDO 
+  ENDIF
   !
   !
   IF (DTI%LDATA_VEGTYPE) THEN
diff --git a/src/SURFEX/snow3l.F90 b/src/SURFEX/snow3l.F90
index 938ae4091ff2ece19d9bd85c4a8893e066237dcc..6d6771c3fd9944015b2b56019e195974d76f3199 100644
--- a/src/SURFEX/snow3l.F90
+++ b/src/SURFEX/snow3l.F90
@@ -545,7 +545,7 @@ ZSNOWTEMPO1(:) = ZSNOWTEMP(:,1) ! save surface snow temperature before update
 !
 ZGRNDFLUXI(:)  = ZGRNDFLUX(:)
 !
- CALL SNOW3LSOLVT(OMEB,PTSTEP,XSNOWDZMIN,PSNOWDZ,ZSCOND,ZSCAP,PTG,              &
+CALL SNOW3LSOLVT(OMEB,PTSTEP,XSNOWDZMIN,PSNOWDZ,ZSCOND,ZSCAP,PTG,              &
                    PSOILCOND,PD_G,ZRADSINK,ZCT,ZTSTERM1,ZTSTERM2,              &
                    ZPET_A_COEF_T,ZPEQ_A_COEF_T,ZPET_B_COEF_T,ZPEQ_B_COEF_T,    &
                    ZTA_IC,ZQA_IC,ZGRNDFLUX,ZGRNDFLUXO,ZSNOWTEMP,PSNOWFLUX      )  
@@ -1117,8 +1117,8 @@ IF(OMEB)THEN
 
 ! Diagnose surface layer coef (should be very close/identical to ZCOEF(:,1) computed above)
 
-   ZCOEF(:,1)           = 1.0 - PSWNETSNOWS(:)/MAX(1.E-4,PSWNETSNOW(:))
-
+   ZCOEF(:,1)           = (PSWNETSNOW(:)-PSWNETSNOWS(:))/MAX(1.E-4,PSW_RAD(:))
+   
 ELSE
 
 ! Consider 3 bands:
diff --git a/src/SURFEX/treat_field.F90 b/src/SURFEX/treat_field.F90
index 8e52f043d0a8c12454420a3a61894a7ce9f2c90f..b2d01b9aa52903e69c9aadd2e24f19f2ac29707a 100644
--- a/src/SURFEX/treat_field.F90
+++ b/src/SURFEX/treat_field.F90
@@ -491,7 +491,7 @@ SELECT CASE (HSUBROUTINE)
 
    ELSE
 
-     ALLOCATE(XSUMVAL(U%NSIZE_FULL,1))
+     ALLOCATE(XSUMVAL(U%NSIZE_FULL,SIZE(NSIZE,2)))
      IF (HFILETYPE=='DIRECT' .AND. NPROC>1) THEN
        CALL ABOR1_SFX("TREAT_FIELD: MAJ is not possible with DIRECT filetype and NPROC>1")
      ELSE
diff --git a/src/SURFEX/write_diag_pgd_isban.F90 b/src/SURFEX/write_diag_pgd_isban.F90
index 70928be6bcfd4929fbbbf754e283a5bfbe6b4e7a..96417bb19f3c71c32751e2ec36ab05c108c3f79f 100644
--- a/src/SURFEX/write_diag_pgd_isban.F90
+++ b/src/SURFEX/write_diag_pgd_isban.F90
@@ -299,9 +299,9 @@ IF(IO%CISBA=='DIF')THEN
   !* Root fraction for each patch
   !
   ALLOCATE(ZWORK1(ILU))
-  DO JP = 1,IO%NPATCH
-    PK => NP%AL(JP)
-    DO JL=1,SIZE(PK%XROOTFRAC,2)
+  DO JL=1,SIZE(PK%XROOTFRAC,2)
+    DO JP = 1,IO%NPATCH  
+      PK => NP%AL(JP)
       IF (JL<10) THEN
         WRITE(YRECFM,FMT='(A8,I1)') 'ROOTFRAC',JL
       ELSE
diff --git a/src/SURFEX/write_field_1d_patch.F90 b/src/SURFEX/write_field_1d_patch.F90
index 2907ca95fbde04fab4b0d09c58ee902d832ee090..b398135ba5fb11c7ebca37617dad971b84b957ae 100644
--- a/src/SURFEX/write_field_1d_patch.F90
+++ b/src/SURFEX/write_field_1d_patch.F90
@@ -29,7 +29,7 @@ REAL, DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: PWORK_WR
 REAL, DIMENSION(KSIZE,1) :: ZWORK
  CHARACTER(LEN=LEN_HREC) :: YRECFM
  CHARACTER(LEN=2) :: YPAT
-INTEGER :: IRESP
+INTEGER :: IRESP, ILEN
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
 IF (LHOOK) CALL DR_HOOK('WRITE_FIELD_1D_PATCH',0,ZHOOK_HANDLE)
@@ -47,6 +47,8 @@ IF (LSPLIT_PATCH) THEN
   !
 ELSE
   !
+  ILEN = LEN_TRIM(YRECFM)
+  IF (YRECFM(ILEN:ILEN)=="_") YRECFM = ADJUSTL(YRECFM(:LEN_TRIM(YRECFM)))//'P'  
   IF (KP/=0) THEN
     PWORK_WR(:,KP) = ZWORK(:,1)
     IF ( KP==SIZE(PWORK_WR,2) ) THEN
diff --git a/src/SURFEX/write_field_2d_patch.F90 b/src/SURFEX/write_field_2d_patch.F90
index f296a11cd0f930e9a90ba39ee54f90f5130b5cf6..830ff1a9a2742a6bb2485dc7913bcaa07fa1cd1a 100644
--- a/src/SURFEX/write_field_2d_patch.F90
+++ b/src/SURFEX/write_field_2d_patch.F90
@@ -33,7 +33,7 @@ REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: PWORK_WR
 REAL, DIMENSION(KSIZE,SIZE(PFIELD_IN,2)) :: ZWORK
  CHARACTER(LEN=LEN_HREC) :: YRECFM
  CHARACTER(LEN=2) :: YPAT
-INTEGER :: IRESP
+INTEGER :: IRESP, ILEN
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
 IF (LHOOK) CALL DR_HOOK('WRITE_FIELD_2D_PATCH',0,ZHOOK_HANDLE)
@@ -51,6 +51,8 @@ IF (LSPLIT_PATCH) THEN
   !
 ELSE
   !
+  ILEN = LEN_TRIM(YRECFM)
+  IF (YRECFM(ILEN:ILEN)=="_") YRECFM = ADJUSTL(YRECFM(:LEN_TRIM(YRECFM)))//'P'  
   IF (KP/=0) THEN
     PWORK_WR(:,:,KP) = ZWORK(:,:)
     IF ( KP==SIZE(PWORK_WR,3) ) THEN