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(KI,KLVL,PTSTEP,PRHOA,             &
                                   pz,pzf,pdz,pdzf,                  &
                                   pforc_e,pdforc_ede,               &
                                   pu,pth,puw,pwth,pwq,pleps,ptke    )  
!     #########################################
!
!!****  *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 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
!              -------------------------
!
INTEGER,                  INTENT(IN)    :: ki        ! number of horizontal points
INTEGER,                  INTENT(IN)    :: klvl      ! number of levels in canopy
REAL,                     INTENT(IN)    :: ptstep    ! time-step                             (s)
REAL, DIMENSION(KI),      INTENT(IN)    :: prhoa     ! Air density                           (kg/m3)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pz        ! Z at full levels                      (m)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pzf       ! Z at half levels                      (m)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdz       ! deltaZ between canopy half levels,
!                                                    ! located at full levels                (m)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdzf      ! deltaZ between canopy (full) levels,
!                                                    ! located at half levels                (m)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pforc_e   ! tendency of wind due to canopy drag   (m2/s3)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdforc_ede! formal derivative of the tendency of
!                                                    ! wind due to canopy drag               (1/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pu        ! wind speed at canopy levels           (m/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pth       ! pot. temp. at canopy levels           (K)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: puw       ! turbulent flux (at half levels)       (m2/s2)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pwth      ! w'Th' flux (at half levels)           (mK/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pwq       ! w'q'  flux (at half levels)           (kg/m2/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pleps     ! dissipative length (full levels)      (m)
REAL, DIMENSION(KI,KLVL), INTENT(INOUT) :: ptke      ! TKE at canopy levels   
!
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER                    :: jlayer   ! loop counter on layers
!
REAL, DIMENSION(KI,KLVL)   :: zdudz    ! dU/dz at mid levels
REAL, DIMENSION(KI,KLVL)   :: zdp      ! dynamical production at full levels
!                                      ! (at full levels)
REAL, DIMENSION(KI,KLVL)   :: ztp      ! thermal   production at full levels
!                                      ! (at full levels)
REAL, DIMENSION(KI,KLVL)   :: zdiss_o_tke ! dissipation/TKE ratio at full levels
REAL, DIMENSION(KI,KLVL)   :: zf       ! turbulent flux at mid levels
REAL, DIMENSION(KI,KLVL)   :: zdfddvdz ! derivative of turbulent flux as a
!                                      ! function of vertical gradient of wind variable
!                                      ! (at mid levels)
REAL, DIMENSION(KI,KLVL)   :: zext     ! external forcing at full levels
REAL, DIMENSION(KI,KLVL)   :: zdextdv  ! derivative of external forcing as a
!                                      ! function of vertical variable
!                                      ! (at full levels)
REAL, DIMENSION(KI,KLVL)   :: 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,klvl
  zdudz(:,jlayer) = (pu(:,jlayer) - pu(:,jlayer-1)) / pdzf(:,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) * (pzf(:,2)/pz(:,1))

! other levels
DO jlayer=2,klvl-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(:,klvl) = - puw(:,klvl) * zdudz(:,klvl) * (pzf(:,klvl)/pz(:,klvl))
!
!* 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,klvl-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(:,klvl) = xg/pth(:,klvl) * pwth(:,klvl)
!
!* 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(ptke(:,:)) / pleps(:,:)
!
!* adds dissipation in non-transport forces
zext    = zext    +       zdiss_o_tke * ptke(:,:)
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(ptke,zf,zdfddvdz,zext,zdextdv,ptstep, &
                    pdzf,pdz,ztke                         )  
!
!-------------------------------------------------------------------------------
!
!*    7. New value of TKE (at full levels)
!        ----------------
!
ptke(:,:) = ztke(:,:)
!
!-------------------------------------------------------------------------------
!
!*    8. Security at all levels : set minimum threshold
!        ----------------------
!
ptke(:,:) = max(ptke,xtkemin)
IF (lhook) CALL dr_hook('CANOPY_EVOL_TKE',1,zhook_handle)
!
!----------------------------------------------------------------
!
END SUBROUTINE canopy_evol_tke