MODULE surf_inlandsis_mod

    IMPLICIT NONE

CONTAINS


SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, &
            rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
            precip_rain, precip_snow, &
            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
            rugos, snow_cont_air, alb_soil, alt, slope, cloudf, &
            radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
            runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, &
            tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf)

        ! |                                                                        |
        ! |   SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow|
        ! |                              (INLANDSIS)                               |
        ! |   SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)         |
        ! |   surface scheme of the Modele Atmospherique Regional (MAR)            |
        ! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
        ! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
        ! |   Updated by Etienne Vignon, Cecile Agosta                             |
        ! |                                                                        |
        ! +------------------------------------------------------------------------+
        ! |
        ! |   In the current setup, SISVAT is used only to model the land ice      |
        ! |   part of the surface; hence it is called with the compressed variables|
        ! |   from pbl_surface, and only by the surf_landice routine.              |
        ! |                                                                        |
        ! |   In this interface it is assumed that the partitioning of the soil,   |
        ! |   and hence the number of grid points is constant during a simulation, |
        ! |   hence eg. snow properties remain stored in the global SISVAT         |
        ! |   variables between the calls and don't need to be handed over as      |
        ! |   arguments. When the partitioning is supposed to change, make sure to |
        ! |   update the variables.                                                |
        ! |                                                                        |
        ! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV ...)                    |
        ! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |
        ! |                                                                        |
        ! +------------------------------------------------------------------------+

        USE dimphy
        USE VAR_SV
        USE VARdSV
        USE VARxSV
        USE VARySV
        USE VARtSV
        USE VARphy
        USE surface_data, ONLY: iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor

        IMPLICIT NONE

        ! +--INTERFACE Variables
        ! +  ===================
        !    include  "dimsoil.h"

        ! +--Global Variables
        ! +  ================
        ! Input Variables for SISVAT
        INTEGER, INTENT(IN) :: knon
        INTEGER, INTENT(IN) :: itime
        REAL, INTENT(IN) :: dtime
        LOGICAL, INTENT(IN) :: debut     ! true if first step
        LOGICAL, INTENT(IN) :: lafin     ! true if last step

        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i     ! Index Decompression
        REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
        REAL, DIMENSION(klon), INTENT(IN) :: rmu0      ! cos sol. zenith angle
        REAL, DIMENSION(klon), INTENT(IN) :: swdown    !
        REAL, DIMENSION(klon), INTENT(IN) :: lwdown    !
        REAL, DIMENSION(klon), INTENT(IN) :: albedo_old
        REAL, DIMENSION(klon), INTENT(IN) :: pexner    ! Exner potential
        REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
        REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo
        REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay
        REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf
        REAL, DIMENSION(klon), INTENT(IN) :: rugos
        REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air
        REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope
        REAL, DIMENSION(klon), INTENT(IN) :: alt       ! surface elevation
        REAL, DIMENSION(klon), INTENT(IN) :: cloudf
        REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
        REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
        REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh
        REAL, DIMENSION(klon), INTENT(IN) :: ustar   ! friction velocity

        ! Variables exchanged between LMDZ and SISVAT
        REAL, DIMENSION(klon), INTENT(IN) :: radsol    ! Surface absorbed rad.
        REAL, DIMENSION(klon), INTENT(INOUT) :: snow      ! Tot snow mass [kg/m2]
        REAL, DIMENSION(klon), INTENT(INOUT) :: zfra      ! snwo surface fraction [0-1]
        REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
        REAL, DIMENSION(klon), INTENT(OUT) :: qsol      ! Soil Water Content
        REAL, DIMENSION(klon), INTENT(INOUT) :: z0m    ! Momentum Roughn Lgt
        REAL, DIMENSION(klon), INTENT(INOUT) :: z0h    ! Momentum Roughn Lgt

        ! Output Variables for LMDZ
        REAL, DIMENSION(klon), INTENT(OUT) :: alb1      ! Albedo SW
        REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW
        REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo
        REAL, DIMENSION(klon), INTENT(OUT) :: emis_new  ! Surface Emissivity
        REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff
        REAL, DIMENSION(klon), INTENT(OUT) :: ffonte    ! enthalpy flux due to surface melting
        REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte   ! water flux due to surface melting
        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s   ! d/dT sens. ht flux
        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l   ! d/dT latent ht flux
        REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens  ! Sensible ht flux
        REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat   ! Latent heat flux
        REAL, DIMENSION(klon), INTENT(OUT) :: evap      ! Evaporation
        REAL, DIMENSION(klon), INTENT(OUT) :: erod      ! Erosion of surface snow (flux)
        REAL, DIMENSION(klon), INTENT(OUT) :: agesno    ! Snow age (top layer)
        REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature
        REAL, DIMENSION(klon), INTENT(OUT) :: qsurf     ! Surface Humidity

        ! Specific INLANDIS outputs
        REAL, DIMENSION(klon), INTENT(OUT) :: qsnow     ! Total H2O snow[kg/m2]
        REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt   ! Snow height (m)
        REAL, DIMENSION(klon), INTENT(OUT) :: to_ice    ! Snow passed to ice
        REAL, DIMENSION(klon), INTENT(OUT) :: sissnow   ! Snow in model (kg/m2)

        ! +--Internal  Variables
        ! +  ===================

        CHARACTER(len = 20) :: fn_outfor ! Name for output file
        CHARACTER (len = 80)              :: abort_message
        CHARACTER (len = 20)              :: modname = 'surf_inlandsis_mod'

        INTEGER :: i, ig, ikl, isl, isn, nt
        INTEGER :: gp_outfor, un_outfor
        REAL, PARAMETER :: f1 = 0.5
        REAL, PARAMETER :: sn_upp = 10000., sn_low = 500.
        REAL, PARAMETER :: sn_add = 400., sn_div = 2.
        ! snow mass upper,lower limit,
        ! added mass/division lowest layer
        REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6
        REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3
        ! Parameters for drainage
        ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
        ! +...        Run Off Parameters
        ! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)
        ! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)

        REAL, DIMENSION(klon) :: eps0SL          ! surface Emissivity
        REAL :: zsigma, Ua_min, Us_min, lati
        REAL, PARAMETER :: cdmax=0.05
        REAL :: lambda          ! Par. soil discret.
        REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2         ! Soil layer thicknesses
        !$OMP THREADPRIVATE(dz1,dz2)
        LOGICAL, SAVE :: firstcall
        !$OMP THREADPRIVATE(firstcall)

        INTEGER :: iso
        LOGICAL :: file_exists
        CHARACTER(len = 20) :: fichnom
        LOGICAL :: is_init_domec
        ! CA initialization
        ! dz_profil_15 : 1 m in 15 layers [m]
        real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, &
                                                0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/)
        ! mean_temp : mean annual surface temperature [K]
        real, dimension(klon) :: mean_temp
        ! mean_dens : mean surface density [kg/m3]
        real, dimension(klon) :: mean_dens
        ! lat_scale : temperature lapse rate against latitude [K degree-1]
        real :: lat_scale
        ! sh_scale : temperature lapse rate against altitude [K km-1]
        real :: sh_scale
        ! variables for density profile
        ! E0, E1 : exponent
        real :: E0, E1
        ! depth at which 550 kg m-3 is reached [m]
        real :: z550
        ! depths of snow layers
        real :: depth, snow_depth, distup
        ! number of initial snow layers
        integer :: nb_snow_layer
        ! For density calc.
        real :: alpha0, alpha1, ln_smb
        ! theoritical densities [kg m-3]
        real :: rho0, rho1, rho1_550
        ! constants for density profile
        ! C0, C1 : constant, 0.07 for z <= 550 kg m-3
        real, parameter :: C0 = 0.07
        real, parameter :: C1 = 0.03
        ! rho_i : ice density [kg m-3]
        real, parameter :: rho_ice = 917.
        ! E_c : activation energy [J mol-1]
        real, parameter :: E_c = 60000.
        ! E_g : activation energy [J mol-1]
        real, parameter :: E_g = 42400.
        ! R : gas constant [J mol-1 K-1]
        real, parameter :: R = 8.3144621

      
     


        ! + PROGRAM START
        ! + -----------------------------------------

        zsigma = 1000.
        dt__SV = dtime

        IF (debut) THEN
            firstcall = .TRUE.
            INI_SV = .FALSE.
        ELSE
            firstcall = .FALSE.
            INI_SV = .TRUE.
        END IF

        IF (ok_outfor) THEN
            un_outfor = 51    ! unit number for point output file
            gp_outfor = 1    ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC
            fn_outfor = 'outfor_SV.dat'
        END IF

        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! + INITIALISATION: BEGIN +++
        ! + -----------------------------------------
        IF (firstcall) THEN

            ! +--Array size
            ! +  -----------------------

            klonv = klon
            knonv = knon
                write(*, *) 'ikl, lon and lat in INLANDSIS'

            DO ikl = 1, knon
                i=ikl2i(ikl)
                write(*, *) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i)
            END DO

            ! +--Variables initizialisation
            ! +  ---------------------------

            CALL INIT_VARtSV
            CALL INIT_VARxSV
            CALL INIT_VARySV



            ! +--Surface Fall Line Slope
            ! +  -----------------------
            IF (SnoMod)  THEN
                DO ikl = 1, knon
                    slopSV(ikl) = slope(ikl)
                    SWf_SV(ikl) = &   ! Normalized Decay of the
                            exp(-dt__SV             &   ! Surficial Water Content
                                    / (c1_zuo                &   !(Zuo and Oerlemans 1996,
                                            + c2_zuo * exp(-c3_zuo * abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
                END DO
            END IF



            ! +--Soil layer thickness . Compute soil discretization (as for LMDZ)
            ! +  ----------------------------------------------------------------
            !        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
            CALL get_soil_levels(dz1, dz2, lambda)

            lambSV = lambda
            dz1_SV(1:knon, 1:) = 0.
            dz2_SV(1:knon, 1:) = 0.

            DO isl = -nsol, 0
                dz_dSV(isl) = 0.5e-3 * dz2(1 - isl)           ! Soil layer thickness
                DO ikl = 1, knon
                    dz1_SV(ikl, isl) = dz1(1 - isl)    !1.e-3*
                    dz2_SV(ikl, isl) = dz2(1 - isl)    !1.e-3*
                END DO
            END DO


            ! Set variables
            ! =============
            DO ikl = 1, knon
                ! LSmask : Land/Sea Mask
                LSmask(ikl) = 1
                ! isotSV : Soil Type -> 12 = ice
                isotSV(ikl) = 12
                ! iWaFSV : Soil Drainage (1,0)=(y,n)
                iWaFSV(ikl) = 1
                ! eps0SL : Surface Emissivity
                eps0SL(ikl) = 1.
                ! alb0SV : Soil Albedo
                alb0SV(ikl) = alb_soil(ikl)
                ! Tsf_SV : Surface Temperature, must be bellow freezing
                Tsf_SV(ikl) = min(temp_air(ikl), TfSnow)
            END DO

            ! +--Initialization of soil and snow variables in case startsis is not read
            ! +  ----------------------------------------------------------------------

            is_init_domec=.FALSE.


            IF (is_init_domec) THEN
            ! Coarse initilization inspired from vertcical profiles at Dome C,
            ! Antarctic Plateaui (10m of snow, 19 levels)

                 DO ikl = 1,knon
! + Soil
                 DO isl =   -nsol,0    
                   TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2)   !temp_air(ikl)
                   !tsoil(ikl,1-isl)   Soil Temperature
                   !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
                   eta_SV(ikl,isl) = epsi           !etasoil(ikl,1-isl) Soil Water[m3/m3]
                   ro__SV(ikl,isl) = rhoIce         !rosoil(ikl,1-isl) volumic mass
                 END DO   


           ! Snow
                 isnoSV(ikl) = 19
                 istoSV(ikl, 1:isnoSV(ikl)) = 100
                 ro__SV(ikl, 1:isnoSV(ikl)) = 350.
                 eta_SV(ikl, 1:isnoSV(ikl)) = epsi
                 TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2)
                 G1snSV(ikl, 1:isnoSV(ikl)) = 99.
                 G2snSV(ikl, 1:isnoSV(ikl)) = 2.
                 agsnSV(ikl, 1:isnoSV(ikl)) = 50.
                 dzsnSV(ikl, 19) = 0.015
                 dzsnSV(ikl, 18) = 0.015
                 dzsnSV(ikl, 17) = 0.020
                 dzsnSV(ikl, 16) = 0.030
                 dzsnSV(ikl, 15) = 0.040
                 dzsnSV(ikl, 14) = 0.060
                 dzsnSV(ikl, 13) = 0.080
                 dzsnSV(ikl, 12) = 0.110
                 dzsnSV(ikl, 11) = 0.150
                 dzsnSV(ikl, 10) = 0.200
                 dzsnSV(ikl, 9) = 0.300
                 dzsnSV(ikl, 8) = 0.420
                 dzsnSV(ikl, 7) = 0.780
                 dzsnSV(ikl, 6) = 1.020
                 dzsnSV(ikl, 5) = 0.980
                 dzsnSV(ikl, 4) = 1.020
                 dzsnSV(ikl, 3) = 3.970
                 dzsnSV(ikl, 2) = 1.020
                 dzsnSV(ikl, 1) = 1.020

                 END DO
            ELSE

            ! Initilialisation with climatological temperature and density
            ! profiles as in MAR. Methodology developed by Cecile Agosta
  
            ! initialize with 0., for unused snow layers
            dzsnSV = 0.
            G1snSV = 0.
            G2snSV = 0.
            istoSV = 0
            TsisSV = 0.


            ! initialize mean variables (unrealistic)
            mean_temp = TfSnow
            mean_dens = 300.
            ! loop on grid cells
            DO ikl = 1, knon
            	lati=rlat(ikl2i(ikl))
                ! approximations for mean_temp and mean_dens
                ! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1)
                ! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1
                ! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica,
                ! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed.
                if (lati > 60.) then
                    ! CA todo : add longitude bounds
                    ! Greenland mean temperature : function of altitude and latitude
                    ! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1
                    lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9
                    lat_scale = max(min(lat_scale, 0.9), 0.75)
                    ! sh_scale equals the environmental lapse rate : 6.5 °C km-1
                    sh_scale = 6.5
                    mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.)
                    ! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051
                    mean_dens(ikl) = 315.
                else if (lati < -60.) then
                    ! Antarctica mean temperature : function of altitude and latitude
                    ! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1
                    lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3
                    lat_scale = max(min(lat_scale, 1.3), 0.6)
                    ! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1
                    sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5
                    sh_scale = max(min(sh_scale, 9.8), 6.5)
                    mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.)
                    ! Antarctica surface density : function of mean annual temperature
                    ! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013)
                    ! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m.
                    ! Weinhart et al 2020  https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201
                    ! 320 kg m-3 is reached at Dome A, 4100 m a.s.l.
                    ! Dome C : st_ant_param(3233, -75.1) = -47.7
                    ! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7
                    mean_dens(ikl) =  (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450.
                    mean_dens(ikl) = min(450., max(320., mean_dens(ikl)))
                else

                !    write(*, *) 'Attention: temperature initialization is only defined for Greenland and Antarctica'

                     mean_dens(ikl) =350.
                     mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2)

                !abort_message='temperature initialization is only defined for Greenland and Antarctica'
                !CALL abort_physic(modname,abort_message,1)

                end if
 
                ! mean_temp is defined for ice ground only
                mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2)

                ! Soil layers
                ! ===========
                DO isl = -nsol, 0
                    ! TsisSV : Temperature [K]
                    TsisSV(ikl, isl) = mean_temp(ikl)
                    ! eta_SV : Soil Water [m3/m3]
                    eta_SV(ikl, isl) = epsi
                    ! ro__SV : Volumic Mass [kg/m3]
                    ro__SV(ikl, isl) = rhoIce
                END DO

                ! Snow layers
                ! ===========
                ! snow_depth : initial snow depth
                snow_depth = 20.
                ! nb_snow_layer : initial nb of snow layers
                nb_snow_layer = 15
                ! isnoSV : total nb of snow layers
                isnoSV(ikl) = nb_snow_layer
                ! depth : depth of each layer
                depth = snow_depth
                do isl = 1, nb_snow_layer
                    ! dzsnSV : snow layer thickness
                    dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1))
                    ! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical
                    G1snSV(ikl, isl) = 99.
                    ! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size
                    G2snSV(ikl, isl) = 3.
                    agsnSV(ikl, isl) = 0.
                    istoSV(ikl, isl) = 0
                    ! eta_SV : Liquid Water Content [m3/m3]
                    eta_SV(ikl, isl) = 0.
                    ! distance to surface
                    depth = depth - dzsnSV(ikl,isl) / 2.
                    distup = min(1., max(0., depth / snow_depth))
                    ! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom)
                    TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup**2) + mean_temp(ikl) * distup**2
                    ! firn density : densification formulas from :
                    ! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/)
                    ! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306)
                    ! Integration of the steady state equation
                    ! ln_smb approximated as a function of temperature
                    ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.)
                    ! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1
                    alpha0 = max(1.435 - 0.151 * ln_smb, 0.25)
                    alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25)
                    E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0
                    E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1
                    z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0
                    rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
                    rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice
                    if (depth <= z550) then
                        ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
                    else
                        ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice
                    end if
                    depth = depth - dzsnSV(ikl,isl) / 2.
                    
                END DO

            END DO

            END IF


            ! + Numerics paramaters, SISVAT_ini 
            ! +  ----------------------
            CALL SISVAT_ini(knon)


            ! +--Read restart file
            ! +  =================================================

            INQUIRE(FILE = "startsis.nc", EXIST = file_exists)
            IF (file_exists) THEN
                CALL sisvatetat0("startsis.nc", ikl2i)
            END IF



            ! +--Output ascii file
            ! +  =================================================

            ! open output file
            IF (ok_outfor) THEN
                open(unit = un_outfor, status = 'replace', file = fn_outfor)
                ikl = gp_outfor     ! index sur la grille land ice
                write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl))
                write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] &
   G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]'
            END IF

        END IF  ! firstcall
        ! +
        ! +  +++  INITIALISATION:  END  +++
        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! + READ FORCINGS
        ! + ------------------------

        ! + Update Forcings for SISVAT given by the LMDZ model.
        ! +
        DO ikl = 1, knon

            ! +--Atmospheric Forcing                                    (INPUT)
            ! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^
            za__SV(ikl) = zsl_height(ikl)               ! surface layer height (fisr model level) [m]
            Ua_min = 0.2 * sqrt(za__SV(ikl))            !
            VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
            TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
            ExnrSV(ikl) = pexner(ikl)                   ! Exner potential
            rhT_SV(ikl) = dens_air(ikl)                 ! Air density
            QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity
            ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
            p1l_SV(ikl) = p1lay(ikl)                    ! lowest atm. layer press[Pa]

            ! +--Surface properties
            ! +  ^^^^^^^^^^^^^^^^^^

            Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
            Z0h_SV(ikl) = z0h(ikl)                      ! Moment.Roughn.L.

            ! +--Energy Fluxes                                          (INPUT)
            ! +  ^^^^^^^^^^^^^                                           ^^^^^
            coszSV(ikl) = max(czemin, rmu0(ikl))         ! cos(zenith.Dist.)
            sol_SV(ikl) = swdown(ikl)                   ! downward Solar
            IRd_SV(ikl) = lwdown(ikl)                   ! downward IR
            rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.

            ! +--Water  Fluxes                                          (INPUT)
            ! +  ^^^^^^^^^^^^^                                           ^^^^^
            drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
            dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]

            ! #BS    dbs_SV(ikl) = blowSN(i,j,n)
            ! dbs_SV = Maximum potential erosion amount [kg/m2]
            ! => Upper bound for eroded snow mass
            !        uss_SV(ikl) = SLussl(i,j,n) ! u*qs* (only for Tv in sisvatesbl.f)
            ! #BS  if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then
            ! #BS    dsnbSV(ikl) =1.0-min(qsHY(i,j,kB)     !BS neglib. at kb ~100 magl)
            ! #BS.                        /max(qshy(i,j,mz),eps9),unun)
            ! #BS    dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl))
            ! #BS    dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl)))
            ! #BS  else
            ! #BS    dsnbSV(ikl) = 0.
            ! #BS  endif
            !      dsnbSV is the drift fraction of deposited snow updated in sisvat.f
            !      will be used for characterizing the Buffer Layer
            !      (see update of  Bros_N, G1same, G2same, zroOLD, zroNEW)
            ! #BS  if(n==1) qbs_HY(i,j) = dsnbSV(ikl)
            qsnoSV(ikl) = snow_cont_air(ikl)



            ! +--Soil/BL                                      (INPUT)
            ! +  ^^^^^^^                                       ^^^^^
            alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo
            AcoHSV(ikl) = AcoefH(ikl)
            BcoHSV(ikl) = BcoefH(ikl)
            AcoQSV(ikl) = AcoefQ(ikl)
            BcoQSV(ikl) = BcoefQ(ikl)
            cdH_SV(ikl) = min(cdragh(ikl),cdmax)
            cdM_SV(ikl) = min(cdragm(ikl),cdmax)
            rcdmSV(ikl) = sqrt(cdM_SV(ikl))
            Us_min = 0.01
            us__SV(ikl) = max(Us_min, ustar(ikl))
            ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6))
            rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6))

            ! +--Energy Fluxes                                          (INPUT/OUTPUT)
            ! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^
            !IF (.not.firstcall) THEN
            Tsrfsv(ikl)  = tsurf(ikl)                     !hj 12 03 2010
            cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness
            !END IF

         END DO

        ! +  +++  READ FORCINGS:  END  +++
        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 

        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! +--SISVAT EXECUTION
        ! +  ----------------

        call  INLANDSIS(SnoMod, BloMod, 1)


       
        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        ! + RETURN RESULTS
        ! + --------------
        ! + Return (compressed) SISVAT variables to LMDZ
        ! +
        DO  ikl = 1, knon                  ! use only 1:knon (actual ice sheet..)
            dflux_s(ikl) = dSdTSV(ikl)         ! Sens.H.Flux T-Der.
            dflux_l(ikl) = dLdTSV(ikl)         ! Latn.H.Flux T-Der.
            fluxsens(ikl) = HSs_sv(ikl)         ! HS
            fluxlat(ikl) = HLs_sv(ikl)         ! HL
            evap(ikl) = -1*HLs_sv(ikl) / LHvH2O  ! Evaporation
            erod(ikl) = 0.

            IF (BloMod) THEN
                ! + Blowing snow

                !       SLussl(i,j,n)= 0.
                ! #BS   SLussl(i,j,n)=                     !Effective erosion
                ! #BS. (- dbs_ER(ikl))/(dt*rhT_SV(ikl))    !~u*qs* from previous time step
                ! #BS   blowSN(i,j,n)=  dt*uss_SV(ikl)     !New max. pot. Erosion [kg/m2]
                ! #BS.                    *rhT_SV(ikl)     !(further bounded in sisvat_bsn.f)
                ! #BS  erprev(i,j,n) =     dbs_Er(ikl)/dt__SV
                erod(ikl) = dbs_Er(ikl) / dt__SV
            ENDIF

            ! +   Check snow thickness,  substract if too thick, add if too thin

            sissnow(ikl) = 0.  !()
            DO  isn = 1, isnoSV(ikl)
                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
            END DO

            IF (sissnow(ikl) <= sn_low) THEN  !add snow
                IF (isnoSV(ikl)>=1) THEN
                    dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi)
                    toicSV(ikl) = toicSV(ikl) - sn_add
                ELSE
                    write(*, *) 'Attention, bare ice... point ', ikl
                    isnoSV(ikl) = 1
                    istoSV(ikl, 1) = 0
                    ro__SV(ikl, 1) = 350.
                    dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi)  ! 1.
                    eta_SV(ikl, 1) = epsi
                    TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2)
                    G1snSV(ikl, 1) = 0.
                    G2snSV(ikl, 1) = 0.3
                    agsnSV(ikl, 1) = 10.
                    toicSV(ikl) = toicSV(ikl) - sn_add
                END IF
            END IF

            IF (sissnow(ikl) >= sn_upp) THEN  !thinnen snow layer below
                dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div
                toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div
            END IF

            sissnow(ikl) = 0.
            qsnow(ikl) = 0.
            snow(ikl) = 0.
            snowhgt(ikl) = 0.

            DO  isn = 1, isnoSV(ikl)
                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
                snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn)
                ! Etienne: check calc qsnow
                qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn)
            END DO

            zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0)
            ! Etienne: comment following line
            ! snow(ikl)    = sissnow(ikl)+toicSV(ikl)
            snow(ikl) = sissnow(ikl)

            to_ice(ikl) = toicSV(ikl)
            runoff_lic(ikl) = RnofSV(ikl)    ! RunOFF: intensity (flux due to melting + liquid precip)
            fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing
            ffonte(ikl)=fqfonte(ikl)*Lf_H2O

            qsol(ikl) = 0.
            DO  isl = -nsol, 0
                tsoil(ikl, 1 - isl) = TsisSV(ikl, isl)       ! Soil Temperature
                ! Etienne: check calc qsol
                qsol(ikl) = qsol(ikl)                      &
                        + eta_SV(ikl, isl) * dz_dSV(isl)
            END DO
            agesno(ikl) = agsnSV(ikl, isnoSV(ikl))        !          [day]

            alb1(ikl) = alb1sv(ikl)             ! Albedo VIS
!            alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl)                   &
!                    & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1
            alb2(ikl)=alb2sv(ikl)
            ! Albedo NIR
            alb3(ikl) = alb3sv(ikl)             ! Albedo FIR
            ! 6 band Albedo
            alb6(ikl,:)=alb6sv(ikl,:)

            tsurf_new(ikl) = Tsrfsv(ikl)

            qsurf(ikl) = QsT_SV(ikl)
            emis_new(ikl) = eps0SL(ikl)
            z0m(ikl) = Z0m_SV(ikl)
            z0h(ikl) = Z0h_SV(ikl)


        END DO

            IF (ok_outfor) THEN
             ikl= gp_outfor
            write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++'
            write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl)
            write(un_outfor, *) dzsnSV(ikl, :)
            write(un_outfor, *) TsisSV(ikl, :)
            write(un_outfor, *) ro__SV(ikl, :)
            write(un_outfor, *) eta_SV(ikl, :)
            write(un_outfor, *) G1snSV(ikl, :)
            write(un_outfor, *) G2snSV(ikl, :)
            write(un_outfor, *) agsnSV(ikl, :)
            write(un_outfor, *) istoSV(ikl, :)
            write(un_outfor, *) DOPsnSV(ikl, :)
        ENDIF



        ! +  -----------------------------
        ! +  END --- RETURN RESULTS
        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        IF (lafin) THEN
            fichnom = "restartsis.nc"
            CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat)

            IF (ok_outfor) THEN
                close(unit = un_outfor)
            END IF
        END IF

    END SUBROUTINE surf_inlandsis


    !=======================================================================

    SUBROUTINE get_soil_levels(dz1, dz2, lambda)
        ! ======================================================================
        ! Routine to compute the vertical discretization of the soil in analogy
        ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
        ! of SISVAT, therefore it's needed here.

        USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
        USE mod_phys_lmdz_para
        USE VAR_SV


        !    INCLUDE "dimsoil.h"

        REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
        REAL, INTENT(OUT) :: lambda


        !-----------------------------------------------------------------------
        !   Depthts:
        !   --------
        REAL fz, rk, fz1, rk1, rk2
        REAL min_period, dalph_soil
        INTEGER ierr, jk

        fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)

        !    write(*,*)'Start soil level computation'
        !-----------------------------------------------------------------------
        ! Calculation of some constants
        ! NB! These constants do not depend on the sub-surfaces
        !-----------------------------------------------------------------------
        !-----------------------------------------------------------------------
        !   ground levels
        !   grnd=z/l where l is the skin depth of the diurnal cycle:
        !-----------------------------------------------------------------------

        min_period = 1800. ! en secondes
        dalph_soil = 2.    ! rapport entre les epaisseurs de 2 couches succ.
        ! !$OMP MASTER
        !     IF (is_mpi_root) THEN
        !        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
        !        IF (ierr == 0) THEN ! Read file only if it exists
        !           READ(99,*) min_period
        !           READ(99,*) dalph_soil
        !           PRINT*,'Discretization for the soil model'
        !           PRINT*,'First level e-folding depth',min_period, &
        !                '   dalph',dalph_soil
        !           CLOSE(99)
        !        END IF
        !     ENDIF
        ! !$OMP END MASTER
        !     CALL bcast(min_period)
        !     CALL bcast(dalph_soil)

        !   la premiere couche represente un dixieme de cycle diurne
        fz1 = SQRT(min_period / 3.14)

        DO jk = 1, nsoilmx
            rk1 = jk
            rk2 = jk - 1
            dz2(jk) = fz(rk1) - fz(rk2)
        ENDDO
        DO jk = 1, nsoilmx - 1
            rk1 = jk + .5
            rk2 = jk - .5
            dz1(jk) = 1. / (fz(rk1) - fz(rk2))
        ENDDO
        lambda = fz(.5) * dz1(1)
        DO jk = 1, nsoilmx
            rk = jk
            rk1 = jk + .5
            rk2 = jk - .5
        ENDDO

    END SUBROUTINE get_soil_levels


    !===========================================================================

    SUBROUTINE SISVAT_ini(knon)

        !C +------------------------------------------------------------------------+
        !C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR |
        !C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters |
        !C +------------------------------------------------------------------------+
        !C |   PARAMETERS:  klonv: Total Number of columns =                        |
        !C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
        !C |                     X       Number of Mosaic Cell per grid box         |
        !C |                                                                        |
        !C |   INPUT:   dt__SV   : Time  Step                                   [s] |
        !C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] |
        !C |                                                                        |
        !C |   OUTPUT:             [-] |
        !C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] |
        !C |            etamSV   : Soil Minimum Humidity                    [m3/m3] |
        !C |                      (based on a prescribed Soil Relative Humidity)    |
        !C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       |
        !C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        |
        !C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   |
        !C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] |
        !C |            dzsnSV(0): Soil first Layer Thickness                   [m] |
        !C |            dzmiSV   : Distance between two contiguous levels       [m] |
        !C |            dz78SV   : 7/8 (Layer Thickness)                        [m] |
        !C |            dz34SV   : 3/4 (Layer Thickness)                        [m] |
        !C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] |
        !C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] |
        !C |            dtz_SV   : dt/dz                                      [s/m] |
        !C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
        !C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      |
        !C |            Explic   : Explicit Parameter = 1.0 - Implic                |
        !C |                                                                        |
        !C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
        !C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
        !C | #          #SH: Hapex-Sahel Values                                     !
        !C |                                                                        |
        !C +------------------------------------------------------------------------+


        !C +--Global Variables
        !C +  ================

        USE dimphy
        USE VARphy
        USE VAR_SV
        USE VARdSV
        USE VAR0SV
        USE VARxSV
        USE VARtSV
        USE VARxSV
        USE VARySV
        IMPLICIT NONE



        !C +--Arguments
        !C +  ==================
        INTEGER, INTENT(IN) :: knon

        !C +--Internal Variables
        !C +  ==================

        INTEGER :: ivt, ist, ikl, isl, isn, ikh
        INTEGER :: misl_2, nisl_2
        REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2
        REAL, PARAMETER :: RHsMin = 0.001        ! Min.Soil Relative Humidity
        REAL :: PsiMax                        ! Max.Soil Water    Potential
        REAL :: a_Khyd, b_Khyd                 ! Water conductivity


        !c #WR REAL    ::  Khyd_x,Khyd_y



        !C +--Non Time Dependant SISVAT parameters
        !C +  ====================================

        !C +--Soil Discretization
        !C +  -------------------

        !C +--Numerical Scheme Parameters
        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
        Implic = 0.75                           ! 0.5  <==> Crank-Nicholson  
        Explic = 1.00 - Implic                  !                            

        !C +--Soil/Snow Layers Indices
        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^
        DO  isl = -nsol, 0
            islpSV(isl) = isl + 1
            islpSV(isl) = min(islpSV(isl), 0)
            islmSV(isl) = isl - 1
            islmSV(isl) = max(-nsol, islmSV(isl))
        END DO

        DO  isn = 1, nsno
            isnpSV(isn) = isn + 1
            isnpSV(isn) = min(isnpSV(isn), nsno)
        END DO

        !C +--Soil      Layers Thicknesses
        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
        !c #kd IF (nsol.gt.4)                                              THEN
        !c #kd   DO isl=-5,-nsol,-1
        !c #kd     dz_dSV(isl)=   1.
        !c #kd   END DO
        !c #kd END IF

        !      IF (nsol.ne.4)                                              THEN
        !        DO isl= 0,-nsol,-1
        !          misl_2 =     -mod(isl,2)
        !          nisl_2 =         -isl/2
        !          dz_dSV(isl)=(((1-misl_2) * 0.001
        !     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.
        !C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm

        !c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001
        !c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.

        !c #05     dz_dSV(isl)=(((1-misl_2) * 0.001
        !c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5
        !        END DO
        !          dz_dSV(0)  =               0.001
        !          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004
        !      END IF

        zz_dSV = 0.
        DO  isl = -nsol, 0
            dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl)))
            dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl)
            dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl)
            dtz_SV(isl) = dt__SV / dz_dSV(isl)
            dtz_SV2(isl) = 1. / dz_dSV(isl)
            dz78SV(isl) = 0.875 * dz_dSV(isl)
            dz34SV(isl) = 0.750 * dz_dSV(isl)
            dz_8SV(isl) = 0.125 * dz_dSV(isl)
            dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl))                        &
   + 0.750 * dz_dSV(isl)                                &
   + 0.125 * dz_dSV(islpSV(isl))
            zz_dSV = zz_dSV + dz_dSV(isl)
        END DO
        DO ikl = 1, knon !v
            dzsnSV(ikl, 0) = dz_dSV(0)
        END DO

        !C +--Conversion to a 50 m Swab Ocean Discretization
        !C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        OcndSV = 0.
        DO isl = -nsol, 0
            OcndSV = OcndSV + dz_dSV(isl)
        END DO
        OcndSV = 50. / OcndSV


        !C +--Secondary Soil       Parameters
        !C +  -------------------------------

        DO  ist = 0, nsot
            rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6   ! Soil Contrib. to (ro c)_s
            s1__SV(ist) = bCHdSV(ist)          & ! Factor of (eta)**(b+2)
   * psidSV(ist) * Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)
   / (etadSV(ist)**(bCHdSV(ist) + 3.))     !
            s2__SV(ist) = Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)
   / (etadSV(ist)**(2. * bCHdSV(ist) + 3.))     !    in DR97, Eqn.(3.35)

            !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
            !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
            Psimax = -(log(RHsMin)) / 7.2E-5        ! DR97, Eqn 3.15 Inversion
            etamSV(ist) = etadSV(ist)                                      &
   * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist)))
        END DO
        etamSV(12) = 0.

        !C +--Piecewise Hydraulic Conductivity Profiles
        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        DO   ist = 0, nsot

            d__eta = etadSV(ist) / nkhy
            eta__1 = 0.
            eta__2 = d__eta
            DO ikh = 0, nkhy
                Khyd_1 = s2__SV(ist)             & ! DR97, Eqn.(3.35)
   * (eta__1      **(2. * bCHdSV(ist) + 3.))        !
                Khyd_2 = s2__SV(ist)             &!
   * (eta__2      **(2. * bCHdSV(ist) + 3.))        !

                a_Khyd = (Khyd_2 - Khyd_1) / d__eta   !
                b_Khyd = Khyd_1 - a_Khyd * eta__1   !
                !c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !
                !c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !
                aKdtSV(ist, ikh) = a_Khyd * dt__SV   !
                bKdtSV(ist, ikh) = b_Khyd * dt__SV   !

                eta__1 = eta__1 + d__eta
                eta__2 = eta__2 + d__eta
            END DO
        END DO

        return

    END SUBROUTINE SISVAT_ini


    !***************************************************************************

    SUBROUTINE sisvatetat0 (fichnom, ikl2i)

        USE dimphy
        USE mod_grid_phy_lmdz
        USE mod_phys_lmdz_para

        USE iostart
        USE VAR_SV
        USE VARdSV
        USE VARxSV
        USE VARtSV
        USE indice_sol_mod

        IMPLICIT none
        !======================================================================
        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
        ! Objet: Lecture du fichier de conditions initiales pour SISVAT
        !======================================================================
        include "clesphys.h"
        include "compbl.h"

        !======================================================================
        CHARACTER(LEN = *) :: fichnom

        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
        REAL, DIMENSION(klon) :: rlon
        REAL, DIMENSION(klon) :: rlat

        ! les variables globales ecrites dans le fichier restart
        REAL, DIMENSION(klon) :: isno
        REAL, DIMENSION(klon) :: ispi
        REAL, DIMENSION(klon) :: iice
        REAL, DIMENSION(klon) :: rusn
        REAL, DIMENSION(klon, nsno) :: isto

        REAL, DIMENSION(klon, nsismx) :: Tsis
        REAL, DIMENSION(klon, nsismx) :: eta
        REAL, DIMENSION(klon, nsismx) :: ro

        REAL, DIMENSION(klon, nsno) :: dzsn
        REAL, DIMENSION(klon, nsno) :: G1sn
        REAL, DIMENSION(klon, nsno) :: G2sn
        REAL, DIMENSION(klon, nsno) :: agsn

        REAL, DIMENSION(klon) :: toic

        INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts
        CHARACTER (len = 2) :: str2
        LOGICAL :: found

        errT = 0
        errro = 0
        erreta = 0
        errdz = 0
        snopts = 0
        ! Ouvrir le fichier contenant l'etat initial:

        CALL open_startphy(fichnom)

        ! Lecture des latitudes, longitudes (coordonnees):

        CALL get_field("latitude", rlat, found)
        CALL get_field("longitude", rlon, found)

        CALL get_field("n_snows", isno, found)
        IF (.NOT. found) THEN
            PRINT*, 'phyetat0: Le champ <n_snows> est absent'
            PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
        ENDIF

        CALL get_field("n_ice_top", ispi, found)
        CALL get_field("n_ice", iice, found)
        CALL get_field("surf_water", rusn, found)


        CALL get_field("to_ice", toic, found)
        IF (.NOT. found) THEN
            PRINT*, 'phyetat0: Le champ <to_ice> est absent'
            toic(:) = 0.
        ENDIF

        DO isn = 1, nsno
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("AGESNOW" // str2, &
                        agsn(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        DO isn = 1, nsno
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("DZSNOW" // str2, &
                        dzsn(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        DO isn = 1, nsno
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("G2SNOW" // str2, &
                        G2sn(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        DO isn = 1, nsno
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("G1SNOW" // str2, &
                        G1sn(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        DO isn = 1, nsismx
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("ETA" // str2, &
                        eta(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        DO isn = 1, nsismx
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("RO" // str2, &
                        ro(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        DO isn = 1, nsismx
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("TSS" // str2, &
                        Tsis(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        DO isn = 1, nsno
            IF (isn<=99) THEN
                WRITE(str2, '(i2.2)') isn
                CALL get_field("HISTORY" // str2, &
                        isto(:, isn), found)
            ELSE
                PRINT*, "Trop de couches"
                CALL abort
            ENDIF
        ENDDO
        write(*, *)'Read ', fichnom, ' finished!!'

        !*********************************************************************************
        ! Compress restart file variables for SISVAT

        DO  ikl = 1, klon
            i = ikl2i(ikl)
            IF (i > 0) THEN
                isnoSV(ikl) = INT(isno(i))          ! Nb Snow/Ice Lay.
                ispiSV(ikl) = INT(ispi(i))          ! Nb Supr.Ice Lay.
                iiceSV(ikl) = INT(iice(i))          ! Nb      Ice Lay.

                DO isl = -nsol, 0
                    ro__SV(ikl, isl) = ro(i, nsno + 1 - isl)       !
                    eta_SV(ikl, isl) = eta(i, nsno + 1 - isl)         ! Soil Humidity
                    !hjp 15/10/2010
                    IF (eta_SV(ikl, isl) <= 1.e-6) THEN          !hj check
                        eta_SV(ikl, isl) = 1.e-6
                    ENDIF
                    TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl)        ! Soil Temperature
                    IF (TsisSV(ikl, isl) <= 1.) THEN             !hj check
                        !                errT=errT+1
                        TsisSV(ikl, isl) = 273.15 - 0.2              ! Etienne: negative temperature since soil is ice
                    ENDIF

                END DO
                write(*, *)'Copy histo', ikl

                DO  isn = 1, isnoSV(ikl) !nsno
                    snopts = snopts + 1
                    IF (isto(i, isn) > 10.) THEN          !hj check
                        write(*, *)'Irregular isto', ikl, i, isn, isto(i, isn)
                        isto(i, isn) = 1.
                    ENDIF

                    istoSV(ikl, isn) = INT(isto(i, isn))     ! Snow     History
                    ro__SV(ikl, isn) = ro(i, isn)            !        [kg/m3]
                    eta_SV(ikl, isn) = eta(i, isn)           !        [m3/m3]
                    TsisSV(ikl, isn) = Tsis(i, isn)          !            [K]

                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
                        errT = errT + 1
                        TsisSV(ikl, isn) = TsisSV(ikl, 0)
                    ENDIF
                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
                        TsisSV(ikl, isn) = 263.15
                    ENDIF
                    IF (eta_SV(ikl, isn) < 1.e-9) THEN          !hj check
                        eta_SV(ikl, isn) = 1.e-6
                        erreta = erreta + 1
                    ENDIF
                    IF (ro__SV(ikl, isn) <= 10.) THEN          !hj check
                        ro__SV(ikl, isn) = 11.
                        errro = errro + 1
                    ENDIF
                    write(*, *)ikl, i, isn, Tsis(i, isn), G1sn(i, isn)
                    G1snSV(ikl, isn) = G1sn(i, isn)          ! [-]        [-]
                    G2snSV(ikl, isn) = G2sn(i, isn)          ! [-] [0.0001 m]
                    dzsnSV(ikl, isn) = dzsn(i, isn)          !            [m]
                    agsnSV(ikl, isn) = agsn(i, isn)          !          [day]
                END DO
                rusnSV(ikl) = rusn(i)              ! Surficial Water
                toicSV(ikl) = toic(i)              ! bilan snow to ice
            END IF
        END DO

    END SUBROUTINE sisvatetat0


    !======================================================================
    SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat)



        !======================================================================
        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
        ! Objet: Ecriture de l'etat de redemarrage pour SISVAT
        !======================================================================
        USE mod_grid_phy_lmdz
        USE mod_phys_lmdz_para
        USE iostart
        USE VAR_SV
        USE VARxSV
        USE VARySV !hj tmp 12 03 2010
        USE VARtSV
        USE indice_sol_mod
        USE dimphy

        IMPLICIT none

        include "clesphys.h"
        include "compbl.h"

        !======================================================================

        CHARACTER(LEN = *) :: fichnom
        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
        REAL, DIMENSION(klon), INTENT(IN) :: rlon
        REAL, DIMENSION(klon), INTENT(IN) :: rlat

        ! les variables globales ecrites dans le fichier restart
        REAL, DIMENSION(klon) :: isno
        REAL, DIMENSION(klon) :: ispi
        REAL, DIMENSION(klon) :: iice
        REAL, DIMENSION(klon, nsnowmx) :: isto

        REAL, DIMENSION(klon, nsismx) :: Tsis
        REAL, DIMENSION(klon, nsismx) :: eta
        REAL, DIMENSION(klon, nsnowmx) :: dzsn
        REAL, DIMENSION(klon, nsismx) :: ro
        REAL, DIMENSION(klon, nsnowmx) :: G1sn
        REAL, DIMENSION(klon, nsnowmx) :: G2sn
        REAL, DIMENSION(klon, nsnowmx) :: agsn
        REAL, DIMENSION(klon) :: IRs
        REAL, DIMENSION(klon) :: LMO
        REAL, DIMENSION(klon) :: rusn
        REAL, DIMENSION(klon) :: toic
        REAL, DIMENSION(klon) :: Bufs
        REAL, DIMENSION(klon) :: alb1, alb2, alb3

        INTEGER isl, ikl, i, isn, ierr
        CHARACTER (len = 2) :: str2
        INTEGER :: pass

        isno(:) = 0
        ispi(:) = 0
        iice(:) = 0
        IRs(:) = 0.
        LMO(:) = 0.
        eta(:, :) = 0.
        Tsis(:, :) = 0.
        isto(:, :) = 0
        ro(:, :) = 0.
        G1sn(:, :) = 0.
        G2sn(:, :) = 0.
        dzsn(:, :) = 0.
        agsn(:, :) = 0.
        rusn(:) = 0.
        toic(:) = 0.
        Bufs(:) = 0.
        alb1(:) = 0.
        alb2(:) = 0.
        alb3(:) = 0.

        !***************************************************************************
        ! Uncompress SISVAT output variables for storage

        DO  ikl = 1, klon
            i = ikl2i(ikl)
            IF (i > 0) THEN
                isno(i) = 1. * isnoSV(ikl)               ! Nb Snow/Ice Lay.
                ispi(i) = 1. * ispiSV(ikl)               ! Nb Supr.Ice Lay.
                iice(i) = 1. * iiceSV(ikl)               ! Nb      Ice Lay.

                !        IRs(i)        = IRs_SV(ikl)
                !        LMO(i)        = LMO_SV(ikl)

                DO isl = -nsol, 0                           !
                    eta(i, nsno + 1 - isl) = eta_SV(ikl, isl)            ! Soil Humidity
                    Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl)            ! Soil Temperature
                    ro(i, nsno + 1 - isl) = ro__SV(ikl, isl)            !        [kg/m3]
                END DO

                DO  isn = 1, nsno
                    isto(i, isn) = 1. * istoSV(ikl, isn)         ! Snow     History
                    ro(i, isn) = ro__SV(ikl, isn)            !        [kg/m3]
                    eta(i, isn) = eta_SV(ikl, isn)            !        [m3/m3]
                    Tsis(i, isn) = TsisSV(ikl, isn)            !            [K]
                    G1sn(i, isn) = G1snSV(ikl, isn)            ! [-]        [-]
                    G2sn(i, isn) = G2snSV(ikl, isn)            ! [-] [0.0001 m]
                    dzsn(i, isn) = dzsnSV(ikl, isn)            !            [m]
                    agsn(i, isn) = agsnSV(ikl, isn)            !          [day]
                END DO
                rusn(i) = rusnSV(ikl)                  ! Surficial Water
                toic(i) = toicSV(ikl)                  ! to ice
                alb1(i) = alb1sv(ikl)
                alb2(i) = alb2sv(ikl)
                alb3(i) = alb3sv(ikl)
                !        Bufs(i)       = BufsSV(ikl)
            END IF
        END DO

        CALL open_restartphy(fichnom)

        DO pass = 1, 2
            CALL put_field(pass, "longitude", &
                    "Longitudes de la grille physique", rlon)
            CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat)

            CALL put_field(pass, "n_snows", "number of snow/ice layers", isno)
            CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi)
            CALL put_field(pass, "n_ice", "number of ice layers", iice)
            CALL put_field(pass, "IR_soil", "Soil IR flux", IRs)
            CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO)
            CALL put_field(pass, "surf_water", "Surficial water", rusn)
            CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs)
            CALL put_field(pass, "alb_1", "albedo sw", alb1)
            CALL put_field(pass, "alb_2", "albedo nIR", alb2)
            CALL put_field(pass, "alb_3", "albedo fIR", alb3)
            CALL put_field(pass, "to_ice", "Snow passed to ice", toic)

            DO isn = 1, nsno
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "AGESNOW" // str2, &
                            "Age de la neige layer No." // str2, &
                            agsn(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO
            DO isn = 1, nsno
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "DZSNOW" // str2, &
                            "Snow/ice thickness layer No." // str2, &
                            dzsn(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO
            DO isn = 1, nsno
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "G2SNOW" // str2, &
                            "Snow Property 2, layer No." // str2, &
                            G2sn(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO
            DO isn = 1, nsno
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "G1SNOW" // str2, &
                            "Snow Property 1, layer No." // str2, &
                            G1sn(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO
            DO isn = 1, nsismx
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "ETA" // str2, &
                            "Soil/snow water content layer No." // str2, &
                            eta(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO
            DO isn = 1, nsismx   !nsno
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "RO" // str2, &
                            "Snow density layer No." // str2, &
                            ro(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO
            DO isn = 1, nsismx
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "TSS" // str2, &
                            "Soil/snow temperature layer No." // str2, &
                            Tsis(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO
            DO isn = 1, nsno
                IF (isn<=99) THEN
                    WRITE(str2, '(i2.2)') isn
                    CALL put_field(pass, "HISTORY" // str2, &
                            "Snow history layer No." // str2, &
                            isto(:, isn))
                ELSE
                    PRINT*, "Trop de couches"
                    CALL abort
                ENDIF
            ENDDO

            CALL enddef_restartphy
        ENDDO
        CALL close_restartphy

    END SUBROUTINE sisvatredem

END MODULE surf_inlandsis_mod
