mode_meb.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_meb
!     ################
!
!!****  *MODE_MEB * - contains multi-energy balance characteristics functions
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!    
!!    EXTERNAL
!!    --------
!!       
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!      P. Samuelsson            * SMHI *
!!      A. Boone                 * Meteo France *       
!!      S. Gollvik               * SMHI *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original        18/01/11
!-----------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!
!
INTERFACE mebpalphan
  MODULE PROCEDURE mebpalphan_3d
  MODULE PROCEDURE mebpalphan_2d
  MODULE PROCEDURE mebpalphan_1d
  MODULE PROCEDURE mebpalphan_0d
END INTERFACE
!
INTERFACE sfc_heatcap_veg
   MODULE PROCEDURE sfc_heatcap_veg_2d
   MODULE PROCEDURE sfc_heatcap_veg_1d
   MODULE PROCEDURE sfc_heatcap_veg_0d
END INTERFACE
!
INTERFACE swdown_diff
   MODULE PROCEDURE swdown_diff_2d
   MODULE PROCEDURE swdown_diff_1d
   MODULE PROCEDURE swdown_diff_0d
END INTERFACE
!
INTERFACE snow_intercept_eff
   MODULE PROCEDURE snow_intercept_eff_2d
   MODULE PROCEDURE snow_intercept_eff_1d
   MODULE PROCEDURE snow_intercept_eff_0d
END INTERFACE
!
INTERFACE meb_shield_factor
   MODULE PROCEDURE meb_shield_factor_2d
   MODULE PROCEDURE meb_shield_factor_1d
   MODULE PROCEDURE meb_shield_factor_0d
END INTERFACE
!
!-------------------------------------------------------------------------------
 CONTAINS
!
!####################################################################
!####################################################################
!####################################################################
      FUNCTION mebpalphan_3d(PSNOWDEPTH,PH_VEG) RESULT(PPALPHAN)
!
!!    PURPOSE
!!    -------
!     Calculation of p_alphan, the snow/canopy transition coefficient
!     0 for snow at canopy base
!     1 for snow at canopy top
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:,:), INTENT(IN)     :: psnowdepth,ph_veg
!
REAL, DIMENSION(SIZE(PSNOWDEPTH,1),SIZE(PSNOWDEPTH,2),SIZE(PSNOWDEPTH,3)) :: ppalphan
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL, DIMENSION(SIZE(PSNOWDEPTH,1),SIZE(PSNOWDEPTH,2),SIZE(PSNOWDEPTH,3))   :: zh_baseveg  ! height of the base of the canopy
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_3D',0,zhook_handle)
!
! multi-energy balance: Maybe the calculation of ZH_BASEVEG should be modified.
!WHERE(PH_VEG>1.)
!  ZH_BASEVEG(:,:,:)=0.3*PH_VEG(:,:,:)
!ELSEWHERE
!  ZH_BASEVEG(:,:,:)=0.
!END WHERE
zh_baseveg(:,:,:)=max(0.2*(ph_veg(:,:,:)-2.0),0.0);
!
ppalphan(:,:,:)=min(1.,max(0., (psnowdepth(:,:,:)-zh_baseveg(:,:,:))/(ph_veg(:,:,:)-zh_baseveg(:,:,:)) ))
!
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_3D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION mebpalphan_3d
!####################################################################
!####################################################################
!####################################################################
      FUNCTION mebpalphan_2d(PSNOWDEPTH,PH_VEG) RESULT(PPALPHAN)
!
!!    PURPOSE
!!    -------
!     Calculation of p_alphan, the snow/canopy transition coefficient
!     0 for snow at canopy base
!     1 for snow at canopy top
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)     :: psnowdepth,ph_veg
!
REAL, DIMENSION(SIZE(PSNOWDEPTH,1),SIZE(PSNOWDEPTH,2)) :: ppalphan
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL, DIMENSION(SIZE(PSNOWDEPTH,1),SIZE(PSNOWDEPTH,2))   :: zh_baseveg  ! height of the base of the canopy
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_2D',0,zhook_handle)
!
!WHERE(PH_VEG>1.)
!  ZH_BASEVEG(:,:)=0.3*PH_VEG(:,:)
!ELSEWHERE
!  ZH_BASEVEG(:,:)=0.
!END WHERE
zh_baseveg(:,:)=max(0.2*(ph_veg(:,:)-2.0),0.0);
!
ppalphan(:,:)=min(1.,max(0., (psnowdepth(:,:)-zh_baseveg(:,:))/(ph_veg(:,:)-zh_baseveg(:,:)) ))
!
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_2D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION mebpalphan_2d
!####################################################################
!####################################################################
!####################################################################
      FUNCTION mebpalphan_1d(PSNOWDEPTH,PH_VEG) RESULT(PPALPHAN)
!
!!    PURPOSE
!!    -------
!     Calculation of p_alphan, the snow/canopy transition coefficient
!     0 for snow at canopy base
!     1 for snow at canopy top
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)     :: psnowdepth,ph_veg
!
REAL, DIMENSION(SIZE(PSNOWDEPTH,1)) :: ppalphan
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL, DIMENSION(SIZE(PSNOWDEPTH,1))   :: zh_baseveg  ! height of the base of the canopy
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_1D',0,zhook_handle)
!
!WHERE(PH_VEG>1.)
!  ZH_BASEVEG(:)=0.3*PH_VEG(:)
!ELSEWHERE
!  ZH_BASEVEG(:)=0.
!END WHERE
zh_baseveg(:)=max(0.2*(ph_veg(:)-2.0),0.0);
!
!
ppalphan(:)=min(1.,max(0., (psnowdepth(:)-zh_baseveg(:))/(ph_veg(:)-zh_baseveg(:)) ))
!
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_1D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION mebpalphan_1d
!####################################################################
!####################################################################
!####################################################################
      FUNCTION mebpalphan_0d(PSNOWDEPTH,PH_VEG) RESULT(PPALPHAN)
!
!!    PURPOSE
!!    -------
!     Calculation of p_alphan, the snow/canopy transition coefficient
!     0 for snow at canopy base
!     1 for snow at canopy top
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)     :: psnowdepth,ph_veg
!
REAL :: ppalphan
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL :: zh_baseveg  ! height of the base of the canopy
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_0D',0,zhook_handle)
!
!IF(PH_VEG>1.)THEN
!  ZH_BASEVEG=0.3*PH_VEG
!ELSE
!  ZH_BASEVEG=0.
!ENDIF
zh_baseveg=max(0.2*(ph_veg-2.0),0.0);
!
ppalphan=min(1.,max(0., (psnowdepth-zh_baseveg)/(ph_veg-zh_baseveg) ))
!
IF (lhook) CALL dr_hook('MODE_MEB:MEBPALPHAN_0D',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION mebpalphan_0d
!####################################################################
!####################################################################
!####################################################################
FUNCTION sfc_heatcap_veg_0d(PWRN,PWR,PCV) RESULT(ZCHEATV)

! Compute the bulk heat capacity of the vegetation canopy

USE modd_csts,     ONLY : xcl, xci
USE modd_isba_par, ONLY : xcvheatf

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!*      0.1    declarations of arguments
!
REAL,               INTENT(IN)  :: pwrn, pwr, pcv
!                                  PWRN      = Liquid water equivalent mass of intercepted snow (kg m-2)
!                                  PCV       = Thermal inertia of the vegetation (m2 K J-1)
!                                  PWR       = Liquid water mass intercepted (kg m-2)
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL                            :: zcheatv
!                                  ZCHEATV = Total bulk Vegetation canopy heat capacity (J m-2 K-1)
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER                 :: zcheatvmin = 1.e+4 ! minimum limit (J m-2 K-1)
!
!------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SFC_HEATCAP_VEG_0D',0,zhook_handle)
!
! Total bulk canopy heat capacity
! Method: we use the ratio of biomass to LAI to get a total biomass estimate,
! next assume that the dry biomass heat capacity is small compared to that
! of the water contained in the vegetation to arrive at the vegetation part, then
! we add the heat capacities of intercepted liquid and frozen water.
! Finally, use a minimum value to avoid numerical jumps, while still ensuring that
! the heat capacity of the vegetation is generally < a typical daily restore for the soil

zcheatv   = max(zcheatvmin,xcvheatf/pcv)           +   & ! stems, branches, trunk...
               xci*pwrn                            +   & ! intercepted snow
               xcl*pwr                                   ! intercepted water

IF (lhook) CALL dr_hook('MODE_MEB:SFC_HEATCAP_VEG_0D',1,zhook_handle)

END FUNCTION sfc_heatcap_veg_0d
!####################################################################
!####################################################################
!####################################################################
FUNCTION sfc_heatcap_veg_1d(PWRN,PWR,PCV) RESULT(ZCHEATV)

! Compute the bulk heat capacity of the vegetation canopy

USE modd_csts,     ONLY : xcl, xci
USE modd_isba_par, ONLY : xcvheatf

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)  :: pwrn, pwr, pcv
!                                  PWRN      = Liquid water equivalent mass of intercepted snow (kg m-2)
!                                  PCV       = Thermal inertia of the vegetation (m2 K J-1)
!                                  PWR       = Liquid water mass intercepted (kg m-2)
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL, DIMENSION(SIZE(PCV))     :: zcheatv
!                                  ZCHEATV = Total bulk Vegetation canopy heat capacity (J m-2 K-1)
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER                 :: zcheatvmin = 1.e+4 ! minimum limit (J m-2 K-1)
!
!------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SFC_HEATCAP_VEG_1D',0,zhook_handle)
!
! Total bulk canopy heat capacity
! Method: we use the ratio of biomass to LAI to get a total biomass estimate,
! next assume that the dry biomass heat capacity is small compared to that
! of the water contained in the vegetation to arrive at the vegetation part, then
! we add the heat capacities of intercepted liquid and frozen water.
! Finally, use a minimum value to avoid numerical jumps, while still ensuring that
! the heat capacity of the vegetation is generally < a typical daily restore for the soil

zcheatv(:)   = max(zcheatvmin,xcvheatf/pcv(:))        +   & ! stems, branches, trunk...
               xci*pwrn(:)                            +   & ! intercepted snow
               xcl*pwr(:)                                   ! intercepted water

IF (lhook) CALL dr_hook('MODE_MEB:SFC_HEATCAP_VEG_1D',1,zhook_handle)

END FUNCTION sfc_heatcap_veg_1d
!####################################################################
!####################################################################
!####################################################################
FUNCTION sfc_heatcap_veg_2d(PWRN,PWR,PCV) RESULT(ZCHEATV)

! Compute the bulk heat capacity of the vegetation canopy

USE modd_csts,     ONLY : xcl, xci
USE modd_isba_par, ONLY : xcvheatf

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN) :: pwrn, pwr, pcv
!                                   PWRN      = Liquid water equivalent mass of intercepted snow (kg m-2)
!                                   PCV       = Thermal inertia of the vegetation (m2 K J-1)
!                                   PWR       = Liquid water mass intercepted (kg m-2)
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL, DIMENSION(SIZE(PCV),SIZE(PCV,2)) :: zcheatv
!                                  ZCHEATV = Total bulk Vegetation canopy heat capacity (J m-2 K-1)
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER                 :: zcheatvmin = 1.e+4 ! minimum limit (J m-2 K-1)
!
!------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SFC_HEATCAP_VEG_2D',0,zhook_handle)
!
! Total bulk canopy heat capacity
! Method: we use the ratio of biomass to LAI to get a total biomass estimate,
! next assume that the dry biomass heat capacity is small compared to that
! of the water contained in the vegetation to arrive at the vegetation part, then
! we add the heat capacities of intercepted liquid and frozen water.
! Finally, use a minimum value to avoid numerical jumps, while still ensuring that
! the heat capacity of the vegetation is generally < a typical daily restore for the soil

zcheatv(:,:) = max(zcheatvmin,xcvheatf/pcv(:,:))      +   & ! stems, branches, trunk...
               xci*pwrn(:,:)                          +   & ! intercepted snow
               xcl*pwr(:,:)                                 ! intercepted water

IF (lhook) CALL dr_hook('MODE_MEB:SFC_HEATCAP_VEG_2D',1,zhook_handle)

END FUNCTION sfc_heatcap_veg_2d
!####################################################################
!####################################################################
!####################################################################
!####################################################################
!####################################################################
!####################################################################
FUNCTION swdown_diff_2d(PSW_RAD,PCOSZENITH) RESULT(ZSWDOWN_DIFF)

! Based on incoming total rad (direct and diffuse) and cosine of the
! solar zenith angle, compute the fraction of that rad which is diffuse.
! D.G. Erbs, S.A. Klein and J.A. Duffie, 1982:
! Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation.
! Solar Energy, 28(4), 293-302.
!
! Author: A. Boone CNRM-GAME

USE modd_csts, ONLY : xi0

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!*      0.1    declarations of arguments
!
REAL,    DIMENSION(:,:), INTENT(IN)                 :: psw_rad 
                                                             
REAL,    DIMENSION(:,:), INTENT(IN)                 :: pcoszenith 

!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL,    DIMENSION(SIZE(PSW_RAD,1),SIZE(PSW_RAD,2)) :: zswdown_diff, zratio

!*      0.3    declarations of local parameters

REAL, PARAMETER                                     :: zswcnt = 1.0 ! 1.000 if all wavelength SWdown 
                                                                    ! 2.083 (=1/0.48) if just visible
!--------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SWDOWN_DIFF_2D',0,zhook_handle)

zratio(:,:)       = psw_rad(:,:)*(zswcnt/xi0)/max(0.01,pcoszenith(:,:))

! Numerical Check:

zratio(:,:)       = min(1.0, zratio(:,:))                                 

zswdown_diff(:,:) = 0.165 ! RATIO >= 0.8

WHERE(zratio(:,:) <  0.22 )                                        &
     zswdown_diff(:,:) = 1.0 - 0.09*zratio(:,:)

WHERE(zratio(:,:) >= 0.22 .AND. zratio(:,:) < 0.80)                &
     zswdown_diff(:,:) = 0.9511 + (-0.1604 + (4.388 + (-16.64 +    &
     12.34*zratio(:,:))*zratio(:,:))*zratio(:,:))*zratio(:,:)

IF (lhook) CALL dr_hook('MODE_MEB:SWDOWN_DIFF_2D',1,zhook_handle)

END FUNCTION swdown_diff_2d
!####################################################################
!####################################################################
!####################################################################
FUNCTION swdown_diff_1d(PSW_RAD,PCOSZENITH) RESULT(ZSWDOWN_DIFF)

! Based on incoming total rad (direct and diffuse) and cosine of the
! solar zenith angle, compute the fraction of that rad which is diffuse.
! D.G. Erbs, S.A. Klein and J.A. Duffie, 1982:
! Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation.
! Solar Energy, 28(4), 293-302.
!
! Author: A. Boone CNRM-GAME

USE modd_csts, ONLY : xi0

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!*      0.1    declarations of arguments
!
REAL,    DIMENSION(:), INTENT(IN)                 :: psw_rad 
                                                             
REAL,    DIMENSION(:), INTENT(IN)                 :: pcoszenith 

!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL,    DIMENSION(SIZE(PSW_RAD))                 :: zswdown_diff, zratio

!*      0.3    declarations of local parameters

REAL, PARAMETER                                   :: zswcnt = 1.0 ! 1.000 if all wavelength SWdown 
                                                                  ! 2.083 (=1/0.48) if just visible
!
!--------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SWDOWN_DIFF_1D',0,zhook_handle)

zratio(:) = psw_rad(:)*(zswcnt/xi0)/max(0.01,pcoszenith(:))

! Numerical Check:

zratio(:) = min(1.0, zratio(:))                                 

zswdown_diff(:) = 0.165 ! RATIO >= 0.8

WHERE(zratio(:) <  0.22 )                                      &
     zswdown_diff(:) = 1.0 - 0.09*zratio(:)

WHERE(zratio(:) >= 0.22 .AND. zratio(:) < 0.80)                &
     zswdown_diff(:) = 0.9511 + (-0.1604 + (4.388 + (-16.64 +  &
     12.34*zratio(:))*zratio(:))*zratio(:))*zratio(:)

IF (lhook) CALL dr_hook('MODE_MEB:SWDOWN_DIFF_1D',1,zhook_handle)

END FUNCTION swdown_diff_1d
!####################################################################
!####################################################################
!####################################################################
FUNCTION swdown_diff_0d(PSW_RAD,PCOSZENITH) RESULT(ZSWDOWN_DIFF)

! Based on incoming total rad (direct and diffuse) and cosine of the
! solar zenith angle, compute the fraction of that rad which is diffuse.
! D.G. Erbs, S.A. Klein and J.A. Duffie, 1982:
! Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation.
! Solar Energy, 28(4), 293-302.
!
! Author: A. Boone CNRM-GAME

USE modd_csts, ONLY : xi0

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                 :: psw_rad 
                                                             
REAL, INTENT(IN)                 :: pcoszenith 

!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
REAL                             :: zswdown_diff, zratio
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER                  :: zswcnt = 1.0 ! 1.000 if all wavelength SWdown 
                                                                  ! 2.083 (=1/0.48) if just visible
!--------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SWDOWN_DIFF_0D',0,zhook_handle)

zratio       = psw_rad*(zswcnt/xi0)/max(0.01,pcoszenith)

! Numerical Check:

zratio       = min(1.0, zratio)                                 

zswdown_diff = 0.165 ! RATIO >= 0.8

IF    (zratio <  0.22 )THEN
     zswdown_diff = 1.0 - 0.09*zratio
ELSEIF(zratio >= 0.22 .AND. zratio < 0.80)THEN
     zswdown_diff = 0.9511 + (-0.1604 + (4.388 + (-16.64 +  &
                    12.34*zratio)*zratio)*zratio)*zratio
ENDIF

IF (lhook) CALL dr_hook('MODE_MEB:SWDOWN_DIFF_0D',1,zhook_handle)

END FUNCTION swdown_diff_0d
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow_intercept_eff_0d(PCHIP,PVELC,PWRVNMAX) RESULT(ZKVN)
!!
!! Calculate snow interception efficiency.
!!
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      P. Samuelsson           * SMHI *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       02/2011
!!      (A. Boone)  08/08/2011 Transform from subroutine to function
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL,               INTENT(IN)   :: pchip, pvelc
!                                     PCHIP      = view factor (for LW) 
!                                     PVELC      = wind speed at top of vegetation
!
REAL,               INTENT(IN)   :: pwrvnmax
!                                     PWRVNMAX   = maximum equivalent snow content
!                                                  in the canopy vegetation
!
!*      0.2    declarations of local variables
!
!
REAL                             :: zkvn
!                                     ZKVN       = Snow interception efficiency
!                                                  coefficient.
!                                                  Note: if this is set=0 it means
!                                                  that snow interception is shut
!                                                  off.
!
REAL                             :: zfcp
!                                     ZFCP       = snow interception factor
REAL(KIND=JPRB) :: zhook_handle
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER             :: zwsnow        = 0.8        ! Snow fall velocity (m/s)
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SNOW_INTERCEPT_EFF_0D',0,zhook_handle)
!
! Initialization:
!
zfcp        = 0.0 
zkvn        = 0.0
!
! Snow interception efficiency
!
IF(pwrvnmax > 0.0)THEN
   zfcp     = min(1.,max(0., pvelc/((2*zwsnow)*pchip) ) )
   zkvn     = (1.-pchip+zfcp*pchip)/pwrvnmax
ENDIF
!
IF (lhook) CALL dr_hook('MODE_MEB:SNOW_INTERCEPT_EFF_0D',1,zhook_handle)
END FUNCTION snow_intercept_eff_0d
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow_intercept_eff_1d(PCHIP,PVELC,PWRVNMAX) RESULT(ZKVN)
!!
!! Calculate snow interception efficiency.
!!
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      P. Samuelsson           * SMHI *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       02/2011
!!      (A. Boone)  08/08/2011 Transform from subroutine to function
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)   :: pchip, pvelc
!                                     PCHIP      = view factor (for LW) 
!                                     PVELC      = wind speed at top of vegetation
!
REAL, DIMENSION(:), INTENT(IN)   :: pwrvnmax
!                                     PWRVNMAX   = maximum equivalent snow content
!                                                  in the canopy vegetation
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PVELC))     :: zkvn
!                                     ZKVN       = Snow interception efficiency
!                                                  coefficient.
!                                                  Note: if this is set=0 it means
!                                                  that snow interception is shut
!                                                  off.
!
REAL, DIMENSION(SIZE(PVELC))     :: zfcp
!                                     ZFCP       = snow interception factor
REAL(KIND=JPRB) :: zhook_handle
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER                  :: zwsnow       = 0.8        ! Snow fall velocity (m/s)
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SNOW_INTERCEPT_EFF_1D',0,zhook_handle)
!
! Initialization:
!
zfcp(:)        = 0.0 
zkvn(:)        = 0.0
!
! Snow interception efficiency
!
WHERE(pwrvnmax(:) > 0.0)
   zfcp(:)     = min(1.,max(0., pvelc(:)/((2*zwsnow)*pchip(:)) ) )
   zkvn(:)     = (1.-pchip(:)+zfcp(:)*pchip(:))/pwrvnmax(:)
END WHERE
!
IF (lhook) CALL dr_hook('MODE_MEB:SNOW_INTERCEPT_EFF_1D',1,zhook_handle)
END FUNCTION snow_intercept_eff_1d
!####################################################################
!####################################################################
!####################################################################
FUNCTION snow_intercept_eff_2d(PCHIP,PVELC,PWRVNMAX) RESULT(ZKVN)
!!
!! Calculate snow interception efficiency.
!!
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      P. Samuelsson           * SMHI *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       02/2011
!!      (A. Boone)  08/08/2011 Transform from subroutine to function
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)   :: pchip, pvelc
!                                     PCHIP      = view factor (for LW) 
!                                     PVELC      = wind speed at top of vegetation
!
REAL, DIMENSION(:,:), INTENT(IN)   :: pwrvnmax
!                                     PWRVNMAX   = maximum equivalent snow content
!                                                  in the canopy vegetation
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PVELC,1),SIZE(PVELC,2)) :: zkvn
!                                     ZKVN       = Snow interception efficiency
!                                                  coefficient.
!                                                  Note: if this is set=0 it means
!                                                  that snow interception is shut
!                                                  off.
!
REAL, DIMENSION(SIZE(PVELC,1),SIZE(PVELC,2)) :: zfcp
!                                     ZFCP       = snow interception factor
REAL(KIND=JPRB) :: zhook_handle
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER                  :: zwsnow       = 0.8        ! Snow fall velocity (m/s)
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:SNOW_INTERCEPT_EFF_2D',0,zhook_handle)
!
! Initialization:
!
zfcp(:,:)     = 0.0 
zkvn(:,:)     = 0.0
!
! Snow interception efficiency
!
WHERE(pwrvnmax(:,:) > 0.0)
   zfcp(:,:)  = min(1.,max(0., pvelc(:,:)/((2*zwsnow)*pchip(:,:)) ) )
   zkvn(:,:)  = (1.-pchip(:,:)+zfcp(:,:)*pchip(:,:))/pwrvnmax(:,:)
END WHERE
!
IF (lhook) CALL dr_hook('MODE_MEB:SNOW_INTERCEPT_EFF_2D',1,zhook_handle)
END FUNCTION snow_intercept_eff_2d
!####################################################################
!####################################################################
!####################################################################
FUNCTION meb_shield_factor_0d(PLAI,PPALPHAN) RESULT(PCHIP)
!!
!! Calculate MEB shielding factor
!!
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!  A. Boone                * Meteo France *       
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       02/2011
!!      P. Samuelsson  07/2014 Transform from subroutine to function
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_meb_par,   ONLY : xtau_lw
!
USE yomhook,        ONLY : lhook,   dr_hook
USE parkind1,       ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL,               INTENT(IN)   :: plai, ppalphan
!                                   PWRVNMAX   = canopy vegetation leaf area index
!                                   PPALPHAN   = snow/canopy transition coefficient
!
REAL                             :: pchip
!                                   PCHIP      = shielding factor (for LW) 
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:MEB_SHIELD_FACTOR_0D',0,zhook_handle)
!
pchip = exp(-xtau_lw*plai*(1.-ppalphan))
!
IF (lhook) CALL dr_hook('MODE_MEB:MEB_SHIELD_FACTOR_0D',1,zhook_handle)
END FUNCTION meb_shield_factor_0d
!####################################################################
!####################################################################
!####################################################################
FUNCTION meb_shield_factor_1d(PLAI,PPALPHAN) RESULT(PCHIP)
!!
!! Calculate MEB shielding factor
!!
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!  A. Boone                * Meteo France *       
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       02/2011
!!      P. Samuelsson  07/2014 Transform from subroutine to function
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_meb_par,   ONLY : xtau_lw
!
USE yomhook,        ONLY : lhook,   dr_hook
USE parkind1,       ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)   :: plai, ppalphan
!                                   PWRVNMAX   = canopy vegetation leaf area index
!                                   PPALPHAN   = snow/canopy transition coefficient
!
REAL, DIMENSION(SIZE(PLAI,1))    :: pchip
!                                   PCHIP      = shielding factor (for LW) 
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:MEB_SHIELD_FACTOR_1D',0,zhook_handle)
!
pchip(:) = exp(-xtau_lw*plai(:)*(1.-ppalphan(:)))
!
IF (lhook) CALL dr_hook('MODE_MEB:MEB_SHIELD_FACTOR_1D',1,zhook_handle)
END FUNCTION meb_shield_factor_1d
!####################################################################
!####################################################################
!####################################################################
FUNCTION meb_shield_factor_2d(PLAI,PPALPHAN) RESULT(PCHIP)
!!
!! Calculate MEB shielding factor
!!
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!  A. Boone                * Meteo France *       
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       02/2011
!!      P. Samuelsson  07/2014 Transform from subroutine to function
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_meb_par,   ONLY : xtau_lw
!
USE yomhook,        ONLY : lhook,   dr_hook
USE parkind1,       ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)   :: plai, ppalphan
!                                     PWRVNMAX   = canopy vegetation leaf area index
!                                     PPALPHAN   = snow/canopy transition coefficient
!
REAL, DIMENSION(SIZE(PLAI,1),SIZE(PLAI,2))    :: pchip
!                                                PCHIP      = shielding factor (for LW) 
!
!*      0.2    declarations of local variables
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_MEB:MEB_SHIELD_FACTOR_2D',0,zhook_handle)
!
pchip(:,:) = exp(-xtau_lw*plai(:,:)*(1.-ppalphan(:,:)))
!
IF (lhook) CALL dr_hook('MODE_MEB:MEB_SHIELD_FACTOR_2D',1,zhook_handle)
END FUNCTION meb_shield_factor_2d
!####################################################################
!####################################################################
!####################################################################
END MODULE mode_meb