      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   !
c #sw real      PorVol,rWater                 !
c #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
c #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
c #SZ real      Z0Sa_N                        ! Regime   Snow Roughness Length
c #SZ real      Z0SaSi                        ! 1.IF Rgm Snow Roughness Length 
c #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*3    qsalt_param              ! Switch for saltation flux param.
      character*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)
c #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
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


! 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

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



        DO ikl=1,knonv
            usthSV(ikl) =                     1.0e+2
        END DO


        IF (BloMod) THEN
 
        if (klonv==1) then
          if(isnoSV(1)>=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
C +                       **********
                     CALL SISVAT_BSn
          endif
         else
                     CALL SISVAT_BSn
C +                       **********
        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 == "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 == "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 == "pom") then
          qSalSV(ikl) = (us__SV(ikl)**2 - usthSV(ikl)**2) *signus
     .               / (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25)
        endif
 
        if (qsalt_param == "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 +  -----------------------------------------------------
 
          DO ikl=1,knonv
c  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      !
 
 
C +--Snow Density
C +  ^^^^^^^^^^^^
            Polair      =      zero
c #NP       Polair      =  max(zero,                    !
c #NP.                         sign(unun,TaPole         !
c #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
c #NP       BufPro      =  max( rosMin,                 ! Fallen Snow Density
c #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.
c #BS       density_kotlyakov = .FALSE.  !C.Amory BS 2018
C + ...     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
C + ...         [ A compromise between
C + ...           Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ]
                if (tt_tmp>=-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
C + ...         [ density derived from observations of the first 50cm of
C + ...           snow - cf. Rajashree Datta - and multiplied by 0.8 ]
C + ...           C. Agosta, 2016-09
cc #SD           BufPro = 149.2 + 6.84*VV10SV(ikl) + 0.48*Tsrfsv(ikl)
cc #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)     !
 
 
C +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
C +     ^^^^^^^^^^^^^^^^^^^^^^^
            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        )) !
C +...      REMARQUE: le  type moyen (dendritique ou non) depend
C +         ^^^^^^^^  de la  comparaison avec le diametre optique
C +                   d'une neige recente de   dendricite nulle
C +...      REMARK:   the mean type  (dendritic   or not) depends
C +         ^^^^^^    on the comparaison with the optical diameter
C +                   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        )) !
C +...      REMARQUE: le  type moyen (dendritique ou non) depend
C +         ^^^^^^^^  de la  comparaison avec le diametre optique
C +                   d'une neige recente de   dendricite nulle
C +...      REMARK:   the mean type  (dendritic   or not) depends
C +         ^^^^^^    on the comparaison with the optical diameter
C +                   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) !
 

C +--Update of Buffer Layer Content & Decision about creating a new snow layer
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
            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==1) then

       if(isnoSV(1)>=1.or.NLaysv(1)>=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
! ====================

          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 == 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 == 0) THEN

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

          ELSE IF (iflag_tsurf_inlandsis > 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) == 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==1) then
      if(isnoSV(1)>=1) then
C +          **********
      CALL SISVAT_GSn
C +          **********
      endif
      else
C +          **********
        CALL SISVAT_GSn
C +          **********
      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


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,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)


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

        
! Calculation of snow roughness length
!=====================================
          IF (iflag_z0m_snow == 0) THEN

          Z0m_Sn=prescribed_z0m_snow

          ELSE IF (iflag_z0m_snow == 1) THEN

          Z0m_Sn=Z0enSV(ikl) 

          ELSE IF (iflag_z0m_snow == 2) THEN

C +--Rough   Snow Surface Roughness Length (Variable Sastrugi Height)
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          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) < ta1) then
            Z0_obs = z01
          else if (TaT_SV(ikl)>=ta1 .and. TaT_SV(ikl)<ta2) then
            Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
          else if (TaT_SV(ikl)>=ta2 .and. TaT_SV(ikl)<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 
c #SZ     Z0Sa_N =                   (us__SV(ikl) -0.2)*param   ! 0.0001=TUNING
c #SZ.           * max(zero,sign(unun,TfSnow-eps9
c #SZ.                               -TsisSV(ikl , isnoSV(ikl))))
!!#SZ     Z0SaSi = max(zero,sign(unun,Z0Sa_N                  ))! 1 if erosion
c #SZ     Z0SaSi = max(zero,sign(unun,zero  -eps9 -uss_SV(ikl)))!
c #SZ     Z0Sa_N = max(zero,          Z0Sa_N)
c #SZ     Z0SaSV(ikl) =
c #SZ.             max(Z0SaSV(ikl)   ,Z0SaSV(ikl)
c #SZ.               + Z0SaSi*(Z0Sa_N-Z0SaSV(ikl))*exp(-dt__SV/43200.))
c #SZ.               -            min(dz0_SV(ikl) ,     Z0SaSV(ikl))
 
c #SZ     A_Fact      =               Z0SaSV(ikl) *  5.0/0.15   ! A=5 if h~10cm
C +...    CAUTION: The influence of the sastrugi direction is not yet included
 
c #SZ     Z0m_Sn =                    Z0SaSV(ikl)               !
c #SZ.                              - Z0m_nu                    !
 
C +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46)  p.222)
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
c #ZN     sqrrZ0 =       usthSV(ikl)/max( us__SV(ikl),0.001)
c #ZN     sqrrZ0 =                   min( sqrrZ0     ,0.999)
c #ZN     Z0mBSn =       0.55 *0.55 *exp(-sqrrZ0     *sqrrZ0)
c #ZN.                  *us__SV(ikl)*     us__SV(ikl)*grvinv*0.5
 
C +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46)  p.222)
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
c #ZN     Z0enSV(ikl) = (Z0m_nu     **    sqrrZ0 )
c #ZN.                * (Z0mBSn     **(1.-sqrrZ0))
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)
c #ZA     Z0m_nu = 0.135*akmol  / max(us__SV(ikl) , epsi)
 
C +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^  ams.confex.com/ams/pdfpapers/68601.pdf)
c #ZA     Z0mBSn = 0.035*u2star      *grvinv
 
C +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
!    (      used by Erosion)     ams.confex.com/ams/pdfpapers/68601.pdf)
!    ^^^^^^^^^^^^^^^^^^^^^^^^^^
c #ZA     Z0enSV(ikl) =  Z0m_nu
c #ZA.                +  Z0mBSn
 
C +--Z0 Rough  Regime over Snow (Andreas etAl., 2004
C +  (.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)
c #ZA     Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
c #ZA.            *exp(-0.4/sqrt(.00275+.00001*max(0.,VVs_SV(ikl)-5.)))
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)
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
          Z0e_SV(ikl) =  Z0enSV(ikl)
 
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
 

C +--GIS  Roughness Length
C +  ^^^^^^^^^^^^^^^^^^^^^
c #GL     Z0mnSV(ikl) =
c #GL.      (1-LSmask(ikl)) *     Z0mnSV(ikl)
c #GL.    +    LSmask(ikl)  * max(Z0mnSV(ikl),max(Z0_GIM,
c #GL.                                            Z0_GIM+
c #GL.      (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.)   !
c #GL.                     /(920.00                 -600.))) !
 
C +--Mom. Roughness Length, Instantaneous
C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          Z0m_SV(ikl) =  Z0mnSV(ikl)                         ! Z0mnSV  instant.
 
 
C +--Roughness Length for Scalars
C +  ----------------------------
 
          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)
 

c #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










 
