! ================================================================================================================================
!  MODULE       : explicitsnow
!
!  CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
!  LICENCE      : IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF  Computes hydrologic snow processes on continental points.
!!
!!\n DESCRIPTION: This module computes hydrologic snow processes on continental points.
!!
!! RECENT CHANGES: 
!!
!! REFERENCES	: None
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/tags/ORCHIDEE_2_0/ORCHIDEE/src_sechiba/explicitsnow.f90 $
!! $Date: 2018-03-13 17:38:36 +0100 (mar. 13 mars 2018) $
!! $Revision: 5091 $
!! \n
!_ ================================================================================================================================

MODULE explicitsnow
  USE ioipsl_para
  USE constantes_soil
  USE constantes
  USE time, ONLY : one_day, dt_sechiba
  USE pft_parameters 
  USE qsat_moisture 
  USE sechiba_io_p
  USE xios_orchidee

  IMPLICIT NONE

  ! Public routines :
  PRIVATE
  PUBLIC :: explicitsnow_main, explicitsnow_initialize, explicitsnow_finalize

CONTAINS

!================================================================================================================================
!! SUBROUTINE   : explicitsnow_initialize
!!
!>\BRIEF        Read variables for explictsnow module from restart file
!!                
!! DESCRIPTION  : Read variables for explictsnow module from restart file
!!
!! \n
!_
!================================================================================================================================
  SUBROUTINE explicitsnow_initialize( kjit,     kjpindex, rest_id,    snowrho,   &
                                      snowtemp, snowdz,   snowheat,   snowgrain)

    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                           :: kjit           !! Time step number 
    INTEGER(i_std), INTENT(in)                           :: kjpindex       !! Domain size
    INTEGER(i_std),INTENT (in)                           :: rest_id        !! Restart file identifier
    
    !! 0.2 Output variables
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out)  :: snowrho        !! Snow density (Kg/m^3)
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out)  :: snowtemp       !! Snow temperature
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out)  :: snowdz         !! Snow layer thickness
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out)  :: snowheat       !! Snow heat content (J/m^2)
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out)  :: snowgrain      !! Snow grainsize (m)
    
    
    !! 1. Read from restart file
    CALL restget_p (rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit,.TRUE.,snowrho, "gather", nbp_glo, index_g)
    CALL setvar_p (snowrho, val_exp, 'Snow Density profile', xrhosmin)
    
    CALL restget_p (rest_id, 'snowtemp', nbp_glo,  nsnow, 1, kjit,.TRUE.,snowtemp, "gather", nbp_glo,  index_g)
    CALL setvar_p (snowtemp, val_exp, 'Snow Temperature profile', tp_00)
    
    CALL restget_p (rest_id, 'snowdz', nbp_glo,  nsnow, 1, kjit,.TRUE.,snowdz, "gather", nbp_glo,  index_g)
    CALL setvar_p (snowdz, val_exp, 'Snow depth profile', 0.0)
    
    CALL restget_p (rest_id, 'snowheat', nbp_glo,  nsnow, 1, kjit,.TRUE.,snowheat, "gather", nbp_glo,  index_g)
    CALL setvar_p (snowheat, val_exp, 'Snow Heat profile', 0.0)
    
    CALL restget_p (rest_id, 'snowgrain', nbp_glo,  nsnow, 1, kjit,.TRUE.,snowgrain, "gather", nbp_glo,  index_g)
    CALL setvar_p (snowgrain, val_exp, 'Snow grain profile', 0.0)
    
  END SUBROUTINE explicitsnow_initialize
  

!================================================================================================================================
!! SUBROUTINE   : explicitsnow_main
!!
!>\BRIEF        Main call for snow calculations
!!                
!! DESCRIPTION  : Main routine for calculation of the snow processes with explictsnow module.
!!
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : 
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================
  SUBROUTINE explicitsnow_main(kjpindex,    precip_rain,  precip_snow,   temp_air,    pb,       & ! in
                               u,           v,            temp_sol_new,  soilcap,     pgflux,   & ! in
                               frac_nobio,  totfrac_nobio,gtemp,                                & ! in
                               lambda_snow, cgrnd_snow,   dgrnd_snow,    contfrac,              & ! in  
                               vevapsno,    snow_age,     snow_nobio_age,snow_nobio,  snowrho,  & ! inout
                               snowgrain,   snowdz,       snowtemp,      snowheat,    snow,     & ! inout
                               temp_sol_add,          					        & ! inout
                               snowliq,     subsnownobio, grndflux,      snowmelt,    tot_melt, & ! output
                               subsinksoil )                                                      ! output
             

    !! 0. Variable and parameter declaration

    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                               :: kjpindex         !! Domain size
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: precip_rain      !! Rainfall
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: precip_snow      !! Snowfall
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: temp_air         !! Air temperature
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: pb               !! Surface pressure
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: u,v              !! Horizontal wind speed
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: temp_sol_new     !! Surface temperature
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: soilcap          !! Soil capacity
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: pgflux           !! Net energy into snowpack
    REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in)     :: frac_nobio       !! Fraction of continental ice, lakes, ...
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: totfrac_nobio    !! Total fraction of continental ice+lakes+ ...
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: gtemp            !! First soil layer temperature 
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: lambda_snow      !! Coefficient of the linear extrapolation of surface temperature 
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in)     :: cgrnd_snow       !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in)     :: dgrnd_snow       !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex), INTENT (in)           :: contfrac         !! Fraction of continent in the grid

    !! 0.2 Output fields
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out)     :: snowliq          !! Snow liquid content (m)
    REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(out)    :: subsnownobio     !! Sublimation of snow on other surface types (ice, lakes, ...)
    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: grndflux         !! Net flux into soil [W/m2]
    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: snowmelt         !! Snow melt [mm/dt_sechiba]
    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: tot_melt         !! Total melt from ice and snow
    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: subsinksoil      !! Excess of sublimation as a sink for the soil

    !! 0.3 Modified fields
    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: vevapsno         !! Snow evaporation  @tex ($kg m^{-2}$) @endtex 
    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: snow_age         !! Snow age
    REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout)  :: snow_nobio       !! Ice water balance
    REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout)  :: snow_nobio_age   !! Snow age on ice, lakes, ...
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)   :: snowrho          !! Snow density (Kg/m^3)
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)   :: snowgrain        !! Snow grainsize (m)
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)   :: snowdz           !! Snow layer thickness [m]
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)   :: snowtemp         !! Snow temperature
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)   :: snowheat         !! Snow heat content (J/m^2)
    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: snow             !! Snow mass [Kg/m^2]
    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: temp_sol_add     !! Additional energy to melt snow for snow ablation case (K)

 
    !! 0.4 Local declaration
    INTEGER(i_std)                                           :: ji, iv, jj,m,jv
    REAL(r_std),DIMENSION  (kjpindex)                        :: snow_depth_tmp
    REAL(r_std),DIMENSION  (kjpindex)                        :: snowmelt_from_maxmass
    REAL(r_std)                                              :: snowdzm1
    REAL(r_std), DIMENSION (kjpindex)                        :: thrufal          !! Water leaving snow pack [kg/m2/s]
    REAL(r_std), DIMENSION (kjpindex)                        :: d_age            !! Snow age change
    REAL(r_std), DIMENSION (kjpindex)                        :: xx               !! Temporary
    REAL(r_std), DIMENSION (kjpindex)                        :: snowmelt_tmp,snowmelt_ice,icemelt,temp_sol_new_old
    REAL(r_std), DIMENSION (kjpindex,nsnow)                  :: snowdz_old
    REAL(r_std), DIMENSION (kjpindex)                        :: ZLIQHEATXS
    REAL(r_std), DIMENSION (kjpindex)                        :: ZSNOWEVAPS, ZSNOWDZ,subsnowveg,ZSNOWMELT_XS
    REAL(r_std)                                              :: maxmass_snowdepth, snow_d1k
    REAL(r_std), DIMENSION (kjpindex,nsnow)                  :: WSNOWDZ,SMASS
    REAL(r_std), DIMENSION (kjpindex)                        :: SMASSC,snowacc,snow_remove
    INTEGER(i_std) :: locjj
    REAL(r_std)                                              :: grndflux_tmp
    REAL(r_std), DIMENSION (nsnow)                           :: snowtemp_tmp
    REAL(r_std)                                              :: s2flux_tmp,fromsoilflux
    REAL(r_std), DIMENSION (kjpindex,nsnow)                  :: pcapa_snow 
    REAL(r_std), DIMENSION (kjpindex)                        :: psnowhmass
    REAL(r_std), PARAMETER                                   :: XP00 = 1.E5
    REAL(r_std), DIMENSION (kjpindex)                        :: snowheattot_begin,snowheattot_end,del_snowheattot !![J/m2]
    REAL(r_std), DIMENSION (kjpindex,nsnow)                  :: snowliq_diag         ! Snow liquid content per snow layer [kg/m2]
    REAL(r_std), DIMENSION (kjpindex)                        :: snowliqtot_diag      ! Snow liquid content integrated over snow depth [kg/m2]
    REAL(r_std), DIMENSION (kjpindex,nsnow)                  :: snowrho_diag
    REAL(r_std), DIMENSION (kjpindex,nsnow)                  :: snowheat_diag
    REAL(r_std), DIMENSION (kjpindex,nsnow)                  :: snowgrain_diag

    !! 1. Initialization
    
    temp_sol_new_old = temp_sol_new
    snowmelt_ice(:) = zero
    icemelt(:) = zero
    tot_melt(:) = zero
    snowmelt(:) = zero
    snowmelt_from_maxmass(:) = zero

    
    !! 2. on Vegetation
    ! 2.1 Snow fall
    snowheattot_begin(:)=SUM(snowheat(:,:),2)
    CALL explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
         & snowheat,snowgrain,snowtemp,psnowhmass)
    
    ! 2.2 calculate the new snow discretization
    snow_depth_tmp(:) = SUM(snowdz(:,:),2)
    
    snowdz_old = snowdz
    
    CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
    
    ! 2.3 Snow heat redistribution
    CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
    
    ! 2.4 Diagonize water portion of the snow from snow heat content:
    DO ji=1, kjpindex
       IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
          snowtemp(ji,:) = snow3ltemp_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:))
          snowliq(ji,:) = snow3lliq_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
       ELSE
          snowliq(ji,:) = zero
          snowtemp(ji,:) = tp_00
       ENDIF
    END DO
    
    ! 2.5 snow compaction
    CALL explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
    ! Update snow heat 
    DO ji = 1, kjpindex
       snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
    ENDDO

    ! 2.6 Calculate the snow temperature profile based on heat diffusion
    CALL explicitsnow_profile (kjpindex,cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
    
    ! 2.7 Test whether snow is existed on the ground or not
    grndflux(:)=0.0
    CALL explicitsnow_gone(kjpindex,pgflux,&
         & snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
    
    ! 2.8 Calculate snow melt/refreezing processes
    CALL explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,&
         & snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air)
    
    
    ! 2.9 Snow sublimation changing snow thickness
    snow(:) = 0.0
    DO ji=1,kjpindex !domain size
       snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
    ENDDO
    
    subsnownobio(:,:) = zero
    subsinksoil(:) = zero
    
    DO ji=1, kjpindex ! domain size
       IF ( snow(ji) > snowcri ) THEN
          ! There is enough snow to split the sublimation demand between the nobio and vegetation parts.
          IF (snow_nobio(ji,iice) .GE. vevapsno(ji)) THEN
             ! There is enough snow on the ice fraction to sustain the sublimation demand.
             subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
          ELSE
             ! There is not enough snow on the ice fraction to sustain the full sublimation demand.
             ! We take all the nobio snow and set the remaining sublimation demand on the vegetation fraction.
             ! We do this because the nobio snow cannot deal with a too high sublimation demand, whereas the 
             ! vegetation fraction can, through the subsinksoil variable (see below 2.8.1), which is taken into 
             ! account in the water budget (see TWBR in hydrol.f90).
             subsnownobio(ji,iice) = snow_nobio(ji,iice)
          ENDIF
          subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
       ELSE
          ! There is a small amount of snow.
          IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
             ! The ice fraction is not too small. In this case, we make the hypothesis that the snow is mainly on the ice fraction.
             IF (snow_nobio(ji,iice) .GE. vevapsno(ji)) THEN
                ! There is enough snow on the ice fraction to sustain the sublimation demand.
                subsnownobio(ji,iice) = vevapsno(ji)
                subsnowveg(ji) = zero
             ELSE
                ! There is not enough snow on the ice fraction to sustain the sublimation demand.
                ! We take all the nobio snow and set the remaining sublimation demand on the vegetation fraction.
                ! We do this beacuse the nobio snow cannot deal with a too high sublimation demand, whereas the 
                ! vegetation fraction can, through the subsinksoil variable (see below 2.8.1), which is taken into 
                ! account in the water budget (see TWBR in hydrol.f90).
                subsnownobio(ji,iice) = snow_nobio(ji,iice)
                subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
             ENDIF
          ELSE
             ! The nobio snow fraction is too small, the vegetation fraction deals with the whole sublimation demand.
             subsnownobio(ji,iice) = zero
             subsnowveg(ji) = vevapsno(ji)
          ENDIF
       ENDIF
       
       !! 2.8.1 Check that sublimation on the vegetated fraction is possible.
       IF (subsnowveg(ji) .GT. snow(ji)) THEN
          IF( (un - totfrac_nobio(ji)).GT.min_sechiba) THEN
             subsinksoil (ji) = (subsnowveg(ji) - snow(ji))/ (un - totfrac_nobio(ji))
          END IF
          ! Sublimation is thus limited to what is available
          subsnowveg(ji) = snow(ji)
          snow(ji) = zero
          snowdz(ji,:)  =  0
          snowliq(ji,:)   =  0
          snowtemp(ji,:) = tp_00 
          vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
       ELSE
          ! Calculating the snow accumulation  
          WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:)
          SMASSC (ji)= 0.0
          DO jj=1,nsnow
             SMASS(ji,jj)   = SMASSC(ji) + WSNOWDZ(ji,jj)
             SMASSC(ji)      = SMASSC(ji) + WSNOWDZ(ji,jj)
          ENDDO
          ! Finding the layer
          locjj=1
          DO jj=1,nsnow-1
             IF ((SMASS(ji,jj) .LE. subsnowveg(ji)) .AND. (SMASS(ji,jj+1) .GE. subsnowveg(ji)) ) THEN
                locjj=jj+1
             ENDIF
          ENDDO
          
          ! Calculating the removal of snow depth
          IF (locjj .EQ. 1) THEN
             ZSNOWEVAPS(ji)  = subsnowveg(ji)/snowrho(ji,1)
             ZSNOWDZ(ji)     = snowdz(ji,1) - ZSNOWEVAPS(ji)
             snowdz(ji,1)     = MAX(0.0, ZSNOWDZ(ji))
          ELSE IF (locjj .GT. 1) THEN
             snowacc(ji)=0
             DO jj=1,locjj-1
                snowacc(ji)=snowacc(ji)+snowdz(ji,jj)*snowrho(ji,jj)
                snowdz(ji,jj)=0
             ENDDO
             ZSNOWEVAPS(ji)  = (subsnowveg(ji)-snowacc(ji))/snowrho(ji,locjj)
             ZSNOWDZ(ji)     = snowdz(ji,locjj) - ZSNOWEVAPS(ji)
             snowdz(ji,locjj)     = MAX(0.0, ZSNOWDZ(ji))
!          ELSE
!             ZSNOWEVAPS(ji)  = subsnowveg(ji)/snowrho(ji,1)
!             ZSNOWDZ(ji)     = snowdz(ji,1) - ZSNOWEVAPS(ji)
!             snowdz(ji,1)     = MAX(0.0, ZSNOWDZ(ji))
          ENDIF
          
       ENDIF
    ENDDO
    
! 2.9.1 Snowmelt_from_maxmass 

     DO ji=1,kjpindex !domain size
       snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
       IF(snow(ji).GT.maxmass_snow)THEN
           IF (snow(ji).GT.(1.2*maxmass_snow)) THEN 
              ! melt faster for points where accumulation is too high
              snow_d1k = trois * soilcap(ji) / chalfu0
           ELSE
              ! standard melting
              snow_d1k = un * soilcap(ji) / chalfu0
           ENDIF
           snowmelt_from_maxmass(ji) = MIN((snow(ji) - maxmass_snow),snow_d1k)
           ! Calculating the snow accumulation  
           WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:) ! WSNOWDZ en kg/m2
           SMASSC (ji)= 0.0
           DO jj=nsnow,1,-1
             SMASS(ji,jj)   = SMASSC(ji) + WSNOWDZ(ji,jj)
             SMASSC(ji)      = SMASSC(ji) + WSNOWDZ(ji,jj)
           ENDDO

          ! Finding the layer
          locjj = nsnow
          DO jj=nsnow,2,-1
             IF ((SMASS(ji,jj) .LE. snowmelt_from_maxmass(ji)) .AND. (SMASS(ji,jj-1) .GE. snowmelt_from_maxmass(ji)) ) THEN
                locjj=jj-1
             ENDIF
          ENDDO
          
          ! Calculating the removal of snow depth
          IF (locjj .EQ. 3) THEN
! ZSNOWMELT_XS : Epaisseur de la couche qui a disparu partiellement 
             ZSNOWMELT_XS(ji)  = snowmelt_from_maxmass(ji)/snowrho(ji,3)
             ZSNOWDZ(ji)     = snowdz(ji,3) - ZSNOWMELT_XS(ji)
             snowdz(ji,3)     = MAX(0.0, ZSNOWDZ(ji))
          ELSE IF (locjj .LT. 3) THEN
             snow_remove(ji)=0   ! masse de neige des couches qui disparaissent totalement
             DO jj=3,locjj+1,-1
                snow_remove(ji)=snow_remove(ji)+snowdz(ji,jj)*snowrho(ji,jj)
                snowdz(ji,jj)=0
             ENDDO
             ZSNOWMELT_XS(ji)  = (snowmelt_from_maxmass(ji) - snow_remove(ji))/snowrho(ji,locjj)
             ZSNOWDZ(ji)     = snowdz(ji,locjj) - ZSNOWMELT_XS(ji)
             snowdz(ji,locjj)     = MAX(0.0, ZSNOWDZ(ji))             
          ENDIF
      ENDIF
      
    ENDDO   
    
    !2.10 Calculate snow grain size using the updated thermal gradient
    
    CALL explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
    
    !2.11  Update snow heat
    ! Update the heat content (variable stored each time step)
    ! using current snow temperature and liquid water content:
    !
    ! First, make check to make sure heat content not too large
    ! (this can result due to signifigant heating of thin snowpacks):
    ! add any excess heat to ground flux:
    !
    DO ji=1,kjpindex
       DO jj=1,nsnow
          ZLIQHEATXS(ji)  = MAX(0.0, snowliq(ji,jj)*ph2o - 0.10*snowdz(ji,jj)*snowrho(ji,jj))*chalfu0/dt_sechiba
          snowliq(ji,jj) = snowliq(ji,jj) - ZLIQHEATXS(ji)*dt_sechiba/(ph2o*chalfu0)
          snowliq(ji,jj) = MAX(0.0, snowliq(ji,jj))
          grndflux(ji)   = grndflux(ji)   + ZLIQHEATXS(ji)
       ENDDO
    ENDDO
    
    snow(:) = 0.0
    DO ji=1,kjpindex !domain size
       snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
    ENDDO
        
    DO ji = 1, kjpindex
       snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
    ENDDO
    snowheattot_end(:)=SUM(snowheat(:,:),2)
    del_snowheattot(:)=snowheattot_end(:)-snowheattot_begin(:)
    
    !3. on land ice (using the default ORCHIDEE snow scheme)
    
    DO ji = 1,kjpindex
       
       !! 3.1. It is snowing
       
       snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
            & frac_nobio(ji,iice)*precip_rain(ji)
       
       !! 3.2. Sublimation - was calculated before, it cannot give us negative snow_nobio (see 2.9).
       
       snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
       
       !! 3.3. ice melt only for continental ice fraction
       
       snowmelt_tmp(ji) = zero
       IF (temp_sol_new_old(ji) .GT. tp_00) THEN
          
          !! 3.3.1 If there is snow on the ice-fraction it can melt
          
          snowmelt_tmp(ji) = frac_nobio(ji,iice)*(temp_sol_new_old(ji) - tp_00) * soilcap(ji) / chalfu0
          
          IF ( snowmelt_tmp(ji) .GT. snow_nobio(ji,iice) ) THEN
             snowmelt_tmp(ji) = MAX( 0., snow_nobio(ji,iice))
          ENDIF
          snowmelt_ice(ji) = snowmelt_ice(ji) + snowmelt_tmp(ji)
          snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp(ji)
          
       ENDIF

       !! Ice melt only if there is more than a given mass : maxmass_snow, i.e. only weight melts glaciers !
       
       IF ( snow_nobio(ji,iice) .GE. maxmass_snow ) THEN
          icemelt(ji) = snow_nobio(ji,iice) - maxmass_snow
          snow_nobio(ji,iice) = maxmass_snow
       ENDIF
       
    END DO
    
    
    !! 4. On other surface types - not done yet
    
    IF ( nnobio .GT. 1 ) THEN
       WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
       WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES'
       CALL ipslerr_p(3,'explicitsnow_main','Cannot treat snow on these surface types.','','')
    ENDIF
    
    !! 5. Computes snow age on land and land ice (for albedo)
    
    DO ji = 1, kjpindex
       
       !! 5.1. Snow age on land 
       
       IF (snow(ji) .LE. zero) THEN
          snow_age(ji) = zero
       ELSE
          snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dt_sechiba/one_day) &
               & * EXP(-precip_snow(ji) / snow_trans)
       ENDIF
       
       !! 5.2. Snow age on land ice
       
       !! Age of snow on ice: a little bit different because in cold regions, we really
       !! cannot negect the effect of cold temperatures on snow metamorphism any more.
       
       IF (snow_nobio(ji,iice) .LE. zero) THEN
          snow_nobio_age(ji,iice) = zero
       ELSE
          
          d_age(ji) = ( snow_nobio_age(ji,iice) + &
               &  (un - snow_nobio_age(ji,iice)/max_snow_age) * dt_sechiba/one_day ) * &
               &  EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
          IF (d_age(ji) .GT. 0. ) THEN
             xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
             xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
             d_age(ji) = d_age(ji) / (un+xx(ji))
          ENDIF
          snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
          
       ENDIF
       
    ENDDO
    
    
    !! 6. Check the snow on land 
    DO ji=1,kjpindex
       IF (snow(ji) .EQ. 0) THEN
          snowrho(ji,:)=50.0
          snowgrain(ji,:)=0.0
          snowdz(ji,:)=0.0
          snowliq(ji,:)=0.0
       ENDIF
    ENDDO
    
    
    ! Snow melt only if there is more than a given mass : maxmass_snow
    ! Here I suggest to remove the snow based on a certain threshold of snow
    ! depth instead of snow mass because it is quite difficult for
    ! explictsnow module to calculate other snow properties following the
    ! removal of snow mass
    ! to define the threshold of snow depth based on old snow density (330
    ! kg/m3)
    !       maxmass_snowdepth=maxmass_snow/sn_dens 
!    snowmelt_from_maxmass(:) = 0.0
    
    !! 7. compute total melt 
    
    DO ji=1,kjpindex 
       tot_melt(ji) = icemelt(ji) + snowmelt(ji) + snowmelt_ice(ji) + snowmelt_from_maxmass(ji)
    ENDDO
    IF (printlev>=3) WRITE(numout,*) 'explicitsnow_main done'
    

    ! Write output variables

    ! Add XIOS default value where no snow
    DO ji=1,kjpindex 
       IF (snow(ji) .GT. zero) THEN
          ! Output for snowliq and snowliqtot change unit from m to kg/m2 and are multiplied by contfrac 
          ! to take into acount the portion form land fraction only
          snowliq_diag(ji,:) = snowliq(ji,:) * 1000 * contfrac(ji)
          snowliqtot_diag(ji) = SUM(snowliq_diag(ji,:)) * 1000 * contfrac(ji)
          snowrho_diag(ji,:) = snowrho(ji,:)
          snowheat_diag(ji,:) = snowheat(ji,:)
          snowgrain_diag(ji,:) = snowgrain(ji,:)
       ELSE
          snowliq_diag(ji,:) = xios_default_val
          snowliqtot_diag(ji) = xios_default_val
          snowrho_diag(ji,:) = xios_default_val
          snowheat_diag(ji,:) = xios_default_val
          snowgrain_diag(ji,:) = xios_default_val
       END IF
    END DO

    CALL xios_orchidee_send_field("snowliq",snowliq_diag)
    CALL xios_orchidee_send_field("snowliqtot", snowliqtot_diag)
    CALL xios_orchidee_send_field("snowrho",snowrho_diag)
    CALL xios_orchidee_send_field("snowheat",snowheat_diag)
    CALL xios_orchidee_send_field("snowgrain",snowgrain_diag)
    CALL xios_orchidee_send_field("snowmelt_from_maxmass",snowmelt_from_maxmass)
    CALL xios_orchidee_send_field("soilcap",soilcap)
    CALL xios_orchidee_send_field("del_snowheattot",del_snowheattot)

  END SUBROUTINE explicitsnow_main
  

!================================================================================================================================
!! SUBROUTINE   : explicitsnow_finalize
!!
!>\BRIEF        Write variables for explictsnow module to restart file
!!                
!! DESCRIPTION  : Write variables for explictsnow module to restart file
!!
!! \n
!_
!================================================================================================================================
  SUBROUTINE explicitsnow_finalize ( kjit,     kjpindex, rest_id,    snowrho,   &
       snowtemp, snowdz,   snowheat,   snowgrain)
    
    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                           :: kjit           !! Time step number 
    INTEGER(i_std), INTENT(in)                           :: kjpindex       !! Domain size
    INTEGER(i_std),INTENT (in)                           :: rest_id        !! Restart file identifier
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowrho        !! Snow density (Kg/m^3)
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowtemp       !! Snow temperature
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowdz         !! Snow layer thickness
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowheat       !! Snow heat content (J/m^2)
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowgrain      !! Snow grainsize (m)
    
    
    !! 1. Write to restart file
    CALL restput_p(rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit, snowrho, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit, snowtemp, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit, snowdz, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit, snowheat, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit, snowgrain, 'scatter', nbp_glo, index_g)
    
  END SUBROUTINE explicitsnow_finalize
  

!================================================================================================================================
!! SUBROUTINE   : explicitsnow_grain
!!
!>\BRIEF        Compute evolution of snow grain size
!!                
!! DESCRIPTION  : 
!!
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : R. Jordan (1991)
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================


 SUBROUTINE explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)

      !! 0.1 Input variables
      INTEGER(i_std),INTENT(in)                                  :: kjpindex         !! Domain size
      REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)           :: snowliq          !! Liquid water content
      REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)           :: snowdz           !! Snow depth (m)
      REAL(r_std),DIMENSION(kjpindex),INTENT(in)                 :: gtemp            !! First soil layer temperature
      REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)           :: snowtemp         !! Snow temperature
      REAL(r_std),DIMENSION (kjpindex),INTENT(in)                :: pb               !! Surface pressure (hpa) 

      !! 0.2 Output variables

      !! 0.3 Modified variables

      REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)        :: snowgrain        !! Snow grain size (m)

      !! 0.4 Local variables
      REAL(r_std),DIMENSION(kjpindex,nsnow)                      :: zsnowdz,zdz,ztheta
      REAL(r_std),DIMENSION(kjpindex,0:nsnow)                    :: ztemp,zdiff,ztgrad,zthetaa,zfrac,&
                                                                    zexpo,zckt_liq,zckt_ice,zckt
      REAL(r_std),DIMENSION(kjpindex,nsnow)                      :: zrhomin,zgrainmin
      INTEGER(i_std) :: ji,jj

      !! 0.5 Local parameters
      REAL(r_std), PARAMETER                                     :: ztheta_crit = 0.02     !! m3 m-3
      REAL(r_std), PARAMETER                                     :: zc1_ice     = 8.047E+9 !! kg m-3 K
      REAL(r_std), PARAMETER                                     :: zc1_liq     = 5.726E+8 !! kg m-3 K
      REAL(r_std), PARAMETER                                     :: zdeos       = 0.92E-4  !! effective diffusion
                                                                                           !! coef for water vapor in snow
                                                                                           !! at 0C and 1000 mb (m2 s-1)
      REAL(r_std), PARAMETER                                     :: zg1         = 5.0E-7   !! m4 kg-1
      REAL(r_std), PARAMETER                                     :: zg2         = 4.0E-12  !! m2 s-1
      REAL(r_std), PARAMETER                                     :: ztheta_w      = 0.05   !! m3 m-3
      REAL(r_std), PARAMETER                                     :: ztheta_crit_w = 0.14   !! m3 m-3
      REAL(r_std), PARAMETER                                     :: zdzmin        = 0.01   !! m : minimum thickness
                                                                                           !! for thermal gradient evaluation:
                                                                                           !! to prevent excessive gradients
                                                                                           !! for vanishingly thin snowpacks.
      REAL(r_std), PARAMETER                                     :: xp00=1.E5 
      !! 1. initialize
 
      DO ji=1,kjpindex 


         zsnowdz(ji,:)  = MAX(xsnowdmin/nsnow, snowdz(ji,:))

         DO jj=1,nsnow-1
            zdz(ji,jj)      = zsnowdz(ji,jj) + zsnowdz(ji,jj+1)
         ENDDO
            zdz(ji,nsnow)     = zsnowdz(ji,nsnow)

         ! compute interface average volumetric water content (m3 m-3):
         ! first, layer avg VWC:
         !
          ztheta(ji,:) = snowliq(ji,:)/MAX(xsnowdmin, zsnowdz(ji,:))

         ! at interfaces:
          zthetaa(ji,0)      = ztheta(ji,1)
          DO jj=1,nsnow-1
             zthetaa(ji,jj)  = (zsnowdz(ji,jj)  *ztheta(ji,jj)   +             &
                                   zsnowdz(ji,jj+1)*ztheta(ji,jj+1))/zdz(ji,jj)
          ENDDO
          zthetaa(ji,nsnow) = ztheta(ji,nsnow)
          ! compute interface average temperatures (K):
          ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
          !
          ztemp(ji,0)      = snowtemp(ji,1)
          DO jj=1,nsnow-1
             ztemp(ji,jj)  = (zsnowdz(ji,jj)  *snowtemp(ji,jj)   +             &
                                 zsnowdz(ji,jj+1)*snowtemp(ji,jj+1))/zdz(ji,jj)
          ENDDO
          ztemp(ji,nsnow) = snowtemp(ji,nsnow)
!
          ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
          ! compute variation of saturation vapor pressure with temperature
          ! for solid and liquid phases:
          ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
          zexpo(ji,:)    = chalsu0/(xrv*ztemp(ji,:))
          zckt_ice(ji,:) = (zc1_ice/ztemp(ji,:)**2)*(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
          !
          zexpo(ji,:)    = chalev0/(xrv*ztemp(ji,:))
          zckt_liq(ji,:) = (zc1_liq/ztemp(ji,:)**2)*(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
          !
          ! compute the weighted ice/liquid total variation (N m-2 K):
          !
          zfrac(ji,:)    = MIN(1.0, zthetaa(ji,:)/ztheta_crit)
          zckt(ji,:)     = zfrac(ji,:)*zckt_liq(ji,:) + (1.0 - zfrac(ji,:))*zckt_ice(ji,:)

          ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
          ! Compute effective diffusion coefficient (m2 s-1):
          ! -diffudivity relative to a reference diffusion at 1000 mb and freezing point
          !  multiplied by phase energy coefficient
          ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
          !
          DO jj=0,nsnow
             zdiff(ji,jj) = zdeos*(xp00/(pb(ji)*100.))*((ztemp(ji,jj)/tp_00)**6)*zckt(ji,jj)
          ENDDO 

          ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
          ! Temperature gradient (K m-1):

          ztgrad(ji,0)      = 0.0 ! uppermost layer-mean and surface T's are assumed to be equal
          DO jj=1,nsnow-1
             ztgrad(ji,jj)  = 2*(snowtemp(ji,jj)-snowtemp(ji,jj+1))/MAX(zdzmin, zdz(ji,jj))
          ENDDO
          !
          ! assume at base of snow, temperature is in equilibrium with soil
          ! (but obviously must be at or below freezing point!)
          !
          ztgrad(ji,nsnow) = 2*(snowtemp(ji,nsnow) - MIN(tp_00, gtemp(ji)))/MAX(zdzmin, zdz(ji,nsnow))
          ! prognostic grain size (m) equation:
          !-------------------------------------------------------------------
          ! first compute the minimum grain size (m):
          !
          zrhomin(ji,:)     = xrhosmin
          zgrainmin(ji,:)   = snow3lgrain_1d(zrhomin(ji,:))

          ! dry snow:
          ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
          !
          DO jj=1,nsnow

             IF(ztheta(ji,jj) == 0.0) THEN

              ! only allow growth due to vapor flux INTO layer: 
              ! aab add sublimation(only condensation) as upper BC...?

                snowgrain(ji,jj)  = snowgrain(ji,jj) +                                      &
                                       (dt_sechiba*zg1/MAX(zgrainmin(ji,jj),snowgrain(ji,jj)))*      &
                                       ( zdiff(ji,jj-1)*MAX(0.0,ztgrad(ji,jj-1)) -                &
                                       zdiff(ji,jj)  *MIN(0.0,ztgrad(ji,jj)) )
             ELSE
 
              ! wet snow
              ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
              !
                snowgrain(ji,jj)  = snowgrain(ji,jj) +                                      &
                                       (dt_sechiba*zg2/MAX(zgrainmin(ji,jj),snowgrain(ji,jj)))*      &
                                       MIN(ztheta_crit_w, ztheta(ji,jj) + ztheta_w)
             END IF

          ENDDO


      ENDDO


    END SUBROUTINE explicitsnow_grain

!!
!================================================================================================================================
!! SUBROUTINE   : explicitsnow_compactn
!!
!>\BRIEF        Compute Compaction/Settling
!!                
!! DESCRIPTION  : 
!!     Snow compaction due to overburden and settling.
!!     Mass is unchanged: layer thickness is reduced
!!     in proportion to density increases. Method
!!     of Anderson (1976): see Loth and Graf, 1993,
!!     J. of Geophys. Res., 98, 10,451-10,464.
!!
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : Loth and Graf (1993), Mellor (1964) and Anderson (1976)
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================


   SUBROUTINE explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)

         !! 0.1 Input variables

         INTEGER(i_std),INTENT(in)                                 :: kjpindex         !! Domain size
         REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)          :: snowtemp         !! Snow temperature

         !! 0.2 Output variables

         !! 0.3 Modified variables

         REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)       :: snowrho          !! Snow density (Kg/m^3)
         REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)       :: snowdz           !! Snow depth

         !! 0.4 Local variables

         REAL(r_std),DIMENSION(kjpindex,nsnow)                     :: zwsnowdz,zsmass,zsnowrho2,zviscocity,zsettle 
         REAL(r_std),DIMENSION(kjpindex)                           :: zsmassc          !! cummulative snow mass (kg/m2)
         REAL(r_std),DIMENSION(kjpindex)                           :: snowdepth_crit
         INTEGER(i_std)                                            :: ji,jj

        !! 1. initialize

        zsnowrho2  = snowrho
        zsettle(:,:)    = ZSNOWCMPCT_ACM
        zviscocity(:,:) = ZSNOWCMPCT_V0

        !! 2. Calculating Cumulative snow mass (kg/m2):

        DO ji=1, kjpindex


            IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN 

              zwsnowdz(ji,:)= snowdz(ji,:)*snowrho(ji,:)

              zsmassc (ji)= 0.0

              DO jj=1,nsnow
                 zsmass(ji,jj)   = zsmassc(ji) + zwsnowdz(ji,jj)
                 zsmassc(ji)      = zsmassc(ji) + zwsnowdz(ji,jj)
              ENDDO


              !! 3. Computing compaction/Settling
              ! ----------------------
              ! Compaction/settling if density below upper limit
              ! (compaction is generally quite small above ~ 500 kg m-3):
              !
              DO jj=1,nsnow
                 IF (snowrho(ji,jj) .LT. xrhosmax) THEN
            
              !
              ! First calculate settling due to freshly fallen snow: (NOTE:bug here for the snow temperature profile)
              !
             
              zsettle(ji,jj)     = ZSNOWCMPCT_ACM*EXP(                                      &
                                     -ZSNOWCMPCT_BCM*(tp_00-MIN(tp_00,snowtemp(ji,jj)))                      &
                                     -ZSNOWCMPCT_CCM*MAX(0.0,                                  &
                                      snowrho(ji,jj)-ZSNOWCMPCT_RHOD))
              !
              ! Snow viscocity:
              !
              zviscocity(ji,jj)   = ZSNOWCMPCT_V0*EXP( ZSNOWCMPCT_VT*(tp_00-MIN(tp_00,snowtemp(ji,jj))) +        &
                                  ZSNOWCMPCT_VR*snowrho(ji,jj) )

              ! Calculate snow density: compaction from weight/over-burden
              ! Anderson 1976 method:
              zsnowrho2(ji,jj)    = snowrho(ji,jj) + snowrho(ji,jj)*dt_sechiba*(          &
                                   (cte_grav*zsmass(ji,jj)/zviscocity(ji,jj))                 &
                                    + zsettle(ji,jj) )
              ! Conserve mass by decreasing grid thicknesses in response
              ! to density increases
              !
              snowdz(ji,jj)  = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj))

                 ENDIF
              ENDDO

              ! Update density (kg m-3):
              snowrho(ji,:) = zsnowrho2(ji,:)

            ENDIF

        ENDDO

      END SUBROUTINE explicitsnow_compactn

!!
!================================================================================================================================
!! SUBROUTINE   : explicitsnow_transf
!!
!>\BRIEF        Computing snow mass and heat redistribution due to grid thickness configuration resetting
!!                
!! DESCRIPTION  : Snow mass and heat redistibution due to grid thickness
!!                configuration resetting. Total mass and heat content
!!                of the overall snowpack unchanged/conserved within this routine. 
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : 
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================

   SUBROUTINE explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)

     !! 0.1 Input variables

     INTEGER(i_std),INTENT(in)                                           :: kjpindex         !! Domain size
     REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)                    :: snowdz_old       !! Snow depth at the previous time step 

     !! 0.2 Output variables

     !! 0.3 Modified variables

     REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)                 :: snowrho          !! Snow density (Kg/m^3)
     REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)                 :: snowgrain        !! Snow grain size (m)
     REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)                 :: snowdz           !! Snow depth (m)
     REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)                 :: snowheat         !! Snow heat content/enthalpy (J/m2)
                                                                        
     !! 0.4 Local varibles                                              
      
     REAL(r_std),DIMENSION(kjpindex,0:nsnow)                             :: zsnowzo
     REAL(r_std),DIMENSION(kjpindex,0:nsnow)                             :: zsnowzn 
     REAL(r_std),DIMENSION(kjpindex,nsnow)                               :: zsnowddz 
     REAL(r_std),DIMENSION(kjpindex,nsnow)                               :: zdelta
     REAL(r_std),DIMENSION(kjpindex,nsnow)                               :: zsnowrhon,zsnowheatn,zsnowgrainn
     REAL(r_std),DIMENSION(kjpindex)                                     :: zsnowmix_delta
     REAL(r_std),DIMENSION(kjpindex)                                     :: zsumheat, zsumswe, zsumgrain
     INTEGER(i_std),DIMENSION(nsnow,2)                                   :: locflag
     REAL(r_std)                                                         :: psnow
     INTEGER(i_std)                                                      :: ji,jj,jjj

     ! Initialization
     zsumheat(:)       = 0.0
     zsumswe(:)        = 0.0
     zsumgrain(:)      = 0.0
     zsnowmix_delta(:) = 0.0
     locflag(:,:)        = 0

     DO ji=1, kjpindex 
 

        psnow = SUM(snowdz(ji,:))

        IF (psnow .GE. xsnowcritd .AND. snowdz_old(ji,1) .NE. 0 .AND. snowdz_old(ji,2) .NE. 0 .AND. snowdz_old(ji,3) .NE. 0) THEN
          !
          zsnowzo(ji,0)     = 0.
          zsnowzn(ji,0)     = 0.
          zsnowzo(ji,1)     = snowdz_old(ji,1)
          zsnowzn(ji,1)     = snowdz(ji,1)
          !
          DO jj=2,nsnow
             zsnowzo(ji,jj) = zsnowzo(ji,jj-1) + snowdz_old(ji,jj)
             zsnowzn(ji,jj) = zsnowzn(ji,jj-1) + snowdz(ji,jj)
          ENDDO

          DO jj=1,nsnow

             IF (jj .EQ. 1) THEN 
                locflag(jj,1)=1 
             ELSE       
                DO jjj=nsnow,1,-1

                   !upper bound of the snow layer
                   IF (zsnowzn(ji,jj-1) .LE. zsnowzo(ji,jjj)) THEN
                      locflag(jj,1) = jjj
                   ENDIF
                ENDDO
	     ENDIF

             IF (jj .EQ. nsnow) THEN 
                locflag(jj,2)=nsnow 
             ELSE       
                DO jjj=nsnow,1,-1

                   !lower bound of the snow layer
                   IF (zsnowzn(ji,jj)   .LE. zsnowzo(ji,jjj)) THEN
                      locflag(jj,2) = jjj
                   ENDIF
                ENDDO
             ENDIF

             !to interpolate
             ! when heavy snow occurred
             IF (locflag(jj,1) .EQ. locflag(jj,2)) THEN 

                zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1))

                zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1))*snowdz(ji,jj)/snowdz_old(ji,locflag(jj,1)) 

                zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) 

             ELSE 

                !snow density
                zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1)) * &
                                      (zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
                                      snowrho(ji,locflag(jj,2)) * &
                                      (zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))

                DO jjj=locflag(jj,1),locflag(jj,2)-1
                   zsnowrhon(ji,jj)=zsnowrhon(ji,jj) + &
                   (jjj-locflag(jj,1))*snowrho(ji,jjj)*snowdz_old(ji,jjj)  
                ENDDO

                zsnowrhon(ji,jj) = zsnowrhon(ji,jj) / snowdz(ji,jj)


                !snow heat
                IF (snowdz_old(ji,locflag(jj,1)) .GT. 0.0) THEN

                   zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1)) * &
                                 (zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1))/snowdz_old(ji,locflag(jj,1)) + &
                                       snowheat(ji,locflag(jj,2)) * &
                                 (zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))/snowdz_old(ji,locflag(jj,2))
	        ELSE
                   zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,2))/snowdz_old(ji,locflag(jj,2)) * &
                        (zsnowzn(ji,locflag(jj,1))-zsnowzo(ji,locflag(jj,1)))
		ENDIF

                DO jjj=locflag(jj,1),locflag(jj,2)-1
                   zsnowheatn(ji,jj)=zsnowheatn(ji,jj) + &
                                        (jjj-locflag(jj,1))*snowheat(ji,jjj)
                ENDDO


                !snow grain
                zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) * &
                                      (zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
                                      snowgrain(ji,locflag(jj,2)) * &
                                      (zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))

                DO jjj=locflag(jj,1),locflag(jj,2)-1
                   zsnowgrainn(ji,jj)=zsnowgrainn(ji,jj) + &
                   (jjj-locflag(jj,1))*snowgrain(ji,jjj)*snowdz_old(ji,jjj)
                ENDDO

                zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) / snowdz(ji,jj)

             ENDIF
                        
          ENDDO

          snowrho(ji,:)    = zsnowrhon(ji,:)
          snowheat(ji,:)   = zsnowheatn(ji,:)
          snowgrain(ji,:)  = zsnowgrainn(ji,:)

        ENDIF

          ! Vanishing or very thin snowpack check:
          ! -----------------------------------------
          !
          ! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous):
          ! this avoids problems related to heat and mass exchange for
          ! thin layers during heavy snowfall or signifigant melt: one
          ! new/old layer can exceed the thickness of several old/new layers.
          ! Therefore, mix (conservative):
          !
          ! modified by Tao Wang
        IF (psnow > 0 .AND. (psnow < xsnowcritd .OR. snowdz_old(ji,1) &
            & .eq. 0 .OR. snowdz_old(ji,2) .eq. 0 .OR. snowdz_old(ji,3) .eq. 0)) THEN
                zsumheat(ji) = SUM(snowheat(ji,:))
                zsumswe(ji)  = SUM(snowrho(ji,:)*snowdz_old(ji,:))
                zsumgrain(ji)= SUM(snowgrain(ji,:)*snowdz_old(ji,:))
                zsnowmix_delta(ji) = 1.0
                DO jj=1,nsnow
                   zsnowheatn(ji,jj)  = zsnowmix_delta(ji)*(zsumheat(ji)/nsnow)
                   snowdz(ji,jj)    = zsnowmix_delta(ji)*(psnow/nsnow) 
                   zsnowrhon(ji,jj)   = zsnowmix_delta(ji)*(zsumswe(ji)/psnow)
                   zsnowgrainn(ji,jj) = zsnowmix_delta(ji)*(zsumgrain(ji)/psnow)
                ENDDO
          ! Update mass (density and thickness), heat and grain size:
          ! ------------------------------------------------------------
          !
          snowrho(ji,:)    = zsnowrhon(ji,:)
          snowheat(ji,:)   = zsnowheatn(ji,:)
          snowgrain(ji,:)  = zsnowgrainn(ji,:)

        ENDIF

     ENDDO


   END SUBROUTINE explicitsnow_transf

  
!!
!================================================================================================================================
!! SUBROUTINE   : explicitsnow_fall
!!
!>\BRIEF    Computes snowfall    
!!                
!! DESCRIPTION  : Computes snowfall    
!routine. 
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : 
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================
 
   SUBROUTINE explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
                        & snowheat,snowgrain,snowtemp,psnowhmass)

       !! 0.1 Input variables

       INTEGER(i_std),INTENT(in)                              :: kjpindex            !! Domain size
       REAL(r_std),DIMENSION(kjpindex),INTENT(in)             :: precip_snow         !! Snow rate (SWE) (kg/m2 per dt_sechiba)
       REAL(r_std),DIMENSION(kjpindex),INTENT(in)             :: temp_air            !! Air temperature
       REAL(r_std),DIMENSION(kjpindex),INTENT(in)             :: u,v                 !! Horizontal wind speed
       REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)       :: snowtemp            !! Snow temperature
       REAL(r_std),DIMENSION(kjpindex),INTENT(in)             :: totfrac_nobio
       !! 0.2 Output variables

       REAL(r_std), DIMENSION(kjpindex),INTENT(out)           :: psnowhmass          !! Heat content of snowfall (J/m2)

       !! 0.3 Modified variables
 
       REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)    :: snowrho             !! Snow density profile (kg/m3)
       REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)    :: snowdz              !! Snow layer thickness profile (m)
       REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)    :: snowheat            !! Snow heat content/enthalpy (J/m2)
       REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout)    :: snowgrain           !! Snow grain size (m)

       !! 0.4 Local variables

       REAL(r_std), DIMENSION(kjpindex)                       :: rhosnew             !! Snowfall density
       REAL(r_std), DIMENSION(kjpindex)                       :: dsnowfall           !! Snowfall thickness (m)
       REAL(r_std), DIMENSION(kjpindex,nsnow)                 :: snowdz_old
       REAL(r_std), DIMENSION(kjpindex)                       :: snow_depth,snow_depth_old,newgrain
       REAL(r_std)                                            :: snowfall_delta,speed
       INTEGER(i_std)                                         :: ji,jj

       !! 1. initialize the variables

       snowdz_old = snowdz
       DO ji=1,kjpindex
                 snow_depth(ji) = SUM(snowdz(ji,:))
       ENDDO

       snow_depth_old = snow_depth
 
       snowfall_delta = 0.0 

       !! 2. incorporate snowfall into snowpack 
       DO ji = 1, kjpindex
   
       speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
        
          ! new snow fall on snowpack 
          ! NOTE: when the surface temperature is zero, it means that snowfall has no
          ! heat content but it can be used for increasing the thickness and changing the density (maybe it is a bug)
          psnowhmass(ji) = 0.0
          IF ( (precip_snow(ji) .GT. 0.0) ) THEN

             !calculate

             psnowhmass(ji) = precip_snow(ji)*(un-totfrac_nobio(ji))* &
                                (xci*(snowtemp(ji,1)-tp_00)-chalfu0)
            
             ! Snowfall density: Following CROCUS (Pahaut 1976)
             !
             rhosnew(ji)   = MAX(xrhosmin, snowfall_a_sn + snowfall_b_sn*(temp_air(ji)-tp_00)+         &
                   snowfall_c_sn*SQRT(speed))

             ! Augment total pack depth:
             ! Note that dsnowfall(ji) can be equal to zero if totfrac_nobio(ji)=1
             dsnowfall(ji) = (precip_snow(ji)*(un-totfrac_nobio(ji)))/rhosnew(ji) !snowfall thickness (m)

             snow_depth(ji) = snow_depth(ji) + dsnowfall(ji) 


             ! Fresh snowfall changes the snowpack density and liquid content in uppermost layer 

             IF (dsnowfall(ji) .GT. zero) THEN
               snowrho(ji,1) = (snowdz(ji,1)*snowrho(ji,1) + dsnowfall(ji)*rhosnew(ji))/     &
                                (snowdz(ji,1)+dsnowfall(ji))

               snowdz(ji,1) = snowdz(ji,1) + dsnowfall(ji)

               ! Add energy of snowfall to snowpack:
               ! Update heat content (J/m2) (therefore the snow temperature
               ! and liquid content):
               !
               snowheat(ji,1)  = snowheat(ji,1) + psnowhmass(ji)
               !
               ! Incorporate snowfall grain size:
               !
               newgrain(ji)    = MIN(dgrain_new_max, snow3lgrain_0d(rhosnew(ji)))
               
               snowgrain(ji,1) = (snowdz_old(ji,1)*snowgrain(ji,1) + dsnowfall(ji)*newgrain(ji))/ &
                    snowdz(ji,1)
            ENDIF
          ELSE
            dsnowfall(ji) = 0.
          ENDIF

          ! new snow fall on snow free surface. 
          ! we use the linearization for the new snow fall on snow-free ground 

          IF ( (dsnowfall(ji) .GT. zero) .AND. (snow_depth_old(ji) .EQ. zero) ) THEN

             snowfall_delta = 1.0
              
             DO jj=1,nsnow

                snowdz(ji,jj)   = snowfall_delta*(dsnowfall(ji)/nsnow) + &
                                  (1.0-snowfall_delta)*snowdz(ji,jj)

                snowheat(ji,jj)   = snowfall_delta*(psnowhmass(ji)/nsnow) + &
                                  (1.0-snowfall_delta)*snowheat(ji,jj)

                snowrho(ji,jj)    = snowfall_delta*rhosnew(ji)            + &
                                  (1.0-snowfall_delta)*snowrho(ji,jj)

                snowgrain(ji,jj)=   snowfall_delta*newgrain(ji)           + &
                                  (1.0-snowfall_delta)*snowgrain(ji,jj)

             ENDDO


          ENDIF


       ENDDO 

     END SUBROUTINE explicitsnow_fall

!!
!================================================================================================================================
!! SUBROUTINE   : explicitsnow_gone
!!
!>\BRIEF        Check whether snow is gone 
!!                
!! DESCRIPTION  : If so, set thickness (and therefore mass and heat) and liquid
!!                content to zero, and adjust fluxes of water, evaporation and
!!                heat into underlying surface. 
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : 
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================

   SUBROUTINE explicitsnow_gone(kjpindex,pgflux,&
                        snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)

     !! 0.1 Input variables

     INTEGER(i_std), INTENT(in)                                 :: kjpindex     !! Domain size
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)              :: pgflux       !! Net energy into snow pack(w/m2)
     REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in)          :: snowheat     !! Snow heat content (J/m^2)

     !! 0.2 Output variables

     !! 0.3 Modified variables

     REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)     :: snowtemp     !! Snow temperature
     REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)     :: snowdz       !! Snow depth
     REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout)      :: snowrho      !! Snow density (Kg/m^3)
     REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout)      :: snowliq      !! Liquid water content
     REAL(r_std),DIMENSION(kjpindex), INTENT(inout)             :: grndflux     !! Soil/snow interface heat flux (W/m2)
     REAL(r_std),DIMENSION(kjpindex),INTENT(inout)              :: snowmelt     !! Snow melt
     REAL(r_std),DIMENSION(kjpindex)                            :: thrufal      !! Water leaving snowpack(kg/m2/s)

     !! 0.4 Local variables

     INTEGER(i_std)                                             :: ji,jj
     REAL(r_std),DIMENSION(kjpindex)                            :: snowgone_delta
     REAL(r_std),DIMENSION (kjpindex)                           :: totsnowheat  !!snow heat content at each layer 
     REAL(r_std),DIMENSION(kjpindex)                            :: snowdepth_crit

     ! first caculate total snowpack snow heat content
     snowgone_delta(:) = un
     thrufal(:)=0.0
     snowmelt(:)=0
     totsnowheat(:)  = SUM(snowheat(:,:),2) 
     
     DO ji = 1, kjpindex 

           IF ( (SUM(snowdz(ji,:)) .GT. min_sechiba)) THEN

              IF( pgflux(ji) >= (-totsnowheat(ji)/dt_sechiba) ) THEN

                 grndflux(ji) = pgflux(ji) + (totsnowheat(ji)/dt_sechiba)

                 thrufal(ji)=SUM(snowrho(ji,:)*snowdz(ji,:))

                 snowgone_delta(ji) = 0.0       

                 snowmelt(ji) = snowmelt(ji)+thrufal(ji)

              ENDIF

              ! update of snow state (either still present or not)

              DO jj=1,nsnow
                 snowdz(ji,jj)  =   snowdz(ji,jj) *snowgone_delta(ji)
                 snowliq(ji,jj)   =   snowliq(ji,jj) *snowgone_delta(ji)
                 snowtemp(ji,jj) = (1.0-snowgone_delta(ji))*tp_00 + snowtemp(ji,jj)*snowgone_delta(ji)
              ENDDO

           ELSE
	      snowdz(ji,:)=0.
              snowliq(ji,:)=0.
	      snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
           ENDIF
     ENDDO 

   END SUBROUTINE explicitsnow_gone

!================================================================================================================================
!! SUBROUTINE   : explicitsnow_melt_refrz
!!
!>\BRIEF        Computes snow melt and refreezing processes within snowpack
!!                
!! DESCRIPTION  : 
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : 
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================

   SUBROUTINE explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,&
                     snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air)

      !! 0.1 Input variables

      INTEGER(i_std), INTENT (in)                             :: kjpindex     !! Domain size
      REAL(r_std),DIMENSION (kjpindex,nsnow)                  :: pcapa_snow   !! Heat capacity for snow
      REAL(r_std),DIMENSION (kjpindex), INTENT(in)            :: precip_rain  !! Rainfall      
      REAL(r_std),DIMENSION (kjpindex), INTENT(in)            :: temp_air     !! Air temperature
      REAL(r_std),DIMENSION (kjpindex),INTENT(in)             :: pgflux       !! Net energy into snowpack(w/m2)
      REAL(r_std),DIMENSION (kjpindex),INTENT(in)             :: soilcap      !! Soil heat capacity

      !! 0.2 Output variables

      !! 0.3 Modified variables

      REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)  :: snowtemp     !! Snow temperature 
      REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)  :: snowdz       !! Snow depth
      REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout)   :: snowrho      !! Snow layer density (Kg/m^3)
      REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout)   :: snowliq      !! Liquid water content
      REAL(r_std),DIMENSION(kjpindex),INTENT(inout)           :: grndflux     !! Net energy input to soil 
      REAL(r_std),DIMENSION (kjpindex), INTENT(inout)         :: snowmelt     !! Snowmelt

      !! 0.4 Local variables

      REAL(r_std),DIMENSION (kjpindex)                        :: meltxs       !! Residual snowmelt energy applied to underlying soil
      REAL(r_std)                                             :: enerin,melttot,hrain 
      REAL(r_std),DIMENSION (nsnow)                           :: zsnowlwe
      REAL(r_std),DIMENSION (nsnow)                           :: flowliq
      REAL(r_std),DIMENSION (kjpindex)                        :: snowmass
      REAL(r_std),DIMENSION (nsnow)                           :: zphase       !! Phase change (from ice to water) (J/m2)
      REAL(r_std),DIMENSION (nsnow)                           :: zphase2      !! Phase change (from water to ice)
      REAL(r_std),DIMENSION (nsnow)                           :: zphase3      !! Phase change related with net energy input to snowpack
      REAL(r_std),DIMENSION (nsnow)                           :: zsnowdz      !! Snow layer depth
      REAL(r_std),DIMENSION (nsnow)                           :: zsnowmelt    !! Snow melt (liquid water) (m)
      REAL(r_std),DIMENSION (nsnow)                           :: zsnowtemp
      REAL(r_std),DIMENSION (nsnow)                           :: zmeltxs      !! Excess melt
      REAL(r_std),DIMENSION (nsnow)                           :: zwholdmax    !! Maximum liquid water holding (m)
      REAL(r_std),DIMENSION (nsnow)                           :: zcmprsfact   !! Compression factor due to densification from melting
      REAL(r_std),DIMENSION (nsnow)                           :: zscap        !! Snow heat capacity (J/m3 K)
      REAL(r_std),DIMENSION (nsnow)                           :: zsnowliq     !! (m)
      REAL(r_std),DIMENSION (nsnow)                           :: snowtemp_old
      REAL(r_std),DIMENSION (0:nsnow)                         :: zflowliqt    !!(m)
      REAL(r_std)                                             :: zrainfall,zpcpxs
      REAL(r_std)                                             :: ztotwcap
      REAL(r_std),DIMENSION(kjpindex,nsnow)                   :: snowdz_old,snowliq_old
      INTEGER(i_std)                                          :: jj,ji, iv
      REAL(r_std),DIMENSION(nsnow)                            :: snowdz_old2
      REAL(r_std),DIMENSION(nsnow)                            :: zsnowrho 
    
      !initialize
      snowdz_old = snowdz
      snowliq_old = snowliq

      DO ji = 1, kjpindex


         snowmass(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
         IF ((snowmass(ji) .GT. min_sechiba)) THEN

           !! 1 snow melting due to positive snowpack snow temperature
 
             !! 1.0 total liquid equivalent water content of each snow layer

               zsnowlwe(:) = snowrho(ji,:) * snowdz(ji,:)/ ph2o 

             !! 1.1 phase change (J/m2)

               pcapa_snow(ji,:) = snowrho(ji,:)*xci

               zphase(:)  = MIN(pcapa_snow(ji,:)*MAX(0.0, snowtemp(ji,:)-tp_00)*      &
                                snowdz(ji,:),                                       &
                            MAX(0.0,zsnowlwe(:)-snowliq(ji,:))*chalfu0*ph2o)

             !! 1.2 update snow liq water and temperature if melting

               zsnowmelt(:) = zphase(:)/(chalfu0*ph2o)

             !! 1.3 cool off the snow layer temperature due to melt

               zsnowtemp(:) = snowtemp(ji,:) - zphase(:)/(pcapa_snow(ji,:)* snowdz(ji,:))

               snowtemp(ji,:) = MIN(tp_00, zsnowtemp(:))

               zmeltxs(:)   = (zsnowtemp(:)-snowtemp(ji,:))*pcapa_snow(ji,:)*snowdz(ji,:)

             !! 1.4 loss of snowpack depth and liquid equivalent water

               zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension

               zcmprsfact(:) = (zsnowlwe(:)-MIN(snowliq(ji,:)+zsnowmelt(:),zwholdmax(:)))/ &
                               (zsnowlwe(:)-MIN(snowliq(ji,:),zwholdmax(:)))

               snowdz(ji,:)    = snowdz(ji,:)*zcmprsfact(:)

               snowrho(ji,:)     = zsnowlwe(:)*ph2o/snowdz(ji,:)

               snowliq(ji,:)   = snowliq(ji,:) + zsnowmelt(:)


           !! 2 snow refreezing process
              !! 2.1 freeze liquid water in any layer
               zscap(:)     = snowrho(ji,:)*xci  !J K-1 m-3
               zphase2(:)    = MIN(zscap(:)*                                          &
                                MAX(0.0, tp_00 - snowtemp(ji,:))*snowdz(ji,:),             &
                                snowliq(ji,:)*chalfu0*ph2o)

               ! warm layer and reduce liquid if freezing occurs
               zsnowdz(:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
               snowtemp_old(:) = snowtemp(ji,:) 
               snowtemp(ji,:) = snowtemp(ji,:) + zphase2(:)/(zscap(:)*zsnowdz(:))
               
               ! Reduce liquid portion if freezing occurs:
               snowliq(ji,:) = snowliq(ji,:) - ( (snowtemp(ji,:)-snowtemp_old(:))*       &
               zscap(:)*zsnowdz(:)/(chalfu0*ph2o) )
               snowliq(ji,:) = MAX(snowliq(ji,:), 0.0)
           !! 3. thickness change due to snowmelt in excess of holding capacity
               zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
               flowliq(:) = MAX(0.,(snowliq(ji,:)-zwholdmax(:)))
               snowliq(ji,:)  = snowliq(ji,:) - flowliq(:)
               snowdz(ji,:) = snowdz(ji,:) - flowliq(:)*ph2o/snowrho(ji,:)
               ! to prevent possible very small negative values (machine
               ! prescision as snow vanishes
               snowdz(ji,:) = MAX(0.0, snowdz(ji,:)) 

           !! 4. liquid water flow
               ztotwcap = SUM(zwholdmax(:)) 
               ! Rain entering snow (m):
               !ORCHIDEE has assumed that all precipitation entering snow has
               !left the snowpack and finally become runoff in hydrolc_soil. Here we just put zrainfall as 0.0
               !!zrainfall = precip_rain(ji)/ph2o ! rainfall (m)
               zrainfall = 0.0

               zflowliqt(0) = MIN(zrainfall,ztotwcap)
               ! Rain assumed to directly pass through the pack to runoff (m):
               zpcpxs = zrainfall - zflowliqt(0)
               !
               DO jj=1,nsnow
                  zflowliqt(jj) = flowliq(jj)
               ENDDO

               ! translated into a density increase:
               flowliq(:) = 0.0                ! clear this array for work
               zsnowliq(:) = snowliq(ji,:)      ! reset liquid water content
               !

               DO jj=1,nsnow
                  snowliq(ji,jj)  = snowliq(ji,jj) + zflowliqt(jj-1)
                  flowliq(jj)   = MAX(0.0, (snowliq(ji,jj)-zwholdmax(jj)))
                  snowliq(ji,jj)  = snowliq(ji,jj) - flowliq(jj)

                  !Modified by Tao Wang based on previous ISBA-ES scheme
                  snowrho(ji,jj)  = snowrho(ji,jj)  + (snowliq(ji,jj) - zsnowliq(jj))*       &
                                        & ph2o/MAX(xsnowdmin/nsnow,snowdz(ji,jj))

                  zflowliqt(jj) = zflowliqt(jj) + flowliq(jj)  
               ENDDO

               snowmelt(ji)  = snowmelt(ji) + (zflowliqt(nsnow) + zpcpxs) *  ph2o

             ! excess heat from melting, using it to warm underlying ground to conserve energy
               meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba  ! (W/m2)

             ! energy flux into the soil 
               grndflux(ji) = grndflux(ji) + meltxs(ji)  

         ELSE
	       snowdz(ji,:)=0.
               snowliq(ji,:)=0.
	       snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
         ENDIF

      ENDDO

    END SUBROUTINE explicitsnow_melt_refrz

!================================================================================================================================
!! SUBROUTINE   : explicitsnow_levels
!!
!>\BRIEF        Computes snow discretization based on given total snow depth
!!                
!! DESCRIPTION  : 
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : 
!!
!! FLOWCHART    : None
!! \n
!_
!================================================================================================================================ 
  SUBROUTINE explicitsnow_levels( kjpindex,snow_thick, snowdz)
    
    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                                     :: kjpindex     !! Domain size
    REAL(r_std),DIMENSION (kjpindex),   INTENT (in)                :: snow_thick   !! Total snow depth

    !! 0.2 Output variables

    !! 0.3 Modified variables

    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout)         :: snowdz                           !! Snow depth

    !! 0.4 Local variables

    INTEGER(i_std)                                                 :: il,it,ji
    INTEGER(i_std)                                                 :: i,j, ix
    REAL(r_std), DIMENSION(kjpindex)                               :: xi, xf
    REAL(r_std), PARAMETER, DIMENSION(3)                           :: ZSGCOEF1  = (/0.25, 0.50, 0.25/) !! Snow grid parameters 
    REAL(r_std), PARAMETER, DIMENSION(2)                           :: ZSGCOEF2  = (/0.05, 0.34/)       !! Snow grid parameters 
    REAL(r_std), PARAMETER                                         :: ZSNOWTRANS = 0.20                !! Minimum total snow depth at which surface layer thickness is constant (m)
    REAL(r_std), PARAMETER                                         :: XSNOWCRITD = 0.03                !! (m)

    DO ji=1,kjpindex
       IF ( snow_thick(ji) .LE. (XSNOWCRITD+0.01)) THEN

        snowdz(ji,1) = MIN(0.01, snow_thick(ji)/nsnow)
        snowdz(ji,3) = MIN(0.01, snow_thick(ji)/nsnow)
        snowdz(ji,2) = snow_thick(ji) - snowdz(ji,1) - snowdz(ji,3)
   
       ENDIF 
    ENDDO

    WHERE ( snow_thick(:) .LE. ZSNOWTRANS .AND. &
            snow_thick(:) .GT. (XSNOWCRITD+0.01) )
       !
        snowdz(:,1) = snow_thick(:)*ZSGCOEF1(1)
        snowdz(:,2) = snow_thick(:)*ZSGCOEF1(2) 
        snowdz(:,3) = snow_thick(:)*ZSGCOEF1(3)
       !
    END WHERE

     DO ji = 1,kjpindex
      IF (snow_thick(ji) .GT. ZSNOWTRANS) THEN
          snowdz(ji,1) = ZSGCOEF2(1)
          snowdz(ji,2) = (snow_thick(ji)-ZSGCOEF2(1))*ZSGCOEF2(2) + ZSGCOEF2(1)
        !When using simple finite differences, limit the thickness
        !factor between the top and 2nd layers to at most 10
          snowdz(ji,2)  = MIN(10*ZSGCOEF2(1),  snowdz(ji,2) )
          snowdz(ji,3)  = snow_thick(ji) - snowdz(ji,2) - snowdz(ji,1)

      ENDIF
    ENDDO
      

  END SUBROUTINE explicitsnow_levels

!!
!================================================================================================================================
!! SUBROUTINE   : explicitsnow_profile
!!
!>\BRIEF        
!!
!! DESCRIPTION  : In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile.
!!
!! RECENT CHANGE(S) : None
!! 
!! MAIN OUTPUT VARIABLE(S): 
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None 
!! \n 
!_
!================================================================================================================================
  SUBROUTINE explicitsnow_profile (kjpindex, cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)

  !! 0. Variables and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                               :: kjpindex     !! Domain size (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_sol_new !! skin temperature
    REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in)      :: cgrnd_snow   !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in)      :: dgrnd_snow   !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex),INTENT(in)             :: lambda_snow  !! Coefficient of the linear extrapolation of surface temperature 
    REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in)       :: snowdz       !! Snow layer thickness
    REAL(r_std), DIMENSION (kjpindex),INTENT(inout)          :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)

    !! 0.3 Modified variable

    !! 0.2 Output variables

    !! 0.3 Modified variables
    REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout)   :: snowtemp

    !! 0.4 Local variables

    INTEGER(i_std)                                           :: ji, jg
!_
!================================================================================================================================
  !! 1. Computes the snow temperatures ptn.
    DO ji = 1,kjpindex
      IF (SUM(snowdz(ji,:)) .GT. 0) THEN
        snowtemp(ji,1) = (lambda_snow(ji) * cgrnd_snow(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / &
             (lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un)
        temp_sol_add(ji) = zero 
        DO jg = 1,nsnow-1
            snowtemp(ji,jg+1) = cgrnd_snow(ji,jg) + dgrnd_snow(ji,jg) * snowtemp(ji,jg)
        ENDDO
      ENDIF
    ENDDO
    IF (printlev>=3) WRITE (numout,*) ' explicitsnow_profile done '

  END SUBROUTINE explicitsnow_profile

END MODULE explicitsnow
