drag_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.
!   ############################################################################
!
SUBROUTINE drag_meb(IO, PEK, DMK, DK, PTG, PTA, PQA, PVMOD, &
                    PWG, PWGI, PWSAT, PWFC, PEXNS, PEXNA, PPS,  PRR, PSR,   &
                    PRHOA, PZ0G_WITHOUT_SNOW, PZ0_MEBV, PZ0H_MEBV,          &
                    PZ0EFF_MEBV, PZ0_MEBN, PZ0H_MEBN, PZ0EFF_MEBN, PSNOWSWE,&
                    PCHIP, PTSTEP, PRS_VG, PRS_VN, PPALPHAN, PZREF, PUREF,  &
                    PDIRCOSZW, PSNCV, PDELTA,  PVELC,                       &
                    PRISNOW, PUSTAR2SNOW, PHUGI, PHVG, PHVN, PLEG_DELTA,    &
                    PLEGI_DELTA, PHSGL, PHSGF, PFLXC_C_A, PFLXC_N_A,        &
                    PFLXC_G_C, PFLXC_N_C, PFLXC_VG_C, PFLXC_VN_C, PFLXC_MOM,&
                    PQSATG, PQSATV, PQSATC, PQSATN, PDELTAVK           )
!                    
!   ############################################################################
!
!!****  *DRAG_MEB*  
!!
!!    PURPOSE
!!    -------
!!
!!     Calculates the coefficients for heat transfers
!!     for the multiple-energy balance, and Halstead coefficients.
!!     Also estimate patch average specific and relative humidity, as well
!!     as exchange coefficients and Richardson's number
!!         
!!     
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    USE MODD_CST 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      P.Samuelsson/S.Gollvik            * SMHI *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    06/2010
!!      For MEB     01/2011
!                   10/2014 (A. Boone) Remove understory compsite vegetation
!                   10/2014 (A. Napoly) Added ground/litter resistance
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_pe_t
USE modd_diag_n, ONLY : diag_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
!
USE modd_csts,     ONLY : xpi
USE modd_snow_par, ONLY : xz0sn
USE modd_isba_par, ONLY : xwgmin
USE modd_surf_atm, ONLY : ldrag_coef_arp, xrimax, lrrgust_arp, xrrscale,   &
                          xrrgamma, xutilgust, lcpl_arp
USE modd_meb_par,  ONLY : xkdelta_wr
!
USE modi_surface_cdch_1darp
USE modi_wind_threshold
USE modi_disph_for_meb
USE modi_preps_for_meb_drag
USE modi_surface_air_meb
!
USE mode_thermos, ONLY : qsat, qsati
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(diag_t), INTENT(INOUT) :: DK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
!
REAL, INTENT(IN)                 ::   PTSTEP 
!                                     PTSTEP = Model time step (s)        
!
REAL, DIMENSION(:), INTENT(IN)   :: PTG, PTA, PQA, PVMOD, PWG, PWGI, PWSAT, PWFC, &
                                    PEXNS, PEXNA, PPS, PSNOWSWE
!                                     PTG     = surface temperature
!                                     PSNOWSWE  = sfc layer snow water equiv (SWE) kg m-2
!                                     PTA     = temperature of the atmosphere/forcing-level
!                                     PQA = specific humidity of the atmosphere/forcing-level
!                                     PVMOD = module of the horizontal wind
!                                             NOTE it should be input as >= 1. (m)
!                                     PWG     = near-surface volumetric water content
!                                     PWGI    = near-surface frozen volumetric water content
!                                     PWSAT   = volumetric water content at saturation
!                                     PWFC = field capacity volumetric water content
!                                     PEXNS   = Exner function at the surface
!                                     PEXNA= Exner function of the atmosphere/forcing-level
!                                     PPS = surface pressure
!
REAL, DIMENSION(:), INTENT(IN)   ::  PRR, PSR, PRHOA
!                                     PRR = rain rate    
!                                     PSR = snow rate             
!                                     PRHOA = near-ground air density
!
REAL, DIMENSION(:), INTENT(IN)   :: PZ0G_WITHOUT_SNOW, &
                                    PZ0_MEBV, PZ0H_MEBV, PZ0EFF_MEBV, &
                                    PZ0_MEBN, PZ0H_MEBN, PZ0EFF_MEBN
!
!                    PZ0G_WITHOUT_SNOW  ! roughness length for momentum at snow-free canopy floor
!                    PZ0_MEBV           ! roughness length for momentum over MEB vegetation part of patch
!                    PZ0H_MEBV          ! roughness length for heat over MEB vegetataion part of path
!                    PZ0_MEBN           ! roughness length for momentum over MEB snow part of patch
!                    PZ0H_MEBN          ! roughness length for heat over MEB snow part of path
!
REAL, DIMENSION(:), INTENT(IN)   :: PCHIP
!                                     PCHIP = view factor (for LW) 
!
REAL, DIMENSION(:), INTENT(IN)   :: PRS_VG, PRS_VN, PSNCV,  &
                                    PPALPHAN, PZREF, PUREF, PDIRCOSZW
!
!                                     PRS_VG   = surface resistance for canopy (1-png)
!                                     PRS_VN   = surface resistance for canopy (png)
!                                     PSNCV   = fraction of the canopy vegetation covered by snow
!                                     PPALPHAN = weight between canopy air flow and direct flow 
!                                                 between snow and atmosphere
!                                     PZREF    = reference height of the first
!                                                atmospheric level 
!                                     PUREF    = reference height of the wind
!                                                NOTE this is different from ZZREF
!                                                ONLY in stand-alone/forced mode,
!                                                NOT when coupled to a model (MesoNH)
!                                     PDIRCOSZW= Cosinus of the angle between 
!                                                the normal to the surface and the vertical
!
REAL, DIMENSION(:), INTENT(IN)    :: PDELTA
!                                     PDELTA   = fraction of the canopy foliage covered
!                                                by intercepted water (-)
!
REAL, DIMENSION(:), INTENT(OUT)  :: PDELTAVK
!                                     PDELTAVK = fraction of the canopy foliage covered
!                                                by intercepted water *including* K-factor (-)
!                                                (i.e. that the intercepted water fraction is not 
!                                                 necessarily reaching 100%)
!
REAL, DIMENSION(:), INTENT(OUT)  :: PVELC
!                                     PVELC = wind speed at top of vegetation
!
REAL, DIMENSION(:), INTENT(OUT)  :: PHUGI, PHVG, PHVN, PLEG_DELTA, PLEGI_DELTA, PHSGL, PHSGF
!
!                                    PHUGI = ground (ice) relative humidity
!                                    PHVN  = Halstead coefficient vegetation canopy above snow
!                                    PHVG  = Halstead coefficient vegetation canopy above ground
!                                    PLEG_DELTA  = soil evaporation delta fn
!                                    PLEGI_DELTA = soil sublimation delta fn
!                                    PHSGL = surface halstead cofficient for bare soil (currently==1)
!                                    PHSGF = surface halstead cofficient for bare soil ice (currently==1)
!
!
REAL, DIMENSION(:), INTENT(OUT)  :: PFLXC_C_A, PFLXC_N_A, PFLXC_G_C, PFLXC_N_C, PFLXC_VG_C, PFLXC_VN_C
!  
!                                     EXCHANGE coefficients i.e. rho/resistance [kg/m2/s]
!  
!                                    PFLXC_C_A   between canopy air and atmosphere
!                                    PFLXC_N_A   between the snow on the ground and atmosphere
!                                    PFLXC_G_C   between snow-free ground and canopy air
!                                    PFLXC_N_C   between snow on the ground and canopy air
!                                    PFLXC_VG_C  between canopy over snow-free ground and canopy air
!                                    PFLXC_VN_C  between canopy over the snow on the ground and canopy air
!
REAL, DIMENSION(:), INTENT(OUT)  :: PFLXC_MOM, PQSATG, PQSATV, PQSATC, PQSATN
!                                    PFLXC_MOM = effective drag coefficient for momentum [kg/m2/s]
!                                    PQSATG = qsat for PTG
!                                    PQSATV = qsat for PEK%XTV(:)
!                                    PQSATC = qsat for PEK%XTC(:)
!                                    PQSATN = qsat for DMK%XSNOWTEMP(:)
!
REAL, DIMENSION(:), INTENT(OUT)  :: PRISNOW, PUSTAR2SNOW 
!                                    PRISNOW     = Richardson number over snow (-)
!                                    PUSTAR2SNOW = Surface friction velocity squared  (m2 s-2)
!                                                   Just a diagnostic, not used in coupling

!                                    
!*      0.2    declarations of local variables
!
!
!
REAL, DIMENSION(SIZE(PTG)) :: ZAC,ZWFC,ZWSAT,ZFP,ZRRCOR 
!                              ZQSATG = specific humidity at saturation at ground
!                              ZQSATN = specific humidity of snow
!                              ZAC    =  aerodynamical conductance
!                              ZWFC   = field capacity in presence of ice
!                              ZWSAT  = saturation in presence of ice
!                              ZFP    = working variable
!                              ZRRCOR = correction of CD, CH, CDN due to moist-gustiness
!
REAL, DIMENSION(SIZE(PTG)) :: ZCHIL, ZLAISN, ZLW, ZDISPH, ZVELC, ZRICN, ZRA_C_A, &
                              ZRA_G_C, ZG_VG_C, ZRA_N_C, ZG_VN_C 
!                              ZLAISN = leaf area index for snow covered ground
!                              ZLW = leaf width
!                              ZVELC = wind speed at top of vegetation
!                              ZRICN = Richardson nr between canopy air and atmosphere
!                              ZRA_C_A = resistance between canopy air and atmosphere
!                              ZRA_G_C = resistance between canopy air and snow free ground
!                              ZG_VG_C = conductance between canopy and canopy air (1-png)
!                              ZRA_N_C = aerodynamic resistance between the partly
!                                       snowcovered (png>0), ground and canopy air
!                              ZG_VN_C = conductance between canopy and canopy air (png)
!
REAL, DIMENSION(SIZE(PTG)) :: ZCHCN,ZCDNCN,ZCDCN !exchange coefficients between 
!                                                 canopy air and atmosphere
!
REAL, DIMENSION(SIZE(PTG)) :: ZRINN,ZRANN,ZCHNN,ZCDNNN,ZCDNN,ZTEFF,ZDELTAMAX,ZDELTAV,ZVMOD 
!
REAL, DIMENSION(SIZE(PTG)) :: ZRSGL,ZRSGF,ZZ0SN
!                             ZRSGL = surface resistance for bare soil (currently==0)
!                             ZRSGF = surface resistance for bare soil ice (currently==0)
!                             ZZ0SN = array of XZ0SN
!
REAL, DIMENSION(SIZE(PTG)) :: ZRSNFRAC, ZDENOM
!                             ZRSNFRAC = fraction to prevent/reduce sublimation of snow if too thin (-)
!                             ZDENOM   = working variable for denominator of an expression (*)
!
REAL, DIMENSION(SIZE(PTG)) :: ZUSTAR2G, ZCDG, ZCHG, ZRIG, ZPSNA 
!                                   ZUSTAR2G = canopy top friction velocity squared (m2/s2)
!                                   ZCDG     = drag coefficient (-)
!                                   ZCHG     = heat transfor coefficient (ground to canopy air) (-)
!                                   ZRIG     = Richardson number (ground to canopy air) (-)!                             
!                                   ZPSNA    = buried (by snow) canopy fraction (-)
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER            :: ZRAEPS       = 1.e-3  ! Safe limit of aerodynamic resistance to avoid dividing
!                                                   ! by zero when calculating exchange coefficients
REAL, PARAMETER            :: ZSNOWSWESMIN = 1.e-4  ! (kg m-2) Reduce sublimation from ground based snowpack as it
!                                                   !          becomes vanishingly thin 
REAL, PARAMETER            :: ZRG_COEF1    = 8.206  ! Ground/litter resistance coefficient 
REAL, PARAMETER            :: ZRG_COEF2    = 4.255  ! Ground/litter resistance coefficient 
!
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!               ---------------
!
IF (lhook) CALL dr_hook('DRAG_MEB',0,zhook_handle)
!
dk%XCH(:)   = 0.
dk%XCD(:)   = 0.
dk%XCDN(:)  = 0.
dk%XRI(:)   = 0.
dk%XHUG(:)  = 0.
!
dmk%XHV(:)   = 0.
dk%XHU(:)   = 0.
!
phugi(:) = 0.
phvn(:)  = 0.
phvg(:)  = 0.
phsgl(:) = 0.
phsgf(:) = 0.
!
pflxc_c_a(:)  = 0.
pflxc_g_c(:)  = 0.
pflxc_n_c(:)  = 0.
pflxc_vg_c(:)  = 0.
pflxc_vn_c(:) = 0.
pflxc_n_a(:)  = 0.
pflxc_mom(:) = 0.
!
zrsgl(:) = 0.
zrsgf(:) = 0.
zz0sn(:) = xz0sn
!
!-------------------------------------------------------------------------------
!
!*       1.     RELATIVE HUMIDITY OF THE GROUND HU
!               ----------------------------------
!
! this relative humidity is related to
! the superficial soil moisture and the
! field capacity of the ground
!
zwsat(:) = pwsat(:)-pwgi(:)
zwfc(:)  = pwfc(:)*zwsat(:)/pwsat(:)
!
dk%XHUG(:) = 0.5 * ( 1.-cos(xpi*min((pwg(:)-xwgmin) /zwfc(:),1.)) )
!
zwsat(:) = max(xwgmin, pwsat(:)-pwg(:))
zwfc(:) = pwfc(:)*zwsat(:)/pwsat(:)
phugi(:) = 0.5 * ( 1.-cos(xpi*min(pwgi(:)/zwfc(:),1.)) )
!
! there is a specific treatment for dew
! (see Mahfouf and Noilhan, jam, 1991)
!
pqsatg(:) = qsat(ptg(:),pps(:)) 
pqsatv(:) = qsat(pek%XTV(:),pps(:)) 
pqsatc(:) = qsat(pek%XTC(:),pps(:)) 
pqsatn(:) = qsati(dmk%XSNOWTEMP(:,1),pps(:)) 
!
!-------------------------------------------------------------------------------
!
!*       2.     GRID-AVERAGED HUMIDITY OF THE SURFACE
!               -------------------------------------
!       average of canopy air humidity and saturation humidity of snowpack
!
!
dk%XQS(:) =(1.-pek%XPSN(:)*ppalphan(:))*pek%XQC(:)          + pek%XPSN(:)*ppalphan(:)*pqsatn(:)
dk%XHU(:) =(1.-pek%XPSN(:)*ppalphan(:))*pek%XQC(:)/pqsatc(:)+ pek%XPSN(:)*ppalphan(:)
!
!-------------------------------------------------------------------------------
!
!*       3.     SPECIFIC HEAT OF THE SURFACE
!               ----------------------------
! Comment S.Gollvik 20110120
!     AARON IN E-BUDGET:
!
!IF (HCPSURF=='DRY') THEN
!    PCPS(:) = XCPD
!ELSEIF(.NOT.LCPL_ARP)THEN
!    PCPS(:) = XCPD + ( XCPV - XCPD ) * PQA(:)   
!ENDIF
!
!-------------------------------------------------------------------------------
!
!*       4.     COMPUTATION OF RESISTANCES/CONDUCTIVITES
!               ----------------------------------------------
!
zchil(:)=0.12
zlw(:)=0.02
!
! Calculate the displacement height
!
 CALL disph_for_meb(zchil,pek%XLAI,zlw,pek%XH_VEG,pzref,pz0_mebv,zdisph)
!
! Here ZRICN,ZRA_C_A,ZCHCN,ZCDNCN,ZCDCN are valid for the snow free part in meb:
!
 CALL preps_for_meb_drag(.true.,io%LFORC_MEASURE, pz0_mebv, pz0h_mebv, pz0eff_mebv,  &
                        pek%XH_VEG, pzref, pek%XTC, pta, pek%XQC, pqa, puref,   &
                        pvmod, pexna, pexns, pdircoszw, zdisph, pvelc, zvmod, zricn,&
                        zra_c_a, zchcn, zcdncn, zcdcn                    )
!
IF (lrrgust_arp) THEN

   zfp(:)=max(0.0,prr(:)+psr(:))
   zrrcor(:)=sqrt(1.0+((((zfp(:)/(zfp(:)+xrrscale))**xrrgamma)*xutilgust)**2) &
    & /(zcdcn(:)*zvmod(:)**2))

   zcdcn(:)  = zcdcn(:)  * zrrcor(:)
   zchcn(:)  = zchcn(:)  * zrrcor(:)
   zcdncn(:) = zcdncn(:) * zrrcor(:)

ENDIF
!       
pflxc_c_a(:)=zchcn(:)*zvmod(:)*prhoa(:)
!
!
! The momentum drag for this area:
!
pflxc_mom(:)=zcdcn(:)*zvmod(:)*prhoa(:)
!
!Calculate the aerodynamic resistance between the snowfree part of the ground
!and canopy air, ZRA_G_C, and the conductance between the canopy and canopy air, ZG_VG_C
!  
 CALL surface_air_meb(pz0_mebv, pz0h_mebv, pz0g_without_snow, pek%XH_VEG, pek%XLAI,   &
                     ptg, pek%XTC, pek%XTV, pvelc, zlw, zdisph,      &
                     zra_g_c, zg_vg_c, zustar2g, zcdg, zchg, zrig         )
!
!Compute the lai of the canopy that is above snow
!
zlaisn(:)= pek%XLAI(:)*(1.-ppalphan(:))
!
! Note that we use the same displacement height also for the snow covered part
!
! The same as ZRA_G_C/ZG_VG_C but for the snow part (png) =>ZRA_N_C/ZG_VN_C
!
 CALL surface_air_meb(pz0_mebn, pz0h_mebn, zz0sn, pek%XH_VEG, zlaisn,     &
                     dmk%XSNOWTEMP(:,1), pek%XTC, pek%XTV,  &
                     pvelc, zlw, zdisph, zra_n_c, zg_vn_c,           &
                     zustar2g, zcdg, zchg, zrig                      )

! save values over snow for diagnostic purposes:

pustar2snow(:) = zustar2g(:)
dmk%XCDSNOW(:)     = zcdg(:)  
dmk%XCHSNOW(:)     = zchg(:)
prisnow(:)     = zrig(:)
!
!------------------------------------------------------------------------------
! Now calculate the aerodynamic resistance for the completely snow covered part,
! i.e. between the snow surface and atmosphere directly
!
 CALL preps_for_meb_drag(.false.,io%LFORC_MEASURE,                        &
                   pz0_mebn, pz0h_mebn, pz0eff_mebn, pek%XH_VEG, pzref,  &
                   dmk%XSNOWTEMP(:,1), pta, pqsatn, pqa, puref, pvmod,   &
                   pexna, pexns, pdircoszw, zdisph,                      &
                   zvelc, zvmod, zrinn, zrann,                           &
                   zchnn,zcdnnn,zcdnn                                    )
!
IF (lrrgust_arp) THEN

   zrrcor(:)=sqrt(1.0+((((zfp(:)/(zfp(:)+xrrscale))**xrrgamma)*xutilgust)**2) &
          /(zcdnn(:)*zvmod(:)**2))

   zcdnn(:)  = zcdnn(:)  * zrrcor(:)
   zchnn(:)  = zchnn(:)  * zrrcor(:)
   zcdnnn(:) = zcdnnn(:) * zrrcor(:)

ENDIF
!
pflxc_n_a(:)=zchnn(:)*zvmod(:)*prhoa(:)
!
! The effective momentum drag: 
!
pflxc_mom(:)=(1.-pek%XPSN(:)*ppalphan(:))*pflxc_mom(:) +  &
                 pek%XPSN(:)*ppalphan(:)*zcdnn(:)*zvmod(:)*prhoa(:)
!
! Now calculate the aerodynamic resistance for  the whole meb-patch
!
! Calculate  the effective temperature
!
zpsna(:)   = pek%XPSN(:)*ppalphan(:)
!
zteff(:)   = (1.-zpsna(:))*pek%XTC(:)+ zpsna(:)*dmk%XSNOWTEMP(:,1)
!
! Some additional diagnostics:
!
pustar2snow(:) = (1.-zpsna(:))*pustar2snow(:) + zpsna(:)*zcdnn(:)*zvmod(:)**2
dmk%XCDSNOW(:) = (1.-zpsna(:))*dmk%XCDSNOW(:) + zpsna(:)*zcdnn(:)
dmk%XCHSNOW(:) = (1.-zpsna(:))*dmk%XCHSNOW(:) + zpsna(:)*zchnn(:)
prisnow(:)     = (1.-zpsna(:))*prisnow(:)     + zpsna(:)*zrinn(:)
!
!-------------------------------------------------------------------------------
!
 CALL preps_for_meb_drag(.false.,io%LFORC_MEASURE, dk%XZ0, dk%XZ0H, dk%XZ0EFF, &
                        pek%XH_VEG, pzref, zteff, pta, dk%XQS, pqa, puref, pvmod, &
                        pexna, pexns, pdircoszw, zdisph, zvelc, zvmod, dk%XRI,    &
                        pek%XRESA(:), dk%XCH,dk%XCDN,dk%XCD               )
!
IF (lrrgust_arp) THEN

   zrrcor(:)=sqrt(1.0+((((zfp(:)/(zfp(:)+xrrscale))**xrrgamma)*xutilgust)**2) &
                 /(dk%XCD(:)*zvmod(:)**2))

   dk%XCD(:)  = dk%XCD(:)  * zrrcor(:)
   dk%XCH(:)  = dk%XCH(:)  * zrrcor(:)
   dk%XCDN(:) = dk%XCDN(:) * zrrcor(:)

ENDIF

!-------------------------------------------------------------------------------
! Now start to define the internal exchange coefficients 
!-------------------------------------------------------------------------------
!
! PFLXC_G_C between snow-free ground and canopy air
! PFLXC_N_C  between snow on the ground and canopy air
! PFLXC_VG_C  between canopy over snow-free ground and canopy air
! PFLXC_VN_C  between canopy over the snow on the ground and canopy air
!
! ZRA_G_C = resistance between canopy air and snow free ground
! ZRA_N_C = aerodynamic resistance between the partly
!          snowcovered (png>0), ground and canopy air
! ZG_VG_C = conductance between canopy and canopy air (1-png)
! ZG_VN_C = conductance between canopy and canopy air (png)
!
! The resistances between ground/snow and canopy air go to zero when lai => 0
!
zra_g_c(:)   = max(zraeps,zra_g_c(:))
pflxc_g_c(:) = prhoa(:)  / zra_g_c(:)
phsgl(:)     = zra_g_c(:)/(zra_g_c(:)+zrsgl(:))
phsgf(:)     = zra_g_c(:)/(zra_g_c(:)+zrsgf(:))

zra_n_c(:)   = max(zra_n_c(:),zraeps)
pflxc_n_c(:) = prhoa(:)/zra_n_c(:)
!
! Reduce sublimation from ground based snowpack as it
! becomes vanishingly thin. This is mainly just for numerical reasons within the snow scheme
! and has little to
! no impact on actual grid box average fluxes since sublimation is multiplied by the snow fraction,
! which in this case is quite small. 
!
pflxc_n_c(:)  = pflxc_n_c(:)*min(1., (psnowswe(:) + psr(:)*ptstep)/zsnowswesmin)
!
! unit conversion:
!
pflxc_vg_c(:) = prhoa(:)*zg_vg_c(:)
pflxc_vn_c(:) = prhoa(:)*zg_vn_c(:)
!
!-------------------------------------------------------------------------------
!
!*       5.     HALSTEAD COEFFICIENT (RELATIVE HUMIDITY OF THE VEGETATION)
!               ----------------------------------------------------------
!
zdeltamax(:) = (1.-pchip(:))*(1.-pek%XPSN(:)*ppalphan(:))*prr(:)+ pek%XWR(:)/ptstep
zdenom(:)    = (1.-psncv(:))*xkdelta_wr *                                               &
               ( pek%XPSN(:)*(1.-ppalphan(:))*pflxc_vn_c(:) + (1.-pek%XPSN(:))*pflxc_vg_c(:) )* &
               ( pqsatv(:)-pek%XQC(:))
zdeltamax(:) = max(0., min(1.0, zdeltamax(:)/max(1.e-10,zdenom(:))))
!
zdeltav(:)   = xkdelta_wr*min(zdeltamax(:),pdelta(:))
pdeltavk(:)  = xkdelta_wr*pdelta(:)                    ! delta including K factor
!
phvg(:) = 1. - max(0.,sign(1.,pqsatv(:)-pek%XQC(:)))             &
          *(1.-zdeltav(:))*prs_vg(:)*zg_vg_c(:) / (1.+prs_vg(:)*zg_vg_c(:))
!
phvn(:) = 1. - max(0.,sign(1.,pqsatv(:)-pek%XQC(:)))             &
          *(1.-zdeltav(:))*prs_vn(:)*zg_vn_c(:) / (1.+prs_vn(:)*zg_vn_c(:))

! Diagnostics:
! Compute an effective canopy stomatal resistance (s m-1) and Halstead Coef: 
!
dmk%XRS(:)  = ppalphan(:)*prs_vn(:) + (1.0-ppalphan(:))*prs_vg(:)
dmk%XHV(:)  = ppalphan(:)*phvn(:)   + (1.0-ppalphan(:))*phvg(:)
!
!
!-------------------------------------------------------------------------
!*       6.    LITTER/GROUND RESISTANCE
!              ------------------------
!
! Inclusion of a ground resistance in the computation of the ground evaporation. 
! We use the existing LEG_DELTA (formerly a delta function) as a Beta-type-function
! (based on Sellers et al., 1992, J Geophys Res)
!
IF (io%LMEB_GNDRES) THEN
  pleg_delta(:) = zra_g_c(:) / ( zra_g_c(:) + exp(zrg_coef1 - zrg_coef2 * pwg(:) / zwsat(:) ) ) 
  plegi_delta(:) = zra_g_c(:) / ( zra_g_c(:) + exp(zrg_coef1 - zrg_coef2 * pwgi(:)/ zwsat(:) ) ) 
ELSE
  pleg_delta(:) = 1.0
  plegi_delta(:) = 1.0
ENDIF
!
! when hu*qsat < qa, there are two
! possibilities
!
! a) low-level air is dry, i.e., 
!    qa < qsat
!
! NOTE the additional delta fn's
! here are needed owing to linearization 
! of Qsat in the surface energy budget.
!
WHERE ( dk%XHUG(:)*pqsatg(:)  < pek%XQC(:) .AND. pqsatg(:) > pek%XQC(:) )
  dk%XHUG(:)        = pek%XQC(:) / pqsatg(:)
  pleg_delta(:)  = 0.0
END WHERE
WHERE ( phugi(:)*pqsatg(:) < pek%XQC(:) .AND. pqsatg(:) > pek%XQC(:) )
  phugi(:)       = pek%XQC(:) / pqsatg(:)
  plegi_delta(:) = 0.0
END WHERE
!
! b) low-level air is humid, i.e., qa >= qsat (condensation)
!
WHERE ( dk%XHUG*pqsatg  < pek%XQC(:) .AND. pqsatg <= pek%XQC(:) )dk%XHUG(:)  = 1.0
WHERE ( phugi*pqsatg < pek%XQC(:) .AND. pqsatg <= pek%XQC(:) )phugi(:) = 1.0
!
!
IF (lhook) CALL dr_hook('DRAG_MEB',1,zhook_handle)
!
END SUBROUTINE drag_meb