surface_air_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 surface_air_meb(PZ0, PZ0H, PZ0G, PH_VEG, PLAI,          &
                               ptg, ptc, ptv, pvelc, plw,              &
                               pdisph,                                 &
                               pragnc, pgvnc,                          &
                               pustar2, pcd, pch, pri                  )
!
! typical values for nordic forest:
!     PZ0 = 0.8 m
!     PZ0G = 0.007m
!     PH_VEG =   15 m
!     PLAI =  5.
!     PCHIL =   0.12
!     PLW   =   0.02 m
!
!
!!****  *SURFACE_AIR_MEB*
!!
!!    PURPOSE
!!    -------
!
!     Calculates the aerodynamic resistance (PRAGNC) between the soil and eventually
!     understory vegetation and canopy air, based on Choudhury and Monteith (1998)
!     and Monteith (1975). 
!     Modification for free convection based on Sellers et.al (1986), leaf drag 
!     coefficient according to Sellers et.al. (1996).
!     Also calculates the aerodynamic conductance (PGVNC) between the canopy itself
!     and canopy air, also based on Choudhury and Monteith (1998), and modified for
!     free conveection by Sellers et.al. (1986)
!     Only used for double energy balance
!         
!     
!!**  METHOD
!!    ------
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      P. Samuelsson/S.Gollvik           * SMHI *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original       11/2010
!!     (A. Boone)   25/06/2014  Use stability fn from ISBA for stable conditions
!!                              (since none existed before). Added to handle extremely 
!!                              stable conditions, such as encountered over a snowpack
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,     ONLY : xg, xkarman
USE modd_surf_atm, ONLY : xrimax
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)   :: pz0, pz0h, pz0g, ph_veg, plai
!                                     PZ0  = roughness length for momentum
!                                     PZ0H = roughness length for heat
!                                     PZ0G = roughness length for soil(understory veg?)/snow
!                                     PH_VEG = height of the vegetation
!                                     PLAI = leaf area index
!
REAL, DIMENSION(:), INTENT(IN)   :: ptg, ptc, ptv, pvelc
!                                     PTG     = surface temperature
!                                     PTC     = canopy air temperature
!                                     PTV     = canopy temperature
!                                     PVELC  =  wind speed at top of vegetation
!
REAL, DIMENSION(:), INTENT(IN)   :: plw, pdisph
!                                     PLW = leaf width
!                                     PDISPH = displacement height
!
REAL, DIMENSION(:), INTENT(OUT)  :: pragnc, pgvnc
!                                     PRAGNC = aerodynamic resistance between
!                                              soil/snow and canopy air
!                                     PGVNC = aerodynamic conductance between
!                                              the canopy and canopy air
!
REAL, DIMENSION(:), INTENT(OUT)  :: pustar2, pcd, pch, pri
!                                   PUSTAR2 = canopy top friction velocity squared (m2/s2)
!                                   PCD     = drag coefficient (-)
!                                   PCH     = heat transfor coefficient (ground to canopy air) (-)
!                                   PRI     = Richardson number (ground to canopy air) (-)
!
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PTG)) :: zdiffh, zk, zpih, zdifft, zrif, zz0hg, zustar
!
REAL(KIND=JPRB) :: zhook_handle
!
!*      0.3    declarations of local parameters
!
REAL, PARAMETER             :: zny          = 0.15e-04 !kinematic viscosity for air
REAL, PARAMETER             :: zalpha       = 4.       !attenuation coef. for mom.
                                                       ! ~2 for Wheat, ~4 for coniferous forest
                                                       ! NOTE Eventually should possibly be made a fn of veg type
REAL, PARAMETER             :: zalphaprim   = 3.       !attenuation coef. for wind
REAL, PARAMETER             :: za           = 0.01
REAL, PARAMETER             :: zrafa        = 9.       !resistance factor for stability correction
                                                       !for unstable conditions

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!               ---------------
!
IF (lhook) CALL dr_hook('SURFACE_AIR_MEB',0,zhook_handle)
!
pragnc(:)  = 0.
pgvnc(:)   = 0.
!
zdiffh(:)  = max(ph_veg(:)-pdisph(:),0.1)
!
zk(:)      = (xkarman*xkarman)*pvelc(:)*zdiffh(:)/log(zdiffh(:)/pz0(:))
!
! Just a diagnostic: Ustar and the equivalent drag coef at the top of the canopy (m/s):
!
zustar(:)  = zk(:)/(xkarman*zdiffh(:))
pcd(:)     = zustar(:)/max(1.,pvelc(:))
pustar2(:) = zustar(:)**2
!
! Aerodynamic resistance, Eq. 25 Choudhury and Monteith, 1988:
!
pragnc(:)=ph_veg(:)*exp(zalpha)/(zalpha*zk(:))*(exp(-zalpha*pz0g(:)/ph_veg(:)) &
          -exp(-zalpha*(pdisph(:)+pz0(:))/ph_veg(:)))
!
! Modify the aerodynamic resistance, with an unstable transfer correction 
! Eq. A15, Sellers et.al. 1986 (RI < 0)
!
! For stable conditions (RI > 0), use a form which asymptotically approaches 0
! as Ri==>infinty. This is a fairly generic form typical of such curves...
! As we have no predefined stable correction, we use ISBAs
! (see Noilhan and Mahfouf, 1996). This stable correction
! is most critical for snow...BOTH sublimation from ground-based snowpack and
! melt conditions which can be characterized by very
! strong stable conditions, thus a limited, constant value of the turbulent
! exchange can permit very large negative heat fluxes (towards the snow)
! and excessive melt rates ...also excessive sublimation can occur in winter.
! Since the ISBA stability fn has a roughness length factor, we use
! a simple Ri-based weight factor to ensure a continuous transition
! between stability regimes. In the limit as Ri==>Ricrit, the 
! transfer function collapses into the ISBA relationship with the full roughness factor.
! Finally, note that we assume the ratio z0/z0h is the same for vegetation and 
! underlying surface (as done for composite ISBA).
!
pri(:)   = -xg*ph_veg(:)*(ptg(:)-ptc(:))/(ptg(:)*pvelc(:)*pvelc(:))
!
zrif(:)  = 0.
zz0hg(:) = pz0g(:)
!
WHERE(pri(:) <= 0.)                                        ! Unstable - Sellers 
   zpih(:) = sqrt(1. - zrafa*pri(:))                       
ELSEWHERE                                                  ! Stable   - Noilhan and Mahfouf
   zz0hg(:)= pz0g(:)*pz0h(:)/pz0(:)                        
   zrif(:) = min(1., pri(:)/xrimax)
   zpih(:) = (1/(1. + 15*pri(:)*sqrt(1.+5*pri(:))))* &  
             ((1.-zrif(:)) + zrif(:)*(log(zdiffh(:)/pz0g(:))/log(zdiffh(:)/zz0hg(:))) )            ! continuous
END WHERE
!
pragnc(:)  = pragnc(:)/zpih(:)
!
! Diagnose CH (transfer coefficient: -) just for diagnostic purposes:
!
pch(:)     = 1/(pragnc(:)*max(1., pvelc(:)))
!
! Aerodynamic resistance within the canopy layer, i.e. between canopy and canopy air
! Eq. 29 and 30, Choudhury and Monteith, 1988:
!
pgvnc(:)  = (2.*za)*plai(:)/zalphaprim*sqrt(pvelc(:)/plw(:))*(1.-exp(-0.5*zalphaprim)) 
!
! Limit the aerodynamic conductance to a small value
!
pgvnc(:)  = max(1.e-06,pgvnc(:))
!
! Free convection correction Eq. A9 Sellers et.al., 1986: 
!
zdifft(:) = max(1.e-06,ptv(:)-ptc(:))
zdifft(:) = sqrt(sqrt(zdifft(:)/plw(:)))
pgvnc(:)  = pgvnc(:)+zdifft(:)*plai(:)/890.
!
IF (lhook) CALL dr_hook('SURFACE_AIR_MEB',1,zhook_handle)
!
END SUBROUTINE surface_air_meb