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(LFORC_MEASURE,                                     &
                    ptg, ptc, ptv, psnowtemp, pta, pqc, 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,                  &
                    pz0_with_snow, pz0h_with_snow, pz0eff,             &
                    psnowswe,                                          &
                    pwr, pchip, ptstep, prs_vg, prs_vn,                &
                    ppsn, ppalphan, pzref, puref, ph_veg, pdircoszw,   &
                    ppsncv, pdelta, plai, omeb_gndres,                 &
                    pch, pcd, pcdn, pri, pra, pvelc,                   &
                    pcdsnow, pchsnow, prisnow, pustar2snow,            &
                    phug, phugi, phv, phvg, phvn, phu, pqs, prs,       &
                    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_csts,     ONLY : xpi
USE modd_snow_par, ONLY : xz0sn
USE modd_isba_par, ONLY : xwgmin, xflxmax
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
!
LOGICAL, INTENT(IN)              ::   lforc_measure
!                                       LFORC_MEASURE = switch for using measured data
!
REAL, INTENT(IN)                 ::   ptstep 
!                                     PTSTEP = Model time step (s)        
!
REAL, DIMENSION(:), INTENT(IN)   :: ptg, ptc, ptv, psnowtemp, pta, pqc, pqa, pvmod, &
                                    pwg, pwgi, pwsat, pwfc, pexns, pexna, pps,      &
                                    psnowswe
!                                     PTG     = surface temperature
!                                     PTC     = canopy air temperature
!                                     PTV     = canopy temperature
!                                     PSNOWTEMP = sfc layer snow temperature (K)
!                                     PSNOWSWE  = sfc layer snow water equiv (SWE) kg m-2
!                                     PTA     = temperature of the atmosphere/forcing-level
!                                     PQC     = canopy air humidity
!                                     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, &
                                    pz0_with_snow, pz0h_with_snow, pz0eff
!
!                    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
!                    PZ0EFF_MEBV        ! roughness length for momentum over MEB vegetation part of patch
!                                         eventuelly including orograhic roughness
!                    PZ0_MEBN           ! roughness length for momentum over MEB snow part of patch
!                    PZ0H_MEBN          ! roughness length for heat over MEB snow part of path
!                    PZ0EFF_MEBN        ! roughness length for momentum over MEB snow part of patch
!                                         eventuelly including orograhic roughness
!                    PZ0_WITH_SNOW      ! roughness length for momentum over MEB total patch
!                    PZ0H_WITH_SNOW     ! roughness length for heat over MEB total path
!                    PZ0EFF             ! roughness length for momentum over MEB total patch
!                                         eventuelly including orograhic roughness
!
REAL, DIMENSION(:), INTENT(IN)   :: pwr,pchip
!                                     PWR  = intercepted water for the canopy
!                                     PCHIP = view factor (for LW) 
!
REAL, DIMENSION(:), INTENT(IN)   :: prs_vg, prs_vn, ppsn, ppsncv,  &
                                    ppalphan, pzref, puref, ph_veg, pdircoszw
!
!                                     PRS_VG   = surface resistance for canopy (1-png)
!                                     PRS_VN   = surface resistance for canopy (png)
!                                     PPSN     = fraction of the patch covered by snow
!                                     PPSNCV   = 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)
!                                     PH_VEG   = height of the vegetation
!                                     PDIRCOSZW= Cosinus of the angle between 
!                                                the normal to the surface and the vertical
!
REAL, DIMENSION(:), INTENT(IN)    :: pdelta, plai
!                                     PDELTA   = fraction of the canopy foliage covered
!                                                by intercepted water (-)
!                                     PLAI     = vegetation LAI (m2 m-2)
!
LOGICAL,              INTENT(IN)    :: omeb_gndres   ! Flag for ground resistance
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)  :: pch, pcd, pcdn, pri, pra, pvelc
!                                     PCH = drag coefficient for heat, averaged over whole patch
!                                     PCD = drag coefficient for momentum, averaged over whole patch
!                                     PCDN= neutral drag coefficient for momentum, 
!                                           averaged over whole patch
!                                     PRI = Richardson number, averaged over whole patch
!                                     PRA = aerodynamic resistance between 
!                                           the whole patch and atmosphere
!                                     PVELC = wind speed at top of vegetation
!
REAL, DIMENSION(:), INTENT(OUT)  :: phug, phugi, phvg, phvn, phv,                        &
                                    phu, pqs, pleg_delta, plegi_delta, phsgl, phsgf, prs
!
!                                    PRS   = total effective stomatal resistance (under and overstory) diagnostic (s m-1)
!                                    PHUG  = ground relative humidity
!                                    PHUGI = ground (ice) relative humidity
!                                    PHV   = Total effective Halstead coefficient 
!                                    PHVN  = Halstead coefficient vegetation canopy above snow
!                                    PHVG  = Halstead coefficient vegetation canopy above ground
!                                    PHU   = effective relative humidity at the surface
!                                    PQS   = effective specific humidity at surface
!                                    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 PTV
!                                    PQSATC = qsat for PTC
!                                    PQSATN = qsat for PSNOWTEMP
!
REAL, DIMENSION(:), INTENT(OUT)  :: pcdsnow, pchsnow, prisnow, pustar2snow 
!                                    PCDSNOW     = drag coefficient over snow (-)
!                                    PCHSNOW     = heat/mass exchange coefficient over snow (-)
!                                    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)
!
pch(:)   = 0.
pcd(:)   = 0.
pcdn(:)  = 0.
pri(:)   = 0.
phug(:)  = 0.
phugi(:) = 0.
phv(:)   = 0.
phvn(:)  = 0.
phvg(:)  = 0.
phu(:)   = 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(:)
phug(:)  = 0.5 * ( 1.-cos(xpi*min((pwg(:)-xwgmin) /zwfc(:),1.)) )
zwsat(:) = max(xwgmin, zwsat(:))
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(ptv(:),pps(:)) 
pqsatc(:) = qsat(ptc(:),pps(:)) 
pqsatn(:) = qsati(psnowtemp(:),pps(:)) 
!
!-------------------------------------------------------------------------------
!
!*       2.     GRID-AVERAGED HUMIDITY OF THE SURFACE
!               -------------------------------------
!       average of canopy air humidity and saturation humidity of snowpack
!
!
   pqs(:) =(1.-ppsn(:)*ppalphan(:))*pqc(:)+ ppsn(:)*ppalphan(:)*pqsatn(:)
   phu(:) =(1.-ppsn(:)*ppalphan(:))*pqc(:)/pqsatc(:)+ ppsn(:)*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,plai,zlw,ph_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.,lforc_measure,            &
                   pz0_mebv, pz0h_mebv, pz0eff_mebv,     &
                   ph_veg, pzref,                        &
                   ptc, pta, pqc, 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, ph_veg, plai,   &
                     ptg, ptc, ptv, pvelc, zlw,                              &
                     zdisph,                                                 &
                     zra_g_c, zg_vg_c,                                       &
                     zustar2g, zcdg, zchg, zrig                              )
!
!Compute the lai of the canopy that is above snow
!
zlaisn(:)=plai(:)*(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, ph_veg, zlaisn,     &
                     psnowtemp, ptc, ptv, pvelc, zlw,                &
                     zdisph,                                         &
                     zra_n_c, zg_vn_c,                               &
                     zustar2g, zcdg, zchg, zrig                      )

! save values over snow for diagnostic purposes:

pustar2snow(:) = zustar2g(:)
pcdsnow(:)     = zcdg(:)  
pchsnow(:)     = 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.,lforc_measure,                           &
                   pz0_mebn, pz0h_mebn, pz0eff_mebn,                     &
                   ph_veg, pzref,                                        &
                   psnowtemp, 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.-ppsn(:)*ppalphan(:))*pflxc_mom(:) +  &
            ppsn(:)*ppalphan(:)*zcdnn(:)*zvmod(:)*prhoa(:)
!
! Now calculate the aerodynamic resistance for  the whole meb-patch
!
! Calculate  the effective temperature
!
zpsna(:)   = ppsn(:)*ppalphan(:)
!
zteff(:)   = (1.-zpsna(:))*ptc(:)+ zpsna(:)*psnowtemp(:)
!
! Some additional diagnostics:
!
pustar2snow(:) = (1.-zpsna(:))*pustar2snow(:) + zpsna(:)*zcdnn(:)*zvmod(:)**2
pcdsnow(:)     = (1.-zpsna(:))*pcdsnow(:)     + zpsna(:)*zcdnn(:)
pchsnow(:)     = (1.-zpsna(:))*pchsnow(:)     + zpsna(:)*zchnn(:)
prisnow(:)     = (1.-zpsna(:))*prisnow(:)     + zpsna(:)*zrinn(:)
!
!-------------------------------------------------------------------------------
!
 CALL preps_for_meb_drag(.false.,lforc_measure,            & 
                   pz0_with_snow, pz0h_with_snow, pz0eff, &
                   ph_veg, pzref,                         &
                   zteff, pta, pqs, pqa, puref, pvmod,    &
                   pexna, pexns, pdircoszw, zdisph,       &
                   zvelc, zvmod, pri, pra,                &
                   pch,pcdn,pcd                           )
!
IF (lrrgust_arp) THEN

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

   pcd(:)  = pcd(:)  * zrrcor(:)
   pch(:)  = pch(:)  * zrrcor(:)
   pcdn(:) = pcdn(:) * 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,
! limit exchange coefficients to 5000.
!
WHERE(zra_g_c(:) > zraeps)
   pflxc_g_c(:) = prhoa(:)  / zra_g_c(:)
   phsgl(:)     = zra_g_c(:)/(zra_g_c(:)+zrsgl(:))
   phsgf(:)     = zra_g_c(:)/(zra_g_c(:)+zrsgf(:))
ELSEWHERE
   pflxc_g_c(:) = xflxmax
   phsgl(:)     = 1.
   phsgf(:)     = 1.
END WHERE
!
WHERE(zra_n_c>zraeps)
   pflxc_n_c(:) = prhoa(:)/zra_n_c(:)
ELSEWHERE
   pflxc_n_c(:) = xflxmax
END WHERE
!
! 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.-ppsn(:)*ppalphan(:))*prr(:)+ pwr(:)/ptstep
zdenom(:)    = (1.-ppsncv(:))*xkdelta_wr*                                               &
               ( ppsn(:)*(1.-ppalphan(:))*pflxc_vn_c(:) + (1.-ppsn(:))*pflxc_vg_c(:) )* &
               ( pqsatv(:)-pqc(:))
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(:)-pqc(:)))             &
          *(1.-zdeltav(:))*prs_vg(:)*zg_vg_c(:) / (1.+prs_vg(:)*zg_vg_c(:))
!
phvn(:) = 1. - max(0.,sign(1.,pqsatv(:)-pqc(:)))             &
          *(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: 
!
prs(:)  = ppalphan(:)*prs_vn(:) + (1.0-ppalphan(:))*prs_vg(:)
phv(:)  = 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 (omeb_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 ( phug(:)*pqsatg(:)  < pqc(:) .AND. pqsatg(:) > pqc(:) )
  phug(:)        = pqc(:) / pqsatg(:)
  pleg_delta(:)  = 0.0
END WHERE
WHERE ( phugi(:)*pqsatg(:) < pqc(:) .AND. pqsatg(:) > pqc(:) )
  phugi(:)       = pqc(:) / pqsatg(:)
  plegi_delta(:) = 0.0
END WHERE
!
! b) low-level air is humid, i.e., qa >= qsat (condensation)
!
WHERE ( phug*pqsatg  < pqc .AND. pqsatg <= pqc )phug(:)  = 1.0
WHERE ( phugi*pqsatg < pqc .AND. pqsatg <= pqc )phugi(:) = 1.0
!
!
IF (lhook) CALL dr_hook('DRAG_MEB',1,zhook_handle)
!
END SUBROUTINE drag_meb