canopy_evol_neutral.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 canopy_evol_neutral(KI,KLVL,PTSTEP,KIMPL,PWIND,PRHOA, &
                             psflux_u,                                 &
                             pforc_u,pdforc_udu,pforc_e,pdforc_ede,    &
                             pz,pzf,pdz,pdzf,pu,ptke,pustar,           &
                             palfau,pbetau                             )
!     #########################################
!
!!****  *CANOPY_EVOL* - evolution of canopy
!!                        
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    05/2010 
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_csts,        ONLY : xkarman
USE modd_surf_par,    ONLY : xundef
USE modd_canopy_turb, ONLY : xcmfs, xalpsbl, xced
!
USE modi_canopy_evol_wind
USE modi_canopy_evol_tke
!
USE mode_sbls
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
INTEGER,                  INTENT(IN)    :: ki        ! number of horizontal points
INTEGER,                  INTENT(IN)    :: klvl      ! number of levels in canopy
REAL,                     INTENT(IN)    :: ptstep    ! atmospheric time-step                 (s)
INTEGER,                  INTENT(IN)    :: kimpl     ! implicitation code: 
!                                                    ! 1 : computes only alfa and beta coupling
!                                                    !     coefficients for all variables
!                                                    ! 2 : computes temporal evolution of the
!                                                    !     variables
REAL, DIMENSION(KI),      INTENT(IN)    :: pwind     ! wind speed                            (m/s)
REAL, DIMENSION(KI),      INTENT(IN)    :: prhoa     ! Air density at forcing level          (kg/m3)
REAL, DIMENSION(KI),      INTENT(IN)    :: psflux_u  ! surface flux u'w'                     (m2/s2)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pforc_u   ! tendency of wind due to canopy drag   (m/s2)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdforc_udu! formal derivative of the tendency of
!                                                    ! wind due to canopy drag               (1/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pforc_e   ! tendency of TKE  due to canopy drag   (m2/s3)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdforc_ede! formal derivative of the tendency of
!                                                    ! TKE  due to canopy drag               (1/s)
!
REAL, DIMENSION(KI,KLVL), INTENT(INOUT) :: pz        ! heights of canopy levels              (m)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pzf       ! heights of bottom of canopy levels    (m)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdz       ! depth   of canopy levels              (m)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdzf      ! depth between canopy levels           (m)
REAL, DIMENSION(KI,KLVL), INTENT(INOUT) :: pu        ! wind speed at canopy levels           (m/s)
REAL, DIMENSION(KI,KLVL), INTENT(INOUT) :: ptke      ! TKE at canopy levels                  (m2/s2)
REAL, DIMENSION(KI),      INTENT(OUT)   :: pustar    ! friction velocity at forcing level    (m/s)
REAL, DIMENSION(KI),      INTENT(OUT)   :: palfau   ! V+(1) = alfa u'w'(1) + beta
REAL, DIMENSION(KI),      INTENT(OUT)   :: pbetau   ! V+(1) = alfa u'w'(1) + beta
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER :: jlayer                 ! loop counter on layers
INTEGER :: ji                     ! loop counter
!
REAL, DIMENSION(KI,KLVL) :: zk       ! mixing coefficient
REAL, DIMENSION(KI,KLVL) :: zdkddvdz ! formal derivative of mixing coefficient according to variable vertical gradient
REAL, DIMENSION(KI,KLVL) :: zexn     ! Exner function        at full levels
REAL, DIMENSION(KI,KLVL) :: zuw      ! friction at mid levels
REAL, DIMENSION(KI,KLVL) :: zth      ! temperature(meaningless here)
REAL, DIMENSION(KI,KLVL) :: zwth     ! heat  flux (supposed equal to zero)
REAL, DIMENSION(KI,KLVL) :: zwq      ! vapor flux (supposed equal to zero)
REAL, DIMENSION(KI,KLVL) :: zleps    ! Dissipative length at full levels
REAL, DIMENSION(KI,KLVL) :: zz       ! height above displacement height
REAL, DIMENSION(KI,KLVL) :: zlm      ! Mixing      length at mid levels
REAL, DIMENSION(KI)      :: zustar   ! friction velocity (estimated from
!                                    ! fluxes at atmospheric forcing level)
!
REAL                     :: zz0      ! a value of z0 just for first time-step init.
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
!*    1. First time step initialization
!        ------------------------------
!
!* first time step (i.e. wind profile is initially zero) : 
!  set a neutral wind profile in relation with forcing wind
IF (lhook) CALL dr_hook('CANOPY_EVOL_NEUTRAL',0,zhook_handle)
DO ji=1,ki
  IF(pwind(ji)>0. .AND. pu(ji,klvl-1)==0.) THEN
    zz0 = pz(ji,1)/2.
    pu(ji,:) = pwind(ji) * log(pz(ji,:)/zz0) / log(pz(ji,klvl)/zz0)
  END IF
END DO
!-------------------------------------------------------------------------------
!
!*    2. mixing and dissipative lengths (at full levels)
!        ------------------------------
!
!* influence of ground surface on mixing length in neutral case
DO jlayer=1,klvl
  zlm(:,jlayer) = xkarman/sqrt(xalpsbl)/xcmfs * pz(:,jlayer)
  zleps(:,jlayer) = xkarman*(xalpsbl**1.5)*xced * pz(:,jlayer)
END DO
!
!-------------------------------------------------------------------------------
!
!*    3. time evolution of wind in canopy
!        --------------------------------
!
!*    3.1 mixing coefficient for wind at mid layers (at half levels)
!         ---------------------------
!
zk = -999.
DO jlayer=2,klvl
  zk(:,jlayer) = 0.5 * xcmfs * zlm(:,jlayer)   * sqrt(ptke(:,jlayer)  ) &
               + 0.5 * xcmfs * zlm(:,jlayer-1) * sqrt(ptke(:,jlayer-1)) 
END DO
!
!
!*    3.2 mixing coefficient vertical derivative at mid layers (at half levels)
!         --------------------------------------
!
!* there is no formal dependency on wind
zdkddvdz = 0.
!
!
!*    3.3  time evolution of wind in canopy
!          --------------------------------
!
 CALL canopy_evol_wind(ki,klvl,ptstep,kimpl,pwind,zk,zdkddvdz,psflux_u,pforc_u,pdforc_udu,pdz,pdzf,pu,zuw,palfau,pbetau)
!
!*    3.4  Friction velocity at top of SBL layers
!          --------------------------------------
!
pustar = sqrt(abs(zuw(:,klvl)))
!
!-------------------------------------------------------------------------------
!
!*    4. time evolution of TKE in canopy
!        -------------------------------
!
!* neutral case, no thermal production
zwth = 0.
zwq  = 0.
zth  = 300.
!
 CALL canopy_evol_tke(ki,klvl,ptstep,prhoa,pz,pzf,pdz,pdzf,pforc_e,pdforc_ede, &
                     pu,zth,zuw,zwth,zwq,zleps,ptke                           )
!
!-------------------------------------------------------------------------------
!
!*    5. Security at atmospheric forcing level
!        -------------------------------------
!
pu(:,klvl) = pwind(:)
IF (lhook) CALL dr_hook('CANOPY_EVOL_NEUTRAL',1,zhook_handle)
!
!-------------------------------------------------------------------------------
END SUBROUTINE canopy_evol_neutral