canopy_evol.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(KI,KLVL,PTSTEP, KIMPL,                                     &
                              pzz,pwind,pta,pqa,ppa,prhoa,                              &
                              psflux_u,psflux_t,psflux_q,                               &
                              pforc_u,pdforc_udu,pforc_e,pdforc_ede,pforc_t,pdforc_tdt, &
                              pforc_q,pdforc_qdq,                                       &
                              pz,pzf,pdz,pdzf,pu,ptke,pt,pq,plmo,plm,pleps,pp,pustar,   &
                              palfau,pbetau,palfath,pbetath,palfaq,pbetaq, oneutral     ) 
!     #########################################
!
!!****  *CANOPY_EVOL* - evolution of canopy
!!                        
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    07/2006 
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_csts,        ONLY : xg, xrd, xcpd, xp00
USE modd_surf_par,    ONLY : xundef
USE modd_canopy_turb, ONLY : xcmfs, xcshf
!
USE modi_rmc01_surf
USE modi_canopy_evol_wind
USE modi_canopy_evol_tke
USE modi_canopy_evol_temp
!
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,KLVL), INTENT(IN)    :: pzz       ! Mixing length generic profile at mid levels (-)
REAL, DIMENSION(KI),      INTENT(IN)    :: pwind     ! wind speed                            (m/s)
REAL, DIMENSION(KI),      INTENT(IN)    :: pta       ! Air temperature                       (K)
REAL, DIMENSION(KI),      INTENT(IN)    :: pqa       ! Air humidity                          (kg/m3)
REAL, DIMENSION(KI),      INTENT(IN)    :: ppa       ! Pressure at forcing level             (Pa)
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),      INTENT(IN)    :: psflux_t  ! surface flux w'T'                     (Km/s)
REAL, DIMENSION(KI),      INTENT(IN)    :: psflux_q  ! surface flux w'q'                     (kg/m2/s)
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(IN)    :: pforc_t   ! tendency of Temp due to canopy drag   (T/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdforc_tdt! formal derivative of the tendency of
!                                                    ! Temp due to canopy drag               (1/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pforc_q   ! tendency of Hum. due to canopy drag   (kg/m3/s)
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: pdforc_qdq! formal derivative of the tendency of
!                                                    ! Hum. due to canopy drag               (1/s)
!
REAL, DIMENSION(KI,KLVL), INTENT(IN)    :: 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,KLVL), INTENT(INOUT) :: pt        ! Temperature at canopy levels          (K)
REAL, DIMENSION(KI,KLVL), INTENT(INOUT) :: pq        ! Humidity    at canopy levels          (kg/m3)
REAL, DIMENSION(KI,KLVL), INTENT(INOUT) :: plmo      ! Monin-Obhukov length                  (m)
REAL, DIMENSION(KI,KLVL), INTENT(OUT)   :: plm       ! mixing length                         (m)
REAL, DIMENSION(KI,KLVL), INTENT(OUT)   :: pleps     ! dissipative length                    (m)
REAL, DIMENSION(KI,KLVL), INTENT(INOUT) :: pp        ! Pressure    at canopy levels          (Pa)
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
REAL, DIMENSION(KI),      INTENT(OUT)   :: palfath  ! Th+(1) = alfa w'th'(1) + beta
REAL, DIMENSION(KI),      INTENT(OUT)   :: pbetath  ! Th+(1) = alfa w'th'(1) + beta
REAL, DIMENSION(KI),      INTENT(OUT)   :: palfaq   ! Q+(1) = alfa w'q'(1) + beta
REAL, DIMENSION(KI),      INTENT(OUT)   :: pbetaq   ! Q+(1) = alfa w'q'(1) + beta
!
LOGICAL, OPTIONAL, INTENT(IN) :: oneutral
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
LOGICAL :: gneutral
!
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) :: zth                ! potential temperature at full levels
REAL, DIMENSION(KI)      :: ztha               ! potential temperature at forcing level
REAL, DIMENSION(KI,KLVL) :: zexn               ! Exner function        at full levels
REAL, DIMENSION(KI,KLVL) :: zuw                ! friction at mid levels
REAL, DIMENSION(KI)      :: zsflux_th          ! Surface flux w'Th'                    (mK/s)
REAL, DIMENSION(KI,KLVL) :: zforc_th           ! tendency of Temp due to canopy drag   (T/s)
REAL, DIMENSION(KI,KLVL) :: zdforc_thdth       ! formal derivative of the tendency of
!                                              ! Temp due to canopy drag               (1/s)
REAL, DIMENSION(KI,KLVL) :: zwth               ! w'Th' at mid levels
REAL, DIMENSION(KI,KLVL) :: zwq                ! w'q'  at mid levels
REAL, DIMENSION(KI,KLVL) :: zsfth              ! heat flux at atmospheric forcing level
REAL, DIMENSION(KI,KLVL) :: zsfrv              ! vapor flux at atmospheric forcing level
REAL                     :: zz0                ! a value of z0 just for first time-step init.
REAL, DIMENSION(KI,KLVL) :: zrhoa              ! air density profile
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('CANOPY_EVOL',0,zhook_handle)
!
gneutral = .false. 
!
IF (present(oneutral)) gneutral = oneutral
!
!*    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
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
!-------------------------------------------------------------------------------
!
!*    5. mixing and dissipative lengths (at full levels)
!        ------------------------------
!
 CALL rmc01_surf(pzz,plmo,plm,pleps,gneutral)
!
!-------------------------------------------------------------------------------
!
!*    6. time evolution of wind in canopy
!        --------------------------------
!
!*    6.1 mixing coefficient for wind at mid layers (at half levels)
!         ---------------------------
!
zk = -999.
DO jlayer=2,klvl
  zk(:,jlayer) = 0.5 * xcmfs * plm(:,jlayer)   * sqrt(ptke(:,jlayer)  ) &
               + 0.5 * xcmfs * plm(:,jlayer-1) * sqrt(ptke(:,jlayer-1))  
END DO
!
!
!*    6.2 mixing coefficient vertical derivative at mid layers (at half levels)
!         --------------------------------------
!
!* there is no formal dependency on wind
zdkddvdz = 0.
!
!
!*    6.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)
!
!*    6.4  Friction velocity at top of SBL layers
!          --------------------------------------
!
pustar = sqrt(abs(zuw(:,klvl)))
!
!-------------------------------------------------------------------------------
!
IF (gneutral) THEN
  !
  zth  = 300.  
  zwth = 0.
  zwq  = 0.
  !
ELSE
  !
  !*    7. time evolution of temperature in canopy
  !        ---------------------------------------
  !
  !*    7.3 convertion into potential temperature (at half levels)
  !         -------------------------------------
  !
  DO jlayer=1,klvl
    pp(:,jlayer) = ppa(:) + xg * prhoa(:) * (pz(:,klvl) - pz(:,jlayer))
  END DO
  zexn = (pp/xp00)**(xrd/xcpd)
  !
  zth  = xundef
  WHERE(pt/=xundef) zth(:,:) = pt(:,:) / zexn(:,:)
  !
  ztha(:) = pta(:) / zexn(:,klvl)
  !
  !*    7.1 mixing coefficient for wind at mid layers (at half levels)
  !         ---------------------------
  !
  zk = -999.
  DO jlayer=2,klvl
    zk(:,jlayer) = 0.5 * xcshf * plm(:,jlayer)   * sqrt(ptke(:,jlayer)  ) &
                 + 0.5 * xcshf * plm(:,jlayer-1) * sqrt(ptke(:,jlayer-1))  
  END DO
  !
  !*    7.2 mixing coefficient vertical derivative at mid layers (at half levels)
  !         --------------------------------------
  !
  !* there is no formal dependency on temperature
  zdkddvdz = 0.
  !
  !*    7.4  conversion of canopy tendency into potential temperature tendency
  !          -----------------------------------------------------------------
  !
  zsflux_th    = psflux_t / zexn(:,1)
  zforc_th     = pforc_t  / zexn
  zdforc_thdth = pdforc_tdt
  !
  !
  !*    7.5  time evolution of temperature in canopy
  !          ---------------------------------------
  !
  CALL canopy_evol_temp(ki,klvl,ptstep,kimpl,ztha,zk,zdkddvdz,zsflux_th,zforc_th,zdforc_thdth,pdz,pdzf,zth,zwth,palfath,pbetath)
  !
  !*    7.6  convertion into absolute temperature
  !          ------------------------------------
  !
  WHERE(pt/=xundef) pt(:,:) = zth(:,:) * zexn(:,:)
  !
  !-------------------------------------------------------------------------------
  !
  !*    8. time evolution of Humidity in canopy
  !        ------------------------------------
  !
  CALL canopy_evol_temp(ki,klvl,ptstep,kimpl,pqa,zk,zdkddvdz,psflux_q,pforc_q,pdforc_qdq,pdz,pdzf,pq,zwq,palfaq,pbetaq)
  !
  !-------------------------------------------------------------------------------
  IF (kimpl==1 .AND. lhook) CALL dr_hook('CANOPY_EVOL',1,zhook_handle)
  IF (kimpl==1) RETURN
  !-------------------------------------------------------------------------------
  !
ENDIF
!
!*    9. time evolution of TKE in canopy
!        -------------------------------
!
 CALL canopy_evol_tke(ki,klvl,ptstep,prhoa,pz,pzf,pdz,pdzf,pforc_e,pdforc_ede, &
                      pu,zth,zuw,zwth,zwq,pleps,ptke                          ) 
!
!-------------------------------------------------------------------------------
!
IF (.NOT.gneutral) THEN
  !
  !*   10. Monin-Obuhkov length
  !        --------------------
  !
  !* MO length is estimated using the heat and vapor turbulent fluxes at atmospheric level
  !  (it avoids the problems of vertical variation of the fluxes in the canopy)
  !  However, friction flux MUST be taken as the maximum flux on the
  !  profile, in order to avoid unrealistically small MO length when using
  !  small time-steps
  !
  zrhoa(:,:) = spread(prhoa(:),2,klvl)
  !
  zsfth(:,:)  = zwth(:,:)
  zsfrv(:,:)  = zwq(:,:) / zrhoa(:,:)
  !
  plmo(:,:) = lmo(sqrt(abs(zuw)),pt,pq,zsfth,zsfrv)
  !
  DO jlayer=1,klvl
    WHERE (plmo(:,jlayer)>0.) plmo(:,jlayer) = max(plmo(:,jlayer),pz(:,klvl))
    WHERE (plmo(:,jlayer)<0.) plmo(:,jlayer) = min(plmo(:,jlayer),-pz(:,klvl))     
  ENDDO
  !
  !-------------------------------------------------------------------------------
  !
  !*   11. Security at atmospheric forcing level
  !        -------------------------------------
  !
  pt(:,klvl) = pta(:)
  !
  pq(:,klvl) = pqa(:)
  !
ENDIF
!
pu(:,klvl) = pwind(:)
!
IF (lhook) CALL dr_hook('CANOPY_EVOL',1,zhook_handle)
!
!-------------------------------------------------------------------------------
END SUBROUTINE canopy_evol