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(SB, KI, 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, PLM, PLEPS, 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_canopy_n, ONLY : canopy_t
!
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
!              -------------------------
!
TYPE(canopy_t), INTENT(INOUT) :: SB
!
INTEGER,                  INTENT(IN)    :: KI        ! number of horizontal points
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,SB%NLVL), 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,SB%NLVL), INTENT(IN)    :: PFORC_U   ! tendency of wind due to canopy drag   (m/s2)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PDFORC_UDU! formal derivative of the tendency of
!                                                    ! wind due to canopy drag               (1/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PFORC_E   ! tendency of TKE  due to canopy drag   (m2/s3)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PDFORC_EDE! formal derivative of the tendency of
!                                                    ! TKE  due to canopy drag               (1/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PFORC_T   ! tendency of Temp due to canopy drag   (T/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PDFORC_TDT! formal derivative of the tendency of
!                                                    ! Temp due to canopy drag               (1/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PFORC_Q   ! tendency of Hum. due to canopy drag   (kg/m3/s)
REAL, DIMENSION(KI,SB%NLVL), INTENT(IN)    :: PDFORC_QDQ! formal derivative of the tendency of
!                                                    ! Hum. due to canopy drag               (1/s)
!
REAL, DIMENSION(KI,SB%NLVL), INTENT(OUT)   :: PLM       ! mixing length                         (m)
REAL, DIMENSION(KI,SB%NLVL), INTENT(OUT)   :: PLEPS     ! dissipative length                    (m)
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,SB%NLVL) :: ZK                 ! mixing coefficient
REAL, DIMENSION(KI,SB%NLVL) :: ZDKDDVDZ           ! formal derivative of mixing coefficient according to variable vertical gradient
REAL, DIMENSION(KI,SB%NLVL) :: ZTH                ! potential temperature at full levels
REAL, DIMENSION(KI)         :: ZTHA               ! potential temperature at forcing level
REAL, DIMENSION(KI,SB%NLVL) :: ZEXN               ! Exner function        at full levels
REAL, DIMENSION(KI,SB%NLVL) :: ZUW                ! friction at mid levels
REAL, DIMENSION(KI)         :: ZSFLUX_TH          ! Surface flux w'Th'                    (mK/s)
REAL, DIMENSION(KI,SB%NLVL) :: ZFORC_TH           ! tendency of Temp due to canopy drag   (T/s)
REAL, DIMENSION(KI,SB%NLVL) :: ZDFORC_THDTH       ! formal derivative of the tendency of
!                                              ! Temp due to canopy drag               (1/s)
REAL, DIMENSION(KI,SB%NLVL) :: ZWTH               ! w'Th' at mid levels
REAL, DIMENSION(KI,SB%NLVL) :: ZWQ                ! w'q'  at mid levels
REAL, DIMENSION(KI,SB%NLVL) :: ZSFTH              ! heat flux at atmospheric forcing level
REAL, DIMENSION(KI,SB%NLVL) :: ZSFRV              ! vapor flux at atmospheric forcing level
REAL                        :: ZZ0                ! a value of z0 just for first time-step init.
REAL, DIMENSION(KI,SB%NLVL) :: 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. sb%XU(ji,sb%NLVL-1)==0.) THEN
    zz0 = sb%XZ(ji,1)/2.
    sb%XU(ji,:) = pwind(ji) * log(sb%XZ(ji,:)/zz0) / log(sb%XZ(ji,sb%NLVL)/zz0)
  END IF
END DO
!-------------------------------------------------------------------------------
!
!*    5. mixing and dissipative lengths (at full levels)
!        ------------------------------
!
 CALL rmc01_surf(pzz,sb%XLMO,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,sb%NLVL
  zk(:,jlayer) = 0.5 * xcmfs * plm(:,jlayer)   * sqrt(sb%XTKE(:,jlayer)  ) &
               + 0.5 * xcmfs * plm(:,jlayer-1) * sqrt(sb%XTKE(:,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(sb, ki, ptstep, kimpl, pwind, zk, zdkddvdz,&
                       psflux_u, pforc_u, pdforc_udu, zuw, palfau, pbetau)
!
!*    6.4  Friction velocity at top of SBL layers
!          --------------------------------------
!
pustar = sqrt(abs(zuw(:,sb%NLVL)))
!
!-------------------------------------------------------------------------------
!
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,sb%NLVL
    sb%XP(:,jlayer) = ppa(:) + xg * prhoa(:) * (sb%XZ(:,sb%NLVL) - sb%XZ(:,jlayer))
  END DO
  zexn = (sb%XP/xp00)**(xrd/xcpd)
  !
  zth  = xundef
  WHERE(sb%XT/=xundef) zth(:,:) = sb%XT(:,:) / zexn(:,:)
  !
  ztha(:) = pta(:) / zexn(:,sb%NLVL)
  !
  !*    7.1 mixing coefficient for wind at mid layers (at half levels)
  !         ---------------------------
  !
  zk = -999.
  DO jlayer=2,sb%NLVL
    zk(:,jlayer) = 0.5 * xcshf * plm(:,jlayer)   * sqrt(sb%XTKE(:,jlayer)  ) &
                 + 0.5 * xcshf * plm(:,jlayer-1) * sqrt(sb%XTKE(:,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(sb, ki, ptstep, kimpl, ztha, zk, zdkddvdz, &
                        zsflux_th, zforc_th, zdforc_thdth, zth, zwth, palfath, pbetath)
  !
  !*    7.6  convertion into absolute temperature
  !          ------------------------------------
  !
  WHERE(sb%XT/=xundef) sb%XT(:,:) = zth(:,:) * zexn(:,:)
  !
  !-------------------------------------------------------------------------------
  !
  !*    8. time evolution of Humidity in canopy
  !        ------------------------------------
  !
  CALL canopy_evol_temp(sb, ki, ptstep, kimpl, pqa, zk, zdkddvdz, &
                        psflux_q, pforc_q, pdforc_qdq, sb%XQ, 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(sb, ki, ptstep, prhoa, pforc_e, pdforc_ede, zth, zuw, zwth, zwq, pleps) 
!
!-------------------------------------------------------------------------------
!
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,sb%NLVL)
  !
  zsfth(:,:)  = zwth(:,:)
  zsfrv(:,:)  = zwq(:,:) / zrhoa(:,:)
  !
  sb%XLMO(:,:) = lmo(sqrt(abs(zuw)),sb%XT,sb%XQ,zsfth,zsfrv)
  !
  DO jlayer=1,sb%NLVL
    WHERE (sb%XLMO(:,jlayer)>0.) sb%XLMO(:,jlayer) = max(sb%XLMO(:,jlayer),sb%XZ(:,sb%NLVL))
    WHERE (sb%XLMO(:,jlayer)<0.) sb%XLMO(:,jlayer) = min(sb%XLMO(:,jlayer),-sb%XZ(:,sb%NLVL))     
  ENDDO
  !
  !-------------------------------------------------------------------------------
  !
  !*   11. Security at atmospheric forcing level
  !        -------------------------------------
  !
  sb%XT(:,sb%NLVL) = pta(:)
  !
  sb%XQ(:,sb%NLVL) = pqa(:)
  !
ENDIF
!
sb%XU(:,sb%NLVL) = pwind(:)
!
IF (lhook) CALL dr_hook('CANOPY_EVOL',1,zhook_handle)
!
!-------------------------------------------------------------------------------
END SUBROUTINE canopy_evol