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 snow3lradabs_sfc
  MODULE PROCEDURE snow3lradabs_sfc
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
!
INTERFACE snow3lfall
  MODULE PROCEDURE snow3lfall
END INTERFACE
!
INTERFACE snow3ltransf
  MODULE PROCEDURE snow3ltransf
END INTERFACE
!
INTERFACE snow3lcompactn
  MODULE PROCEDURE snow3lcompactn
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
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow3lradabs_sfc(PSNOWRHO,PSNOWDZ,PSPECTRALALBEDO,PZENITH,PPERMSNOWFRAC,PDSGRAIN) RESULT(PCOEF)
!
!!    PURPOSE
!!    -------
!     Calculate the transmission of shortwave radiation within the snowpack
!     (with depth) for 3 spectral bands
!     A. Boone 02/2011
!     A. Boone 11/2014 Updated to use spectral dependence
!                      NOTE, assumes 3 spectral bands
!     A. Boone 06/2017 ONLY considers albedo of uppermost layer  
!
USE modd_surf_par, ONLY : xundef
USE modd_snow_par, ONLY : xvspec1,xvspec2,xvspec3,xvbeta1,xvbeta2, &
                          xvbeta4,xvbeta3,xvbeta5, xmincoszen
USE modd_meb_par,  ONLY : xsw_wght_vis, xsw_wght_nir
!
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
!
INTEGER                                            :: JJ, JI, INLVLS, INI
REAL, DIMENSION(SIZE(PSNOWRHO,1))                  :: ZPROJLAT,  ZOPTICALPATH1,         &
                                                      ZOPTICALPATH2, ZOPTICALPATH3
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZWORK, ZBETA1, ZBETA2, ZBETA3
!
REAL(KIND=JPRB)                                    :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_SFC',0,zhook_handle)
!
ini    = SIZE(psnowdz(:,:),1)
inlvls = SIZE(psnowdz(:,:),2)
!
! 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(:)      = 0.0
zopticalpath2(:)      = 0.0
zopticalpath3(:)      = 0.0
!
DO jj=1,inlvls
   DO ji=1,ini
      zopticalpath1(ji) = zopticalpath1(ji) + zbeta1(ji,jj)*psnowdz(ji,jj)
      zopticalpath2(ji) = zopticalpath2(ji) + zbeta2(ji,jj)*psnowdz(ji,jj)

      IF(pspectralalbedo(ji,3)==xundef)THEN 

         pcoef(ji,jj)     = xsw_wght_vis*(1.0-pspectralalbedo(ji,1))*exp(-zopticalpath1(ji)*zprojlat(ji)) &
                           + xsw_wght_nir*(1.0-pspectralalbedo(ji,2))*exp(-zopticalpath2(ji)*zprojlat(ji))  
      ELSE
      
         zopticalpath3(ji) = zopticalpath3(ji) + zbeta3(ji,jj)*psnowdz(ji,jj)

         pcoef(ji,jj)     = xvspec1*(1.0-pspectralalbedo(ji,1))*exp(-zopticalpath1(ji)*zprojlat(ji)) &
                           + xvspec2*(1.0-pspectralalbedo(ji,2))*exp(-zopticalpath2(ji)*zprojlat(ji)) &
                           + xvspec3*(1.0-pspectralalbedo(ji,3))*exp(-zopticalpath3(ji)*zprojlat(ji))

      ENDIF
      
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LRADABS_SFC',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION snow3lradabs_sfc
!####################################################################
!####################################################################
!####################################################################
      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
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snow3lfall(PTSTEP,PSR,PTA,PVMOD,PSNOW,PSNOWRHO,PSNOWDZ,        &
                      PSNOWHEAT,PSNOWHMASS,PSNOWAGE,PPERMSNOWFRAC)  
!
!!    PURPOSE
!!    -------
!     Calculate changes to snowpack resulting from snowfall.
!     Update mass and heat content of uppermost layer.
!
!
USE modd_csts,     ONLY : xlmtt, xtt, xci
USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin, &
                          xsnowfall_a_sn,         &
                          xsnowfall_b_sn,         &
                          xsnowfall_c_sn
!                   
USE yomhook   ,    ONLY : lhook,   dr_hook
USE parkind1  ,    ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: PTSTEP
!
REAL, DIMENSION(:), INTENT(IN)      :: PSR, PTA, PVMOD, PPERMSNOWFRAC
!
REAL, DIMENSION(:), INTENT(INOUT)   :: PSNOW
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ, PSNOWHEAT, PSNOWAGE
!
REAL, DIMENSION(:), INTENT(OUT)     :: PSNOWHMASS
!
!
!*      0.2    declarations of local variables
!
INTEGER                             :: JJ, JI
!
INTEGER                             :: INI
INTEGER                             :: INLVLS
!
REAL, DIMENSION(SIZE(PTA))          :: ZSNOWFALL, ZRHOSNEW,        &
                                       ZSNOW, ZSNOWTEMP,           &
                                       ZSNOWFALL_DELTA, ZSCAP,     &
                                       ZAGENEW
!                               
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
! 0. Initialize:
! ------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LFALL',0,zhook_handle)
!
ini             = SIZE(psnowdz(:,:),1)
inlvls          = SIZE(psnowdz(:,:),2)
!
zrhosnew(:)     = xrhosmin_es
zagenew(:)     = 0.0
zsnowfall(:)    = 0.0
zscap(:)        = 0.0
zsnow(:)        = psnow(:)
!
psnowhmass(:)   = 0.0
!
! 1. Incorporate snowfall into snowpack:
! --------------------------------------
!
!
! Heat content of newly fallen snow (J/m2):
! NOTE for now we assume the snowfall has
! the temperature of the snow surface upon reaching the snow.
! This is done as opposed to using the air temperature since
! this flux is quite small and has little to no impact
! on the time scales of interest. If we use the above assumption
! then, then the snowfall advective heat flux is zero.
!
zsnowtemp(:)  = xtt
zscap(:)  = snow3lscap(psnowrho(:,1))
!
WHERE (psr(:) > 0.0 .AND. psnowdz(:,1)>0.)
  zsnowtemp(:)  = xtt + (psnowheat(:,1) +                              &
                    xlmtt*psnowrho(:,1)*psnowdz(:,1))/                   &
                    (zscap(:)*max(xsnowdmin/inlvls,psnowdz(:,1)))  
  zsnowtemp(:)  = min(xtt, zsnowtemp(:))
END WHERE
!
WHERE (psr(:) > 0.0)
!
  psnowhmass(:) = psr(:)*(xci*(zsnowtemp(:)-xtt)-xlmtt)*ptstep
!
! Snowfall density: Following CROCUS (Pahaut 1976)
!
   zrhosnew(:)   = max(xrhosmin_es, xsnowfall_a_sn + xsnowfall_b_sn*(pta(:)-xtt)+         &
                     xsnowfall_c_sn*sqrt(pvmod(:)))  
!
!
! Fresh snowfall changes the snowpack age,
! decreasing in uppermost snow layer (mass weighted average):
!
   psnowage(:,1) = (psnowage(:,1)*psnowdz(:,1)*psnowrho(:,1)+zagenew(:)*psr(:)*ptstep) / &
                   (psnowdz(:,1)*psnowrho(:,1)+psr(:)*ptstep)
!
! Augment total pack depth:
!
   zsnowfall(:)  = psr(:)*ptstep/zrhosnew(:)    ! snowfall thickness (m)
!
   psnow(:)      = psnow(:) + zsnowfall(:)
!
! Fresh snowfall changes the snowpack
! density, increases the total liquid water
! equivalent: in uppermost snow layer:
!
   psnowrho(:,1) = (psnowdz(:,1)*psnowrho(:,1) + zsnowfall(:)*zrhosnew(:))/     &
                   (psnowdz(:,1)+zsnowfall(:))  
!
   psnowdz(:,1)  = psnowdz(:,1) + zsnowfall(:)
!
! Add energy of snowfall to snowpack:
! Update heat content (J/m2) (therefore the snow temperature
! and liquid content):
!
   psnowheat(:,1)  = psnowheat(:,1) + psnowhmass(:)
!
END WHERE
!
!
! 2. Case of new snowfall on a previously snow-free surface:
! ----------------------------------------------------------
!
! When snow first falls on a surface devoid of snow,
! redistribute the snow mass throughout the 3 layers:
! (temperature already set in the calling routine
! for this case)
!
zsnowfall_delta(:)    = 0.0
WHERE(zsnow(:) == 0.0 .AND. psr(:) > 0.0)
   zsnowfall_delta(:) = 1.0
END WHERE
!
DO jj=1,inlvls
   DO ji=1,ini
!
      psnowdz(ji,jj)   = zsnowfall_delta(ji)*(zsnowfall(ji) /inlvls) + &
                        (1.0-zsnowfall_delta(ji))*psnowdz(ji,jj)  
!
      psnowheat(ji,jj) = zsnowfall_delta(ji)*(psnowhmass(ji)/inlvls) + &
                       (1.0-zsnowfall_delta(ji))*psnowheat(ji,jj)  
!
      psnowrho(ji,jj)  = zsnowfall_delta(ji)*zrhosnew(ji)            + &
                       (1.0-zsnowfall_delta(ji))*psnowrho(ji,jj)
!
      psnowage(ji,jj)  = zsnowfall_delta(ji)*(zagenew(ji)/inlvls)    + &
                       (1.0-zsnowfall_delta(ji))*psnowage(ji,jj)  
!
   ENDDO
ENDDO
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LFALL',1,zhook_handle)
!
!
END SUBROUTINE snow3lfall
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snow3lcompactn(PTSTEP,PSNOWDZMIN,PSNOWRHO,PSNOWDZ,PSNOWTEMP,PSNOW,PSNOWLIQ)  
!
!!    PURPOSE
!!    -------
!     Snow compaction due to overburden and settling.
!     Mass is unchanged: layer thickness is reduced
!     in proportion to density increases. Method
!     of Brun et al (1989) and Vionnet et al. (2012)
!
!     
USE modd_surf_par, ONLY : xundef
!
USE modd_csts,     ONLY : xtt, xg
USE modd_snow_par, ONLY : xrhosmax_es
!
USE modd_snow_metamo, ONLY : xvvisc1,xvvisc3,xvvisc4, &
                             xvvisc5,xvvisc6,xvro11
!
USE yomhook   ,    ONLY : lhook,   dr_hook
USE parkind1  ,    ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: PTSTEP
REAL, INTENT(IN)                    :: PSNOWDZMIN
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PSNOWTEMP, PSNOWLIQ
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ
!
REAL, DIMENSION(:), INTENT(OUT)     :: PSNOW
!
!
!*      0.2    declarations of local variables
!
INTEGER                             :: JJ, JI
!
INTEGER                             :: INI
INTEGER                             :: INLVLS
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHO2, ZVISCOCITY, ZF1, &
                                                      ZTEMP, ZSMASS, ZSNOWDZ,     &
                                                      ZWSNOWDZ, ZWHOLDMAX
!
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
! 0. Initialization:
! ------------------
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LCOMPACTN',0,zhook_handle)
!
ini             = SIZE(psnowdz(:,:),1)
inlvls          = SIZE(psnowdz(:,:),2)
!
zsnowrho2(:,:) = psnowrho(:,:)
zsnowdz(:,:) = max(psnowdzmin,psnowdz(:,:))
zviscocity(:,:) = 0.0
ztemp(:,:) = 0.0
!
! 1. Cumulative snow mass (kg/m2):
! --------------------------------
!
zsmass(:,:) = 0.0
DO jj=2,inlvls
   DO ji=1,ini
      zsmass(ji,jj) = zsmass(ji,jj-1) + psnowdz(ji,jj-1)*psnowrho(ji,jj-1)
   ENDDO
ENDDO
! overburden of half the mass of the uppermost layer applied to itself
zsmass(:,1) = 0.5 * psnowdz(:,1) * psnowrho(:,1)
!
! 2. Compaction
! -------------
!
!Liquid water effect
!
zwholdmax(:,:) = snow3lhold(psnowrho,psnowdz)
zwholdmax(:,:) = max(1.e-10, zwholdmax(:,:))
zf1(:,:) = 1.0/(xvvisc5+10.*min(1.0,psnowliq(:,:)/zwholdmax(:,:)))
!
!Snow viscocity, density and grid thicknesses
!
DO jj=1,inlvls
   DO ji=1,ini
!   
      IF(psnowrho(ji,jj) < xrhosmax_es)THEN
!   
!       temperature dependence limited to 5K: Schleef et al. (2014)
        ztemp(ji,jj) = xvvisc4*min(5.0,abs(xtt-psnowtemp(ji,jj)))
!        
!       Calculate snow viscocity: Brun et al. (1989), Vionnet et al. (2012)
        zviscocity(ji,jj) = xvvisc1*zf1(ji,jj)*exp(xvvisc3*psnowrho(ji,jj)+ztemp(ji,jj))*psnowrho(ji,jj)/xvro11
!
!       Calculate snow density:
        zsnowrho2(ji,jj) = psnowrho(ji,jj) + psnowrho(ji,jj)*ptstep &
                         * ( (xg*zsmass(ji,jj)/zviscocity(ji,jj)) )
!         
!       Conserve mass by decreasing grid thicknesses in response to density increases
        psnowdz(ji,jj) = psnowdz(ji,jj)*(psnowrho(ji,jj)/zsnowrho2(ji,jj))  
!        
      ENDIF
!
   ENDDO
ENDDO
!
! 3. Update total snow depth and density profile:
! -----------------------------------------------
!
! Compaction/augmentation of total snowpack depth
!
psnow(:) = 0.
DO jj=1,inlvls
   DO ji=1,ini
      psnow(ji) = psnow(ji) + psnowdz(ji,jj)
   ENDDO
ENDDO
!
! Update density (kg m-3):
!
psnowrho(:,:)  = zsnowrho2(:,:)
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LCOMPACTN',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE snow3lcompactn
!####################################################################
!####################################################################
!####################################################################
        SUBROUTINE snow3ltransf(PSNOW,PSNOWDZ,PSNOWDZN,    &
                                PSNOWRHO,PSNOWHEAT,PSNOWAGE)  
!
!!    PURPOSE
!!    -------
!     Snow mass,heat and characteristics redistibution in case of
!     grid resizing. Total mass and heat content of the overall snowpack
!     unchanged/conserved within this routine.
!     Same method as in Crocus
!
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(INOUT) :: PSNOWDZN  
REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHEAT
REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO
REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ
REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWAGE
!
!*      0.2    declarations of local variables
!
INTEGER                             :: JI, JL, JLO
!
INTEGER                             :: INI
INTEGER                             :: INLVLS
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHON
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEATN
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWAGEN
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZTOP_NEW
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZBOT_NEW                                       
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHOO
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEATO
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWAGEO
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWDZO
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZTOP_OLD
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWZBOT_OLD  
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEAN
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWAGN
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZMASTOTN
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZMASSDZO
!
REAL, DIMENSION(SIZE(PSNOW)) :: ZPSNOW_OLD, ZPSNOW_NEW
REAL, DIMENSION(SIZE(PSNOW)) :: ZSUMHEAT, ZSUMSWE, ZSUMAGE, ZSNOWMIX_DELTA
!
REAL :: ZPROPOR
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
! 0. Initialization:
! ------------------
!
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LTRANSF',0,zhook_handle)
!
ini        = SIZE(psnowrho,1)
inlvls     = SIZE(psnowrho,2)
!
zpsnow_new(:) = 0.0
zpsnow_old(:) = psnow(:)
!
DO jl=1,inlvls
   DO ji=1,ini
      zpsnow_new(ji)=zpsnow_new(ji)+psnowdzn(ji,jl)
   ENDDO
ENDDO
!
! initialization of variables describing the initial snowpack 
!
zsnowdzo(:,:) = psnowdz(:,:)
zsnowrhoo(:,:) = psnowrho(:,:)
zsnowheato(:,:) = psnowheat(:,:)
zsnowageo(:,:) = psnowage(:,:)
zmassdzo(:,:) = xundef
!
! 1. Calculate vertical grid limits (m):
! --------------------------------------
!
zsnowztop_old(:,1) = zpsnow_old(:)
zsnowztop_new(:,1) = zpsnow_new(:)
zsnowzbot_old(:,1) = zsnowztop_old(:,1)-zsnowdzo(:,1)
zsnowzbot_new(:,1) = zsnowztop_new(:,1)-psnowdzn(:,1)
!
DO jl=2,inlvls
   DO ji=1,ini
      zsnowztop_old(ji,jl) = zsnowzbot_old(ji,jl-1)
      zsnowztop_new(ji,jl) = zsnowzbot_new(ji,jl-1)
      zsnowzbot_old(ji,jl) = zsnowztop_old(ji,jl  )-zsnowdzo(ji,jl)
      zsnowzbot_new(ji,jl) = zsnowztop_new(ji,jl  )-psnowdzn(ji,jl)
   ENDDO
ENDDO
zsnowzbot_old(:,inlvls)=0.0
zsnowzbot_new(:,inlvls)=0.0
!
! 3. Calculate mass, heat, charcateristics mixing due to vertical grid resizing:
! --------------------------------------------------------------------
!
! loop over the new snow layers
! Summ or avergage of the constituting quantities of the old snow layers
! which are totally or partially inserted in the new snow layer
! For snow age, mass weighted average is used.
!
zsnowhean(:,:)=0.0
zmastotn(:,:)=0.0
zsnowagn(:,:)=0.0
!
DO jl=1,inlvls
   DO jlo=1, inlvls   
      DO ji=1,ini
        IF((zsnowztop_old(ji,jlo)>zsnowzbot_new(ji,jl)).AND.(zsnowzbot_old(ji,jlo)<zsnowztop_new(ji,jl)))THEN
!                
          zpropor = (min(zsnowztop_old(ji,jlo), zsnowztop_new(ji,jl)) &
                  -  max(zsnowzbot_old(ji,jlo), zsnowzbot_new(ji,jl)))&
                  / zsnowdzo(ji,jlo) 
!
          zmassdzo(ji,jlo)=zsnowrhoo(ji,jlo)*zsnowdzo(ji,jlo)*zpropor
!
          zmastotn(ji,jl)=zmastotn(ji,jl)+zmassdzo(ji,jlo)
          zsnowagn(ji,jl)=zsnowagn(ji,jl)+zsnowageo(ji,jlo)*zmassdzo(ji,jlo)
!
          zsnowhean(ji,jl)=zsnowhean(ji,jl)+zsnowheato(ji,jlo)*zpropor
!          
        ENDIF
      ENDDO 
    ENDDO 
ENDDO  
!
! the new layer inherits from the weighted average properties of the old ones
! heat and mass
!
zsnowheatn(:,:)= zsnowhean(:,:)
zsnowagen(:,:)= zsnowagn(:,:)/zmastotn(:,:)
zsnowrhon(:,:)= zmastotn(:,:)/psnowdzn(:,:)
!
!
! 4. 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):
!
zsumheat(:)       = 0.0
zsumswe(:)        = 0.0
zsumage(:)        = 0.0
zsnowmix_delta(:) = 0.0
!
DO jl=1,inlvls
   DO ji=1,ini
      IF(psnow(ji) < xsnowcritd)THEN
         zsumheat(ji) = zsumheat(ji) + psnowheat(ji,jl)
         zsumswe(ji) = zsumswe(ji) + psnowrho(ji,jl)*psnowdz(ji,jl)
         zsumage(ji) = zsumage(ji) + psnowage(ji,jl)
         zsnowmix_delta(ji) = 1.0
      ENDIF
   ENDDO
ENDDO
!
! Heat and mass are evenly distributed vertically:
! heat and mass (density and thickness) are constant
! in profile:
!
DO jl=1,inlvls
   DO ji=1,ini
!         
      zsnowheatn(ji,jl) = zsnowmix_delta(ji)*(zsumheat(ji)/inlvls)  + &
                         (1.0-zsnowmix_delta(ji))*zsnowheatn(ji,jl)  
!
      psnowdzn(ji,jl)   = zsnowmix_delta(ji)*(psnow(ji)/inlvls)     + &
                         (1.0-zsnowmix_delta(ji))*psnowdzn(ji,jl)  
!
      zsnowrhon(ji,jl)  = zsnowmix_delta(ji)*(zsumswe(ji)/psnow(ji)) + &
                         (1.0-zsnowmix_delta(ji))*zsnowrhon(ji,jl)
!
      zsnowagen(ji,jl)  = zsnowmix_delta(ji)*(zsumage(ji)/inlvls)  + &
                         (1.0-zsnowmix_delta(ji))*zsnowagen(ji,jl)  
!
   ENDDO
ENDDO
!
! 5. Update mass (density and thickness) and heat:
! ------------------------------------------------
!
psnowdz(:,:) = psnowdzn(:,:)
psnowrho(:,:) = zsnowrhon(:,:)
psnowheat(:,:) = zsnowheatn(:,:)
psnowage(:,:) = zsnowagen(:,:)
!
IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LTRANSF',1,zhook_handle)
!
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE snow3ltransf
!####################################################################
!####################################################################
!####################################################################


END MODULE mode_snow3l