coupling_isba_canopyn.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 coupling_isba_canopy_n (DTCO, UG, U, USS, SB, NAG, CHI, NCHI, DTV, ID, NGB, GB, &
                                   ISS, NISS, IG, NIG, IO, S, K, NK, NP, NPE, NDST, SLT,   &
                                   HPROGRAM, HCOUPLING, PTSTEP,                            &
                                   KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN,        &
                                   PZENITH, PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV,    &
                                   PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW,     &
                                   PDIR_SW, PSCA_SW, PSW_BANDS, PPS, PPA, PSFTQ, PSFTH,    &
                                   PSFTS, PSFCO2, PSFU, PSFV, PTRAD, PDIR_ALB, PSCA_ALB,   &
                                   PEMIS, PTSURF, PZ0,PZ0H, PQSURF, PPEW_A_COEF,           &
                                   PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF,     &
                                   PPEQ_B_COEF, HTEST  )
!     ###############################################################################
!
!!****  *COUPLING_ISBA_CANOPY_n * - Adds a SBL into ISBA
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    09/2007
!!      S. Riette   06/2009 Initialisation of XT, XQ, XU and XTKE on canopy levels
!!      S. Riette   01/2010 Use of interpol_sbl to compute 10m wind diagnostic
!!      Modified    09/2012  : J. Escobar , SIZE(PTA) not allowed without-interface , replace by KI
!!      B. Decharme  04/2013 new coupling variables
!----------------------------------------------------------------
!
USE modd_agri_n, ONLY : agri_np_t
USE modd_ch_isba_n, ONLY : ch_isba_t, ch_isba_np_t
USE modd_data_isba_n, ONLY : data_isba_t
USE modd_surfex_n, ONLY : isba_diag_t
USE modd_gr_biog_n, ONLY : gr_biog_t, gr_biog_np_t
USE modd_sso_n, ONLY : sso_t, sso_np_t
USE modd_sfx_grid_n, ONLY : grid_t, grid_np_t
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_t, isba_p_t, isba_pe_t, isba_nk_t, isba_np_t, isba_npe_t
!
USE modd_dst_n, ONLY : dst_np_t

USE modd_canopy_n, ONLY : canopy_t
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_sso_n, ONLY : sso_t
USE modd_data_isba_n, ONLY : data_isba_t
USE modd_slt_n, ONLY : slt_t
!
!
USE modd_csts,          ONLY : xcpd
USE modd_surf_par,      ONLY : xundef
USE modd_canopy_turb,   ONLY : xalpsbl
!
USE mode_coupling_canopy
!
USE modi_init_isba_sbl
!
USE modi_canopy_evol
USE modi_canopy_grid_update
!
USE modi_coupling_isba_n
!
USE modi_isba_canopy
USE modi_sso_beljaars04
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(agri_np_t), INTENT(INOUT) :: NAG
TYPE(ch_isba_t), INTENT(INOUT) :: CHI
TYPE(ch_isba_np_t), INTENT(INOUT) :: NCHI
TYPE(data_isba_t), INTENT(INOUT) :: DTV
TYPE(isba_diag_t), INTENT(INOUT) :: ID
TYPE(gr_biog_np_t), INTENT(INOUT) :: NGB
TYPE(gr_biog_t), INTENT(INOUT) :: GB
TYPE(sso_t), INTENT(INOUT) :: ISS 
TYPE(sso_np_t), INTENT(INOUT) :: NISS
TYPE(grid_t), INTENT(INOUT) :: IG
TYPE(grid_np_t), INTENT(INOUT) :: NIG
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_nk_t), INTENT(INOUT) :: NK
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) ::NPE
!
TYPE(dst_np_t), INTENT(INOUT) :: NDST
!
TYPE(canopy_t), INTENT(INOUT) :: SB
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
TYPE(slt_t), INTENT(INOUT) :: SLT
!
 CHARACTER(LEN=6),    INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=1),    INTENT(IN)  :: HCOUPLING ! type of coupling
                                              ! 'E' : explicit
                                              ! 'I' : implicit
INTEGER,             INTENT(IN)  :: KYEAR     ! current year (UTC)
INTEGER,             INTENT(IN)  :: KMONTH    ! current month (UTC)
INTEGER,             INTENT(IN)  :: KDAY      ! current day (UTC)
REAL,                INTENT(IN)  :: PTIME     ! current time since midnight (UTC, s)
INTEGER,             INTENT(IN)  :: KI        ! number of points
INTEGER,             INTENT(IN)  :: KSV       ! number of scalars
INTEGER,             INTENT(IN)  :: KSW       ! number of short-wave spectral bands
REAL, DIMENSION(KI), INTENT(IN)  :: PTSUN     ! solar time                    (s from midnight)
REAL,                INTENT(IN)  :: PTSTEP    ! atmospheric time-step                 (s)
REAL, DIMENSION(KI), INTENT(IN)  :: PZREF     ! height of T,q forcing                 (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PUREF     ! height of wind forcing                (m)
!
REAL, DIMENSION(KI), INTENT(IN)  :: PTA       ! air temperature forcing               (K)
REAL, DIMENSION(KI), INTENT(IN)  :: PQA       ! air humidity forcing                  (kg/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PRHOA     ! air density                           (kg/m3)
REAL, DIMENSION(KI,KSV),INTENT(IN) :: PSV     ! scalar variables
!                                             ! chemistry:   first char. in HSV: '#'  (molecule/m3)
!                                             !
 CHARACTER(LEN=6), DIMENSION(KSV),INTENT(IN):: HSV  ! name of all scalar variables
REAL, DIMENSION(KI), INTENT(IN)  :: PU        ! zonal wind                            (m/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PV        ! meridian wind                         (m/s)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PDIR_SW ! direct  solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PSCA_SW ! diffuse solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KSW),INTENT(IN)  :: PSW_BANDS ! mean wavelength of each shortwave band (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH   ! zenithal angle at t       (radian from the vertical)
REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH2  ! zenithal angle at t+1     (radian from the vertical)
REAL, DIMENSION(KI), INTENT(IN)  :: PAZIM     ! azimuthal angle      (radian from North, clockwise)
REAL, DIMENSION(KI), INTENT(IN)  :: PLW       ! longwave radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI), INTENT(IN)  :: PPS       ! pressure at atmospheric model surface (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PPA       ! pressure at forcing level             (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PZS       ! atmospheric model CANOPY           (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PCO2      ! CO2 concentration in the air          (kg/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PSNOW     ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PRAIN     ! liquid precipitation                  (kg/m2/s)
!
!
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTH     ! flux of heat                          (W/m2)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTQ     ! flux of water vapor                   (kg/m2/s)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFU      ! zonal momentum flux                   (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFV      ! meridian momentum flux                (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFCO2    ! flux of CO2                           (m/s*kg_CO2/kg_air)
REAL, DIMENSION(KI,KSV),INTENT(OUT):: PSFTS   ! flux of scalar var.                   (kg/m2/s)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTRAD     ! radiative temperature                 (K)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PDIR_ALB! direct albedo for each spectral band  (-)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PSCA_ALB! diffuse albedo for each spectral band (-)
REAL, DIMENSION(KI), INTENT(OUT) :: PEMIS     ! emissivity                            (-)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTSURF    ! surface effective temperature         (K)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0       ! roughness length for momentum         (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0H      ! roughness length for heat             (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PQSURF    ! specific humidity at surface          (kg/kg)
!
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_A_COEF! implicit coefficients
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_B_COEF! needed if HCOUPLING='I'
REAL, DIMENSION(KI), INTENT(IN) :: PPET_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPET_B_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_B_COEF
 CHARACTER(LEN=2),    INTENT(IN) :: HTEST ! must be equal to 'OK'
!
!*      0.2    declarations of local variables
!
!* forcing variables
!
TYPE(isba_p_t), POINTER :: PK
TYPE(isba_pe_t), POINTER :: PEK
!
REAL, DIMENSION(KI)     :: ZWIND    ! lowest atmospheric level wind speed           (m/s)
REAL, DIMENSION(KI)     :: ZEXNA    ! Exner function at lowest SBL scheme level     (-)
REAL, DIMENSION(KI)     :: ZTA      ! temperature                                   (K)
REAL, DIMENSION(KI)     :: ZPA      ! pressure                                      (Pa)
REAL, DIMENSION(KI)     :: ZZREF    ! temperature forcing level                     (m)
REAL, DIMENSION(KI)     :: ZUREF    ! wind        forcing level                     (m)
REAL, DIMENSION(KI)     :: ZU       ! zonal wind                                    (m/s)
REAL, DIMENSION(KI)     :: ZV       ! meridian wind                                 (m/s)
REAL, DIMENSION(KI)     :: ZQA      ! specific humidity                             (kg/m3)
REAL, DIMENSION(KI)     :: ZPEQ_A_COEF ! specific humidity implicit
REAL, DIMENSION(KI)     :: ZPEQ_B_COEF ! coefficients (hum. in kg/kg)
!
!
! canopy turbulence scheme
!
REAL, DIMENSION(KI)        :: ZCANOPY   ! height of canopy   (m)
REAL, DIMENSION(KI)        :: ZSFLUX_U  ! Surface flux u'w' (m2/s2)
REAL, DIMENSION(KI)        :: ZSFLUX_T  ! Surface flux w'T' (mK/s)
REAL, DIMENSION(KI)        :: ZSFLUX_Q  ! Surface flux w'q' (kgm2/s)
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_U   ! tendency due to drag force for wind
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_UDU! formal derivative of
!                                              ! tendency due to drag force for wind
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_E   ! tendency due to drag force for TKE
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_EDE! formal derivative of
!                                              ! tendency due to drag force for TKE
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_T   ! tendency due to drag force for Temp
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_TDT! formal derivative of
!                                              ! tendency due to drag force for Temp
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_Q   ! tendency due to drag force for Temp
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_QDQ! formal derivative of
!                                              ! tendency due to drag force for hum.
REAL, DIMENSION(KI,SB%NLVL)   :: ZLM       ! mixing length
REAL, DIMENSION(KI,SB%NLVL)   :: ZLEPS     ! dissipative length
REAL, DIMENSION(KI)     :: ZH           ! canopy height (m)
REAL, DIMENSION(KI)     :: ZUSTAR       ! friction velocity including drag effect (m/s)
REAL, DIMENSION(KI)     :: ZUSTAR_GROUND! friction velocity at ground only (ISBA) (m/s)
!
REAL, DIMENSION(KI)     :: ZPET_A_COEF ! temperature implicit
REAL, DIMENSION(KI)     :: ZPET_B_COEF ! coefficients (K)
REAL, DIMENSION(KI)     :: ZPEW_A_COEF ! wind implicit
REAL, DIMENSION(KI)     :: ZPEW_B_COEF ! coefficients (m/s)
!
REAL, DIMENSION(KI)   :: ZALFAU   ! V+(1) = - alfa rho u'w'(1) + beta
REAL, DIMENSION(KI)   :: ZBETAU   ! V+(1) = - alfa rho u'w'(1) + beta
REAL, DIMENSION(KI)   :: ZALFATH  ! Th+(1) = - alfa rho w'th'(1) + beta
REAL, DIMENSION(KI)   :: ZBETATH  ! Th+(1) = - alfa rho w'th'(1) + beta
REAL, DIMENSION(KI)   :: ZALFAQ   ! Q+(1) = - alfa rho w'q'(1) + beta
REAL, DIMENSION(KI)   :: ZBETAQ   ! Q+(1) = - alfa rho w'q'(1) + beta
!
 CHARACTER(LEN=1) :: GCOUPLING
!
REAL, DIMENSION(KI)   ::ZCANOPY_DENSITY
REAL, DIMENSION(KI)   ::ZUW_GROUND
REAL, DIMENSION(KI)   ::ZDUWDU_GROUND
!
INTEGER                      :: JJ, JLAYER, IMASK, JP, JI
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------------------
!
!
!*      1.     Preliminary computations of the SBL scheme
!              ------------------------------------------
!
IF (lhook) CALL dr_hook('COUPLING_ISBA_CANOPY_N',0,zhook_handle)
!
IF (io%LCANOPY) THEN
!
!*      1.1    Updates canopy vertical grid as a function of forcing height
!              ------------------------------------------------------------
!
!* determines where is the forcing level and modifies the upper levels of the canopy grid
!
  zcanopy = 0.
  CALL canopy_grid_update(ki,zcanopy,puref,sb)
!
!
!
!*      1.2    Allocations and initialisations
!              -------------------------------
!
!
!       1.2.1  First time step canopy initialisation
!
  IF(any(sb%XT(:,:) == xundef)) THEN
    CALL init_isba_sbl(io, k, np, npe, sb, &
                       ptstep, ppa, pps, pta, pqa, prhoa, pu, pv, pdir_sw,  &
                       psca_sw, psw_bands, prain, psnow, pzref, puref, iss%XSSO_SLOPE )
  ENDIF
!
!*      1.3    Allocations
!              -----------
!
  CALL init_forc(zforc_u, zdforc_udu, zforc_e, zdforc_ede, zforc_t, zdforc_tdt, zforc_q, zdforc_qdq )
!
  zsflux_u = 0.
  zsflux_t = 0.
  zsflux_q = 0.
!
!* default :
!* no canopy in ISBA scheme
!
  zh = 0.
  zcanopy_density(:) = 0.
!
!
!
!*      1.4   canopy for wind drag only
!             -------------------------
!
  IF (io%LCANOPY_DRAG) THEN
!* mean canopy height
!
!* in ecoclimap, height is set retrieved from roughness length (z0/0.13)
    DO jp = 1,io%NPATCH
      pk => np%AL(jp)
      pek => npe%AL(jp)
      DO jj = 1,pk%NSIZE_P
        imask = pk%NR_P(jj)
        zh(imask) = zh(imask) + pk%XPATCH(jj) * pek%XZ0(jj) / 0.13
        zcanopy_density(imask) = zcanopy_density(imask) + pk%XPATCH(jj) *  pek%XLAI(jj)
      ENDDO
    ENDDO
    !
    DO jj=1,ki
      zh(jj) = min(zh(jj), sb%XZF(jj,sb%NLVL))
      IF (zh(jj)<=sb%XDZ(jj,1)) zh(jj) = 0.
!
!* canopy for wind drag only
      zuw_ground(jj)      = 0.
      zduwdu_ground(jj)   = 0.
      !
    ENDDO
!
!* computes tendencies on wind and Tke due to canopy
    CALL isba_canopy(io%XCDRAG, ki, sb, zh,zcanopy_density, zuw_ground, &
                     zduwdu_ground, zforc_u, zdforc_udu, zforc_e, zdforc_ede )
!
  ENDIF
!
!*      1.5   Computes coefficients for implicitation
!             ---------------------------------------
!
  zwind = sqrt(pu**2+pv**2)
  CALL canopy_evol(sb, ki, ptstep, 1, sb%XZ, zwind, pta, pqa, ppa, prhoa, &
                   zsflux_u, zsflux_t, zsflux_q, zforc_u, zdforc_udu,       &
                   zforc_e, zdforc_ede, zforc_t, zdforc_tdt, zforc_q,       &
                   zdforc_qdq, zlm, zleps, zustar, zalfau, zbetau, zalfath, &
                   zbetath, zalfaq, zbetaq   )
!
!*     1.6     Goes from atmospheric forcing to canopy forcing height
!              ------------------------------------------------------
!
  gcoupling ='I'
!
  CALL init_coupling_canopy(sb, ppa, pu, pv, prhoa, zalfau, &
                            zbetau, zalfath, zbetath, zalfaq,&
                            zbetaq, zpa, zta, zqa, zu, zv,   &
                            zuref, zzref, zexna, zpew_a_coef,&
                            zpew_b_coef, zpet_a_coef,        &
                            zpet_b_coef, zpeq_a_coef,        &
                            zpeq_b_coef          )
!
!-------------------------------------------------------------------------------------
ELSE
!-------------------------------------------------------------------------------------
!
!*      2.     If no canopy scheme is used, forcing is not modified
!              ----------------------------------------------------
!
  gcoupling = hcoupling
!
  CALL init_coupling(hcoupling, pps, ppa, pta, pqa, &
                     pu, pv, puref, pzref,          &
                     ppew_a_coef, ppew_b_coef,      &
                     ppet_a_coef, ppet_b_coef,      &
                     ppeq_a_coef, ppeq_b_coef,      &
                     zpa, zta, zqa, zu, zv, zuref,  &
                     zzref, zpew_a_coef,            &
                     zpew_b_coef, zpet_a_coef,      &
                     zpet_b_coef, zpeq_a_coef,      &
                     zpeq_b_coef    ) 
! 
END IF
!
!-------------------------------------------------------------------------------------
!
!*      2.     Call of ISBA
!              ------------
!
 CALL coupling_isba_n(dtco, ug, u, uss, nag, chi, nchi, dtv, id, ngb, gb, iss,niss, ig, &
                      nig, io, s, k, nk, np, npe, ndst, slt, hprogram, gcoupling,       &
                      ptstep, kyear, kmonth, kday, ptime, ki, ksv, ksw, ptsun, pzenith, &
                      pzenith2, zzref, zuref, pzs, zu, zv, zqa, zta, prhoa, psv, pco2,  &
                      hsv, prain, psnow, plw, pdir_sw, psca_sw, psw_bands, pps, zpa,    &
                      psftq, psfth, psfts, psfco2, psfu, psfv, ptrad, pdir_alb,         &
                      psca_alb, pemis, ptsurf, pz0, pz0h, pqsurf, zpew_a_coef,          &
                      zpew_b_coef, zpet_a_coef, zpeq_a_coef, zpet_b_coef, zpeq_b_coef,  &
                      'OK' )
!
!-------------------------------------------------------------------------------------
!
!*      3.     End if no canopy is used
!              ------------------------
!
IF (.NOT. io%LCANOPY .AND. lhook) CALL dr_hook('COUPLING_ISBA_CANOPY_N',1,zhook_handle)
IF (.NOT. io%LCANOPY) RETURN
!
!-------------------------------------------------------------------------------------
!
!*      4.     Computes the impact of surface on air
!              -------------------------------------
!
 CALL init_forc(zforc_u, zdforc_udu, zforc_e, zdforc_ede, &
                zforc_t, zdforc_tdt, zforc_q, zdforc_qdq )
!
zsflux_u = - sqrt(psfu(:)**2+psfv(:)**2) / prhoa(:)
zsflux_t(:) = psfth(:) / xcpd * zexna(:) / prhoa(:)
zsflux_q(:) = psftq(:)
!
!-------------------------------------------------------------------------------------
!
!*      5.     Computes the impact of canopy on air
!              ------------------------------------
!
IF (io%LCANOPY_DRAG) THEN
!
  DO jj=1,ki
    zuw_ground(jj) = -sqrt(psfu(jj)**2+psfv(jj)**2)/ prhoa(jj)
    zduwdu_ground(jj) = 0.
    IF (sb%XU(jj,1)/=0.) zduwdu_ground(jj) = 2. * zuw_ground(jj) / sb%XU(jj,1)
  ENDDO

!* computes tendencies on wind and Tke due to canopy and surface
  CALL isba_canopy(io%XCDRAG, ki, sb, zh, zcanopy_density, zuw_ground, &
                   zduwdu_ground, zforc_u, zdforc_udu, zforc_e, zdforc_ede  )

  zsflux_u = 0.  ! surface friction is incorporated in ZFORC_U by ISBA_CANOPY routine
!
END IF
!
!-------------------------------------------------------------------------------------
!
!*      6.    Evolution of canopy air due to these impacts
!             --------------------------------------------
!
zwind = sqrt(pu**2+pv**2)
 CALL canopy_evol(sb, ki, ptstep, 2, sb%XZ, zwind, pta, pqa, ppa, prhoa,  &
                  zsflux_u, zsflux_t, zsflux_q, zforc_u, zdforc_udu,zforc_e,&
                  zdforc_ede, zforc_t, zdforc_tdt, zforc_q, zdforc_qdq, zlm,&
                  zleps, zustar, zalfau, zbetau, zalfath, zbetath, zalfaq,  &
                  zbetaq                  )
!
DO jlayer=1,sb%NLVL-1
  sb%XLMO(:,jlayer) = sb%XLMO(:,sb%NLVL)
ENDDO
!
! Momentum fluxes if canopy is used
!
!* Total friction due to surface averaged friction and averaged canopy drag
IF (io%LCANOPY_DRAG) THEN
  DO jj=1,ki
    zustar_ground(jj) = sqrt(sqrt(psfu(jj)**2+psfv(jj)**2)/prhoa(jj))
    IF (zustar_ground(jj)>0.) THEN
      psfu(jj) = psfu(jj) * zustar(jj)**2/zustar_ground(jj)**2
      psfv(jj) = psfv(jj) * zustar(jj)**2/zustar_ground(jj)**2
    ENDIF
  ENDDO
!* Total friction due to surface averaged friction and averaged canopy drag
  IF (id%O%LSURF_BUDGET) THEN
    id%D%XFMU = psfu
    id%D%XFMV = psfv          
  ENDIF
END IF
!
!-------------------------------------------------------------------------------------
!
!*      7.    2m and 10m diagnostics if canopy is used
!             ----------------------------------------
!
IF (id%O%N2M>=1) CALL init_2m_10m(sb, id%D, pu, pv, zwind, prhoa )
!
IF (lhook) CALL dr_hook('COUPLING_ISBA_CANOPY_N',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
END SUBROUTINE coupling_isba_canopy_n