drag.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(HISBA, HSNOW_ISBA, HCPSURF, PTSTEP,                  &
                      ptg, pwg, pwgi,                                    &
                      pexns, pexna, pta, pvmod, pqa, prr, psr,           &
                      pps, prs, pveg, pz0, pz0eff, pz0h,                 &
                      pwfc, pwsat, ppsng, ppsnv, pzref, puref,           &
                      pdircoszw, pdelta, pf5, pra,                       &
                      pch, pcd, pcdn, pri, phug, phugi,                  &
                      phv, phu, pcps, pqs, pffg, pffv, pff,              &
                      pffg_nosnow, pffv_nosnow,                          &
                      pleg_delta, plegi_delta, pwr, prhoa, plvtt, pqsat  )  
!   ############################################################################
!
!!****  *DRAG*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the drag coefficients for heat and momentum transfers
!     over ground (i.e., Ch and Cd).
!         
!     
!!**  METHOD
!!    ------
!
!     1) computes hu, hv, and delta
!
!     2) use this to find qsoil, the grid-averaged relative humidity
!        of the soil
!
!     3) calculates the Richardson's number
!
!     4) find the transfer and resistance coefficients Ch, Cd, and Ra
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    USE MODD_CST 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Mascart et al. (1995)
!!    Belair (1995)
!!      
!!    AUTHOR
!!    ------
!!
!!      S. Belair           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    13/03/95 
!!      (J.Stein)   15/11/95  use the potential temperature to compute Ri
!!                            and PVMOD instead of ZVMOD
!!      (P.Lacarrere)15/03/96 replace * PEXNS by / PEXNS
!!      (V.Masson)   22/12/97 computation of z0eff after snow treatment
!!      (V.Masson)   05/10/98 clear routine
!!      (A.Boone)    11/26/98 Option for PDELTA: forested vs default surface
!!      (V.Masson)   01/03/03 Puts PDELTA in another routine
!!      (P.LeMoigne) 28/07/05 dependance on qs for cp
!!      (P.LeMoigne) 09/02/06 z0h and snow  
!!      (P.LeMoigne) 09/01/08 limitation of Ri
!!      (B.Decharme) 01/05/08 optional limitation of Ri
!!      (B.Decharme) 03/06/08 flooding scheme
!!      (A.Boone)    03/15/10 Add delta fnctions to force LEG ans LEGI=0
!!                            when hug(i)Qsat < Qa and Qsat > Qa
!!      (A.Boone)    21/11/11 Add Rs_max limit for dry conditions with Etv
!!      (B. Decharme)   09/12 limitation of Ri in surface_ri.F90
!!      (C. Ardilouze)  09/13 Halstead coef set to 0 for very low values
!!      (B. Decharme)   03/16 Bug in limitation of Er for Interception reservoir
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,     ONLY : xpi, xcpd, xcpv
USE modd_isba_par, ONLY : xwgmin, xrs_max
USE modd_surf_atm, ONLY : ldrag_coef_arp, lrrgust_arp, xrrscale,   &
                            xrrgamma, xutilgust, lcpl_arp  
!
USE modi_surface_ri
USE modi_surface_aero_cond
USE modi_surface_cd
USE modi_surface_cdch_1darp
USE modi_wind_threshold
!
USE mode_thermos
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
 CHARACTER(LEN=*),     INTENT(IN)  :: hisba      ! type of ISBA version:
!                                               ! '2-L' (default)
!                                               ! '3-L'
!                                               ! 'DIF'
 CHARACTER(LEN=*),     INTENT(IN)  :: hsnow_isba ! 'DEF' = Default F-R snow scheme
!                                               !         (Douville et al. 1995)
!                                               ! '3-L' = 3-L snow scheme (option)
!                                               !         (Boone and Etchevers 2000)
 CHARACTER(LEN=*),     INTENT(IN)  :: hcpsurf    ! option for specific heat Cp:
!                                               ! 'DRY' = dry Cp
!                                               ! 'HUM' = Cp as a function of qs
!
REAL,                 INTENT(IN) :: ptstep      ! timestep of ISBA (not split time step)
!
REAL, DIMENSION(:), INTENT(IN)   :: ptg, pwg, pwgi, pexns
!                                     PTG     = surface temperature
!                                     PWG     = near-surface volumetric water content
!                                     PWGI    = near-surface frozen volumetric water content
!                                     PEXNS   = Exner function at the surface
!
REAL, DIMENSION(:), INTENT(IN)   :: pexna, pta, pvmod, pqa, prr, psr, pps
!                                     PEXNA= Exner function
!                                     near surface atmospheric variables
!                                     PTA = 2m temperature
!                                     PVMOD = module of the horizontal wind
!                                             NOTE it should be input as >= 1. (m)
!                                     PQA = specific humidity
!                                     PPS = surface pressure
!                                     PRR = rain rate    
!                                     PSR = snow rate             
!
REAL, DIMENSION(:), INTENT(IN)   :: prs, pveg, pz0, pz0h, pz0eff
REAL, DIMENSION(:), INTENT(IN)   :: pwfc, pwsat, ppsng, ppsnv, pzref, puref
!                                     PRS = stomatal resistance
!                                     PVEG = vegetation fraction
!                                     PZ0  = roughness length for momentum
!                                     PZ0H = roughness length for heat
!                                     PZ0EFF= effective roughness length
!                                     PWFC = field capacity volumetric water content
!                                     PWSAT = volumetric water content at saturation
!                                     PPSNG = fraction of the bare ground covered
!                                             by snow
!                                     PPSNV = fraction of the vegetation covered
!                                             by snow
!                                     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)
!
REAL, DIMENSION(:), INTENT(INOUT) :: pdelta
!                                     PDELTA = fraction of the foliage covered
!                                              by intercepted water
REAL, DIMENSION(:), INTENT(IN)    :: pf5
!                                     PF5 = water stress function for Hv
REAL, DIMENSION(:), INTENT(IN)    :: pdircoszw 
! Cosinus of the angle between the normal to the surface and the vertical
REAL, DIMENSION(:), INTENT(INOUT) :: pra
!                                      aerodynamic surface resistance
REAL, DIMENSION(:), INTENT(INOUT) :: pcps
!                                      specific heat at surface
!
REAL, DIMENSION(:), INTENT(OUT)  :: pch, pcd, pcdn, pri
!                                     PCH = drag coefficient for heat
!                                     PCD = drag coefficient for momentum
!                                     PCDN= neutral drag coefficient for momentum
!                                     PRI = Richardson number
!
REAL, DIMENSION(:), INTENT(OUT)  :: phug, phugi, phv, phu, pqs, pleg_delta, plegi_delta
!                                     PHUG  = ground relative humidity
!                                     PHUGI = ground (ice) relative humidity
!                                     PHV = Halstead coefficient
!                                     PHU = relative humidity at the surface
!                                     PQS = humidity at surface
!                                     PLEG_DELTA  = soil evaporation delta fn
!                                     PLEGI_DELTA = soil sublimation delta fn
!
REAL, DIMENSION(:), INTENT(IN)   :: pffv, pff, pffg, pffg_nosnow, pffv_nosnow
!                                   PFFG = Floodplain fraction over ground
!                                   PFFV = Floodplain fraction over vegetation
!                                   PFF  = Floodplain fraction at the surface
!
REAL, DIMENSION(:), INTENT(IN)   :: pwr, prhoa, plvtt
!                                   PWR = liquid water retained on the foliage
!                                   PRHOA = near-ground air density
!                                   PLVTT = Vaporization heat
!
REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pqsat
!                                            PQSAT = specific humidity at saturation
!
!*      0.2    declarations of local variables
!
REAL, PARAMETER            :: zhvlim = 1.e-12 ! set halstead coef to 0 at this limit
!
REAL, DIMENSION(SIZE(PTG)) :: zqsat,           &
!                                              ZQSAT = specific humidity at saturation
                                 zac,            &
!                                              ZAC = aerodynamical conductance
                                 zwfc,           &
!                                              ZWFC = field capacity in presence of ice
                                 zvmod,          &
!                                              ZVMOD = wind modulus with minimum threshold
                                  zfp,           &
!                                              ZFP = working variable                               
                                 zzhv,           &
!                                              ZZHV = condensation delta fn for Hv
                                 zrrcor
!                                              ZRRCOR = correction of CD, CH, CDN due to moist-gustiness
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!               ---------------
!
IF (lhook) CALL dr_hook('DRAG',0,zhook_handle)
pch(:)   =0.
pcd(:)   =0.
phug(:)  =0.
phugi(:) =0.
phv(:)   =0.
phu(:)   =0.
!
zqsat(:)    =0.
zwfc(:)     =0.
pri(:)      =0.
!
zvmod = wind_threshold(pvmod,puref)
!
zqsat(:) = qsat(ptg(:),pps(:)) 
!
IF(present(pqsat))pqsat(:)=zqsat(:)
!
!-------------------------------------------------------------------------------
!
!*       1.     RELATIVE HUMIDITY OF THE GROUND HU
!               ----------------------------------
!
!                                         this relative humidity is related to
!                                         the superficial soil moisture and the
!                                         field capacity of the ground
!
!
zwfc(:)  = pwfc(:)*(pwsat(:)-pwgi(:))/pwsat(:)
!
phug(:)  = 0.5 * ( 1.-cos(xpi*min((pwg(:)-xwgmin) /zwfc(:),1.)) )

zwfc(:)  = pwfc(:)*max(xwgmin, pwsat(:)-pwg(:))/pwsat(:)
!
phugi(:) = 0.5 * ( 1.-cos(xpi*min(pwgi(:)/zwfc(:),1.)) )
!
!
!                                         there is a specific treatment for dew
!                                         (see Mahfouf and Noilhan, jam, 1991)
!                                         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.
!
pleg_delta(:)    = 1.0
plegi_delta(:)   = 1.0
WHERE ( phug*zqsat < pqa .AND. zqsat > pqa )
  phug(:)  = pqa(:) / zqsat(:)
  pleg_delta(:)  = 0.0
END WHERE
WHERE ( phugi*zqsat < pqa .AND. zqsat > pqa )
  phugi(:) = pqa(:) / zqsat(:)
  plegi_delta(:) = 0.0
END WHERE
!
!                                         b) low-level air is humid, i.e.,
!                                            qa >= qsat
!
WHERE ( phug*zqsat  < pqa .AND. zqsat <= pqa )
  phug(:)  = 1.0
END WHERE
WHERE ( phugi*zqsat < pqa .AND. zqsat <= pqa )
  phugi(:) = 1.0
END WHERE
!
!-------------------------------------------------------------------------------
!
!*       3.     HALSTEAD COEFFICIENT (RELATIVE HUMIDITY OF THE VEGETATION)
!               ----------------------------------------------------------
!
zzhv(:)       = max(0.,sign(1.,zqsat(:)-pqa(:))) ! condensation on foilage if = 1
!
phv(:)        = pdelta(:) + (1.- pdelta(:))*                                    &
                ( pra(:) + prs(:)*(1.0 - zzhv(:)) )*                            &
                ( (1/(pra(:)+prs(:))) - (zzhv(:)*(1.-pf5(:))/(pra(:)+xrs_max)) )
!
WHERE(phv(:)<zhvlim)
      phv(:)=0.0
ENDWHERE
!
!-------------------------------------------------------------------------------
!
!*       4.     GRID-AVERAGED HUMIDITY OF THE SURFACE
!               -------------------------------------
!
IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
!
! For explicit snow, humidity HERE only over snow-free region:
!
   pqs(:) =( (1.-pveg(:))*(1.-pffg_nosnow(:))*(pwg(:) /(pwg(:)+pwgi(:)))*phug(:)   &
             + (1.-pveg(:))*(1.-pffg_nosnow(:))*(pwgi(:)/(pwg(:)+pwgi(:)))*phugi(:)   &
             +     pveg(:) *(1.-pffv_nosnow(:))*    phv(:)                            &
             + ((1.-pveg(:))*pffg_nosnow(:)+pveg(:)*pffv_nosnow(:)))*zqsat(:)         &
             +     pveg(:) *(1.-pffv_nosnow(:))*(1.-phv(:))*pqa(:)  
!
ELSE
!
   pqs(:) =( (1.-pveg(:))*(1.-ppsng(:)-pffg(:))*(pwg(:) /(pwg(:)+pwgi(:)))*phug(:)  &
             + (1.-pveg(:))*(1.-ppsng(:)-pffg(:))*(pwgi(:)/(pwg(:)+pwgi(:)))*phugi(:)  &
             + (1.-pveg(:))*    ppsng(:)                                               &
             +     pveg(:) *(1.-ppsnv(:)-pffv(:))*phv(:)                               &
             +     pveg(:) *    ppsnv(:)                                               &
             +     pff(:)                            )*zqsat(:)                        &
             +     pveg(:) *(1.-ppsnv(:)-pffv(:))*(1.-phv(:))*pqa(:)  
!
ENDIF
!
phu(:) = pqs(:)/zqsat(:)
!
!-------------------------------------------------------------------------------
!
!*       5.     SPECIFIC HEAT OF THE SURFACE
!               ----------------------------
!
IF (hcpsurf=='DRY') THEN
    pcps(:) = xcpd
ELSEIF(.NOT.lcpl_arp)THEN
    pcps(:) = xcpd + ( xcpv - xcpd ) * pqa(:)   
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       6.     COMPUTATION OF RESISTANCE AND DRAG COEFFICIENT
!               ----------------------------------------------
!
 CALL surface_ri(ptg, pqs, pexns, pexna, pta, pqa,                    &
                pzref, puref, pdircoszw, pvmod, pri                  )  
!
!-------------------------------------------------------------------------------
!
!*       6.5    LIMITATION OF RICHARDSON NUMBER
!               -------------------------------
!
!Now done directly in SURFACE_RI subroutine
!
!-------------------------------------------------------------------------------
!*       7.0     DRAG COEFFICIENT FOR HEAT AND MOMENTUM TRANSFERS
!               -------------------------------------------------
!
IF (ldrag_coef_arp) THEN

   CALL surface_cdch_1darp(pzref, pz0eff, pz0h, zvmod, pta, ptg, &
                             pqa, pqs, pcd, pcdn, pch              )  
   pra(:) = 1. / ( pch(:) * zvmod(:) )

ELSE

!
!*       7.1    SURFACE AERODYNAMIC RESISTANCE FOR HEAT TRANSFERS
!               -------------------------------------------------
!
   CALL surface_aero_cond(pri, pzref, puref, zvmod, pz0, pz0h, zac, pra, pch)
!
!-------------------------------------------------------------------------------
!
!*       7.2    DRAG COEFFICIENT FOR MOMENTUM TRANSFERS
!               ---------------------------------------
!
!
   CALL surface_cd(pri, pzref, puref, pz0eff, pz0h, pcd, pcdn)
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       8.     CORRECTION DUE TO MOIST GUSTINESS
!               ---------------------------------
!
IF (lrrgust_arp) THEN

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

   pcd  = pcd  * zrrcor
   pch  = pch  * zrrcor
   pcdn = pcdn * zrrcor

ENDIF
!
!-------------------------------------------------------------------------------
!
!
!*       9.     HALSTEAD COEFFICIENT (WITH THE NEW VALUES OF PRA)
!               ----------------------------------------------------------
!
!
phv(:)        = pdelta(:) + (1.- pdelta(:))*                                       &
                ( pra(:) + prs(:)*(1.0 - zzhv(:)) )*                               &
                ( (1/(pra(:)+prs(:))) - (zzhv(:)*(1.-pf5(:))/(pra(:)+xrs_max)) )
!
 CALL limit_ler(pdelta)
!
WHERE(phv(:)<zhvlim)
      phv(:)=0.0
ENDWHERE
!                
IF (lhook) CALL dr_hook('DRAG',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
 CONTAINS
!
!-------------------------------------------------------------------------------
!
SUBROUTINE limit_ler(PDELTA)
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(INOUT) :: pdelta
!
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PDELTA)) ::  zpsnv
REAL, DIMENSION(SIZE(PDELTA)) ::  zlev
REAL, DIMENSION(SIZE(PDELTA)) ::  zletr
REAL, DIMENSION(SIZE(PDELTA)) ::  zlecoef
REAL, DIMENSION(SIZE(PDELTA)) ::  zer
REAL, DIMENSION(SIZE(PDELTA)) ::  zrrveg
REAL, DIMENSION(SIZE(PDELTA)) ::  zwr_delta
!
REAL(KIND=JPRB)   :: zhook_handle
!
IF (lhook) CALL dr_hook('DRAG:LIMIT_LER',0,zhook_handle)
!
!-------------------------------------------------------------------------------
!
!*       1.     Initialization:
!               ---------------
!
IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
   zlecoef(:) = (1.0-ppsnv(:)-pffv(:))
   zpsnv(:) = 0.0
ELSE
   zlecoef(:) = 1.0
   zpsnv(:) = ppsnv(:)+pffv(:)
ENDIF
!
!
!*       2.     Interception reservoir consistency:
!               -----------------------------------
!
!  In DRAG, we use the timestep of ISBA (PTSTEP) and not the split time step (ZTSTEP)
!  because diagnostic canopy evaporation (Er) must be consistent with PWR to
!  limit negative dripping in hydro_veg

!
zlev(:)  = prhoa(:) * plvtt(:) * pveg(:) * (1-zpsnv(:)) * phv(:) * (zqsat(:) - pqa(:)) / pra(:)
!
zletr(:) = zzhv(:) * prhoa(:) * (1. - pdelta(:)) * plvtt(:) * pveg(:) *(1-zpsnv(:))         &
         * (zqsat(:) - pqa(:)) * ( (1/(pra(:) + prs(:))) - ((1.-pf5(:))/(pra(:) + xrs_max)) )
!
zer(:)=ptstep*(zlev(:)-zletr(:))*zlecoef(:)/plvtt(:)
!
zrrveg(:) = ptstep*pveg(:)*(1.-ppsnv(:))*prr(:)
!
zwr_delta(:)=1.0
!
WHERE( zzhv(:)>0.0 .AND. zer(:)/=0.0 .AND. (pwr(:)+zrrveg(:))<zer(:) )
!
       zwr_delta(:) = max(0.01,min(1.0,(pwr(:)+zrrveg(:))/zer(:)))
!       
       pdelta(:) = pdelta(:) * zwr_delta(:)
!
       phv(:)    = pdelta(:) + (1.- pdelta(:))*( pra(:) + prs(:)*(1.0 - zzhv(:)) )* &
                  ( (1/(pra(:)+prs(:))) - (zzhv(:)*(1.-pf5(:))/(pra(:)+xrs_max)) )
!
ENDWHERE
!
IF (lhook) CALL dr_hook('DRAG:LIMIT_LER',1,zhook_handle)
!
END SUBROUTINE limit_ler
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE drag