canopy_evol_tke.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_tke(SB, KI, PTSTEP, PRHOA, PFORC_E, PDFORC_EDE, &
                                 PTH, PUW, PWTH, PWQ, PLEPS )  
!     #########################################
!
!!****  *CANOPY_EVOL_TKE* - evolution of wind in canopy 
!!                        
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    07/2006 
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_canopy_n, ONLY : canopy_t
!
USE modi_tridiag_surf
USE modd_canopy_turb, ONLY : xced, xtkemin
USE modd_csts,        ONLY : xg, xrd, xrv
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
TYPE(canopy_t), INTENT(INOUT) :: SB
!
INTEGER,                  INTENT(IN)    :: KI        ! number of horizontal points
REAL,                     INTENT(IN)    :: PTSTEP    ! time-step                             (s)
REAL, DIMENSION(KI),      INTENT(IN)    :: PRHOA     ! Air density                           (kg/m3)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PFORC_E   ! tendency of wind due to canopy drag   (m2/s3)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PDFORC_EDE! formal derivative of the tendency of
!                                                    ! wind due to canopy drag               (1/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PTH       ! pot. temp. at canopy levels           (K)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PUW       ! turbulent flux (at half levels)       (m2/s2)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PWTH      ! w'Th' flux (at half levels)           (mK/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PWQ       ! w'q'  flux (at half levels)           (kg/m2/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PLEPS     ! dissipative length (full levels)      (m)
!
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER                    :: JLAYER   ! loop counter on layers
!
REAL, DIMENSION(KI,SB%NLVL)   :: ZDUDZ    ! dU/dz at mid levels
REAL, DIMENSION(KI,SB%NLVL)   :: ZDP      ! dynamical production at full levels
!                                      ! (at full levels)
REAL, DIMENSION(KI,SB%NLVL)   :: ZTP      ! thermal   production at full levels
!                                      ! (at full levels)
REAL, DIMENSION(KI,SB%NLVL)   :: ZDISS_O_TKE ! dissipation/TKE ratio at full levels
REAL, DIMENSION(KI,SB%NLVL)   :: ZF       ! turbulent flux at mid levels
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFDDVDZ ! derivative of turbulent flux as a
!                                      ! function of vertical gradient of wind variable
!                                      ! (at mid levels)
REAL, DIMENSION(KI,SB%NLVL)   :: ZEXT     ! external forcing at full levels
REAL, DIMENSION(KI,SB%NLVL)   :: ZDEXTDV  ! derivative of external forcing as a
!                                      ! function of vertical variable
!                                      ! (at full levels)
REAL, DIMENSION(KI,SB%NLVL)   :: ZTKE     ! TKE     at canopy levels (work var.)
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!
!*    1. initializations
!        ---------------
!
!* external forces
!
IF (lhook) CALL dr_hook('CANOPY_EVOL_TKE',0,zhook_handle)
zext = 0.
zdextdv = 0.
!
!* turbulent flux
!
zf       = 0.
zdfddvdz = 0.
!
!-------------------------------------------------------------------------------
!
!*    2. wind vertical derivative mid layers (at half levels below full levels)
!        ------------------------
!
zdudz(:,:) = -999.
DO jlayer=2,sb%NLVL
  zdudz(:,jlayer) = (sb%XU(:,jlayer) - sb%XU(:,jlayer-1)) / sb%XDZF(:,jlayer)
END DO
!-------------------------------------------------------------------------------
!
!*    3. Dynamical production of TKE (at full levels)
!        ---------------------------
!
zdp(:,:) = -999.
!
!* first level using an extrapolation using a 1/z law
zdp(:,1) = - puw(:,2) * zdudz(:,2) * (sb%XZF(:,2)/sb%XZ(:,1))

! other levels
DO jlayer=2,sb%NLVL-1
  zdp(:,jlayer) = - 0.5 * (puw(:,jlayer)   * zdudz(:,jlayer)  ) &
                    - 0.5 * (puw(:,jlayer+1) * zdudz(:,jlayer+1))  
END DO
!
!* upper level using an extrapolation using a 1/z law
zdp(:,sb%NLVL) = - puw(:,sb%NLVL) * zdudz(:,sb%NLVL) * (sb%XZF(:,sb%NLVL)/sb%XZ(:,sb%NLVL))
!
!* adds dynamical production in non-transport forces
zext    = zext    + zdp
zdextdv = zdextdv + 0.
!
!-------------------------------------------------------------------------------
!
!*    4. Thermal production of TKE (at full levels)
!        -------------------------
!
ztp(:,:) = -999.
!
! other levels
DO jlayer=1,sb%NLVL-1
  ztp(:,jlayer) = xg/pth(:,jlayer) * 0.5 * ( (pwth(:,jlayer) + pwth(:,jlayer+1) )          &
                             + (xrv/xrd-1) * (pwq(:,jlayer) + pwq(:,jlayer)   )/prhoa(:) )  
END DO
!
!* upper level
ztp(:,sb%NLVL) = xg/pth(:,sb%NLVL) * pwth(:,sb%NLVL)
!
!* adds dynamical production in non-transport forces
zext    = zext    + ztp
zdextdv = zdextdv + 0.
!
!-------------------------------------------------------------------------------
!
!*    4. Dissipation/TKE ratio (to prepare implicitation of dissipation)
!        ---------------------
!
zdiss_o_tke = - xced * sqrt(sb%XTKE(:,:)) / pleps(:,:)
!
!* adds dissipation in non-transport forces
zext    = zext    +       zdiss_o_tke * sb%XTKE(:,:)
zdextdv = zdextdv + 1.5 * zdiss_o_tke
!
!-------------------------------------------------------------------------------
!
!*    6. Adds Creation force due to canopy
!        ---------------------------------
!
!
zext    = zext    + pforc_e(:,:)
zdextdv = zdextdv + pdforc_ede(:,:)


!-------------------------------------------------------------------------------
!
!*    7. Evolution of TKE due to Dyn prod, Dissipation and Drag production
!        -----------------------------------------------------------------
!
!* note that dissipation is implicited
!
 CALL tridiag_surf(sb%XTKE,zf,zdfddvdz,zext,zdextdv,ptstep,sb%XDZF,sb%XDZ,ztke)  
!
!-------------------------------------------------------------------------------
!
!*    7. New value of TKE (at full levels)
!        ----------------
!
sb%XTKE(:,:) = ztke(:,:)
!
!-------------------------------------------------------------------------------
!
!*    8. Security at all levels : set minimum threshold
!        ----------------------
!
sb%XTKE(:,:) = max(sb%XTKE,xtkemin)
!
IF (lhook) CALL dr_hook('CANOPY_EVOL_TKE',1,zhook_handle)
!
!----------------------------------------------------------------
!
END SUBROUTINE canopy_evol_tke