Changeset 3900

Timestamp:

May 17, 2021, 3:35:58 PM (3 years ago)

Author:

evignon

Message:

Commission de la nouvelle interface LMDZ-SISVAT
la prochaine commission consistera a supprimer l'ancien repertoire sisvat
et a faire un peu de nettoyage.

Location:

LMDZ6/trunk/libf/phylmd

Files:

• conf_phys_m.F90 (modified) (15 diffs)
• create_etat0_unstruct.F90 (modified) (1 diff)
• dimsoil.h (modified) (1 diff)
• fonte_neige_mod.F90 (modified) (1 diff)
• inlandsis/VARphy.F90 (modified) (2 diffs)
• inlandsis/VARtSV.F90 (modified) (2 diffs)
• inlandsis/VARxSV.F90 (modified) (6 diffs)
• inlandsis/VARySV.F90 (modified) (3 diffs)
• inlandsis/inlandsis.F (modified) (34 diffs)
• inlandsis/sisvat_bsn.F (modified) (1 diff)
• inlandsis/sisvat_qsn.F (modified) (6 diffs)
• inlandsis/sisvat_sop.F (modified) (18 diffs)
• inlandsis/sisvat_tso.F (modified) (11 diffs)
• inlandsis/sisvat_zsn.F (modified) (2 diffs)
• inlandsis/surf_inlandsis_mod.F90 (modified) (1 diff)
• pbl_surface_mod.F90 (modified) (14 diffs)
• phys_output_write_mod.F90 (modified) (1 diff)
• physiq_mod.F90 (modified) (1 diff)
• surf_landice_mod.F90 (modified) (10 diffs)
• surface_data.F90 (modified) (1 diff)

LMDZ6/trunk/libf/phylmd/conf_phys_m.F90

@@ -27 +27 @@
     USE phys_cal_mod
     USE carbon_cycle_mod,  ONLY: carbon_cycle_tr, carbon_cycle_cpl, carbon_cycle_rad, level_coupling_esm
-    USE carbon_cycle_mod,  ONLY: read_fco2_ocean_cor,var_fco2_ocean_cor
-    USE carbon_cycle_mod,  ONLY: read_fco2_land_cor,var_fco2_land_cor
     USE mod_grid_phy_lmdz, ONLY: klon_glo
     USE print_control_mod, ONLY: lunout
@@ -90 +88 @@
     CHARACTER (len = 8), SAVE  :: aer_type_omp
     INTEGER, SAVE       :: landice_opt_omp
-    INTEGER, SAVE       :: n_dtis_omp
-    INTEGER, SAVE       :: iflag_tsurf_inlandsis_omp
-    INTEGER, SAVE       :: iflag_albzenith_omp    
-    LOGICAL, SAVE       :: SnoMod_omp,BloMod_omp,ok_outfor_omp
+    INTEGER, SAVE       :: iflag_tsurf_inlandsis_omp,iflag_temp_inlandsis_omp
+    INTEGER, SAVE       :: iflag_albcalc_omp,iflag_z0m_snow_omp   
+    LOGICAL, SAVE       :: SnoMod_omp,BloMod_omp,ok_outfor_omp,ok_zsn_ii_omp
+    LOGICAL, SAVE       :: discret_xf_omp,opt_runoff_ac_omp  
+    LOGICAL, SAVE       :: is_ok_slush_omp,is_ok_z0h_rn_omp,is_ok_density_kotlyakov_omp
+    REAL, SAVE          :: prescribed_z0m_snow_omp,correc_alb_omp
+    REAL, SAVE          :: buf_sph_pol_omp,buf_siz_pol_omp
     LOGICAL, SAVE       :: ok_newmicro_omp
     LOGICAL, SAVE       :: ok_all_xml_omp
@@ -239 +240 @@
     LOGICAL, SAVE :: carbon_cycle_rad_omp
     INTEGER, SAVE :: level_coupling_esm_omp
-    LOGICAL, SAVE :: read_fco2_ocean_cor_omp
-    REAL, SAVE    :: var_fco2_ocean_cor_omp
-    LOGICAL, SAVE :: read_fco2_land_cor_omp
-    REAL, SAVE    :: var_fco2_land_cor_omp
     LOGICAL, SAVE :: adjust_tropopause_omp
     LOGICAL, SAVE :: ok_daily_climoz_omp
@@ -340 +337 @@
 
     !Etienne
+    !Config Key  = iflag_temp_inlandsis
+    !Config Desc = which method to calculate temp within the soil in INLANDSIS
+    !Config Def  = 0
+    iflag_temp_inlandsis_omp = 0
+    CALL getin('iflag_temp_inlandsis', iflag_temp_inlandsis_omp)
+
+    !Etienne
     !Config Key  = iflag_tsurf_inlandsis
     !Config Desc = which method to calculate tsurf in INLANDSIS
     !Config Def  = 0
-    iflag_tsurf_inlandsis_omp = 0
+    iflag_tsurf_inlandsis_omp = 1
     CALL getin('iflag_tsurf_inlandsis', iflag_tsurf_inlandsis_omp)
 
+
     !Etienne
-    !Config Key  = iflag_albzenith
-    !Config Desc = method to account for albedo sensitivity to solar zenith angle
-    !Config Def  = 0
-    iflag_albzenith_omp = 0
-    CALL getin('iflag_albzenith', iflag_albzenith_omp)
-
-    !Etienne
-    !Config Key  = n_dtis
-    !Config Desc = number of subtimesteps for INLANDSIS
-    !Config Def  = 1
-    n_dtis_omp = 1
-    CALL getin('n_dtis', n_dtis_omp)
+    !Config Key  = iflag_albcalc
+    !Config Desc = method to calculate snow albedo in INLANDSIS
+    !Config Def  = 0
+    iflag_albcalc_omp = 0
+    CALL getin('iflag_albcalc', iflag_albcalc_omp)
+
 
     !Etienne
     !Config Key  = SnoMod
     !Config Desc = activation of snow modules in inlandsis
-    !Config Def  = 1
+    !Config Def  = .TRUE.
     SnoMod_omp = .TRUE.
     CALL getin('SnoMod', SnoMod_omp)
@@ -370 +369 @@
     !Config Key  = BloMod
     !Config Desc = activation of blowing snow in inlandsis
-    !Config Def  = 1
+    !Config Def  = .FALSE.
     BloMod_omp = .FALSE.
     CALL getin('BloMod', BloMod_omp)
@@ -377 +376 @@
     !Config Key  = ok_outfor
     !Config Desc = activation of output ascii file in inlandsis
-    !Config Def  = 1
+    !Config Def  = .FALSE.
     ok_outfor_omp = .FALSE.
     CALL getin('ok_outfor', ok_outfor_omp)
+
+
+    !Etienne
+    !Config Key  = ok_sn_ii
+    !Config Desc = activation of ice/snow detection
+    !Config Def  = .TRUE.
+    ok_zsn_ii = .TRUE.
+    CALL getin('ok_zsn_ii', ok_zsn_ii_omp)
+
+
+    !Etienne
+    !Config Key  = discret_xf
+    !Config Desc = snow discretization following XF
+    !Config Def  = .TRUE.
+    discret_xf = .TRUE.
+    CALL getin('discret_xf', discret_xf_omp)
+
+
+    !Etienne
+    !Config Key  = is_ok_slush
+    !Config Desc = activation of the slush option
+    !Config Def  = .TRUE.
+    is_ok_slush_omp = .TRUE.
+    CALL getin('is_ok_slush', is_ok_slush_omp)
+
+    !Etienne
+    !Config Key  = opt_runoff_ac
+    !Config Desc = option runoff AC
+    !Config Def  = .TRUE.
+    opt_runoff_ac_omp = .TRUE.
+    CALL getin('opt_runoff_ac', opt_runoff_ac_omp)
+
+    !Etienne
+    !Config Key  = is_ok_z0h_rn
+    !Config Desc = z0h calculation following RN method
+    !Config Def  = .TRUE.
+    is_ok_z0h_rn = .TRUE.
+    CALL getin('is_ok_z0h_rn', is_ok_z0h_rn_omp)
+
+
+    !Etienne
+    !Config Key  = is_ok_density_kotlyakov
+    !Config Desc = snow density calculation following kotlyakov
+    !Config Def  = .FALSE.
+    is_ok_density_kotlyakov = .FALSE.
+    CALL getin('is_ok_density_kotlyakov', is_ok_density_kotlyakov_omp)
+
+
+    !Etienne
+    !Config Key  = prescribed_z0m_snow
+    !Config Desc = prescribed snow z0m
+    !Config Def  = 0.005
+    prescribed_z0m_snow = 0.005
+    CALL getin('prescribed_z0m_snow', prescribed_z0m_snow_omp)
+
+
+    !Etienne
+    !Config Key  = iflag_z0m_snow
+    !Config Desc = method to calculate snow z0m
+    !Config Def  = 0
+    iflag_z0m_snow = 0
+    CALL getin('iflag_z0m_snow', iflag_z0m_snow_omp)
+
+
+    !Etienne
+    !Config Key  = correc_alb
+    !Config Desc = correction term for albedo
+    !Config Def  = 1.01
+    correc_alb_omp=1.01
+    CALL getin('correc_alb', correc_alb_omp)
+
+
+    !Etienne
+    !Config Key  = buf_sph_pol
+    !Config Desc = sphericity of buffer layer in polar regions
+    !Config Def  = 99.
+    buf_sph_pol_omp=99.
+    CALL getin('buf_sph_pol', buf_sph_pol_omp)
+
+    !Etienne
+    !Config Key  = buf_siz_pol
+    !Config Desc = grain size of buffer layer in polar regions in e-4m
+    !Config Def  = 4.
+    buf_siz_pol_omp=4.
+    CALL getin('buf_siz_pol', buf_siz_pol_omp)
 
     !==================================================================
@@ -1140 +1224 @@
     CALL getin('coefw_cld_cv',coefw_cld_cv_omp)
 
+
+
+
     !
     !Config Key  = iflag_pdf 
@@ -1148 +1235 @@
     iflag_pdf_omp = 0
     CALL getin('iflag_pdf',iflag_pdf_omp)
-
     !
     !Config Key  = fact_cldcon
@@ -1175 +1261 @@
     ok_newmicro_omp = .TRUE.
     CALL getin('ok_newmicro',ok_newmicro_omp)
-
     !
     !Config Key  = ratqsbas
@@ -1184 +1269 @@
     ratqsbas_omp = 0.01
     CALL getin('ratqsbas',ratqsbas_omp)
-
     !
     !Config Key  = ratqshaut
@@ -2218 +2302 @@
     CALL getin('carbon_cycle_rad',carbon_cycle_rad_omp)
 
-    read_fco2_ocean_cor_omp=.FALSE.
-    CALL getin('read_fco2_ocean_cor',read_fco2_ocean_cor_omp)
-
-    var_fco2_ocean_cor_omp=0. ! default value
-    CALL getin('var_fco2_ocean_cor',var_fco2_ocean_cor_omp)
-
-    read_fco2_land_cor_omp=.FALSE.
-    CALL getin('read_fco2_land_cor',read_fco2_land_cor_omp)
-
-    var_fco2_land_cor_omp=0. ! default value
-    CALL getin('var_fco2_land_cor',var_fco2_land_cor_omp)
-
+    ! >> PC
     ! level_coupling_esm : level of coupling of the biogeochemical fields between LMDZ, ORCHIDEE and NEMO
     ! Definitions of level_coupling_esm in physiq.def
@@ -2242 +2315 @@
     level_coupling_esm_omp=0 ! default value
     CALL getin('level_coupling_esm',level_coupling_esm_omp)
+    ! << PC
 
     !$OMP END MASTER
@@ -2358 +2432 @@
        ok_veget=.FALSE.
     ENDIF
-    ! SISVAT and INLANDSIS
+    ! INLANDSIS
     !=================================================
     landice_opt = landice_opt_omp
     iflag_tsurf_inlandsis = iflag_tsurf_inlandsis_omp
-    iflag_albzenith = iflag_albzenith_omp
-    n_dtis=n_dtis_omp
+    iflag_temp_inlandsis = iflag_temp_inlandsis_omp
+    iflag_albcalc = iflag_albcalc_omp
     SnoMod=SnoMod_omp
     BloMod=BloMod_omp
     ok_outfor=ok_outfor_omp
+    is_ok_slush=is_ok_slush_omp
+    opt_runoff_ac=opt_runoff_ac_omp
+    is_ok_z0h_rn=is_ok_z0h_rn_omp
+    is_ok_density_kotlyakov=is_ok_density_kotlyakov_omp
+    prescribed_z0m_snow=prescribed_z0m_snow_omp
+    correc_alb=correc_alb_omp
+    iflag_z0m_snow=iflag_z0m_snow_omp
+    ok_zsn_ii=ok_zsn_ii_omp
+    discret_xf=discret_xf_omp
+    buf_sph_pol=buf_sph_pol_omp
+    buf_siz_pol=buf_siz_pol_omp
     !=================================================
     ok_all_xml = ok_all_xml_omp
@@ -2507 +2592 @@
     carbon_cycle_rad = carbon_cycle_rad_omp
     level_coupling_esm = level_coupling_esm_omp
-    read_fco2_ocean_cor = read_fco2_ocean_cor_omp
-    var_fco2_ocean_cor = var_fco2_ocean_cor_omp
-    read_fco2_land_cor = read_fco2_land_cor_omp
-    var_fco2_land_cor = var_fco2_land_cor_omp
 
     ! Test of coherence between type_ocean and version_ocean
@@ -2831 +2912 @@
     WRITE(lunout,*) ' level_coupling_esm = ', level_coupling_esm
     WRITE(lunout,*) ' iflag_tsurf_inlandsis = ', iflag_tsurf_inlandsis
-    WRITE(lunout,*) ' iflag_albzenith = ', iflag_albzenith
-    WRITE(lunout,*) ' n_dtis = ', n_dtis
+    WRITE(lunout,*) ' iflag_temp_inlandsis = ', iflag_temp_inlandsis
+    WRITE(lunout,*) ' iflag_albcalc = ', iflag_albcalc
     WRITE(lunout,*) ' SnoMod = ', SnoMod
     WRITE(lunout,*) ' BloMod = ', BloMod
     WRITE(lunout,*) ' ok_outfor = ', ok_outfor
-
+    WRITE(lunout,*) ' is_ok_slush = ', is_ok_slush
+    WRITE(lunout,*) ' opt_runoff_ac = ', opt_runoff_ac
+    WRITE(lunout,*) ' is_ok_z0h_rn = ', is_ok_z0h_rn
+    WRITE(lunout,*) ' is_ok_density_kotlyakov = ', is_ok_density_kotlyakov
+    WRITE(lunout,*) ' prescribed_z0m_snow = ', prescribed_z0m_snow
+    WRITE(lunout,*) ' iflag_z0m_snow = ', iflag_z0m_snow
+    WRITE(lunout,*) ' ok_zsn_ii = ', ok_zsn_ii
+    WRITE(lunout,*) ' discret_xf = ', discret_xf
+    WRITE(lunout,*) ' correc_alb= ', correc_alb
+    WRITE(lunout,*) ' buf_sph_pol = ', buf_sph_pol
+    WRITE(lunout,*) ' buf_siz_pol= ', buf_siz_pol
 
     !$OMP END MASTER

LMDZ6/trunk/libf/phylmd/create_etat0_unstruct.F90

@@ -209 +209 @@
     z0m(:,is_oce) = rugmer(:)
 
-   z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-   z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
+   z0m(:,is_ter) = 0.01 ! MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
+   z0m(:,is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
 
    z0m(:,is_sic) = 0.001

LMDZ6/trunk/libf/phylmd/dimsoil.h

@@ -8 +8 @@
 
       INTEGER nsnowmx
-      PARAMETER (nsnowmx=35)
+      PARAMETER (nsnowmx=30)
       
       INTEGER nsismx
-      PARAMETER (nsismx=46) 
+      PARAMETER (nsismx=41) 
 
 ! nsismx should be equal to nsoilmx+nsnowmx 

LMDZ6/trunk/libf/phylmd/fonte_neige_mod.F90

@@ -28 +28 @@
   REAL, PRIVATE                               :: tau_calv  
   !$OMP THREADPRIVATE(tau_calv)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: ffonte_global
+  REAL, ALLOCATABLE, DIMENSION(:,:)           :: ffonte_global
   !$OMP THREADPRIVATE(ffonte_global)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: fqfonte_global
+  REAL, ALLOCATABLE, DIMENSION(:,:)           :: fqfonte_global
   !$OMP THREADPRIVATE(fqfonte_global)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: fqcalving_global
+  REAL, ALLOCATABLE, DIMENSION(:,:)           :: fqcalving_global
   !$OMP THREADPRIVATE(fqcalving_global)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE  :: runofflic_global
+  REAL, ALLOCATABLE, DIMENSION(:)             :: runofflic_global
   !$OMP THREADPRIVATE(runofflic_global)
 

LMDZ6/trunk/libf/phylmd/inlandsis/VARphy.F90

@@ -26 +26 @@
       INTEGER, PARAMETER ::  iun=1                                             
       REAL, PARAMETER    ::  zer0 = 0.0e+0, half = 0.5e+0, un_1 = 1.0e+0,     &
-     &                       eps6 = 1.0e-6, R_1000=1.e3       
+     &                       eps6 = 1.0e-6, R_1000=1.e3   
       REAL, PARAMETER    ::  zero = 0.0e+0, demi = 0.5e+0, unun = 1.0e+0,     &
      &                       epsi = 1.0e-6, eps9 = 1.0e-9         
@@ -91 +91 @@
 ! A1.6 Turbulent and molecular diffusion
 !----------------------------------------
-      REAL, PARAMETER    ::  A_MolV = 1.35e-5, vonKrm = 0.40e0
+      REAL, PARAMETER    ::  A_MolV = 1.35e-5, vonKrm = 0.40e0, r_turb=3.0
+      REAL, PARAMETER    ::  A_turb=5.8, akmol=1.35e-5
 !C +...                A_MolV: Air Viscosity                 = 1.35d-5 m2/s   
 !C +                   vonKrm: von Karman constant           = 0.4            
-                                                                          
+!C +                   r_turb:   Turbulent Diffusivities Ratio K*/Km        
+!C +                   A_turb:    Stability  Coefficient Moment                                 
+!C +                   Air Viscosity                 = 1.35d-5 m2/s                 
+
+
 
 END MODULE VARphy

LMDZ6/trunk/libf/phylmd/inlandsis/VARtSV.F90

@@ -41 +41 @@
 
   SUBROUTINE INIT_VARtSV
+  
   IMPLICIT NONE
-  
+
+  INTEGER ikl
+
+
+
+
+
+
+
       ALLOCATE(toicSV(klonv))
 
@@ -59 +68 @@
       ALLOCATE(rsolSV(klonv))                ! Radiation balance surface
 
+      DO ikl=1,klonv        
+           
+         toicSV(ikl)   = 0.
+         dz1_SV(ikl,:) = 0.
+         dz2_SV(ikl,:) = 0.
+         Tsf_SV(ikl)   = 0.
+         TsfnSV(ikl)   = 0.
+         AcoHSV(ikl)   = 0.
+         BcoHSV(ikl)   = 0.
+         AcoQSV(ikl)   = 0.
+         ps__SV(ikl)   = 0.
+         p1l_SV(ikl)   = 0.
+         cdH_SV(ikl)   = 0.
+         cdM_SV(ikl)   = 0.
+         rsolSV(ikl)   = 0.
+      END DO
+
+
+
   END SUBROUTINE INIT_VARtSV
 

LMDZ6/trunk/libf/phylmd/inlandsis/VARxSV.F90

@@ -67 +67 @@
       REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QaT_SV  ! SBL Top   Specific Humidity     
 !$OMP THREADPRIVATE(QaT_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QsT_SV  ! SBL Top   Specific Humidity
+!$OMP THREADPRIVATE(QsT_SV)
       REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dQa_SV  ! SBL Flux  Limitation of Qa      
 !$OMP THREADPRIVATE(dQa_SV)
@@ -78 +80 @@
 
                                                               
-      REAL,SAVE                      ::   zSBLSV  ! SBL Height (Initial Value)      
-!$OMP THREADPRIVATE(zSBLSV)
       REAL,SAVE                      ::   dt__SV  ! Time Step                       
 !$OMP THREADPRIVATE(dt__SV)
@@ -160 +160 @@
       REAL,ALLOCATABLE,SAVE    ::   agsnSV(:,:)  ! Snow Age                        
 !$OMP THREADPRIVATE(agsnSV)
+      REAL,ALLOCATABLE,SAVE    ::   DOPsnSV(:,:)  ! Snow optical diameter [m]
+!$OMP THREADPRIVATE(DOPsnSV)
       REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BufsSV  ! Snow Buffer Layer               
 !$OMP THREADPRIVATE(BufsSV)
@@ -260 +262 @@
       ALLOCATE(dLdTSV(klonv))  ! Latent   Heat Flux T Derivat.   
       ALLOCATE(rhT_SV(klonv))  ! SBL Top   Air  Density          
-      ALLOCATE(QaT_SV(klonv))  ! SBL Top   Specific Humidity     
+      ALLOCATE(QaT_SV(klonv))  ! SBL Top   Specific Humidity    
+      ALLOCATE(QsT_SV(klonv))  ! surface   Specific Humidity 
       ALLOCATE(dQa_SV(klonv))  ! SBL Flux  Limitation of Qa      
       ALLOCATE(qsnoSV(klonv))  ! SBL Mean  Snow       Content    
@@ -309 +312 @@
       ALLOCATE(dzsnSV(klonv,    0:nsno))  ! Snow Layer  Thickness           
       ALLOCATE(agsnSV(klonv,    0:nsno))  ! Snow Age                        
+      ALLOCATE(DOPsnSV(klonv,    0:nsno))  ! Snow Optical diameter                        
       ALLOCATE(BufsSV(klonv))  ! Snow Buffer Layer               
       ALLOCATE(rusnSV(klonv))  ! Surficial   Water               
@@ -339 +343 @@
 
       DO ikl=1,klonv        
-        LSmask(ikl)   = 0
-        isotSV(ikl)   = 0                
-        iWaFSV(ikl)   = 0 
-        isnoSV(ikl)   = 0  
-        ispiSV(ikl)   = 0
-        iiceSV(ikl)   = 0          
-        istoSV(ikl,:) = 0
-        ii__SV(ikl)   = 0
-        jj__SV(ikl)   = 0
-        nn__SV(ikl)   = 0 
+
+
+      isnoSV(ikl)  =0.        
+      ispiSV(ikl)  =0.   
+      iiceSV(ikl)  =0.        
+      istoSV(ikl,:)=0.                                                                        
+      alb_SV(ikl)  =0.      
+      emi_SV(ikl)  =0.    
+      IRs_SV(ikl)  =0.      
+      LMO_SV(ikl)  =0.
+      us__SV(ikl)  =0.
+      uts_SV(ikl)  =0.
+      cutsSV(ikl)  =0.
+      uqs_SV(ikl)  =0.
+      uss_SV(ikl)  =0.
+      usthSV(ikl)  =0.
+      rCDmSV(ikl)  =0.
+      rCDhSV(ikl)  =0.
+      Z0m_SV(ikl)  =0.
+      Z0mmSV(ikl)  =0.
+      Z0mnSV(ikl)  =0.
+      Z0roSV(ikl)  =0.
+      Z0SaSV(ikl)  =0.
+      Z0e_SV(ikl)  =0.
+      Z0emSV(ikl)  =0.
+      Z0enSV(ikl)  =0.
+      Z0h_SV(ikl)  =0.
+      Z0hmSV(ikl)  =0.
+      Z0hnSV(ikl)  =0.
+                                                                          
+                                                                        
+      TsisSV(ikl,:)  =0.
+      ro__SV(ikl,:)  =0.
+      eta_SV(ikl,:)  =0.
+      G1snSV(ikl,:)  =0.  
+      G2snSV(ikl,:)  =0.          
+      dzsnSV(ikl,:)  =0.      
+      agsnSV(ikl,:)  =0.              
+      DOPsnSV(ikl,:) =0.                  
+      BufsSV(ikl)  =0.         
+      rusnSV(ikl)  =0.     
+      SWf_SV(ikl)  =0.        
+      SWS_SV(ikl)  =0.
+      HFraSV(ikl)  =0.       
+                                                                            
+      zWE_SV(ikl)  =0.
+      zWEcSV(ikl)  =0.
+      wem_SV(ikl)  =0.
+      wer_SV(ikl)  =0.
+      wes_SV(ikl)  =0.
+      zn4_SV(ikl)  =0.
+      zn5_SV(ikl)  =0.                                                        
+                                       
+                                                                                                                                                          
+      ii__SV(ikl)  =0. 
+      jj__SV(ikl)  =0.
+      nn__SV(ikl)  =0.
+                                                                               
+      IRu_SV(ikl)  =0. 
+      hSalSV(ikl)  =0.     
+      qSalSV(ikl)  =0.
+      RnofSV(ikl)  =0.    
+      RuofSV(ikl,:)  =0.
+
+
+
+
       END DO
   END SUBROUTINE INIT_VARxSV

LMDZ6/trunk/libf/phylmd/inlandsis/VARySV.F90

@@ -22 +22 @@
       REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb3sv  ! Surface Albedo FIR      
 !$OMP THREADPRIVATE(alb3sv)
-
+      REAL, DIMENSION(:,:),SAVE,ALLOCATABLE  ::   alb6sv  ! 6 band-albedo    
+!$OMP THREADPRIVATE(alb6sv)
       REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   albssv  ! Soil               Albedo [-]   
 !$OMP THREADPRIVATE(albssv)
@@ -83 +84 @@
       ALLOCATE(alb2sv(klonv))  ! Surface Albedo NIR      
       ALLOCATE(alb3sv(klonv))  ! Surface Albedo FIR       
+      ALLOCATE(alb6sv(klonv,6))! 6-band  Albedo      
 
       !
@@ -110 +112 @@
 
       DO ikl=1,klonv
-        NLaysv(ikl)   = 0
-        i_thin(ikl)   = 0
-        LIndsv(ikl)   = 0
+
+      NLaysv(ikl) =0.
+      i_thin(ikl) =0.
+      LIndsv(ikl) =0.
+      albisv(ikl) =0.
+      alb1sv(ikl) =0. 
+      alb2sv(ikl) =0.   
+      alb3sv(ikl) =0.
+      alb6sv(ikl,:)=0.
+      albssv(ikl) =0.
+      SoSosv(ikl) =0.
+      Eso_sv(ikl) =0.
+      HSv_sv(ikl) =0.   
+      HLv_sv(ikl) =0.
+      HSs_sv(ikl) =0.       
+      HLs_sv(ikl) =0.
+      sqrCm0(ikl) =0.   
+      sqrCh0(ikl) =0.    
+      Lx_H2O(ikl) =0.
+      ram_sv(ikl) =0.
+      rah_sv(ikl) =0.
+      Fh__sv(ikl) =0.      
+      dFh_sv(ikl) =0.
+      Evp_sv(ikl) =0.
+      EvT_sv(ikl) =0.
+      LSdzsv(ikl) =0.
+      Tsrfsv(ikl) =0.
+      sEX_sv(ikl,:)  =0.
+      zzsnsv(ikl,:)  =0.
+      psi_sv(ikl,:)  =0. 
+      Khydsv(ikl,:)  =0. 
+      EExcsv(ikl)  =0.    
+
+
       END DO
 

LMDZ6/trunk/libf/phylmd/inlandsis/inlandsis.F

@@ -1 @@
-      subroutine INLANDSIS(SnoMod,BloMod,jjtime)
+      subroutine INLANDSIS(SnoMod,BloMod,jjtime,debut)
 
       USE dimphy
@@ -173 +173 @@
       USE VARySV
       USE VARtSV
-      USE surface_data, only: iflag_tsurf_inlandsis
-
+      USE surface_data, ONLY: is_ok_z0h_rn,
+     .                        is_ok_density_kotlyakov,
+     .                        prescribed_z0m_snow,
+     .                        iflag_z0m_snow,
+     .                        iflag_tsurf_inlandsis,
+     .                        iflag_temp_inlandsis,
+     .                        discret_xf, buf_sph_pol,buf_siz_pol      
 
       IMPLICIT NONE
@@ -180 +185 @@
       logical   SnoMod
       logical   BloMod
+      logical   debut
       integer   jjtime
 
@@ -213 +219 @@
       integer   IceMsk,IcIndx(klonv)          !      Ice / No      Ice Mask
       integer   SnoMsk                        ! Snow     / No Snow     Mask
-
       real      roSMin,roSMax,roSn_1,roSn_2,roSn_3   ! Fallen Snow Density (PAHAUT)
       real      Dendr1,Dendr2,Dendr3          ! Fallen Snow Dendric.(GIRAUD)
       real      Spher1,Spher2,Spher3,Spher4   ! Fallen Snow Spheric.(GIRAUD)
       real      Polair                        ! Polar  Snow Switch
-      real      PorSno,Por_BS,Salt_f,PorRef   !
+      real      PorSno,Salt_f,PorRef   !
 c #sw real      PorVol,rWater                 !
 c #sw real      rusNEW,rdzNEW,etaNEW          !
@@ -244 +249 @@
       real      Z0m_Sn,Z0m_90                 ! Snow  Surface Roughness Length
       real      SnoWat                        ! Snow Layer    Switch
-c #RN real      rstar,alors                   !
-c #RN real      rstar0,rstar1,rstar2          !
+      real      rstar,alors                   !
+      real      rstar0,rstar1,rstar2          !
       real      SameOK                        ! 1. => Same Type of Grains
       real      G1same                        ! Averaged G1,  same Grains
@@ -263 +268 @@
       real      Sph_av                        ! Averaged    Grain Spher.
       real      Den_av                        ! Averaged    Grain Dendr.
-      real      DendOK                        ! 1. => Average is  Dendr.
       real      G1diff                        ! Averaged G1, diff. Grains
       real      G2diff                        ! Averaged G2, diff. Grains
@@ -277 +281 @@
       real      tt_c,vv_c                     ! Critical param.
       real      tt_tmp,vv_tmp,vv_virt         ! Temporary variables
-      logical   density_kotlyakov             ! .true. if Kotlyakov 1961
       real      e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations
       real      zm1, zm2, coefslope                    ! variables for surface temperature extrapolation
-
+! for Aeolian erosion and blowing snow
+      integer   nit   ,iit
+      real      Fac                           ! Correc. factor for drift ratio
+      real      dusuth,signus
+      real      sss__F,sss__N
+      real      sss__K,sss__G
+      real      us_127,us_227,us_327,us_427,us_527
+      real      VVa_OK, usuth0
+      real      ssstar
+      real      SblPom
+      real      rCd10n                        ! Square root of drag coefficient
+      real      DendOK                        ! Dendricity Switch
+      real      SaltOK                        ! Saltation  Switch
+      real      MeltOK                        ! Saltation  Switch (Melting Snow)
+      real      SnowOK                        ! Pack Top   Switch
+      real      SaltM1,SaltM2,SaltMo,SaltMx   ! Saltation  Parameters
+      real      ShearX, ShearS                ! Arg. Max Shear Stress
+      real      Por_BS                        ! Snow Porosity
+      real      Salt_us                       ! New thresh.friction velocity u*t
+      real      Fac_Mo,ArguSi,FacRho          ! Numerical factors for u*t
+      real      SaltSI(klonv,0:nsno)          ! Snow Drift Index              !
+      real      MIN_Mo                        ! Minimum Mobility Fresh Fallen *
+      character*3    qsalt_param              ! Switch for saltation flux param.
+      character*3    usth_param               ! Switch for u*t param
 
 
@@ -287 +313 @@
 
       data      T__Min / 200.00/              ! Minimum realistic Temperature
-      data      TaPole / 263.15/              ! Maximum Polar     Temperature
-      data      roSMin /  30.  /              ! Minimum Snow  Density
+      data      TaPole / 268.15/              ! Maximum Polar Temperature (value from C. Agosta)
+      data      roSMin / 300.  /              ! Minimum Snow  Density
       data      roSMax / 400.  /              ! Max Fresh Snow Density
       data      tt_c   / -2.0  /              ! Critical Temp. (degC)
@@ -305 +331 @@
       data      EmiWat /   0.99999999/        ! Emissivity of a Water Area
       data      EmiSno /   0.99999999/        ! Emissivity of Snow
+
       
 !     DATA      Emissivities                  ! Pielke, 1984, pp. 383,409
@@ -321 +348 @@
       data      Z0_ICE/    0.0010/            ! Sea-Ice Z0 = 0.0010 m (Andreas)
 !                                             !    (Ice Station Weddel -- ISW)
+! for aerolian erosion
+      data      SblPom/ 1.27/   ! Lower Boundary Height Parameter
+C +                             ! for Suspension
+C +                             ! Pommeroy, Gray and Landine 1993,
+C +                             ! J. Hydrology, 144(8) p.169
+      data      nit   / 5   /   ! us(is0,uth) recursivity: Nb Iterations
+cc#AE data      qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p
+      data      qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray
+cc#AE data      usth_param/"lis"/  ! u*t from Liston et al. 2007
+      data      usth_param/"gal"/  ! u*t from Gallee et al. 2001
+      data      SaltMx/-5.83e-2/
+
       vk2    =  vonKrm  *  vonKrm             ! Square of Von Karman Constant
 
@@ -352 +391 @@
 
 
-! Blowing Particles Threshold Friction velocity
-! =============================================
-
-c #AE       usthSV(ikl) =                     1.0e+2
-!          END DO
-!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
-
-
-
-
-! Contribution of Snow to the Surface Snow Pack
-! =============================================
-
-      IF (SnoMod)                                                 THEN
-
-
- 
-C +--Blowing Snow
-C +  ------------
- 
-        IF (BloMod) then
-         if (klonv.eq.1) then
+
+
+
+      IF (SnoMod)                            THEN
+
+ 
+C +--Aeolian erosion and Blowing Snow
+C +==================================
+
+
+
+        DO ikl=1,knonv
+            usthSV(ikl) =                     1.0e+2
+        END DO
+
+
+        IF (BloMod) THEN
+ 
+        if (klonv.eq.1) then
           if(isnoSV(1).ge.2                   .and.
-     .       TsisSV(1,max(1,isnoSV(1)))<273.  .and.
-     .       ro__SV(1,max(1,isnoSV(1)))<500.  .and.
-     .       eta_SV(1,max(1,isnoSV(1)))<epsi) then
+     .         TsisSV(1,max(1,isnoSV(1)))<273.  .and.
+     .         ro__SV(1,max(1,isnoSV(1)))<500.  .and.
+     .         eta_SV(1,max(1,isnoSV(1)))<epsi) then
 C +                       **********
                      call SISVAT_BSn
@@ -384 +420 @@
                      call SISVAT_BSn
 C +                       **********
-         endif
-        ENDIF
- 
-
-
+        endif
+
+
+
+
+
+
+
+! Calculate threshold erosion velocity for next time step
+! Unlike in sisvat, computation is of threshold velocity made here (instead of sisvaesbl)
+! since we do not use sisvatesbl for the coupling with LMDZ
+
+C +--Computation of threshold friction velocity for snow erosion
+C --------------------------------------------------------------- 
+
+        rCd10n =  1. / 26.5 ! Vt / u*t = 26.5
+                     ! Budd et al. 1965, Antarct. Res. Series Fig.13
+                     ! ratio developped during assumed neutral conditions
+ 
+
+C +--Snow Properties
+C +  ~~~~~~~~~~~~~~~
+
+        DO ikl = 1,knonv
+
+          isn      =  isnoSV(ikl)
+
+
+ 
+          DendOK   =  max(zero,sign(unun,epsi-G1snSV(ikl,isn)  ))  !
+          SaltOK   =  min(1   , max(istdSV(2)-istoSV(ikl,isn),0))  !
+          MeltOK   =     (unun                                     !
+     .             -max(zero,sign(unun,TfSnow-epsi                 !
+     .             -TsisSV(ikl,isn)  )))                           ! Melting Snow
+     .             *  min(unun,DendOK                              !
+     .                  +(1.-DendOK)                               !
+     .                      *sign(unun,     G2snSV(ikl,isn)-1.0))  ! 1.0 for 1mm
+          SnowOK   =  min(1   , max(isnoSV(ikl)      +1 -isn ,0))  ! Snow Switch
+ 
+          G1snSV(ikl,isn) =      SnowOK *    G1snSV(ikl,isn)
+     .                  + (1.- SnowOK)*min(G1snSV(ikl,isn),G1_dSV)
+          G2snSV(ikl,isn) =      SnowOK *    G2snSV(ikl,isn)
+     .                  + (1.- SnowOK)*min(G2snSV(ikl,isn),G1_dSV)
+ 
+          SaltOK   =  min(unun, SaltOK + MeltOK) * SnowOK
+ 
+ 
+C +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
+C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          SaltM1   = -0.750e-2 * G1snSV(ikl,isn)
+     .             -0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case
+C +     CAUTION:  Guyomarc'h & Merindol Dendricity Sign is +
+C +     ^^^^^^^^                    MAR Dendricity Sign is -
+          SaltM2   = -0.833d-2 * G1snSV(ikl,isn)
+     .             -0.583d-2 * G2snSV(ikl,isn)+ 0.833d00 !non-dendritic case
+ 
+c       SaltMo   = (DendOK   * SaltM1 + (1.-DendOK) *     SaltM2       )
+          SaltMo   = 0.625 !SaltMo pour d=s=0.5
+ 
+!weighting SaltMo with surface snow density (Vionnet et al. 2012)
+cc#AE   FacRho   = 1.25 - 0.0042 * ro__SV(ikl,isn)
+cc#AE   SaltMo   = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow
+          MIN_Mo   =  0.
+c       SaltMo   =  max(SaltMo,MIN_Mo)
+c       SaltMo   =  SaltOK   * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx)
+c #TUNE SaltMo   =  SaltOK   * SaltMo - (1.-SaltOK) *     0.9500
+          SaltMo   =  max(SaltMo,epsi-unun)
+ 
+C +--Influence of Density on Threshold Shear Stress
+C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          Por_BS =  1. - 300. / ro_Ice
+          ShearS = Por_BS / (1.-Por_BS)
+C +...         SheaBS =  Arg(sqrt(shear = max shear stress in snow)):
+C +            shear  =  3.420d00 * exp(-(Por_BS      +Por_BS)
+C +  .                                  /(unun        -Por_BS))
+C +            SheaBS :  see de Montmollin         (1978),
+C +                      These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124
+ 
+C +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
+C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          ArguSi      =     -0.085 *us__SV(ikl)/rCd10n
+!V=u*/sqrt(CD) eqs 2 to 4 Gallee et al. 2001
+ 
+          SaltSI(ikl,isn) = -2.868 * exp(ArguSi) + 1 + SaltMo
+ 
+
+C +--Threshold Friction Velocity
+C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          if(ro__SV(ikl,isn)>300.) then
+             Por_BS      =  1.000       - ro__SV(ikl,isn)     /ro_Ice
+          else
+             Por_BS      =  1.000  - 300. /ro_Ice
+          endif
+ 
+          ShearX =  Por_BS/max(epsi,1.-Por_BS)
+          Fac_Mo = exp(-ShearX+ShearS)
+C +     Gallee et al., 2001    eq 5, p5
+ 
+          if (usth_param .eq. "gal") then
+            Salt_us   =   (log(2.868) - log(1 + SaltMo)) * rCd10n/0.085
+            Salt_us   = Salt_us * Fac_Mo
+C +...  Salt_us   :  Extension of  Guyomarc'h & Merindol 1998 with
+C +...              de Montmollin (1978). Gallee et al. 2001
+          endif
+ 
+          if (usth_param .eq. "lis") then !Liston et al. 2007
+            if(ro__SV(ikl,isn)>300.) then
+              Salt_us   = 0.005*exp(0.013*ro__SV(ikl,isn))
+            else
+              Salt_us   = 0.01*exp(0.003*ro__SV(ikl,isn))
+            endif
+          endif
+ 
+          SnowOK   =  1 -min(1,iabs(isn-isnoSV(ikl))) !Switch new vs old snow
+ 
+          usthSV(ikl) =     SnowOK *   (Salt_us)
+     .                + (1.-SnowOK)*    usthSV(ikl)
+ 
+        END DO
+
+
+ 
+!  Feeback between blowing snow turbulent Scale  u* (commented here
+!  since ustar is an input variable (not in/out) of inlandsis)
+!  -----------------------------------------------------------------
+
+
+!           VVa_OK      =  max(0.000001,       VVaSBL(ikl))
+!           sss__N      =  vonkar      *       VVa_OK
+!           sss__F      = (sqrCm0(ikl) - psim_z + psim_0)
+!           usuth0      =  sss__N /sss__F                ! u* if NO Blow. Snow
+ 
+!           sss__G      =  0.27417     * gravit
+ 
+! !  ______________               _____
+! !  Newton-Raphson (! Iteration, BEGIN)
+! !  ~~~~~~~~~~~~~~               ~~~~~
+!           DO iit=1,nit
+!           sss__K      =  gravit      * r_Turb * A_Turb *za__SV(ikl)
+!      .                                     *rCDmSV(ikl)*rCDmSV(ikl)
+!      .                           /(1.+0.608*QaT_SV(ikl)-qsnoSV(ikl))
+!           us_127      =  exp(    SblPom *log(us__SV(ikl)))
+!           us_227      =  us_127         *    us__SV(ikl)
+!           us_327      =  us_227         *    us__SV(ikl)
+!           us_427      =  us_327         *    us__SV(ikl)
+!           us_527      =  us_427         *    us__SV(ikl)
+ 
+!           us__SV(ikl) =  us__SV(ikl)
+!      .    - (  us_527     *sss__F     /sss__N
+!      .      -  us_427
+!      .      -  us_227     *qsnoSV(ikl)*sss__K
+!      .      + (us__SV(ikl)*us__SV(ikl)-usthSV(ikl)*usthSV(ikl))/sss__G)
+!      .     /(  us_427*5.27*sss__F     /sss__N
+!      .      -  us_327*4.27
+!      .      -  us_127*2.27*qsnoSV(ikl)*sss__K
+!      .      +  us__SV(ikl)*2.0                                 /sss__G)
+ 
+!           us__SV(ikl)= min(us__SV(ikl),usuth0)
+!           us__SV(ikl)= max(us__SV(ikl),epsi  )
+!           rCDmSV(ikl)=     us__SV(ikl)/VVa_OK
+! ! #AE     sss__F     =     vonkar     /rCDmSV(ikl)
+!           ENDDO
+ 
+! !  ______________               ___
+! !  Newton-Raphson (! Iteration, END  )
+! !  ~~~~~~~~~~~~~~               ~~~
+ 
+!           us_127      =  exp(    SblPom *log(us__SV(ikl)))
+!           us_227      =  us_127         *    us__SV(ikl)
+ 
+! !  Momentum            Turbulent Scale  u*: 0-Limit in case of no Blow. Snow
+! !  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!           dusuth      =  us__SV(ikl) - usthSV(ikl)       ! u* - uth*
+!           signus      =  max(sign(unun,dusuth),zero)     ! 1 <=> u* - uth* > 0
+!           us__SV(ikl) =                                  !
+!      .                   us__SV(ikl)  *signus  +         ! u* (_BS)
+!      .                   usuth0                          ! u* (nBS)
+!      .                            *(1.-signus)           !        
+
+
+
+
+!  Blowing Snow        Turbulent Scale ss*
+!  ---------------------------------------
+ 
+        hSalSV(ikl) = 8.436e-2  * us__SV(ikl)**SblPom
+ 
+        if (qsalt_param .eq. "pom") then
+          qSalSV(ikl) = (us__SV(ikl)**2 - usthSV(ikl)**2) *signus
+     .               / (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25)
+        endif
+ 
+        if (qsalt_param .eq. "bin") then
+          qSalSV(ikl) = (us__SV(ikl) * us__SV(ikl)
+     .                -usthSV(ikl) * usthSV(ikl))*signus
+     .                * 0.535 / (hSalSV(ikl) * gravit)
+        endif
+ 
+        qSalSV(ikl) = qSalSV(ikl)/rht_SV(ikl) ! conversion kg/m3 to kg/kg
+ 
+        ssstar      = rCDmSV(ikl) * (qsnoSV(ikl) - qSalSV(ikl))
+     .              * r_Turb !Bintanja 2000, BLM
+!r_Turb compensates for an overestim. of the blown snow part. fall velocity
+ 
+        uss_SV(ikl) = min(zero    , us__SV(ikl) *ssstar)
+        uss_SV(ikl) = max(-0.0001 , uss_SV(ikl))   
+
+
+
+
+        ENDIF   ! BloMod
+ 
+C + ------------------------------------------------------
 C +--Buffer Layer
-C +  ------------
+C +  -----------------------------------------------------
  
           DO ikl=1,knonv
@@ -414 +658 @@
 c #NP.         104. *sqrt( max( VV10SV(ikl)-6.0,0.0)))  ! Kotlyakov (1961)
  
-            density_kotlyakov = .true.
-c #AC       density_kotlyakov = .false.  !C.Agosta snow densisty as if BS is on b
+!          C.Agosta option for snow density, same as for BS i.e.
+!          is_ok_density_kotlyakov=.false.
 c #BS       density_kotlyakov = .false.  !C.Amory BS 2018
 C + ...     Fallen Snow Density, Adapted for Antarctica
-            if (density_kotlyakov) then
+            if (is_ok_density_kotlyakov) then
                 tt_tmp = TaT_SV(ikl)-TfSnow
                 !vv_tmp = VV10SV(ikl)
@@ -452 +696 @@
 !    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  
-c #BS       Bros_N      = frsno
-c #BS       ro_new      = ro__SV(ikl,max(1,isnoSV(ikl)))
-c #BS       ro_new      = max(Bros_N,min(roBdSV,ro_new))
-c #BS       Fac         = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
-c #BS.                     -roBdSV)/(500.-roBdSV))
-c #BS       Fac         = max(0.,min(1.,Fac))
-c #BS       dsnbSV(ikl) = Fac*dsnbSV(ikl)
-c #BS       Bros_N      = Bros_N     * (1.0-dsnbSV(ikl))
-c #BS.                  + ro_new     *      dsnbSV(ikl)
-
+         if (BloMod) then
+         Bros_N      = frsno
+         ro_new      = ro__SV(ikl,max(1,isnoSV(ikl)))
+         ro_new      = max(Bros_N,min(roBdSV,ro_new))
+         Fac         = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
+     .               -roBdSV)/(500.-roBdSV))
+         Fac         = max(0.,min(1.,Fac))
+         dsnbSV(ikl) = Fac*dsnbSV(ikl)
+         Bros_N      = Bros_N     * (1.0-dsnbSV(ikl))
+     .               + ro_new     *      dsnbSV(ikl)
+         endif
 
  
@@ -480 +725 @@
      .               max(Spher1*VV__SV(ikl)+Spher2,     !     Sphericity
      .                   Spher3                   ))    !
+! EV: now control buf_sph_pol and bug_siz_pol in physiq.def
             Buf_G1      = (1. - Polair) *   Buf_G1      ! Temperate Snow
-     .                        + Polair  *   G1_dSV      ! Polar     Snow
+     .                        + Polair  *   buf_sph_pol ! Polar Snow 
             Buf_G2      = (1. - Polair) *   Buf_G2      ! Temperate Snow
-     .                        + Polair  *   ADSdSV      ! Polar     Snow
+     .                        + Polair  *   buf_siz_pol ! Polar Snow
                 G1      =                   Buf_G1      ! NO  Blown Snow
                 G2      =                   Buf_G2      ! NO  Blown Snow
 
- 
+
+
+            IF (BloMod) THEN 
+
 !     S.1. Meme  Type  de Neige  / same Grain Type
 !          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-c #BS       SameOK  =  max(zero,
-c #BS.                     sign(unun,    Buf_G1             *G1_dSV
-c #BS.                                 - eps_21                    ))
-c #BS       G1same  = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV)
-c #BS       G2same  = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV)
+
+           SameOK  =  max(zero,
+     .         sign(unun,    Buf_G1             *G1_dSV
+     .                            - eps_21                    ))
+           G1same  = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV)
+           G2same  = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV)
 !           Blowing Snow Properties:                         G1_dSV, ADSdSV
  
 !     S.2. Types differents / differents Types
 !          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-c #BS       typ__1  =  max(zero,sign(unun,epsi-Buf_G1))   ! =1.=> Dendritic
-c #BS       zroNEW  =     typ__1  *(1.0-dsnbSV(ikl))      ! fract.Dendr.Lay.
-c #BS.              + (1.-typ__1) *     dsnbSV(ikl)       !
-c #BS       G1_NEW  =     typ__1  *Buf_G1                 ! G1 of Dendr.Lay.
-c #BS.              + (1.-typ__1) *G1_dSV                 !
-c #BS       G2_NEW  =     typ__1  *Buf_G2                 ! G2 of Dendr.Lay.
-c #BS.              + (1.-typ__1) *ADSdSV                 !
-c #BS       zroOLD  = (1.-typ__1) *(1.0-dsnbSV(ikl))      ! fract.Spher.Lay.
-c #BS.              +     typ__1  *     dsnbSV(ikl)       !
-c #BS       G1_OLD  = (1.-typ__1) *Buf_G1                 ! G1 of Spher.Lay.
-c #BS.              +     typ__1  *G1_dSV                 !
-c #BS       G2_OLD  = (1.-typ__1) *Buf_G2                 ! G2 of Spher.Lay.
-c #BS.              +     typ__1  *ADSdSV                 !
-c #BS       SizNEW  =    -G1_NEW  *DDcdSV/G1_dSV          ! Size  Dendr.Lay.
-c #BS.               +(1.+G1_NEW         /G1_dSV)         !
-c #BS.                  *(G2_NEW  *DScdSV/G1_dSV          !
-c #BS.               +(1.-G2_NEW         /G1_dSV)*DFcdSV) !
-c #BS       SphNEW  =     G2_NEW         /G1_dSV          ! Spher.Dendr.Lay.
-c #BS       SizOLD  =     G2_OLD                          ! Size  Spher.Lay.
-c #BS       SphOLD  =     G1_OLD         /G1_dSV          ! Spher.Spher.Lay.
-c #BS       Siz_av =     (zroNEW*SizNEW+zroOLD*SizOLD)    ! Averaged Size
-c #BS       Sph_av = min( zroNEW*SphNEW+zroOLD*SphOLD     !
-c #BS.                   ,  unun)                         ! Averaged Sphericity
-c #BS       Den_av = min((Siz_av -(    Sph_av *DScdSV     !
-c #BS.                            +(1.-Sph_av)*DFcdSV))   !
-c #BS.                 / (DDcdSV -(    Sph_av *DScdSV     !
-c #BS.                            +(1.-Sph_av)*DFcdSV))   !
-c #BS.                   ,  unun)                         !
-c #BS       DendOK  = max(zero,                           !
-c #BS.                    sign(unun,     Sph_av *DScdSV   ! Small   Grains
-c #BS.                              +(1.-Sph_av)*DFcdSV   ! Faceted Grains
-c #BS.                              -    Siz_av        )) !
+           typ__1  =  max(zero,sign(unun,epsi-Buf_G1))   ! =1.=> Dendritic
+           zroNEW  =     typ__1  *(1.0-dsnbSV(ikl))      ! fract.Dendr.Lay.
+     .            + (1.-typ__1) *     dsnbSV(ikl)       !
+           G1_NEW  =     typ__1  *Buf_G1                 ! G1 of Dendr.Lay.
+     .            + (1.-typ__1) *G1_dSV                 !
+           G2_NEW  =     typ__1  *Buf_G2                 ! G2 of Dendr.Lay.
+     .            + (1.-typ__1) *ADSdSV                 !
+           zroOLD  = (1.-typ__1) *(1.0-dsnbSV(ikl))      ! fract.Spher.Lay.
+     .            +     typ__1  *     dsnbSV(ikl)       !
+           G1_OLD  = (1.-typ__1) *Buf_G1                 ! G1 of Spher.Lay.
+     .            +     typ__1  *G1_dSV                 !
+           G2_OLD  = (1.-typ__1) *Buf_G2                 ! G2 of Spher.Lay.
+     .            +     typ__1  *ADSdSV                 !
+           SizNEW  =    -G1_NEW  *DDcdSV/G1_dSV          ! Size  Dendr.Lay.
+     .            +(1.+G1_NEW         /G1_dSV)          !
+     .                  *(G2_NEW  *DScdSV/G1_dSV        !
+     .            +(1.-G2_NEW         /G1_dSV)*DFcdSV)  !
+           SphNEW  =     G2_NEW         /G1_dSV          ! Spher.Dendr.Lay.
+           SizOLD  =     G2_OLD                          ! Size  Spher.Lay.
+           SphOLD  =     G1_OLD         /G1_dSV          ! Spher.Spher.Lay.
+           Siz_av  =     (zroNEW*SizNEW+zroOLD*SizOLD)   ! Averaged Size
+           Sph_av  = min( zroNEW*SphNEW+zroOLD*SphOLD    !
+     .                 ,  unun)                         ! Averaged Sphericity
+           Den_av  = min((Siz_av -(    Sph_av *DScdSV    !
+     .            +(1.-Sph_av)*DFcdSV))                 !
+     .            / (DDcdSV -(    Sph_av *DScdSV        !
+     .            +(1.-Sph_av)*DFcdSV))                 !
+     .                   ,  unun)                       !
+           DendOK  = max(zero,                           !
+     .                    sign(unun,     Sph_av *DScdSV   ! Small   Grains
+     .                              +(1.-Sph_av)*DFcdSV   ! Faceted Grains
+     .                              -    Siz_av        )) !
 C +...      REMARQUE: le  type moyen (dendritique ou non) depend
 C +         ^^^^^^^^  de la  comparaison avec le diametre optique
@@ -538 +788 @@
 C +                   of a recent snow    having zero dendricity
  
-c #BS       G1diff  =(   -DendOK *Den_av
-c #BS.               +(1.-DendOK)*Sph_av) *G1_dSV
-c #BS       G2diff  =     DendOK *Sph_av  *G1_dSV
-c #BS.               +(1.-DendOK)*Siz_av
-c #BS       G1      =     SameOK *G1same
-c #BS.               +(1.-SameOK)*G1diff
-c #BS       G2      =     SameOK *G2same
-c #BS.               +(1.-SameOK)*G2diff
- 
+           G1diff  =(   -DendOK *Den_av
+     .            +(1.-DendOK)*Sph_av) *G1_dSV
+           G2diff  =     DendOK *Sph_av  *G1_dSV
+     .            +(1.-DendOK)*Siz_av
+           G1      =     SameOK *G1same
+     .            +(1.-SameOK)*G1diff
+           G2      =     SameOK *G2same
+     .            +(1.-SameOK)*G2diff
+           ENDIF
+
 
  
@@ -634 +885 @@
      .                            /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m]
  
-
  
           END DO
@@ -640 +890 @@
 
 
-! Snow Pack Discretization
-! ========================
-
-!            **********
+! Snow Pack Discretization(option XF in MAR)
+! ==========================================
+
+          
+      if (discret_xf.AND.klonv.eq.1) then
+
+       if(isnoSV(1).ge.1.or.NLaysv(1).ge.1) then
+C +          **********
+         call SISVAT_zSn
+C +          **********
+       endif
+      else
+C +          **********
         call SISVAT_zSn
-!            **********
-
-!            **********
+C +          **********
+      endif
+ 
+C +          **********
 ! #ve   call SISVAT_wEq('_zSn  ',0)
-!            **********
-
-
+C +          **********
 
 ! Add a new Snow Layer
@@ -664 +922 @@
             TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl))
      .                  + min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl) 
-
             ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl))
      .                      + Brossv(ikl)     *    NLaysv(ikl) 
@@ -699 +956 @@
 
 
-      END IF
+      END IF  ! SnoMod
 
 
@@ -740 +997 @@
 ! =============================
 !Etienne: as in inlandis we do not call vgopt, we need to define
-!the albedo  alb_SV and to calculate the
+!the albedo alb_SV and to calculate the
 !absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv
 
@@ -810 +1067 @@
 
 
-! Aerodynamic Resistance
-! ^^^^^^^^^^^^^^^^^^^^^^
-
-
-       DO ikl=1,knonv
-          ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
-          rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
-        END DO
+! Aerodynamic Resistance (calculated from drags given by LMDZ)
+! Commented because already calculated in surf_inlandsis_mod
+! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+!       DO ikl=1,knonv
+!          ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
+!          rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
+!        END DO
 
 
@@ -825 +1081 @@
 
 
-      if (iflag_tsurf_inlandsis .eq. 0) then
+      if (iflag_temp_inlandsis .eq. 0) then
 
        call SISVAT_TSo
 
       else
+        DO ikl=1,knonv
+        Tsf_SV(ikl)=Tsrfsv(ikl)
+        END DO
 
        call SISVAT_TS2
@@ -938 +1197 @@
 ! Surface Temperature
 ! ^^^^^^^^^^^^^^^^^^^^
-!           Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
-
+
+          IF (iflag_tsurf_inlandsis .EQ. 0) THEN   
+
+            Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
+
+          ELSE IF (iflag_tsurf_inlandsis .GT. 0) THEN
 ! Etienne: extrapolation from the two uppermost levels:
 
@@ -959 +1222 @@
 
 
-        END DO
-
+         ELSE !(default)
+
+           Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
+
+         END IF
+
+
+         END DO
 
 ! Snow Pack Properties (sphericity, dendricity, size)
@@ -967 +1236 @@
       IF (SnoMod)                                                 THEN
 
-!            **********
+      if (discret_xf .AND. klonv.eq.1) then
+      if(isnoSV(1).ge.1) then
+C +          **********
+      call SISVAT_GSn
+C +          **********
+      endif
+      else
+C +          **********
         call SISVAT_GSn
-!            **********
-
-!            **********
-! #ve   call SISVAT_wEq('_GSn  ',0)
-!            **********
-
+C +          **********
+      endif
 
 
@@ -990 +1262 @@
 C +--Roughness Length for Momentum
 C +  -----------------------------
+
+! ETIENNE WARNING: changes have been made wrt original SISVAT
  
 C +--Land+Sea-Ice / Ice-free Sea Mask
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-        DO ikl=1,klonv
+        DO ikl=1,knonv
           IcIndx(ikl) = 0
         ENDDO
         DO isn=1,nsno
-        DO ikl=1,klonv
+        DO ikl=1,knonv
+
           IcIndx(ikl) = max(IcIndx(ikl),
-     .                      isn*max(0,
-     .                              sign(1,
-     .                                   int(ro__SV(ikl,isn)-900.))))
+     .                  isn*max(0,
+     .                  sign(1,
+     .                  int(ro__SV(ikl,isn)-900.))))
         ENDDO
         ENDDO
  
-        DO ikl=1,klonv
+        DO ikl=1,knonv
           LISmsk    =     1. ! in inlandsis, land only
           IceMsk    =     max(0,sign(1   ,IcIndx(ikl)-1)  )
           SnoMsk    = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
 
- 
-
-          Z0mLnd      =max( Z0_ICE    ,    5.e-5  )  ! Min set := Z0 on *
 
 C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
           Z0m_nu =       5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
- 
+
 C +--Z0 Saltat.Regime over Snow (Gallee  et al., 2001, BLM 99 (19) p.11)
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
+
           u2star =       us__SV(ikl) *us__SV(ikl)
           Z0mBSn =       u2star      *0.536e-3   -  61.8e-6
           Z0mBSn =   max(Z0mBS0      ,Z0mBSn)
- 
+
 C +--Z0 Smooth + Saltat. Regime
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
           Z0enSV(ikl) =  Z0m_nu
      .                +  Z0mBSn
- 
-C +--Rough   Snow Surface Roughness Length (Typical Value)
-C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-c #tz     Z0m_Sn =    0.250e-3 ! Andreas 1995, CRREL Report 95-16, fig.1&p.2
-                               ! z0r~(10-d)*exp(-vonkar/sqrt(1.5e-03))-5.e-5
-          Z0m_Sn =    2.000e-3 ! Calibration    of MAR
-c #TZ     Z0m_Sn =    1.000e-3 ! Exemple Tuning in RACMO
-c #TZ     Z0m_Sn =    0.500e-3 ! Exemple Tuning in MAR
- 
+
+        
+! Calculation of snow roughness length
+!=====================================
+          IF (iflag_z0m_snow .EQ. 0) THEN
+
+          Z0m_Sn=prescribed_z0m_snow
+
+          ELSE IF (iflag_z0m_snow .EQ. 1) THEN
+
+          Z0m_Sn=Z0enSV(ikl) 
+
+          ELSE IF (iflag_z0m_snow .EQ. 2) THEN                             
+
 C +--Rough   Snow Surface Roughness Length (Variable Sastrugi Height)
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -1045 +1323 @@
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013
+
+          
           coefa = 0.1658 !0.1862 !Ant
           coefb = -50.3869 !-55.7718 !Ant
@@ -1055 +1335 @@
           coefc = log(z03/z02)/(ta3-ta2)
           coefd = log(z03)-coefc*ta3
+
           if (TaT_SV(ikl) .lt. ta1) then
             Z0_obs = z01
@@ -1066 +1347 @@
           endif
  
-
-! pour le moment, on choisit une valeur fixe
-          Z0_obs = 1.000e-3
- 
-cCA       Snow roughness lenght deduced from observations
-cCA       (parametrization if no Blowing Snow)
-cCA       ----------------------------------- C. Agosta 09-2016 -----
-cCA       Substract Z0enSV(ikl) because re-added later in Z0mnSV(ikl)
-          Z0m_Sn = Z0_obs - Z0enSV(ikl)
-cCA       -----------------------------------------------------------
- 
-          param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING
- 
+          Z0m_Sn=Z0_obs
+
+
+          ELSE
+
+          Z0m_Sn=0.500e-3  ! default=0.500e-3m (tuning of MAR) 
+
+          ENDIF
+ 
+
+
+!          param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING 
 c #SZ     Z0Sa_N =                   (us__SV(ikl) -0.2)*param   ! 0.0001=TUNING
 c #SZ.           * max(zero,sign(unun,TfSnow-eps9
@@ -1109 +1389 @@
 c #ZN     Z0enSV(ikl) =  max(Z0enSV(ikl), Z0m_nu)
  
+
 C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^  ams.confex.com/ams/pdfpapers/68601.pdf)
@@ -1132 +1413 @@
 c #ZA     Z0m_Sn =           DDs_SV(ikl)* Z0m_90 / 45.
 c #ZA.         - DDs_SV(ikl)*DDs_SV(ikl)* Z0m_90 /(90.*90.)
- 
-C +--Z0  (Erosion)    over Snow (instantaneous or time average)
+
+
+
+
+C +--Z0  (Erosion)    over Snow (instantaneous)
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
           Z0e_SV(ikl) =  Z0enSV(ikl)
-          Z0e_SV(ikl) =  Z0emSV(ikl)
- 
-C +--Momentum  Roughness Length
-C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^                              ! Contribution of
-          Z0mnSV(ikl) =  Z0mLnd                              ! land Form
-     .                + (Z0m_Sn                              ! Sastrugi   Form
-     .                +  Z0enSV(ikl))   *SnoMsk              ! Snow    Erosion
+ 
+C +--Momentum  Roughness Length (Etienne: changes wrt original SISVAT)
+C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                              
+          Z0mnSV(ikl) =  Z0m_nu *(1-SnoMsk)                     ! Ice z0
+     .                + (Z0m_Sn)*SnoMsk                         ! Snow Sastrugi Form and Snow Erosion
  
 
@@ -1154 +1436 @@
 c #GL.                     /(920.00                 -600.))) !
  
-C +--Mom. Roughness Length, Instantaneous OR Box Moving Average in Time
-C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+C +--Mom. Roughness Length, Instantaneous
+C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
           Z0m_SV(ikl) =  Z0mnSV(ikl)                         ! Z0mnSV  instant.
-!          Z0m_SV(ikl) =  Z0mmSV(ikl)                         ! Z0mnSV  Average
- 
-C +--Corrected Threshold Friction Velocity before Erosion    ! Marticorena and
-C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^    ! Bergametti 1995
-! not used anymore since Marticorena and Bergametti disabled !CK 18/07/2018
-cc #BS     Z0e_SV(ikl) =   min(Z0m_SV(ikl),Z0e_SV(ikl))       !
-cc #MB     f_eff=    log(0.35*(0.1        /Z0e_SV(ikl))**0.8) ! JGR 100
-cc #MB     f_eff=1.-(log(      Z0m_SV(ikl)/Z0e_SV(ikl)      ))! (20) p. 16420
-cc #MB.            /(max(      f_eff      ,epsi             ))! p.16426 2nd ?
-cc #MB     f_eff=    max(      f_eff      ,epsi              )! CONTROL
-cc #MB     f_eff=1.0   -(1.0 - f_eff)     /5.00               ! TUNING
-cc #MB     f_eff=    min(      f_eff      ,1.00              )!
-cc #MB    usthSV(ikl) =       usthSV(ikl)/f_eff              !
- 
  
  
@@ -1177 +1445 @@
  
           Z0hnSV(ikl) =     Z0mnSV(ikl)/  7.4
-c #SH     Z0hnSV(ikl) =     Z0mnSV(ikl)/100.0
-C +                         Z0h = Z0m  /100.0   over the Sahel
-C +                                            (Taylor & Clark, QJRMS 127,p864)
- 
-c #RN     rstar       =     Z0mnSV(ikl) * us__SV(ikl) / akmol
-c #RN     rstar       = max(epsi,min(rstar,thous))
-c #RN     alors       =          log(rstar)
-c #RN     rstar0      = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar))
-c #RN.                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
-c #RN.                *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar))
-c #RN.                + 0.317e0
-c #RN.                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
-c #RN     rstar1      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
-c #RN.                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
-c #RN.                *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar))
-c #RN.                - 0.565
-c #RN.                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
-c #RN     rstar2      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
-c #RN.                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
-c #RN.                *(0.      * max(zero,sign(unun,2.500e0 - rstar))
-c #RN.                - 0.183
-c #RN.                *(unun    - max(zero,sign(unun,2.500e0 - rstar))))
- 
-cXF    #RN does not work over bare ice
-cXF    MAR is then too warm and not enough melt
- 
-c #RN     if(ro__SV(ikl,isnoSV(ikl))>50
-c #RN.  .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
- 
-c #RN     Z0hnSV(ikl) = max(zero
-c #RN.                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))
-c #RN.                * exp(rstar0+rstar1*alors+rstar2*alors*alors)
-c #RN.                * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero
-c #RN.                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
- 
-c #RN     endif
+ 
+          IF (is_ok_z0h_rn) THEN
+
+          rstar       =     Z0mnSV(ikl) * us__SV(ikl) / akmol
+          rstar       = max(epsi,min(rstar,R_1000))
+          alors       =          log(rstar)
+          rstar0      = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar))
+     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
+     .                *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar))
+     .                + 0.317e0
+     .                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
+          rstar1      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
+     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
+     .                *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar))
+     .                - 0.565
+     .                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
+          rstar2      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
+     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
+     .                *(0.      * max(zero,sign(unun,2.500e0 - rstar))
+     .                - 0.183
+     .                *(unun    - max(zero,sign(unun,2.500e0 - rstar))))
+ 
+          
+
+!XF    #RN (is_ok_z0h_rn) does not work well over bare ice
+!XF    MAR is then too warm and not enough melt
+ 
+         if(ro__SV(ikl,isnoSV(ikl))>50
+     .  .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
+ 
+             Z0hnSV(ikl) = max(zero
+     .                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))
+     .                * exp(rstar0+rstar1*alors+rstar2*alors*alors)
+     .                * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero
+     .                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
+ 
+          endif
+
+
+          ENDIF
   
           Z0h_SV(ikl) =     Z0hnSV(ikl)
-!          Z0h_SV(ikl) =     Z0hmSV(ikl)
  
 

LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_bsn.F

@@ -9 +9 @@
 C |                                                                        |
 C |   SISVAT_bsn computes the snow erosion mass according to both the      |
-C |   theoretical maximum erosion amount computed in SISVATesbl and the    |
+C |   theoretical maximum erosion amount computed in inlandsis and the     |
 C |   availability of snow (currently in the uppermost snow layer only)    |
 C |                                                                        |
-C |   Preprocessing  Option: SISVAT IO (not always a standard preprocess.) |
-C |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^                                     |
-C |   FILE                 |      CONTENT                                  |
-C |   ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
-C | # stdout               | #sb: OUTPUT of Snow Erosion                   |
-C |                        |      unit  6, SubRoutine  SISVAT_BSn **ONLY** |
 C +------------------------------------------------------------------------+
  

LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_qsn.F

@@ -61 +61 @@
       use VARxSV
       use VARySV
+      use surface_data, only: is_ok_slush,opt_runoff_ac
+
 
       IMPLICIT NONE
@@ -235 +237 @@
  
       DO ikl=1,knonv
-       DO isn=min(nsno,isnoSV(ikl)+1),1,-1
+
+      DO isn=min(nsno,isnoSV(ikl)+1),1,-1
 ! EV          DO isn=nsno,1,-1
 C +--Energy, store Previous Content
@@ -243 +246 @@
      .                + ro__SV(ikl,isn) * Cn_dSV * dTSnow
      .                                           * dzsnSV(ikl,isn)
-
-          Tsave       = TsisSV(ikl,isn)
-
           TsisSV(ikl,isn) =                        TfSnow
  
@@ -312 +312 @@
           rdzNEW      = WaFrez + rdzsno
           ro__SV(ikl,isn) =      rdzNEW /max(epsi, dzsnSV(ikl,isn))
-
-! EV: condition on Enfrez
-!          if (EnFrez .eq. 0.) then
-          
-          TsisSV(ikl,isn) = Tsave
-!          else
           TsisSV(ikl,isn) =      TfSnow
      .                + EnFrez /(Cn_dSV *max(epsi, rdzNEW)        )
-!          end if
           EExcsv(ikl) =          EExcsv(ikl)     - EnFrez
           wer_SV(ikl) = WaFrez
@@ -499 +492 @@
           rusnew      = rusnSV(ikl) * SWf_SV(ikl)
  
-          if(isnoSV(ikl)<=1) rusnew = 0.
+          if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0.
           !if(ivgtSV(ikl)>=1) rusnew = 0.
  
 c #EU                        rusnew = 0.
-c #AC                        rusnew = 0.
+c #AC               rusnew = 0.
+
           RnofSV(ikl) = RnofSV(ikl)
      .                +(rusnSV(ikl) - rusnew     ) / dt__SV
@@ -545 +539 @@
         ENDDO
  
-C +--Slush Formation (CAUTION: ADD RunOff Possibility before Activation)
+C +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
 C +  ---------------  ^^^^^^^  ^^^
  
- 
-c #SU DO  ikl=1,knonv
-c #SU  DO isn=1,isnoSV(ikl)
-c #SU     kSlush = min(1,max(0,isn+1-ispiSV(ikl)))        ! Slush Switch
+      IF (is_ok_slush) THEN 
+
+      DO  ikl=1,knonv
+       DO isn=1,isnoSV(ikl)
+          kSlush = min(1,max(0,isn+1-ispiSV(ikl)))        ! Slush Switch
  
 C +--Available Additional Pore   Volume [-]
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-c #SU     PorVol = 1. - ro__SV(ikl,isn)                    ! [--]
-c #SU.           *(1. - eta_SV(ikl,isn))/ ro_Ice           !
-c #SU.           -      eta_SV(ikl,isn)                    !
-c #SU.                 *ro__SV(ikl,isn) / ro_Wat           !
-c #SU     PorVol =  max(PorVol          , zero  )          !
-c #SU     zWater =      dzsnSV(ikl,isn) * PorVol * 1000.   ! [mm] OR [kg/m2]
-c #SU.           * (1. -SWS_SV(ikl)                        ! 0 <=> freezing
-c #SU.                *(1 -min(1,iabs(isn-isnoSV(ikl)))))  ! 1 <=> isn=isnoSV
-c #SU     zSlush =  min(rusnSV(ikl)     , zWater)          ! [mm] OR [kg/m2]
-c #SU     ro_new      =(dzsnSV(ikl,isn) * ro__SV(ikl,isn)  !
-c #SU.                 +zSlush                           ) !
-c #SU.            / max(dzsnSV(ikl,isn) , epsi           ) !
-c #SU     if(ro_new<ro_Ice+20) then ! MAX 940kg/m3         !
-c #SU      rusnSV(ikl)  = rusnSV(ikl)          - zSlush    ! [mm] OR [kg/m2]
-c #SU      RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)
-c #SU      eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn)      !
-c #SU.                     *(1.     - eta_SV(ikl,isn)))    !
-c #SU.                / max (ro_new , epsi            )    !
-c #SU      ro__SV(ikl,isn) =      ro_new                   !
-c #SU     endif
-c #SU   END DO
-c #SU END DO
- 
+          PorVol = 1. - ro__SV(ikl,isn)                    ! [--]
+     .           *(1. - eta_SV(ikl,isn))/ ro_Ice           !
+     .           -      eta_SV(ikl,isn)                    !
+     .                 *ro__SV(ikl,isn) / ro_Wat           !
+          PorVol =  max(PorVol          , zero  )          !
+          zWater =      dzsnSV(ikl,isn) * PorVol * 1000.   ! [mm] OR [kg/m2]
+     .           * (1. -SWS_SV(ikl)                        ! 0 <=> freezing
+     .                *(1 -min(1,iabs(isn-isnoSV(ikl)))))  ! 1 <=> isn=isnoSV
+          zSlush =  min(rusnSV(ikl)     , zWater)          ! [mm] OR [kg/m2]
+          ro_new      =(dzsnSV(ikl,isn) * ro__SV(ikl,isn)  !
+     .                 +zSlush                           ) !
+     .            / max(dzsnSV(ikl,isn) , epsi           ) !
+          if(ro_new<ro_Ice+20) then ! MAX 940kg/m3         !
+           rusnSV(ikl)  = rusnSV(ikl)          - zSlush    ! [mm] OR [kg/m2]
+           RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)
+           eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn)      !
+     .                     *(1.     - eta_SV(ikl,isn)))    !
+     .                / max (ro_new , epsi            )    !
+           ro__SV(ikl,isn) =      ro_new                   !
+          endif
+        END DO
+      END DO
+      END IF
  
 C +--Impact of the Sublimation/Deposition on the Surface Mass Balance

LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F

@@ -1 @@
- 
- 
       subroutine SnOptP(jjtime)
  
@@ -30 +28 @@
 C |   OUTPUT:  albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
 C |   ^^^^^^   sEX_sv   : Verticaly Integrated Extinction Coefficient      |
+C |            DOPsnSV   : Snow Optical diameter [m]                       |
 C |                                                                        |
 C |   Internal Variables:                                                  |
@@ -68 +67 @@
       use VARySV
       use VARtSV
-      USE surface_data, only: iflag_albzenith
+      USE surface_data, only: iflag_albcalc,correc_alb
 
       IMPLICIT NONE
@@ -89 +88 @@
 c #AG real      agesno
  
-      integer   isn   ,ikl   ,isn1  
+      integer   isn   ,ikl   ,isn1, i  
       real      sbeta1,sbeta2,sbeta3,sbeta4,sbeta5
       real      AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK
       real      dalbed,dalbeS,dalbeW
-      real      bsegal,czemax,csegal
+      real      bsegal,czemax,csegal,csza
       real      RoFrez,DiffRo,SignRo,SnowOK,OpSqrt
       real      albSn1,albIc1,a_SnI1,a_SII1
@@ -103 +102 @@
       real      albedo_old,albCor
       real      ro_ave,dz_ave,minalb
-      real      thetazs,thetazs0,aap,bbp,ccp,alb0      ! parameters for zenoth angle effect at Dome C
-
-
+      real      l1min,l1max,l2min,l2max,l3min,l3max
+      real      l6min(6), l6max(6), albSn6(6), a_SII6(6)
+      real      lmintmp,lmaxtmp,albtmp
  
 C +--Local   DATA
@@ -138 +137 @@
       data      albmax  /0.99       /  ! Albedo max
 
-C +-- For the computation of solar zenoth angle effect from Dome C data
-C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      data      aap/0.00016/,bbp/-0.017/,ccp/1.2/
-
+C +-- wavelength limits [m] for each broad band
+
+      data      l1min/400.0e-9/,l1max/800.0e-9/
+      data      l2min/800.0e-9/,l2max/1500.0e-9/
+      data      l3min/1500.0e-9/,l3max/2800.0e-9/
+
+      data      l6min/185.0e-9,250.0e-9,400.0e-9,
+     .               690.0e-9,1190.0e-9,2380.0e-9/
+      data      l6max/250.0e-9,400.0e-9,
+     .          690.0e-9,1190.0e-9,2380.0e-9,4000.0e-9/
  
  
@@ -177 +182 @@
      .                                       *DFcdSV                 )))
           SnOpSV(ikl,isn) =  max(zero,SnOpSV(ikl,isn))
+
+C + --For outputs (Etienne)
+C + ------------------------
+          DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
         END DO
       END DO
  
+
+
  
 C +--Snow/Ice Albedo
@@ -191 +202 @@
         DO ikl=1,knonv
  
-           isn   =  max(iun,isnoSV(ikl))
+          isn   =  max(iun,isnoSV(ikl))
  
           SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
@@ -200 +211 @@
 cCA    +--Correction of snow albedo for Antarctica/Greenland
 cCA       --------------------------------------------------
-          albCor = 1.
+
+          
+          albCor = correc_alb 
 c #GL     albCor = 1.01
-c #AC     albCor = 1.01
- 
+c #AC    albCor = 1.01
+
+
+          IF (iflag_albcalc .GE. 1) THEN  ! Albedo calculation according to Kokhanovsky and Zege 2004
+
+          dalbed = 0.0
+          doptic=SnOpSV(ikl,isn)
+          csza=coszSV(ikl)
+
+          CALL albedo_kokhanovsky(l1min,l1max,csza,doptic,albSn1)
+          CALL albedo_kokhanovsky(l2min,l2max,csza,doptic,albSn2)
+          CALL albedo_kokhanovsky(l3min,l3max,csza,doptic,albSn3)
+
+          DO i=1,6
+             lmintmp=l6min(i)
+             lmaxtmp=l6max(i)
+             CALL albedo_kokhanovsky(lmintmp,lmaxtmp,csza,doptic,albtmp)
+             albSn6(i)=albtmp
+          ENDDO
+
+
+          ELSE ! Default calculation in SISVAT
+
+!    Zenith Angle Correction (Segal et al.,         1991, JAS 48, p.1025)
+!--------------------------- (Wiscombe & Warren, dec1980, JAS   , p.2723)
+!                            (Warren,               1982,  RG   , p.  81)
+!                            -------------------------------------------------
+          
+          dalbed = 0.0
+
+          csegal = max(czemax                   ,coszSV(ikl))
+c #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
+c #cz.                -       unun                          )*0.32
+c #cz.              /  bsegal
+c #cz     dalbeS = max(dalbeS,zero)
+c #cz     dalbed =     dalbeS      *       min(1,isnoSV(ikl))
+ 
+          dalbeW =(0.64 - csegal  )*0.0625  ! Warren 1982, RevGeo, fig.12b
+                                            ! 0.0625 = 5% * 1/0.8,   p.81
+                                            ! 0.64   = cos(50)
+          dalbed =     dalbeW      *       min(1,isnoSV(ikl))
+!-------------------------------------------------------------------------
+
           albSn1 =  0.96-1.580*OpSqrt
           albSn1 =  max(albSn1,AlbMin)
@@ -219 +273 @@
           albSn3 =  min(albSn3*albCor,unun)
  
+          albSn6(1:3)=albSn1
+          albSn6(4:6)=albSn2
+
           !snow albedo corection if wetsnow
 c #GL     albSn1 =  albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
 c #GL     albSn2 =  albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
 c #GL     albSn3 =  albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
+
+          ENDIF
+
  
           albSno =  So1dSV*albSn1
@@ -237 +297 @@
           albSn2 =  SnowOK*albSn2+(1.0-SnowOK)*max(albSno,minalb)
           albSn3 =  SnowOK*albSn3+(1.0-SnowOK)*max(albSno,minalb)
+          albSn6(:) =  SnowOK*albSn6(:)+(1.0-SnowOK)*max(albSno,minalb)
+
  
 c +           ro < roSdSV |          min al > aI3dSV
@@ -304 +366 @@
           a_SII3      =     albWIc      +(albSn3-albWIc)     *SnownH
           a_SII3      = min(a_SII3       ,albSn3)
- 
+
+          DO i=1,6
+          a_SII6(i)   = albWIc      +(albSn6(i)-albWIc)     *SnownH
+          a_SII6(i)   = min(a_SII6(i)       ,albSn6(i))
+          ENDDO
+
 cc #AG     agesno =      min(agsnSV(ikl,isn)          ,AgeMax)
 cc #AG     a_SII1      =     a_SII1      -0.175*agesno/AgeMax
 C +...                                   Impurities: Col de Porte Parameter.
  
- 
-!    Zenith Angle Correction (Segal et al.,         1991, JAS 48, p.1025)
-!    ----------------------- (Wiscombe & Warren, dec1980, JAS   , p.2723)
-!                            (Warren,               1982,  RG   , p.  81)
-!                            (Vignon, PhD Thesis, pp 150, conditions at Dome C)
-!                            -------------------------------------------------
- 
-         
-          dalbed = 0.0
-
-          if (iflag_albzenith .GE. 0) then
-          csegal = max(czemax                   ,coszSV(ikl))
-c #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
-c #cz.                -       unun                          )*0.32
-c #cz.              /  bsegal
-c #cz     dalbeS = max(dalbeS,zero)
-c #cz     dalbed =     dalbeS      *       min(1,isnoSV(ikl))
- 
-          dalbeW =(0.64 - csegal  )*0.0625  ! Warren 1982, RevGeo, fig.12b
-                                            ! 0.0625 = 5% * 1/0.8,   p.81
-                                            ! 0.64   = cos(50)
-          dalbed =     dalbeW      *       min(1,isnoSV(ikl))
-
-
-! Vignon PhD thesis, Dome C conditions
-
-            if (iflag_albzenith .EQ. 1) then
-            thetazs0=-bbp/(2.0*aap)
-            alb0=(bbp**2)/4./aap-(bbp**2)/2./aap+ccp
-            thetazs=max(acos(coszSV(ikl))/pi*180.0,thetazs0) ! in degrees
-
-
-            dalbeW = aap*(thetazs**2)+bbp*thetazs+ccp-alb0       
-            dalbed =     dalbeW      *       min(1,isnoSV(ikl))
-            end if
-
-
-            end if
+
  
 C +--Elsewhere    Integrated Snow/Ice Albedo
@@ -368 +398 @@
             alb3sv(ikl) =     albssv(ikl) +(a_SII3-albssv(ikl))*SIcenH
             alb3sv(ikl) = min(alb3sv(ikl)  ,a_SII3)
- 
+
             albisv(ikl) =     albssv(ikl) +(albSII-albssv(ikl))*SIcenH
             albisv(ikl) = min(albisv(ikl)  ,albSII)
- 
+
+            DO i=1,6
+            alb6sv(ikl,i) = albssv(ikl) +(a_SII6(i)-albssv(ikl))*SIcenH
+            alb6sv(ikl,i) = min(alb6sv(ikl,i)  ,a_SII6(i))
+            ENDDO
+
  
 C +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
@@ -382 +417 @@
             alb3sv(ikl) = alb3sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH
      .                  + dalbed      *    (1.-cld_SV(ikl))
+            alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 *(cld_SV(ikl)-0.5)*SIcenH
+     .                  + dalbed      *    (1.-cld_SV(ikl))
             albisv(ikl) = albisv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH
      .                  + dalbed      *    (1.-cld_SV(ikl))
@@ -397 +434 @@
             alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0             ! 33 %
      .                  * (albedo_old-albisv(ikl)) / So3dSV
- 
- 
-C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
+            alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0     ! 16 %
+     .                  * (albedo_old-albisv(ikl)) / (So1dSV/3)
+            alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0     ! 16 %
+     .                  * (albedo_old-albisv(ikl)) / (So2dSV/3)
+
+
+C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95% 
 C +  -----------------------------------------------------
  
@@ -410 +451 @@
             alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0             ! 33 %
      .                  * (albedo_old-albisv(ikl)) / So3dSV
+            alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0     ! 16 %
+     .                  * (albedo_old-albisv(ikl)) / (So1dSV/3)
+            alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0     ! 16 %
+     .                  * (albedo_old-albisv(ikl)) / (So2dSV/3)
 
 
@@ -436 +481 @@
             alb2sv(ikl) = min(max(zero,alb2sv(ikl)),albmax)
             alb3sv(ikl) = min(max(zero,alb3sv(ikl)),albmax)
- 
+            
+            DO i=1,6
+                alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax)
+            ENDDO 
         END DO
  
@@ -519 +567 @@
  
       return
+
+
       end
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+      SUBROUTINE albedo_kokhanovsky(lambdamin,lambdamax,
+     .                              cossza,dopt,albint)
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Authors: Hajar El Habchi El Fenniri, Etienne Vignon, Cecile Agosta
+!          Ghislain Picard
+! Routine that calculates the  integrated snow spectral albedo between 
+! lambdamin and lambdamax following Kokhanisky and Zege 2004,
+! Scattering optics of snow, Applied Optics, Vol 43, No7 
+! Code inspired from the snowoptics package of Ghislain Picard:
+! https://github.com/ghislainp/snowoptics
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+       
+      USE VARphy
+
+      IMPLICIT NONE
+
+! Inputs
+!--------
+      REAL lambdamin   ! minimum wavelength for integration [m]
+      REAL lambdamax   ! maximum wavelength for integration [m] 
+      REAL cossza     ! solar zenith angle cosinus
+      REAL dopt        ! optical diameter [m]
+
+!Outputs
+!-------
+      REAL albint
+
+! Local Variables
+!-----------------
+
+      REAL ropt,cosalb,norm,Pas
+      REAL SSA,alpha,gamm,R,cos30,alb30
+      INTEGER i
+
+
+      REAL B_amp                       ! amplification factor
+      PARAMETER (B_amp=1.6)  
+
+      REAL g_asy                       ! asymetry factor 
+      PARAMETER (g_asy=0.845)
+
+      INTEGER nlambda                  ! length of wavelength vector 
+      PARAMETER(nlambda=200)
+
+      REAL lmin
+      PARAMETER(lmin=185.0E-9)
+
+      REAL lmax
+      PARAMETER(lmax=4000.0E-9)
+
+      REAL albmax
+      PARAMETER(albmax=0.99)
+
+      REAL wavelengths(nlambda)
+      REAL ni(nlambda)
+
+      DATA wavelengths / 1.85000000e-07, 2.04170854e-07,
+     . 2.23341709e-07, 2.42512563e-07,
+     . 2.61683417e-07, 2.80854271e-07, 3.00025126e-07, 3.19195980e-07,
+     . 3.38366834e-07, 3.57537688e-07, 3.76708543e-07, 3.95879397e-07,
+     . 4.15050251e-07, 4.34221106e-07, 4.53391960e-07, 4.72562814e-07,
+     . 4.91733668e-07, 5.10904523e-07, 5.30075377e-07, 5.49246231e-07,
+     . 5.68417085e-07, 5.87587940e-07, 6.06758794e-07, 6.25929648e-07,
+     . 6.45100503e-07, 6.64271357e-07, 6.83442211e-07, 7.02613065e-07,
+     . 7.21783920e-07, 7.40954774e-07, 7.60125628e-07, 7.79296482e-07,
+     . 7.98467337e-07, 8.17638191e-07, 8.36809045e-07, 8.55979899e-07,
+     . 8.75150754e-07, 8.94321608e-07, 9.13492462e-07, 9.32663317e-07,
+     . 9.51834171e-07, 9.71005025e-07, 9.90175879e-07, 1.00934673e-06,
+     . 1.02851759e-06, 1.04768844e-06, 1.06685930e-06, 1.08603015e-06,
+     . 1.10520101e-06, 1.12437186e-06, 1.14354271e-06, 1.16271357e-06,
+     . 1.18188442e-06, 1.20105528e-06, 1.22022613e-06, 1.23939698e-06,
+     . 1.25856784e-06, 1.27773869e-06, 1.29690955e-06, 1.31608040e-06,
+     . 1.33525126e-06, 1.35442211e-06, 1.37359296e-06, 1.39276382e-06,
+     . 1.41193467e-06, 1.43110553e-06, 1.45027638e-06, 1.46944724e-06,
+     . 1.48861809e-06, 1.50778894e-06, 1.52695980e-06, 1.54613065e-06,
+     . 1.56530151e-06, 1.58447236e-06, 1.60364322e-06, 1.62281407e-06,
+     . 1.64198492e-06, 1.66115578e-06, 1.68032663e-06, 1.69949749e-06,
+     . 1.71866834e-06, 1.73783920e-06, 1.75701005e-06, 1.77618090e-06,
+     . 1.79535176e-06, 1.81452261e-06, 1.83369347e-06, 1.85286432e-06,
+     . 1.87203518e-06, 1.89120603e-06, 1.91037688e-06, 1.92954774e-06,
+     . 1.94871859e-06, 1.96788945e-06, 1.98706030e-06, 2.00623116e-06,
+     . 2.02540201e-06, 2.04457286e-06, 2.06374372e-06, 2.08291457e-06,
+     . 2.10208543e-06, 2.12125628e-06, 2.14042714e-06, 2.15959799e-06,
+     . 2.17876884e-06, 2.19793970e-06, 2.21711055e-06, 2.23628141e-06,
+     . 2.25545226e-06, 2.27462312e-06, 2.29379397e-06, 2.31296482e-06,
+     . 2.33213568e-06, 2.35130653e-06, 2.37047739e-06, 2.38964824e-06,
+     . 2.40881910e-06, 2.42798995e-06, 2.44716080e-06, 2.46633166e-06,
+     . 2.48550251e-06, 2.50467337e-06, 2.52384422e-06, 2.54301508e-06,
+     . 2.56218593e-06, 2.58135678e-06, 2.60052764e-06, 2.61969849e-06,
+     . 2.63886935e-06, 2.65804020e-06, 2.67721106e-06, 2.69638191e-06,
+     . 2.71555276e-06, 2.73472362e-06, 2.75389447e-06, 2.77306533e-06,
+     . 2.79223618e-06, 2.81140704e-06, 2.83057789e-06, 2.84974874e-06,
+     . 2.86891960e-06, 2.88809045e-06, 2.90726131e-06, 2.92643216e-06,
+     . 2.94560302e-06, 2.96477387e-06, 2.98394472e-06, 3.00311558e-06,
+     . 3.02228643e-06, 3.04145729e-06, 3.06062814e-06, 3.07979899e-06,
+     . 3.09896985e-06, 3.11814070e-06, 3.13731156e-06, 3.15648241e-06,
+     . 3.17565327e-06, 3.19482412e-06, 3.21399497e-06, 3.23316583e-06,
+     . 3.25233668e-06, 3.27150754e-06, 3.29067839e-06, 3.30984925e-06,
+     . 3.32902010e-06, 3.34819095e-06, 3.36736181e-06, 3.38653266e-06,
+     . 3.40570352e-06, 3.42487437e-06, 3.44404523e-06, 3.46321608e-06,
+     . 3.48238693e-06, 3.50155779e-06, 3.52072864e-06, 3.53989950e-06,
+     . 3.55907035e-06, 3.57824121e-06, 3.59741206e-06, 3.61658291e-06,
+     . 3.63575377e-06, 3.65492462e-06, 3.67409548e-06, 3.69326633e-06,
+     . 3.71243719e-06, 3.73160804e-06, 3.75077889e-06, 3.76994975e-06,
+     . 3.78912060e-06, 3.80829146e-06, 3.82746231e-06, 3.84663317e-06,
+     . 3.86580402e-06, 3.88497487e-06, 3.90414573e-06, 3.92331658e-06,
+     . 3.94248744e-06, 3.96165829e-06, 3.98082915e-06, 4.00000000e-06/
+
+
+      DATA ni /7.74508407e-10, 7.74508407e-10, 
+     .  7.74508407e-10, 7.74508407e-10,
+     .  7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 7.74508407e-10,
+     .  6.98381122e-10, 6.23170274e-10, 5.97655992e-10, 5.84106004e-10,
+     .  5.44327597e-10, 5.71923510e-10, 6.59723827e-10, 8.05183870e-10,
+     .  1.03110161e-09, 1.36680386e-09, 1.85161253e-09, 2.56487751e-09,
+     .  3.56462855e-09, 4.89450926e-09, 6.49252022e-09, 9.62029335e-09,
+     .  1.32335015e-08, 1.75502184e-08, 2.19240625e-08, 3.03304156e-08,
+     .  4.07715972e-08, 5.00414911e-08, 7.09722331e-08, 1.00773751e-07,
+     .  1.31427409e-07, 1.42289041e-07, 1.49066787e-07, 2.01558515e-07,
+     .  2.99106105e-07, 4.03902086e-07, 4.54292169e-07, 5.21906983e-07,
+     .  6.27643362e-07, 9.43955678e-07, 1.33464494e-06, 1.97278315e-06,
+     .  2.31801329e-06, 2.20584676e-06, 1.85568138e-06, 1.73395193e-06,
+     .  1.73101406e-06, 1.91333936e-06, 2.26413019e-06, 3.23959718e-06,
+     .  4.94316963e-06, 6.89378896e-06, 1.02237444e-05, 1.21439656e-05,
+     .  1.31567585e-05, 1.33448288e-05, 1.32000000e-05, 1.31608040e-05,
+     .  1.33048369e-05, 1.40321464e-05, 1.51526244e-05, 1.80342858e-05,
+     .  3.82875736e-05, 1.07325259e-04, 2.11961637e-04, 3.82008054e-04,
+     .  5.30897470e-04, 5.29244735e-04, 4.90876605e-04, 4.33905427e-04,
+     .  3.77795349e-04, 3.17633099e-04, 2.81078564e-04, 2.57579485e-04,
+     .  2.42203100e-04, 2.23789060e-04, 2.04306870e-04, 1.87909255e-04,
+     .  1.73117146e-04, 1.61533186e-04, 1.53420328e-04, 1.47578033e-04,
+     .  1.42334776e-04, 1.35691466e-04, 1.30495414e-04, 1.36065123e-04,
+     .  1.70928821e-04, 2.66389730e-04, 4.80957955e-04, 8.25041961e-04,
+     .  1.21654792e-03, 1.50232875e-03, 1.62316078e-03, 1.61649750e-03,
+     .  1.53736801e-03, 1.42343711e-03, 1.24459117e-03, 1.02388611e-03,
+     .  7.89112523e-04, 5.97204264e-04, 4.57152413e-04, 3.62341259e-04,
+     .  2.99128332e-04, 2.57035569e-04, 2.26992203e-04, 2.07110355e-04,
+     .  2.05835688e-04, 2.25108810e-04, 2.64262893e-04, 3.23594011e-04,
+     .  3.93061117e-04, 4.62789970e-04, 5.19664416e-04, 5.59739628e-04,
+     .  5.93476084e-04, 6.22797885e-04, 6.57484833e-04, 6.92849600e-04,
+     .  7.26584901e-04, 7.56604648e-04, 7.68009488e-04, 7.65579073e-04,
+     .  7.50526164e-04, 7.39809972e-04, 7.55622847e-04, 8.05099514e-04,
+     .  9.67279246e-04, 1.16281559e-03, 1.42570247e-03, 2.04986585e-03,
+     .  2.93971170e-03, 4.49827711e-03, 7.32537532e-03, 1.18889734e-02,
+     .  1.85851805e-02, 2.86242532e-02, 4.34131035e-02, 6.37828307e-02,
+     .  9.24145850e-02, 1.35547945e-01, 1.94143245e-01, 2.54542814e-01,
+     .  3.02282024e-01, 3.42214181e-01, 3.85475065e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
+     .  4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01/
+
+
+      Pas     =(lmax-lmin)/nlambda
+      ropt    = dopt/2.0
+      SSA     = 3.0/(rhoIce*ropt)
+      cos30   = cos(30.0/180.0*pi)
+
+
+      albint=0.
+      norm=0.
+        
+      DO i = 1,nlambda
+          gamm = 4.0 * pi * ni(i) / wavelengths(i)
+          cosalb = 2.0 / (SSA * rhoice) * B_amp * gamm
+          alpha = 16. / 3 * cosalb / (1.0 - g_asy)
+          R = exp(-(alpha**0.5) * 3.0 / 7.0 * (1.0 + 2.0 * cossza))
+          alb30 = exp(-(alpha**0.5)* 3.0 / 7.0 * (1.0 + 20 * cos30))
+
+          IF ((wavelengths(i).GE.lambdamin).AND.
+     .       (wavelengths(i).LT.lambdamax)) THEN
+          albint=albint+R*Pas  ! rectangle integration -> can be improved with trapezintegration
+          norm=norm+Pas
+          ENDIF
+
+      END DO
+
+      albint=max(0.,min(albint/max(norm,1E-10),albmax))
+
+      END
+  

LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_tso.F

@@ -132 +132 @@
 
       integer   nt_srf,it_srf,itEuBk          ! HL: Surface Scheme
-      parameter(nt_srf=10)                    !
+      parameter(nt_srf=10)                     ! 10 before
       real      agpsrf,xgpsrf,dt_srf,dt_ver   !
       real      etaBAK(knonv)                 !
@@ -153 +153 @@
 C +                                           ! including   Snow Melt  Energy
  
- 
+C +-- Initilialisation of local arrays
+C +   ================================
+        DO ikl=1,knonv
+
+          mu_sno(ikl)=0.
+          mu__dz(ikl,:)=0.    
+          dtC_sv(ikl,:)=0.      
+          IRs__D(ikl)=0.                 
+          dIRsdT(ikl)=0.                 
+          f_HSHL(ikl)=0.                 
+          dRidTs(ikl)=0.                 
+          HS___D(ikl)=0.                
+          f___HL(ikl)=0.             
+          HL___D(ikl)=0.            
+          TSurf0(ikl)=0.      
+          qsatsg(ikl)=0. 
+          dqs_dT(ikl)=0.                 
+          Psi(ikl)=0.               
+          RHuSol(ikl)=0.                 
+          Diag_A(ikl,:)=0.      
+          Diag_B(ikl,:)=0.      
+          Diag_C(ikl,:)=0.      
+          Term_D(ikl,:)=0.      
+          Aux__P(ikl,:)=0.      
+          Aux__Q(ikl,:)=0.      
+          etaBAK(ikl)=0.                
+          etaNEW(ikl)=0.               
+          etEuBk(ikl)=0.                 
+          fac_dt(ikl)=0. 
+          faceta(ikl)=0.  
+          PsiArg(ikl)=0. 
+          SHuSol(ikl)=0.  
+
+        END DO
+
+        
 
 C +--Heat Conduction Coefficient (zero in the Layers over the highest one)
@@ -336 +371 @@
 C +--Snow highest Layer (dummy!)
 C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
-        isl=  min(isnoSV(1)+1,nsno)
-        DO ikl=1,knonv
+
+        !EV!isl=  min(isnoSV(1)+1,nsno)
+
+        DO ikl=1,knonv
+! EV try to calculate isl at the ikl grid point
+          isl=  min(isnoSV(ikl)+1,nsno)
+
           Elem_A          =  dtC_sv(ikl,isl)  *mu__dz(ikl,isl)
           Elem_C          =  0.
@@ -384 +424 @@
 c    .                                   / den_qs              !
 c         qsatsg(ikl) = .0038 *        exp(arg_qs)             !
-
 !          sp = (pst_SV(ikl) + ptopSV) * 10.
 
-          sp=ps__SV(ikl)
+          !sp=ps__SV(ikl)
+          ! Etienne: in the formula herebelow sp should be in hPa, not
+          ! in Pa so I divide by 100.
+          sp=ps__SV(ikl)/100.
           psat_ice = 6.1070 * exp(6150. *(1./273.16 -
      .                                              1./TsisSV(ikl,isl)))
@@ -399 +441 @@
             qsatsg(ikl) = 0.622 * psat_wat / (sp - 0.378 * psat_wat)
           endif
+          QsT_SV(ikl)=qsatsg(ikl)
 
 c         dqs_dT(ikl) = qsatsg(ikl)* 4099.2   /(den_qs *den_qs)!
           fac_dt(ikl) = f_HSHL(ikl)/(ro_Wat   * dz_dSV(0))     !
         END DO
+
+
  
 C +--Surface: Latent    Heat Flux: Surface    Relative Humidity
@@ -410 +455 @@
      .                             /(   1.0-xgpsrf**nt_srf)    !
               dt_srf       = agpsrf                            !
-              dt_ver       = 0.                                !
+              dt_ver       = 0.                
+
             DO ikl=1,knonv
-              isl          =          isnoSV(ikl)              !
+              isl          =          isnoSV(ikl)              
+              ist          = max(0,isotSV(ikl)-100*isnoSV(ikl))! 0 if    H2O                          
+              ist__s       = min(1,ist)                                                                             
               etaBAK(ikl)  = max(epsi,eta_SV(ikl ,isl))        !
               etaNEW(ikl)  =          etaBAK(ikl)              !
               etEuBk(ikl)  =          etaNEW(ikl)              !
-            END DO                                             !
+            END DO      
+
+        if(ist__s==1) then ! to reduce computer time                                                  
+                                          !
         DO it_srf=1,nt_srf                                     !
               dt_ver       = dt_ver     +dt_srf                !
@@ -458 +509 @@
             END DO                                             !
               dt_srf      =      dt_srf         * xgpsrf       !
-        END DO                                                 !
+        END DO          
+
+
+        endif                                       !
  
 C +--Surface: Latent    Heat Flux: Soil/Water Surface Contributions
@@ -579 +633 @@
        
       END DO
+
+
  
 C +--Temperature Limits (avoids problems in case of no Snow Layers)
@@ -584 +640 @@
         DO ikl=     1,knonv
            isl              = isnoSV(ikl)
-          dTSurf            = TsisSV(ikl,isl) -     TSurf0(ikl)
+
+           dTSurf            = TsisSV(ikl,isl) -     TSurf0(ikl)
           TsisSV(ikl,isl)   = TSurf0(ikl) + sign(1.,dTSurf) ! 180.0 dgC/hr
      .              * min(abs(dTSurf),5.e-2*dt__SV)         ! =0.05 dgC/s
@@ -602 +659 @@
 C +--Update Surface    Fluxes
 C +  ========================
- 
+        
+
+
         DO ikl=      1,knonv
           isl         = isnoSV(ikl)
@@ -613 +672 @@
         END DO
 
- 
- 
       return
       end

LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_zsn.F

@@ -52 +52 @@
       use VARxSV
       use VARySV
+      use surface_data, only: ok_zsn_ii
 
       IMPLICIT NONE
@@ -716 +717 @@
         END DO
 
+
+C +--Search new Ice/Snow Interface (option II in MAR)
+C +  ===============================================
+
+        IF (ok_zsn_ii) THEN 
+        
+        DO ikl=1,knonv
+          iiceSV(ikl) =  0
+        END DO
+ 
+        DO ikl=1,knonv
+        DO   isn=1,isnoSV(ikl)
+          OK_ICE      = max(zero,sign(unun,ro__SV(ikl,isn)-ro_ice+20.))
+     .                * max(zero,sign(unun,dzsnSV(ikl,isn)-epsi))
+          iiceSV(ikl) = (1.-OK_ICE)       *iiceSV(ikl)
+     .                +     OK_ICE        *isn
+        END DO
+        END DO
+
+        END IF
  
       return

LMDZ6/trunk/libf/phylmd/inlandsis/surf_inlandsis_mod.F90

@@ -1 +1 @@
 MODULE surf_inlandsis_mod
 
-  IMPLICIT NONE
-
+    IMPLICIT NONE
+
+CONTAINS
+
+
+SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, &
+            rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
+            precip_rain, precip_snow, &
+            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
+            rugos, snow_cont_air, alb_soil, alt, slope, cloudf, &
+            radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
+            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
+            runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, &
+            tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf)
+
+        ! |                                                                        |
+        ! |   SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow|
+        ! |                              (INLANDSIS)                               |
+        ! |   SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)         |
+        ! |   surface scheme of the Modele Atmospherique Regional (MAR)            |
+        ! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
+        ! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
+        ! |   Updated by Etienne Vignon, Cecile Agosta                             |
+        ! |                                                                        |
+        ! +------------------------------------------------------------------------+
+        ! |
+        ! |   In the current setup, SISVAT is used only to model the land ice      |
+        ! |   part of the surface; hence it is called with the compressed variables|
+        ! |   from pbl_surface, and only by the surf_landice routine.              |
+        ! |                                                                        |
+        ! |   In this interface it is assumed that the partitioning of the soil,   |
+        ! |   and hence the number of grid points is constant during a simulation, |
+        ! |   hence eg. snow properties remain stored in the global SISVAT         |
+        ! |   variables between the calls and don't need to be handed over as      |
+        ! |   arguments. When the partitioning is supposed to change, make sure to |
+        ! |   update the variables.                                                |
+        ! |                                                                        |
+        ! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV ...)                    |
+        ! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |
+        ! |                                                                        |
+        ! +------------------------------------------------------------------------+
+
+        USE dimphy
+        USE VAR_SV
+        USE VARdSV
+        USE VARxSV
+        USE VARySV
+        USE VARtSV
+        USE VARphy
+        USE surface_data, only : iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor
+
+        IMPLICIT NONE
+
+        ! +--INTERFACE Variables
+        ! +  ===================
+        !    include  "dimsoil.h"
+
+        ! +--Global Variables
+        ! +  ================
+        ! Input Variables for SISVAT
+        INTEGER, INTENT(IN) :: knon
+        INTEGER, INTENT(IN) :: itime
+        REAL, INTENT(IN) :: dtime
+        LOGICAL, INTENT(IN) :: debut     ! true if first step
+        LOGICAL, INTENT(IN) :: lafin     ! true if last step
+
+        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i     ! Index Decompression
+        REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
+        REAL, DIMENSION(klon), INTENT(IN) :: rmu0      ! cos sol. zenith angle
+        REAL, DIMENSION(klon), INTENT(IN) :: swdown    !
+        REAL, DIMENSION(klon), INTENT(IN) :: lwdown    !
+        REAL, DIMENSION(klon), INTENT(IN) :: albedo_old
+        REAL, DIMENSION(klon), INTENT(IN) :: pexner    ! Exner potential
+        REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
+        REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo
+        REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay
+        REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf
+        REAL, DIMENSION(klon), INTENT(IN) :: rugos
+        REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air
+        REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope
+        REAL, DIMENSION(klon), INTENT(IN) :: alt       ! surface elevation
+        REAL, DIMENSION(klon), INTENT(IN) :: cloudf
+        REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
+        REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
+        REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh
+        REAL, DIMENSION(klon), INTENT(IN) :: ustar   ! friction velocity
+
+        ! Variables exchanged between LMDZ and SISVAT
+        REAL, DIMENSION(klon), INTENT(IN) :: radsol    ! Surface absorbed rad.
+        REAL, DIMENSION(klon), INTENT(INOUT) :: snow      ! Tot snow mass [kg/m2]
+        REAL, DIMENSION(klon), INTENT(INOUT) :: zfra      ! snwo surface fraction [0-1]
+        REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
+        REAL, DIMENSION(klon), INTENT(OUT) :: qsol      ! Soil Water Content
+        REAL, DIMENSION(klon), INTENT(INOUT) :: z0m    ! Momentum Roughn Lgt
+        REAL, DIMENSION(klon), INTENT(INOUT) :: z0h    ! Momentum Roughn Lgt
+
+        ! Output Variables for LMDZ
+        REAL, DIMENSION(klon), INTENT(OUT) :: alb1      ! Albedo SW
+        REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW
+        REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo
+        REAL, DIMENSION(klon), INTENT(OUT) :: emis_new  ! Surface Emissivity
+        REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff
+        REAL, DIMENSION(klon), INTENT(OUT) :: ffonte    ! enthalpy flux due to surface melting
+        REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte   ! water flux due to surface melting
+        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s   ! d/dT sens. ht flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l   ! d/dT latent ht flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens  ! Sensible ht flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat   ! Latent heat flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: evap      ! Evaporation
+        REAL, DIMENSION(klon), INTENT(OUT) :: erod      ! Erosion of surface snow (flux)
+        REAL, DIMENSION(klon), INTENT(OUT) :: agesno    ! Snow age (top layer)
+        REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature
+        REAL, DIMENSION(klon), INTENT(OUT) :: qsurf     ! Surface Humidity
+
+        ! Specific INLANDIS outputs
+        REAL, DIMENSION(klon), INTENT(OUT) :: qsnow     ! Total H2O snow[kg/m2]
+        REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt   ! Snow height (m)
+        REAL, DIMENSION(klon), INTENT(OUT) :: to_ice    ! Snow passed to ice
+        REAL, DIMENSION(klon), INTENT(OUT) :: sissnow   ! Snow in model (kg/m2)
+
+        ! +--Internal  Variables
+        ! +  ===================
+
+        CHARACTER(len = 20) :: fn_outfor ! Name for output file
+        CHARACTER (len = 80)              :: abort_message
+        CHARACTER (len = 20)              :: modname = 'surf_inlandsis_mod'
+
+        INTEGER :: i, ig, ikl, isl, isn, nt
+        INTEGER :: gp_outfor, un_outfor
+        REAL, PARAMETER :: f1 = 0.5
+        REAL, PARAMETER :: sn_upp = 10000., sn_low = 500.
+        REAL, PARAMETER :: sn_add = 400., sn_div = 2.
+        ! snow mass upper,lower limit,
+        ! added mass/division lowest layer
+        REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6
+        REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3
+        ! Parameters for drainage
+        ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
+        ! +...        Run Off Parameters
+        ! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)
+        ! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)
+
+        REAL, DIMENSION(klon) :: eps0SL          ! surface Emissivity
+        REAL :: zsigma, Ua_min, Us_min, lati
+        REAL, PARAMETER :: cdmax=0.05
+        REAL :: lambda          ! Par. soil discret.
+        REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2         ! Soil layer thicknesses
+        !$OMP THREADPRIVATE(dz1,dz2)
+        LOGICAL, SAVE :: firstcall
+        !$OMP THREADPRIVATE(firstcall)
+
+        INTEGER :: iso
+        LOGICAL :: file_exists
+        CHARACTER(len = 20) :: fichnom
+        LOGICAL :: is_init_domec
+        ! CA initialization
+        ! dz_profil_15 : 1 m in 15 layers [m]
+        real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, &
+                                                0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/)
+        ! mean_temp : mean annual surface temperature [K]
+        real, dimension(klon) :: mean_temp
+        ! mean_dens : mean surface density [kg/m3]
+        real, dimension(klon) :: mean_dens
+        ! lat_scale : temperature lapse rate against latitude [K degree-1]
+        real :: lat_scale
+        ! sh_scale : temperature lapse rate against altitude [K km-1]
+        real :: sh_scale
+        ! variables for density profile
+        ! E0, E1 : exponent
+        real :: E0, E1
+        ! depth at which 550 kg m-3 is reached [m]
+        real :: z550
+        ! depths of snow layers
+        real :: depth, snow_depth, distup
+        ! number of initial snow layers
+        integer :: nb_snow_layer
+        ! For density calc.
+        real :: alpha0, alpha1, ln_smb
+        ! theoritical densities [kg m-3]
+        real :: rho0, rho1, rho1_550
+        ! constants for density profile
+        ! C0, C1 : constant, 0.07 for z <= 550 kg m-3
+        real, parameter :: C0 = 0.07
+        real, parameter :: C1 = 0.03
+        ! rho_i : ice density [kg m-3]
+        real, parameter :: rho_ice = 917.
+        ! E_c : activation energy [J mol-1]
+        real, parameter :: E_c = 60000.
+        ! E_g : activation energy [J mol-1]
+        real, parameter :: E_g = 42400.
+        ! R : gas constant [J mol-1 K-1]
+        real, parameter :: R = 8.3144621
+
+      
+     
+
+
+        ! + PROGRAM START
+        ! + -----------------------------------------
+
+        zsigma = 1000.
+        dt__SV = dtime
+
+        IF (debut) THEN
+            firstcall = .TRUE.
+            INI_SV = .false.
+        ELSE
+            firstcall = .false.
+            INI_SV = .true.
+        END IF
+
+        IF (ok_outfor) THEN
+            un_outfor = 51    ! unit number for point output file
+            gp_outfor = 1    ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC
+            fn_outfor = 'outfor_SV.dat'
+        END IF
+
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! + INITIALISATION: BEGIN +++
+        ! + -----------------------------------------
+        IF (firstcall) THEN
+
+            ! +--Array size
+            ! +  -----------------------
+
+            klonv = klon
+            knonv = knon
+                write(*, *) 'ikl, lon and lat in INLANDSIS'
+
+            DO ikl = 1, knon
+                i=ikl2i(ikl)
+                write(*, *) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i)
+            END DO
+
+            ! +--Variables initizialisation
+            ! +  ---------------------------
+
+            CALL INIT_VARtSV
+            CALL INIT_VARxSV
+            CALL INIT_VARySV
+
+
+
+            ! +--Surface Fall Line Slope
+            ! +  -----------------------
+            IF (SnoMod)  THEN
+                DO ikl = 1, knon
+                    slopSV(ikl) = slope(ikl)
+                    SWf_SV(ikl) = &   ! Normalized Decay of the
+                            exp(-dt__SV             &   ! Surficial Water Content
+                                    / (c1_zuo                &   !(Zuo and Oerlemans 1996,
+                                            + c2_zuo * exp(-c3_zuo * abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
+                END DO
+            END IF
+
+
+
+            ! +--Soil layer thickness . Compute soil discretization (as for LMDZ)
+            ! +  ----------------------------------------------------------------
+            !        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
+            CALL get_soil_levels(dz1, dz2, lambda)
+
+            lambSV = lambda
+            dz1_SV(1:knon, 1:) = 0.
+            dz2_SV(1:knon, 1:) = 0.
+
+            DO isl = -nsol, 0
+                dz_dSV(isl) = 0.5e-3 * dz2(1 - isl)           ! Soil layer thickness
+                DO ikl = 1, knon
+                    dz1_SV(ikl, isl) = dz1(1 - isl)    !1.e-3*
+                    dz2_SV(ikl, isl) = dz2(1 - isl)    !1.e-3*
+                END DO
+            END DO
+
+
+            ! Set variables
+            ! =============
+            DO ikl = 1, knon
+                ! LSmask : Land/Sea Mask
+                LSmask(ikl) = 1
+                ! isotSV : Soil Type -> 12 = ice
+                isotSV(ikl) = 12
+                ! iWaFSV : Soil Drainage (1,0)=(y,n)
+                iWaFSV(ikl) = 1
+                ! eps0SL : Surface Emissivity
+                eps0SL(ikl) = 1.
+                ! alb0SV : Soil Albedo
+                alb0SV(ikl) = alb_soil(ikl)
+                ! Tsf_SV : Surface Temperature, must be bellow freezing
+                Tsf_SV(ikl) = min(temp_air(ikl), TfSnow)
+            END DO
+
+            ! +--Initialization of soil and snow variables in case startsis is not read
+            ! +  ----------------------------------------------------------------------
+
+            is_init_domec=.FALSE.
+
+
+            IF (is_init_domec) THEN
+            ! Coarse initilization inspired from vertcical profiles at Dome C,
+            ! Antarctic Plateaui (10m of snow, 19 levels)
+
+                 DO ikl = 1,knon
+! + Soil
+                 DO isl =   -nsol,0    
+                   TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2)   !temp_air(ikl)
+                   !tsoil(ikl,1-isl)   Soil Temperature
+                   !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
+                   eta_SV(ikl,isl) = epsi           !etasoil(ikl,1-isl) Soil Water[m3/m3]
+                   ro__SV(ikl,isl) = rhoIce         !rosoil(ikl,1-isl) volumic mass
+                 END DO   
+
+
+           ! Snow
+                 isnoSV(ikl) = 19
+                 istoSV(ikl, 1:isnoSV(ikl)) = 100
+                 ro__SV(ikl, 1:isnoSV(ikl)) = 350.
+                 eta_SV(ikl, 1:isnoSV(ikl)) = epsi
+                 TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2)
+                 G1snSV(ikl, 1:isnoSV(ikl)) = 99.
+                 G2snSV(ikl, 1:isnoSV(ikl)) = 2.
+                 agsnSV(ikl, 1:isnoSV(ikl)) = 50.
+                 dzsnSV(ikl, 19) = 0.015
+                 dzsnSV(ikl, 18) = 0.015
+                 dzsnSV(ikl, 17) = 0.020
+                 dzsnSV(ikl, 16) = 0.030
+                 dzsnSV(ikl, 15) = 0.040
+                 dzsnSV(ikl, 14) = 0.060
+                 dzsnSV(ikl, 13) = 0.080
+                 dzsnSV(ikl, 12) = 0.110
+                 dzsnSV(ikl, 11) = 0.150
+                 dzsnSV(ikl, 10) = 0.200
+                 dzsnSV(ikl, 9) = 0.300
+                 dzsnSV(ikl, 8) = 0.420
+                 dzsnSV(ikl, 7) = 0.780
+                 dzsnSV(ikl, 6) = 1.020
+                 dzsnSV(ikl, 5) = 0.980
+                 dzsnSV(ikl, 4) = 1.020
+                 dzsnSV(ikl, 3) = 3.970
+                 dzsnSV(ikl, 2) = 1.020
+                 dzsnSV(ikl, 1) = 1.020
+
+                 END DO
+            ELSE
+
+            ! Initilialisation with climatological temperature and density
+            ! profiles as in MAR. Methodology developed by Cecile Agosta
+  
+            ! initialize with 0., for unused snow layers
+            dzsnSV = 0.
+            G1snSV = 0.
+            G2snSV = 0.
+            istoSV = 0
+            TsisSV = 0.
+
+
+            ! initialize mean variables (unrealistic)
+            mean_temp = TfSnow
+            mean_dens = 300.
+            ! loop on grid cells
+            DO ikl = 1, knon
+                lati=rlat(ikl2i(ikl))
+                ! approximations for mean_temp and mean_dens
+                ! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1)
+                ! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1
+                ! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica,
+                ! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed.
+                if (lati > 60.) then
+                    ! CA todo : add longitude bounds
+                    ! Greenland mean temperature : function of altitude and latitude
+                    ! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1
+                    lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9
+                    lat_scale = max(min(lat_scale, 0.9), 0.75)
+                    ! sh_scale equals the environmental lapse rate : 6.5 °C km-1
+                    sh_scale = 6.5
+                    mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.)
+                    ! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051
+                    mean_dens(ikl) = 315.
+                else if (lati < -60.) then
+                    ! Antarctica mean temperature : function of altitude and latitude
+                    ! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1
+                    lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3
+                    lat_scale = max(min(lat_scale, 1.3), 0.6)
+                    ! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1
+                    sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5
+                    sh_scale = max(min(sh_scale, 9.8), 6.5)
+                    mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.)
+                    ! Antarctica surface density : function of mean annual temperature
+                    ! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013)
+                    ! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m.
+                    ! Weinhart et al 2020  https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201
+                    ! 320 kg m-3 is reached at Dome A, 4100 m a.s.l.
+                    ! Dome C : st_ant_param(3233, -75.1) = -47.7
+                    ! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7
+                    mean_dens(ikl) =  (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450.
+                    mean_dens(ikl) = min(450., max(320., mean_dens(ikl)))
+                else
+
+                !    write(*, *) 'Attention: temperature initialization is only defined for Greenland and Antarctica'
+
+                     mean_dens(ikl) =350.
+                     mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2)
+
+                !abort_message='temperature initialization is only defined for Greenland and Antarctica'
+                !CALL abort_physic(modname,abort_message,1)
+
+                end if
  
-  CONTAINS
-
-                                       
-
-  SUBROUTINE surf_inlandsis(knon,rlon,rlat, ikl2i, itime, dtime, debut, lafin, &
-             rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
-             precip_rain, precip_snow, precip_snow_adv, snow_adv, &
-             zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
-             rugos, snow_cont_air, alb_soil, slope, cloudf, &
-             radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, & 
-             AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
-             runoff_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &      
-             tsurf_new, alb1, alb2, alb3, &
-             emis_new, z0m, z0h, qsurf)       
-                                                                             
-! +------------------------------------------------------------------------+   
-! |                                                                        |   
-! |   SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow|
-! |                              (INLANDSIS)                               |
-! |   SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)         |   
-! |   surface scheme of the Modele Atmospherique Regional (MAR)            |
-! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
-! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
-! |           Update Etienne Vignon, LMD, Novembre 2020                    |
-! |                                                                        |   
-! +------------------------------------------------------------------------+   
-! |   
-! |   In the current setup, SISVAT is used only to model the land ice      |
-! |   part of the surface; hence it is called with the compressed variables| 
-! |   from pbl_surface, and only by the surf_landice routine.              |
-! |                                                                        |   
-! |   In this interface it is assumed that the partitioning of the soil,   |
-! |   and hence the number of grid points is constant during a simulation, |
-! |   hence eg. snow properties remain stored in the global SISVAT         |
-! |   variables between the calls and don't need to be handed over as      |
-! |   arguments. When the partitioning is supposed to change, make sure to |
-! |   update the variables.                                                | 
-! |                                                                        |  
-! |   INPUT                                                                |  
-! |             SnoMod: Snow Pack is set up when .T.                       |  
-! |             reaLBC: Update Bound.Condit.when .T.                       |   
-! |                                                                        |  
-! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV)                        |  
-! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |   
-! |                                                                        |  
-! |   Preprocessing  Option: SISVAT PHYSICS                                |   
-! |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^^^^^^                                |   
-! | #                       #HY                                            |   
-! | #                       #SN: Snow         Model                        |   
-! | #                       #BS: Blowing Snow Parameterization             |   
-! +------------------------------------------------------------------------+ 
-        
-    USE dimphy                                            
-    USE VAR_SV 
-    USE VARdSV 
-    USE VARxSV 
-    USE VARySV
-    USE VARtSV
-    USE VARphy
-    USE surface_data, only: iflag_tsurf_inlandsis,SnoMod,BloMod,ok_outfor                                       
-
-    IMPLICIT NONE    
-                                                                
-! +--INTERFACE Variables                                                     
-! +  =================== 
-                                            
-!    include  "dimsoil.h"                                     
-                                                                    
-
-! +--Global Variables                                                          
-! +  ================  
-! Input Variables for SISVAT
-    INTEGER,               INTENT(IN)      :: knon
-    INTEGER,               INTENT(IN)      :: itime   
-    REAL,                  INTENT(IN)      :: dtime 
-    LOGICAL,               INTENT(IN)      :: debut     ! true if first step
-    LOGICAL,               INTENT(IN)      :: lafin     ! true if last step
-
-    INTEGER, DIMENSION(klon), INTENT(IN)   :: ikl2i     ! Index Decompression
-    REAL, DIMENSION(klon), INTENT(IN)      :: rlon, rlat
-    REAL, DIMENSION(klon), INTENT(IN)      :: rmu0      ! cos sol. zenith angle
-    REAL, DIMENSION(klon), INTENT(IN)      :: swdown    !
-    REAL, DIMENSION(klon), INTENT(IN)      :: lwdown    ! 
-    REAL, DIMENSION(klon), INTENT(IN)      :: albedo_old    
-    REAL, DIMENSION(klon), INTENT(IN)      :: pexner    ! Exner potential
-    REAL, DIMENSION(klon), INTENT(IN)      :: precip_rain, precip_snow
-    REAL, DIMENSION(klon), INTENT(IN)      :: precip_snow_adv, snow_adv 
-                                                        !Snow Drift
-    REAL, DIMENSION(klon), INTENT(IN)      :: zsl_height, wind_velo
-    REAL, DIMENSION(klon), INTENT(IN)      :: temp_air, spechum, ps,p1lay
-    REAL, DIMENSION(klon), INTENT(IN)      :: dens_air, tsurf            
-    REAL, DIMENSION(klon), INTENT(IN)      :: rugos,snow_cont_air
-    REAL, DIMENSION(klon), INTENT(IN)      :: alb_soil, slope 
-    REAL, DIMENSION(klon), INTENT(IN)      :: cloudf   
-    REAL, DIMENSION(klon), INTENT(IN)      :: AcoefH, AcoefQ
-    REAL, DIMENSION(klon), INTENT(IN)      :: BcoefH, BcoefQ
-    REAL, DIMENSION(klon), INTENT(IN)      :: cdragm, cdragh
-    REAL, DIMENSION(klon), INTENT(IN)      :: ustar   ! friction velocity
-
-! Variables exchanged between LMDZ and SISVAT
-    REAL, DIMENSION(klon), INTENT(IN)      :: radsol    ! Surface absorbed rad.
-    REAL, DIMENSION(klon), INTENT(INOUT)   :: snow      ! Tot snow mass [kg/m2]
-    REAL, DIMENSION(klon), INTENT(INOUT)   :: zfra      ! snwo surface fraction [0-1]
-    REAL, DIMENSION(klon,nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature 
-    REAL, DIMENSION(klon), INTENT(OUT)       :: qsol      ! Soil Water Content  
-    REAL, DIMENSION(klon), INTENT(INOUT)     :: z0m    ! Momentum Roughn Lgt
-    REAL, DIMENSION(klon), INTENT(INOUT)     :: z0h    ! Momentum Roughn Lgt 
-
-
-! Output Variables for LMDZ
-    REAL, DIMENSION(klon), INTENT(OUT)     :: alb1      ! Albedo SW 
-    REAL, DIMENSION(klon), INTENT(OUT)     :: alb2,alb3 ! Albedo NIR and LW
-    REAL, DIMENSION(klon), INTENT(OUT)     :: emis_new  ! Surface Emissivity
-    REAL, DIMENSION(klon), INTENT(OUT)     :: runoff_lic ! Runoff            
-    REAL, DIMENSION(klon), INTENT(OUT)     :: dflux_s   ! d/dT sens. ht flux  
-    REAL, DIMENSION(klon), INTENT(OUT)     :: dflux_l   ! d/dT latent ht flux 
-    REAL, DIMENSION(klon), INTENT(OUT)     :: fluxsens  ! Sensible ht flux    
-    REAL, DIMENSION(klon), INTENT(OUT)     :: fluxlat   ! Latent heat flux
-    REAL, DIMENSION(klon), INTENT(OUT)     :: evap      ! Evaporation
-    REAL, DIMENSION(klon), INTENT(OUT)     :: agesno    ! Snow age (top layer)
-    REAL, DIMENSION(klon), INTENT(OUT)     :: tsurf_new ! Surface Temperature 
-    REAL, DIMENSION(klon), INTENT(OUT)     :: qsurf     ! Surface Humidity 
-
-! Specific INLANDIS outputs
-
-    REAL, DIMENSION(klon), INTENT(OUT)     :: qsnow     ! Total H2O snow[kg/m2]
-    REAL, DIMENSION(klon), INTENT(OUT)     :: snowhgt   ! Snow height (m)
-    REAL, DIMENSION(klon), INTENT(OUT)     :: to_ice    ! Snow passed to ice 
-    REAL, DIMENSION(klon), INTENT(OUT)     :: sissnow   ! Snow in model (kg/m2)
-
-                                                                          
-
-                                                                    
-! +--Internal  Variables                                                       
-! +  ===================                                          
-
-    CHARACTER(len=20)               :: fn_outfor ! Name for output file
-    INTEGER                         :: i, ig, ikl, isl, isn, nt
-    INTEGER                         :: gp_outfor, un_outfor
-    REAL, PARAMETER                 :: f1=0.5
-    REAL, PARAMETER                 :: sn_upp=5000.,sn_low=500.
-    REAL, PARAMETER                 :: sn_add=400.,sn_div=2.
-                                             ! snow mass upper,lower limit, 
-                                             ! added mass/division lowest layer
-    REAL, PARAMETER                 :: c1_zuo=12.960e+4, c2_zuo=2.160e+6
-    REAL, PARAMETER                 :: c3_zuo=1.400e+2,  czemin=1.e-3  
-                                             ! Parameters for drainage
-! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
-! +...        Run Off Parameters                                              
-! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)    
-! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)             
-
-    REAL, DIMENSION(klon)           :: eps0SL          ! surface Emissivity 
-    REAL                            :: zsigma, Ua_min, Us_min
-    REAL                            :: lambda          ! Par. soil discret.
-    REAL, DIMENSION(nsoilmx), SAVE  :: dz1,dz2         ! Soil layer thicknesses
-!$OMP THREADPRIVATE(dz1,dz2)
-    LOGICAL, SAVE                   :: firstcall
-!$OMP THREADPRIVATE(firstcall)                
-
-               
-                                        
-! +--Internal Variables
-! +  ==================
-
-    INTEGER                         ::  iso
-    LOGICAL                         ::  file_exists
-    CHARACTER(len=20)               ::  fichnom
-!========================================================================
-
-      PRINT*, 'je rentre dans inlandsis'
-
-      zsigma=1000.
-      dt__SV=dtime
-      
-    
-
-!     write(*,*)'Start of simulation? ',debut        !hj
-
-      IF (debut) THEN 
-        firstcall=.TRUE. 
-        INI_SV=.false.
-
-      ELSE
-        firstcall=.false.
-        INI_SV=.true.
-      END IF
-
-
-
-
-       IF (ok_outfor) THEN
-        un_outfor=51                 ! unit number for point output file
-        gp_outfor= 1        ! grid point number for point output
-        fn_outfor='outfor_SV.dat'
-       END IF
-
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-
-
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! + INITIALISATION: BEGIN +++ 
-! + -------------------------
-! +
-! + Compute soil discretization (as for LMDZ)
-! + -----------------------------------------  
-      IF (firstcall) THEN
-
-! +--Array size
-     klonv=klon
-     knonv=knon
-
-
-        write(*,*)'klon',klon,'klonv',klonv,'knon',knon,'nsol',nsol,'nsno',nsno
-
-
-        CALL INIT_VARtSV
-        CALL INIT_VARxSV
-        CALL INIT_VARySV
-
-        eps0SL(:)=0.
-        
-        
-! +--Soil layer thickness                                                   
-! +  -----------------------  
-!        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
-        CALL get_soil_levels(dz1,dz2,lambda)
-
-        
-        lambSV=lambda
-        dz1_SV(1:knon,1:) = 0.     
-        dz2_SV(1:knon,1:) = 0.
-                      
-        DO isl =   -nsol,0   
-          dz_dSV(isl) = 0.5e-3*dz2(1-isl)           ! Soil layer thickness 
-          DO ikl=1,knon
-            dz1_SV(ikl,isl) = dz1(1-isl)    !1.e-3* 
-            dz2_SV(ikl,isl) = dz2(1-isl)    !1.e-3*
-          END DO
+                ! mean_temp is defined for ice ground only
+                mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2)
+
+                ! Soil layers
+                ! ===========
+                DO isl = -nsol, 0
+                    ! TsisSV : Temperature [K]
+                    TsisSV(ikl, isl) = mean_temp(ikl)
+                    ! eta_SV : Soil Water [m3/m3]
+                    eta_SV(ikl, isl) = epsi
+                    ! ro__SV : Volumic Mass [kg/m3]
+                    ro__SV(ikl, isl) = rhoIce
+                END DO
+
+                ! Snow layers
+                ! ===========
+                ! snow_depth : initial snow depth
+                snow_depth = 20.
+                ! nb_snow_layer : initial nb of snow layers
+                nb_snow_layer = 15
+                ! isnoSV : total nb of snow layers
+                isnoSV(ikl) = nb_snow_layer
+                ! depth : depth of each layer
+                depth = snow_depth
+                do isl = 1, nb_snow_layer
+                    ! dzsnSV : snow layer thickness
+                    dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1))
+                    ! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical
+                    G1snSV(ikl, isl) = 99.
+                    ! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size
+                    G2snSV(ikl, isl) = 3.
+                    agsnSV(ikl, isl) = 0.
+                    istoSV(ikl, isl) = 0
+                    ! eta_SV : Liquid Water Content [m3/m3]
+                    eta_SV(ikl, isl) = 0.
+                    ! distance to surface
+                    depth = depth - dzsnSV(ikl,isl) / 2.
+                    distup = min(1., max(0., depth / snow_depth))
+                    ! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom)
+                    TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup**2) + mean_temp(ikl) * distup**2
+                    ! firn density : densification formulas from :
+                    ! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/)
+                    ! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306)
+                    ! Integration of the steady state equation
+                    ! ln_smb approximated as a function of temperature
+                    ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.)
+                    ! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1
+                    alpha0 = max(1.435 - 0.151 * ln_smb, 0.25)
+                    alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25)
+                    E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0
+                    E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1
+                    z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0
+                    rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
+                    rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice
+                    if (depth <= z550) then
+                        ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
+                    else
+                        ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice
+                    end if
+                    depth = depth - dzsnSV(ikl,isl) / 2.
+                    
+                end do
+
+            END DO
+
+            END IF
+
+
+            ! + Numerics paramaters, SISVAT_ini 
+            ! +  ----------------------
+            CALL SISVAT_ini(knon)
+
+
+            ! +--Read restart file
+            ! +  =================================================
+
+            INQUIRE(FILE = "startsis.nc", EXIST = file_exists)
+            IF (file_exists) THEN
+                CALL sisvatetat0("startsis.nc", ikl2i)
+            END IF
+
+
+
+            ! +--Output ascii file
+            ! +  =================================================
+
+            ! open output file
+            IF (ok_outfor) THEN
+                open(unit = un_outfor, status = 'replace', file = fn_outfor)
+                ikl = gp_outfor     ! index sur la grille land ice
+                write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl))
+                write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] &
+                        & G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]'
+            END IF
+
+        END IF  ! firstcall
+        ! +
+        ! +  +++  INITIALISATION:  END  +++
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! + READ FORCINGS
+        ! + ------------------------
+
+        ! + Update Forcings for SISVAT given by the LMDZ model.
+        ! +
+        DO ikl = 1, knon
+
+            ! +--Atmospheric Forcing                                    (INPUT)
+            ! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^
+            za__SV(ikl) = zsl_height(ikl)               ! surface layer height (fisr model level) [m]
+            Ua_min = 0.2 * sqrt(za__SV(ikl))            !
+            VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
+            TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
+            ExnrSV(ikl) = pexner(ikl)                   ! Exner potential
+            rhT_SV(ikl) = dens_air(ikl)                 ! Air density
+            QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity
+            ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
+            p1l_SV(ikl) = p1lay(ikl)                    ! lowest atm. layer press[Pa]
+
+            ! +--Surface properties
+            ! +  ^^^^^^^^^^^^^^^^^^
+
+            Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
+            Z0h_SV(ikl) = z0h(ikl)                      ! Moment.Roughn.L.
+
+            ! +--Energy Fluxes                                          (INPUT)
+            ! +  ^^^^^^^^^^^^^                                           ^^^^^
+            coszSV(ikl) = max(czemin, rmu0(ikl))         ! cos(zenith.Dist.)
+            sol_SV(ikl) = swdown(ikl)                   ! downward Solar
+            IRd_SV(ikl) = lwdown(ikl)                   ! downward IR
+            rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.
+
+            ! +--Water  Fluxes                                          (INPUT)
+            ! +  ^^^^^^^^^^^^^                                           ^^^^^
+            drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
+            dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]
+
+            ! #BS    dbs_SV(ikl) = blowSN(i,j,n)
+            ! dbs_SV = Maximum potential erosion amount [kg/m2]
+            ! => Upper bound for eroded snow mass
+            !        uss_SV(ikl) = SLussl(i,j,n) ! u*qs* (only for Tv in sisvatesbl.f)
+            ! #BS  if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then
+            ! #BS    dsnbSV(ikl) =1.0-min(qsHY(i,j,kB)     !BS neglib. at kb ~100 magl)
+            ! #BS.                        /max(qshy(i,j,mz),eps9),unun)
+            ! #BS    dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl))
+            ! #BS    dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl)))
+            ! #BS  else
+            ! #BS    dsnbSV(ikl) = 0.
+            ! #BS  endif
+            !      dsnbSV is the drift fraction of deposited snow updated in sisvat.f
+            !      will be used for characterizing the Buffer Layer
+            !      (see update of  Bros_N, G1same, G2same, zroOLD, zroNEW)
+            ! #BS  if(n==1) qbs_HY(i,j) = dsnbSV(ikl)
+            qsnoSV(ikl) = snow_cont_air(ikl)
+
+
+
+            ! +--Soil/BL                                      (INPUT)
+            ! +  ^^^^^^^                                       ^^^^^
+            alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo
+            AcoHSV(ikl) = AcoefH(ikl)
+            BcoHSV(ikl) = BcoefH(ikl)
+            AcoQSV(ikl) = AcoefQ(ikl)
+            BcoQSV(ikl) = BcoefQ(ikl)
+            cdH_SV(ikl) = min(cdragh(ikl),cdmax)
+            cdM_SV(ikl) = min(cdragm(ikl),cdmax)
+            rcdmSV(ikl) = sqrt(cdM_SV(ikl))
+            Us_min = 0.01
+            us__SV(ikl) = max(Us_min, ustar(ikl))
+            ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6))
+            rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6))
+
+            ! +--Energy Fluxes                                          (INPUT/OUTPUT)
+            ! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^
+            !IF (.not.firstcall) THEN
+            Tsrfsv(ikl)  = tsurf(ikl)                     !hj 12 03 2010
+            cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness
+            !END IF
+
+         END DO
+
+        !
+        ! +  +++  READ FORCINGS:  END  +++
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 
+
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! +--SISVAT EXECUTION
+        ! +  ----------------
+
+        call  INLANDSIS(SnoMod, BloMod, 1)
+
+
+       
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! + RETURN RESULTS
+        ! + --------------
+        ! + Return (compressed) SISVAT variables to LMDZ
+        ! +
+        DO  ikl = 1, knon                  ! use only 1:knon (actual ice sheet..)
+            dflux_s(ikl) = dSdTSV(ikl)         ! Sens.H.Flux T-Der.
+            dflux_l(ikl) = dLdTSV(ikl)         ! Latn.H.Flux T-Der.
+            fluxsens(ikl) = HSs_sv(ikl)         ! HS
+            fluxlat(ikl) = HLs_sv(ikl)         ! HL
+            evap(ikl) = -1*HLs_sv(ikl) / LHvH2O  ! Evaporation
+            erod(ikl) = 0.
+
+            IF (BloMod) THEN
+                ! + Blowing snow
+
+                !       SLussl(i,j,n)= 0.
+                ! #BS   SLussl(i,j,n)=                     !Effective erosion
+                ! #BS. (- dbs_ER(ikl))/(dt*rhT_SV(ikl))    !~u*qs* from previous time step
+                ! #BS   blowSN(i,j,n)=  dt*uss_SV(ikl)     !New max. pot. Erosion [kg/m2]
+                ! #BS.                    *rhT_SV(ikl)     !(further bounded in sisvat_bsn.f)
+                ! #BS  erprev(i,j,n) =     dbs_Er(ikl)/dt__SV
+                erod(ikl) = dbs_Er(ikl) / dt__SV
+            ENDIF
+
+            ! +   Check snow thickness,  substract if too thick, add if too thin
+
+            sissnow(ikl) = 0.  !()
+            DO  isn = 1, isnoSV(ikl)
+                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
+            END DO
+
+            IF (sissnow(ikl) .LE. sn_low) THEN  !add snow
+                IF (isnoSV(ikl).GE.1) THEN
+                    dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi)
+                    toicSV(ikl) = toicSV(ikl) - sn_add
+                ELSE
+                    write(*, *) 'Attention, bare ice... point ', ikl
+                    isnoSV(ikl) = 1
+                    istoSV(ikl, 1) = 0
+                    ro__SV(ikl, 1) = 350.
+                    dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi)  ! 1.
+                    eta_SV(ikl, 1) = epsi
+                    TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2)
+                    G1snSV(ikl, 1) = 0.
+                    G2snSV(ikl, 1) = 0.3
+                    agsnSV(ikl, 1) = 10.
+                    toicSV(ikl) = toicSV(ikl) - sn_add
+                END IF
+            END IF
+
+            IF (sissnow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
+                dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div
+                toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div
+            END IF
+
+            sissnow(ikl) = 0.
+            qsnow(ikl) = 0.
+            snow(ikl) = 0.
+            snowhgt(ikl) = 0.
+
+            DO  isn = 1, isnoSV(ikl)
+                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
+                snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn)
+                ! Etienne: check calc qsnow
+                qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn)
+            END DO
+
+            zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0)
+            ! Etienne: comment following line
+            ! snow(ikl)    = sissnow(ikl)+toicSV(ikl)
+            snow(ikl) = sissnow(ikl)
+
+            to_ice(ikl) = toicSV(ikl)
+            runoff_lic(ikl) = RnofSV(ikl)    ! RunOFF: intensity (flux due to melting + liquid precip)
+            fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing
+            ffonte(ikl)=fqfonte(ikl)*Lf_H2O
+
+            qsol(ikl) = 0.
+            DO  isl = -nsol, 0
+                tsoil(ikl, 1 - isl) = TsisSV(ikl, isl)       ! Soil Temperature
+                ! Etienne: check calc qsol
+                qsol(ikl) = qsol(ikl)                      &
+                        + eta_SV(ikl, isl) * dz_dSV(isl)
+            END DO
+            agesno(ikl) = agsnSV(ikl, isnoSV(ikl))        !          [day]
+
+            alb1(ikl) = alb1sv(ikl)             ! Albedo VIS
+!            alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl)                   &
+!                    & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1
+            alb2(ikl)=alb2sv(ikl)
+            ! Albedo NIR
+            alb3(ikl) = alb3sv(ikl)             ! Albedo FIR
+            ! 6 band Albedo
+            alb6(ikl,:)=alb6sv(ikl,:)
+
+            tsurf_new(ikl) = Tsrfsv(ikl)
+
+            qsurf(ikl) = QsT_SV(ikl)
+            emis_new(ikl) = eps0SL(ikl)
+            z0m(ikl) = Z0m_SV(ikl)
+            z0h(ikl) = Z0h_SV(ikl)
+
 
         END DO
 
-
-
-
-        DO ikl=1,knon     
-
-
-! Initialise variables
-          
-          ispiSV(ikl)             = 0 
-          iiceSV(ikl)             = 0  
-          rusnSV(ikl)             = 0.   
-          toicSV(ikl)             = 0.     
-          isnoSV(ikl)             = 0.       ! # snow layers                            
-          istoSV(ikl,:)           = 0. 
-          eta_SV(ikl,:)           = 0.     
-          TsisSV(ikl,:)           = 0.
-          ro__SV(ikl,:)           = 0.        
-          G1snSV(ikl,:)           = 0.  
-          G2snSV(ikl,:)           = 0. 
-          agsnSV(ikl,:)           = 0. 
-          dzsnSV(ikl,:)           = 0. 
-          zzsnsv(ikl,:)           = 0.                                           
-          BufsSV(ikl)             = 0.    
-          qsnoSV(ikl)             = 0.     ! BL snow content  
-          zWEcSV(ikl)             = 0.
-          dbs_SV(ikl)             = 0.
-          dsnbSV(ikl)             = 0.
-          esnbSV(ikl)             = 0.
-          BrosSV(ikl)             = 0.
-          BG1sSV(ikl)             = 0.          
-          BG2sSV(ikl)             = 0.
-          SWS_SV(ikl)             = 0.
-          RnofSV(ikl)             = 0.     ! RunOFF Intensity 
-          RRs_SV(ikl)             = 0.
-          DDs_SV(ikl)             = 0.
-          VVs_SV(ikl)             = 0.
-          cld_SV(ikl)             = 0.      
-          uts_SV(ikl)             = 0.     ! u*T*  arbitrary  
-          uqs_SV(ikl)             = 0.     ! u*q*    " 
-          uss_SV(ikl)             = 0.     ! u*s*    "
-          LMO_SV(ikl)             = 0.
-
-
-! Set variables 
-
-          LSmask(ikl) = 1                          ! Land/Sea   Mask    
-          isotSV(ikl) = 12                         ! Soil       Type  -> 12= ice  
-          iWaFSV(ikl) = 1                          ! Soil Drainage                                     
-          eps0SL(ikl )= 1.
-          alb0SV(ikl) = alb_soil(ikl)                 ! Soil Albedo        
-          Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
-          Z0h_SV(ikl) = z0h(ikl)                      ! heat Roughn.L.
-
-! + Soil Upward IR Flux, Water Fluxes, roughness length  
-          IRs_SV(ikl) =                               &
-              -eps0SL(ikl)* StefBo*(temp_air(ikl)**4)   ! Upward IR Flux    
-          Tsf_SV(ikl) = min(temp_air(ikl),TfSnow)       
-          
-! + Soil
-        DO isl =   -nsol,0    
-          TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2)   !temp_air(ikl)  !tsoil(ikl,1-isl)   Soil Temperature
-          !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
-          eta_SV(ikl,isl) = epsi                        !etasoil(ikl,1-isl) Soil Water[m3/m3]
-          ro__SV(ikl,isl) = rhoIce                         !rosoil(ikl,1-isl)  volumic mass
-        END DO      
-
-
-
-!! Initialise with snow
-          !  G1snSV(ikl,0)          = 0.                   !      [-]     
-          !  G2snSV(ikl,0)          = 1.6                  ! [-] [0.0001 m] 
-          !  dzsnSV(ikl,0)          = dz_dSV(0)            !           [m] 
-
-
-          ! if (snow(ikl) .GT. 0.) then
-          !   isnoSV(ikl)             = 1       ! snow layers                            
-          !   istoSV(ikl,1:nsno)      = 0     ! 0,...,5 :   Snow     History (see istdSV data)
-          !   eta_SV(ikl,1:nsno)      = epsi     
-          !   TsisSV(ikl,1:nsno)      = tsoil(ikl,1)         
-          !   ro__SV(ikl,1:nsno)      = 350.0       
-          !   G1snSV(ikl,1:nsno)      = 0.   !      [-]  
-          !   G2snSV(ikl,1:nsno)      = 1.6   !     [-] [0.0001 m] 
-          !   agsnSV(ikl,1:nsno)      = 50.   !          [day] 
-          !   dzsnSV(ikl,1)           = snow(ikl)/max(ro__SV(ikl,1),epsi) ![m] 
-! ! ecrete si trop de neige:
-!              IF (snow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
-!               dzsnSV(ikl,1)      = dzsnSV(ikl,1)/sn_div
-!               toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div
-!              END IF
-!            zzsnsv(ikl,1)      =  dzsnSV(ikl,1)                    ! Total snow pack thickness
-!          endif
-
-
-! Initialise la neige avec un profil de densité prochde des conditions de Dôme C (~10m de neige avec 19 niveaux) (Etienne):
-          isnoSV(ikl)                    = 19
-          istoSV(ikl,1:isnoSV(ikl))      = 100
-          ro__SV(ikl,1:isnoSV(ikl))      = 350.      
-          eta_SV(ikl,1:isnoSV(ikl))      = epsi
-          TsisSV(ikl,1:isnoSV(ikl))      = min(tsoil(ikl,1),TfSnow-0.2)
-          G1snSV(ikl,1:isnoSV(ikl))      = 0
-          G2snSV(ikl,1:isnoSV(ikl))      = 1.6 
-          agsnSV(ikl,1:isnoSV(ikl))      = 50.    
-          dzsnSV(ikl,19)                  = 0.015
-          dzsnSV(ikl,18)                  =0.015
-          dzsnSV(ikl,17)                  =0.020
-          dzsnSV(ikl,16)                  =0.030
-          dzsnSV(ikl,15)                  =0.040
-          dzsnSV(ikl,14)                  =0.060
-          dzsnSV(ikl,13)                  =0.080
-          dzsnSV(ikl,12)                  =0.110
-          dzsnSV(ikl,11)                  =0.150
-          dzsnSV(ikl,10)                  =0.200
-          dzsnSV(ikl,9)                   =0.300
-          dzsnSV(ikl,8)                   =0.420
-          dzsnSV(ikl,7)                   =0.780
-          dzsnSV(ikl,6)                   =1.020
-          dzsnSV(ikl,5)                   =0.980
-          dzsnSV(ikl,4)                   =1.020
-          dzsnSV(ikl,3)                   =3.970
-          dzsnSV(ikl,2)                   =1.020
-          dzsnSV(ikl,1)                   =0.100
-
-
-        END DO                                            
-
-! +--Surface Fall Line Slope                                                   
-! +  -----------------------  
-        IF (SnoMod)  THEN                
-          DO ikl=1,knon  
-            slopSV(ikl) = slope(ikl)
-            SWf_SV(ikl) =             &   ! Normalized Decay of the  
-              exp(-dt__SV             &   ! Surficial Water Content  
-              /(c1_zuo                &   !(Zuo and Oerlemans 1996,  
-            +c2_zuo*exp(-c3_zuo*abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
-          END DO                                     
-        END IF                            
-
-! + SISVAT_ini (as for use with MAR, but not computing soil layers)
-! + -------------------------------------------------------------
-!        write(*,'(/a)') 'Start SISVAT initialization: SISVAT_ini' 
-        CALL SISVAT_ini(knon)
-
-
-! +--Read restart file
-! +  =================================================   
-        
-        INQUIRE(FILE="startsis.nc", EXIST=file_exists)
-        IF (file_exists) THEN
-        CALL sisvatetat0("startsis.nc",ikl2i)
+            IF (ok_outfor) THEN
+             ikl= gp_outfor
+            write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++'
+            write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl)
+            write(un_outfor, *) dzsnSV(ikl, :)
+            write(un_outfor, *) TsisSV(ikl, :)
+            write(un_outfor, *) ro__SV(ikl, :)
+            write(un_outfor, *) eta_SV(ikl, :)
+            write(un_outfor, *) G1snSV(ikl, :)
+            write(un_outfor, *) G2snSV(ikl, :)
+            write(un_outfor, *) agsnSV(ikl, :)
+            write(un_outfor, *) istoSV(ikl, :)
+            write(un_outfor, *) DOPsnSV(ikl, :)
+        ENDIF
+
+
+
+        ! +  -----------------------------
+        ! +  END --- RETURN RESULTS
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        IF (lafin) THEN
+            fichnom = "restartsis.nc"
+            CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat)
+
+            IF (ok_outfor) THEN
+                close(unit = un_outfor)
+            END IF
         END IF
-        
-        
-        
-! +--Output ascii file
-! +  =================================================   
-               
-        
-        
-        ! open output file 
-        IF (ok_outfor) THEN
-          open(unit=un_outfor,status='replace',file=fn_outfor)         
-          ikl=gp_outfor     ! index sur la grille land ice
-          write(un_outfor,*) fn_outfor, ikl, dt__SV   
-          write(un_outfor,*) 'nsnow - albedo - z0m - z0h , dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46] &
-           & G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]'
-
-        END IF
- 
-      END IF  ! firstcall                        
-! +                                
-! +  +++  INITIALISATION:  END  +++                               
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! + READ FORCINGS  
-! + ------------------------ 
-
-! + Update Forcings for SISVAT given by the LMDZ model.
-! +
-      DO ikl=1,knon
-
-! +--Atmospheric Forcing                                    (INPUT)            
-! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^ 
-        zSBLSV      = 1000.                         ! [m]                
-        za__SV(ikl) = zsl_height(ikl)               ! surface layer height (fisr model level) [m]
-        Ua_min      = epsi                          !
-        Ua_min      = 0.2 * sqrt(za__SV(ikl)   )    !                    
-        VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
-        TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
-        ExnrSV(ikl) = pexner(ikl)                   ! Exner potential          
-        rhT_SV(ikl) = dens_air(ikl)                 ! Air density           
-        QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity 
-        ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
-        p1l_SV(ikl) = p1lay(ikl)                    ! lowest atm. layer press[Pa]
-
-! +--Surface properties
-! +  ^^^^^^^^^^^^^^^^^^
-
-        Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
-        Z0h_SV(ikl) = z0h(ikl)                      ! Moment.Roughn.L.
-
-! +--Energy Fluxes                                          (INPUT)           
-! +  ^^^^^^^^^^^^^                                           ^^^^^             
-        coszSV(ikl) = max(czemin,rmu0(ikl))         ! cos(zenith.Dist.)  
-        sol_SV(ikl) = swdown(ikl)                   ! downward Solar  
-        IRd_SV(ikl) = lwdown(ikl)                   ! downward IR    
-        rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.   
-
-! +--Water  Fluxes                                          (INPUT)           
-! +  ^^^^^^^^^^^^^                                           ^^^^^             
-        drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
-        dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]
-!c #BS  dbsnow      = -SLussl(i,j,n)                ! Erosion    
-!c #BS.               *dtPhys     *rhT_SV(ikl) /ro_Wat                    
-!c #BS  dsnbSV(ikl) = snow_adv(ikl)  ! min(max(zero,dbsnow)              
-!c #BS.                    /    max(epsi,d_snow),unun)                    
-!c #BS  dbs_SV(ikl) = snow_cont_air(ikl)
-!c #BS                  blowSN(i,j,n)               !          [kg/m2]  
-                                                                              
-! +--Soil/BL                                      (INPUT)            
-! +  ^^^^^^^                                       ^^^^^            
-        alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo 
-        AcoHSV(ikl) = AcoefH(ikl)  
-        BcoHSV(ikl) = BcoefH(ikl)                     
-        AcoQSV(ikl) = AcoefQ(ikl)  
-        BcoQSV(ikl) = BcoefQ(ikl)             
-        cdH_SV(ikl) = cdragh(ikl)     
-        cdM_SV(ikl) = cdragm(ikl)     
-        Us_min      = 0.01
-        us__SV(ikl) = max(Us_min, ustar(ikl)) 
-        ram_sv(ikl) = 1./(cdragm(ikl)*max(VV__SV(ikl),eps6))
-        rah_sv(ikl) = 1./(cdragh(ikl)*max(VV__SV(ikl),eps6))   
-
-! +--Energy Fluxes                                          (INPUT/OUTPUT)    
-! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^    
-        IF (.not.firstcall) THEN  
-        Tsf_SV(ikl) = tsurf(ikl)                     !hj 12 03 2010
-        cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness         
-        END IF
- 
-
-      END DO
-
-!                            
-! +  +++  READ FORCINGS:  END  +++    
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! +--SISVAT EXECUTION                                                          
-! +  ----------------                                                          
-
-      call  INLANDSIS(SnoMod,BloMod,1) 
-
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! + RETURN RESULTS  
-! + -------------- 
-! + Return (compressed) SISVAT variables to LMDZ             
-! + 
-      DO  ikl=1,knon                  ! use only 1:knon (actual ice sheet..) 
-        runoff_lic(ikl)    = RnofSV(ikl)*dtime   ! RunOFF: intensity* time step
-        dflux_s(ikl)       = dSdTSV(ikl)         ! Sens.H.Flux T-Der. 
-        dflux_l(ikl)       = dLdTSV(ikl)         ! Latn.H.Flux T-Der. 
-        fluxsens(ikl)      = HSs_sv(ikl)         ! HS                 
-        fluxlat(ikl)       = HLs_sv(ikl)         ! HL                 
-        evap(ikl)          = HLs_sv(ikl)/LHvH2O  ! Evaporation  
-        snow(ikl)          = 0.
-        snowhgt(ikl)       = 0.
-        qsnow(ikl)         = 0.
-        qsol(ikl)          = 0. 
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! +
-! +   Check snow thickness,  substract if too thick   (commended by etienne: add if too thin) 
-
-        sissnow(ikl)       = 0.  !()
-      DO  isn = 1,isnoSV(ikl)                                               
-        sissnow(ikl)       = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)   
-      END DO
-
-       IF (sissnow(ikl) .LE. sn_low) THEN  !add snow 
-       IF (isnoSV(ikl).GE.1) THEN
-         dzsnSV(ikl,1)      = dzsnSV(ikl,1) + sn_add/max(ro__SV(ikl,1),epsi)  
-         toicSV(ikl)        = toicSV(ikl)   - sn_add
-!       ELSE
-!         write(*,*) 'Attention, bare ice... point ',ikl
-!         isnoSV(ikl)        = 1
-!         istoSV(ikl,1)      = 0
-!         ro__SV(ikl,1)      = 350.      
-!         dzsnSV(ikl,1)      = sn_add/max(ro__SV(ikl,1),epsi)  ! 1.
-!         eta_SV(ikl,1)      = epsi
-!         TsisSV(ikl,1)      = min(TsisSV(ikl,0),TfSnow-0.2)   
-!         G1snSV(ikl,1)      = 0.
-!         G2snSV(ikl,1)      = 0.3 
-!         agsnSV(ikl,1)      = 10.    
-!         toicSV(ikl)        = toicSV(ikl)   - sn_add
-       END IF
-       END IF
-
-      IF (sissnow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
-        dzsnSV(ikl,1)      = dzsnSV(ikl,1)/sn_div
-        toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div
-      END IF
-
-        sissnow(ikl)       = 0.  !()
-
-        DO  isn = 1,isnoSV(ikl)                                               
-        sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)           
-        snowhgt(ikl) = snowhgt(ikl)+dzsnSV(ikl,isn)              
-        qsnow(ikl)   = qsnow(ikl)+1e03*eta_SV(ikl,isn)*dzsnSV(ikl,isn)   
-        END DO
-
-        ! Etienne: pourquoi ajouter toicSV ici? Pour bilan d'eau?
-        snow(ikl)    = sissnow(ikl)+toicSV(ikl)
-        to_ice(ikl)  = toicSV(ikl)
-
-
-      DO  isl =   -nsol,0    
-        tsoil(ikl,1-isl)   = TsisSV(ikl,isl)       ! Soil Temperature  
-        qsol(ikl)          = qsol(ikl)                      &    
-                            +eta_SV(ikl,isl) * dz_dSV(isl)  
-      END DO                                               
-        agesno(ikl)        = agsnSV(ikl,isnoSV(ikl))        !          [day] 
-
-        alb1(ikl)          = alb1sv(ikl)             ! Albedo VIS  
-        alb2(ikl)          = ((So1dSV-f1)*alb1sv(ikl)                   &
-     &                       +So2dSV*alb2sv(ikl)+So3dSV*alb3sv(ikl))/f1    
-                                                     ! Albedo NIR 
-        alb3(ikl)          = alb3sv(ikl)             ! Albedo FIR  
-
-        tsurf_new(ikl)     =Tsrfsv(ikl)
-
-        zfra(ikl)          = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
-        qsurf(ikl)         = QaT_SV(ikl) 
-        emis_new(ikl)      = eps0SL(ikl)  
-        z0m(ikl)           = Z0m_SV(ikl) 
-        z0h(ikl)           = Z0h_SV(ikl) 
-
-      END DO  ! ikl
-      
-
-
-
-
-
-! write variables in output file
-
-      IF (ok_outfor) THEN
-        ikl=gp_outfor
-
-!        write(un_outfor,*) 'nsnow [-,1], dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46]'
-!        write(un_outfor,*) 'G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]'
-        write(un_outfor,*) '+++++++++++++++++++++++++++++++++++++++++++++++'
-        write(un_outfor,*) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl)
-        write(un_outfor,*) dzsnSV(ikl,:)
-        write(un_outfor,*) TsisSV(ikl,:)
-        write(un_outfor,*) ro__SV(ikl,:)
-        write(un_outfor,*) eta_SV(ikl,:)
-        write(un_outfor,*) G1snSV(ikl,:)
-        write(un_outfor,*) G2snSV(ikl,:)
-        write(un_outfor,*) agsnSV(ikl,:)
-        write(un_outfor,*) istoSV(ikl,:)
-       
-      ENDIF
-
-
-
-
-! +  -----------------------------                              
-! +  END --- RETURN RESULTS    
-! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-      IF (lafin) THEN
-        fichnom = "restartsis.nc"
-        CALL sisvatredem("restartsis.nc",ikl2i,rlon,rlat)     
-        
-        IF (ok_outfor) THEN
-          close(unit=un_outfor)                    
-        END IF 
-      END IF                                                          
-
-
-  END SUBROUTINE surf_inlandsis
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-!=======================================================================
-
-  SUBROUTINE get_soil_levels(dz1, dz2, lambda)
-! ======================================================================
-! Routine to compute the vertical discretization of the soil in analogy
-! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case 
-! of SISVAT, therefore it's needed here.
-!
-    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
-    USE mod_phys_lmdz_para
-    USE VAR_SV 
-
-
-!    INCLUDE "dimsoil.h"
-
-    REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
-    REAL, INTENT(OUT)                     :: lambda
-
-
-!-----------------------------------------------------------------------
-!   Depthts:
-!   --------
-    REAL fz,rk,fz1,rk1,rk2
-    REAL min_period, dalph_soil
-    INTEGER ierr,jk
-
-    fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
-
-!    write(*,*)'Start soil level computation' 
-!-----------------------------------------------------------------------
-! Calculation of some constants
-! NB! These constants do not depend on the sub-surfaces
-!-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-!   ground levels 
-!   grnd=z/l where l is the skin depth of the diurnal cycle:
-!-----------------------------------------------------------------------
-
-     min_period=1800. ! en secondes
-     dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
-! !$OMP MASTER
-!     IF (is_mpi_root) THEN
-!        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
-!        IF (ierr == 0) THEN ! Read file only if it exists
-!           READ(99,*) min_period
-!           READ(99,*) dalph_soil
-!           PRINT*,'Discretization for the soil model'
-!           PRINT*,'First level e-folding depth',min_period, &
-!                '   dalph',dalph_soil
-!           CLOSE(99)
-!        END IF
-!     ENDIF
-! !$OMP END MASTER
-!     CALL bcast(min_period)
-!     CALL bcast(dalph_soil)
-
-!   la premiere couche represente un dixieme de cycle diurne
-     fz1=SQRT(min_period/3.14)
-     
-     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)
-     PRINT*,'full layers, intermediate layers (seconds)'
-     DO jk=1,nsoilmx
-        rk=jk
-        rk1=jk+.5
-        rk2=jk-.5
-        PRINT *,'fz=', &
-             fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
-     ENDDO
-
-  END SUBROUTINE get_soil_levels
-  
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-!===========================================================================
-
-  SUBROUTINE SISVAT_ini(knon)                                                     
-                                                                               
-!C +------------------------------------------------------------------------+  
-!C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR | 
-!C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters |  
-!C +------------------------------------------------------------------------+  
-!C |   PARAMETERS:  klonv: Total Number of columns =                        |  
-!C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |  
-!C |                     X       Number of Mosaic Cell per grid box         |  
-!C |                                                                        |  
-!C |   INPUT:   dt__SV   : Time  Step                                   [s] | 
-!C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] |  
-!C |                                                                        |  
-!C |   OUTPUT:             [-] |  
-!C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] |  
-!C |            etamSV   : Soil Minimum Humidity                    [m3/m3] |  
-!C |                      (based on a prescribed Soil Relative Humidity)    |  
-!C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       |  
-!C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        |  
-!C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   |  
-!C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] |  
-!C |            dzsnSV(0): Soil first Layer Thickness                   [m] |  
-!C |            dzmiSV   : Distance between two contiguous levels       [m] |  
-!C |            dz78SV   : 7/8 (Layer Thickness)                        [m] | 
-!C |            dz34SV   : 3/4 (Layer Thickness)                        [m] | 
-!C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] | 
-!C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] | 
-!C |            dtz_SV   : dt/dz                                      [s/m] | 
-!C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
-!C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      |  
-!C |            Explic   : Explicit Parameter = 1.0 - Implic                |  
-!C |                                                                        | 
-!C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
-!C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
-!C | #          #SH: Hapex-Sahel Values                                     !  
-!C |                                                                        | 
-!C +------------------------------------------------------------------------+  
-!                                                                              
-!                                                                             
-                                                                              
-!C +--Global Variables                                                         
-!C +  ================         
-
-      USE dimphy
-      USE VARphy                                           
-      USE VAR_SV                                                     
-      USE VARdSV                                                         
-      USE VAR0SV                                                           
-      USE VARxSV
-      USE VARtSV
-      USE VARxSV
-      USE VARySV
-      IMPLICIT NONE                                                           
-                                                                              
-                                                                              
-                                                                              
-!C +--Arguments                                                     
-!C +  ==================                                                      
-       INTEGER,INTENT(IN) ::  knon                                       
-
-!C +--Internal Variables                                                     
-!C +  ==================                                                      
-                                                                               
-      INTEGER ::  ivt   ,ist   ,ikl   ,isl   ,isn   ,ikh               
-      INTEGER ::  misl_2,nisl_2                                              
-      REAL    ::  d__eta,eta__1,eta__2,Khyd_1,Khyd_2                           
-      REAL,PARAMETER  ::  RHsMin=  0.001        ! Min.Soil Relative Humidity   
-      REAL    ::  PsiMax                        ! Max.Soil Water    Potential 
-      REAL    ::  a_Khyd,b_Khyd                 ! Piecewis.https://www.lequipe.fr/Water Conductivity 
-
-
-!c #WR REAL    ::  Khyd_x,Khyd_y                                               
-                                                                              
-                                                                              
-                                                                 
-!C +--Non Time Dependant SISVAT parameters                                    
-!C +  ====================================                                
-                                                                              
-!C +--Soil Discretization                                                     
-!C +  -------------------                                                     
-                                                                              
-!C +--Numerical Scheme Parameters                                             
-!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^                                             
+
+    END SUBROUTINE surf_inlandsis
+
+
+    !=======================================================================
+
+    SUBROUTINE get_soil_levels(dz1, dz2, lambda)
+        ! ======================================================================
+        ! Routine to compute the vertical discretization of the soil in analogy
+        ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
+        ! of SISVAT, therefore it's needed here.
+        !
+        USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
+        USE mod_phys_lmdz_para
+        USE VAR_SV
+
+
+        !    INCLUDE "dimsoil.h"
+
+        REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
+        REAL, INTENT(OUT) :: lambda
+
+
+        !-----------------------------------------------------------------------
+        !   Depthts:
+        !   --------
+        REAL fz, rk, fz1, rk1, rk2
+        REAL min_period, dalph_soil
+        INTEGER ierr, jk
+
+        fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
+
+        !    write(*,*)'Start soil level computation'
+        !-----------------------------------------------------------------------
+        ! Calculation of some constants
+        ! NB! These constants do not depend on the sub-surfaces
+        !-----------------------------------------------------------------------
+        !-----------------------------------------------------------------------
+        !   ground levels
+        !   grnd=z/l where l is the skin depth of the diurnal cycle:
+        !-----------------------------------------------------------------------
+
+        min_period = 1800. ! en secondes
+        dalph_soil = 2.    ! rapport entre les epaisseurs de 2 couches succ.
+        ! !$OMP MASTER
+        !     IF (is_mpi_root) THEN
+        !        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
+        !        IF (ierr == 0) THEN ! Read file only if it exists
+        !           READ(99,*) min_period
+        !           READ(99,*) dalph_soil
+        !           PRINT*,'Discretization for the soil model'
+        !           PRINT*,'First level e-folding depth',min_period, &
+        !                '   dalph',dalph_soil
+        !           CLOSE(99)
+        !        END IF
+        !     ENDIF
+        ! !$OMP END MASTER
+        !     CALL bcast(min_period)
+        !     CALL bcast(dalph_soil)
+
+        !   la premiere couche represente un dixieme de cycle diurne
+        fz1 = SQRT(min_period / 3.14)
+
+        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)
+        DO jk = 1, nsoilmx
+            rk = jk
+            rk1 = jk + .5
+            rk2 = jk - .5
+        ENDDO
+
+    END SUBROUTINE get_soil_levels
+
+
+    !===========================================================================
+
+    SUBROUTINE SISVAT_ini(knon)
+
+        !C +------------------------------------------------------------------------+
+        !C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR |
+        !C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters |
+        !C +------------------------------------------------------------------------+
+        !C |   PARAMETERS:  klonv: Total Number of columns =                        |
+        !C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
+        !C |                     X       Number of Mosaic Cell per grid box         |
+        !C |                                                                        |
+        !C |   INPUT:   dt__SV   : Time  Step                                   [s] |
+        !C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] |
+        !C |                                                                        |
+        !C |   OUTPUT:             [-] |
+        !C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] |
+        !C |            etamSV   : Soil Minimum Humidity                    [m3/m3] |
+        !C |                      (based on a prescribed Soil Relative Humidity)    |
+        !C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       |
+        !C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        |
+        !C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   |
+        !C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] |
+        !C |            dzsnSV(0): Soil first Layer Thickness                   [m] |
+        !C |            dzmiSV   : Distance between two contiguous levels       [m] |
+        !C |            dz78SV   : 7/8 (Layer Thickness)                        [m] |
+        !C |            dz34SV   : 3/4 (Layer Thickness)                        [m] |
+        !C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] |
+        !C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] |
+        !C |            dtz_SV   : dt/dz                                      [s/m] |
+        !C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
+        !C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      |
+        !C |            Explic   : Explicit Parameter = 1.0 - Implic                |
+        !C |                                                                        |
+        !C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
+        !C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
+        !C | #          #SH: Hapex-Sahel Values                                     !
+        !C |                                                                        |
+        !C +------------------------------------------------------------------------+
+        !
+        !
+
+        !C +--Global Variables
+        !C +  ================
+
+        USE dimphy
+        USE VARphy
+        USE VAR_SV
+        USE VARdSV
+        USE VAR0SV
+        USE VARxSV
+        USE VARtSV
+        USE VARxSV
+        USE VARySV
+        IMPLICIT NONE
+
+
+
+        !C +--Arguments
+        !C +  ==================
+        INTEGER, INTENT(IN) :: knon
+
+        !C +--Internal Variables
+        !C +  ==================
+
+        INTEGER :: ivt, ist, ikl, isl, isn, ikh
+        INTEGER :: misl_2, nisl_2
+        REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2
+        REAL, PARAMETER :: RHsMin = 0.001        ! Min.Soil Relative Humidity
+        REAL :: PsiMax                        ! Max.Soil Water    Potential
+        REAL :: a_Khyd, b_Khyd                 ! Water conductivity
+
+
+        !c #WR REAL    ::  Khyd_x,Khyd_y
+
+
+
+        !C +--Non Time Dependant SISVAT parameters
+        !C +  ====================================
+
+        !C +--Soil Discretization
+        !C +  -------------------
+
+        !C +--Numerical Scheme Parameters
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
         Implic = 0.75                           ! 0.5  <==> Crank-Nicholson  
         Explic = 1.00 - Implic                  !                            
-                                                                             
-!C +--Soil/Snow Layers Indices                                                
-!C +  ^^^^^^^^^^^^^^^^^^^^^^^^                                                
-      DO  isl=-nsol,0                                                          
-        islpSV(isl) =           isl+1                                         
-        islpSV(isl) = min(      islpSV(isl),0)                                
-        islmSV(isl) =           isl-1                                         
-        islmSV(isl) = max(-nsol,islmSV(isl))                                 
-      END DO                                                                  
-                                                                               
-      DO  isn=1,nsno                                                           
-        isnpSV(isn) =           isn+1                                          
-        isnpSV(isn) = min(      isnpSV(isn),nsno)                            
-      END DO                                                                 
-                                                                              
-!C +--Soil      Layers Thicknesses                                             
-!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^  
-! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!    
-!c #kd IF (nsol.gt.4)                                              THEN        
-!c #kd   DO isl=-5,-nsol,-1                                                    
-!c #kd     dz_dSV(isl)=   1.                                                  
-!c #kd   END DO                                                                
-!c #kd END IF                                                                 
-!                                                                             
-!      IF (nsol.ne.4)                                              THEN        
-!        DO isl= 0,-nsol,-1                                                   
-!          misl_2 =     -mod(isl,2)                                          
-!          nisl_2 =         -isl/2                                            
-!          dz_dSV(isl)=(((1-misl_2) * 0.001                                   
-!     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.             
-!C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm            
-!                                                                             
-!c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001                                   
-!c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.              
-!                                                                             
-!c #05     dz_dSV(isl)=(((1-misl_2) * 0.001                                    
-!c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5             
-!        END DO                                                                
-!          dz_dSV(0)  =               0.001                                    
-!          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004         
-!      END IF           
-
-                                                                               
-        zz_dSV      = 0.                                                      
-      DO  isl=-nsol,0                                                        
-        dzmiSV(isl) = 0.500*(dz_dSV(isl)        +dz_dSV(islmSV(isl)))        
-        dziiSV(isl) = 0.500* dz_dSV(isl)        /dzmiSV(isl)                  
-        dzi_SV(isl) = 0.500* dz_dSV(islmSV(isl))/dzmiSV(isl)                   
-        dtz_SV(isl) =        dt__SV             /dz_dSV(isl)  
-        dtz_SV2(isl) =        1.            /dz_dSV(isl)                          
-        dz78SV(isl) = 0.875* dz_dSV(isl)                                      
-        dz34SV(isl) = 0.750* dz_dSV(isl)                                      
-        dz_8SV(isl) = 0.125* dz_dSV(isl)                                       
-        dzAvSV(isl) = 0.125* dz_dSV(islmSV(isl))                        &
-     &              + 0.750* dz_dSV(isl)                                &
-     &              + 0.125* dz_dSV(islpSV(isl))                                                              
-        zz_dSV      = zz_dSV+dz_dSV(isl)                                      
-      END DO                                                                  
-      DO ikl=1,knon !v                                                          
-        dzsnSV(ikl,0) =      dz_dSV(0)                                        
-      END DO                                                                  
-                                                                              
-!C +--Conversion to a 50 m Swab Ocean Discretization                           
-!C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                           
-        OcndSV = 0.                                                           
-      DO isl=-nsol,0                                                          
-        OcndSV = OcndSV +dz_dSV(isl)                                           
-      END DO                                                                   
-        OcndSV = 50.    /OcndSV                                                
-                                                                                                            
-                                                                                                                                                                                                                                
-!C +--Secondary Soil       Parameters                                          
-!C +  -------------------------------                                         
-                                                                               
-      DO  ist=0,nsot                                                           
-         rocsSV(ist)=(1.0-etadSV(ist))*1.2E+6   ! Soil Contrib. to (ro c)_s    
-         s1__SV(ist)=     bCHdSV(ist)          & ! Factor of (eta)**(b+2)     
-     &  *psidSV(ist)     *Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)     
-     & /(etadSV(ist)**(   bCHdSV(ist)+3.))     !                              
-         s2__SV(ist)=     Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)     
-     & /(etadSV(ist)**(2.*bCHdSV(ist)+3.))     !    in DR97, Eqn.(3.35)       
-                                                                             
-!C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)     
-!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^    
-         Psimax = -(log(RHsMin))/7.2E-5        ! DR97, Eqn 3.15 Inversion   
-         etamSV(ist) =  etadSV(ist)                                      &
-     &         *(PsiMax/psidSV(ist))**(-min(10.,1./bCHdSV(ist)))             
-      END DO                                                                  
-         etamSV(12)  =  0.                                                   
-                                                                             
-!C +--Piecewise Hydraulic Conductivity Profiles                               
-!C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
-      DO   ist=0,nsot                                                          
- 
-                                                                             
-          d__eta          =  etadSV(ist)/nkhy                                 
-          eta__1          =  0.                                              
-          eta__2          =  d__eta                                           
-        DO ikh=0,nkhy                                                         
-          Khyd_1          =  s2__SV(ist)             & ! DR97, Eqn.(3.35)      
-     &  *(eta__1      **(2. *bCHdSV(ist)+3.))        !                       
-          Khyd_2          =  s2__SV(ist)             &!                       
-     &  *(eta__2      **(2. *bCHdSV(ist)+3.))        !                       
-                                                                             
-          a_Khyd          = (Khyd_2-Khyd_1)/d__eta   !                       
-          b_Khyd          =  Khyd_1-a_Khyd *eta__1   !                       
-!c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !                       
-!c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !                       
-          aKdtSV(ist,ikh) =  a_Khyd       * dt__SV   !                       
-          bKdtSV(ist,ikh) =  b_Khyd       * dt__SV   !                         
-               
-          eta__1          = eta__1  + d__eta                              
-          eta__2          = eta__2  + d__eta                              
-        END DO                                                             
-      END DO                                                               
-
-                                                                           
-      return                                                               
-
-  END SUBROUTINE SISVAT_ini
-
-
-
-
-
-
-
-!***************************************************************************
-
-    SUBROUTINE sisvatetat0 (fichnom,ikl2i)
-
-    USE dimphy
-    USE mod_grid_phy_lmdz
-    USE mod_phys_lmdz_para
-
-    USE iostart
-    USE VAR_SV
-    USE VARdSV
-    USE VARxSV          
-    USE VARtSV
-    USE indice_sol_mod
-
-      IMPLICIT none
-!======================================================================
-! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
-! Objet: Lecture du fichier de conditions initiales pour SISVAT
-!======================================================================
-    include "netcdf.inc"
-!    include "indicesol.h"
-
-!    include "dimsoil.h"
-    include "clesphys.h"
-    include "thermcell.h"
-    include "compbl.h"
-
-!======================================================================
-    CHARACTER(LEN=*) :: fichnom
-
-
-    INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
-    REAL, DIMENSION(klon) :: rlon
-    REAL, DIMENSION(klon) :: rlat
-
-! les variables globales ecrites dans le fichier restart
-    REAL, DIMENSION(klon) :: isno
-    REAL, DIMENSION(klon) :: ispi
-    REAL, DIMENSION(klon) :: iice
-    REAL, DIMENSION(klon) :: rusn 
-    REAL, DIMENSION(klon, nsno) :: isto
-
-    REAL, DIMENSION(klon, nsismx) :: Tsis
-    REAL, DIMENSION(klon, nsismx) :: eta
-    REAL, DIMENSION(klon, nsismx) :: ro 
-
-    REAL, DIMENSION(klon, nsno) :: dzsn      
-    REAL, DIMENSION(klon, nsno) :: G1sn
-    REAL, DIMENSION(klon, nsno) :: G2sn
-    REAL, DIMENSION(klon, nsno) :: agsn
-
-    REAL, DIMENSION(klon) :: toic
-
-
-    INTEGER  :: isl, ikl, i, isn , errT, erreta, errro, errdz, snopts
-    CHARACTER (len=2) :: str2
-    LOGICAL :: found
-  
-    errT=0
-    errro=0
-    erreta=0
-    errdz=0
-    snopts=0
-! Ouvrir le fichier contenant l'etat initial:
-
-      CALL open_startphy(fichnom)
-
-! Lecture des latitudes, longitudes (coordonnees):
-
-      CALL get_field("latitude",rlat,found)
-      CALL get_field("longitude",rlon,found)
-
-      CALL get_field("n_snows", isno,found)
-      IF (.NOT. found) THEN
-        PRINT*, 'phyetat0: Le champ <n_snows> est absent'
-        PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
-      ENDIF
-
-      CALL get_field("n_ice_top",ispi,found)
-      CALL get_field("n_ice",iice,found)
-      CALL get_field("surf_water",rusn,found)
-!      IF (.NOT. found) THEN
-!        PRINT*, 'phyetat0: Le champ <surf_water> est absent'
-!        rusn(:)=0.  
-!      ENDIF
-
-
-      CALL get_field("to_ice",toic,found)
-      IF (.NOT. found) THEN
-        PRINT*, 'phyetat0: Le champ <to_ice> est absent'
-        toic(:)=0.  
-      ENDIF
-
-
-
-      DO isn = 1,nsno
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("AGESNOW"//str2, &
-                          agsn(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
+
+        !C +--Soil/Snow Layers Indices
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^
+        DO  isl = -nsol, 0
+            islpSV(isl) = isl + 1
+            islpSV(isl) = min(islpSV(isl), 0)
+            islmSV(isl) = isl - 1
+            islmSV(isl) = max(-nsol, islmSV(isl))
+        END DO
+
+        DO  isn = 1, nsno
+            isnpSV(isn) = isn + 1
+            isnpSV(isn) = min(isnpSV(isn), nsno)
+        END DO
+
+        !C +--Soil      Layers Thicknesses
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+        ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
+        !c #kd IF (nsol.gt.4)                                              THEN
+        !c #kd   DO isl=-5,-nsol,-1
+        !c #kd     dz_dSV(isl)=   1.
+        !c #kd   END DO
+        !c #kd END IF
+        !
+        !      IF (nsol.ne.4)                                              THEN
+        !        DO isl= 0,-nsol,-1
+        !          misl_2 =     -mod(isl,2)
+        !          nisl_2 =         -isl/2
+        !          dz_dSV(isl)=(((1-misl_2) * 0.001
+        !     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.
+        !C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm
+        !
+        !c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001
+        !c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.
+        !
+        !c #05     dz_dSV(isl)=(((1-misl_2) * 0.001
+        !c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5
+        !        END DO
+        !          dz_dSV(0)  =               0.001
+        !          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004
+        !      END IF
+
+        zz_dSV = 0.
+        DO  isl = -nsol, 0
+            dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl)))
+            dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl)
+            dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl)
+            dtz_SV(isl) = dt__SV / dz_dSV(isl)
+            dtz_SV2(isl) = 1. / dz_dSV(isl)
+            dz78SV(isl) = 0.875 * dz_dSV(isl)
+            dz34SV(isl) = 0.750 * dz_dSV(isl)
+            dz_8SV(isl) = 0.125 * dz_dSV(isl)
+            dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl))                        &
+                    & + 0.750 * dz_dSV(isl)                                &
+                    & + 0.125 * dz_dSV(islpSV(isl))
+            zz_dSV = zz_dSV + dz_dSV(isl)
+        END DO
+        DO ikl = 1, knon !v
+            dzsnSV(ikl, 0) = dz_dSV(0)
+        END DO
+
+        !C +--Conversion to a 50 m Swab Ocean Discretization
+        !C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        OcndSV = 0.
+        DO isl = -nsol, 0
+            OcndSV = OcndSV + dz_dSV(isl)
+        END DO
+        OcndSV = 50. / OcndSV
+
+
+        !C +--Secondary Soil       Parameters
+        !C +  -------------------------------
+
+        DO  ist = 0, nsot
+            rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6   ! Soil Contrib. to (ro c)_s
+            s1__SV(ist) = bCHdSV(ist)          & ! Factor of (eta)**(b+2)
+                    & * psidSV(ist) * Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)
+                    & / (etadSV(ist)**(bCHdSV(ist) + 3.))     !
+            s2__SV(ist) = Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)
+                    & / (etadSV(ist)**(2. * bCHdSV(ist) + 3.))     !    in DR97, Eqn.(3.35)
+
+            !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
+            !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+            Psimax = -(log(RHsMin)) / 7.2E-5        ! DR97, Eqn 3.15 Inversion
+            etamSV(ist) = etadSV(ist)                                      &
+                    & * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist)))
+        END DO
+        etamSV(12) = 0.
+
+        !C +--Piecewise Hydraulic Conductivity Profiles
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+        DO   ist = 0, nsot
+
+            d__eta = etadSV(ist) / nkhy
+            eta__1 = 0.
+            eta__2 = d__eta
+            DO ikh = 0, nkhy
+                Khyd_1 = s2__SV(ist)             & ! DR97, Eqn.(3.35)
+                        & * (eta__1      **(2. * bCHdSV(ist) + 3.))        !
+                Khyd_2 = s2__SV(ist)             &!
+                        & * (eta__2      **(2. * bCHdSV(ist) + 3.))        !
+
+                a_Khyd = (Khyd_2 - Khyd_1) / d__eta   !
+                b_Khyd = Khyd_1 - a_Khyd * eta__1   !
+                !c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !
+                !c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !
+                aKdtSV(ist, ikh) = a_Khyd * dt__SV   !
+                bKdtSV(ist, ikh) = b_Khyd * dt__SV   !
+
+                eta__1 = eta__1 + d__eta
+                eta__2 = eta__2 + d__eta
+            END DO
+        END DO
+
+        return
+
+    END SUBROUTINE SISVAT_ini
+
+
+    !***************************************************************************
+
+    SUBROUTINE sisvatetat0 (fichnom, ikl2i)
+
+        USE dimphy
+        USE mod_grid_phy_lmdz
+        USE mod_phys_lmdz_para
+
+        USE iostart
+        USE VAR_SV
+        USE VARdSV
+        USE VARxSV
+        USE VARtSV
+        USE indice_sol_mod
+
+        IMPLICIT none
+        !======================================================================
+        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
+        ! Objet: Lecture du fichier de conditions initiales pour SISVAT
+        !======================================================================
+        include "netcdf.inc"
+        !    include "indicesol.h"
+
+        !    include "dimsoil.h"
+        include "clesphys.h"
+        include "thermcell.h"
+        include "compbl.h"
+
+        !======================================================================
+        CHARACTER(LEN = *) :: fichnom
+
+        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
+        REAL, DIMENSION(klon) :: rlon
+        REAL, DIMENSION(klon) :: rlat
+
+        ! les variables globales ecrites dans le fichier restart
+        REAL, DIMENSION(klon) :: isno
+        REAL, DIMENSION(klon) :: ispi
+        REAL, DIMENSION(klon) :: iice
+        REAL, DIMENSION(klon) :: rusn
+        REAL, DIMENSION(klon, nsno) :: isto
+
+        REAL, DIMENSION(klon, nsismx) :: Tsis
+        REAL, DIMENSION(klon, nsismx) :: eta
+        REAL, DIMENSION(klon, nsismx) :: ro
+
+        REAL, DIMENSION(klon, nsno) :: dzsn
+        REAL, DIMENSION(klon, nsno) :: G1sn
+        REAL, DIMENSION(klon, nsno) :: G2sn
+        REAL, DIMENSION(klon, nsno) :: agsn
+
+        REAL, DIMENSION(klon) :: toic
+
+        INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts
+        CHARACTER (len = 2) :: str2
+        LOGICAL :: found
+
+        errT = 0
+        errro = 0
+        erreta = 0
+        errdz = 0
+        snopts = 0
+        ! Ouvrir le fichier contenant l'etat initial:
+
+        CALL open_startphy(fichnom)
+
+        ! Lecture des latitudes, longitudes (coordonnees):
+
+        CALL get_field("latitude", rlat, found)
+        CALL get_field("longitude", rlon, found)
+
+        CALL get_field("n_snows", isno, found)
+        IF (.NOT. found) THEN
+            PRINT*, 'phyetat0: Le champ <n_snows> est absent'
+            PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
         ENDIF
-      ENDDO
-      DO isn = 1,nsno
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("DZSNOW"//str2, &
-                          dzsn(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
+
+        CALL get_field("n_ice_top", ispi, found)
+        CALL get_field("n_ice", iice, found)
+        CALL get_field("surf_water", rusn, found)
+
+
+        CALL get_field("to_ice", toic, found)
+        IF (.NOT. found) THEN
+            PRINT*, 'phyetat0: Le champ <to_ice> est absent'
+            toic(:) = 0.
         ENDIF
-      ENDDO
-      DO isn = 1,nsno
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("G2SNOW"//str2, &
-                          G2sn(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-        ENDIF
-      ENDDO
-      DO isn = 1,nsno
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("G1SNOW"//str2, &
-                          G1sn(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-        ENDIF
-      ENDDO
-      DO isn = 1,nsismx
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("ETA"//str2, &
-                          eta(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-        ENDIF
-      ENDDO
-      DO isn = 1,nsismx
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("RO"//str2, &
-                          ro(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-        ENDIF
-      ENDDO
-      DO isn = 1,nsismx
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("TSS"//str2, &
-                          Tsis(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-        ENDIF
-      ENDDO
-      DO isn = 1,nsno
-        IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL get_field("HISTORY"//str2, &
-                          isto(:,isn),found)
-        ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-        ENDIF
-      ENDDO
-      write(*,*)'Read ',fichnom,' finished!!'
-
-!*********************************************************************************
-! Compress restart file variables for SISVAT
-
-
-      DO  ikl = 1,klon                                                   
-          i   = ikl2i(ikl)   
-          IF (i > 0) THEN
-              isnoSV(ikl)     = INT(isno(i))          ! Nb Snow/Ice Lay.   
-              ispiSV(ikl)     = INT(ispi(i))          ! Nb Supr.Ice Lay.  
-              iiceSV(ikl)     = INT(iice(i))          ! Nb      Ice Lay.   
-                                                                              
-
-            DO isl =   -nsol,0    
-              ro__SV(ikl,isl) = ro(i,nsno+1-isl)       !                    
-              eta_SV(ikl,isl) = eta(i,nsno+1-isl)         ! Soil Humidity      
-!hjp 15/10/2010
-              IF (eta_SV(ikl,isl) <= 1.e-6) THEN          !hj check
-                eta_SV(ikl,isl) = 1.e-6
-              ENDIF
-              TsisSV(ikl,isl) = Tsis(i,nsno+1-isl)        ! Soil Temperature  
-              IF (TsisSV(ikl,isl) <= 1.) THEN             !hj check
-!                errT=errT+1
-                TsisSV(ikl,isl) = 273.15-0.2              ! Etienne: negative temperature since soil is ice
-              ENDIF
-
-            END DO        
-            write(*,*)'Copy histo', ikl
-    
-   
-            DO  isn = 1,isnoSV(ikl) !nsno      
-              snopts=snopts+1 
-              IF (isto(i,isn) > 10.) THEN          !hj check
-                write(*,*)'Irregular isto',ikl,i,isn,isto(i,isn)
-                isto(i,isn) = 1.
-              ENDIF
-
-              istoSV(ikl,isn) = INT(isto(i,isn))     ! Snow     History 
-              ro__SV(ikl,isn) = ro(i,isn)            !        [kg/m3]     
-              eta_SV(ikl,isn) = eta(i,isn)           !        [m3/m3]  
-              TsisSV(ikl,isn) = Tsis(i,isn)          !            [K]  
-
-             IF (TsisSV(ikl,isn) <= 1.) THEN          !hj check
-              errT=errT+1
-              TsisSV(ikl,isn) = TsisSV(ikl,0)
-             ENDIF  
-             IF (TsisSV(ikl,isn) <= 1.) THEN          !hj check
-              TsisSV(ikl,isn) = 263.15
-             ENDIF
-             IF (eta_SV(ikl,isn) < 1.e-9) THEN          !hj check
-              eta_SV(ikl,isn) = 1.e-6  
-              erreta=erreta+1
-             ENDIF  
-             IF (ro__SV(ikl,isn) <= 10.) THEN          !hj check
-              ro__SV(ikl,isn) = 11.
-              errro=errro+1
-             ENDIF 
-              write(*,*)ikl,i,isn,Tsis(i,isn),G1sn(i,isn)
-              G1snSV(ikl,isn) = G1sn(i,isn)          ! [-]        [-]     
-              G2snSV(ikl,isn) = G2sn(i,isn)          ! [-] [0.0001 m] 
-              dzsnSV(ikl,isn) = dzsn(i,isn)          !            [m]         
-              agsnSV(ikl,isn) = agsn(i,isn)          !          [day]     
-            END DO  
-              rusnSV(ikl)     = rusn(i)              ! Surficial Water    
-              toicSV(ikl)     = toic(i)              ! bilan snow to ice   
-          END IF                       
-        END DO    
+
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("AGESNOW" // str2, &
+                        agsn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("DZSNOW" // str2, &
+                        dzsn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("G2SNOW" // str2, &
+                        G2sn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("G1SNOW" // str2, &
+                        G1sn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsismx
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("ETA" // str2, &
+                        eta(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsismx
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("RO" // str2, &
+                        ro(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsismx
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("TSS" // str2, &
+                        Tsis(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("HISTORY" // str2, &
+                        isto(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        write(*, *)'Read ', fichnom, ' finished!!'
+
+        !*********************************************************************************
+        ! Compress restart file variables for SISVAT
+
+        DO  ikl = 1, klon
+            i = ikl2i(ikl)
+            IF (i > 0) THEN
+                isnoSV(ikl) = INT(isno(i))          ! Nb Snow/Ice Lay.
+                ispiSV(ikl) = INT(ispi(i))          ! Nb Supr.Ice Lay.
+                iiceSV(ikl) = INT(iice(i))          ! Nb      Ice Lay.
+
+                DO isl = -nsol, 0
+                    ro__SV(ikl, isl) = ro(i, nsno + 1 - isl)       !
+                    eta_SV(ikl, isl) = eta(i, nsno + 1 - isl)         ! Soil Humidity
+                    !hjp 15/10/2010
+                    IF (eta_SV(ikl, isl) <= 1.e-6) THEN          !hj check
+                        eta_SV(ikl, isl) = 1.e-6
+                    ENDIF
+                    TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl)        ! Soil Temperature
+                    IF (TsisSV(ikl, isl) <= 1.) THEN             !hj check
+                        !                errT=errT+1
+                        TsisSV(ikl, isl) = 273.15 - 0.2              ! Etienne: negative temperature since soil is ice
+                    ENDIF
+
+                END DO
+                write(*, *)'Copy histo', ikl
+
+                DO  isn = 1, isnoSV(ikl) !nsno
+                    snopts = snopts + 1
+                    IF (isto(i, isn) > 10.) THEN          !hj check
+                        write(*, *)'Irregular isto', ikl, i, isn, isto(i, isn)
+                        isto(i, isn) = 1.
+                    ENDIF
+
+                    istoSV(ikl, isn) = INT(isto(i, isn))     ! Snow     History
+                    ro__SV(ikl, isn) = ro(i, isn)            !        [kg/m3]
+                    eta_SV(ikl, isn) = eta(i, isn)           !        [m3/m3]
+                    TsisSV(ikl, isn) = Tsis(i, isn)          !            [K]
+
+                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
+                        errT = errT + 1
+                        TsisSV(ikl, isn) = TsisSV(ikl, 0)
+                    ENDIF
+                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
+                        TsisSV(ikl, isn) = 263.15
+                    ENDIF
+                    IF (eta_SV(ikl, isn) < 1.e-9) THEN          !hj check
+                        eta_SV(ikl, isn) = 1.e-6
+                        erreta = erreta + 1
+                    ENDIF
+                    IF (ro__SV(ikl, isn) <= 10.) THEN          !hj check
+                        ro__SV(ikl, isn) = 11.
+                        errro = errro + 1
+                    ENDIF
+                    write(*, *)ikl, i, isn, Tsis(i, isn), G1sn(i, isn)
+                    G1snSV(ikl, isn) = G1sn(i, isn)          ! [-]        [-]
+                    G2snSV(ikl, isn) = G2sn(i, isn)          ! [-] [0.0001 m]
+                    dzsnSV(ikl, isn) = dzsn(i, isn)          !            [m]
+                    agsnSV(ikl, isn) = agsn(i, isn)          !          [day]
+                END DO
+                rusnSV(ikl) = rusn(i)              ! Surficial Water
+                toicSV(ikl) = toic(i)              ! bilan snow to ice
+            END IF
+        END DO
 
     END SUBROUTINE sisvatetat0
 
 
-
-
-!======================================================================
-    SUBROUTINE sisvatredem (fichnom,ikl2i,rlon,rlat)
-    
-     
-    
-!======================================================================
-! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
-! Objet: Ecriture de l'etat de redemarrage pour SISVAT
-!======================================================================
-    USE mod_grid_phy_lmdz
-    USE mod_phys_lmdz_para
-    USE iostart
-    USE VAR_SV
-    USE VARxSV          
-    USE VARySV !hj tmp 12 03 2010
-    USE VARtSV
-    USE indice_sol_mod
-    USE dimphy
-
-
-    IMPLICIT none
-
-    include "netcdf.inc"
-!    include "indicesol.h" 
-!    include "dimsoil.h"
-    include "clesphys.h" 
-    include "thermcell.h"
-    include "compbl.h"
-
-!======================================================================
-
-    CHARACTER(LEN=*) :: fichnom
-    INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
-    REAL, DIMENSION(klon), INTENT(IN) :: rlon
-    REAL, DIMENSION(klon), INTENT(IN) :: rlat
-
-! les variables globales ecrites dans le fichier restart
-    REAL, DIMENSION(klon) :: isno
-    REAL, DIMENSION(klon) :: ispi
-    REAL, DIMENSION(klon) :: iice
-    REAL, DIMENSION(klon, nsnowmx) :: isto
-
-    REAL, DIMENSION(klon, nsismx) :: Tsis
-    REAL, DIMENSION(klon, nsismx) :: eta
-    REAL, DIMENSION(klon, nsnowmx) :: dzsn
-    REAL, DIMENSION(klon, nsismx) :: ro        
-    REAL, DIMENSION(klon, nsnowmx) :: G1sn
-    REAL, DIMENSION(klon, nsnowmx) :: G2sn
-    REAL, DIMENSION(klon, nsnowmx) :: agsn
-    REAL, DIMENSION(klon) :: IRs
-    REAL, DIMENSION(klon) :: LMO
-    REAL, DIMENSION(klon) :: rusn 
-    REAL, DIMENSION(klon) :: toic
-    REAL, DIMENSION(klon) :: Bufs
-    REAL, DIMENSION(klon) :: alb1,alb2,alb3
-
-    INTEGER isl, ikl, i, isn, ierr
-    CHARACTER (len=2) :: str2
-    INTEGER           :: pass
-
-      isno(:)       = 0       
-      ispi(:)       = 0 
-      iice(:)       = 0                                
-      IRs(:)        = 0.
-      LMO(:)        = 0.                             
-      eta(:,:)      = 0.     
-      Tsis(:,:)     = 0.         
-      isto(:,:)     = 0 
-      ro(:,:)       = 0.        
-      G1sn(:,:)     = 0.    
-      G2sn(:,:)     = 0. 
-      dzsn(:,:)     = 0.        
-      agsn(:,:)     = 0. 
-      rusn(:)       = 0.   
-      toic(:)       = 0.   
-      Bufs(:)       = 0.   
-      alb1(:)       = 0.
-      alb2(:)       = 0.
-      alb3(:)       = 0.
-
-!***************************************************************************
-! Uncompress SISVAT output variables for storage
-           
-
-      print*, 'je rentre dans restart inlandsis'     
-      DO  ikl = 1,klon  
-           i   = ikl2i(ikl)
-      IF (i > 0) THEN
-        isno(i)       = 1.*isnoSV(ikl)               ! Nb Snow/Ice Lay.   
-        ispi(i)       = 1.*ispiSV(ikl)               ! Nb Supr.Ice Lay.   
-        iice(i)       = 1.*iiceSV(ikl)               ! Nb      Ice Lay.       
-   
-!        IRs(i)        = IRs_SV(ikl)
-!        LMO(i)        = LMO_SV(ikl)                             
-
-
-        DO isl =   -nsol,0                           !                    
-          eta(i,nsno+1-isl)  = eta_SV(ikl,isl)            ! Soil Humidity      
-          Tsis(i,nsno+1-isl) = TsisSV(ikl,isl)            ! Soil Temperature   
-          ro(i,nsno+1-isl)   = ro__SV(ikl,isl)            !        [kg/m3]    
-        END DO        
- 
-  
-        DO  isn = 1,nsno              
-          isto(i,isn)   = 1.*istoSV(ikl,isn)         ! Snow     History 
-          ro(i,isn)     = ro__SV(ikl,isn)            !        [kg/m3]     
-          eta(i,isn)    = eta_SV(ikl,isn)            !        [m3/m3]  
-          Tsis(i,isn)   = TsisSV(ikl,isn)            !            [K]    
-          G1sn(i,isn)   = G1snSV(ikl,isn)            ! [-]        [-]     
-          G2sn(i,isn)   = G2snSV(ikl,isn)            ! [-] [0.0001 m] 
-          dzsn(i,isn)   = dzsnSV(ikl,isn)            !            [m]         
-          agsn(i,isn)   = agsnSV(ikl,isn)            !          [day]     
-        END DO  
-        rusn(i)       = rusnSV(ikl)                  ! Surficial Water  
-        toic(i)       = toicSV(ikl)                  ! to ice 
-        alb1(i)       = alb1sv(ikl) 
-        alb2(i)       = alb2sv(ikl)
-        alb3(i)       = alb3sv(ikl)
-!        Bufs(i)       = BufsSV(ikl)      
-      END IF                     
-      END DO                                                
-
-
-      print*, 'je call open_restart'     
-
-      CALL open_restartphy(fichnom)
-
-      print*, 'je sors open_restart'     
-
-
-      DO pass = 1, 2
-        CALL put_field(pass,"longitude", &
-                    "Longitudes de la grille physique",rlon)     
-        CALL put_field(pass,"latitude","Latitudes de la grille physique",rlat)
-  
-        CALL put_field(pass,"n_snows", "number of snow/ice layers",isno)
-        CALL put_field(pass,"n_ice_top", "number of top ice layers",ispi)
-        CALL put_field(pass,"n_ice", "number of ice layers",iice)
-        CALL put_field(pass,"IR_soil", "Soil IR flux",IRs)
-        CALL put_field(pass,"LMO", "Monin-Obukhov Scale",LMO)
-        CALL put_field(pass,"surf_water", "Surficial water",rusn)
-        CALL put_field(pass,"snow_buffer", "Snow buffer layer",Bufs)
-        CALL put_field(pass,"alb_1", "albedo sw",alb1)
-        CALL put_field(pass,"alb_2", "albedo nIR",alb2)
-        CALL put_field(pass,"alb_3", "albedo fIR",alb3)
-        CALL put_field(pass,"to_ice", "Snow passed to ice",toic)
-
-
-
-        DO isn = 1,nsno
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"AGESNOW"//str2, &
-                         "Age de la neige layer No."//str2, &
-                         agsn(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
+    !======================================================================
+    SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat)
+
+
+
+        !======================================================================
+        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
+        ! Objet: Ecriture de l'etat de redemarrage pour SISVAT
+        !======================================================================
+        USE mod_grid_phy_lmdz
+        USE mod_phys_lmdz_para
+        USE iostart
+        USE VAR_SV
+        USE VARxSV
+        USE VARySV !hj tmp 12 03 2010
+        USE VARtSV
+        USE indice_sol_mod
+        USE dimphy
+
+        IMPLICIT none
+
+        include "netcdf.inc"
+        !    include "indicesol.h"
+        !    include "dimsoil.h"
+        include "clesphys.h"
+        include "thermcell.h"
+        include "compbl.h"
+
+        !======================================================================
+
+        CHARACTER(LEN = *) :: fichnom
+        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
+        REAL, DIMENSION(klon), INTENT(IN) :: rlon
+        REAL, DIMENSION(klon), INTENT(IN) :: rlat
+
+        ! les variables globales ecrites dans le fichier restart
+        REAL, DIMENSION(klon) :: isno
+        REAL, DIMENSION(klon) :: ispi
+        REAL, DIMENSION(klon) :: iice
+        REAL, DIMENSION(klon, nsnowmx) :: isto
+
+        REAL, DIMENSION(klon, nsismx) :: Tsis
+        REAL, DIMENSION(klon, nsismx) :: eta
+        REAL, DIMENSION(klon, nsnowmx) :: dzsn
+        REAL, DIMENSION(klon, nsismx) :: ro
+        REAL, DIMENSION(klon, nsnowmx) :: G1sn
+        REAL, DIMENSION(klon, nsnowmx) :: G2sn
+        REAL, DIMENSION(klon, nsnowmx) :: agsn
+        REAL, DIMENSION(klon) :: IRs
+        REAL, DIMENSION(klon) :: LMO
+        REAL, DIMENSION(klon) :: rusn
+        REAL, DIMENSION(klon) :: toic
+        REAL, DIMENSION(klon) :: Bufs
+        REAL, DIMENSION(klon) :: alb1, alb2, alb3
+
+        INTEGER isl, ikl, i, isn, ierr
+        CHARACTER (len = 2) :: str2
+        INTEGER :: pass
+
+        isno(:) = 0
+        ispi(:) = 0
+        iice(:) = 0
+        IRs(:) = 0.
+        LMO(:) = 0.
+        eta(:, :) = 0.
+        Tsis(:, :) = 0.
+        isto(:, :) = 0
+        ro(:, :) = 0.
+        G1sn(:, :) = 0.
+        G2sn(:, :) = 0.
+        dzsn(:, :) = 0.
+        agsn(:, :) = 0.
+        rusn(:) = 0.
+        toic(:) = 0.
+        Bufs(:) = 0.
+        alb1(:) = 0.
+        alb2(:) = 0.
+        alb3(:) = 0.
+
+        !***************************************************************************
+        ! Uncompress SISVAT output variables for storage
+
+        DO  ikl = 1, klon
+            i = ikl2i(ikl)
+            IF (i > 0) THEN
+                isno(i) = 1. * isnoSV(ikl)               ! Nb Snow/Ice Lay.
+                ispi(i) = 1. * ispiSV(ikl)               ! Nb Supr.Ice Lay.
+                iice(i) = 1. * iiceSV(ikl)               ! Nb      Ice Lay.
+
+                !        IRs(i)        = IRs_SV(ikl)
+                !        LMO(i)        = LMO_SV(ikl)
+
+                DO isl = -nsol, 0                           !
+                    eta(i, nsno + 1 - isl) = eta_SV(ikl, isl)            ! Soil Humidity
+                    Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl)            ! Soil Temperature
+                    ro(i, nsno + 1 - isl) = ro__SV(ikl, isl)            !        [kg/m3]
+                END DO
+
+                DO  isn = 1, nsno
+                    isto(i, isn) = 1. * istoSV(ikl, isn)         ! Snow     History
+                    ro(i, isn) = ro__SV(ikl, isn)            !        [kg/m3]
+                    eta(i, isn) = eta_SV(ikl, isn)            !        [m3/m3]
+                    Tsis(i, isn) = TsisSV(ikl, isn)            !            [K]
+                    G1sn(i, isn) = G1snSV(ikl, isn)            ! [-]        [-]
+                    G2sn(i, isn) = G2snSV(ikl, isn)            ! [-] [0.0001 m]
+                    dzsn(i, isn) = dzsnSV(ikl, isn)            !            [m]
+                    agsn(i, isn) = agsnSV(ikl, isn)            !          [day]
+                END DO
+                rusn(i) = rusnSV(ikl)                  ! Surficial Water
+                toic(i) = toicSV(ikl)                  ! to ice
+                alb1(i) = alb1sv(ikl)
+                alb2(i) = alb2sv(ikl)
+                alb3(i) = alb3sv(ikl)
+                !        Bufs(i)       = BufsSV(ikl)
+            END IF
+        END DO
+
+        CALL open_restartphy(fichnom)
+
+        DO pass = 1, 2
+            CALL put_field(pass, "longitude", &
+                    "Longitudes de la grille physique", rlon)
+            CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat)
+
+            CALL put_field(pass, "n_snows", "number of snow/ice layers", isno)
+            CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi)
+            CALL put_field(pass, "n_ice", "number of ice layers", iice)
+            CALL put_field(pass, "IR_soil", "Soil IR flux", IRs)
+            CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO)
+            CALL put_field(pass, "surf_water", "Surficial water", rusn)
+            CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs)
+            CALL put_field(pass, "alb_1", "albedo sw", alb1)
+            CALL put_field(pass, "alb_2", "albedo nIR", alb2)
+            CALL put_field(pass, "alb_3", "albedo fIR", alb3)
+            CALL put_field(pass, "to_ice", "Snow passed to ice", toic)
+
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "AGESNOW" // str2, &
+                            "Age de la neige layer No." // str2, &
+                            agsn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "DZSNOW" // str2, &
+                            "Snow/ice thickness layer No." // str2, &
+                            dzsn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "G2SNOW" // str2, &
+                            "Snow Property 2, layer No." // str2, &
+                            G2sn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "G1SNOW" // str2, &
+                            "Snow Property 1, layer No." // str2, &
+                            G1sn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsismx
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "ETA" // str2, &
+                            "Soil/snow water content layer No." // str2, &
+                            eta(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsismx   !nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "RO" // str2, &
+                            "Snow density layer No." // str2, &
+                            ro(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsismx
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "TSS" // str2, &
+                            "Soil/snow temperature layer No." // str2, &
+                            Tsis(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "HISTORY" // str2, &
+                            "Snow history layer No." // str2, &
+                            isto(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+
+            CALL enddef_restartphy
         ENDDO
-        DO isn = 1,nsno
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"DZSNOW"//str2, &
-                         "Snow/ice thickness layer No."//str2, &
-                         dzsn(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
-        ENDDO
-        DO isn = 1,nsno
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"G2SNOW"//str2, &
-                         "Snow Property 2, layer No."//str2, &
-                         G2sn(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
-        ENDDO
-        DO isn = 1,nsno
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"G1SNOW"//str2, &
-                         "Snow Property 1, layer No."//str2, &
-                         G1sn(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
-        ENDDO
-        DO isn = 1,nsismx
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"ETA"//str2, &
-                         "Soil/snow water content layer No."//str2, &
-                         eta(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
-        ENDDO
-        DO isn = 1,nsismx   !nsno
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"RO"//str2, &
-                           "Snow density layer No."//str2, &
-                           ro(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
-        ENDDO
-        DO isn = 1,nsismx
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"TSS"//str2, &
-                           "Soil/snow temperature layer No."//str2, &
-                           Tsis(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
-        ENDDO
-        DO isn = 1,nsno
-          IF (isn.LE.99) THEN
-            WRITE(str2,'(i2.2)') isn
-            CALL put_field(pass,"HISTORY"//str2, &
-                           "Snow history layer No."//str2, &
-                           isto(:,isn))
-          ELSE
-            PRINT*, "Trop de couches"
-            CALL abort
-          ENDIF
-        ENDDO
-
-      CALL enddef_restartphy
-      ENDDO
-      CALL close_restartphy
-
-
-  END SUBROUTINE sisvatredem
+        CALL close_restartphy
+
+    END SUBROUTINE sisvatredem
 
 END MODULE surf_inlandsis_mod

LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

@@ -182 +182 @@
        debut,     lafin,                              &
        rlon,      rlat,      rugoro,   rmu0,          &
-       zsig,      lwdown_m,  pphi,     cldt,          &
+       lwdown_m,  cldt,          &
        rain_f,    snow_f,    solsw_m,  solswfdiff_m, sollw_m,       &
        gustiness,                                     &
@@ -277 +277 @@
 ! z0m, z0h ----input-R- longeur de rugosite (en m)
 ! Martin
-! zsig-----input-R- slope
 ! cldt-----input-R- total cloud fraction
-! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)
 ! Martin
 !
@@ -315 +313 @@
     USE print_control_mod,  ONLY : prt_level,lunout
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal
+    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, zsig, zmea
     use phys_output_var_mod, only: dter, dser, tkt, tks, taur, sss
 #ifdef CPP_XIOS
@@ -322 +320 @@
     use netcdf, only: missing_val => nf90_fill_real
 #endif
+
+     
+
 
     IMPLICIT NONE
@@ -359 +360 @@
     REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf  ! sub-surface fraction
 ! Martin
-    REAL, DIMENSION(klon),        INTENT(IN)        :: zsig    ! slope
     REAL, DIMENSION(klon),        INTENT(IN)        :: lwdown_m ! downward longwave radiation at mean s    
     REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
 
     REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud fraction
-    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pphi    ! geopotential (m2/s2)
-! Martin
 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -873 +871 @@
     REAL, DIMENSION(klon)              :: ytoice
     REAL, DIMENSION(klon)              :: ysnowhgt, yqsnow, ysissnow, yrunoff
+    REAL, DIMENSION(klon)              :: yzmea
     REAL, DIMENSION(klon)              :: yzsig
-    REAL, DIMENSION(klon,klev)         :: ypphi
     REAL, DIMENSION(klon)              :: ycldt
     REAL, DIMENSION(klon)              :: yrmu0
@@ -1054 +1052 @@
     ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
     yalb3_new = 0.0  ; ysissnow = 0.0 
-    ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
+    ycldt = 0.0      ; yrmu0 = 0.0
     ! Martin
 
@@ -1394 +1392 @@
           ywindsp(j) = windsp(i,nsrf)
 !>jyg
-          ! Martin
+          ! Martin and Etienne
+          yzmea(j)   = zmea(i)
           yzsig(j)   = zsig(i)
           ycldt(j)   = cldt(i)
@@ -1529 +1528 @@
 !
 !****************************************************************************************
+
 
 !!! jyg le 07/02/2012
@@ -1579 +1579 @@
            speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2)
         ENDDO
-        CALL cdrag(knon, nsrf, &
+
+
+            CALL cdrag(knon, nsrf, &
             speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
             yts_x, yqsurf_x, yz0m, yz0h, &
@@ -2040 +2042 @@
           CALL surf_landice(itap, dtime, knon, ni, &
                rlon, rlat, debut, lafin, &
-               yrmu0, ylwdown, yalb, ypphi(:,1), &
+               yrmu0, ylwdown, yalb, zgeo1, &
                ysolsw, ysollw, yts, ypplay(:,1), &
 !!jyg               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
@@ -2050 +2052 @@
                ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
                ytsurf_new, y_dflux_t, y_dflux_q, &
-               yzsig, ycldt, &
+               yzmea, yzsig, ycldt, &
                ysnowhgt, yqsnow, ytoice, ysissnow, &
                yalb3_new, yrunoff, &
@@ -2840 +2842 @@
           sss(ni(:knon)) = ysss(:knon)
        end if
+
 
 !****************************************************************************************
@@ -3253 +3256 @@
     ENDDO
 !!!
-   
+
 !
 ! Incrementer la temperature du sol

LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90

@@ -1054 +1054 @@
        CALL histwrite_phy(o_tauy, zx_tmp_fi2d)
 
-       IF (landice_opt .GE. 1) THEN
-          CALL histwrite_phy(o_snowsrf, snow_o)
-          CALL histwrite_phy(o_qsnow, qsnow)
-          CALL histwrite_phy(o_snowhgt,snowhgt)
-          CALL histwrite_phy(o_toice,to_ice)
-          CALL histwrite_phy(o_sissnow,sissnow)
-          CALL histwrite_phy(o_runoff,runoff)
-          CALL histwrite_phy(o_albslw3,albsol3_lic)
-       ENDIF
+       ! Etienne: test sorties pour compil sur JZ
+!       IF (landice_opt .GE. 1) THEN
+!          CALL histwrite_phy(o_snowsrf, snow_o)
+!          CALL histwrite_phy(o_qsnow, qsnow)
+!          CALL histwrite_phy(o_snowhgt,snowhgt)
+!          CALL histwrite_phy(o_toice,to_ice)
+!          CALL histwrite_phy(o_sissnow,sissnow)
+!          CALL histwrite_phy(o_runoff,runoff)
+!          CALL histwrite_phy(o_albslw3,albsol3_lic)
+!       ENDIF
 
        DO nsrf = 1, nbsrf

LMDZ6/trunk/libf/phylmd/physiq_mod.F90

@@ -2537 +2537 @@
             debut,     lafin, &
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
-            zsig,      sollwdown, pphi,    cldt,      &
+             sollwdown,    cldt,      &
             rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
             gustiness,                                &

LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90

@@ -19 +19 @@
        tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
        tsurf_new, dflux_s, dflux_l, &
-       slope, cloudf, &
+       alt, slope, cloudf, &
        snowhgt, qsnow, to_ice, sissnow, &
        alb3, runoff, &
@@ -25 +25 @@
 
     USE dimphy
-    USE surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, n_dtis
-    USE fonte_neige_mod,  ONLY : fonte_neige, run_off_lic
+    USE surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, iflag_albcalc
+    USE fonte_neige_mod,  ONLY : fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
@@ -75 +75 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: albedo  !mean albedo
     REAL, DIMENSION(klon), INTENT(IN)             :: pphi1    
+    REAL, DIMENSION(klon), INTENT(IN)             :: alt   !mean altitude of the grid box  
     REAL, DIMENSION(klon), INTENT(IN)             :: slope   !mean slope in grid box  
     REAL, DIMENSION(klon), INTENT(IN)             :: cloudf  !total cloud fraction
@@ -115 +116 @@
     REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
     INTEGER                  :: i,j,nt
-
+    REAL, DIMENSION(klon)    :: fqfonte,ffonte
     REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
     REAL, DIMENSION(klon)    :: swdown,lwdown
-    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection
-    REAL, DIMENSION(klon)    :: zsl_height, wind_velo      !surface layer height, wind spd
+    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
+    REAL, DIMENSION(klon)    :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
+    REAL, DIMENSION(klon)    :: zsl_height, wind_velo     !surface layer height, wind spd
     REAL, DIMENSION(klon)    :: dens_air,  snow_cont_air  !air density; snow content air
     REAL, DIMENSION(klon)    :: alb_soil                  !albedo of underlying ice
@@ -132 +134 @@
 
 
-!albedo SB >>>
-    real,dimension(klon) :: alb1,alb2
-!albedo SB <<<
-
+    REAL,DIMENSION(klon) :: alb1,alb2
+    REAL, DIMENSION (klon,6) :: alb6
 ! End definition
 !****************************************************************************************
@@ -230 +230 @@
             AcoefH, AcoefQ, BcoefH, BcoefQ, cdragh, &
             run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &        
-            tsurf_new, alb1, alb2, alb3, &
-            emis_new, z0m, qsurf)
+            tsurf_new, alb1, alb2, alb3, emis_new, z0m, qsurf)
        z0h(1:knon)=z0m(1:knon) ! en attendant mieux
        
@@ -278 +277 @@
        swdown(:)        = 0.0
        lwdown(:)        = 0.0
-       snow_adv(:)      = 0.  ! no snow blown in for now
-       snow_cont_air(:) = 0.       
+       snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model      
        alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
        ustar(:)         = 0.
@@ -296 +294 @@
        
 
-       ! Subtimestepping
-       
-        dtis=dtime/n_dtis
-
-        DO nt=1,n_dtis
-
-          IF (lafin .and. nt.eq.n_dtis) THEN
+
+        dtis=dtime
+
+          IF (lafin) THEN
             lafin_is=.true.
           END IF
 
-        !PRINT*,'RENTRE DANS INLANDSIS','itime',itime,'dtime',dtime,'dtis',dtis
-        CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is, &
-            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, &
-            precip_rain, precip_snow, precip_snow_adv, snow_adv, &
-            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
-            rugoro, snow_cont_air, alb_soil, slope, cloudf, &
-            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice,sissnow, agesno, &
+          CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
+            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
+            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
+            rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
+            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
             AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
-            run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &        
-            tsurf_new, alb1, alb2, alb3, &
-            emis_new, z0m, z0h, qsurf)      
-
-         debut_is=.false.
-
-        END DO
+            run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
+            tsurf_new, alb1, alb2, alb3, alb6, &
+            emis_new, z0m, z0h, qsurf)
+
+          debut_is=.false.
+
+
+        ! Treatment of snow melting and calving
+
+        ! for consistency with standard LMDZ, add calving to run_off_lic
+        run_off_lic(:)=run_off_lic(:) + to_ice(:)
+
+        DO i = 1, knon
+           ffonte_global(knindex(i),is_lic)    = ffonte(i)
+           fqfonte_global(knindex(i),is_lic)   = fqfonte(i)! net melting= melting - refreezing
+           fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
+           runofflic_global(knindex(i)) = run_off_lic(i)
+        ENDDO
+        ! Here, we assume that the calving term is equal to the to_ice term
+        ! (no ice accumulation)
 
 
@@ -419 +425 @@
 
 
-
-
-    
 
     END IF ! landice_opt
@@ -476 +479 @@
        alb_dir(1:knon,5)=alb2(1:knon)
        alb_dir(1:knon,6)=alb2(1:knon)
+
+       IF ((landice_opt .EQ.2) .AND. (iflag_albcalc .EQ. 2)) THEN
+       alb_dir(1:knon,1)=alb6(1:knon,1)
+       alb_dir(1:knon,2)=alb6(1:knon,2)
+       alb_dir(1:knon,3)=alb6(1:knon,3)
+       alb_dir(1:knon,4)=alb6(1:knon,4)
+       alb_dir(1:knon,5)=alb6(1:knon,5)
+       alb_dir(1:knon,6)=alb6(1:knon,6)
+       ENDIF
+
      end select
 alb_dif=alb_dir
 !albedo SB <<<
 
-
-
+  
+  
 
   END SUBROUTINE surf_landice

LMDZ6/trunk/libf/phylmd/surface_data.F90

@@ -29 +29 @@
   ! FOR INLANDSIS:
   !===============
-
-   INTEGER, SAVE          :: landice_opt       ! 1 for coupling with SISVAT, 2 for coupling with INLANDSIS
+ 
+   ! 1 for coupling with SISVAT, 2 for coupling with INLANDSIS and number of subtimesteps for INLANDSIS
+   INTEGER, SAVE          :: landice_opt      ! 1 for coupling with SISVAT, 2 for coupling with INLANDSIS
   !$OMP THREADPRIVATE(landice_opt)
 
-   INTEGER, SAVE          :: iflag_tsurf_inlandsis ! 0 SISVAT method, 1 LMDZ method
-  !$OMP THREADPRIVATE(iflag_tsurf_inlandsis)
+  ! temperature calculation options within the soil and at the surface
+   INTEGER, SAVE          :: iflag_tsurf_inlandsis,iflag_temp_inlandsis 
+   !$OMP THREADPRIVATE(iflag_tsurf_inlandsis,iflag_temp_inlandsis)
 
-   INTEGER, SAVE          :: iflag_albzenith ! dependency of albedo to zenith angle
-  !$OMP THREADPRIVATE(iflag_albzenith)  
-
-   INTEGER, SAVE          :: n_dtis           ! number of subtimesteps for INLANDSIS
-  !$OMP THREADPRIVATE(n_dtis)
+  ! flags for albedo and roughness calc.
+   INTEGER, SAVE          :: iflag_albcalc,iflag_z0m_snow
+  !$OMP THREADPRIVATE(iflag_albcalc,iflag_z0m_snow)  
 
   ! with or without snow module/ blowing snow, ascii outfile
-   LOGICAL, SAVE           :: SnoMod,BloMod,ok_outfor 
+   LOGICAL, SAVE          :: SnoMod,BloMod,ok_outfor 
   !$OMP THREADPRIVATE(SnoMod,BloMod,ok_outfor)    
 
+  ! activate slush, korlyakov snow density, RN z0h calc.
+   LOGICAL, SAVE          :: is_ok_slush,is_ok_density_kotlyakov,is_ok_z0h_rn
+  !$OMP THREADPRIVATE(is_ok_slush,is_ok_density_kotlyakov,is_ok_z0h_rn) 
+
+  ! activate detection snow/ice layers and option XF discrtet/option runoff AC
+   LOGICAL, SAVE          :: ok_zsn_ii,discret_xf,opt_runoff_ac
+  !$OMP THREADPRIVATE(ok_zsn_ii,discret_xf, opt_runoff_ac)
+
+  ! value of z0m snow when prescribed and albedo correction term
+   REAL, SAVE             :: prescribed_z0m_snow,correc_alb
+  !$OMP THREADPRIVATE(prescribed_z0m_snow, correc_alb)
+
+  ! value of sphericity [0-99] and snow grain size [e-4m] for polar buffer snow
+  ! layer
+   REAL, SAVE             :: buf_sph_pol,buf_siz_pol
+  !$OMP THREADPRIVATE(buf_sph_pol,buf_siz_pol)
+
+
+
 END MODULE surface_data