Index: /LMDZ6/trunk/libf/phylmd/compbl_mod_h.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/compbl_mod_h.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/compbl_mod_h.f90	(revision 5627)
@@ -3,5 +3,5 @@
 MODULE compbl_mod_h
   IMPLICIT NONE; PRIVATE
-  PUBLIC iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
+  PUBLIC iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree, iflag_hetero_surf
 
   !!      integer iflag_pbl,iflag_pbl_split
@@ -9,6 +9,6 @@
   !!FC      integer iflag_pbl, iflag_pbl_split, iflag_order2_sollw
   !FC      common/compbl/iflag_pbl, iflag_pbl_split, iflag_order2_sollw
-  INTEGER iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
-  !$OMP THREADPRIVATE(iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree)
+  INTEGER iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree, iflag_hetero_surf
+  !$OMP THREADPRIVATE(iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree, iflag_hetero_surf)
 
   !>jyg+al1
Index: /LMDZ6/trunk/libf/phylmd/conf_phys_m.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/conf_phys_m.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/conf_phys_m.f90	(revision 5627)
@@ -193,4 +193,6 @@
     REAL,SAVE :: Cd_frein_omp
 !FC
+!AM
+    INTEGER,SAVE :: iflag_hetero_surf_omp
     INTEGER,SAVE :: iflag_order2_sollw_omp
     INTEGER, SAVE :: lev_histins_omp, lev_histLES_omp 
@@ -1531,5 +1533,12 @@
     Cd_frein_omp = 7.5E-02
     CALL getin('Cd_frein',Cd_frein_omp)
-
+!AM
+    !Config Key  = iflag_hetero_surf
+    !Config Desc = type of treatment for heterogeneous continental sub-surfaces
+    !Config Def  = 0
+    !Config Help = 0: homo. surface; 1: heteo. surface with parameter aggregation; 2: heteo surface with flux aggregation
+    !
+    iflag_hetero_surf_omp = 0
+    CALL getin('iflag_hetero_surf',iflag_hetero_surf_omp)
     !
     !Config Key  = iflag_pbl_split
@@ -2390,4 +2399,6 @@
     ifl_pbltree = ifl_pbltree_omp
     Cd_frein    =Cd_frein_omp
+!AM
+    iflag_hetero_surf = iflag_hetero_surf_omp
     iflag_order2_sollw = iflag_order2_sollw_omp
     lev_histhf = lev_histhf_omp
@@ -2839,4 +2850,6 @@
     WRITE(lunout,*) ' ifl_pbltree = ', ifl_pbltree
     WRITE(lunout,*) ' Cd_frein = ', Cd_frein
+!AM
+    WRITE(lunout,*) ' iflag_hetero_surf = ', iflag_hetero_surf
     WRITE(lunout,*) ' iflag_pbl_split = ', iflag_pbl_split
     WRITE(lunout,*) ' iflag_order2_sollw = ', iflag_order2_sollw
Index: /LMDZ6/trunk/libf/phylmd/dimphy.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dimphy.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/dimphy.f90	(revision 5627)
@@ -12,4 +12,6 @@
   INTEGER,SAVE :: klevm1
   INTEGER,SAVE :: kflev
+  INTEGER,SAVE :: nbtersrf !AM
+  INTEGER,SAVE :: nbtsoildepths !AM
 
 !$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon)
@@ -44,5 +46,5 @@
   END SUBROUTINE Init_dimphy
 
-  SUBROUTINE Init_dimphy1D(klon0,klev0)
+  SUBROUTINE Init_dimphy1D(klon0,klev0,nbtersrf0,nbtsoildepths0)
 ! 1D special version of dimphy without ALLOCATE(zmasq)
 ! which will be allocated in iniphysiq
@@ -51,5 +53,7 @@
     INTEGER, INTENT(in) :: klon0
     INTEGER, INTENT(in) :: klev0
-    
+    INTEGER, INTENT(in), OPTIONAL :: nbtersrf0
+    INTEGER, INTENT(in), OPTIONAL :: nbtsoildepths0
+
     klon=klon0
     kdlon=klon
@@ -60,5 +64,7 @@
     klevm1=klev-1
     kflev=klev
-    
+    nbtersrf=nbtersrf0
+    nbtsoildepths=nbtsoildepths0
+
   END SUBROUTINE Init_dimphy1D
 
Index: /LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h	(revision 5627)
@@ -415,4 +415,102 @@
        CALL getin('tau_soil_nudge',tau_soil_nudge)
 
+!Config  Key  = nb_ter_srf
+!Config  Desc = nb_ter_srf
+!Config  Def  = 0
+!Config  Help =
+       nb_ter_srf = 0
+       CALL getin('nb_ter_srf',nb_ter_srf)
+
+!Config  Key  = alpha_soil_ter_srf
+!Config  Desc = alpha_soil_ter_srf
+!Config  Def  = 2.
+!Config  Help =
+       alpha_soil_ter_srf = 2.
+       CALL getin('alpha_soil_ter_srf',alpha_soil_ter_srf)
+
+!Config  Key  = period_ter_srf
+!Config  Desc = period_ter_srf
+!Config  Def  = 1800.
+!Config  Help =
+       period_ter_srf = 1800.
+       CALL getin('period_ter_srf',period_ter_srf)
+
+!Config  Key  = frac_ter_srf
+!Config  Desc = frac_ter_srf
+!Config  Def  = 0.
+!Config  Help = 
+       frac_ter_srf = 0.
+       CALL getin('frac_ter_srf',frac_ter_srf)
+
+!Config  Key  = rugos_ter_srf
+!Config  Desc = rugos_ter_srf
+!Config  Def  = 0.
+!Config  Help = 
+       rugos_ter_srf = 0.
+       CALL getin('rugos_ter_srf',rugos_ter_srf)
+
+!Config  Key  = ratio_z0m_z0h_ter_srf
+!Config  Desc = ratio_z0m_z0h_ter_srf
+!Config  Def  = 10.
+!Config  Help = 
+       ratio_z0m_z0h_ter_srf = 10.
+       CALL getin('ratio_z0m_z0h_ter_srf',ratio_z0m_z0h_ter_srf)
+
+!Config  Key  = albedo_ter_srf
+!Config  Desc = albedo_ter_srf
+!Config  Def  = 0.
+!Config  Help = 
+       albedo_ter_srf = 0.
+       CALL getin('albedo_ter_srf',albedo_ter_srf)
+
+!Config  Key  = beta_ter_srf
+!Config  Desc = beta_ter_srf
+!Config  Def  = 0.
+!Config  Help = 
+       beta_ter_srf = 0.
+       CALL getin('beta_ter_srf',beta_ter_srf)
+
+!Config  Key  = inertie_ter_srf
+!Config  Desc = inertie_ter_srf
+!Config  Def  = 0.
+!Config  Help = 
+       inertie_ter_srf = 0.
+       CALL getin('inertie_ter_srf',inertie_ter_srf)
+
+!Config  Key  = hcond_ter_srf
+!Config  Desc = hcond_ter_srf
+!Config  Def  = 0.
+!Config  Help =
+       hcond_ter_srf = 0.
+       CALL getin('hcond_ter_srf',hcond_ter_srf)
+
+!Config  Key  = tsurf_ter_srf
+!Config  Desc = tsurf_ter_srf
+!Config  Def  = 283.
+!Config  Help =
+       tsurf_ter_srf = 283.
+       CALL getin('tsurf_ter_srf',tsurf_ter_srf)
+
+!Config  Key  = tsoil_ter_srf
+!Config  Desc = tsoil_ter_srf
+!Config  Def  = 283.
+!Config  Help =
+       tsoil_ter_srf = 283.
+       CALL getin('tsoil_ter_srf',tsoil_ter_srf)
+
+!Config  Key  = tsoil_depths
+!Config  Desc = tsoil_depths
+!Config  Def  = 0.
+!Config  Help =
+       tsoil_depths = 0.
+       CALL getin('tsoil_depths',tsoil_depths)
+
+!Config  Key  = nb_tsoil_depths
+!Config  Desc = nb_tsoil_depths
+!Config  Def  = 0
+!Config  Help =
+       nb_tsoil_depths = 0
+       CALL getin('nb_tsoil_depths',nb_tsoil_depths)
+
 !----------------------------------------------------------
 ! Param??tres de for??age pour les forcages communs:
@@ -631,4 +729,17 @@
       write(lunout,*)' nudging_t  = ', nudging_t
       write(lunout,*)' nudging_qv  = ', nudging_qv
+      write(lunout,*)' nb_ter_srf = ', nb_ter_srf
+      write(lunout,*)' alpha_soil_ter_srf = ', alpha_soil_ter_srf
+      write(lunout,*)' period_ter_srf = ', period_ter_srf 
+      write(lunout,*)' frac_ter_srf = ', frac_ter_srf
+      write(lunout,*)' rugos_ter_srf = ', rugos_ter_srf
+      write(lunout,*)' ratio_z0m_z0h_ter_srf = ', ratio_z0m_z0h_ter_srf
+      write(lunout,*)' albedo_ter_srf = ', albedo_ter_srf
+      write(lunout,*)' beta_ter_srf = ', beta_ter_srf
+      write(lunout,*)' inertie_ter_srf = ', inertie_ter_srf
+      write(lunout,*)' hcond_ter_srf = ', hcond_ter_srf
+      write(lunout,*)' tsurf_ter_srf = ', tsurf_ter_srf
+      write(lunout,*)' tsoil_ter_srf = ', tsoil_ter_srf
+
       IF (forcing_type .eq.40) THEN
         write(lunout,*) '--- Forcing type GCSS Old --- with:'
Index: /LMDZ6/trunk/libf/phylmd/dyn1d/compar1d_mod_h.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dyn1d/compar1d_mod_h.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/dyn1d/compar1d_mod_h.f90	(revision 5627)
@@ -8,4 +8,8 @@
           iflag_nudge, snowmass, &
           restart, ok_old_disvert, &
+          nb_ter_srf, alpha_soil_ter_srf, period_ter_srf, frac_ter_srf, &
+          rugos_ter_srf, ratio_z0m_z0h_ter_srf, albedo_ter_srf, beta_ter_srf, &
+          inertie_ter_srf, hcond_ter_srf, tsurf_ter_srf, tsoil_ter_srf, &
+          tsoil_depths, nb_tsoil_depths, &
           tadv, tadvv, tadvh, qadv, qadvv, qadvh, thadv, thadvv, thadvh, &
           trad, forc_omega, forc_u, forc_v, forc_w, forc_geo, forc_ustar, &
@@ -44,4 +48,19 @@
   LOGICAL :: ok_old_disvert
 
+  INTEGER :: nb_ter_srf
+  REAL :: alpha_soil_ter_srf
+  REAL :: period_ter_srf
+  REAL, DIMENSION(5) :: frac_ter_srf
+  REAL, DIMENSION(5) :: rugos_ter_srf
+  REAL, DIMENSION(5) :: ratio_z0m_z0h_ter_srf
+  REAL, DIMENSION(5) :: albedo_ter_srf
+  REAL, DIMENSION(5) :: beta_ter_srf
+  REAL, DIMENSION(5) :: inertie_ter_srf
+  REAL, DIMENSION(5) :: hcond_ter_srf
+  REAL, DIMENSION(5) :: tsurf_ter_srf
+  REAL, DIMENSION(5*5) :: tsoil_ter_srf
+  REAL, DIMENSION(5*5) :: tsoil_depths
+  INTEGER :: nb_tsoil_depths
+
   ! Pour les forcages communs: ces entiers valent 0 ou 1
   ! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale
@@ -65,4 +84,8 @@
   !$OMP      iflag_nudge, snowmass, &
   !$OMP      restart, ok_old_disvert, &
+  !$OMP      nb_ter_srf, frac_ter_srf, rugos_ter_srf, albedo_ter_srf,         &
+  !$OMP      beta_ter_srf, inertie_ter_srf, alpha_soil_ter_srf,               &
+  !$OMP      period_ter_srf, hcond_ter_srf, ratio_z0m_z0h_ter_srf,            &
+  !$OMP      tsurf_ter_srf, tsoil_ter_srf, tsoil_depths, nb_tsoil_depths,     &
   !$OMP      tadv, tadvv, tadvh, qadv, qadvv, qadvh, thadv, thadvv, thadvh, &
   !$OMP      trad, forc_omega, forc_u, forc_v, forc_w, forc_geo, forc_ustar, &
Index: /LMDZ6/trunk/libf/phylmd/dyn1d/scm.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dyn1d/scm.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/dyn1d/scm.f90	(revision 5627)
@@ -8,4 +8,8 @@
        clwcon, detr_therm, &
        qsol, fevap, z0m, z0h, agesno, &
+       frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, &
+       albedo_tersrf, beta_tersrf, inertie_tersrf, &
+       alpha_soil_tersrf, period_tersrf, hcond_tersrf, &
+       tsurfi_tersrf, tsoili_tersrf, tsoil_depth, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
        falb_dir, falb_dif, &
@@ -178,4 +182,7 @@
       real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
       real :: tsoil(1,nsoilmx,nbsrf)
+      ! AM
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: tsoil_ter_srf2 ! resized initial soil temperature on vertical levels (K) 
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: tsoil_depths2  ! resized soil depth at which inititial temperature is given (m)
 
 !---------------------------------------------------------------------
@@ -222,5 +229,6 @@
 !                           <> 0, tendencies of forcing are not added
       INTEGER :: flag_inhib_forcing = 0
-
+      CHARACTER(len=80) :: abort_message
+      CHARACTER(len=20) :: modname = 'scm'
 
       print*,'VOUS ENTREZ DANS LE 1D FORMAT STANDARD'
@@ -385,5 +393,5 @@
 !     call init_phys_lmdz(1,1,llm,1,(/1/)) ! job now done via iniphysiq
 !     but we still need to initialize dimphy module (klon,klev,etc.)  here.
-      call init_dimphy1D(1,llm)
+      call init_dimphy1D(1,llm,nb_ter_srf,nb_tsoil_depths)
       call suphel
       call init_infotrac
@@ -560,4 +568,48 @@
         agesno  = xagesno
         tsoil(:,:,:)=tsurf
+
+        iflag_hetero_surf = 0
+        CALL getin('iflag_hetero_surf',iflag_hetero_surf)
+
+        IF (iflag_hetero_surf .GT. 0) THEN
+          PRINT*, 'scm iflag_hetero_surf', iflag_hetero_surf
+          IF ((nbtersrf .LT. 2) .OR. (nbtersrf .GT. max_nbtersrf)) THEN
+            abort_message='The number of continental sub-surfaces (nb_ter_srf) must be between 2 and 5'
+            CALL abort_physic(modname,abort_message,1)
+          ENDIF
+          ! resized initial soil temperature on vertical levels and soil depth at which inititial temperature is given
+          ALLOCATE(tsoil_ter_srf2(nbtsoildepths,nbtersrf))
+          ALLOCATE(tsoil_depths2(nbtsoildepths,nbtersrf))
+          tsoil_ter_srf2(:,:) = 0.
+          tsoil_depths2(:,:) = 0.
+          DO i=1, nbtersrf
+            DO l=1, nbtsoildepths
+              k = nbtsoildepths*(i-1)+l
+              tsoil_ter_srf2(l,i) = tsoil_ter_srf(k)
+              tsoil_depths2(l,i) = tsoil_depths(k)
+            ENDDO
+          ENDDO
+          !
+          DO i=1, nbtersrf
+            frac_tersrf(:,i) = frac_ter_srf(i)                   ! fraction of land surface heterogeneity (-)
+            z0m_tersrf(:,i) = rugos_ter_srf(i)                   ! roughness length for momentum of land sub-surfaces (m)
+            ratio_z0m_z0h_tersrf(:,i) = ratio_z0m_z0h_ter_srf(i) ! ratio of heat to momentum roughness length of land sub-surfaces (-)
+            albedo_tersrf(:,i) = albedo_ter_srf(i)               ! albedo of land sub-surfaces (-)
+            beta_tersrf(:,i) = beta_ter_srf(i)                   ! evapotranspiration coef of land sub-surfaces (-)
+            inertie_tersrf(:,i) = inertie_ter_srf(i)             ! soil thermal inertia of land sub-surfaces (J/m2/K/s1/2)
+            hcond_tersrf(:,i) = hcond_ter_srf(i)                 ! soil heat conductivity (W/(m.K))
+            tsurfi_tersrf(:,i) = tsurf_ter_srf(i)                ! initial surface temperature (K)
+            DO l=1, nbtsoildepths
+              tsoili_tersrf(:,l,i) = tsoil_ter_srf2(l,i)         ! initial soil temperature on vertical levels (K)
+              tsoil_depth(:,l,i) = tsoil_depths2(l,i)
+            ENDDO
+          ENDDO
+          alpha_soil_tersrf = alpha_soil_ter_srf               ! ratio between the thicknesses of 2 successive layers (-)
+          period_tersrf = period_ter_srf                       ! temperature oscillation amplitude period
+          !
+          DEALLOCATE(tsoil_ter_srf2)
+          DEALLOCATE(tsoil_depths2)
+        ENDIF
+
 !-----------------------------------------------------------------------
         call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
Index: /LMDZ6/trunk/libf/phylmd/indice_sol_mod.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/indice_sol_mod.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/indice_sol_mod.f90	(revision 5627)
@@ -14,4 +14,12 @@
            INTEGER, SAVE    :: nvm_orch ! Nombre de type de vegetation ds ORCHIDEE                 
            !$OMP THREADPRIVATE(nvm_orch)
+!
+!AM heterogeneous continental sub-surfaces
+            !!! If max_nbtersrf is modified, please change also the output number in phys_output_ctrlout_mod.F90
+            INTEGER, PARAMETER :: max_nbtersrf = 5  ! maximal number of continental sub-surfaces
+            CHARACTER(len=1), DIMENSION(max_nbtersrf), PARAMETER :: nb_tersrf = (/'1', '2', '3', '4', '5'/)
+            !!! nsoilout must be lower than nsoilmx
+            INTEGER, PARAMETER :: nsoilout = 10  ! number of soil layers for output
+            CHARACTER(len=2), DIMENSION(nsoilout), PARAMETER :: nb_soil = (/'01','02','03','04','05','06','07','08','09','10'/)!,'11','12','13','14'/)
 
       END MODULE indice_sol_mod
Index: /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5627)
@@ -313,6 +313,10 @@
 !!        tke_x,     tke_w                              &
        wake_dltke,                                   &
-        treedrg                                      &
+        treedrg,                                      &
 !FC
+!AM heterogeneous continental sub-surfaces
+       tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
+       cdragm_tersrf, cdragh_tersrf, &
+       swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf &
 !!!
 #ifdef ISO
@@ -390,5 +394,13 @@
 ! pblT-----output-R- T au nveau HCL
 ! treedrg--output-R- tree drag (m)               
-!
+! qsurf_tersrf--output-R- surface specific humidity of continental sub-surfaces
+! cdragm_tersrf--output-R- momentum drag coefficient of continental sub-surfaces
+! cdragh_tersrf--output-R- heat drag coefficient of continental sub-surfaces
+! tsurf_new_tersrf--output-R- surface temperature of continental sub-surfaces
+! swnet_tersrf--output-R- net shortwave radiation of continental sub-surfaces
+! lwnet_tersrf--output-R- net longwave radiation of continental sub-surfaces
+! fluxsens_tersrf--output-R- sensible heat flux of continental sub-surfaces
+! fluxlat_tersrf--output-R- latent heat flux of continental sub-surfaces
+
     USE carbon_cycle_mod,   ONLY : carbon_cycle_cpl, carbon_cycle_tr, level_coupling_esm 
     USE carbon_cycle_mod,   ONLY : co2_send, nbcf_out, fields_out, yfields_out, cfname_out
@@ -414,5 +426,6 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
     use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &
-         dser, dt_ds, zsig, zmea
+         dser, dt_ds, zsig, zmea, &
+         frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, albedo_tersrf !AM 
     use phys_output_var_mod, only: tkt, tks, taur, sss
     use lmdz_blowing_snow_ini, only : zeta_bs
@@ -420,4 +433,5 @@
     USE netcdf, only: missing_val_netcdf => nf90_fill_real
     USE dimsoil_mod_h, ONLY: nsoilmx
+    USE surf_param_mod, ONLY: eff_surf_param  !AM
 
     USE yomcst_mod_h
@@ -620,5 +634,16 @@
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
 !FC
-    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg  ! tree drag (m)               
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg  ! tree drag (m)     
+!AM heterogeneous continental sub-surfaces
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_tersrf     ! surface temperature of continental sub-surfaces (K)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: qsurf_tersrf     ! surface specific humidity of continental sub-surfaces (kg/kg)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_new_tersrf ! surface temperature of continental sub-surfaces (K)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragm_tersrf    ! momentum drag coefficient of continental sub-surfaces (-)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragh_tersrf    ! heat drag coefficient of continental sub-surfaces (-)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: swnet_tersrf     ! net shortwave radiation of continental sub-surfaces (W/m2)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: lwnet_tersrf     ! net longwave radiation of continental sub-surfaces (W/m2)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxsens_tersrf  ! sensible heat flux of continental sub-surfaces (W/m2)               
+    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxlat_tersrf   ! latent heat flux of continental sub-surfaces (W/m2)               
+    REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf ! soil temperature of continental sub-surfaces (K)               
 #ifdef ISO        
     REAL, DIMENSION(niso,klon),   INTENT(OUT)       :: xtsol      ! water height in the soil (mm)
@@ -1050,4 +1075,7 @@
     ! dt_ds, tkt, tks, taur, sss on ocean points
     REAL :: missing_val
+    ! AM !
+    REAL, DIMENSION(klon) :: z0m_eff, z0h_eff, ratio_z0m_z0h_eff, albedo_eff
+    REAL, DIMENSION(klon, nbtersrf) :: z0h_tersrf
 #ifdef ISO
     REAL, DIMENSION(klon)       :: h1
@@ -1474,4 +1502,10 @@
     ENDDO
 
+    ! AM heterogeneous continental subsurfaces
+    ! compute time-independent effective surface parameters
+    IF (iflag_hetero_surf .GT. 0) THEN
+      albedo_eff = eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI')
+    ENDIF
+
 ! Mean calculations of albedo
 !
@@ -1486,4 +1520,10 @@
      DO nsrf = 1, nbsrf
        DO i = 1, klon
+          ! AM heterogeneous continental sub-surfaces
+          IF (nsrf .EQ. is_ter .AND. iflag_hetero_surf .GT. 0) THEN
+            alb_dir(i,k,nsrf) = albedo_eff(i)
+            alb_dif(i,k,nsrf) = albedo_eff(i)
+          ENDIF
+          !
           alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)
           alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)
@@ -1883,4 +1923,27 @@
            speed(i) = SQRT(yu(i,1)**2+yv(i,1)**2)
         ENDDO
+!       
+        !!! AM heterogeneous continental subsurfaces
+        IF (nsrf .EQ. is_ter) THEN
+          ! compute time-dependent effective surface parameters (function of zgeo1) !! AM
+          IF (iflag_hetero_surf .GT. 0) THEN
+            DO i=1,klon
+              DO j=1,nbtersrf
+                IF (ratio_z0m_z0h_tersrf(i,j) .NE. 0.) THEN
+                  z0h_tersrf(i,j) = z0m_tersrf(i,j) / ratio_z0m_z0h_tersrf(i,j)
+                ELSE
+                  z0h_tersrf(i,j) = 0.
+                ENDIF
+              ENDDO
+            ENDDO
+            !
+            z0m_eff = eff_surf_param(klon, nbtersrf, z0m_tersrf, frac_tersrf, 'CDN', zgeo1/RG)
+            z0h_eff = eff_surf_param(klon, nbtersrf, z0h_tersrf, frac_tersrf, 'CDN', zgeo1/RG)
+            yz0m = z0m_eff
+            yz0h = z0h_eff
+            !
+          ENDIF
+        ENDIF
+!
         CALL cdrag(knon, nsrf, &
             speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, &
@@ -2425,5 +2488,7 @@
                yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
                y_flux_u1, y_flux_v1, &
-               yveget,ylai,yheight   &
+               yveget,ylai,yheight, tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
+               cdragm_tersrf, cdragh_tersrf, &
+               swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf  &
 #ifdef ISO
          &      ,yxtrain_f, yxtsnow_f,yxt1, &
@@ -2432,5 +2497,7 @@
 #endif               
          &      )
- 
+
+          tsurf_tersrf(:,:) =  tsurf_new_tersrf(:,:) ! for next time step
+
 !FC quid qd yveget ylai yheight ne sont pas definit
 !FC  yveget,ylai,yheight, &
Index: /LMDZ6/trunk/libf/phylmd/phyaqua_mod.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyaqua_mod.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/phyaqua_mod.f90	(revision 5627)
@@ -355,4 +355,12 @@
     alp_bl =0.
     treedrg(:,:,:)=0.
+    tsurf_tersrf(:,:) = 0.
+    qsurf_tersrf(:,:) = 0.
+    cdragm_tersrf(:,:) = 0.
+    cdragh_tersrf(:,:) = 0.
+    swnet_tersrf(:,:) = 0.
+    lwnet_tersrf(:,:) = 0.
+    fluxsens_tersrf(:,:) = 0.
+    fluxlat_tersrf(:,:) = 0.
 
     u10m = 0.
Index: /LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90	(revision 5627)
@@ -11,5 +11,5 @@
 
   USE clesphys_mod_h
-  USE dimphy, only: klon, zmasq, klev
+  USE dimphy, only: klon, zmasq, klev, nbtersrf, nbtsoildepths
   USE iophy, ONLY : init_iophy_new
   USE ocean_cpl_mod,    ONLY : ocean_cpl_init
@@ -30,5 +30,9 @@
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
        ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, dter, dser, &
-       dt_ds, ratqs_inter_
+       dt_ds, ratqs_inter_, frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, &
+       albedo_tersrf, beta_tersrf, inertie_tersrf, alpha_soil_tersrf, &
+       period_tersrf, hcond_tersrf, tsurfi_tersrf, tsoili_tersrf, tsoil_depth, &
+       qsurf_tersrf, tsurf_tersrf, tsoil_tersrf, tsurf_new_tersrf, cdragm_tersrf, &
+       cdragh_tersrf, swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf
 !FC
   USE geometry_mod,     ONLY: longitude_deg, latitude_deg
@@ -44,4 +48,5 @@
   use netcdf, only: missing_val_netcdf => nf90_fill_real
   use config_ocean_skin_m, only: activate_ocean_skin
+  USE surf_param_mod, ONLY: average_surf_var, interpol_tsoil !AM
   USE dimsoil_mod_h, ONLY: nsoilmx
   USE yomcst_mod_h
@@ -386,4 +391,56 @@
   ENDIF
 
+  IF (iflag_hetero_surf .GT. 0) THEN
+    found=phyetat0_srf(frac_tersrf,"frac_tersrf","fraction of continental sub-surfaces",0.)
+    found=phyetat0_srf(z0m_tersrf,"z0m_tersrf","roughness length for momentum of continental sub-surfaces",0.)
+    found=phyetat0_srf(ratio_z0m_z0h_tersrf,"ratio_z0m_z0h_tersrf","ratio of heat to momentum roughness length of continental sub-surfaces",0.)
+    found=phyetat0_srf(albedo_tersrf,"albedo_tersrf","albedo of continental sub-surfaces",0.)
+    found=phyetat0_srf(beta_tersrf,"beta_tersrf","evapotranspiration coef of continental sub-surfaces",0.)
+    found=phyetat0_srf(inertie_tersrf,"inertie_tersrf","soil thermal inertia of continental sub-surfaces",0.)
+    found=phyetat0_srf(hcond_tersrf,"hcond_tersrf","heat conductivity of continental sub-surfaces",0.)
+    found=phyetat0_srf(tsurfi_tersrf,"tsurfi_tersrf","initial surface temperature of continental sub-surfaces",0.)
+    !
+    ! Check if the sum of the sub-surface fractions is equal to 1
+    DO it=1,klon
+      IF (SUM(frac_tersrf(it,:)) .NE. 1.) THEN
+        PRINT*, 'SUM(frac_tersrf) = ', SUM(frac_tersrf(it,:))
+        CALL abort_physic('conf_phys', 'the sum of fractions of heterogeneous land subsurfaces must be equal &
+                          & to 1 for iflag_hetero_surf = 1 and 2',1)
+      ENDIF
+    ENDDO
+    !
+    ! Initialisation of surface and soil temperatures (potentially different initial temperatures between sub-surfaces)
+    DO iq=1,nbtersrf
+      DO it=1,klon
+        tsurf_tersrf(it,iq) = tsurfi_tersrf(it,iq)
+      ENDDO
+    ENDDO
+    !
+    DO isoil=1, nbtsoildepths 
+      IF (isoil.GT.99) THEN
+        PRINT*, "Trop de couches "
+        CALL abort_physic("phyetat0", "", 1)
+      ENDIF
+      WRITE(str2,'(i2.2)') isoil
+      found=phyetat0_srf(tsoil_depth(:,isoil,:),"tsoil_depth"//str2//"srf","soil depth of continental sub-surfaces",0.)
+      found=phyetat0_srf(tsoili_tersrf(:,isoil,:),"Tsoili"//str2//"srf","initial soil temperature of continental sub-surfaces",0.)
+      IF (.NOT. found) THEN
+        PRINT*, "phyetat0: Le champ <Tsoili"//str2//"> est absent"
+        PRINT*, "          Il prend donc la valeur de surface"
+        tsoili_tersrf(:, isoil, :) = tsurfi_tersrf(:, :)
+      ENDIF
+    ENDDO
+    !
+    tsoil_tersrf = interpol_tsoil(klon, nbtersrf, nsoilmx, nbtsoildepths, alpha_soil_tersrf, period_tersrf, &
+                   inertie_tersrf, hcond_tersrf, tsoil_depth, tsurf_tersrf, tsoili_tersrf)
+    !
+    ! initialise also average surface and soil temperatures
+    ftsol(:,is_ter) = average_surf_var(klon, nbtersrf, tsurf_tersrf, frac_tersrf, 'ARI')
+    DO k=1, nsoilmx
+      tsoil(:,k,is_ter) = average_surf_var(klon, nbtersrf, tsoil_tersrf(:,k,:), frac_tersrf, 'ARI')
+    ENDDO
+    !
+  ENDIF ! iflag_hetero_surf > 0
+
   endif ! iflag_physiq <= 1
 
Index: /LMDZ6/trunk/libf/phylmd/phyredem.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyredem.f90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/phyredem.f90	(revision 5627)
@@ -34,5 +34,11 @@
                                 du_gwd_rando, du_gwd_front, u10m, v10m, &
                                 treedrg, solswfdiff, delta_sal, ds_ns, dt_ns, &
-                                delta_sst, ratqs_inter_, dter, dser, dt_ds
+                                delta_sst, ratqs_inter_, dter, dser, dt_ds,  &
+                                frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, &
+                                albedo_tersrf, beta_tersrf, inertie_tersrf,  &
+                                hcond_tersrf, tsurfi_tersrf, tsoili_tersrf, tsoil_depth, &
+                                qsurf_tersrf, tsurf_tersrf, tsoil_tersrf, tsurf_new_tersrf, &
+                                cdragm_tersrf, cdragh_tersrf, &
+                                swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf
 
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
@@ -189,5 +195,25 @@
 !  CALL put_field_srf2("treedrg","freinage arbres",treedrg(:,:,:))
     CALL put_field(pass,"treedrg_ter","freinage arbres",treedrg(:,:,is_ter))
-
+!AM
+    CALL put_field_srf1(pass,"frac_tersrf","fraction sous surface", frac_tersrf(:,:))
+    CALL put_field_srf1(pass,"z0m_tersrf","rugosite sous surface", z0m_tersrf(:,:))
+    CALL put_field_srf1(pass,"ratio_z0m_z0h_tersrf","ratio rugosites sous surface", ratio_z0m_z0h_tersrf(:,:))
+    CALL put_field_srf1(pass,"albedo_tersrf","albedo sous surface", albedo_tersrf(:,:))
+    CALL put_field_srf1(pass,"beta_tersrf","beta sous surface", beta_tersrf(:,:))
+    CALL put_field_srf1(pass,"inertie_tersrf","inertie sous surface", inertie_tersrf(:,:))
+    CALL put_field_srf1(pass,"hcond_tersrf","conductivité thermique sous surface", hcond_tersrf(:,:))
+    CALL put_field_srf1(pass,"tsurfi_tersrf","temperature surface sous surface initiale", tsurfi_tersrf(:,:))
+    CALL put_field_srf2(pass,"Tsoili","temperature sol sous surface initiale", tsoili_tersrf(:,:,:))
+    CALL put_field_srf2(pass,"tsoil_depth","profondeur temperature sol sous surface", tsoil_depth(:,:,:))
+    CALL put_field_srf1(pass,"qsurf_tersrf","humidite surface sous surface", qsurf_tersrf(:,:))
+    CALL put_field_srf1(pass,"tsurf_tersrf","temperature surface sous surface", tsurf_tersrf(:,:))
+    CALL put_field_srf1(pass,"tsurf_new_tersrf","temperature surface sous surface", tsurf_new_tersrf(:,:))
+    CALL put_field_srf1(pass,"cdragm_tersrf","coeff trainee quantite mouvement sous surface", cdragm_tersrf(:,:))
+    CALL put_field_srf1(pass,"cdragh_tersrf","coeff trainee chaleur sous surface", cdragh_tersrf(:,:))
+    CALL put_field_srf1(pass,"swnet_tersrf","shortwave net sous surface", swnet_tersrf(:,:))
+    CALL put_field_srf1(pass,"lwnet_tersrf","longwave net sous surface", lwnet_tersrf(:,:))
+    CALL put_field_srf1(pass,"fluxsens_tersrf","flux sensible sous surface", fluxsens_tersrf(:,:))
+    CALL put_field_srf1(pass,"fluxlat_tersrf","flux latent sous surface", fluxlat_tersrf(:,:))
+    CALL put_field_srf2(pass,"tsoil_tersrf","temperature sol sous surface", tsoil_tersrf(:,:,:))
 
     CALL put_field_srf1(pass,"QS"  , "Humidite",qsurf(:,:))
Index: /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 5627)
@@ -538,4 +538,289 @@
   TYPE(ctrl_out), SAVE :: o_tauy = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'tauy', 'Meridional wind stress', 'Pa', (/ ('', i=1, 10) /))
+
+! AM
+  !!! The number of continental sub-surfaces (max_nbtersrf) is defined in indice_sol_mod
+
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_frac_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'frac_tersrf'//nb_tersrf(1),  &
+    "Fraction of each continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'frac_tersrf'//nb_tersrf(2),  &
+    "Fraction of each continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'frac_tersrf'//nb_tersrf(3),  &
+    "Fraction of each continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'frac_tersrf'//nb_tersrf(4),  &
+    "Fraction of each continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'frac_tersrf'//nb_tersrf(5),  &
+    "Fraction of each continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_qsurf_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'qsurf_tersrf'//nb_tersrf(1),  &
+    "Surface humidity of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'qsurf_tersrf'//nb_tersrf(2),  &
+    "Surface humidity of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'qsurf_tersrf'//nb_tersrf(3),  &
+    "Surface humidity of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'qsurf_tersrf'//nb_tersrf(4),  &
+    "Surface humidity of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'qsurf_tersrf'//nb_tersrf(5),  &
+    "Surface humidity of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_tsurf_new_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsol_tersrf'//nb_tersrf(1),  &
+    "Surface temperature of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsol_tersrf'//nb_tersrf(2),  &
+    "Surface temperature of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsol_tersrf'//nb_tersrf(3),  &
+    "Surface temperature of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsol_tersrf'//nb_tersrf(4),  &
+    "Surface temperature of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsol_tersrf'//nb_tersrf(5),  &
+    "Surface temperature of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_swnet_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sols_tersrf'//nb_tersrf(1),  &
+    "Net SW radiation of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sols_tersrf'//nb_tersrf(2),  &
+    "Net SW radiation of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sols_tersrf'//nb_tersrf(3),  &
+    "Net SW radiation of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sols_tersrf'//nb_tersrf(4),  &
+    "Net SW radiation of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sols_tersrf'//nb_tersrf(5),  &
+    "Net SW radiation of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_lwnet_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'soll_tersrf'//nb_tersrf(1),  &
+    "Net LW radiation of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'soll_tersrf'//nb_tersrf(2),  &
+    "Net LW radiation of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'soll_tersrf'//nb_tersrf(3),  &
+    "Net LW radiation of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'soll_tersrf'//nb_tersrf(4),  &
+    "Net LW radiation of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'soll_tersrf'//nb_tersrf(5),  &
+    "Net LW radiation of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_fluxsens_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sens_tersrf'//nb_tersrf(1),  &
+    "Sensible heat flux of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sens_tersrf'//nb_tersrf(2),  &
+    "Sensible heat flux of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sens_tersrf'//nb_tersrf(3),  &
+    "Sensible heat flux of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sens_tersrf'//nb_tersrf(4),  &
+    "Sensible heat flux of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'sens_tersrf'//nb_tersrf(5),  &
+    "Sensible heat flux of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+ TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_fluxlat_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'flat_tersrf'//nb_tersrf(1),  &
+    "Latent heat flux of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'flat_tersrf'//nb_tersrf(2),  &
+    "Latent heat flux of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'flat_tersrf'//nb_tersrf(3),  &
+    "Latent heat flux of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'flat_tersrf'//nb_tersrf(4),  &
+    "Latent heat flux of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'flat_tersrf'//nb_tersrf(5),  &
+    "Latent heat flux of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+ TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_cdragm_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrm_tersrf'//nb_tersrf(1),  &
+    "Momentum drag coefficient of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrm_tersrf'//nb_tersrf(2),  &
+    "Momentum drag coefficient of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrm_tersrf'//nb_tersrf(3),  &
+    "Momentum drag coefficient of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrm_tersrf'//nb_tersrf(4),  &
+    "Momentum drag coefficient of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrm_tersrf'//nb_tersrf(5),  &
+    "Momentum drag coefficient of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf) :: o_cdragh_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrh_tersrf'//nb_tersrf(1),  &
+    "Heat drag coefficient of continental sub-surface "//nb_tersrf(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrh_tersrf'//nb_tersrf(2),  &
+    "Heat drag coefficient of continental sub-surface "//nb_tersrf(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrh_tersrf'//nb_tersrf(3),  &
+    "Heat drag coefficient of continental sub-surface "//nb_tersrf(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrh_tersrf'//nb_tersrf(4),  &
+    "Heat drag coefficient of continental sub-surface "//nb_tersrf(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'cdrh_tersrf'//nb_tersrf(5),  &
+    "Heat drag coefficient of continental sub-surface "//nb_tersrf(5),"K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(max_nbtersrf*nsoilout) :: o_tsoil_tersrf     = (/ &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(1),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(2),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(3),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(4),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(5),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(5),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(6),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(6),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(7),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(7),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(8),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(8),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(9),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(9),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(10),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(10),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(11),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(11),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(12),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(12),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(13),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(13),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(14),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(1)//" layer "//nb_soil(14),"K", (/ ('', i=1, 10) /)),  &
+    !
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(1),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(2),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(3),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(4),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(5),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(5),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(6),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(6),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(7),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(7),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(8),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(8),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(9),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(9),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(10),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(10),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(11),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(11),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(12),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(12),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(13),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(13),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(2)//"_l"//nb_soil(14),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(2)//" layer "//nb_soil(14),"K", (/ ('', i=1, 10) /)),  &
+    !
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(1),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(2),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(3),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(4),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(5),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(5),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(6),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(6),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(1)//"_l"//nb_soil(7),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(7),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(8),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(8),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(9),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(9),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(10),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(10),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(11),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(11),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(12),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(12),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(13),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(13),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(3)//"_l"//nb_soil(14),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(3)//" layer "//nb_soil(14),"K", (/ ('', i=1, 10) /)),  &
+    !
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(1),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(2),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(3),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(4),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(5),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(5),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(6),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(6),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(7),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(7),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(8),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(8),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(9),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(9),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(10),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(10),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(11),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(11),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(12),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(12),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(13),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(13),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(4)//"_l"//nb_soil(14),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(4)//" layer "//nb_soil(14),"K", (/ ('', i=1, 10) /)),  &
+    !
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(1),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(1),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(2),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(2),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(3),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(3),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(4),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(4),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(5),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(5),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(6),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(6),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(7),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(7),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(8),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(8),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(9),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(9),"K", (/ ('', i=1, 10) /)),  &
+    ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(10),  &
+    "Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(10),"K", (/ ('', i=1, 10) /)) /)
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(11),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(11),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(12),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(12),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(13),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(13),"K", (/ ('', i=1, 10) /)),  &
+    !ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/),'tsoil_tersrf'//nb_tersrf(5)//"_l"//nb_soil(14),  &
+    !"Soil temperature of continental sub-surface "//nb_tersrf(5)//" layer "//nb_soil(14),"K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(nsoilout) :: o_ftsoil = (/ &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(1), &
+    'Continental soil temperature layer '//nb_soil(1), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(2), &
+    'Continental soil temperature layer '//nb_soil(2), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(3), &
+    'Continental soil temperature layer '//nb_soil(3), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(4), &
+    'Continental soil temperature layer '//nb_soil(4), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(5), &
+    'Continental soil temperature layer '//nb_soil(5), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(6), &
+    'Continental soil temperature layer '//nb_soil(6), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(7), &
+    'Continental soil temperature layer '//nb_soil(7), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(8), &
+    'Continental soil temperature layer '//nb_soil(8), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(9), &
+    'Continental soil temperature layer '//nb_soil(9), 'K', (/ ('', i=1, 10) /)), &
+    ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(10), &
+    'Continental soil temperature layer '//nb_soil(10), 'K', (/ ('', i=1, 10) /)) /)
+    !ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(11), &
+    !'Continental soil temperature layer '//nb_soil(11), 'K', (/ ('', i=1, 10) /)), &
+    !ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(12), &
+    !'Continental soil temperature layer '//nb_soil(12), 'K', (/ ('', i=1, 10) /)), &
+    !ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(13), &
+    !'Continental soil temperature layer '//nb_soil(13), 'K', (/ ('', i=1, 10) /)), &
+    !ctrl_out((/ 1, 2, 10, 10, 10, 10, 11, 11, 11, 11/), 'tsoil'//nb_soil(14), &
+    !'Continental soil temperature layer '//nb_soil(14), 'K', (/ ('', i=1, 10) /)) /)
+! AM
 
 !AI Ecrad 3Deffect
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 5627)
@@ -25,5 +25,5 @@
     ! defined and initialised in phys_output_mod.F90
 
-    USE dimphy, ONLY: klon, klev, klevp1
+    USE dimphy, ONLY: klon, klev, klevp1, nbtersrf
     USE infotrac_phy, ONLY: nbtr, nqtot, nqo, type_trac, tracers, niso, ntiso
     USE strings_mod,  ONLY: maxlen
@@ -65,5 +65,5 @@
          o_bils_ec,o_bils_ech, o_bils_tke, o_bils_kinetic, &
          o_bils_latent, o_bils_enthalp, o_sens, &
-         o_fder, o_ffonte, o_fqcalving, o_fqfonte, o_mrroli, o_runofflic, &
+         o_fder, o_ftsoil, o_ffonte, o_fqcalving, o_fqfonte, o_mrroli, o_runofflic, &
          o_taux, o_tauy, o_snowsrf, o_qsnow, &
 ! SN runoff_diag
@@ -203,4 +203,9 @@
          o_lat_prec_sol_oce, o_lat_prec_sol_sic, &
          o_sza, &
+! AM
+         o_frac_tersrf, o_qsurf_tersrf, o_tsurf_new_tersrf, &
+         o_cdragm_tersrf, o_cdragh_tersrf, &
+         o_swnet_tersrf, o_lwnet_tersrf, o_fluxsens_tersrf, o_fluxlat_tersrf, &
+         o_tsoil_tersrf, &
 ! Marine
          o_map_prop_hc, o_map_prop_hist, o_map_emis_hc, o_map_iwp_hc, &
@@ -303,4 +308,6 @@
          rhlevSTD, O3STD, O3daySTD, uvSTD, vqSTD, vTSTD, wqSTD, vphiSTD, &
          wTSTD, u2STD, v2STD, T2STD, missing_val_nf90, delta_sal, ds_ns, &
+         frac_tersrf, qsurf_tersrf, tsurf_new_tersrf, cdragm_tersrf, cdragh_tersrf, &
+         swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf, tsoil_tersrf, &
 #ifdef ISO
          xtrain_con, xtsnow_con, xtrain_fall, xtsnow_fall, fxtevap, &
@@ -459,6 +466,6 @@
     USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, &
         slab_ekman,slab_hdiff,slab_gm,dt_ekman, dt_hdiff, dt_gm, dt_qflux
-    USE pbl_surface_mod, ONLY: snow
-    USE indice_sol_mod, ONLY: nbsrf
+    USE pbl_surface_mod, ONLY: snow, ftsoil
+    USE indice_sol_mod, ONLY: nbsrf, nsoilout
 #ifdef ISO
     USE isotopes_mod, ONLY: iso_HTO, isoName
@@ -518,5 +525,5 @@
     ! Local
     INTEGER :: itau_w
-    INTEGER :: i, iinit, iinitend=1, iff, iq, nsrf, k, ll, naero
+    INTEGER :: i, iinit, iinitend=1, iff, iq, nsrf, k, ll, naero, j
     REAL, DIMENSION (klon) :: zx_tmp_fi2d, zpt_conv2d, wind100m
     REAL, DIMENSION (klon,klev) :: zx_tmp_fi3d, zpt_conv
@@ -944,4 +951,43 @@
        CALL histwrite_phy(o_topl, toplw)
        CALL histwrite_phy(o_topl0, toplw0)
+
+       !AM heterogeneous continental sub-surfaces
+       IF (iflag_hetero_surf .EQ. 2) THEN
+         iq = 0
+         DO j = 1, nbtersrf
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = frac_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_frac_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = qsurf_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_qsurf_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = tsurf_new_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_tsurf_new_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = cdragm_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_cdragm_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = cdragh_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_cdragh_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = swnet_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_swnet_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = lwnet_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_lwnet_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = fluxsens_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_fluxsens_tersrf(j), zx_tmp_fi2d)
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = fluxlat_tersrf( 1 : klon, j)
+           CALL histwrite_phy(o_fluxlat_tersrf(j), zx_tmp_fi2d)
+           !
+           DO k = 1, nsoilout
+             iq = iq + 1
+             IF (vars_defined) zx_tmp_fi2d(1 : klon) = tsoil_tersrf( 1 : klon, k, j)
+             CALL histwrite_phy(o_tsoil_tersrf(iq), zx_tmp_fi2d)
+           ENDDO
+         ENDDO
+       ENDIF
+       ! add tsoil as output
+       IF (iflag_hetero_surf .GT. 0) THEN
+         DO k = 1, nsoilout
+           IF (vars_defined) zx_tmp_fi2d(1 : klon) = ftsoil( 1 : klon, k, is_ter)
+           CALL histwrite_phy(o_ftsoil(k), zx_tmp_fi2d)
+         ENDDO
+       ENDIF
+       !AM
 
 ! offline 
Index: /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 5627)
@@ -10,4 +10,5 @@
 ! Declaration des variables
       USE dimphy
+      USE dimsoil_mod_h, ONLY: nsoilmx
       USE netcdf, only: nf90_fill_real
       INTEGER, PARAMETER :: nlevSTD=17
@@ -39,4 +40,49 @@
       REAL, ALLOCATABLE, SAVE :: treedrg(:,:,:)
 !$OMP THREADPRIVATE(treedrg)
+!AM land surface heterogeneities 
+      REAL, SAVE :: alpha_soil_tersrf
+!$OMP THREADPRIVATE(alpha_soil_tersrf)
+      REAL, SAVE :: period_tersrf
+!$OMP THREADPRIVATE(period_tersrf)
+      REAL, ALLOCATABLE, SAVE :: frac_tersrf(:,:)
+!$OMP THREADPRIVATE(frac_tersrf)
+      REAL, ALLOCATABLE, SAVE :: z0m_tersrf(:,:)
+!$OMP THREADPRIVATE(z0m_tersrf)
+      REAL, ALLOCATABLE, SAVE :: ratio_z0m_z0h_tersrf(:,:)
+!$OMP THREADPRIVATE(ratio_z0m_z0h_tersrf)
+      REAL, ALLOCATABLE, SAVE :: albedo_tersrf(:,:)
+!$OMP THREADPRIVATE(albedo_tersrf)
+      REAL, ALLOCATABLE, SAVE :: beta_tersrf(:,:)
+!$OMP THREADPRIVATE(beta_tersrf)
+      REAL, ALLOCATABLE, SAVE :: inertie_tersrf(:,:)
+!$OMP THREADPRIVATE(inertie_tersrf)
+      REAL, ALLOCATABLE, SAVE :: hcond_tersrf(:,:)
+!$OMP THREADPRIVATE(hcond_tersrf)
+      REAL, ALLOCATABLE, SAVE :: tsurfi_tersrf(:,:)
+!$OMP THREADPRIVATE(tsurfi_tersrf)
+      REAL, ALLOCATABLE, SAVE :: tsoili_tersrf(:,:,:)
+!$OMP THREADPRIVATE(tsoili_tersrf)
+      REAL, ALLOCATABLE, SAVE :: tsoil_depth(:,:,:)
+!$OMP THREADPRIVATE(tsoil_depth)
+      REAL, ALLOCATABLE, SAVE :: tsurf_tersrf(:,:)
+!$OMP THREADPRIVATE(tsurf_tersrf)
+     REAL, ALLOCATABLE, SAVE :: tsoil_tersrf(:,:,:)
+!$OMP THREADPRIVATE(tsoil_tersrf)
+      REAL, ALLOCATABLE, SAVE :: qsurf_tersrf(:,:)
+!$OMP THREADPRIVATE(qsurf_tersrf)
+      REAL, ALLOCATABLE, SAVE :: tsurf_new_tersrf(:,:)
+!$OMP THREADPRIVATE(tsurf_new_tersrf)
+      REAL, ALLOCATABLE, SAVE :: cdragm_tersrf(:,:)
+!$OMP THREADPRIVATE(cdragm_tersrf)
+      REAL, ALLOCATABLE, SAVE :: cdragh_tersrf(:,:)
+!$OMP THREADPRIVATE(cdragh_tersrf)
+      REAL, ALLOCATABLE, SAVE :: swnet_tersrf(:,:)
+!$OMP THREADPRIVATE(swnet_tersrf)
+      REAL, ALLOCATABLE, SAVE :: lwnet_tersrf(:,:)
+!$OMP THREADPRIVATE(lwnet_tersrf)
+      REAL, ALLOCATABLE, SAVE :: fluxsens_tersrf(:,:)
+!$OMP THREADPRIVATE(fluxsens_tersrf)
+      REAL, ALLOCATABLE, SAVE :: fluxlat_tersrf(:,:)
+!$OMP THREADPRIVATE(fluxlat_tersrf)
 
 !      character(len=6), SAVE :: ocean
@@ -563,4 +609,26 @@
 !FC
       ALLOCATE(treedrg(klon,klev,nbsrf))
+!AM
+      ALLOCATE(frac_tersrf(klon,nbtersrf))
+      ALLOCATE(z0m_tersrf(klon,nbtersrf))
+      ALLOCATE(ratio_z0m_z0h_tersrf(klon,nbtersrf))
+      ALLOCATE(albedo_tersrf(klon,nbtersrf))
+      ALLOCATE(beta_tersrf(klon,nbtersrf))
+      ALLOCATE(inertie_tersrf(klon,nbtersrf))
+      ALLOCATE(hcond_tersrf(klon,nbtersrf))
+      ALLOCATE(tsurfi_tersrf(klon,nbtersrf))
+      ALLOCATE(tsoili_tersrf(klon,nbtsoildepths,nbtersrf))
+      ALLOCATE(tsoil_depth(klon,nbtsoildepths,nbtersrf))
+      ALLOCATE(tsurf_tersrf(klon,nbtersrf))
+      ALLOCATE(tsoil_tersrf(klon,nsoilmx,nbtersrf))
+      ALLOCATE(qsurf_tersrf(klon,nbtersrf))
+      ALLOCATE(tsurf_new_tersrf(klon,nbtersrf))
+      ALLOCATE(cdragm_tersrf(klon,nbtersrf))
+      ALLOCATE(cdragh_tersrf(klon,nbtersrf))
+      ALLOCATE(swnet_tersrf(klon,nbtersrf))
+      ALLOCATE(lwnet_tersrf(klon,nbtersrf))
+      ALLOCATE(fluxsens_tersrf(klon,nbtersrf))
+      ALLOCATE(fluxlat_tersrf(klon,nbtersrf))
+
       ALLOCATE(falb1(klon,nbsrf))
       ALLOCATE(falb2(klon,nbsrf))
@@ -806,4 +874,25 @@
 !FC
       DEALLOCATE(treedrg)
+!AM
+      DEALLOCATE(frac_tersrf)
+      DEALLOCATE(z0m_tersrf)
+      DEALLOCATE(ratio_z0m_z0h_tersrf)
+      DEALLOCATE(albedo_tersrf)
+      DEALLOCATE(beta_tersrf)
+      DEALLOCATE(inertie_tersrf)
+      DEALLOCATE(hcond_tersrf)
+      DEALLOCATE(tsurfi_tersrf)
+      DEALLOCATE(tsoili_tersrf)
+      DEALLOCATE(tsoil_depth)
+      DEALLOCATE(tsurf_tersrf)
+      DEALLOCATE(tsoil_tersrf)
+      DEALLOCATE(qsurf_tersrf)
+      DEALLOCATE(tsurf_new_tersrf)
+      DEALLOCATE(cdragm_tersrf)
+      DEALLOCATE(cdragh_tersrf)
+      DEALLOCATE(swnet_tersrf)
+      DEALLOCATE(lwnet_tersrf)
+      DEALLOCATE(fluxsens_tersrf)
+      DEALLOCATE(fluxlat_tersrf)
       DEALLOCATE(rain_fall, snow_fall, bs_fall,solsw, solswfdiff, sollw, radsol, swradcorr)
       DEALLOCATE(zmea, zstd, zsig, zgam)
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5627)
@@ -2932,6 +2932,10 @@
             wake_delta_pbl_TKE, &
                                 !>nrlmd+jyg
-             treedrg )
+             treedrg, &
 !FC
+!AM 
+            tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
+            cdragm_tersrf, cdragh_tersrf, &
+            swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)
        !
        !  Add turbulent diffusion tendency to the wake difference variables
Index: /LMDZ6/trunk/libf/phylmd/soil_hetero.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/soil_hetero.F90	(revision 5627)
+++ /LMDZ6/trunk/libf/phylmd/soil_hetero.F90	(revision 5627)
@@ -0,0 +1,251 @@
+!
+! $Header$
+!
+SUBROUTINE soil_hetero(ptimestep, indice, knon, snow, ptsrf, qsol, &
+     lon, lat, ptsoil, pcapcal, pfluxgrd, ztherm_i, conv_ratio)
+
+  USE YOMCST, ONLY: RTT, RPI
+  USE dimsoil_mod_h, ONLY: nsoilmx
+  USE comsoil_mod_h
+  USE compbl_mod_h
+  USE dimpft_mod_h
+  USE dimphy
+  USE mod_phys_lmdz_para
+  USE indice_sol_mod
+  USE print_control_mod, ONLY: lunout
+  USE phys_state_var_mod, ONLY: alpha_soil_tersrf, period_tersrf 
+  USE surf_param_mod, ONLY: eff_surf_param
+
+  IMPLICIT NONE
+
+!=======================================================================
+!
+!   Auteur:  Frederic Hourdin     30/01/92
+!   -------
+!
+!   Object:  Computation of : the soil temperature evolution
+!   -------                   the surfacic heat capacity "Capcal"
+!                            the surface conduction flux pcapcal
+!
+!   Update: 2021/07 : soil thermal inertia, formerly a constant value,
+!   ------   can also be now a function of soil moisture (F Cheruy's idea)
+!            depending on iflag_inertie, read from physiq.def via conf_phys_m.F90
+!            ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy)
+!           2025/04 : A. Maison, adapting the routine for heterogeneous continental sub-surfaces
+!
+!   Method: Implicit time integration
+!   -------
+!   Consecutive ground temperatures are related by:
+!           T(k+1) = C(k) + D(k)*T(k)  (*)
+!   The coefficients C and D are computed at the t-dt time-step.
+!   Routine structure:
+!   1) C and D coefficients are computed from the old temperature
+!   2) new temperatures are computed using (*)
+!   3) C and D coefficients are computed from the new temperature
+!      profile for the t+dt time-step
+!   4) the coefficients A and B are computed where the diffusive
+!      fluxes at the t+dt time-step is given by
+!             Fdiff = A + B Ts(t+dt)
+!      or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
+!             with F0 = A + B (Ts(t))
+!                 Capcal = B*dt
+!           
+!   Interface:
+!   ----------
+!
+!   Arguments:
+!   ----------
+!   ptimestep            physical timestep (s)
+!   indice               sub-surface index
+!   snow(klon)           snow
+!   ptsrf(klon)          surface temperature at time-step t (K)
+!   qsol(klon)           soil moisture (kg/m2 or mm)
+!   lon(klon)            longitude in radian
+!   lat(klon)            latitude in radian
+!   ptsoil(klon,nsoilmx) temperature inside the ground (K)
+!   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
+!   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
+!   ztherm_i(klon)       soil thermal inertia (J.m-2.K.s-1/2)
+!   conv_ratio(klon)     ratio to convert soil depths in meters (-)
+!   
+!=======================================================================
+!-----------------------------------------------------------------------
+! Arguments
+! ---------
+  REAL, INTENT(IN)                     :: ptimestep
+  INTEGER, INTENT(IN)                  :: indice, knon !, knindex
+  REAL, DIMENSION(klon), INTENT(IN)    :: snow
+  REAL, DIMENSION(klon), INTENT(IN)    :: ptsrf
+  REAL, DIMENSION(klon), INTENT(IN)    :: qsol
+  REAL, DIMENSION(klon), INTENT(IN)    :: lon
+  REAL, DIMENSION(klon), INTENT(IN)    :: lat
+  REAL, DIMENSION(klon), INTENT(IN)    :: ztherm_i
+  REAL, DIMENSION(klon), INTENT(IN)    :: conv_ratio
+
+  REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil
+  REAL, DIMENSION(klon), INTENT(OUT)           :: pcapcal
+  REAL, DIMENSION(klon), INTENT(OUT)           :: pfluxgrd
+
+!-----------------------------------------------------------------------
+! Local variables
+! ---------------
+  INTEGER                             :: ig, jk, ierr
+  REAL, DIMENSION(nsoilmx)            :: zdz2
+  REAL                                :: z1s
+  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef
+
+! Local saved variables
+! ---------------------
+  REAL, SAVE                     :: lambda
+!$OMP THREADPRIVATE(lambda)
+  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
+!$OMP THREADPRIVATE(dz1,dz2)
+  LOGICAL, SAVE                  :: firstcall=.TRUE.
+!$OMP THREADPRIVATE(firstcall)
+    
+!-----------------------------------------------------------------------
+!   Depthts:
+!   --------
+  REAL fz,rk,fz1,rk1,rk2
+  fz(rk)=fz1*(alpha_soil_tersrf**rk-1.)/(alpha_soil_tersrf-1.)
+!
+!-----------------------------------------------------------------------
+! Calculation of some constants
+! NB! These constants do not depend on the sub-surfaces
+!-----------------------------------------------------------------------
+
+  IF (firstcall) THEN
+!-----------------------------------------------------------------------
+!   ground levels 
+!   grnd=z/l where l is the skin depth of the diurnal cycle:
+!-----------------------------------------------------------------------
+!   la premiere couche represente un dixieme de cycle diurne
+     fz1=SQRT(period_tersrf/RPI)
+     
+     DO jk=1,nsoilmx
+        rk1=jk
+        rk2=jk-1
+        dz2(jk)=fz(rk1)-fz(rk2)
+     ENDDO
+     DO jk=1,nsoilmx-1
+        rk1=jk+.5
+        rk2=jk-.5
+        dz1(jk)=1./(fz(rk1)-fz(rk2))
+     ENDDO
+     lambda=fz(.5)*dz1(1)
+     WRITE(lunout,*) 'surface index:', indice
+     WRITE(lunout,*)'full layers, intermediate layers (seconds)'
+     DO jk=1,nsoilmx
+        rk=jk
+        rk1=jk+.5
+        rk2=jk-.5
+        WRITE(lunout,*)'fz=', &
+             fz(rk1)*fz(rk2)*RPI,fz(rk)*fz(rk)*RPI
+     ENDDO
+     WRITE(lunout,*)'full layers, intermediate layers (meters)'
+     DO jk=1,nsoilmx
+        rk=jk
+        rk2=jk-.5
+        WRITE(lunout,*)'fz=', &
+             fz(rk2)*conv_ratio, fz(rk)*conv_ratio
+     ENDDO
+     firstcall =.FALSE.
+  END IF
+
+!-----------------------------------------------------------------------
+! 1)
+! Calculation of Cgrf and Dgrd coefficients using soil temperature from 
+! previous time step.
+!
+! These variables are recalculated on the local compressed grid instead 
+! of saved in restart file.
+!-----------------------------------------------------------------------
+  DO jk=1,nsoilmx
+     zdz2(jk)=dz2(jk)/ptimestep
+  ENDDO
+  
+  DO ig=1,knon
+     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
+     C_coef(ig,nsoilmx-1,indice)= &
+          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
+     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
+  ENDDO
+  
+  DO jk=nsoilmx-1,2,-1
+     DO ig=1,knon
+        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
+             *(1.-D_coef(ig,jk,indice)))
+        C_coef(ig,jk-1,indice)= &
+             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
+        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
+     ENDDO
+  ENDDO
+
+!-----------------------------------------------------------------------
+! 2)
+! Computation of the soil temperatures using the Cgrd and Dgrd
+! coefficient computed above
+!
+!-----------------------------------------------------------------------
+
+!    Surface temperature
+  DO ig=1,knon
+     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
+          (lambda*(1.-D_coef(ig,1,indice))+1.)
+  ENDDO
+  
+!   Other temperatures
+  DO jk=1,nsoilmx-1
+     DO ig=1,knon
+        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
+             *ptsoil(ig,jk)
+     ENDDO
+  ENDDO
+
+  IF (indice == is_sic) THEN
+     DO ig = 1 , knon
+        ptsoil(ig,nsoilmx) = RTT - 1.8
+     END DO
+  ENDIF
+
+!-----------------------------------------------------------------------
+! 3)
+! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil 
+! temperature
+!-----------------------------------------------------------------------
+  DO ig=1,knon
+     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
+     C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
+     D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
+  ENDDO
+  
+  DO jk=nsoilmx-1,2,-1
+     DO ig=1,knon
+        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
+             *(1.-D_coef(ig,jk,indice)))
+        C_coef(ig,jk-1,indice) = &
+             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
+        D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s
+     ENDDO
+  ENDDO
+
+!-----------------------------------------------------------------------
+! 4)
+! Computation of the surface diffusive flux from ground and
+! calorific capacity of the ground
+!-----------------------------------------------------------------------
+  DO ig=1,knon
+     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
+          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
+     pcapcal(ig)  = ztherm_i(ig)* &
+          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
+     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
+     pcapcal(ig)  = pcapcal(ig)/z1s
+     pfluxgrd(ig) = pfluxgrd(ig) &
+          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
+          - lambda * C_coef(ig,1,indice) &
+          - ptsrf(ig)) &
+          /ptimestep
+  ENDDO
+    
+END SUBROUTINE soil_hetero
Index: /LMDZ6/trunk/libf/phylmd/surf_land_bucket_hetero_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/surf_land_bucket_hetero_mod.F90	(revision 5627)
+++ /LMDZ6/trunk/libf/phylmd/surf_land_bucket_hetero_mod.F90	(revision 5627)
@@ -0,0 +1,386 @@
+!
+MODULE surf_land_bucket_hetero_mod
+!
+! 2025/04 A. Maison (adapted from surf_land_bucket_mod)
+! Surface land bucket module with heterogeneous continental sub-surfaces
+! This module is used when no external land model is choosen and iflag_hetero_surf = 1 or 2.
+!
+  IMPLICIT NONE
+
+CONTAINS
+
+  SUBROUTINE surf_land_bucket_hetero(itime, jour, knon, knindex, debut, dtime,&
+       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, &
+       spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, pref, plev, &
+       u1, v1, gustiness, rugoro, swnet, lwnet, &
+       snow, qsol, agesno, tsoil, &
+       qsurf, z0m, z0h, alb1_new, alb2_new, evap, &
+       fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l, &
+       tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
+       cdragm_tersrf, cdragh_tersrf, &
+       swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)
+
+    USE clesphys_mod_h
+    USE dimsoil_mod_h, ONLY: nsoilmx
+    USE YOMCST, ONLY: RD, RG, RCPD, RSIGMA
+    USE compbl_mod_h
+    USE dimpft_mod_h
+    USE limit_read_mod
+    USE surface_data
+    USE fonte_neige_mod
+    USE calcul_fluxs_mod
+    USE cpl_mod
+    USE dimphy
+    USE geometry_mod, ONLY: longitude,latitude 
+    USE mod_grid_phy_lmdz
+    USE mod_phys_lmdz_para
+    USE indice_sol_mod
+    USE phys_state_var_mod, ONLY: frac_tersrf, ratio_z0m_z0h_tersrf, z0m_tersrf, &
+                                  albedo_tersrf, beta_tersrf, inertie_tersrf, &
+                                  hcond_tersrf
+    USE surf_param_mod, ONLY: eff_surf_param, average_surf_var
+    USE cdrag_mod
+
+!****************************************************************************************
+! Bucket calculations for surface. 
+!****************************************************************************************
+!
+! Input variables  
+!****************************************************************************************
+    INTEGER, INTENT(IN)                     :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
+    LOGICAL, INTENT(IN)                     :: debut
+    REAL, INTENT(IN)                        :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
+    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)       :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)       :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)       :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)       :: pref
+    REAL, DIMENSION(klon), INTENT(IN)       :: u1, v1, gustiness
+    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
+    REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
+    REAL, DIMENSION(klon, nbtersrf), INTENT(IN) :: tsurf_tersrf
+
+
+! In/Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: plev
+    REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf
+
+
+! Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new, alb2_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
+    !
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: tsurf_new_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: qsurf_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragm_tersrf, cdragh_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: swnet_tersrf, lwnet_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxsens_tersrf, fluxlat_tersrf
+
+! Local variables
+!****************************************************************************************
+    REAL, DIMENSION(klon) :: soilcap, soilflux
+    REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
+    REAL, DIMENSION(klon) :: alb_neig, alb_lim, icesub
+    REAL, DIMENSION(klon) :: zfra
+    REAL, DIMENSION(klon) :: radsol
+    REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay
+    REAL, DIMENSION(klon) :: dummy_riverflow,dummy_coastalflow 
+    INTEGER               :: i, j, k
+    !
+    REAL, DIMENSION(klon)           :: zlev, geop1, speed, pblh, ri_in, sst
+    REAL, DIMENSION(klon)           :: beta_eff, inertie_eff, conv_ratio_eff
+    REAL, DIMENSION(klon)           :: meansqT
+    REAL, DIMENSION(klon, nbtersrf) :: z0h_tersrf, emis_tersrf, conv_ratio_tersrf
+    REAL, DIMENSION(klon, nbtersrf) :: evap_tersrf
+    REAL, DIMENSION(klon, nbtersrf) :: dflux_s_tersrf, dflux_l_tersrf
+    REAL, DIMENSION(klon, nbtersrf) :: radsol_tersrf
+    REAL, DIMENSION(klon, nbtersrf) :: zri_tersrf
+    REAL, PARAMETER                 :: klon_1D = 1
+
+!
+!****************************************************************************************
+
+    ! *** Calculations common to the two flag values ***
+
+    ! average albedo
+    alb_lim = eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI')
+
+    ! suppose zero surface speed
+    u0(:)=0.0
+    v0(:)=0.0
+    u1_lay(:) = u1(:) - u0(:)
+    v1_lay(:) = v1(:) - v0(:)
+    speed(:) = (u1_lay(:)**2 + v1_lay(:)**2)**0.5
+    !
+    geop1(1:knon) = RD * temp_air(1:knon) / (0.5*(pref(1:knon)+p1lay(1:knon))) &
+            * (pref(1:knon)-p1lay(1:knon))
+    zlev(1:knon) = (plev(1:knon)*RD*temp_air(1:knon)/(pref(1:knon)*RG))/2.
+    !
+    ! compute roughness lengths
+    DO i=1, knon
+      DO j=1, nbtersrf
+        IF (ratio_z0m_z0h_tersrf(i,j) .NE. 0.) THEN
+          z0h_tersrf(i,j) = z0m_tersrf(i,j) / ratio_z0m_z0h_tersrf(i,j)
+        ELSE
+          z0h_tersrf(i,j) = 1.E-12
+        ENDIF
+      ENDDO
+    ENDDO
+
+    z0m = eff_surf_param(klon, nbtersrf, z0m_tersrf, frac_tersrf, 'CDN', zlev)
+    z0h = eff_surf_param(klon, nbtersrf, z0h_tersrf, frac_tersrf, 'CDN', zlev)
+
+    DO i=1, knon
+       z0m(i) = MAX(1.5e-05,SQRT(z0m(i)**2 + rugoro(i)**2))
+    END DO
+
+    ! compute the ratio to convert and print soil depths in meters (conv_ratio = (cond/cap)^0.5 and cap = I^2/cond)
+    DO j=1, nbtersrf
+       conv_ratio_tersrf(:,j) = hcond_tersrf(:,j)/inertie_tersrf(:,j)
+    ENDDO
+
+    !
+    ! *** Surface parameter aggregation ***
+    !
+    IF (iflag_hetero_surf == 1) THEN
+      !* Calcultaion of fluxes 
+      
+      ! calculate total absorbed radiance at surface
+      radsol(:) = 0.0
+      radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
+      ! calculate constants (needeed for capsol and dif_grnd)
+      CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
+      IF (type_veget=='betaclim') THEN
+        CALL calbeta_clim(knon,jour,latitude(knindex(1:knon)),beta)
+      ENDIF
+
+      ! mean evapotranspiration coefficient
+      beta_eff = eff_surf_param(klon, nbtersrf, beta_tersrf, frac_tersrf, 'ARI')
+      beta = beta_eff
+
+      ! calculate temperature, heat capacity and conduction flux in soil
+      IF (soil_model) THEN
+         inertie_eff = eff_surf_param(klon, nbtersrf, inertie_tersrf, frac_tersrf, 'ARI')
+         conv_ratio_eff = eff_surf_param(klon, nbtersrf, conv_ratio_tersrf, frac_tersrf, 'ARI')
+         !
+         CALL soil_hetero(dtime, is_ter, knon, snow, tsurf, qsol,  & 
+           & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux, &
+           & inertie_eff, conv_ratio_eff)
+         !
+         DO i=1, knon
+            cal(i) = RCPD / soilcap(i)
+            radsol(i) = radsol(i)  + soilflux(i)
+         END DO
+      ELSE 
+         cal(:) = RCPD * capsol(:)
+         IF (klon_glo .EQ. 1) THEN
+           cal(:) = 0.
+         ENDIF
+      ENDIF
+
+      ! calculate fluxes
+      CALL calcul_fluxs(knon, is_ter, dtime, &
+           tsurf, p1lay, cal, beta, tq_cdrag, tq_cdrag, pref, &
+           precip_rain, precip_snow, snow, qsurf,  &
+           radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
+           1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
+           tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+    
+      
+      !* Calculate snow height, run_off, age of snow
+
+      CALL fonte_neige( knon, is_ter, knindex, dtime, &
+           tsurf, precip_rain, precip_snow, &
+           snow, qsol, tsurf_new, evap, icesub)
+    
+      ! calculate the age of snow
+      CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
+    
+      WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
+    
+      DO i=1, knon
+         zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
+         alb_lim(i)  = alb_neig(i) *zfra(i) + alb_lim(i)*(1.0-zfra(i))
+      END DO
+
+      
+      !* Return albedo : 
+      !    alb1_new and alb2_new are here given the same values
+      
+      alb1_new(:) = 0.0
+      alb2_new(:) = 0.0
+      alb1_new(1:knon) = alb_lim(1:knon)
+      alb2_new(1:knon) = alb_lim(1:knon)
+       
+      !* Send to coupler
+      !  The run-off from river and coast are not calculated in the bucket modele.
+      !  For testing purpose of the coupled modele we put the run-off to zero.
+      IF (type_ocean=='couple') THEN
+         dummy_riverflow(:)   = 0.0
+         dummy_coastalflow(:) = 0.0
+         CALL cpl_send_land_fields(itime, knon, knindex, &
+              dummy_riverflow, dummy_coastalflow)
+      ENDIF
+
+    !
+    ! ***  Flux aggregation ***
+    !
+    ELSE IF (iflag_hetero_surf == 2) THEN
+      ! initialize output tables
+      evap_tersrf(:,:)      = 0.
+      fluxsens_tersrf(:,:)  = 0.
+      fluxlat_tersrf(:,:)   = 0.
+      tsurf_new_tersrf(:,:) = 0.
+      dflux_s_tersrf(:,:)   = 0.
+      dflux_l_tersrf(:,:)   = 0.
+      radsol_tersrf(:,:)    = 0.
+      swnet_tersrf(:,:)     = 0.
+      lwnet_tersrf(:,:)     = 0.
+      ! hyp: surface emissivity = 1
+      emis_tersrf(:,:)      = 1.
+ 
+      ! * calculate total absorbed radiance at surface
+      DO j=1, nbtersrf
+        ! SW
+        swnet_tersrf(klon_1D,j) = (1. - albedo_tersrf(klon_1D,j)) / (1. - alb_lim(klon_1D)) * swnet(klon_1D)
+        ! LW
+        ! first order
+        lwnet_tersrf(klon_1D,j) = lwnet(klon_1D) + 4. * emis_tersrf(klon_1D,j) * RSIGMA * tsurf(klon_1D)**3 * &
+                                  (tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))
+      ENDDO
+      ! LW second order corrections
+      !- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
+      IF (iflag_order2_sollw == 1) THEN
+        meansqT(:) = 0. ! as working buffer
+        !
+        DO j=1, nbtersrf
+          meansqT(klon_1D) = meansqT(klon_1D) + (tsurf_tersrf(klon_1D,j) - tsurf(klon_1D))**2 * frac_tersrf(klon_1D,j)
+        ENDDO
+        DO j=1, nbtersrf
+          lwnet_tersrf(klon_1D,j) = lwnet_tersrf(klon_1D,j) + 6. * RSIGMA * tsurf(klon_1D)**2 * (meansqT(klon_1D) - &
+                              (tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))**2)
+        ENDDO
+      ENDIF
+      ! net radiation
+      radsol_tersrf(:,:) = swnet_tersrf(:,:) + lwnet_tersrf(:,:)
+ 
+      ! * compute evapotranspiration coefficient
+      capsol(:) = 1.0/(2.5578E+06*0.15)
+      dif_grnd(:) = 0.
+
+      ! unused variables in cdrag routine
+      pblh(:) = 0.
+      ri_in(:) = 0.
+      sst(:) = 0.
+
+      ! Loop on sub-surfaces
+      DO j=1, nbtersrf
+        ! * drag coefficients
+        CALL cdrag(knon, is_ter, speed, temp_air, spechum, geop1, pref, pblh, &
+          tsurf_tersrf(:,j), qsurf_tersrf(:,j), z0m_tersrf(:,j), z0h_tersrf(:,j), ri_in, 0, &
+          cdragm_tersrf(:,j), cdragh_tersrf(:,j), zri_tersrf(:,j), plev, precip_rain, sst, p1lay)
+
+        ! * calculate temperature, heat capacity and conduction flux in soil
+        IF (soil_model) THEN
+          !
+          CALL soil_hetero(dtime, is_ter, knon, snow, tsurf_tersrf(:,j), qsol,  &
+            longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil_tersrf(:,:,j), soilcap, soilflux, &
+            inertie_tersrf(:,j), conv_ratio_tersrf(:,j))
+          !
+          cal(:) = RCPD / soilcap(:)
+          radsol_tersrf(:,j) = radsol_tersrf(:,j)  + soilflux(:)
+          !
+        ELSE
+          cal = RCPD * capsol
+          IF (klon_glo .EQ. 1) THEN
+            cal = 0.
+          ENDIF
+        ENDIF
+
+        ! * calcultaion of fluxes 
+        CALL calcul_fluxs(knon, is_ter, dtime, &
+             tsurf_tersrf(klon_1D,j), p1lay, cal, beta_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), pref, &
+             precip_rain, precip_snow, snow, qsurf_tersrf(klon_1D,j),  &
+             radsol_tersrf(klon_1D,j), dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
+             1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
+             tsurf_new_tersrf(klon_1D,j), evap_tersrf(klon_1D,j), fluxlat_tersrf(klon_1D,j), fluxsens_tersrf(klon_1D,j), &
+             dflux_s_tersrf(klon_1D,j), dflux_l_tersrf(klon_1D,j))
+
+        ! if snow > 0
+        ! calculate snow height, run_off, age of snow
+        CALL fonte_neige( knon, is_ter, knindex, dtime, &
+             tsurf_tersrf(:,j), precip_rain, precip_snow, &
+             snow, qsol, tsurf_new_tersrf(:,j), evap_tersrf(:,j), icesub)
+
+      ENDDO ! loop on sub-surfaces
+
+      ! calculate the age of snow
+      CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
+      WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
+    
+      DO i=1, knon
+        zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
+        alb_lim(i)  = alb_neig(i) *zfra(i) + alb_lim(i)*(1.0-zfra(i))
+      END DO
+
+      ! return albedo : 
+      !    alb1_new and alb2_new are here given the same values
+      alb1_new(:) = 0.0
+      alb2_new(:) = 0.0
+      alb1_new(1:knon) = alb_lim(1:knon)
+      alb2_new(1:knon) = alb_lim(1:knon)
+       
+      ! send to coupler
+      !  the run-off from river and coast are not calculated in the bucket modele.
+      !  for testing purpose of the coupled modele we put the run-off to zero.
+      IF (type_ocean=='couple') THEN
+         dummy_riverflow(:)   = 0.0
+         dummy_coastalflow(:) = 0.0
+         CALL cpl_send_land_fields(itime, knon, knindex, &
+              dummy_riverflow, dummy_coastalflow)
+      ENDIF
+
+      ! * average of fluxes and surface variables
+      qsurf = average_surf_var(klon, nbtersrf, qsurf_tersrf, frac_tersrf, 'ARI')
+      tsurf_new = average_surf_var(klon, nbtersrf, tsurf_new_tersrf, frac_tersrf, 'ARI')
+      evap = average_surf_var(klon, nbtersrf, evap_tersrf, frac_tersrf, 'ARI')
+      fluxlat = average_surf_var(klon, nbtersrf, fluxlat_tersrf, frac_tersrf, 'ARI')
+      fluxsens = average_surf_var(klon, nbtersrf, fluxsens_tersrf, frac_tersrf, 'ARI')
+      dflux_l = average_surf_var(klon, nbtersrf, dflux_l_tersrf, frac_tersrf, 'ARI')
+      dflux_s = average_surf_var(klon, nbtersrf, dflux_s_tersrf, frac_tersrf, 'ARI')
+      DO k=1, nsoilmx
+        tsoil(:,k) = average_surf_var(klon, nbtersrf, tsoil_tersrf(:,k,:), frac_tersrf, 'ARI')
+      ENDDO
+
+      ! order 2 correction to tsurf_new, for radiation computations (main atm effect of Ts)
+      IF (iflag_order2_sollw == 1) THEN
+        meansqT(:) = 0. ! as working buffer
+        DO j=1, nbtersrf
+          meansqT(klon_1D) = meansqT(klon_1D)+(tsurf_tersrf(klon_1D,j)-tsurf_new(klon_1D))**2 *frac_tersrf(klon_1D,j)
+        ENDDO
+        tsurf_new(:) = tsurf_new(:)+1.5*meansqT(:)/tsurf_new(:)
+      ENDIF   ! iflag_order2_sollw == 1
+
+    ENDIF ! iflag_hetero_surf
+!
+!* End
+!
+  END SUBROUTINE surf_land_bucket_hetero
+!
+!****************************************************************************************
+!
+END MODULE surf_land_bucket_hetero_mod
Index: /LMDZ6/trunk/libf/phylmd/surf_land_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/surf_land_mod.F90	(revision 5626)
+++ /LMDZ6/trunk/libf/phylmd/surf_land_mod.F90	(revision 5627)
@@ -20,5 +20,7 @@
        qsurf, tsurf_new, dflux_s, dflux_l, &
        flux_u1, flux_v1 , & 
-       veget,lai,height &
+       veget,lai,height, tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
+       cdragm_tersrf, cdragh_tersrf, &
+       swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf &
 #ifdef ISO
        ,xtprecip_rain, xtprecip_snow,xtspechum, &
@@ -63,4 +65,5 @@
     
     USE surf_land_bucket_mod
+    USE surf_land_bucket_hetero_mod
     USE calcul_fluxs_mod
     USE indice_sol_mod
@@ -78,5 +81,5 @@
 USE print_control_mod, ONLY: lunout
     USE dimsoil_mod_h, ONLY: nsoilmx
-
+    USE compbl_mod_h
 
 ! Input variables  
@@ -89,5 +92,5 @@
     LOGICAL, INTENT(IN)                     :: debut, lafin
     REAL, INTENT(IN)                        :: dtime
-    REAL, DIMENSION(klon), INTENT(IN)       :: zlev, ccanopy
+    REAL, DIMENSION(klon), INTENT(IN)       :: ccanopy
     REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
     REAL, DIMENSION(klon), INTENT(IN)       :: albedo  ! albedo for whole short-wave interval
@@ -106,4 +109,5 @@
                                                          ! corresponds to previous sollwdown
     REAL, DIMENSION(klon), INTENT(IN)       :: q2m, t2m
+    REAL, DIMENSION(klon, nbtersrf), INTENT(IN) :: tsurf_tersrf
 #ifdef ISO
     REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow
@@ -115,4 +119,6 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: zlev
+    REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf
 #ifdef ISO
     REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: xtsnow, xtsol
@@ -136,4 +142,13 @@
     REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: veget,lai
     REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: height
+! AM
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: tsurf_new_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: qsurf_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragm_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragh_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: swnet_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: lwnet_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxsens_tersrf
+    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxlat_tersrf
 #ifdef ISO
     REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
@@ -153,5 +168,7 @@
     REAL, DIMENSION(klon) :: u0, v0     ! surface speed
     REAL, DIMENSION(klon) :: precip_totsnow     ! total solid precip
-    INTEGER               :: i
+    INTEGER               :: i,j
+    CHARACTER (len = 20)  :: modname = 'surf_land'
+    CHARACTER (len = 100) :: abort_message
 
 !albedo SB >>>
@@ -285,4 +302,23 @@
         !write(*,*) 'surf_land 258'
 #endif
+      IF (iflag_hetero_surf .GT. 0) THEN
+        IF (klon .EQ. 1) THEN
+          !
+          CALL surf_land_bucket_hetero(itime, jour, knon, knindex, debut, dtime,&
+              tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
+              spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, pref, zlev, &
+              u1, v1, gustiness, rugoro, swnet, lwnet, &
+              snow, qsol, agesno, tsoil, &
+              qsurf, z0m, z0h, alb1_new, alb2_new, evap, &
+              fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l, &
+              tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
+              cdragm_tersrf, cdragh_tersrf, &
+              swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)
+        ELSE
+          abort_message = 'Heterogeneous continental subsurfaces (iflag_hetero_surf > 0) are only compatible in 1D cases.'
+          CALL abort_physic(modname,abort_message,1)
+        ENDIF
+      !
+      ELSE
        CALL surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
             tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
@@ -300,4 +336,5 @@
         z0h(1:knon)=z0m(1:knon) ! En attendant mieux
 
+      ENDIF ! iflag_hetero_surf
 
     ENDIF ! ok_veget
Index: /LMDZ6/trunk/libf/phylmd/surf_param_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/surf_param_mod.F90	(revision 5627)
+++ /LMDZ6/trunk/libf/phylmd/surf_param_mod.F90	(revision 5627)
@@ -0,0 +1,168 @@
+MODULE surf_param_mod
+!
+  IMPLICIT NONE
+!
+CONTAINS
+!
+!-------------------------------------------------------------------------------
+!
+FUNCTION eff_surf_param(klon, nbtersurf, x, frac, hatype, Zref) RESULT(eff_param)
+!
+!-------------------------------------------------------------------------------
+!
+! Arguments:
+  INTEGER, INTENT(IN)                          :: klon       ! grid point
+  INTEGER, INTENT(IN)                          :: nbtersurf  ! number of land hetero. subsurfaces
+  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: x          ! variable or parameter to integrate
+  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: frac       ! fraction of each land hetero. subsurface
+  CHARACTER(LEN=3), INTENT(IN)                 :: hatype     ! method to integrate the parameter
+  REAL, OPTIONAL, DIMENSION(klon), INTENT(IN)  :: Zref       ! reference height for CDN averaging (m)
+  REAL, DIMENSION(klon)                        :: eff_param  ! effective parameter
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER ik, is
+  REAL :: zrefd = 10.  ! default reference height for CDN averaging (m)
+  REAL :: Cdref        ! reference height for CDN averaging (m)
+!-------------------------------------------------------------------------------
+!
+  eff_param(:) = 0.
+  DO ik = 1, klon
+    DO is = 1, nbtersurf
+      !
+      ! arithmetic averaging
+      IF (hatype == 'ARI') THEN
+        eff_param(ik) = eff_param(ik) + frac(ik,is) * x(ik,is)
+      !
+      ! inverse averaging
+      ELSEIF (hatype == 'INV') THEN
+        IF (x(ik,is) .NE. 0.) THEN
+          eff_param(ik) = eff_param(ik) + frac(ik,is) * 1./x(ik,is)
+        ENDIF
+      !
+      ! inverse of square logarithm averaging
+      ELSEIF (hatype == 'CDN') THEN
+        IF (PRESENT(Zref)) THEN
+          Cdref = Zref(ik)
+        ELSE
+          Cdref = zrefd
+        ENDIF
+      !
+        IF (x(ik,is) .NE. 0.) THEN
+          eff_param(ik) = eff_param(ik) + frac(ik,is) * 1./(LOG(Cdref/x(ik,is)))**2
+        ENDIF
+      !
+      ELSE
+        PRINT*, 'eff_surf_param: invalid averaging type: ', hatype
+      ENDIF
+    ENDDO
+    !
+    IF (hatype == 'CDN') THEN
+      eff_param(ik) = Cdref * exp(-sqrt(1./eff_param(ik)))
+    ENDIF
+  !
+  ENDDO
+!
+END FUNCTION eff_surf_param
+!
+!
+!-------------------------------------------------------------------------------
+!
+FUNCTION average_surf_var(klon, nbtersurf, x, frac, hatype) RESULT(x_avg)
+!
+!-------------------------------------------------------------------------------
+!
+! Arguments:
+  INTEGER, INTENT(IN)                          :: klon       ! grid point
+  INTEGER, INTENT(IN)                          :: nbtersurf  ! number of land hetero. subsurfaces
+  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: x          ! variable or parameter to integrate
+  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: frac       ! fraction of each land hetero. subsurface
+  CHARACTER(LEN=3), INTENT(IN)                 :: hatype     ! method to integrate the parameter
+  REAL, DIMENSION(klon)                        :: x_avg      ! average variable
+!
+! Local variables:
+  INTEGER ik, is
+!
+!-------------------------------------------------------------------------------
+!
+  x_avg(:) = 0.
+  DO ik = 1, klon
+    DO is = 1, nbtersurf
+      !
+      ! arithmetic averaging
+      IF (hatype == 'ARI') THEN
+        x_avg(ik) = x_avg(ik) + frac(ik,is) * x(ik,is)
+      ELSE
+        PRINT*, 'average_surf_var: invalid averaging type: ', hatype
+      ENDIF
+    ENDDO
+  ENDDO
+!
+END FUNCTION average_surf_var
+!
+!-------------------------------------------------------------------------------
+!
+FUNCTION interpol_tsoil(klon, nbtersurf, nsoilmx, nbtsoildepths, alpha, period, inertie, hcond, tsoil_depth, tsurf, tsoil_) RESULT(tsoil)
+!
+!-------------------------------------------------------------------------------
+!
+! Arguments:
+  INTEGER, INTENT(IN)                                         :: klon          ! grid point
+  INTEGER, INTENT(IN)                                         :: nbtersurf     ! number of land hetero. subsurfaces
+  INTEGER, INTENT(IN)                                         :: nsoilmx       ! number of soil layers in the model grid
+  INTEGER, INTENT(IN)                                         :: nbtsoildepths ! number of soil depths for soil temperature initialization
+  REAL, INTENT(IN)                                            :: alpha         ! parameter for soil discretization
+  REAL, INTENT(IN)                                            :: period        ! parameter for soil discretization
+  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: inertie       ! soil thermal inertia
+  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: hcond         ! soil heat conductivity
+  REAL, DIMENSION(klon, nbtsoildepths, nbtersurf), INTENT(IN) :: tsoil_depth   ! soil depth at which temperature is given (m)
+  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: tsurf         ! surface temperature
+  REAL, DIMENSION(klon, nbtsoildepths, nbtersurf), INTENT(IN) :: tsoil_        ! soil temperature given at tsoil_depths
+  REAL, DIMENSION(klon, nsoilmx, nbtersurf)                   :: tsoil         ! soil temperature interpolated in the model grid
+!
+! Local variables:
+  INTEGER ik, is, iq, it
+  REAL pi, slope, inter
+  REAL, DIMENSION(klon, nbtersurf)          :: z1, hcap  ! first layer depth and soil heat capacity
+  REAL, DIMENSION(klon, nsoilmx, nbtersurf) :: z         ! depth of the middle of the soil layer in the model grid (m)
+!
+!-------------------------------------------------------------------------------
+!
+  pi = ACOS(-1.)
+  !
+  DO ik = 1, klon
+    DO is = 1, nbtersurf
+      !
+      hcap(ik,is) = inertie(ik,is)*inertie(ik,is)/hcond(ik,is)
+      z1(ik,is) = SQRT(period*hcond(ik,is)/(pi*hcap(ik,is)))
+      !
+      DO iq = 1, nsoilmx
+        ! compute depth of middle soil layer (in m)
+        z(ik,iq,is) = z1(ik,is) * (alpha**(iq-0.5) - 1.) / (alpha - 1.)
+        ! if z is between the surface and first tsoil_depth
+        IF ((z(ik,iq,is) .GT. 0.) .AND. (z(ik,iq,is) .LE. tsoil_depth(ik,1,is))) THEN
+          slope = (tsoil_(ik,1,is) - tsurf(ik,is)) / (tsoil_depth(ik,1,is) - 0.)
+          inter = tsurf(ik,is)
+          tsoil(ik,iq,is) = slope * z(ik,iq,is) + inter
+        ENDIF
+        ! other levels
+        DO it = 1, nbtsoildepths-1
+          IF ((z(ik,iq,is) .GT. tsoil_depth(ik,it,is)) .AND. (z(ik,iq,is) .LE. tsoil_depth(ik,it+1,is))) THEN
+            slope = (tsoil_(ik,it+1,is) - tsoil_(ik,it,is)) / (tsoil_depth(ik,it+1,is) - tsoil_depth(ik,it,is))
+            inter = tsoil_(ik,it,is) - slope * tsoil_depth(ik,it,is)
+            tsoil(ik,iq,is) = slope * z(ik,iq,is) + inter
+          ENDIF
+        ENDDO
+        ! for layers below
+        IF (z(ik,iq,is) .GT. tsoil_depth(ik,nbtsoildepths,is)) THEN
+          tsoil(ik,iq,is) = tsoil_(ik,nbtsoildepths,is)
+        ENDIF
+      ENDDO
+      !
+    ENDDO
+  ENDDO
+!
+END FUNCTION interpol_tsoil
+!
+!-------------------------------------------------------------------------------
+!
+END MODULE surf_param_mod
