preps_for_meb_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 preps_for_meb_drag(LCVEL, LFORC_MEASURE,                   &
                             pz0, pz0h, pz0eff, ph_veg, pzref,           & 
                             ptc, pta, pqc, pqa, puref, pvmod,            &
                             pexna, pexns, pdircoszw, pdisph,             &
                             pvelc, pzvmod, pri, pra,                     &
                             pch,pcdn,pcd                                 )
!
! typical values for nordic forest:
!     PZ0 = 0.8 m
!     PH_VEG =   15 m
!
!
!!****  *PREPS_FOR_MEB*
!!
!!    PURPOSE
!!    -------
!
!     Calculates the winds near top of vegetation (PVELC)
!     vegetation (PVELC). Used for input to surface_air_meb.
!     Only used for double energy balance
!     Also calculates PRA,PCH,PCDN and PCD
!         
!     
!!**  METHOD
!!    ------
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      P. Samuelsson/S.Gollvik           * SMHI *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    11/2010
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
! LFORC_MEASURE= .TRUE. => heigh is given over ground
!                        .FALSE. => heigh is given over displacement height
!                         example from a numerical model PZREF=10 means that
!                         the level is 10 m above the displacement height
!
!
USE modd_csts,              ONLY : xpi, xkarman, xg, xcpd, xrd
USE modd_surf_atm,          ONLY : ldrag_coef_arp, xrimax
USE modd_isba_par,          ONLY : xlimh
!
USE modi_surface_aero_cond
USE modi_surface_cd
USE modi_surface_ri
USE modi_wind_threshold
!RJ: missing modi
USE modi_surface_cdch_1darp
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)   ::  pz0, pz0h, pz0eff,ph_veg, pzref, puref, pvmod
!                                     PZ0  = roughness length for momentum
!                                     PZ0H  = roughness length for heat
!                                     PZ0EFF  = roughness length for momentum (orographic)
!                                     PH_VEG = height of the vegetation
!                                     PZREF = height of the lowest model layer
!                                     PUREF = reference height of the wind 
!                                             NOTE this is different from PZREF
!                                             ONLY in stand-alone/forced mode,
!                                             NOT when coupled to a model (MesoNH)
!                                     PVMOD = module of the horizontal wind
!
REAL, DIMENSION(:), INTENT(IN)   ::  ptc, pta, pqc, pqa 
!                                     PTC = temperature of the canopy air
!                                     PTA = temperature of the lowest model layer
!                                     PQC = specific humidity of the canopy air
!                                     PQA = specific humidity of the lowest model level
!
REAL, DIMENSION(:), INTENT(IN)   ::  pexna, pexns, pdircoszw, pdisph
!                                     PEXNA = exner function at the lowest model level
!                                     PEXNS = exner function at the surface
!                                     PDIRCOSZW = Cosine of the angle between the normal
!                                                 to the surface and the vertical
!                                     PDISPH = displacement height
!
REAL, DIMENSION(:), INTENT(OUT)  ::  pvelc, pzvmod, pri, pra, pch, pcdn, pcd
!                                     PVELC  =  wind speed atr top of vegetation
!                                     PZVMOD = PMOD after wind threshold
!                                     PRI = Richardson number between canopy air and atm
!                                     PRA = aerodynamic resistance between canopy air
!                                           and lowest model layer
!                                     PCH = exchange coefficient for heat
!                                     PCDN = neutral exchange coefficient for momentum
!                                     PCD = exchange coefficient for momentum
!
LOGICAL, INTENT(IN)              ::   lcvel, lforc_measure
!                                       LCVEL = switch for calculating PVELC
!                                       LFORC_MEASURE = switch for using measured data
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PZ0)) :: zreveg, zcd, zcur, zucur, zratvv, zac
REAL, DIMENSION(SIZE(PZ0)) :: zbn, zbd, zcbnvv, zcbsvv, zcbuvv, zredvv
!
REAL(KIND=JPRB) :: zhook_handle
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER             :: zul      = 1.        ! typical windspeed within the foliage
REAL, PARAMETER             :: zny      = 0.15e-04  ! kinematic viscosity for air
REAL, PARAMETER             :: zvelclim = 0.1       ! safe minimum value ov PVELC
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!
IF (lhook) CALL dr_hook('PREPS_FOR_MEB_DRAG',0,zhook_handle)
!
pvelc(:) = zvelclim
pra(:)   = 0.
zcur(:)  = max(pzref(:),xlimh)
zucur(:) = max(puref(:),xlimh)
!
!Calculate height compared to dispacement height
!
IF(lforc_measure) THEN
!
   zcur(:)  = pzref(:)-pdisph(:)
   zucur(:) = puref(:)-pdisph(:)
!
   IF (any(zcur<0.0 .OR. zucur<0.0))THEN
     print *,'MAXVAL(PH_VEG)=',maxval(ph_veg)
     print *,'MAXVAL(PDISPH)=',maxval(pdisph)
     print *,'MINVAL(PZREF)=',minval(pzref)
     print *,'MINVAL(PUREF)=',minval(puref)
     stop "Forcing height for wind or temperature too low!!"
   ENDIF
!
ELSE
!
! When running a numerical model it is assumed that PZREF is the distance between the 
! lowest model layer and displacement height, thus the "trees are below surface"
!
   zcur(:)  = max(pzref(:),ph_veg(:)-pdisph(:)+xlimh)
   zucur(:) = max(puref(:),ph_veg(:)-pdisph(:)+xlimh)
!
ENDIF
!
! Exner function at displacement height
! For consistancy displacement height has the same pressure as the surface
!
 CALL surface_ri(ptc, pqc, pexns, pexna, pta, pqa,          &
                   zcur, zucur, pdircoszw, pvmod, pri )
!
pri(:) = min(pri(:),xrimax)
!
pzvmod = wind_threshold(pvmod,zucur)
!
IF (ldrag_coef_arp) THEN

   CALL surface_cdch_1darp(zcur, pz0eff, pz0h, pzvmod, pta, ptc, &
                           pqa, pqc, pcd, pcdn, pch              )
   pra(:) = 1. / ( pch(:) * pzvmod(:) )

ELSE
!
!               -------------------------------------------------
!
   CALL surface_aero_cond(pri, zcur, zucur, pzvmod, pz0, pz0h, zac, pra, pch)
!
!-------------------------------------------------------------------------------
!
!*       7.2    DRAG COEFFICIENT FOR MOMENTUM TRANSFERS
!               ---------------------------------------
!
!
   CALL surface_cd(pri, zcur, zucur, pz0eff, pz0h, pcd, pcdn)
!
ENDIF
!
IF(lcvel) THEN
!
! In cases with very low model layer, we estimate the wind at treetop as that of
! the lowest layer
!
   zratvv(:) = min((ph_veg(:)-pdisph(:))/zcur(:),1.)
!
   zbn(:)    = xkarman/sqrt(pcdn(:))
   zbd(:)    = xkarman/sqrt(pcd(:))
   zcbnvv(:) = alog(1.+(exp(zbn(:))-1.)*zratvv(:))
   zcbsvv(:) = -(zbn(:)-zbd(:))*zratvv(:)
   zcbuvv(:) = -alog(1.+(exp(zbn(:)-zbd(:))-1.)*zratvv(:))
!
   WHERE(pri>=0)
      zredvv(:) = (zcbnvv(:) + zcbsvv(:))/zbd(:)
   ELSEWHERE
      zredvv(:) = (zcbnvv(:) + zcbuvv(:))/zbd(:)
   END WHERE
!
   pvelc(:)     = max(zredvv(:)*pzvmod(:),zvelclim)
!
ENDIF
!
IF (lhook) CALL dr_hook('PREPS_FOR_MEB_DRAG',1,zhook_handle)
!
END SUBROUTINE preps_for_meb_drag