subroutine INLANDSIS(SnoMod,BloMod,jjtime,debut)

  USE dimphy

  !--------------------------------------------------------------------------+
  ! INLANDSIS module                                                     |
  ! Simplified SISVAT module, containing ice and snow processes for      |
  ! ice-covered surfaces                                                 |
  ! version MARv3, november 2020                                         |
  ! SubRoutine INLANDSIS contains the fortran 77 code of the             |
  !            Soil/Ice Snow Vegetation Atmosphere Transfer Scheme       |
  !                                                                      |
  !--------------------------------------------------------------------------+
  ! PARAMETERS:  klonv: Total Number of columns =                        |
  ! ^^^^^^^^^^        = Total Number of continental     grid boxes       |
  !                   X       Number of Mosaic Cell per grid box         |
  !                                                                      |
  ! INPUT:   daHost   : Date Host Model                                  |
  ! ^^^^^                                                                |
  !                                                                      |
  ! INPUT:   LSmask   : 1:          Land       MASK                      |
  ! ^^^^^               0:          Sea        MASK                      |
  !          isotSV   = 0,...,12:   Soil       Type                      |
  !                     0:          Water,          Liquid (Sea, Lake)   |
  !                    12:          Water, Solid           (Ice)         |
  !                                                                      |
  ! INPUT:   coszSV   : Cosine of the Sun Zenithal Distance          [-] |
  ! ^^^^^    sol_SV   : Surface Downward  Solar      Radiation    [W/m2] |
  !          IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |
  !          drr_SV   : Rain  Intensity                        [kg/m2/s] |
  !          dsn_SV   : Snow  Intensity                      [mm w.e./s] |
  !          dsnbSV   : Snow  Intensity,  Drift Fraction             [-] |
  !          dbs_SV   : Drift Amount                           [mm w.e.] |
  !          za__SV   : Surface Boundary Layer (SBL) Height          [m] |
  !          VV__SV   :(SBL Top)   Wind Velocity                   [m/s] |
  !          TaT_SV   : SBL Top    Temperature                       [K] |
  !          rhT_SV   : SBL Top    Air  Density                  [kg/m3] |
  !          QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |
  !          qsnoSV   : SBL Mean   Snow      Content             [kg/kg] |
  !          alb0SV   : Soil Basic Albedo                            [-] |
  !          slopSV   : Surface    Slope                             [-] |
  !          dt__SV   : Time  Step                                   [s] |
  !                                                                      |
  ! INPUT /  isnoSV   = total Nb of Ice/Snow Layers                      |
  ! OUTPUT:  ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
  ! ^^^^^^   iiceSV   = total Nb of Ice      Layers                      |
  !          istoSV   = 0,...,5 :   Snow     History (see istdSV data)   |
  !                                                                      |
  ! INPUT /  alb_SV   : Surface        Albedo                        [-] |
  ! OUTPUT:  emi_SV   : Surface        Emissivity                    [-] |
  ! ^^^^^^   IRs_SV   : Soil           IR Flux  (negative)        [W/m2] |
  !          LMO_SV   : Monin-Obukhov               Scale            [m] |
  !          us__SV   : Friction          Velocity                 [m/s] |
  !          uts_SV   : Temperature       Turbulent Scale          [m/s] |
  !          uqs_SV   : Specific Humidity Velocity                 [m/s] |
  !          uss_SV   : Blowing Snow      Turbulent Scale          [m/s] |
  !          usthSV   : Blowing Snow      Erosion   Threshold      [m/s] |
  !          Z0m_SV   : Momentum     Roughness Length                [m] |
  !          Z0mmSV   : Momentum     Roughness Length (time mean)    [m] |
  !          Z0mnSV   : Momentum     Roughness Length (instantaneous)[m] |
  !          Z0SaSV   : Sastrugi     Roughness Length                [m] |
  !          Z0e_SV   : Erosion Snow Roughness Length                [m] |
  !          Z0emSV   : Erosion Snow Roughness Length (time mean)    [m] |
  !          Z0enSV   : Erosion Snow Roughness Length (instantaneous)[m] |
  !          Z0roSV   : Subgrid Topo Roughness Length                [m] |
  !          Z0h_SV   : Heat         Roughness Length                [m] |
  !          TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
  !                   & Snow     Temperatures (layers  1,2,...,nsno) [K] |
  !          ro__SV   : Soil/Snow Volumic Mass                   [kg/m3] |
  !          eta_SV   : Soil/Snow Water   Content                [m3/m3] |
  !          G1snSV   : snow dendricity/sphericity                       |
  !          G2snSV   : snow sphericity/grain size                       |
  !          dzsnSV   : Snow Layer        Thickness                  [m] |
  !          agsnSV   : Snow       Age                             [day] |
  !          BufsSV   : Snow Buffer Layer              [kg/m2] .OR. [mm] |
  !          BrosSV   : Snow Buffer Layer Density      [kg/m3]           |
  !          BG1sSV   : Snow Buffer Layer Dendricity / Sphericity    [-] |
  !          BG2sSV   : Snow Buffer Layer Sphericity / Size [-] [0.1 mm] |
  !          rusnSV   : Surficial   Water              [kg/m2] .OR. [mm] |
  !                                                                      |
  ! OUTPUT:  no__SV   : OUTPUT file Unit Number                      [-] |
  ! ^^^^^^   i___SV   : OUTPUT point   i Coordinate                  [-] |
  !          j___SV   : OUTPUT point   j Coordinate                  [-] |
  !          n___SV   : OUTPUT point   n Coordinate                  [-] |
  !          lwriSV   : OUTPUT point vec Index                       [-] |
  !                                                                      |
  ! OUTPUT:  IRu_SV   : Upward     IR Flux (+, upw., effective)      [K] |
  ! ^^^^^^   hSalSV   : Saltating Layer Height                       [m] |
  !          qSalSV   : Saltating Snow  Concentration            [kg/kg] |
  !          RnofSV   : RunOFF Intensity                       [kg/m2/s] |
  !                                                                      |
  ! Internal Variables:                                                  |
  ! ^^^^^^^^^^^^^^^^^^                                                   |
  !          NLaysv   = New            Snow Layer Switch             [-] |
  !          albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
  !          SoSosv   : Absorbed Solar Radiation by Surfac.(Normaliz)[-] |
  !          TBr_sv   : Brightness Temperature                       [K] |
  !          IRupsv   : Upward     IR Flux (-, upw.)              [W/m2] |
  !          ram_sv   : Aerodynamic Resistance for Momentum        [s/m] |
  !          rah_sv   : Aerodynamic Resistance for Heat            [s/m] |
  !          Evp_sv   : Evaporation                              [kg/m2] |
  !          EvT_sv   : Evapotranspiration                       [kg/m2] |
  !          HSs_sv   : Surface    Sensible Heat Flux + => absorb.[W/m2] |
  !          HLs_sv   : Surface    Latent   Heat Flux + => absorb.[W/m2] |
  !          Lx_H2O   : Latent Heat of Vaporization/Sublimation   [J/kg] |
  !          Tsrfsv   : Surface    Temperature                       [K] |
  !          sEX_sv   : Verticaly Integrated Extinction Coefficient  [-] |
  !          LSdzsv   : Vertical   Discretization Factor             [-] |
  !                   =    1. Soil                                       |
  !                   = 1000. Ocean                                      |
  !          z_snsv   : Snow Pack  Thickness                         [m] |
  !          zzsnsv   : Snow Pack  Thickness                         [m] |
  !          albssv   : Soil       Albedo                            [-] |
  !          Eso_sv   : Soil+Snow       Emissivity                   [-] |
  !          Khydsv   : Soil   Hydraulic    Conductivity           [m/s] |
  !                                                                      |
  !          ETSo_0   : Snow/Soil Energy Power, before Forcing    [W/m2] |
  !          ETSo_1   : Snow/Soil Energy Power, after  Forcing    [W/m2] |
  !          ETSo_d   : Snow/Soil Energy Power         Forcing    [W/m2] |
  !          EqSn_0   : Snow      Energy, before Phase Change     [J/m2] |
  !          EqSn_1   : Snow      Energy, after  Phase Change     [J/m2] |
  !          EqSn_d   : Snow      Energy,       net    Forcing    [J/m2] |
  !          Enrsvd   : SVAT      Energy Power         Forcing    [W/m2] |
  !          Enrbal   : SVAT      Energy Balance                  [W/m2] |
  !          Wats_0   : Soil Water,  before Forcing                 [mm] |
  !          Wats_1   : Soil Water,  after  Forcing                 [mm] |
  !          Wats_d   : Soil Water          Forcing                 [mm] |
  !          SIWm_0   : Snow        initial Mass               [mm w.e.] |
  !          SIWm_1   : Snow        final   Mass               [mm w.e.] |
  !          SIWa_i   : Snow Atmos. initial Forcing            [mm w.e.] |
  !          SIWa_f   : Snow Atmos. final   Forcing(noConsumed)[mm w.e.] |
  !          SIWe_i   : SnowErosion initial Forcing            [mm w.e.] |
  !          SIWe_f   : SnowErosion final   Forcing(noConsumed)[mm w.e.] |
  !          SIsubl   : Snow sublimed/deposed  Mass            [mm w.e.] |
  !          SImelt   : Snow Melted            Mass            [mm w.e.] |
  !          SIrnof   : Surficial Water + Run OFF Change       [mm w.e.] |
  !          SIvAcr   : Sea-Ice    vertical Acretion           [mm w.e.] |
  !          Watsvd   : SVAT Water          Forcing                 [mm] |
  !          Watbal   : SVAT Water  Balance                       [W/m2] |
  !                                                                      |
  !          vk2      : Square of Von Karman Constant                [-] |
  !          sqrCm0   : Factor of   Neutral Drag Coeffic.Momentum  [s/m] |
  !          sqrCh0   : Factor of   Neutral Drag Coeffic.Heat      [s/m] |
  !          EmiSol   : Soil          Emissivity                     [-] |
  !          EmiSno   : Snow          Emissivity                     [-] |
  !          EmiWat   : Water         Emissivity                     [-] |
  !          Z0mLnd   :          Land Roughness Length               [m] |
  !          sqrrZ0   : u*t/u*                                           |
  !          f_eff    : Marticorena & B. 1995 JGR (20)                   |
  !          A_Fact   : Fundamental * Roughness                          |
  !          Z0mBSn   :         BSnow Roughness Length               [m] |
  !          Z0mBS0   : Mimimum BSnow Roughness Length (blown* )     [m] |
  !          Z0m_Sn   :          Snow Roughness Length (surface)     [m] |
  !          Z0m_S0   : Mimimum  Snow Roughness Length               [m] |
  !          Z0m_S1   : Maximum  Snow Roughness Length               [m] |
  !          Z0_GIM   : Minimum GIMEX Roughness Length               [m] |
  !          Z0_ICE   : Sea Ice ISW   Roughness Length               [m] |
  !                                                                      |
  !                                                                      |
  !--------------------------------------------------------------------------+



  ! Global Variables
  ! ================


  USE VARphy
  USE VAR_SV
  USE VARdSV
  USE VAR0SV
  USE VARxSV
  USE VARySV
  USE VARtSV
  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

  logical :: SnoMod
  logical :: BloMod
  logical :: debut
  integer :: jjtime


  ! Internal Variables
  ! ==================

  ! Non Local
  ! ---------

  real :: TBr_sv(klonv)                 ! Brightness Temperature
  real :: IRdwsv(klonv)                 ! DOWNward   IR Flux
  real :: IRupsv(klonv)                 ! UPward     IR Flux
  real :: d_Bufs,Bufs_N                 ! Buffer Snow Layer Increment
  real :: Buf_ro,Bros_N                 ! Buffer Snow Layer Density
  real :: BufPro                        ! Buffer Snow Layer Density
  real :: Buf_G1,BG1__N                 ! Buffer Snow Layer Dendr/Sphe[-]
  real :: Buf_G2,BG2__N                 ! Buffer Snow Layer Spher/Size[-]
  real :: Bdzssv(klonv)                 ! Buffer Snow Layer Thickness
  real :: z_snsv(klonv)                 ! Snow-Ice, current Thickness



  ! Local
  ! -----

  integer :: iwr
  integer :: ikl   ,isn   ,isl   ,ist      !
  integer :: ist__s,ist__w                 ! Soil/Water Body Identifier
  integer :: growth                        ! Seasonal               Mask
  integer :: LISmsk                        ! Land+Ice / Open    Sea Mask
  integer :: LSnMsk                        ! Snow-Ice / No Snow-Ice Mask
  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,Salt_f,PorRef   !
  ! #sw real      PorVol,rWater                 !
  ! #sw real      rusNEW,rdzNEW,etaNEW          !
  real :: ro_new                        !
  real :: TaPole                        ! Maximum     Polar Temperature
  real :: T__Min                        ! Minimum realistic Temperature
  real :: EmiSol                        ! Emissivity of       Soil
  real :: EmiSno                        ! Emissivity of            Snow
  real :: EmiWat                        ! Emissivity of a Water Area
  real :: vk2                           ! Square of Von Karman Constant
  real :: u2star                        !(u*)**2
  real :: Z0mLnd                        !          Land Roughness Length
  ! #ZN real      sqrrZ0                        ! u*t/u*
  real :: f_eff                         ! Marticorena & B. 1995 JGR (20)
  real :: A_Fact                        ! Fundamental * Roughness
  real :: Z0m_nu                        ! Smooth R Snow Roughness Length
  real :: Z0mBSn                        !         BSnow Roughness Length
  real :: Z0mBS0                        ! Mimimum BSnow Roughness Length
  real :: Z0m_S0                        ! Mimimum  Snow Roughness Length
  real :: Z0m_S1                        ! Maximum  Snow Roughness Length
  ! #SZ real      Z0Sa_N                        ! Regime   Snow Roughness Length
  ! #SZ real      Z0SaSi                        ! 1.IF Rgm Snow Roughness Length
  ! #GL real      Z0_GIM                        ! Mimimum GIMEX Roughness Length
  real :: Z0_ICE                        ! Ice ISW   Roughness Length
  real :: Z0m_Sn,Z0m_90                 ! Snow  Surface Roughness Length
  real :: SnoWat                        ! Snow Layer    Switch
  real :: rstar,alors                   !
  real :: rstar0,rstar1,rstar2          !
  real :: SameOK                        ! 1. => Same Type of Grains
  real :: G1same                        ! Averaged G1,  same Grains
  real :: G2same                        ! Averaged G2,  same Grains
  real :: typ__1                        ! 1. => Lay1 Type: Dendritic
  real :: zroNEW                        ! dz X ro, if fresh Snow
  real :: G1_NEW                        ! G1,      if fresh Snow
  real :: G2_NEW                        ! G2,      if fresh Snow
  real :: zroOLD                        ! dz X ro, if old   Snow
  real :: G1_OLD                        ! G1,      if old   Snow
  real :: G2_OLD                        ! G2,      if old   Snow
  real :: SizNEW                        ! Size,    if fresh Snow
  real :: SphNEW                        ! Spheric.,if fresh Snow
  real :: SizOLD                        ! Size,    if old   Snow
  real :: SphOLD                        ! Spheric.,if old   Snow
  real :: Siz_av                        ! Averaged    Grain Size
  real :: Sph_av                        ! Averaged    Grain Spher.
  real :: Den_av                        ! Averaged    Grain Dendr.
  real :: G1diff                        ! Averaged G1, diff. Grains
  real :: G2diff                        ! Averaged G2, diff. Grains
  real :: G1                            ! Averaged G1
  real :: G2                            ! Averaged G2
  real :: param                         ! Polynomial   fit z0=f(T)
  real :: Z0_obs                        ! Fit Z0_obs=f(T) (m)
  real :: tamin                         ! min T of linear fit (K)
  real :: tamax                         ! max T of linear fit (K)
  real :: coefa,coefb,coefc,coefd       ! Coefs for z0=f(T)
  real :: ta1,ta2,ta3                   ! Air temperature thresholds
  real :: z01,z02,z03                   ! z0 thresholds
  real :: tt_c,vv_c                     ! Critical param.
  real :: tt_tmp,vv_tmp,vv_virt         ! Temporary variables
  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(len=3) :: qsalt_param              ! Switch for saltation flux param.
  character(len=3) :: usth_param               ! Switch for u*t param


  ! Internal DATA
  ! =============

  data      T__Min / 200.00/              ! Minimum realistic Temperature
  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)
  data      vv_c   / 14.3  /              ! Critical Wind speed (m/s)
  data      roSn_1 / 109.  /              ! Fall.Sno.Density, Indep. Param.
  data      roSn_2 /   6.  /              ! Fall.Sno.Density, Temper.Param.
  data      roSn_3 /  26.  /              ! Fall.Sno.Density, Wind   Param.
  data      Dendr1 /  17.12/              ! Fall.Sno.Dendric.,Wind 1/Param.
  data      Dendr2 / 128.  /              ! Fall.Sno.Dendric.,Wind 2/Param.
  data      Dendr3 / -20.  /              ! Fall.Sno.Dendric.,Indep. Param.
  data      Spher1 /   7.87/              ! Fall.Sno.Spheric.,Wind 1/Param.
  data      Spher2 /  38.  /              ! Fall.Sno.Spheric.,Wind 2/Param.
  data      Spher3 /  50.  /              ! Fall.Sno.Spheric.,Wind 3/Param.
  data      Spher4 /  90.  /              ! Fall.Sno.Spheric.,Indep. Param.
  data      EmiSol /   0.99999999/        ! 0.94Emissivity of Soil
  data      EmiWat /   0.99999999/        ! Emissivity of a Water Area
  data      EmiSno /   0.99999999/        ! Emissivity of Snow


  ! DATA      Emissivities                  ! Pielke, 1984, pp. 383,409

  data      Z0mBS0 /   0.5e-6/            ! MINimum Snow Roughness Length
                                          ! ! for Momentum if Blowing Snow
                                          ! ! Gallee et al. 2001 BLM 99 (19)
  data      Z0m_S0/    0.00005/           ! MINimum Snow Roughness Length
                                          ! ! MegaDunes    included
  data      Z0m_S1/    0.030  /           ! MAXimum Snow Roughness Length
                                          ! !        (Sastrugis)
  ! #GL data      Z0_GIM/    0.0013/            ! Ice Min Z0 = 0.0013 m (Broeke)
                                          ! ! Old Ice Z0 = 0.0500 m (Bruce)
                                          ! !              0.0500 m (Smeets)
                                          ! !              0.1200 m (Broeke)
  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
  ! +                             ! for Suspension
  ! +                             ! Pommeroy, Gray and Landine 1993,
  ! +                             ! J. Hydrology, 144(8) p.169
  data      nit   / 5   /   ! us(is0,uth) recursivity: Nb Iterations
  !c#AE data      qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p
  data      qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray
  !c#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


  ! BEGIN.main.
  ! ===========================




  ! "Soil" Humidity of Water Bodies
  ! ===============================

  DO ikl=1,knonv

      ist    =      isotSV(ikl)                       ! Soil Type
      ist__s =  min(ist, 1)                           ! 1 => Soil
      ist__w =  1 - ist__s                            ! 1 => Water Body
    DO isl=-nsol,0
      eta_SV(ikl,isl) = eta_SV(ikl,isl) * ist__s & ! Soil
            + etadSV(ist)     * ist__w      ! Water Body
    END DO


  ! Vertical Discretization Factor
  ! ==============================

      LSdzsv(ikl)     =                   ist__s & ! Soil
            + OcndSV          * ist__w      ! Water Body
  END DO





  IF (SnoMod)                            THEN


  ! +--Aeolian erosion and Blowing Snow
  ! +==================================



    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
  ! +                       **********
                 call SISVAT_BSn
      endif
     else
                 call SISVAT_BSn
  ! +                       **********
    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

  ! +--Computation of threshold friction velocity for snow erosion
  ! ---------------------------------------------------------------

    rCd10n =  1. / 26.5 ! Vt / u*t = 26.5
                 ! ! Budd et al. 1965, Antarct. Res. Series Fig.13
                 ! ! ratio developped during assumed neutral conditions


  ! +--Snow Properties
  ! +  ~~~~~~~~~~~~~~~

    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


  ! +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SaltM1   = -0.750e-2 * G1snSV(ikl,isn) &
            -0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case
  ! +     CAUTION:  Guyomarc'h & Merindol Dendricity Sign is +
  ! +     ^^^^^^^^                    MAR Dendricity Sign is -
      SaltM2   = -0.833d-2 * G1snSV(ikl,isn) &
            -0.583d-2 * G2snSV(ikl,isn)+ 0.833d00 !non-dendritic case

    ! 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)
  !c#AE   FacRho   = 1.25 - 0.0042 * ro__SV(ikl,isn)
  !c#AE   SaltMo   = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow
      MIN_Mo   =  0.
    ! SaltMo   =  max(SaltMo,MIN_Mo)
    ! SaltMo   =  SaltOK   * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx)
  ! #TUNE SaltMo   =  SaltOK   * SaltMo - (1.-SaltOK) *     0.9500
      SaltMo   =  max(SaltMo,epsi-unun)

  ! +--Influence of Density on Threshold Shear Stress
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      Por_BS =  1. - 300. / ro_Ice
      ShearS = Por_BS / (1.-Por_BS)
  ! +...         SheaBS =  Arg(sqrt(shear = max shear stress in snow)):
  ! +            shear  =  3.420d00 * exp(-(Por_BS      +Por_BS)
  ! +  .                                  /(unun        -Por_BS))
  ! +            SheaBS :  see de Montmollin         (1978),
  ! +                      These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124

  ! +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      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


  ! +--Threshold Friction Velocity
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
      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)
  ! +     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
  ! +...  Salt_us   :  Extension of  Guyomarc'h & Merindol 1998 with
  ! +...              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

  ! + ------------------------------------------------------
  ! +--Buffer Layer
  ! +  -----------------------------------------------------

      DO ikl=1,knonv
  !  BufsSV(ikl) [mm w.e.] i.e, i.e., [kg/m2]
        d_Bufs      =  max(dsn_SV(ikl) *dt__SV,0.)  !
        dsn_SV(ikl) =      0.                       !
        Bufs_N      =      BufsSV(ikl) +d_Bufs      !


  ! +--Snow Density
  ! +  ^^^^^^^^^^^^
        Polair      =      zero
  ! #NP       Polair      =  max(zero,                    !
  ! #NP.                         sign(unun,TaPole         !
  ! #NP.                                  -TaT_SV(ikl)))  !
        Polair      =  max(zero, & !
              sign(unun,TaPole & !
              -TaT_SV(ikl)))  !
        Buf_ro      =  max( rosMin, & ! Fallen Snow Density
              roSn_1+roSn_2*     (TaT_SV(ikl)-TfSnow) & ! [kg/m3]
              +roSn_3*sqrt( VV__SV(ikl)))           ! Pahaut    (CEN), Etienne: use wind speed at first model level instead of 10m wind
  ! #NP       BufPro      =  max( rosMin,                 ! Fallen Snow Density
  ! #NP.         104. *sqrt( max( VV10SV(ikl)-6.0,0.0)))  ! Kotlyakov (1961)

       ! C.Agosta option for snow density, same as for BS i.e.
       ! is_ok_density_kotlyakov=.false.
  ! #BS       density_kotlyakov = .false.  !C.Amory BS 2018
  ! + ...     Fallen Snow Density, Adapted for Antarctica
        if (is_ok_density_kotlyakov) then
            tt_tmp = TaT_SV(ikl)-TfSnow
            ! !vv_tmp = VV10SV(ikl)
            vv_tmp=VV__SV(ikl) ! Etienne: use wind speed at first model level instead of 10m wind
  ! + ...         [ A compromise between
  ! + ...           Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ]
            if (tt_tmp.ge.-10) then
              BufPro   =  max( rosMin, &
                    104. *sqrt( max( vv_tmp-6.0,0.0))) ! Kotlyakov (1961)
            else
              vv_virt = (tt_c*vv_tmp+vv_c*(tt_tmp+10)) &
                    /(tt_c+tt_tmp+10)
              BufPro  = 104. *sqrt( max( vv_virt-6.0,0.0))
            endif
        else
  ! + ...         [ density derived from observations of the first 50cm of
  ! + ...           snow - cf. Rajashree Datta - and multiplied by 0.8 ]
  ! + ...           C. Agosta, 2016-09
  !c #SD           BufPro = 149.2 + 6.84*VV10SV(ikl) + 0.48*Tsrfsv(ikl)
  !c #SD           BufPro = 125 + 14*VV10SV(ikl) + 0.6*Tsrfsv(ikl) !MAJ CK and CAm
             ! BufPro = 200 + 21 * VV10SV(ikl)!CK 29/07/19
             BufPro = 200 + 21 * VV__SV(ikl)!Etienne: use wind speed at first model level instead of 10m wind
        endif

        Bros_N      = (1. - Polair) *   Buf_ro & ! Temperate Snow
              + Polair  *   BufPro      ! Polar     Snow

        Bros_N = max( 20.,max(rosMin,  Bros_N))
        Bros_N = min(400.,min(rosMax-1,Bros_N)) ! for dz_min in SISVAT_zSn


  !    Density of deposited blown snow
  !    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

     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


  !    Time averaged Density of deposited blown Snow
  !    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        BrosSV(ikl) =(Bros_N     *      d_Bufs & !
              +BrosSV(ikl)*      BufsSV(ikl)) & !
              /         max(epsi,Bufs_N)     !


  ! +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
  ! +     ^^^^^^^^^^^^^^^^^^^^^^^
        Buf_G1      =  max(-G1_dSV, & ! Temperate Snow
              min(Dendr1*VV__SV(ikl)-Dendr2, & !     Dendricity
              Dendr3                   ))    !
        Buf_G2      =  min( Spher4, & ! Temperate Snow
              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  *   buf_sph_pol ! Polar Snow
        Buf_G2      = (1. - Polair) *   Buf_G2 & ! Temperate 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
  !      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

       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
  !      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
       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        )) !
  ! +...      REMARQUE: le  type moyen (dendritique ou non) depend
  ! +         ^^^^^^^^  de la  comparaison avec le diametre optique
  ! +                   d'une neige recente de   dendricite nulle
  ! +...      REMARK:   the mean type  (dendritic   or not) depends
  ! +         ^^^^^^    on the comparaison with the optical diameter
  ! +                   of a recent snow    having zero dendricity

       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



  ! S.1. Meme  Type  de Neige  / same Grain Type
  !      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        SameOK  =  max(zero, &
              sign(unun,    Buf_G1 *BG1sSV(ikl) &
              - eps_21                    ))
        G1same  = (d_Bufs*Buf_G1+BufsSV(ikl)*BG1sSV(ikl)) &
              /max(epsi,Bufs_N)
        G2same  = (d_Bufs*Buf_G2+BufsSV(ikl)*BG2sSV(ikl)) &
              /max(epsi,Bufs_N)

  ! S.2. Types differents / differents Types
  !      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

        typ__1  =  max(zero,sign(unun,epsi-Buf_G1))   ! =1.=> Dendritic
        zroNEW  =(    typ__1  *d_Bufs & ! fract.Dendr.Lay.
              + (1.-typ__1) *BufsSV(ikl)) & !
              /max(epsi,Bufs_N)                !
        G1_NEW  =     typ__1  *Buf_G1 & ! G1 of Dendr.Lay.
              + (1.-typ__1) *BG1sSV(ikl)            !
        G2_NEW  =     typ__1  *Buf_G2 & ! G2 of Dendr.Lay.
              + (1.-typ__1) *BG2sSV(ikl)            !
        zroOLD  =((1.-typ__1) *d_Bufs & ! fract.Spher.Lay.
              +     typ__1  *BufsSV(ikl)) & !
              /max(epsi,Bufs_N)                !
        G1_OLD  = (1.-typ__1) *Buf_G1 & ! G1 of Spher.Lay.
              +     typ__1  *BG1sSV(ikl)            !
        G2_OLD  = (1.-typ__1) *Buf_G2 & ! G2 of Spher.Lay.
              +     typ__1  *BG2sSV(ikl)            !
        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        )) !
  ! +...      REMARQUE: le  type moyen (dendritique ou non) depend
  ! +         ^^^^^^^^  de la  comparaison avec le diametre optique
  ! +                   d'une neige recente de   dendricite nulle
  ! +...      REMARK:   the mean type  (dendritic   or not) depends
  ! +         ^^^^^^    on the comparaison with the optical diameter
  ! +                   of a recent snow    having zero dendricity

        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

        BG1sSV(ikl) =                       G1 & !
              *       Bufs_N/max(epsi,Bufs_N) !
        BG2sSV(ikl) =                       G2 & !
              *       Bufs_N/max(epsi,Bufs_N) !


  ! +--Update of Buffer Layer Content & Decision about creating a new snow layer
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        BufsSV(ikl) =       Bufs_N                  !     [mm w.e.]
        NLaysv(ikl) = min(unun, & !
              max(zero, & ! Allows to create
              sign(unun,BufsSV(ikl) & ! a new snow Layer
              -SMndSV     )) & ! if Buffer > SMndSV
              *max(zero, & ! Except if * Erosion
              sign(unun,0.50 & ! dominates
              -dsnbSV(ikl))) & !
              +max(zero, & ! Allows to create
              sign(unun,BufsSV(ikl) & ! a new snow Layer
              -SMndSV*3.00)))  ! is Buffer > SMndSV*3
        Bdzssv(ikl) = 1.e-3*BufsSV(ikl)*ro_Wat & ! [mm w.e.] -> [m w.e.]
              /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m]


      END DO



  ! 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
  ! +          **********
     call SISVAT_zSn
  ! +          **********
   endif
  else
  ! +          **********
    call SISVAT_zSn
  ! +          **********
  endif

  ! +          **********
  ! #ve   call SISVAT_wEq('_zSn  ',0)
  ! +          **********

  ! Add a new Snow Layer
  ! ====================

      DO ikl=1,knonv

        isnoSV(ikl)     = isnoSV(ikl)         +NLaysv(ikl)
        isn             = isnoSV(ikl)
        dzsnSV(ikl,isn) = dzsnSV(ikl,isn) * (1-NLaysv(ikl)) &
              + Bdzssv(ikl)     *    NLaysv(ikl)
        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)
        eta_SV(ikl,isn) = eta_SV(ikl,isn) * (1-NLaysv(ikl))   ! + 0.
        agsnSV(ikl,isn) = agsnSV(ikl,isn) * (1-NLaysv(ikl))   ! + 0.
        G1snSV(ikl,isn) = G1snSV(ikl,isn) * (1-NLaysv(ikl)) &
              + BG1ssv(ikl)     *    NLaysv(ikl)
        G2snSV(ikl,isn) = G2snSV(ikl,isn) * (1-NLaysv(ikl)) &
              + BG2ssv(ikl)     *    NLaysv(ikl)
        istoSV(ikl,isn) = istoSV(ikl,isn) * (1-NLaysv(ikl)) &
              + max(zer0,sign(un_1,TaT_SV(ikl) &
              -Tf_Sno-eps_21)) *    istdSV(2) &
              *    NLaysv(ikl)
        BufsSV(ikl)     = BufsSV(ikl)     * (1-NLaysv(ikl))
        NLaysv(ikl)     = 0


      END DO


  ! Snow Pack Thickness
  ! -------------------

    DO ikl=1,knonv
        z_snsv(ikl)     = 0.0
    END DO
    DO   isn=1,nsno
      DO ikl=1,knonv
        z_snsv(ikl)     = z_snsv(ikl) + dzsnSV(ikl,isn)
        zzsnsv(ikl,isn) = z_snsv(ikl)
      END DO
    END DO



  END IF  ! SnoMod



  ! Soil Albedo: Soil Humidity Correction
  ! ==========================================

      ! REFERENCE: McCumber and Pielke (1981), Pielke (1984)
      ! ^^^^^^^^^
      DO ikl=1,knonv
        albssv(ikl) = &
              alb0SV(ikl) *(1.0-min(half,eta_SV(       ikl,0) &
              /etadSV(isotSV(ikl))))
      ! REMARK:    Albedo of Water Surfaces (isotSV=0):
      ! ^^^^^^     alb0SV := 2  X  effective value, while
      !            eta_SV :=          etadSV
      END DO


  ! Snow Pack Optical Properties
  ! ============================

  IF (SnoMod)                                                 THEN

         ! ******
    call SnOptP(jjtime)
         ! ******

  ELSE
    DO ikl=1,knonv
      sEX_sv(ikl,1) = 1.0
      sEX_sv(ikl,0) = 0.0
      albisv(ikl)   = albssv(ikl)
    END DO
  END IF



  ! Soil optical properties
  ! =============================
  !Etienne: as in inlandis we do not call vgopt, we need to define
  !the albedo alb_SV and to calculate the
  !absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv


  DO ikl=1,klonv

        e_pRad = 2.5   *  coszSV(ikl)       ! exponential argument,
                                            ! ! V/nIR radiation partitioning,
                                            ! ! DR97, 2, eqn (2.53) & (2.54)
        e1pRad = 1.-exp(-e_pRad)            ! exponential, V/nIR Rad. Part.
        exdRad= 1.

  ! Visible Part of the Solar Radiation Spectrum (V,   0.4--0.7mi.m)
        A_Rad0 =      0.25 + 0.697 * e1pRad ! Absorbed    Vis. Radiation
        absg_V = (1.-albisv(ikl))*(A_Rad0*exdRad)  !

  ! Near-IR Part of the Solar Radiation Spectrum (nIR, 0.7--2.8mi.m)

        A_Rad0 =      0.80 + 0.185 * e1pRad ! Absorbed    nIR. Radiation
        absgnI = (1.-albisv(ikl))*(A_Rad0*exdRad)  !

        SoSosv(ikl) = (absg_V+absgnI)*0.5d0

        alb_SV(ikl) = albisv(ikl)

  END DO

         ! **********
  ! #ve   call SISVAT_wEq('SnOptP',0)
         ! **********


  ! Surface Emissivity (Etienne: simplified calculation for landice)
  ! =============================================================
  !
   DO ikl=1,knonv
        LSnMsk     =     min( 1,isnoSV(ikl))
        Eso_sv(ikl)=  EmiSol*(1-LSnMsk)+EmiSno*LSnMsk  ! Sol+Sno Emissivity
        emi_SV(ikl)= EmiSol*(1-LSnMsk) + EmiSno*LSnMsk
    END DO




  !  Upward IR (INPUT, from previous time step)
  ! ===================================================================

    DO ikl=1,knonv
  ! #e1     Enrsvd(ikl) =    - IRs_SV(ikl)
       IRupsv(ikl) =      IRs_SV(ikl)
    END DO


  ! Turbulence
  ! ==========

  ! Latent Heat of Vaporization/Sublimation
  ! ---------------------------------------

    DO ikl=1,knonv
      SnoWat      =                     min(isnoSV(ikl),0)
      Lx_H2O(ikl) = &
            (1.-SnoWat) * LhvH2O &
            +     SnoWat  *(LhsH2O * (1.-eta_SV(ikl,isnoSV(ikl))) &
            +LhvH2O *     eta_SV(ikl,isnoSV(ikl)) )
    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



  ! Soil   Energy Balance
  ! =====================


  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

  end if


         ! **********
  ! #ve   call SISVAT_wEq('_TSo  ',0)
         ! **********



  ! Soil Water     Potential
  ! ------------------------

  DO   isl=-nsol,0
    DO ikl=1,knonv
      ist             =     isotSV(ikl)        ! Soil Type
      psi_sv(ikl,isl) =     psidSV(ist) & ! DR97, Eqn.(3.34)
            *(etadSV(ist) /max(eps6,eta_SV(ikl,isl))) & !
            **bCHdSV(ist)                              !


  ! Soil Hydraulic Conductivity
  ! ---------------------------

      Khydsv(ikl,isl) =    s2__SV(ist) & ! DR97, Eqn.(3.35)
            *(eta_SV(ikl,isl)**(2.*bCHdSV(ist)+3.))    !
    END DO
  END DO


  ! Melting / Refreezing in the Snow Pack
  ! =====================================

  IF (SnoMod)                                                 THEN

         ! **********
    call SISVAT_qSn
         ! **********

         ! **********
  ! #ve   call SISVAT_wEq('_qSn  ',0)
         ! **********


  ! Snow Pack Thickness
  ! -------------------

      DO ikl=1,knonv
        z_snsv(ikl)     = 0.0
      END DO
    DO   isn=1,nsno
      DO ikl=1,knonv
        z_snsv(ikl)     = z_snsv(ikl) + dzsnSV(ikl,isn)
        zzsnsv(ikl,isn) = z_snsv(ikl)
      END DO
    END DO


  ! Energy in Excess is added to the first Soil Layer
  ! -------------------------------------------------
    DO ikl=1,knonv
        z_snsv(ikl)   = max(zer0, &
              sign(un_1,eps6-z_snsv(ikl)))
        TsisSV(ikl,0) = TsisSV(ikl,0)    + EExcsv(ikl) &
              /(rocsSV(isotSV(ikl)) &
              +rcwdSV*eta_SV(ikl,0))
        EExcsv(ikl)   = 0.
    END DO


  END IF


  ! Soil   Water  Balance
  ! =====================

         ! **********
    call SISVAT_qSo
  ! #m0.                 (Wats_0,Wats_1,Wats_d)
         ! **********


  ! Surface Fluxes
  ! =====================

    DO ikl=1,knonv
     IRdwsv(ikl)=IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
      ! IRdwsv(ikl)=tau_sv(ikl) *IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
  !    .          +(1.0-tau_sv(ikl))*IRd_SV(ikl)*Evg_sv(ikl) !  ! Etienne, remove vegetation component
      IRupsv(ikl) =      IRupsv(ikl)                   ! Upward   IR
      IRu_SV(ikl) =     -IRupsv(ikl) & ! Upward   IR
            +IRd_SV(ikl) & ! (effective)
            -IRdwsv(ikl)                   ! (positive)

      TBr_sv(ikl) =sqrt(sqrt(IRu_SV(ikl)/StefBo))      ! Brightness
                                                       ! ! Temperature
      uts_SV(ikl) =     (HSv_sv(ikl) +HSs_sv(ikl)) & ! u*T*
            /(rhT_SV(ikl) *cp)          !
      uqs_SV(ikl) =     (HLv_sv(ikl) +HLs_sv(ikl)) & ! u*q*
            /(rhT_SV(ikl) *LhvH2O)          !
      LMO_SV(ikl) = TaT_SV(ikl)*(us__SV(ikl)**3) &
            /gravit/uts_SV(ikl)/vonKrm      ! MO length


  ! Surface Temperature
  ! ^^^^^^^^^^^^^^^^^^^^

      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:

     if (isnoSV(ikl) >=2) then
       zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
       zm2=-(dzsnSV(ikl,isnoSV(ikl)) + dzsnSV(ikl,isnoSV(ikl)-1)/2.)
     else if (isnoSV(ikl) .EQ. 1) then
       zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
       zm2=-(dzsnSV(ikl,isnoSV(ikl))+dz_dSV(0)/2.)
     else
       zm1=-dz_dSV(0)/2.
       zm2=-(dz_dSV(0)+dz_dSV(-1)/2.)

     end if

       coefslope=(TsisSV(ikl,isnoSV(ikl))-TsisSV(ikl,isnoSV(ikl)-1)) &
             /(zm1-zm2)
       Tsrfsv(ikl)=TsisSV(ikl,isnoSV(ikl))+coefslope*(0. - zm1)


     ELSE !(default)

       Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))

     END IF


     END DO

  ! Snow Pack Properties (sphericity, dendricity, size)
  ! ===================================================

  IF (SnoMod)                                                 THEN

  if (discret_xf .AND. klonv.eq.1) then
  if(isnoSV(1).ge.1) then
  ! +          **********
  call SISVAT_GSn
  ! +          **********
  endif
  else
  ! +          **********
    call SISVAT_GSn
  ! +          **********
  endif


  END IF


  ! Roughness Length for next time step
  !====================================

  ! Note that in INLANDSIS, we treat only ice covered surfaces so calculation
  ! of z0 is much simpler (no subgrid fraction of ocean or land)
  ! old calculations are commented below


  ! +--Roughness Length for Momentum
  ! +  -----------------------------

  ! ETIENNE WARNING: changes have been made wrt original SISVAT

  ! +--Land+Sea-Ice / Ice-free Sea Mask
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    DO ikl=1,knonv
      IcIndx(ikl) = 0
    ENDDO
    DO isn=1,nsno
    DO ikl=1,knonv

      IcIndx(ikl) = max(IcIndx(ikl), &
            isn*max(0, &
            sign(1, &
            int(ro__SV(ikl,isn)-900.))))
    ENDDO
    ENDDO

    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)


  ! +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
      Z0m_nu =       5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))

  ! +--Z0 Saltat.Regime over Snow (Gallee  et al., 2001, BLM 99 (19) p.11)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^

      u2star =       us__SV(ikl) *us__SV(ikl)
      Z0mBSn =       u2star      *0.536e-3   -  61.8e-6
      Z0mBSn =   max(Z0mBS0      ,Z0mBSn)

  ! +--Z0 Smooth + Saltat. Regime
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
      Z0enSV(ikl) =  Z0m_nu &
            +  Z0mBSn


  ! 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

  ! +--Rough   Snow Surface Roughness Length (Variable Sastrugi Height)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      A_Fact      =  1.0000        ! Andreas et al., 2004, p.4
                                   ! ! ams.confex.com/ams/pdfpapers/68601.pdf

  ! Parameterization of z0 dependance on Temperature (C. Amory, 2017)
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013


      coefa = 0.1658 !0.1862 !Ant
      coefb = -50.3869 !-55.7718 !Ant
      ta1 = 253.15 !255. Ant
      ta2 = 273.15
      ta3 = 273.15+3
      z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
      z02 = exp(coefa*ta2 + coefb) !~6  !~7 mm
      z03 = z01
      coefc = log(z03/z02)/(ta3-ta2)
      coefd = log(z03)-coefc*ta3

      if (TaT_SV(ikl) .lt. ta1) then
        Z0_obs = z01
      else if (TaT_SV(ikl).ge.ta1 .and. TaT_SV(ikl).lt.ta2) then
        Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
      else if (TaT_SV(ikl).ge.ta2 .and. TaT_SV(ikl).lt.ta3) then
        ! ! if st > 0, melting induce smooth surface
        Z0_obs = exp(coefc*TaT_SV(ikl) + coefd)
      else
        Z0_obs = z03
      endif

      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
  ! #SZ     Z0Sa_N =                   (us__SV(ikl) -0.2)*param   ! 0.0001=TUNING
  ! #SZ.           * max(zero,sign(unun,TfSnow-eps9
  ! #SZ.                               -TsisSV(ikl , isnoSV(ikl))))
  !!#SZ     Z0SaSi = max(zero,sign(unun,Z0Sa_N                  ))! 1 if erosion
  ! #SZ     Z0SaSi = max(zero,sign(unun,zero  -eps9 -uss_SV(ikl)))!
  ! #SZ     Z0Sa_N = max(zero,          Z0Sa_N)
  ! #SZ     Z0SaSV(ikl) =
  ! #SZ.             max(Z0SaSV(ikl)   ,Z0SaSV(ikl)
  ! #SZ.               + Z0SaSi*(Z0Sa_N-Z0SaSV(ikl))*exp(-dt__SV/43200.))
  ! #SZ.               -            min(dz0_SV(ikl) ,     Z0SaSV(ikl))

  ! #SZ     A_Fact      =               Z0SaSV(ikl) *  5.0/0.15   ! A=5 if h~10cm
  ! +...    CAUTION: The influence of the sastrugi direction is not yet included

  ! #SZ     Z0m_Sn =                    Z0SaSV(ikl)               !
  ! #SZ.                              - Z0m_nu                    !

  ! +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46)  p.222)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
  ! #ZN     sqrrZ0 =       usthSV(ikl)/max( us__SV(ikl),0.001)
  ! #ZN     sqrrZ0 =                   min( sqrrZ0     ,0.999)
  ! #ZN     Z0mBSn =       0.55 *0.55 *exp(-sqrrZ0     *sqrrZ0)
  ! #ZN.                  *us__SV(ikl)*     us__SV(ikl)*grvinv*0.5

  ! +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46)  p.222)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
  ! #ZN     Z0enSV(ikl) = (Z0m_nu     **    sqrrZ0 )
  ! #ZN.                * (Z0mBSn     **(1.-sqrrZ0))
  ! #ZN     Z0enSV(ikl) =  max(Z0enSV(ikl), Z0m_nu)


  ! +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^  ams.confex.com/ams/pdfpapers/68601.pdf)
  ! #ZA     Z0m_nu = 0.135*akmol  / max(us__SV(ikl) , epsi)

  ! +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^  ams.confex.com/ams/pdfpapers/68601.pdf)
  ! #ZA     Z0mBSn = 0.035*u2star      *grvinv

  ! +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
  !    (      used by Erosion)     ams.confex.com/ams/pdfpapers/68601.pdf)
  !    ^^^^^^^^^^^^^^^^^^^^^^^^^^
  ! #ZA     Z0enSV(ikl) =  Z0m_nu
  ! #ZA.                +  Z0mBSn

  ! +--Z0 Rough  Regime over Snow (Andreas etAl., 2004
  ! +  (.NOT. used by Erosion)     ams.confex.com/ams/pdfpapers/68601.pdf)
  !    ^^^^^^^^^^^^^^^^^^^^^^^^^^
  !!#ZA     u2star =      (us__SV(ikl) -0.1800)     / 0.1
  !!#ZA     Z0m_Sn =A_Fact*Z0mBSn *exp(-u2star*u2star)
  ! #ZA     Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
  ! #ZA.            *exp(-0.4/sqrt(.00275+.00001*max(0.,VVs_SV(ikl)-5.)))
  ! #ZA     Z0m_Sn =           DDs_SV(ikl)* Z0m_90 / 45.
  ! #ZA.         - DDs_SV(ikl)*DDs_SV(ikl)* Z0m_90 /(90.*90.)




  ! +--Z0  (Erosion)    over Snow (instantaneous)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
      Z0e_SV(ikl) =  Z0enSV(ikl)

  ! +--Momentum  Roughness Length (Etienne: changes wrt original SISVAT)
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      Z0mnSV(ikl) =  Z0m_nu *(1-SnoMsk) & ! Ice z0
            + (Z0m_Sn)*SnoMsk                         ! Snow Sastrugi Form and Snow Erosion


  ! +--GIS  Roughness Length
  ! +  ^^^^^^^^^^^^^^^^^^^^^
  ! #GL     Z0mnSV(ikl) =
  ! #GL.      (1-LSmask(ikl)) *     Z0mnSV(ikl)
  ! #GL.    +    LSmask(ikl)  * max(Z0mnSV(ikl),max(Z0_GIM,
  ! #GL.                                            Z0_GIM+
  ! #GL.      (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.)   !
  ! #GL.                     /(920.00                 -600.))) !

  ! +--Mom. Roughness Length, Instantaneous
  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      Z0m_SV(ikl) =  Z0mnSV(ikl)                         ! Z0mnSV  instant.


  ! +--Roughness Length for Scalars
  ! +  ----------------------------

      Z0hnSV(ikl) =     Z0mnSV(ikl)/  7.4

      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)


  ! #MT     Z0m_SV(ikl) = max(2.0e-6     ,Z0m_SV(ikl)) ! Min Z0_m (Garrat Scheme)
       ! Z0m_SV(ikl) = min(Z0m_SV(ikl),za__SV(ikl)*0.3333)


   END DO


   return
end subroutine inlandsis











