mode_snow3l.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ##################
      MODULE mode_snow3l
!     ##################
!
!!****  *MODE_SNOW * - contains explicit snow (ISBA-ES) characteristics functions
!!                     for total liquid water holding capacity of a snow layer (m)
!!                     and the thermal heat capacity of a layer (J K-1 m-3)
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!    direct calculation
!!
!!    EXTERNAL
!!    --------
!!       
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!    REFERENCE
!!    ---------
!!    Boone and Etchevers, J. HydroMeteor., 2001
!!      
!!
!!    AUTHOR
!!    ------
!!      A. Boone       * Meteo France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original        01/08/02
!!      V. Masson       01/2004  add snow grid computations
!!      V. Vionnet      06/2008 -Introduction of Crocus formulation to
!                       calculate maximum liquid water holding capacity
!!                              - New algorithm to compute snow grid :
!                       10 layers
!!                              - Routine to aggregate snow grain type
!                       from 2 layers    
!!                              _ Routine to compute average grain
!                       type when snow depth< 3 cm. 
!     S. Morin          02/2011 - Add routines for Crocus
!     A. Boone          02/2012 - Add optimization of do-loops.
!     C. Carmagnola     12/2012 - Add the case CSNOWMETAMO!='B92' in subroutine SNOW3LAVGRAIN and in function SNOW3LDIFTYP
!     M. Lafaysse       01/2013 - Remove SNOWCROWLIQMAX routines (not used)
!     M. Lafaysse       08/2013 - simplification of routine SNOW3LAVGRAIN (logical GDENDRITIC)
!     B. Decharme       07/2013 - SNOW3LGRID cleanning 
!                                 New algorithm to compute snow grid for 6-L or 12-L
!     A. Boone          10/2014 - Added snow thermal conductivity routines
!     B. Decharme       01/2015 - Added optical snow grain size diameter
!----------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!
!
INTERFACE snow3lwliqmax
  MODULE PROCEDURE snow3lwliqmax_3d
  MODULE PROCEDURE snow3lwliqmax_2d
  MODULE PROCEDURE snow3lwliqmax_1d
END INTERFACE
!
INTERFACE snow3lhold
  MODULE PROCEDURE snow3lhold_3d
  MODULE PROCEDURE snow3lhold_2d
  MODULE PROCEDURE snow3lhold_1d
  MODULE PROCEDURE snow3lhold_0d
END INTERFACE
INTERFACE snowcrohold
  MODULE PROCEDURE snowcrohold_3d
  MODULE PROCEDURE snowcrohold_2d
  MODULE PROCEDURE snowcrohold_1d
  MODULE PROCEDURE snowcrohold_0d
END INTERFACE
!
INTERFACE snow3lscap
  MODULE PROCEDURE snow3lscap_3d
  MODULE PROCEDURE snow3lscap_2d
  MODULE PROCEDURE snow3lscap_1d
  MODULE PROCEDURE snow3lscap_0d
END INTERFACE
!
INTERFACE snow3l_marbouty
  MODULE PROCEDURE snow3l_marbouty
END INTERFACE
!
INTERFACE snow3lgrid
  MODULE PROCEDURE snow3lgrid_2d
  MODULE PROCEDURE snow3lgrid_1d
END INTERFACE
!
INTERFACE snow3lagreg
  MODULE PROCEDURE snow3lagreg
END INTERFACE
!
INTERFACE snow3lavgrain
  MODULE PROCEDURE snow3lavgrain
END INTERFACE
!
INTERFACE snow3ldiftyp
  MODULE PROCEDURE snow3ldiftyp
END INTERFACE
!
INTERFACE get_mass_heat
  MODULE PROCEDURE get_mass_heat
END INTERFACE
!
INTERFACE get_diam
  MODULE PROCEDURE get_diam
END INTERFACE
!
INTERFACE snow3lradabs
  MODULE PROCEDURE snow3lradabs_2d
  MODULE PROCEDURE snow3lradabs_1d
  MODULE PROCEDURE snow3lradabs_0d
END INTERFACE
!
INTERFACE snow3lthrm
  MODULE PROCEDURE snow3lthrm
END INTERFACE
!
INTERFACE snow3ldopt
  MODULE PROCEDURE snow3ldopt_2d
  MODULE PROCEDURE snow3ldopt_1d
  MODULE PROCEDURE snow3ldopt_0d
END INTERFACE
!
INTERFACE snow3lalb
  MODULE PROCEDURE snow3lalb
END INTERFACE
!
!-------------------------------------------------------------------------------
 CONTAINS
!
!####################################################################
      FUNCTION snow3lwliqmax_3d(PSNOWRHO) RESULT(PWLIQMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_snow_par, ONLY : xrhosmax_es, xsnowrhohold,      &
                          xwsnowholdmax2, xwsnowholdmax1
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:,:), INTENT(IN)                                  :: psnowrho ! (kg/m3)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),SIZE(PSNOWRHO,3)) :: pwliqmax ! (kg/m3)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),SIZE(PSNOWRHO,3)) :: zholdmaxr, zsnowrho
REAL(KIND=JPRB) :: zhook_handle
!-----------------------------------------------------------------------
! Evaluate capacity using upper density limit:
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LWLIQMAX_3D',0,zhook_handle)
zsnowrho(:,:,:) = min(xrhosmax_es, psnowrho(:,:,:))
!
! Maximum ratio of liquid to SWE:
!
zholdmaxr(:,:,:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)*    &
                  max(0.,xsnowrhohold-zsnowrho(:,:,:))/xsnowrhohold 
!
! Maximum liquid water holding capacity of the snow (kg/m3):
!
pwliqmax(:,:,:) = zholdmaxr(:,:,:)*zsnowrho(:,:,:)           
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LWLIQMAX_3D',1,zhook_handle)
!
END FUNCTION snow3lwliqmax_3d
!####################################################################
      FUNCTION snow3lwliqmax_2d(PSNOWRHO) RESULT(PWLIQMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_snow_par, ONLY : xrhosmax_es, xsnowrhohold,      &
                          xwsnowholdmax2, xwsnowholdmax1
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)                   :: psnowrho ! (kg/m3)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: pwliqmax ! (kg/m3)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zholdmaxr, zsnowrho
REAL(KIND=JPRB) :: zhook_handle
!-----------------------------------------------------------------------
! Evaluate capacity using upper density limit:
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LWLIQMAX_2D',0,zhook_handle)
zsnowrho(:,:) = min(xrhosmax_es, psnowrho(:,:))
!
! Maximum ratio of liquid to SWE:
!
zholdmaxr(:,:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)*    &
                  max(0.,xsnowrhohold-zsnowrho(:,:))/xsnowrhohold 
!
! Maximum liquid water holding capacity of the snow (kg/m3):
!
pwliqmax(:,:) = zholdmaxr(:,:)*zsnowrho(:,:)               
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LWLIQMAX_2D',1,zhook_handle)
!
END FUNCTION snow3lwliqmax_2d
!####################################################################
      FUNCTION snow3lwliqmax_1d(PSNOWRHO) RESULT(PWLIQMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_snow_par, ONLY : xrhosmax_es, xsnowrhohold,      &
                           xwsnowholdmax2, xwsnowholdmax1
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)  :: psnowrho ! (kg/m3)
!
REAL, DIMENSION(SIZE(PSNOWRHO)) :: pwliqmax ! (kg/m3)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO)) :: zholdmaxr, zsnowrho
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------
! Evaluate capacity using upper density limit:
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LWLIQMAX_1D',0,zhook_handle)
zsnowrho(:) = min(xrhosmax_es, psnowrho(:))
!
! Maximum ratio of liquid to SWE:
!
zholdmaxr(:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)*    &
                  max(0.,xsnowrhohold-zsnowrho(:))/xsnowrhohold 
!
! Maximum liquid water holding capacity of the snow (kg/m3):
!
pwliqmax(:) = zholdmaxr(:)*zsnowrho(:)                
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LWLIQMAX_1D',1,zhook_handle)
!
END FUNCTION snow3lwliqmax_1d

!####################################################################
!####################################################################
!####################################################################
      FUNCTION snow3lhold_3d(PSNOWRHO,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw
USE modd_snow_par, ONLY : xrhosmax_es, xsnowrhohold,      &
                          xwsnowholdmax2, xwsnowholdmax1
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:,:), INTENT(IN)                 :: psnowdz, psnowrho
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),SIZE(PSNOWRHO,3)) :: pwholdmax
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),SIZE(PSNOWRHO,3)) :: zholdmaxr, zsnowrho
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------
! Evaluate capacity using upper density limit:
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_3D',0,zhook_handle)
zsnowrho(:,:,:) = min(xrhosmax_es, psnowrho(:,:,:))
!
! Maximum ratio of liquid to SWE:
!
zholdmaxr(:,:,:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)*    &
                  max(0.,xsnowrhohold-zsnowrho(:,:,:))/xsnowrhohold 
!
! Maximum liquid water holding capacity of the snow (m):
!
pwholdmax(:,:,:) = zholdmaxr(:,:,:)*psnowdz(:,:,:)*zsnowrho(:,:,:)/xrholw
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_3D',1,zhook_handle)
!
END FUNCTION snow3lhold_3d
!####################################################################
      FUNCTION snow3lhold_2d(PSNOWRHO,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw
USE modd_snow_par, ONLY : xrhosmax_es, xsnowrhohold,      &
                          xwsnowholdmax2, xwsnowholdmax1
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)                   :: psnowdz, psnowrho
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: pwholdmax
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zholdmaxr, zsnowrho
REAL(KIND=JPRB) :: zhook_handle
!-----------------------------------------------------------------------
! Evaluate capacity using upper density limit:
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_2D',0,zhook_handle)
zsnowrho(:,:) = min(xrhosmax_es, psnowrho(:,:))
!
! Maximum ratio of liquid to SWE:
!
zholdmaxr(:,:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)*    &
                  max(0.,xsnowrhohold-zsnowrho(:,:))/xsnowrhohold 
!
! Maximum liquid water holding capacity of the snow (m):
!
pwholdmax(:,:) = zholdmaxr(:,:)*psnowdz(:,:)*zsnowrho(:,:)/xrholw
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_2D',1,zhook_handle)
!
END FUNCTION snow3lhold_2d
!####################################################################
      FUNCTION snow3lhold_1d(PSNOWRHO,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw
USE modd_snow_par, ONLY : xrhosmax_es, xsnowrhohold,     &
                          xwsnowholdmax2, xwsnowholdmax1 
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)                     :: psnowdz, psnowrho
!
REAL, DIMENSION(SIZE(PSNOWRHO))                    :: pwholdmax
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO))                    :: zholdmaxr, zsnowrho
REAL(KIND=JPRB) :: zhook_handle
!-----------------------------------------------------------------------
! Evaluate capacity using upper density limit:
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_1D',0,zhook_handle)
zsnowrho(:) = min(xrhosmax_es, psnowrho(:))
!
! Maximum ratio of liquid to SWE:
!
zholdmaxr(:) = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1)*     &
                  max(0.,xsnowrhohold-zsnowrho(:))/xsnowrhohold 
!
! Maximum liquid water holding capacity of the snow (m):
!
pwholdmax(:) = zholdmaxr(:)*psnowdz(:)*zsnowrho(:)/xrholw               
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_1D',1,zhook_handle)
!
END FUNCTION snow3lhold_1d
!####################################################################
      FUNCTION snow3lhold_0d(PSNOWRHO,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw
USE modd_snow_par, ONLY : xrhosmax_es, xsnowrhohold,     &
                          xwsnowholdmax2, xwsnowholdmax1
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)        :: psnowdz, psnowrho
!
REAL                    :: pwholdmax
!
!*      0.2    declarations of local variables
!
REAL                    :: zholdmaxr, zsnowrho
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
! Evaluate capacity using upper density limit:
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_0D',0,zhook_handle)
zsnowrho = min(xrhosmax_es, psnowrho)
!
! Maximum ratio of liquid to SWE:
!
zholdmaxr = xwsnowholdmax1 + (xwsnowholdmax2-xwsnowholdmax1) *     &
                             max(0.,xsnowrhohold-zsnowrho)/xsnowrhohold 
!
! Maximum liquid water holding capacity of the snow (m):
!
pwholdmax = zholdmaxr*psnowdz*zsnowrho/xrholw              
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LHOLD_0D',1,zhook_handle)
!
END FUNCTION snow3lhold_0d
!####################################################################
      FUNCTION snowcrohold_3d(PSNOWRHO,PSNOWLIQ,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw,xrholi
USE modd_snow_par, ONLY : xpercentagepore
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:,:), INTENT(IN)                 :: psnowdz, psnowliq, psnowrho
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),SIZE(PSNOWRHO,3)) :: pwholdmax
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
! computation of water holding capacity based on Crocus, 
!taking into account the conversion between wet and dry density - 
!S. Morin/V. Vionnet 2010 12 09

! PWHOLDMAX is expressed in m water for each layer
! In short, PWHOLDMAX = XPERCENTAGEPORE * porosity * PSNOWDZ .
! The porosity is computed as (rho_ice - (rho_snow - lwc))/(rho_ice)
! where everything has to be in kg m-3 units. In practice, since
! PSNOWLIQ is expressed in m water, expressing the lwc in kg m-3
! is achieved as PSNOWLIQ*XRHOLW/PSNOWDZ. After some rearranging one
! obtains the equation given above.
! Note that equation (19) in Vionnet et al., GMD 2012, is wrong,
! because it does not take into account the fact that liquid water
! content has to be substracted from total density to compute the
! porosity.


IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_3D',0,zhook_handle)
pwholdmax(:,:,:) = xpercentagepore/xrholi * (psnowdz * (xrholi-psnowrho) + psnowliq*xrholw)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_3D',1,zhook_handle)
!
END FUNCTION snowcrohold_3d
!####################################################################
      FUNCTION snowcrohold_2d(PSNOWRHO,PSNOWLIQ,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw,xrholi
USE modd_snow_par, ONLY : xpercentagepore
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)                   :: psnowdz, psnowrho, psnowliq
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: pwholdmax
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
! computation of water holding capacity based on Crocus, 
!taking into account the conversion between wet and dry density - 
!S. Morin/V. Vionnet 2010 12 09

! PWHOLDMAX is expressed in m water for each layer
! In short, PWHOLDMAX = XPERCENTAGEPORE * porosity * PSNOWDZ .
! The porosity is computed as (rho_ice - (rho_snow - lwc))/(rho_ice)
! where everything has to be in kg m-3 units. In practice, since
! PSNOWLIQ is expressed in m water, expressing the lwc in kg m-3
! is achieved as PSNOWLIQ*XRHOLW/PSNOWDZ. After some rearranging one
! obtains the equation given above.
! Note that equation (19) in Vionnet et al., GMD 2012, is wrong,
! because it does not take into account the fact that liquid water
! content has to be substracted from total density to compute the
! porosity.

IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_2D',0,zhook_handle)
pwholdmax(:,:) = xpercentagepore/xrholi * (psnowdz * (xrholi-psnowrho) + psnowliq*xrholw)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_2D',1,zhook_handle)
!
END FUNCTION snowcrohold_2d
!####################################################################
!####################################################################
!####################################################################
      FUNCTION snowcrohold_1d(PSNOWRHO,PSNOWLIQ,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw,xrholi
USE modd_snow_par, ONLY : xpercentagepore
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)                     :: psnowdz, psnowrho, psnowliq
!
REAL, DIMENSION(SIZE(PSNOWRHO))                    :: pwholdmax
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
! computation of water holding capacity based on Crocus, 
!taking into account the conversion between wet and dry density -
!S. Morin/V. Vionnet 2010 12 09

! PWHOLDMAX is expressed in m water for each layer
! In short, PWHOLDMAX = XPERCENTAGEPORE * porosity * PSNOWDZ .
! The porosity is computed as (rho_ice - (rho_snow - lwc))/(rho_ice)
! where everything has to be in kg m-3 units. In practice, since
! PSNOWLIQ is expressed in m water, expressing the lwc in kg m-3
! is achieved as PSNOWLIQ*XRHOLW/PSNOWDZ. After some rearranging one
! obtains the equation given above.
! Note that equation (19) in Vionnet et al., GMD 2012, is wrong,
! because it does not take into account the fact that liquid water
! content has to be substracted from total density to compute the
! porosity.

IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_1D',0,zhook_handle)
pwholdmax(:) = xpercentagepore/xrholi * (psnowdz * (xrholi-psnowrho) + psnowliq*xrholw)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_1D',1,zhook_handle)
!
END FUNCTION snowcrohold_1d
!####################################################################
      FUNCTION snowcrohold_0d(PSNOWRHO,PSNOWLIQ,PSNOWDZ) RESULT(PWHOLDMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the maximum liquid water holding capacity of
!     snow layer(s).
!
USE modd_csts,     ONLY : xrholw,xrholi
USE modd_snow_par, ONLY : xpercentagepore
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)        :: psnowdz, psnowrho, psnowliq
!
REAL                    :: pwholdmax
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
! computation of water holding capacity based on Crocus, 
!taking into account the conversion between wet and dry density - 
!S. Morin/V. Vionnet 2010 12 09

! PWHOLDMAX is expressed in m water for each layer
! In short, PWHOLDMAX = XPERCENTAGEPORE * porosity * PSNOWDZ .
! The porosity is computed as (rho_ice - (rho_snow - lwc))/(rho_ice)
! where everything has to be in kg m-3 units. In practice, since
! PSNOWLIQ is expressed in m water, expressing the lwc in kg m-3
! is achieved as PSNOWLIQ*XRHOLW/PSNOWDZ. After some rearranging one
! obtains the equation given above.
! Note that equation (19) in Vionnet et al., GMD 2012, is wrong,
! because it does not take into account the fact that liquid water
! content has to be substracted from total density to compute the
! porosity.

IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_0D',0,zhook_handle)
pwholdmax = xpercentagepore/xrholi * (psnowdz * (xrholi-psnowrho) + psnowliq*xrholw)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOWCROHOLD_0D',1,zhook_handle)
!
END FUNCTION snowcrohold_0d
!####################################################################
!####################################################################
!####################################################################
      FUNCTION snow3lscap_3d(PSNOWRHO) RESULT(PSCAP)
!
!!    PURPOSE
!!    -------
!     Calculate the heat capacity of a snow layer.
!
USE modd_csts,ONLY : xci
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:,:), INTENT(IN)                                  :: psnowrho
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),SIZE(PSNOWRHO,3)) :: pscap
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!     The method of Verseghy (1991), Int. J. Climat., 11, 111-133.
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_3D',0,zhook_handle)
pscap(:,:,:) = psnowrho(:,:,:)*xci      ! (J K-1 m-3)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_3D',1,zhook_handle)
!
END FUNCTION snow3lscap_3d
!####################################################################
      FUNCTION snow3lscap_2d(PSNOWRHO) RESULT(PSCAP)
!
!!    PURPOSE
!!    -------
!     Calculate the heat capacity of a snow layer.
!
USE modd_csts,ONLY : xci
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)                   :: psnowrho
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: pscap
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!     The method of Verseghy (1991), Int. J. Climat., 11, 111-133.
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_2D',0,zhook_handle)
pscap(:,:) = psnowrho(:,:)*xci      ! (J K-1 m-3)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_2D',1,zhook_handle)
!
END FUNCTION snow3lscap_2d
!####################################################################
      FUNCTION snow3lscap_1d(PSNOWRHO) RESULT(PSCAP)
!
!!    PURPOSE
!!    -------
!     Calculate the heat capacity of a snow layer.
!
USE modd_csts,ONLY : xci
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)                   :: psnowrho
!
REAL, DIMENSION(SIZE(PSNOWRHO))                  :: pscap
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!     The method of Verseghy (1991), Int. J. Climat., 11, 111-133.
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_1D',0,zhook_handle)
pscap(:) = psnowrho(:)*xci      ! (J K-1 m-3)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_1D',1,zhook_handle)
!
END FUNCTION snow3lscap_1d
!####################################################################
      FUNCTION snow3lscap_0d(PSNOWRHO) RESULT(PSCAP)
!
!!    PURPOSE
!!    -------
!     Calculate the heat capacity of a snow layer.
!
USE modd_csts,ONLY : xci
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)       :: psnowrho
!
REAL                   :: pscap
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!     The method of Verseghy (1991), Int. J. Climat., 11, 111-133.
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_0D',0,zhook_handle)
pscap = psnowrho*xci      ! (J K-1 m-3)
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LSCAP_0D',1,zhook_handle)
!
END FUNCTION snow3lscap_0d
!
!####################################################################
!####################################################################
!####################################################################
      FUNCTION snow3l_marbouty(PSNOWRHO,PSNOWTEMP,PGRADT) RESULT(PDANGL) 
!**** *ZDANGL* - CROISSANCE DES GRAINS NON DENDRITIQUES ET ANGULEUX .
!              - GROWTH RATES FOR NON DENDRITIC GRAINS WITH SPHERICITY=0 


!     OBJET.
!     ------

!**   INTERFACE.
!     ----------

!     *ZDANGL*(PST,PSRO,PGRADT)*

!        *PST*     -  TEMPERATURE DE LA STRATE DE NEIGE.
!        *PSRO*    -  MASSE VOLUMIQUE DE LA STRATE DE NEIGE.
!        *PGRADT*  -  GRADIENT DE TEMPERATURE AFFECTANT LA STRATE DE NEIGE.

!     METHODE.
!     --------
!     THE FUNCTION RETURN A VALUE BETWEEN 0 AND 1 WHICH IS USED IN THE DETERMINATION OF THE 
!     GROWTH RATE FOR THE CONSIDERED LAYER.
!     THIS VALUE (WITHOUT UNIT) IS THE PRODUCT OF 3 FUNCTIONS (WHICH HAVE THEIR SOLUTIONS BETWEEN 0 AND 1) :
!     F(TEMPERATURE) * H(DENSITY) * G(TEMPERATURE GRADIENT)

!     EXTERNES.
!     ---------

!     REFERENCES.
!     -----------
!        MARBOUTY D (1980) AN EXPERIMENTAL STUDY OF TEMPERATURE GRADIENT 
!                          METAMORPHISM J GLACIOLOGY

!     AUTEURS.
!     --------
!        DOMINIQUE MARBOUTY (FMARBO/GMARBO/HMARBO).

!     MODIFICATIONS.
!     --------------
!        08/95: YANNICK DANIELOU - CODAGE A LA NORME DOCTOR.
!        03/06: JM WILLEMET      - F90 AND SI UNITS
!        01/08: JM WILLEMET      - ERROR ON THE FORMULATION OF G FUNCTION (WITH GRADIENT) IS CORRECTED 

USE modd_csts, ONLY : xtt
USE modd_snow_metamo  
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!     DECLARATIONS.
!     -------------
!
REAL ,INTENT(IN) :: psnowtemp, psnowrho, pgradt
!
REAL             :: pdangl
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3L_MARBOUTY',0,zhook_handle)
!
pdangl = 0.0
!
! INFLUENCE DE LA TEMPERATURE /TEMPERATURE INFLUENCE.
IF( psnowtemp>=xtt-xvtang1 ) THEN
  !
  IF ( psnowtemp>=xtt-xvtang2 ) THEN
    pdangl = xvtang4 + xvtang5 * (xtt-psnowtemp) / xvtang6
  ELSEIF( psnowtemp>=xtt-xvtang3 ) THEN
    pdangl = xvtang7 - xvtang8 * (xtt-xvtang2-psnowtemp) / xvtang9
  ELSE
    pdangl = xvtanga - xvtangb * (xtt-xvtang3-psnowtemp) / xvtangc
  ENDIF
  !
  ! INFLUENCE DE LA MASSE VOLUMIQUE / DENSITY INFLUENCE.
  IF ( psnowrho<=xvrang1 ) THEN
    !
    IF ( psnowrho>xvrang2 ) THEN
      pdangl = pdangl * ( 1. - (psnowrho-xvrang2)/(xvrang1-xvrang2) )
    ENDIF
    !
    ! INFLUENCE DU GRADIENT DE TEMPERATURE / TEMPERATURE GRADIENT INFLUENCE.
    IF ( pgradt<=xvgang1 ) THEN
      !
      IF ( pgradt<=xvgang2 ) THEN
        pdangl = pdangl * xvgang5 * (pgradt-xvgang6)/(xvgang2-xvgang6)
      ELSEIF( pgradt<=xvgang3 ) THEN
        pdangl = pdangl * ( xvgang7 + xvgang8 * (pgradt-xvgang2)/(xvgang3-xvgang2) )
      ELSEIF( pgradt<=xvgang4 )THEN
        pdangl = pdangl * ( xvgang9 + xvganga * (pgradt-xvgang3)/(xvgang4-xvgang3) )
      ELSE
        pdangl = pdangl * ( xvgangb + xvgangc * (pgradt-xvgang4)/(xvgang1-xvgang4) )
      ENDIF
      !
    ENDIF
    !
  ELSE
    !
    pdangl = 0.
    !
  ENDIF
  !
ELSE
  !
  pdangl = 0.
  !
ENDIF
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3L_MARBOUTY',1,zhook_handle)
!
END FUNCTION snow3l_marbouty     
!       
!####################################################################
!####################################################################
!####################################################################
!
      SUBROUTINE snow3lgrid_2d(PSNOWDZ,PSNOW,PSNOWDZ_OLD)
!
!!    PURPOSE
!!    -------
!     Once during each time step, update grid to maintain
!     grid proportions. Similar to approach of Lynch-Steiglitz,
!     1994, J. Clim., 7, 1842-1855. Corresponding mass and
!     heat adjustments made directly after the call to this
!     routine. 3 grid configurations:
!     1) for very thin snow, constant grid spacing
!     2) for intermediate thicknesses, highest resolution at soil/snow
!        interface and at the snow/atmosphere interface
!     3) for deep snow, vertical resoution finest at snow/atmosphere
!        interface (set to a constant value) and increases with snow depth.
!        Second layer can't be more than an order of magnitude thicker
!        than surface layer.
!
!
USE modd_surf_par,   ONLY : xundef
USE modd_snow_par,   ONLY : xsnowcritd
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:  ), INTENT(IN )           :: psnow
REAL, DIMENSION(:,:), INTENT(OUT)           :: psnowdz
REAL, DIMENSION(:,:), INTENT(IN ), OPTIONAL :: psnowdz_old
!
!*      0.1    declarations of local variables
!
INTEGER                           :: jj, ji
!
INTEGER                           :: inlvls, ini
!   
REAL,     DIMENSION(SIZE(PSNOW))  :: zwork
!
LOGICAL , DIMENSION(SIZE(PSNOW))  :: gregrid

! ISBA-ES snow grid parameters
!
REAL, PARAMETER, DIMENSION(3)     :: zsgcoef1  = (/0.25, 0.50, 0.25/) 
REAL, PARAMETER, DIMENSION(2)     :: zsgcoef2  = (/0.05, 0.34/)       
!      
REAL, PARAMETER, DIMENSION(3)     :: zsgcoef   = (/0.3, 0.4, 0.3/) 
!
! Minimum total snow depth at which surface layer thickness is constant:
!
REAL, PARAMETER                   :: zsnowtrans = 0.20                ! (m)
!      
! Minimum snow depth by layer for 6-L or 12-L configuration :
!
REAL, PARAMETER                   ::  zdz1=0.01
REAL, PARAMETER                   ::  zdz2=0.05
REAL, PARAMETER                   ::  zdz3=0.15
REAL, PARAMETER                   ::  zdz4=0.50
REAL, PARAMETER                   ::  zdz5=1.00
REAL, PARAMETER                   ::  zdzn0=0.02
REAL, PARAMETER                   ::  zdzn1=0.1
REAL, PARAMETER                   ::  zdzn2=0.5
REAL, PARAMETER                   ::  zdzn3=1.0
!
REAL, PARAMETER                   ::  zcoef1 = 0.5
REAL, PARAMETER                   ::  zcoef2 = 1.5
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
! 0. Initialization:
! ------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LGRID_2D',0,zhook_handle)
!
inlvls = SIZE(psnowdz(:,:),2)
ini    = SIZE(psnowdz(:,:),1)
!
zwork(:) = 0.0
gregrid(:) = .true.
!
! 1. Calculate current grid for 3-layer (default) configuration):
! ---------------------------------------------------------------
! Based on formulation of Lynch-Stieglitz (1994)
! except for 3 modifications: 
! i) smooth transition here at ZSNOWTRANS
! ii) constant ratio for very thin snow:
! iii) ratio of layer 2 to surface layer <= 10
!
IF(inlvls == 1)THEN
!
  DO ji=1,ini
     psnowdz(ji,1)  = psnow(ji)
  ENDDO
!
ELSEIF(inlvls == 3)THEN
!
   WHERE(psnow <= xsnowcritd+0.01)
      psnowdz(:,1) = min(0.01, psnow(:)/inlvls)
      psnowdz(:,3) = min(0.01, psnow(:)/inlvls)
      psnowdz(:,2) = psnow(:) - psnowdz(:,1) - psnowdz(:,3)
   END WHERE
!
   WHERE(psnow <= zsnowtrans .AND. psnow > xsnowcritd+0.01)
      psnowdz(:,1) = psnow(:)*zsgcoef1(1)
      psnowdz(:,2) = psnow(:)*zsgcoef1(2)
      psnowdz(:,3) = psnow(:)*zsgcoef1(3)
   END WHERE
!
   WHERE(psnow > zsnowtrans)
      psnowdz(:,1) = zsgcoef2(1)
      psnowdz(:,2) = (psnow(:)-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
! 
      psnowdz(:,2) = min(10*zsgcoef2(1),  psnowdz(:,2))
      psnowdz(:,3) = psnow(:) - psnowdz(:,2) - psnowdz(:,1)
   END WHERE
!
!
! 2. Calculate current grid for 6-layer :
! ---------------------------------------------------------------
!
ELSEIF(inlvls == 6)THEN
!
! critere a satisfaire pour remaillage
!
  IF(present(psnowdz_old))THEN
    gregrid(:) = psnowdz_old(:,1) < zcoef1 * min(zdz1 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,1) > zcoef2 * min(zdz1 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,2) < zcoef1 * min(zdz2 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,2) > zcoef2 * min(zdz2 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,6) < zcoef1 * min(zdzn1,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,6) > zcoef2 * min(zdzn1,psnow(:)/inlvls)
  ENDIF
!
  WHERE(gregrid(:))  
!      top layers 
       psnowdz(:,1) = min(zdz1,psnow(:)/inlvls) 
       psnowdz(:,2) = min(zdz2,psnow(:)/inlvls) 
!      last layers 
       psnowdz(:,6) = min(zdzn1,psnow(:)/inlvls)
!      remaining snow for remaining layers
       zwork(:)     = psnow(:) - psnowdz(:,1) - psnowdz(:,2) - psnowdz(:,6)
       psnowdz(:,3) = zwork(:)*zsgcoef(1)
       psnowdz(:,4) = zwork(:)*zsgcoef(2)
       psnowdz(:,5) = zwork(:)*zsgcoef(3)
!      layer 3 tickness >= layer 2 tickness
       zwork(:)=min(0.0,psnowdz(:,3)-psnowdz(:,2))
       psnowdz(:,3)=psnowdz(:,3)-zwork(:)
       psnowdz(:,4)=psnowdz(:,4)+zwork(:) 
!      layer 5 tickness >= layer 6 tickness  
       zwork(:)=min(0.0,psnowdz(:,5)-psnowdz(:,6))
       psnowdz(:,5)=psnowdz(:,5)-zwork(:)
       psnowdz(:,4)=psnowdz(:,4)+zwork(:)
  ENDWHERE
!
! 3. Calculate current grid for 9-layer :
! ---------------------------------------------------------------
!
ELSEIF(inlvls == 9)THEN
!
! critere a satisfaire pour remaillage
!
  IF(present(psnowdz_old))THEN
    gregrid(:) = psnowdz_old(:,1) < zcoef1 * min(zdz1 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,1) > zcoef2 * min(zdz1 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,2) < zcoef1 * min(zdz2 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,2) > zcoef2 * min(zdz2 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,9) < zcoef1 * min(zdzn0,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,9) > zcoef2 * min(zdzn0,psnow(:)/inlvls) 
  ENDIF
!             
  WHERE(gregrid(:))  
!      top layers 
       psnowdz(:,1) = min(zdz1,psnow(:)/inlvls) 
       psnowdz(:,2) = min(zdz2,psnow(:)/inlvls) 
       psnowdz(:,3) = min(zdz3,psnow(:)/inlvls)
!      last layers 
       psnowdz(:,9)= min(zdzn0,psnow(:)/inlvls)
       psnowdz(:,8)= min(zdzn1,psnow(:)/inlvls)
       psnowdz(:,7)= min(zdzn2,psnow(:)/inlvls)
!      remaining snow for remaining layers
       zwork(:) = psnow(:) - psnowdz(:, 1) - psnowdz(:, 2) - psnowdz(:, 3) &
                           - psnowdz(:, 7) - psnowdz(:, 8) - psnowdz(:, 9)
       psnowdz(:,4) = zwork(:)*zsgcoef(1)
       psnowdz(:,5) = zwork(:)*zsgcoef(2)
       psnowdz(:,6) = zwork(:)*zsgcoef(3)
!      layer 4 tickness >= layer 3 tickness
       zwork(:)=min(0.0,psnowdz(:,4)-psnowdz(:,3))
       psnowdz(:,4)=psnowdz(:,4)-zwork(:)
       psnowdz(:,5)=psnowdz(:,5)+zwork(:) 
!      layer 6 tickness >= layer 7 tickness  
       zwork(:)=min(0.0,psnowdz(:,6)-psnowdz(:,7))
       psnowdz(:,6)=psnowdz(:,6)-zwork(:)
       psnowdz(:,5)=psnowdz(:,5)+zwork(:)
  ENDWHERE
!
! 4. Calculate current grid for 12-layer :
! ---------------------------------------------------------------
!
ELSEIF(inlvls == 12)THEN
!
! critere a satisfaire pour remaillage
!
  IF(present(psnowdz_old))THEN
    gregrid(:) = psnowdz_old(:, 1) < zcoef1 * min(zdz1 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:, 1) > zcoef2 * min(zdz1 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:, 2) < zcoef1 * min(zdz2 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:, 2) > zcoef2 * min(zdz2 ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,12) < zcoef1 * min(zdzn0,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,12) > zcoef2 * min(zdzn0,psnow(:)/inlvls) 
  ENDIF
!             
  WHERE(gregrid(:))  
!      top layers 
       psnowdz(:,1) = min(zdz1,psnow(:)/inlvls) 
       psnowdz(:,2) = min(zdz2,psnow(:)/inlvls) 
       psnowdz(:,3) = min(zdz3,psnow(:)/inlvls)
       psnowdz(:,4) = min(zdz4,psnow(:)/inlvls)
       psnowdz(:,5) = min(zdz5,psnow(:)/inlvls)
!      last layers 
       psnowdz(:,12)= min(zdzn0,psnow(:)/inlvls)
       psnowdz(:,11)= min(zdzn1,psnow(:)/inlvls)
       psnowdz(:,10)= min(zdzn2,psnow(:)/inlvls)
       psnowdz(:, 9)= min(zdzn3,psnow(:)/inlvls)
!      remaining snow for remaining layers
       zwork(:) = psnow(:) - psnowdz(:, 1) - psnowdz(:, 2) - psnowdz(:, 3) &
                           - psnowdz(:, 4) - psnowdz(:, 5) - psnowdz(:, 9) &
                           - psnowdz(:,10) - psnowdz(:,11) - psnowdz(:,12)
       psnowdz(:,6) = zwork(:)*zsgcoef(1)
       psnowdz(:,7) = zwork(:)*zsgcoef(2)
       psnowdz(:,8) = zwork(:)*zsgcoef(3)
!      layer 6 tickness >= layer 5 tickness
       zwork(:)=min(0.0,psnowdz(:,6)-psnowdz(:,5))
       psnowdz(:,6)=psnowdz(:,6)-zwork(:)
       psnowdz(:,7)=psnowdz(:,7)+zwork(:) 
!      layer 8 tickness >= layer 9 tickness  
       zwork(:)=min(0.0,psnowdz(:,8)-psnowdz(:,9))
       psnowdz(:,8)=psnowdz(:,8)-zwork(:)
       psnowdz(:,7)=psnowdz(:,7)+zwork(:)
  ENDWHERE
!
! 4. Calculate other non-optimized grid :
! ---------------------------------------
! 
ELSEIF(inlvls<10.AND.inlvls/=3.AND.inlvls/=6.AND.inlvls/=9) THEN
!
  DO jj=1,inlvls
     DO ji=1,ini
        psnowdz(ji,jj)  = psnow(ji)/inlvls
     ENDDO
  ENDDO
!
  psnowdz(:,inlvls) = psnowdz(:,inlvls) + (psnowdz(:,1) - min(0.05, psnowdz(:,1)))
  psnowdz(:,1)      = min(0.05, psnowdz(:,1))
!
ELSE !(inlvls>=10 and /=12)  
!
! critere a satisfaire pour remaillage
!
  IF(present(psnowdz_old))THEN
    gregrid(:) = psnowdz_old(:,     1) < zcoef1 * min(zdz1         ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,     1) > zcoef2 * min(zdz1         ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,     2) < zcoef1 * min(zdz2         ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,     2) > zcoef2 * min(zdz2         ,psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,inlvls) < zcoef1 * min(0.05*psnow(:),psnow(:)/inlvls) .OR. &
               & psnowdz_old(:,inlvls) > zcoef2 * min(0.05*psnow(:),psnow(:)/inlvls) 
  ENDIF
!
  WHERE(gregrid(:))  
       psnowdz(:,1     ) = min(zdz1         ,psnow(:)/inlvls) 
       psnowdz(:,2     ) = min(zdz2         ,psnow(:)/inlvls) 
       psnowdz(:,3     ) = min(zdz3         ,psnow(:)/inlvls)
       psnowdz(:,4     ) = min(zdz4         ,psnow(:)/inlvls)
       psnowdz(:,5     ) = min(zdz5         ,psnow(:)/inlvls)          
       psnowdz(:,inlvls) = min(0.05*psnow(:),psnow(:)/inlvls) 
  ENDWHERE
!
  DO jj=6,inlvls-1,1
     DO ji=1,ini
        IF(gregrid(ji))THEN           
          zwork(ji) = psnowdz(ji,1)+psnowdz(ji,2)+psnowdz(ji,3)+psnowdz(ji,4)+psnowdz(ji,5)
          psnowdz(ji,jj) = (psnow(ji)-zwork(ji)-psnowdz(ji,inlvls))/(inlvls-6) 
        ENDIF
     ENDDO
  ENDDO
!
ENDIF
!
DO jj=1,inlvls
   DO ji=1,ini
      IF(psnow(ji)==xundef)THEN
         psnowdz(ji,jj) = xundef
      ELSEIF(.NOT.gregrid(ji))THEN           
         psnowdz(ji,jj)=psnowdz_old(ji,jj)
      ENDIF      
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LGRID_2D',1,zhook_handle)
!
END SUBROUTINE snow3lgrid_2d
!####################################################################
!####################################################################
!####################################################################
!
      SUBROUTINE snow3lgrid_1d(PSNOWDZ,PSNOW,PSNOWDZ_OLD)
!
!!    PURPOSE
!!    -------
!     Once during each time step, update grid to maintain
!     grid proportions. Similar to approach of Lynch-Steiglitz,
!     1994, J. Clim., 7, 1842-1855. Corresponding mass and
!     heat adjustments made directly after the call to this
!     routine. 3 grid configurations:
!     1) for very thin snow, constant grid spacing
!     2) for intermediate thicknesses, highest resolution at soil/snow
!        interface and at the snow/atmosphere interface
!     3) for deep snow, vertical resoution finest at snow/atmosphere
!        interface (set to a constant value) and increases with snow depth.
!        Second layer can't be more than an order of magnitude thicker
!        than surface layer.
!
!
USE modd_surf_par,   ONLY : xundef
USE modd_snow_par,   ONLY : xsnowcritd
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL,               INTENT(IN )           :: psnow
REAL, DIMENSION(:), INTENT(OUT)           :: psnowdz
REAL, DIMENSION(:), INTENT(IN ), OPTIONAL :: psnowdz_old
!
!*      0.1    declarations of local variables
!
INTEGER jj
!
INTEGER                           :: inlvls
!
REAL                              :: zwork
!
! modif_EB pour maillage
LOGICAL                           :: gregrid

! ISBA-ES snow grid parameters
!
REAL, PARAMETER, DIMENSION(3)     :: zsgcoef1  = (/0.25, 0.50, 0.25/) 
REAL, PARAMETER, DIMENSION(2)     :: zsgcoef2  = (/0.05, 0.34/)       
!      
REAL, PARAMETER, DIMENSION(3)     :: zsgcoef   = (/0.3, 0.4, 0.3/) 
!
! Minimum total snow depth at which surface layer thickness is constant:
!
REAL, PARAMETER                   :: zsnowtrans  = 0.20                ! (m)
!      
! Minimum snow depth by layer for 6-L or 12-L configuration :
!
REAL, PARAMETER                   ::  zdz1=0.01
REAL, PARAMETER                   ::  zdz2=0.05
REAL, PARAMETER                   ::  zdz3=0.15
REAL, PARAMETER                   ::  zdz4=0.50
REAL, PARAMETER                   ::  zdz5=1.00
REAL, PARAMETER                   ::  zdzn0=0.02
REAL, PARAMETER                   ::  zdzn1=0.1
REAL, PARAMETER                   ::  zdzn2=0.5
REAL, PARAMETER                   ::  zdzn3=1.0
!
REAL, PARAMETER                   ::  zcoef1 = 0.5
REAL, PARAMETER                   ::  zcoef2 = 1.5
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
! 0. Initialization:
! ------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LGRID_1D',0,zhook_handle)
!
inlvls = SIZE(psnowdz(:),1)
!
gregrid = .true.
!
! 1. Calculate current grid for 3-layer (default) configuration):
! ---------------------------------------------------------------
! Based on formulation of Lynch-Stieglitz (1994)
! except for 3 modifications: 
! i) smooth transition here at ZSNOWTRANS
! ii) constant ratio for very thin snow:
! iii) ratio of layer 2 to surface layer <= 10
!
IF(inlvls == 1)THEN
!
  psnowdz(1) = psnow
!
ELSEIF(inlvls == 3)THEN
!
   IF(psnow <= xsnowcritd+0.01)THEN
      psnowdz(1) = min(0.01, psnow/inlvls)
      psnowdz(3) = min(0.01, psnow/inlvls)
      psnowdz(2) = psnow - psnowdz(1) - psnowdz(3)
   ENDIF
!
   IF(psnow <= zsnowtrans .AND. psnow > xsnowcritd+0.01)THEN
      psnowdz(1) = psnow*zsgcoef1(1)
      psnowdz(2) = psnow*zsgcoef1(2)
      psnowdz(3) = psnow*zsgcoef1(3)
   ENDIF
!
   IF(psnow > zsnowtrans)THEN
      psnowdz(1) = zsgcoef2(1)
      psnowdz(2) = (psnow-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
! 
      psnowdz(2) = min(10*zsgcoef2(1),  psnowdz(2))
      psnowdz(3) = psnow - psnowdz(2) - psnowdz(1)
   END IF
!
!
! 2. Calculate current grid for 6-layer :
! ---------------------------------------------------------------
!
ELSEIF(inlvls == 6)THEN
!
! critere a satisfaire pour remaillage
!
  IF(present(psnowdz_old))THEN
    gregrid    = psnowdz_old(1) < zcoef1 * min(zdz1 ,psnow/inlvls) .OR. &
               & psnowdz_old(1) > zcoef2 * min(zdz1 ,psnow/inlvls) .OR. &
               & psnowdz_old(2) < zcoef1 * min(zdz2 ,psnow/inlvls) .OR. &
               & psnowdz_old(2) > zcoef2 * min(zdz2 ,psnow/inlvls) .OR. &
               & psnowdz_old(6) < zcoef1 * min(zdzn1,psnow/inlvls) .OR. &
               & psnowdz_old(6) > zcoef2 * min(zdzn1,psnow/inlvls)
  ENDIF
!
  IF(gregrid)THEN
!      top layers 
       psnowdz(1) = min(zdz1,psnow/inlvls) 
       psnowdz(2) = min(zdz2,psnow/inlvls) 
!      last layers 
       psnowdz(6) = min(zdzn1,psnow/inlvls)
!      remaining snow for remaining layers
       zwork      = psnow - psnowdz(1) - psnowdz(2) - psnowdz(6)
       psnowdz(3) = zwork*zsgcoef(1)
       psnowdz(4) = zwork*zsgcoef(2)
       psnowdz(5) = zwork*zsgcoef(3)
!      layer 3 tickness >= layer 2 tickness
       zwork=min(0.0,psnowdz(3)-psnowdz(2))
       psnowdz(3)=psnowdz(3)-zwork
       psnowdz(4)=psnowdz(4)+zwork
!      layer 5 tickness >= layer 6 tickness  
       zwork=min(0.0,psnowdz(5)-psnowdz(6))
       psnowdz(5)=psnowdz(5)-zwork
       psnowdz(4)=psnowdz(4)+zwork
  ENDIF
!
! 3. Calculate current grid for 9-layer :
! ---------------------------------------------------------------
!
ELSEIF(inlvls == 9)THEN
!
! critere a satisfaire pour remaillage
!
  IF(present(psnowdz_old))THEN
    gregrid    = psnowdz_old(1) < zcoef1 * min(zdz1 ,psnow/inlvls) .OR. &
               & psnowdz_old(1) > zcoef2 * min(zdz1 ,psnow/inlvls) .OR. &
               & psnowdz_old(2) < zcoef1 * min(zdz2 ,psnow/inlvls) .OR. &
               & psnowdz_old(2) > zcoef2 * min(zdz2 ,psnow/inlvls) .OR. &
               & psnowdz_old(9) < zcoef1 * min(zdzn0,psnow/inlvls) .OR. &
               & psnowdz_old(9) > zcoef2 * min(zdzn0,psnow/inlvls) 
  ENDIF
!             
  IF(gregrid)THEN
!      top layers 
       psnowdz(1) = min(zdz1,psnow/inlvls) 
       psnowdz(2) = min(zdz2,psnow/inlvls) 
       psnowdz(3) = min(zdz3,psnow/inlvls)
!      last layers 
       psnowdz(9)= min(zdzn0,psnow/inlvls)
       psnowdz(8)= min(zdzn1,psnow/inlvls)
       psnowdz(7)= min(zdzn2,psnow/inlvls)
!      remaining snow for remaining layers
       zwork = psnow - psnowdz( 1) - psnowdz( 2) - psnowdz( 3) &
                     - psnowdz( 7) - psnowdz( 8) - psnowdz( 9)
       psnowdz(4) = zwork*zsgcoef(1)
       psnowdz(5) = zwork*zsgcoef(2)
       psnowdz(6) = zwork*zsgcoef(3)
!      layer 4 tickness >= layer 3 tickness
       zwork=min(0.0,psnowdz(4)-psnowdz(3))
       psnowdz(4)=psnowdz(4)-zwork
       psnowdz(5)=psnowdz(5)+zwork
!      layer 6 tickness >= layer 7 tickness  
       zwork=min(0.0,psnowdz(6)-psnowdz(7))
       psnowdz(6)=psnowdz(6)-zwork
       psnowdz(5)=psnowdz(5)+zwork
  ENDIF
!
! 4. Calculate current grid for 12-layer :
! ---------------------------------------------------------------
!
ELSEIF(inlvls == 12)THEN
!
! modif_EB pour maillage
! critere a satisfaire pour remaillage
  IF(present(psnowdz_old))THEN
    gregrid    = psnowdz_old(1)  < zcoef1 * min(zdz1 ,psnow/inlvls) .OR. &
               & psnowdz_old(1)  > zcoef2 * min(zdz1 ,psnow/inlvls) .OR. &
               & psnowdz_old(2)  < zcoef1 * min(zdz2 ,psnow/inlvls) .OR. &
               & psnowdz_old(2)  > zcoef2 * min(zdz2 ,psnow/inlvls) .OR. &
               & psnowdz_old(12) < zcoef1 * min(zdzn0,psnow/inlvls) .OR. &
               & psnowdz_old(12) > zcoef2 * min(zdzn0,psnow/inlvls) 
  ENDIF
!
  IF (gregrid)THEN
!    top layers 
     psnowdz(1) = min(zdz1,psnow/inlvls) 
     psnowdz(2) = min(zdz2,psnow/inlvls) 
     psnowdz(3) = min(zdz3,psnow/inlvls)
     psnowdz(4) = min(zdz4,psnow/inlvls)
     psnowdz(5) = min(zdz5,psnow/inlvls)
!    last layers 
     psnowdz(12)= min(zdzn0,psnow/inlvls)
     psnowdz(11)= min(zdzn1,psnow/inlvls)
     psnowdz(10)= min(zdzn2,psnow/inlvls)
     psnowdz( 9)= min(zdzn3,psnow/inlvls)
!    remaining snow for remaining layers
     zwork = psnow - psnowdz( 1) - psnowdz( 2) - psnowdz( 3) &
                   - psnowdz( 4) - psnowdz( 5) - psnowdz( 9) &
                   - psnowdz(10) - psnowdz(11) - psnowdz(12)
     psnowdz(6) = zwork*zsgcoef(1)
     psnowdz(7) = zwork*zsgcoef(2)
     psnowdz(8) = zwork*zsgcoef(3)
!    layer 6 tickness >= layer 5 tickness
     zwork=min(0.0,psnowdz(6)-psnowdz(5))
     psnowdz(6)=psnowdz(6)-zwork
     psnowdz(7)=psnowdz(7)+zwork 
!    layer 8 tickness >= layer 9 tickness  
     zwork=min(0.0,psnowdz(8)-psnowdz(9))
     psnowdz(8)=psnowdz(8)-zwork
     psnowdz(7)=psnowdz(7)+zwork
  ENDIF
!
! 4. Calculate other non-optimized grid to allow CROCUS PREP :
! ------------------------------------------------------------
! 
ELSE IF(inlvls<10.AND.inlvls/=3.AND.inlvls/=6.AND.inlvls/=9) THEN
!        
  DO jj=1,inlvls
     psnowdz(jj)  = psnow/inlvls
  ENDDO
!
  psnowdz(inlvls) = psnowdz(inlvls) + (psnowdz(1) - min(0.05, psnowdz(1)))
  psnowdz(1)      = min(0.05, psnowdz(1))
!
ELSE   
!
  IF(present(psnowdz_old))THEN
    gregrid    = psnowdz_old(     1) < zcoef1 * min(zdz1      ,psnow/inlvls) .OR. &
               & psnowdz_old(     1) > zcoef2 * min(zdz1      ,psnow/inlvls) .OR. &
               & psnowdz_old(     2) < zcoef1 * min(zdz2      ,psnow/inlvls) .OR. &
               & psnowdz_old(     2) > zcoef2 * min(zdz2      ,psnow/inlvls) .OR. &
               & psnowdz_old(inlvls) < zcoef1 * min(0.05*psnow,psnow/inlvls) .OR. &
               & psnowdz_old(inlvls) > zcoef2 * min(0.05*psnow,psnow/inlvls) 
  ENDIF
!
  IF (gregrid)THEN
     psnowdz(     1) = min(zdz1      ,psnow/inlvls) 
     psnowdz(     2) = min(zdz2      ,psnow/inlvls)
     psnowdz(     3) = min(zdz3      ,psnow/inlvls)
     psnowdz(     4) = min(zdz4      ,psnow/inlvls)
     psnowdz(     5) = min(zdz5      ,psnow/inlvls)
     psnowdz(inlvls) = min(0.05*psnow,psnow/inlvls) 
     zwork = sum(psnowdz(1:5))          
     DO jj=6,inlvls-1,1
        psnowdz(jj) = (psnow - zwork -psnowdz(inlvls))/(inlvls-6) 
     END DO
  ENDIF
!
ENDIF
!
DO jj=1,inlvls
  IF(psnow==xundef)THEN
    psnowdz(jj) = xundef
  ELSEIF(.NOT.gregrid)THEN
    psnowdz(jj) = psnowdz_old(jj)
  ENDIF
END DO
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LGRID_1D',1,zhook_handle)
!
END SUBROUTINE snow3lgrid_1d
!
!###################################################################################
!###################################################################################
SUBROUTINE snow3lagreg(PSNOWDZN,PSNOWDZ,PSNOWRHO,PSNOWGRAN1, PSNOWGRAN2, &
                       psnowhist,psnowgran1n,psnowgran2n,psnowhistn,     &
                       kl1,kl2,psnowddz                        ) 
!
USE modd_snow_metamo
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!       0.1 declarations of arguments        
!        
REAL, DIMENSION(:), INTENT(IN)  :: psnowdzn,psnowdz,psnowrho,psnowddz  
!                                                    
REAL, DIMENSION(:), INTENT(IN)  :: psnowgran1,psnowgran2,psnowhist 
REAL, DIMENSION(:), INTENT(OUT) :: psnowgran1n,psnowgran2n,psnowhistn                                              
!
INTEGER, INTENT(IN) :: kl1  ! Indice couche de reference (i)
INTEGER, INTENT(IN) :: kl2 ! Indice de la couche (i-1 ou i+1) dont une 
                                ! partie est aggregee à la couche (i)
!
!       0.2 declaration of local variables
!        
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zsnowrho
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zdiamd,zdiamv,zspherd,zspherv,&
                                     zdiamn,zsphern,zdent 
!
REAL :: zdelta, zcomp
!
INTEGER :: ident, ivieu, il
!
REAL(KIND=JPRB) :: zhook_handle 
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LAGREG',0,zhook_handle)
!
IF( kl1<kl2 ) THEN
  zdelta = 0.0
  il = kl1
ELSE
  zdelta = 1.0
  il = kl2
ENDIF
!
! Mean Properties
!
!       1. History
!
IF ( psnowhist(kl1)/=psnowhist(kl2) ) THEN
  psnowhistn(kl1) = 0.0
ENDIF
!
!       2. New grain types
!
!       2.1 Same grain type
!
IF (  psnowgran1(kl1)*psnowgran1(kl2)>0.0 .OR. &
    ( psnowgran1(kl1)==0.0 .AND. psnowgran1(kl2)>=0.0 ) .OR. &
    ( psnowgran1(kl2)==0.0 .AND. psnowgran1(kl1)>=0.0 ) ) THEN 
  !
  !code original vincent          PSNOWGRAN1N(KL1)=(PSNOWGRAN1(KL1)*PSNOWRHO(KL1)&
  !code original vincent        *(PSNOWDZN(KL1)-(1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))-ZDELTA*&
  !code original vincent        ABS(PSNOWDDZ(KL2)))+PSNOWGRAN1(KL2)*                   &
  !code original vincent        PSNOWRHO(KL2)*((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+        &   
  !code original vincent        ZDELTA*ABS(PSNOWDDZ(KL2))))/((PSNOWDZN(KL1)-(1.0-ZDELTA)&
  !code original vincent        *ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2)))*        &
  !code original vincent        PSNOWRHO(KL1)+((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+        &   
  !code original vincent        ZDELTA*ABS(PSNOWDDZ(KL2)))*PSNOWRHO(KL2))
  !code original vincent !
  !code original vincent          PSNOWGRAN2N(KL1)=(PSNOWGRAN2(KL1)*PSNOWRHO(KL1) &
  !code original vincent        *(PSNOWDZN(KL1)-(1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))-ZDELTA* &
  !code original vincent        ABS(PSNOWDDZ(KL2)))+PSNOWGRAN2(KL2)*                   &
  !code original vincent        PSNOWRHO(KL2)*((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))        &
  !code original vincent        +ZDELTA*ABS(PSNOWDDZ(KL2))))/((PSNOWDZN(KL1)-(1.0-ZDELTA)&
  !code original vincent        *ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2)))*        &
  !code original vincent        PSNOWRHO(KL1)+((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+        &   
  !code original vincent        ZDELTA*ABS(PSNOWDDZ(KL2)))*PSNOWRHO(KL2))
  !     
  !plm
  CALL get_agreg(kl1,kl2,psnowgran1(kl1),psnowgran1(kl2),psnowgran1n(kl1))
  !
  CALL get_agreg(kl1,kl2,psnowgran2(kl1),psnowgran2(kl2),psnowgran2n(kl1))
  !
  !plm
  !     
ELSE
  !
  !       2.2 Different types
  !        
  IF ( psnowgran1(kl1)<0.0 ) THEN
    ident = kl1
    ivieu = kl2
  ELSE
    ident = kl2
    ivieu = kl1
  ENDIF  
  !                        
  zdiamd(kl1) = - psnowgran1(ident)/xgran * xdiaet + ( 1.0 + psnowgran1(ident)/xgran ) * &
                 ( psnowgran2(ident)/xgran * xdiagf + ( 1.0 - psnowgran2(ident)/xgran ) * xdiafp )
  !
  zspherd(kl1) = psnowgran2(ident)/xgran                
  zdiamv(kl1) = psnowgran2(ivieu)
  zspherv(kl1) = psnowgran1(ivieu)/xgran
  !IF(KL1==1)THEN
  !write(*,*) 'ZDD1',ZDIAMD(1),'ZSD1',ZSPHERD(1)
  !write(*,*) 'ZDV1',ZDIAMV(1),'ZSV1',ZSPHERV(1)
  !ENDIF       
  !
  IF ( ident==kl1 ) THEN
    !code original vincent        ZDIAMN(KL1)= (ZDIAMD(KL1)*PSNOWRHO(IDENT)*&
    !code original vincent            (PSNOWDZN(IDENT)-(1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))-ZDELTA*      &
    !code original vincent            ABS(PSNOWDDZ(KL2)))+ZDIAMV(KL1)*PSNOWRHO(IVIEU)*(       &
    !code original vincent            (1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2))))/&
    !code original vincent            ((PSNOWDZN(KL1)-(1.0-ZDELTA)*                                    &
    !code original vincent            ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2)))*            &
    !code original vincent            PSNOWRHO(KL1)+((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+          &   
    !code original vincent            ZDELTA*ABS(PSNOWDDZ(KL2)))*PSNOWRHO(KL2))
    !
    !plm
    CALL get_agreg(ident,ivieu,zdiamd(kl1),zdiamv(kl1),zdiamn(kl1))
    !
    !plm
    !         
    !code original vincent        ZSPHERN(KL1)= (ZSPHERD(KL1)*PSNOWRHO(IDENT)*&
    !code original vincent            (PSNOWDZN(IDENT)-(1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))-ZDELTA*      &
    !code original vincent            ABS(PSNOWDDZ(KL2)))+ZSPHERV(KL1)*PSNOWRHO(IVIEU)*(       &
    !code original vincent            (1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2))))/&
    !code original vincent            ((PSNOWDZN(KL1)-(1.0-ZDELTA)*                                    &
    !code original vincent            ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2)))*            &
    !code original vincent            PSNOWRHO(KL1)+((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+          &   
    !code original vincent            ZDELTA*ABS(PSNOWDDZ(KL2)))*PSNOWRHO(KL2))

    !plm
    CALL get_agreg(ident,ivieu,zspherd(kl1),zspherv(kl1),zsphern(kl1))   
    !plm
    !
  ELSE
    !code original vincent        ZDIAMN(KL1)= (ZDIAMD(KL1)*PSNOWRHO(IDENT)*&
    !code original vincent            ((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+ZDELTA*ABS(PSNOWDDZ(KL2)))&
    !code original vincent            +ZDIAMV(KL1)*PSNOWRHO(IVIEU)*(PSNOWDZN(IVIEU)-(1.0-ZDELTA)*  & 
    !code original vincent            ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2))))/&
    !code original vincent            ((PSNOWDZN(KL1)-(1.0-ZDELTA)*                          &
    !code original vincent            ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2)))*    &
    !code original vincent            PSNOWRHO(KL1)+((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+   &
    !code original vincent            ZDELTA*ABS(PSNOWDDZ(KL2)))*PSNOWRHO(KL2))
    !code original vincent!            
    !code original vincent         ZSPHERN(KL1)= (ZSPHERD(KL1)*PSNOWRHO(IDENT)*&
    !code original vincent            ((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+ZDELTA*ABS(PSNOWDDZ(KL2)))&
    !code original vincent            +ZSPHERV(KL1)*PSNOWRHO(IVIEU)*(PSNOWDZN(IVIEU)-(1.0-ZDELTA)* & 
    !code original vincent            ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2))))/&
    !code original vincent            ((PSNOWDZN(KL1)-(1.0-ZDELTA)*                                    &
    !code original vincent            ABS(PSNOWDDZ(KL1))-ZDELTA*ABS(PSNOWDDZ(KL2)))*            &
    !code original vincent            PSNOWRHO(KL1)+((1.0-ZDELTA)*ABS(PSNOWDDZ(KL1))+          &   
    !code original vincent            ZDELTA*ABS(PSNOWDDZ(KL2)))*PSNOWRHO(KL2))
    !plm
    !
    CALL get_agreg(ivieu,ident,zdiamv(kl1),zdiamd(kl1),zdiamn(kl1))
    !           
    CALL get_agreg(ivieu,ident,zspherv(kl1),zspherd(kl1),zsphern(kl1))
    !plm
    !
  ENDIF
  !       
  zcomp = zsphern(kl1) * xdiagf + ( 1.-zsphern(kl1) ) * xdiafp
  IF( zdiamn(kl1) < zcomp ) THEN
    !
    zdent(kl1) = ( zdiamn(kl1) - zcomp ) / ( xdiaet - zcomp ) 
    !IF(KL1==1) write(*,*) 'ZDENT',ZDENT(1)   
    psnowgran1n(kl1) = - xgran * zdent(kl1)
    psnowgran2n(kl1) =   xgran * zsphern(kl1)
    !
  ELSE
    !
    psnowgran1n(kl1) = xgran * zsphern(kl1)
    psnowgran2n(kl1) = zdiamn(kl1)
    !
  ENDIF
  !  
ENDIF
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LAGREG',1,zhook_handle)
!
!       3. Update snow grains parameters : GRAN1, GRAN2        
!        PSNOWGRAN1(KL1)=ZSNOWGRAN1(KL1)
!        PSNOWGRAN2(KL1)=ZSNOWGRAN2(KL1)
!
 CONTAINS
!
SUBROUTINE get_agreg(KID1,KID2,PFIELD1,PFIELD2,PFIELD)
!
INTEGER, INTENT(IN) :: kid1, kid2
REAL, INTENT(IN) :: pfield1, pfield2
REAL, INTENT(OUT) :: pfield
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LAGREG:GET_AGREG',0,zhook_handle)
!
pfield = ( pfield1 * psnowrho(kid1) * ( psnowdzn(kid1) - abs(psnowddz(il)) ) &
         + pfield2 * psnowrho(kid2) * abs(psnowddz(il))                        ) / &
         (           psnowrho(kl1) * ( psnowdzn(kl1) - abs(psnowddz(il)) ) + &
                     psnowrho(kl2) * abs(psnowddz(il))                        ) 
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LAGREG:GET_AGREG',1,zhook_handle)
!
END SUBROUTINE get_agreg
!
END SUBROUTINE snow3lagreg    
!###############################################################################    
!###############################################################################
!
!       
!ajout EB : ajout des arguments "N" pour faire idem variables d'origine
SUBROUTINE snow3lavgrain(PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,                &
                         psnowgran1n,psnowgran2n,psnowhistn,pndent,pnvieu,&
                         hsnowmetamo) 
!
USE modd_snow_metamo, ONLY : xvdiam6, xuepsi
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!       0.1 declarations of arguments        
!        
REAL, DIMENSION(:,:), INTENT(INOUT)  :: psnowgran1,psnowgran2,psnowhist 
! ajout EB
REAL, DIMENSION(:,:), INTENT(INOUT)  :: psnowgran1n,psnowgran2n,psnowhistn 
!
REAL, DIMENSION(:), INTENT(IN)    :: pndent, pnvieu          
!
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo
!       0.2 declaration of local variables
!
REAL, DIMENSION(SIZE(PSNOWGRAN1,1)) :: zgran1, zgran2, zhist 
!
LOGICAL, DIMENSION(SIZE(PSNOWGRAN1,1),SIZE(PSNOWGRAN1,2)) :: gdendritic
!
INTEGER :: ji, jl
INTEGER :: inlvls, ini
!
REAL(KIND=JPRB) :: zhook_handle
!         
!       0.3 initialization         
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LAVGRAIN',0,zhook_handle)
!
inlvls    = SIZE(psnowgran1,2) 
ini       = SIZE(psnowgran1,1)
!
zgran1(:) = 0.0
zgran2(:) = 0.0
zhist(:) = 0.0
!     
DO ji = 1,ini
  !
  IF ( pndent(ji)==0.0 .AND. pnvieu(ji)==0.0 ) THEN
    !
    zgran1(ji) = 1.0
    zgran2(ji) = 1.0 
    zhist(ji) = 1.0
    !
  ELSE
    !
    DO jl = 1,inlvls
      IF ( hsnowmetamo=='B92' ) THEN
        gdendritic(ji,jl) = ( psnowgran1(ji,jl) < 0.0 )
      ELSE
        gdendritic(ji,jl) = ( psnowgran1(ji,jl) < xvdiam6*(4.-psnowgran2(ji,jl)) - xuepsi )
      ENDIF
    ENDDO
    !
    IF ( pndent(ji)>=pnvieu(ji) ) THEN      ! more dendritic than non dendritic snow layer
      !
      DO jl = 1,inlvls
        IF ( gdendritic(ji,jl) ) THEN
          zgran1(ji) = zgran1(ji) + psnowgran1(ji,jl)
          zgran2(ji) = zgran2(ji) + psnowgran2(ji,jl)
        ENDIF
      ENDDO
      !
      psnowgran1n(ji,:) = zgran1(ji) / pndent(ji)   
      psnowgran2n(ji,:) = zgran2(ji) / pndent(ji)
      psnowhistn(ji,:) = 0.0
      !
    ELSE                              ! more non dendritic than dendritic snow layers  
      !
      DO jl = 1,inlvls
        IF ( .NOT.gdendritic(ji,jl) ) THEN
          zgran1(ji) = zgran1(ji) + psnowgran1(ji,jl)
          zgran2(ji) = zgran2(ji) + psnowgran2(ji,jl)
          zhist(ji) = zhist(ji) + psnowhist(ji,jl) 
        ENDIF
      ENDDO
      !
      psnowgran1n(ji,:) = zgran1(ji) / pnvieu(ji)
      psnowgran2n(ji,:) = zgran2(ji) / pnvieu(ji)
      psnowhistn(ji,:) = zhist(ji) / pnvieu(ji)
      !    
    ENDIF
    !
  ENDIF
  !
ENDDO


 
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LAVGRAIN',1,zhook_handle)
!
END SUBROUTINE snow3lavgrain         
!        
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow3ldiftyp(PGRAIN1,PGRAIN2,PGRAIN3,PGRAIN4,HSNOWMETAMO) RESULT(ZDIFTYPE)
!
! à remplacer sans doute par une routine equivalente du nouveau crocus
!*    CALCUL DE LA DIFFERENCE ENTRE DEUX TYPES DE GRAINS
!     VALEUR ENTRE 200 ET 0
!
USE modd_snow_metamo, ONLY : xgran, xvdiam6, xuepsi
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!*      0.1    declarations of arguments
REAL, INTENT(IN) :: pgrain1, pgrain2, pgrain3, pgrain4
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo
REAL :: zdiftype, zcoef3, zcoef4
REAL(KIND=JPRB) :: zhook_handle

!*      0.2    calcul de la difference entre type de grains
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDIFTYP',0,zhook_handle)
!
IF ( hsnowmetamo=='B92' ) THEN 
  !
  IF ( ( pgrain1<0.  .AND. pgrain2>=0.) .OR. ( pgrain1>=0. .AND. pgrain2<0. ) ) THEN
    zdiftype = 200.
  ELSEIF ( pgrain1<0. ) THEN
    zdiftype = abs( pgrain1-pgrain2 ) * .5 + abs( pgrain3-pgrain4 ) * .5
  ELSE
    zdiftype = abs( pgrain1-pgrain2 )      + abs( pgrain3-pgrain4 ) * 5. * 10000.
  ENDIF
  !
ELSE
  !
  zcoef3 = xvdiam6 * (4.-pgrain3) - xuepsi
  zcoef4 = xvdiam6 * (4.-pgrain4) - xuepsi 
  IF ( ( pgrain1<zcoef3 .AND. pgrain2>=zcoef4 ) .OR. ( pgrain1>=zcoef3 .AND. pgrain2<zcoef4 ) ) THEN
    zdiftype = 200.
  ELSEIF ( pgrain1<zcoef3 ) THEN
    zdiftype = abs( (pgrain3-pgrain4)*xgran ) * .5 + &
               abs( ( (pgrain1/xvdiam6 - 4. + pgrain3) / (pgrain3 - 3.) - &
                      (pgrain2/xvdiam6 - 4. + pgrain4) / (pgrain4 - 3.) ) * xgran ) * .5
             
  ELSE
    zdiftype = abs( (pgrain3-pgrain4)*xgran )      + abs( zcoef3-zcoef4 ) * 5. * 10000.
  ENDIF  
  !
ENDIF
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDIFTYP',1,zhook_handle)
!
END FUNCTION snow3ldiftyp

!####################################################################
SUBROUTINE get_mass_heat(KJ,KNLVLS_NEW,KNLVLS_OLD,                                &
                         psnowztop_old,psnowztop_new,psnowzbot_old,psnowzbot_new, &
                         psnowrhoo,psnowdzo,psnowgran1o,psnowgran2o,psnowhisto,   &
                         psnowageo,psnowheato,                                    &
                         psnowrhon,psnowdzn,psnowgran1n,psnowgran2n,psnowhistn,   &
                         psnowagen, psnowheatn,hsnowmetamo                        )
!
USE modd_snow_par, ONLY : xsnowcritd, xd1, xd2, xd3, xx, xvalb5, xvalb6
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
INTEGER, INTENt(IN) :: kj
INTEGER, INTENT(IN) :: knlvls_new, knlvls_old
REAL, DIMENSION(:), INTENT(IN) :: psnowztop_old, psnowzbot_old
REAL, DIMENSION(:), INTENT(IN) :: psnowztop_new, psnowzbot_new
REAL, DIMENSION(:), INTENT(IN) :: psnowrhoo, psnowdzo, psnowgran1o, psnowgran2o, &
                                  psnowhisto, psnowageo, psnowheato
REAL, DIMENSION(:), INTENT(IN) :: psnowdzn
 CHARACTER(3), INTENT(IN)       :: hsnowmetamo
REAL, DIMENSION(:), INTENT(OUT) :: psnowrhon, psnowgran1n, psnowgran2n, &
                                   psnowhistn, psnowagen, psnowheatn
!
REAL :: zpropor, zmasdz_old, zdiam, zmastot_t07
REAL :: zsnowhean, zmastotn, zdentmoyn, zsphermoyn, zalbmoyn, zhistmoyn
REAL :: zagemoyn
!
INTEGER :: jst_new, jst_old
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('GET_MASS_HEAT',0,zhook_handle)
!
psnowrhon(:) = 0.
psnowgran1n(:) = 0.
psnowgran2n(:) = 0.
psnowhistn(:) = 0.
psnowagen(:) = 0.
psnowheatn(:) = 0.
!
DO jst_new = 1,knlvls_new
  !
  zsnowhean   = 0.
  zmastotn    = 0.
  zmastot_t07 = 0.
  zdentmoyn   = 0.
  zsphermoyn  = 0.
  zalbmoyn    = 0.
  zdiam       = 0.
  zhistmoyn   = 0.
  zagemoyn    = 0.
  !
  ! lopp over the ols snow layers 
  DO jst_old = 1,knlvls_old
    !
    IF( psnowztop_old(jst_old)<=psnowzbot_new(jst_new) ) THEN
      ! JST_OLD lower than JJ_NEW ==> no contribution
    ELSEIF ( psnowzbot_old(jst_old)>=psnowztop_new(jst_new) ) THEN
      ! JST_OLD higher than JJ_NEW ==> no contribution
    ELSE
      ! old layer contributes to the new one
      ! computation of its contributing ratio and mass/heat 
      zpropor = ( min( psnowztop_old(jst_old), psnowztop_new(jst_new) )   &
                - max( psnowzbot_old(jst_old), psnowzbot_new(jst_new) ) ) &
                 / psnowdzo(jst_old) 
      zmasdz_old = zpropor * psnowrhoo(jst_old) * psnowdzo(jst_old)
      !
      zmastotn    = zmastotn + zmasdz_old
      zmastot_t07 = zmastot_t07 + 1.
      !
      zsnowhean = zsnowhean + zpropor * psnowheato(jst_old)
      !
      IF ( hsnowmetamo=='B92' ) THEN
        !
        ! contribution to the grain optical size and then to the albedo
        IF ( psnowgran1o(jst_old)<0. ) THEN
          zdiam = -psnowgran1o(jst_old)*xd1/xx + (1.+psnowgran1o(jst_old)/xx) * &
                 ( psnowgran2o(jst_old)*xd2/xx + (1.-psnowgran2o(jst_old)/xx)*xd3 ) 
          zdiam = zdiam/10000.      
          zdentmoyn  = zdentmoyn  - zmasdz_old * psnowgran1o(jst_old) / xx
          zsphermoyn = zsphermoyn + zmasdz_old * psnowgran2o(jst_old) / xx
        ELSE
          zdiam = psnowgran2o(jst_old)
          zdentmoyn  = zdentmoyn  + zmasdz_old * 0.
          zsphermoyn = zsphermoyn + zmasdz_old * psnowgran1o(jst_old) / xx
        ENDIF
        !
      ELSE
        !
        zdiam = psnowgran1o(jst_old)
        zsphermoyn = zsphermoyn + zmasdz_old * psnowgran2o(jst_old)
        !
      ENDIF
      !
      zalbmoyn  = zalbmoyn  + max( 0., zmasdz_old * (xvalb5-xvalb6*sqrt(zdiam)) )
      zhistmoyn = zhistmoyn + zmasdz_old * psnowhisto(jst_old)
      zagemoyn  = zagemoyn  + zmasdz_old * psnowageo(jst_old)
      !
    ENDIF
    !
  ENDDO
  ! 
  ! the new layer inherits from the weihted average properties of the old ones
  ! heat and mass
  psnowheatn(jst_new) = zsnowhean
  psnowrhon(jst_new) = zmastotn / psnowdzn(jst_new)
  ! grain type and size decuced from the average albedo
  zalbmoyn   = zalbmoyn / zmastotn
  zsphermoyn = max( 0., zsphermoyn/zmastotn )
  zdentmoyn  = max( 0., zdentmoyn /zmastotn )
  zdiam = ( (xvalb5-zalbmoyn)/xvalb6 )**2
  !
  IF ( hsnowmetamo=='B92' ) THEN
    !
    ! size between D2 and D3 and dendricity < 0         
    ! sphericity is firts preserved, if possible. If not,
    ! denditricity =0
    psnowgran1n(jst_new) = -xx * zdentmoyn
    !
    IF ( zdentmoyn/=1.) THEN
      psnowgran2n(jst_new) = xx * ( ( zdiam*10000. + psnowgran1n(jst_new)*xd1/xx ) &
                                 / ( 1. + psnowgran1n(jst_new)/xx ) - xd3 )        &
                             / ( xd2-xd3 )
    ENDIF
    !
    ! dendricity is preserved if possible and sphericity is adjusted
    IF ( zdiam < xd2/10000. - 0.0000001 ) THEN
      !
      IF ( abs( psnowgran1n(jst_new)+xx ) < 0.01 ) THEN
        !
        psnowgran2n(jst_new) = xx * zsphermoyn
        !
      ELSEIF ( abs( psnowgran1n(jst_new) ) < 0.0001 ) THEN ! dendritic snow
        !
        psnowgran1n(jst_new) = xx * zsphermoyn
        psnowgran2n(jst_new) = zdiam
        !
      ELSEIF ( psnowgran2n(jst_new) < 0. ) THEN ! non dendritic
        !
        psnowgran2n(jst_new) = 0.
        !
      ELSEIF ( psnowgran2n(jst_new) > xx + 0.0000001 ) THEN ! non dendritic
        !
        psnowgran2n(jst_new) = xx
        !
      ENDIF
      !
    ELSEIF ( zdiam > xd3/10000. .OR. zdentmoyn <= 0. + 0.0000001 .OR. &
             psnowgran2n(jst_new) < 0. .OR. psnowgran2n(jst_new) > xx ) THEN
      !
      ! dendritic snow
      ! inconsistency with ZDIAM ==>  dendricity = 0
      ! size between D2 and D3 and dendricity == 0          
      psnowgran1n(jst_new) = xx * zsphermoyn
      psnowgran2n(jst_new) = zdiam
      !
    ENDIF
    !
  ELSE
    !
    psnowgran1n(jst_new) = zdiam
    psnowgran2n(jst_new) = min( 1., zsphermoyn )
    !
  ENDIF
  !
  psnowhistn(jst_new) = nint( zhistmoyn/zmastotn )
  psnowagen(jst_new) = zagemoyn / zmastotn
  !
ENDDO     
!
IF (lhook) CALL dr_hook('GET_MASS_HEAT',1,zhook_handle)
!
END SUBROUTINE get_mass_heat
!####################################################################
SUBROUTINE get_diam(PSNOWGRAN1,PSNOWGRAN2,PDIAM,HSNOWMETAMO)
!
USE modd_snow_par, ONLY : xd1, xd2, xd3, xx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: psnowgran1
REAL, INTENT(IN) :: psnowgran2
REAL, INTENT(OUT) :: pdiam
!
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('GET_DIAM',0,zhook_handle)
!
IF ( hsnowmetamo=='B92' ) THEN
  !
  IF( psnowgran1<0. ) THEN
    pdiam = -psnowgran1*xd1/xx + (1.+psnowgran1/xx) * &
           ( psnowgran2*xd2/xx + (1.-psnowgran2/xx) * xd3 ) 
    pdiam = pdiam/10000.      
  ELSE 
    pdiam = psnowgran2*psnowgran1/xx + &
            max( 0.0004, 0.5*psnowgran2 ) * ( 1.-psnowgran1/xx )      
  ENDIF
  !
ELSE
  !
  pdiam = psnowgran1
  !
ENDIF
!
IF (lhook) CALL dr_hook('GET_DIAM',1,zhook_handle)
!
END SUBROUTINE get_diam
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow3lradabs_0d(PSNOWRHO,PSNOWDZ,PSPECTRALALBEDO,PZENITH,PPERMSNOWFRAC,PDSGRAIN) RESULT(PCOEF)
!
!!    PURPOSE
!!    -------
!     Calculate the transmission of shortwave radiation within the snowpack
!     (with depth)
!     A. Boone 02/2011
!     A. Boone 11/2014 Updated to use spectral dependence.
!                      NOTE, assumes 3 spectral bands
!
USE modd_surf_par, ONLY : xundef
USE modd_meb_par,  ONLY : xsw_wght_vis, xsw_wght_nir
USE modd_snow_par, ONLY : xvspec1,xvspec2,xvspec3,xvbeta1,xvbeta2, &
                          xvbeta4,xvbeta3,xvbeta5, xmincoszen
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL,               INTENT(IN)                   :: psnowrho        ! snow density    (kg m-3)
REAL,               INTENT(IN)                   :: psnowdz         ! layer thickness (m)
REAL,               INTENT(IN)                   :: pzenith         ! zenith angle    (rad)
REAL,               INTENT(IN)                   :: ppermsnowfrac   ! permanent snow fraction (-)
REAL, DIMENSION(:), INTENT(IN)                   :: pspectralalbedo ! spectral albedo (-)
REAL,               INTENT(IN)                   :: pdsgrain        ! Snow optical grain diameter (m)
!
REAL                                             :: pcoef           ! -
!
!*      0.2    declarations of local variables
!
REAL                                             :: zwork, zprojlat,                  &
                                                    zbeta1, zbeta2, zbeta3,           &
                                                    zopticalpath1, zopticalpath2,     &
                                                    zopticalpath3
!
REAL(KIND=JPRB)                                  :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_0D',0,zhook_handle)
!
! Coefficient for taking into account the increase of path length of rays
! in snow due to zenithal angle
!
zprojlat         = (1.0-ppermsnowfrac)+ppermsnowfrac/ &
                   max(xmincoszen,cos(pzenith))
!
! Extinction coefficient:
!
zwork            = sqrt(pdsgrain)
zbeta1           = max(xvbeta1*psnowrho/zwork,xvbeta2)
zbeta2           = max(xvbeta3*psnowrho/zwork,xvbeta4)
zbeta3           = xvbeta5
!
zopticalpath1    = zbeta1*psnowdz
zopticalpath2    = zbeta2*psnowdz
zopticalpath3    = xundef
!
IF(pspectralalbedo(3)==xundef)THEN 
   pcoef         = xsw_wght_vis*(1.0-pspectralalbedo(1))*exp(-zopticalpath1*zprojlat) &
                 + xsw_wght_nir*(1.0-pspectralalbedo(2))*exp(-zopticalpath2*zprojlat) 
ELSE
   zopticalpath3 = zbeta3*psnowdz
   pcoef         = xvspec1*(1.0-pspectralalbedo(1))*exp(-zopticalpath1*zprojlat) &
                 + xvspec2*(1.0-pspectralalbedo(2))*exp(-zopticalpath2*zprojlat) &
                 + xvspec3*(1.0-pspectralalbedo(3))*exp(-zopticalpath3*zprojlat)
ENDIF
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_0D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION snow3lradabs_0d
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow3lradabs_1d(PSNOWRHO,PSNOWDZ,PSPECTRALALBEDO,PZENITH,PPERMSNOWFRAC,PDSGRAIN) RESULT(PCOEF)
!
!!    PURPOSE
!!    -------
!     Calculate the transmission of shortwave radiation within the snowpack
!     (with depth)
!     A. Boone 02/2011
!     A. Boone 11/2014 Updated to use spectral dependence
!                      NOTE, assumes 3 spectral bands
!
USE modd_surf_par, ONLY : xundef
USE modd_meb_par,  ONLY : xsw_wght_vis, xsw_wght_nir
USE modd_snow_par, ONLY : xvspec1,xvspec2,xvspec3,xvbeta1,xvbeta2, &
                          xvbeta4,xvbeta3,xvbeta5, xmincoszen
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:),   INTENT(IN)                 :: psnowrho        ! snow density    (kg m-3)
REAL, DIMENSION(:),   INTENT(IN)                 :: psnowdz         ! layer thickness (m)
REAL, DIMENSION(:),   INTENT(IN)                 :: pzenith         ! zenith angle    (rad)
REAL, DIMENSION(:),   INTENT(IN)                 :: ppermsnowfrac   ! permanent snow fraction (-)
REAL, DIMENSION(:,:), INTENT(IN)                 :: pspectralalbedo ! spectral albedo (-)
REAL, DIMENSION(:),   INTENT(IN)                 :: pdsgrain        ! Snow optical grain diameter (m)
!
REAL, DIMENSION(SIZE(PSNOWRHO))                  :: pcoef           ! -
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO))                  :: zwork, zprojlat,                  &
                                                    zbeta1, zbeta2, zbeta3,           &
                                                    zopticalpath1, zopticalpath2,     &
                                                    zopticalpath3
!
REAL(KIND=JPRB)                                  :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_1D',0,zhook_handle)
!
! Coefficient for taking into account the increase of path length of rays
! in snow due to zenithal angle
!
zprojlat(:)         = (1.0-ppermsnowfrac(:))+ppermsnowfrac(:)/ &
                      max(xmincoszen,cos(pzenith(:)))
!
! Extinction coefficient:
!
zwork(:)            = sqrt(pdsgrain(:))
zbeta1(:)           = max(xvbeta1*psnowrho(:)/zwork(:),xvbeta2)
zbeta2(:)           = max(xvbeta3*psnowrho(:)/zwork(:),xvbeta4)
zbeta3(:)           = xvbeta5
!
zopticalpath1(:)    = zbeta1(:)*psnowdz(:)
zopticalpath2(:)    = zbeta2(:)*psnowdz(:)
zopticalpath3(:)    = xundef
!
WHERE(pspectralalbedo(:,3)==xundef)
   pcoef(:)         = xsw_wght_vis*(1.0-pspectralalbedo(:,1))*exp(-zopticalpath1(:)*zprojlat(:)) &
                    + xsw_wght_nir*(1.0-pspectralalbedo(:,2))*exp(-zopticalpath2(:)*zprojlat(:)) 
ELSEWHERE
   zopticalpath3(:) = zbeta3(:)*psnowdz(:)
   pcoef(:)         = xvspec1*(1.0-pspectralalbedo(:,1))*exp(-zopticalpath1(:)*zprojlat(:)) &
                    + xvspec2*(1.0-pspectralalbedo(:,2))*exp(-zopticalpath2(:)*zprojlat(:)) &
                    + xvspec3*(1.0-pspectralalbedo(:,3))*exp(-zopticalpath3(:)*zprojlat(:))
END WHERE
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_1D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION snow3lradabs_1d
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow3lradabs_2d(PSNOWRHO,PSNOWDZ,PSPECTRALALBEDO,PZENITH,PPERMSNOWFRAC,PDSGRAIN) RESULT(PCOEF)
!
!!    PURPOSE
!!    -------
!     Calculate the transmission of shortwave radiation within the snowpack
!     (with depth)
!     A. Boone 02/2011
!     A. Boone 11/2014 Updated to use spectral dependence
!                      NOTE, assumes 3 spectral bands
!
USE modd_surf_par, ONLY : xundef
USE modd_meb_par,  ONLY : xsw_wght_vis, xsw_wght_nir
USE modd_snow_par, ONLY : xvspec1,xvspec2,xvspec3,xvbeta1,xvbeta2, &
                          xvbeta4,xvbeta3,xvbeta5, xmincoszen
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:),   INTENT(IN)                 :: psnowrho        ! snow density    (kg m-3)
REAL, DIMENSION(:,:),   INTENT(IN)                 :: psnowdz         ! layer thickness (m)
REAL, DIMENSION(:,:),   INTENT(IN)                 :: pzenith         ! zenith angle    (rad)
REAL, DIMENSION(:,:),   INTENT(IN)                 :: ppermsnowfrac   ! permanent snow fraction (-)
REAL, DIMENSION(:,:,:), INTENT(IN)                 :: pspectralalbedo ! spectral albedo (-)
REAL, DIMENSION(:,:),   INTENT(IN)                 :: pdsgrain        ! Snow optical grain diameter (m)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: pcoef           ! -
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zwork, zprojlat,                  &
                                                      zbeta1, zbeta2, zbeta3,           &
                                                      zopticalpath1, zopticalpath2,     &
                                                      zopticalpath3
!
REAL(KIND=JPRB)                                    :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_2D',0,zhook_handle)
!
! Coefficient for taking into account the increase of path length of rays
! in snow due to zenithal angle
!
zprojlat(:,:)         = (1.0-ppermsnowfrac(:,:))+ppermsnowfrac(:,:)/ &
                        max(xmincoszen,cos(pzenith(:,:)))
!
! Extinction coefficient:
!
zwork(:,:)            = sqrt(pdsgrain(:,:))
zbeta1(:,:)           = max(xvbeta1*psnowrho(:,:)/zwork(:,:),xvbeta2)
zbeta2(:,:)           = max(xvbeta3*psnowrho(:,:)/zwork(:,:),xvbeta4)
zbeta3(:,:)           = xvbeta5
!
zopticalpath1(:,:)    = zbeta1(:,:)*psnowdz(:,:)
zopticalpath2(:,:)    = zbeta2(:,:)*psnowdz(:,:)
zopticalpath3(:,:)    = xundef
!
WHERE(pspectralalbedo(:,:,3)==xundef)
   pcoef(:,:)         = xsw_wght_vis*(1.0-pspectralalbedo(:,:,1))*exp(-zopticalpath1(:,:)*zprojlat(:,:)) &
                      + xsw_wght_nir*(1.0-pspectralalbedo(:,:,2))*exp(-zopticalpath2(:,:)*zprojlat(:,:)) 
ELSEWHERE
   zopticalpath3(:,:) = zbeta3(:,:)*psnowdz(:,:)
   pcoef(:,:)         = xvspec1*(1.0-pspectralalbedo(:,:,1))*exp(-zopticalpath1(:,:)*zprojlat(:,:)) &
                      + xvspec2*(1.0-pspectralalbedo(:,:,2))*exp(-zopticalpath2(:,:)*zprojlat(:,:)) &
                      + xvspec3*(1.0-pspectralalbedo(:,:,3))*exp(-zopticalpath3(:,:)*zprojlat(:,:))
END WHERE
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_2D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION snow3lradabs_2d
!####################################################################
!####################################################################
!####################################################################
      SUBROUTINE snow3lthrm(PSNOWRHO,PSCOND,PSNOWTEMP,PPS)
!
!!    PURPOSE
!!    -------
!     Calculate snow thermal conductivity from
!     Sun et al. 1999, J. of Geophys. Res., 104, 19587-19579 (vapor) 
!     and Yen, 1981, CRREL Rep 81-10 (snow)
!     or Anderson, 1976, NOAA Tech. Rep. NWS 19 (snow).
!
!
USE modd_csts,     ONLY : xp00, xcondi, xrholw
!
USE modd_snow_par, ONLY : xvrkz6, xsnowthrmcond1, &
                          xsnowthrmcond2,         &
                          xsnowthrmcond_avap,     &
                          xsnowthrmcond_bvap,     &
                          xsnowthrmcond_cvap 
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)      :: pps
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowtemp, psnowrho
!
REAL, DIMENSION(:,:), INTENT(OUT)   :: pscond
!
!
!*      0.2    declarations of local variables
!
INTEGER                              :: jj, ji
!
INTEGER                              :: ini
INTEGER                              :: inlvls
!
 CHARACTER(LEN=5)                     :: ysnowcond !should be in namelist
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LTHRM',0,zhook_handle)
!
ini    = SIZE(psnowrho(:,:),1)
inlvls = SIZE(psnowrho(:,:),2)
!
! 1. Snow thermal conductivity
! ----------------------------
!
ysnowcond='YEN81' !should be in namelist
!
IF(ysnowcond=='AND76')THEN
!  Thermal conductivity coefficients from Anderson (1976)
  pscond(:,:) = (xsnowthrmcond1 + xsnowthrmcond2*psnowrho(:,:)*psnowrho(:,:))
ELSE
! Thermal conductivity coefficients from Yen (1981)
  pscond(:,:) = xcondi * exp(xvrkz6*log(psnowrho(:,:)/xrholw))
ENDIF
!
! 2. Implicit vapor diffn effects
! -------------------------------
!
DO jj=1,inlvls
   DO ji=1,ini
    pscond(ji,jj) = pscond(ji,jj) + max(0.0,(xsnowthrmcond_avap+(xsnowthrmcond_bvap/(psnowtemp(ji,jj) &
                                  + xsnowthrmcond_cvap)))*(xp00/pps(ji)))
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LTHRM',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE snow3lthrm
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow3ldopt_2d(PSNOWRHO,PSNOWAGE) RESULT(PDOPT)
!
!!    PURPOSE
!!    -------
!     Calculate the optical grain diameter.
!
USE modd_snow_par, ONLY : xdsgrain_max,xsnow_agrain, & 
                          xsnow_bgrain,xsnow_cgrain
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)                   :: psnowrho,psnowage
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: pdopt
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zage
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsrho4
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDOPT_2D',0,zhook_handle)
!
zage(:,:) = min(15.,psnowage(:,:))
!
zsrho4(:,:) = psnowrho(:,:)*psnowrho(:,:)*psnowrho(:,:)*psnowrho(:,:)
!
pdopt(:,:) = min(xdsgrain_max,xsnow_agrain+xsnow_bgrain*zsrho4(:,:)+xsnow_cgrain*zage(:,:))
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDOPT_2D',1,zhook_handle)
!
END FUNCTION snow3ldopt_2d
!####################################################################
FUNCTION snow3ldopt_1d(PSNOWRHO,PSNOWAGE) RESULT(PDOPT)
!
!!    PURPOSE
!!    -------
!     Calculate the optical grain diameter.
!
USE modd_snow_par, ONLY : xdsgrain_max,xsnow_agrain, & 
                          xsnow_bgrain,xsnow_cgrain
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)  :: psnowrho,psnowage
!
REAL, DIMENSION(SIZE(PSNOWRHO)) :: pdopt
REAL, DIMENSION(SIZE(PSNOWRHO)) :: zage
REAL, DIMENSION(SIZE(PSNOWRHO)) :: zsrho4
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDOPT_1D',0,zhook_handle)
!
zage(:) = min(15.,psnowage(:))
!
zsrho4(:) = psnowrho(:)*psnowrho(:)*psnowrho(:)*psnowrho(:)
!
pdopt(:) = min(xdsgrain_max,xsnow_agrain+xsnow_bgrain*zsrho4(:)+xsnow_cgrain*zage(:))
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDOPT_1D',1,zhook_handle)
!
END FUNCTION snow3ldopt_1d
!####################################################################
FUNCTION snow3ldopt_0d(PSNOWRHO,PSNOWAGE) RESULT(PDOPT)
!
!!    PURPOSE
!!    -------
!     Calculate the optical grain diameter.
!
USE modd_snow_par, ONLY : xdsgrain_max,xsnow_agrain, & 
                          xsnow_bgrain,xsnow_cgrain
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)  :: psnowrho,psnowage
!
REAL :: pdopt
REAL :: zage
REAL :: zsrho4
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDOPT_0D',0,zhook_handle)
!
zage = min(15.,psnowage)
!
zsrho4 = psnowrho*psnowrho*psnowrho*psnowrho
!
pdopt = min(xdsgrain_max,xsnow_agrain+xsnow_bgrain*zsrho4+xsnow_cgrain*zage)
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LDOPT_0D',1,zhook_handle)
!
END FUNCTION snow3ldopt_0d
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snow3lalb(PALBEDOSC,PSPECTRALALBEDO,PSNOWRHO,PSNOWAGE,   &
                     ppermsnowfrac,pps)  
!
!!    PURPOSE
!!    -------
!     Calculate the snow surface albedo. Use the method of
!     CROCUS with 3 spectral albedo depending on snow density 
!     and age
!
!
USE modd_snow_par, ONLY : xvaging_glacier, xvaging_noglacier,     &
                          xvalb2,xvalb3,xvalb4,xvalb5,xvalb6,     &
                          xvalb7,xvalb8,xvalb9,xvalb10,xvalb11,   &
                          xvdiop1,xvrpre1,xvrpre2,xvpres1,        &
                          xvw1,xvw2,xvspec1,xvspec2,xvspec3
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)      :: psnowrho
REAL, DIMENSION(:), INTENT(IN)      :: psnowage
REAL, DIMENSION(:), INTENT(IN)      :: ppermsnowfrac
REAL, DIMENSION(:), INTENT(IN)      :: pps
!
REAL, DIMENSION(:),   INTENT(INOUT) :: palbedosc
REAL, DIMENSION(:,:), INTENT(INOUT) :: pspectralalbedo
!
!*      0.2    declarations of local variables
!
REAL, PARAMETER                 :: zalbnir1 = 0.3
REAL, PARAMETER                 :: zalbnir2 = 0.0
!
REAL, DIMENSION(SIZE(PSNOWRHO)) :: zvaging, zdiam, zage,  &
                                   zwork, zpres_effect
!
REAL, DIMENSION(SIZE(PSNOWRHO)) :: zalb1, zalb2, zalb3
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LALB',0,zhook_handle)
!
! 0. Initialize:
! ------------------
!
!Snow age effect parameter for Visible (small over glacier)
zvaging(:)=xvaging_glacier*ppermsnowfrac(:) + xvaging_noglacier*(1.0-ppermsnowfrac(:))
!
!Atm pression effect parameter on albedo
zpres_effect(:) = xvalb10*min(max(pps(:)/xvpres1,xvrpre1),xvrpre2)
!
! 1. Snow optical grain diameter :
! --------------------------------
!
!Snow optical diameter do not depend on snow age over glacier or polar regions
zage(:) = (1.0-ppermsnowfrac(:))*psnowage(:)
!
zdiam(:) = snow3ldopt(psnowrho(:),zage(:))
!
! 2. spectral albedo over 3 bands :
! ---------------------------------
!
!Snow age effect limited to 1 year
zage(:) = min(365.,psnowage(:))
!
zwork(:)=sqrt(zdiam(:))
!
! Visible
zalb1(:)=min(xvalb4,xvalb2-xvalb3*zwork(:))
zalb1(:)=max(xvalb11,zalb1(:)-zpres_effect(:)*zage(:)/zvaging(:))
!
! near Infra-red 1
zalb2(:)=xvalb5-xvalb6*zwork(:)
zalb2(:)=max(zalbnir1,zalb2(:))
!
! near Infra-red 2
zdiam(:)=min(xvdiop1,zdiam(:))
zwork(:)=sqrt(zdiam(:))
zalb3(:)=xvalb7*zdiam(:)-xvalb8*zwork(:)+xvalb9
zalb3(:)=max(zalbnir2,zalb3(:))
!
pspectralalbedo(:,1)=zalb1(:)
pspectralalbedo(:,2)=zalb2(:)
pspectralalbedo(:,3)=zalb3(:)
!
! 3. total albedo :
! -----------------
!
palbedosc(:)=xvspec1*zalb1(:)+xvspec2*zalb2(:)+xvspec3*zalb3(:)
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LALB',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE snow3lalb
!####################################################################
!####################################################################
!####################################################################
END MODULE mode_snow3l