canopy_levs.F90

   module canopy_levs_mod
   contains
 
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
   subroutine canopy_levs_init(im, km, nkc, nkt,   &
              ntrac,  ntqv, ntke, & ! ndtend,          &
              zi, zl, zm,                          & ! in: 3D meters
              prsl, prsi,                          & ! in: 3D (Pa)
              dv, du, tdt, rtg,                    & ! in: 3D
              u1, v1, t1, q1,                      & ! in: 3D " 4D q1(ix,km,ntrac1) kg kg-1
              dens, dkt, dku,                      & ! in: 3D
!              dtend,                               & ! in: 4D
              zmom_can3, zmid_can3,                & !out: 3D
              sigmom_can, sigmid_can,              & !out
              zl_can, zm_can,                      & !out
              prsl_can, prsi_can,                  & !out
              dv_can, du_can, tdt_can, rtg_can,    & !out: 3D
              t1_can, qv_can, dens_can,            & !out
              ws_can, dkt_can, dku_can,            & !out
              q1_can, q1_2m,                       & !out
!              DTEND_CAN,                           &
              errmsg, errflg )                       !out
 
   use machine , only : kind_phys
 
   IMPLICIT NONE
 
!...Arguments:
   integer, intent(in)  :: im, km, nkc, nkt, ntrac, ntqv, ntke !, ndtend
 
   real(kind=kind_phys), intent(in) ::    zi(im, km+1),   zl(im, km), &
                                                          zm(im, km), &
                                        prsi(im, km+1), prsl(im, km)
   real(kind=kind_phys), intent(in) ::    dv(im, km),     du(im, km), &
                                         tdt(im, km),    rtg(im, km,ntrac)
   real(kind=kind_phys), intent(in) ::    u1(im, km),     v1(im, km),  t1(im,km)
   real(kind=kind_phys), intent(in) ::  dens(im, km),    dkt(im, km), dku(im,km)
!   real(kind=kind_phys), intent(in) :: dtend(im, km , ndtend)
 
! ** Q1 is concentration field (including gas and aerosol variables) mass mixing ratio kg kg-1
   real(kind=kind_phys), intent(in) ::    q1(im, km, ntrac)
 
   real(kind=kind_phys), intent(out) ::  &
! tendencies
!             DTEND_CAN (im, km , ndtend), &
             dv_can(im, km)        , &
             du_can(im, km)        , &
             tdt_can(im, km)        , &
! tendencies all gas-phase species & TKE
             rtg_can(im, km, ntrac) , &
 
! met3d arrays
             zl_can(im, nkt)        , & ! dim(nkt)
             zm_can(im, nkt)        , & ! dim(nkt)
             t1_can(im, nkt)        , & ! dim(nkt)
             qv_can(im, nkt)        , & ! dim(nkt)
             ws_can(im, nkt)        , & ! dim(nkt)
             prsl_can(im, nkt)        , & ! dim(nkt)
             prsi_can(im, nkt+1)      , & ! dim(nkt+1)
             dens_can(im, nkt)        , & ! dim(nkt)
             dkt_can(im, nkt)        , & ! dim(nkt)
             dku_can(im, nkt)        , & ! dim(nkt)
! all gas-phase species array
             q1_can(im, nkt, ntrac) , & ! dim(nkt)
             q1_2m(im,      ntrac) , & !
! canopy layers height arrays
             zmom_can3(im, nkt+1)      , & ! dim(nkt+1)    ! Paul's sigmcan(:,nkt+1)
             zmid_can3(im, nkt)        , & ! dim(nkt)      ! Paul's sigtcan(:,nkt+1)
             sigmom_can(im, nkt+1)      , & ! dim(nkt) ~ prsi(:,km+1)
             sigmid_can(im, nkt)            ! dim(nkt) ~ prsl(:,km)
 
   character(len=*), intent(out) :: errmsg
   integer,          intent(out) :: errflg
 
!...local variables
 
   integer        :: k, kc
 
! Initialize CCPP error handling variables
   errmsg = ''
   errflg = 0
 
! Initialize with values before in-canopy diffusion
 
! Layers height
   zmom_can3(:,:) = 0.
   zmid_can3(:,:) = 0.
   sigmom_can(:,:) = 0.
   sigmid_can(:,:) = 0.
 
! Zero in-canopy tendencies
!   dtend_can(:, :, : ) = 0.0
 
! Tracers
   q1_2m(:, :) = q1(:,1, :)       ! kg kg-1
 
! Subset (km combined layers minus top nkc layers)
   do k = 1, km-nkc
 
     ! km     is top combined subset
     ! nkc+1  is bot combined
     kc= nkc+k     ! 4th from top (nkt) to nkc+1 combined canopy plus resolved model layer
 
! Tendencies
     du_can(:,kc)    = du(:,k)       ! m s-2
     dv_can(:,kc)    = dv(:,k)       ! m s-2
     tdt_can(:,kc)    = tdt(:,k)       ! K s-1
 
     rtg_can(:,kc, 1:ntrac-1) = rtg(:,k, 1:ntrac-1) ! kg kg-1 s-1
     rtg_can(:,kc,   ntke   ) = rtg(:,k,   ntke   ) ! J   s-1 s-1
 
   end do
 
! All combined canopy plus resolved layers
   do k = 1, km
 
     ! nkc+km is top (nkt) combined
     ! nkc+1  is bot       combined
     kc= nkc+k     ! top (nkt) to nkc+1 combined canopy plus resolved model layer
 
! Height
     zl_can(:,kc) = zl(:,k)
     zm_can(:,kc) = zm(:,k)
 
! Pressure & temperature
     prsl_can(:,kc) = prsl(:,k)  ! km  combined canopy plus resolved layers
     prsi_can(:,kc) = prsi(:,k)  ! km  combined canopy plus resolved layers
     t1_can(:,kc) = t1(:,k)
     dens_can(:,kc) = dens(:,k)
 
! Diffusivities
     dkt_can(:,kc) = dkt(:,k)       ! m2 s-1
     dku_can(:,kc) = dku(:,k)       ! m2 s-1
 
! Wind
     ws_can(:,kc) = sqrt(u1(:,k)**2+v1(:,k)**2)       ! m s-1
 
! Mass tracers
     q1_can(:,kc, 1:ntrac-1) = q1(:,k, 1:ntrac-1) ! all tracers ntrac1
 
! TKE tracer
     q1_can(:,kc,   ntke  ) = q1(:,k,   ntke  ) ! ntke=198  TKE tracer
 
! Humidity
     qv_can(:,kc) = q1(:,k, ntqv) ! ntqv=1
 
   end do
   prsi_can(:,km + nkc +1 ) = prsi(:,km+1)  ! nkt combined canopy plus resolved layers
 
! Canopy layers
   do kc = 1, nkc ! 3-nkc canopy layers
 
! Tendencies
     du_can(:,kc)    = du(:,1)       ! m s-2
     dv_can(:,kc)    = dv(:,1)       ! m s-2
     tdt_can(:,kc)    = tdt(:,1)       ! K s-1
 
     rtg_can(:,kc, 1:ntrac-1) = rtg(:,1, 1:ntrac-1) ! kg kg-1 s-1
     rtg_can(:,kc,   ntke   ) = rtg(:,1,   ntke   ) ! J   s-1 s-1
 
! Height
     zl_can(:,kc) = zl(:,1)
     zm_can(:,kc) = zm(:,1)
 
! Pressure & temperature
     prsl_can(:,kc) = prsl(:,1)  ! km  combined canopy plus resolved layers
     prsi_can(:,kc) = prsi(:,1)  ! km  combined canopy plus resolved layers
     t1_can(:,kc) = t1(:,1)
     dens_can(:,kc) = dens(:,1)
 
! Diffusivities
     dkt_can(:,kc) = dkt(:,1)       ! m2 s-1
     dku_can(:,kc) = dku(:,1)       ! m2 s-1
 
! Wind
     ws_can(:,kc) = sqrt(u1(:,1)**2+v1(:,1)**2)       ! m s-1
 
! Mass tracers
     q1_can(:,kc, 1:ntrac-1) = q1(:,1, 1:ntrac-1) ! all tracers ntrac1
 
! TKE tracer
     q1_can(:,kc,   ntke  ) = q1(:,1,   ntke  ) ! ntke=198  TKE tracer
 
! Water vapor
     qv_can(:,kc) = q1(:,1, ntqv) ! ntqv=1
 
   end do
 
   end subroutine canopy_levs_init
 
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
   subroutine canopy_levs_run(im, km, nkc, nkt, &
              ntrac, ntqv, ntke,                    & ! in
              RDGAS, PI,                            & ! in ?? units ??
              zi, zl, zm,                           & ! in: 1D    zm(i,k) = zi(i,k+1)
              prsl, prsi, psfc,                     & ! in: 3D 3D 2D  (Pa)
              cfch,                                 & ! in: 2D
              garea, u10m, v10m, fm,fh,             & ! in: 2D
              rbsoil,                               & ! in: 2D
              t2m, q2m,                             & ! in: 2D
              stress, spd1,                         & ! in: 2D
              dv, du, tdt, rtg,                     & ! in: 3D
              u1, v1, t1, q1,                       & ! in: 3D " 4D
              dens, dkt, dku,                       & ! in  3D
              frt_mask,                             & ! in  2D
              kmod, kcan3,                          & ! out
              zmom_can3, zmid_can3,                 & ! out      zmom_can3(:, nkt+1)    zmid_can3(im, nkt)
              sigmom_can, sigmid_can,               & ! out 3D sigmom_can (:, nkt+1)  sigmid_can (im, nkt)
              zl_can, zm_can,                       & ! out 3D
              prsl_can, prsi_can,                   & ! out 3D prsi_can (:, nkt+1)
              dv_can, du_can, tdt_can, rtg_can,     & ! out: 3D
              t1_can,  qv_can, dens_can,            & ! out 3D
              ws_can, dkt_can,  dku_can,            & ! out 3D
              q1_can, q1_2m,                        &  !out
              errmsg,errflg)
 
   use machine , only : kind_phys
 
   IMPLICIT NONE
 
! Includes:
 
!...Arguments:
 
   integer, intent(in)  :: im, km, nkc, nkt, ntrac, ntqv, ntke
   real(kind=kind_phys), intent(in) :: rdgas, pi
 
   real(kind=kind_phys), intent(in) ::    zi(im,km+1),   zl(im,km), &
                                                         zm(im,km), &
                                        prsi(im,km+1), prsl(im,km)
   real(kind=kind_phys), intent(in) ::    dv(im,km),     du(im,km),         &
                                         tdt(im,km),    rtg(im,km,ntrac)
 
   real(kind=kind_phys), intent(in) ::    u1(im,km),     v1(im,km),  t1(im,km)
   real(kind=kind_phys), intent(in) ::  dens(im,km),    dkt(im,km), dku(im,km)
 
   real(kind=kind_phys), intent(in) :: psfc(im)   ! Pa
   real(kind=kind_phys), intent(in) :: cfch(im), garea(im), u10m(im), v10m(im),  &
                                       spd1(im),stress(im),                      &
                                        t2m(im),   q2m(im),   fm(im),   fh(im),  &
                                     rbsoil(im)
 
!     ** Q1 is concentration field (including gas and aerosol variables) kg kg-1
   real(kind=kind_phys), intent(in) :: q1(im, km, ntrac)
 
   real(kind=kind_phys), intent(in) :: frt_mask(im)
 
   integer, intent(out) ::  kmod  (im, km) , kcan3 (im, nkc)
 
   real(kind=kind_phys), intent(out) ::  &
! tendencies
             dv_can(im, km)         , &
             du_can(im, km)         , &
             tdt_can(im, km)         , &
! tendencies all gas-phase species & TKE
             rtg_can(im, km, ntrac)  , &
! met3d arrays
             zl_can(im, nkt)        , &
             zm_can(im, nkt)        , &
             t1_can(im, nkt)        , &
             qv_can(im, nkt)        , &
             ws_can(im, nkt)        , &
             prsl_can(im, nkt)        , &
             prsi_can(im, nkt+1)      , &
             dens_can(im, nkt)        , &
             dkt_can(im, nkt)        , &
             dku_can(im, nkt)        , &
! all gas-phase species array
             q1_can(im, nkt, ntrac) , &
             q1_2m(im,      ntrac) , &
! canopy layers height arrays
             zmom_can3(im, nkt+1)      , & ! Paul's sigmcan(:,nkt)
             zmid_can3(im, nkt)        , & ! Paul's sigtcan(:,nkt)
             sigmom_can(im, nkt+1)      , & ! ~ prsi(:,km+1)
             sigmid_can(im, nkt)            ! ~ prsl(:,km)
 
   character(len=*), intent(out) :: errmsg
   integer,          intent(out) :: errflg
 
!...Local arrays:
 
   integer(kind=4) :: kcan_top
   real   (kind=kind_phys) :: hcan
 
   logical         :: sfcflg(im)
 
   integer(kind=4) :: ka    (im)             , &
                      kl    (im)
 
   real(kind=kind_phys) ::      zmid3(km)  , &
                                zmom3(km)  , & ! Paul's zfull
                                sigmom3(km+1), &
                                z2(km+1), & ! Paul's    z2(:,chm_nk+1)
                                sigmid2(km+1), & ! Paul's sigt2(:,chm_nk+1)
             zcan3(nkc),                    &
             ta_can3(nkt),   ta3(km)  , &
             qv_can3(nkt),   qv3(km)  , &
             ws_can3(nkt),   ws3(km)  , &
             dkt_can3(nkt),   dkt3(km)  , &
             dku_can3(nkt),   dku3(km)  , &
             prsl_can3(nkt),   prsl3(km)  , &
             prsi_can3(nkt+1), prsi3(km+1), &
             dens_can3(nkt),   dens3(km)  , &
                                mol3(km)  , &
             klower_can(nkc)
 
   real(kind=kind_phys) ::                  &
             dxdy(im),  ustar(im), &
             ws10m(im),                &
             zol(im),  ilmo(im), &
                    safe_inv_mo_length(im)
 
!...local variables
 
   INTEGER          :: i,L
 
   logical(kind=4) :: flag_error
   integer(kind=4) :: k, kk, kc, k2, II, npass
 
   real(kind=kind_phys) :: tmp
   real(kind=kind_phys) :: hol, a1, b1, c1, rat
 
! del:  Minimum allowable distance between a resolved model layer and a canopy layer
!       (fraction of canopy layer height)
   real(kind=kind_phys),    parameter :: del = 0.2
   real(kind=kind_phys),    parameter :: min_kt = 0.1
   real(kind=kind_phys),    parameter :: zfmin=1.e-8
   real(kind=kind_phys),    parameter :: rimin=-100.
   real(kind=kind_phys),    parameter :: karman=0.4            ! von karman constant
   real(kind=kind_phys),    parameter :: threshold = 1.e06 ! MOL threshold, similar to mach_plumerise
   real(kind=kind_phys),    parameter :: epsilon = 1.e-10
 
   real(kind=kind_phys)    :: zm2, zr, td, hd, ddel
   real(kind=kind_phys)    :: uh, uspr, wndr, sigw, tl, ktr, kur
 
!  Assign the fractional heights of the canopy layers (fraction of canopy height)
   real(kind=kind_phys), dimension(3), parameter   :: can_frac = (/1.0, 0.5, 0.2/)
 
   logical(kind=4)                         :: local_dbg
 
! Initialize CCPP error handling variables
   errmsg = ''
   errflg = 0
 
   local_dbg = (.false.)
 
   kmod(:,:) = -999
   kcan3(:,:) = -999
 
! Initializations
 
! NB. mfpbltq_mod: q1(ix,km,ntrac1) kg kg-1
   q1_2m(:,   :) = q1(:,1,:)       ! kg kg-1
 
! Subset (km combined layers minus top nkc layers)
   do k = 1, km-nkc
     ! km     is top combined subset
     ! nkc+1  is bot combined
     kc= nkc+k     ! 4th from top (nkt) to nkc+1 combined canopy plus resolved model layer
 
! PBL Tendencies are declared in CCPP_typdefs as dim(im,km) instead (im, nkt)
     du_can(:,kc)    = du(:,k)       ! m s-2
     dv_can(:,kc)    = dv(:,k)       ! m s-2
     tdt_can(:,kc)    = tdt(:,k)       ! K s-1
 
     rtg_can(:,kc, 1:ntrac-1) = rtg(:,k, 1:ntrac-1) ! kg kg-1 s-1
     rtg_can(:,kc,   ntke  ) = rtg(:,k,   ntke  ) ! J   s-1 s-1
 
   end do
 
! All combined canopy plus resolved layers
   do k = km, 1, -1 ! top to 1hy model layer
     ! nkc+km is top (nkt) combined
     ! nkc+1  is bot       combined
     kc= nkc+k     ! top (nkt) to nkc+1 combined canopy plus resolved model layer
 
! Pressure & Temperature
     prsl_can(:,kc) = prsl(:,k)  ! km  combined canopy plus resolved layers
     prsi_can(:,kc) = prsi(:,k)  ! km  combined canopy plus resolved layers
     t1_can(:,kc) = t1(:,k)
     dens_can(:,kc) = dens(:,k)
 
! Diffusivities
     dkt_can(:,kc) = dkt(:,k)       ! m2 s-1
     dku_can(:,kc) = dku(:,k)       ! m2 s-1
 
! Wind
     ws_can(:,kc) = sqrt(u1(:,k)**2+v1(:,k)**2)       ! m s-1
 
! Mass tracers
     q1_can(:,kc, 1:ntrac-1) = q1(:,k, 1:ntrac-1) ! all tracers ntrac1
 
! TKE tracer
     q1_can(:,kc,   ntke  ) = q1(:,k,   ntke  ) ! ntke=198  TKE tracer
 
! Water vapor
     qv_can(:,kc) = q1(:,k, ntqv) ! ntqv=1
 
   end do
   prsi_can(:,nkt+1 ) = prsi(:,km+1)  ! km  combined canopy plus resolved layers
 
! Canopy layers
   do kc = 1, nkc ! 3-nkc canopy layers
 
     du_can(:,kc)    = du(:,1)       ! m s-2
     dv_can(:,kc)    = dv(:,1)       ! m s-2
     tdt_can(:,kc)    = tdt(:,1)       ! K s-1
 
     rtg_can(:,kc, 1:ntrac-1) = rtg(:,1, 1:ntrac-1) ! kg kg-1 s-1
     rtg_can(:,kc,   ntke  ) = rtg(:,1,   ntke  ) ! J   s-1 s-1
 
     prsl_can(:,kc) = prsl(:,1)  ! km  combined canopy plus resolved layers
     prsi_can(:,kc) = prsi(:,1)  ! km  combined canopy plus resolved layers
     t1_can(:,kc) = t1(:,1)
     dens_can(:,kc) = dens(:,1)
 
     dkt_can(:,kc) = dkt(:,1)       ! m2 s-1
     dku_can(:,kc) = dku(:,1)       ! m2 s-1
     ws_can(:,kc) = sqrt(u1(:,1)**2+v1(:,1)**2)       ! m s-1
 
! Mass tracers
     q1_can(:,kc, 1:ntrac-1) = q1(:,1, 1:ntrac-1) ! all tracers ntrac-1
 
! TKE tracer
     q1_can(:,kc,   ntke  ) = q1(:,1,   ntke  ) ! ntke=198  TKE tracer
 
! Water vapor
     qv_can(:,kc) = q1(:,1, ntqv) ! ntqv=1
 
   end do
 
 
   DO i = 1, im
 
      sfcflg(i)= .true.
      if(rbsoil(i) > 0.) sfcflg(i) = .false.
 
      dxdy(i) = garea( i ) !  dx*dy ~1.6E+8 m2
 
      ustar(i) = sqrt(stress(i))
!      ws10m(i)  = sqrt(u10m(i)**2+v10m(i)**2)
 
      zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
      if(sfcflg(i)) then
         zol(i) = min(zol(i),-zfmin)
      else
         zol(i) = max(zol(i),zfmin)
      endif
! Inverse of Monin-Obukhov length
      ilmo(i) = 1./zol(i)
 
!!!! Non-Canopy columns
   IF (frt_mask(i) <= 0.) THEN
 
      do k = 1, km     ! from bottom to top
         ii = km + 1 - k  ! from top to bottom of resolved model layers
!!! Paul's zmom is our zmom
! zmom(1)     = ZFULL(km) is top    model  layer height
! zmom(km) = ZFULL(1)     is bottom model  layer height
         ! NB. zm(:,k) = zi(:,k+1)
         zmom3(ii) = zm(i,k) ! ZFULL(i,k)
!! Heights of the original model layers for the canopy columns are extracted to the zmom array.
 
         ! Create temperature & humidity array on reversed layer order for interpolation
         ta3(ii) = t1(i,k)               ! K
         qv3(ii) = q1(i,k,1)      ! 1=water vapor kg kg-1
         prsl3(ii) = prsl(i,k)  ! Pa  mean layer pressure
         dens3(ii) = dens(i,k)  ! kg m-3
         ws3(ii) = sqrt(u1(i,k)**2+v1(i,k)**2)  ! rename wspd3 ???
         dkt3(ii) = dkt(i,k)  ! m2 s-2
         dku3(ii) = dku(i,k)  ! m2 s-2
      end do ! k = 1, km     ! from bottom to top
 
      do k = 1, km+1     ! from bottom to top
         ii = (km + 1) + 1 - k  ! from top to bottom of resolved model layers
 
         prsi3(ii) = prsi(i,k)  ! Pa  air pressure at model layer interfaces
! ! [pgr] surface air pressure meta var
         sigmom3(ii) = prsi(i, k) / psfc(i) ! PRES_FULL
      end do ! k = 1, km+1
 
!  First, carry over original model values for the matching layers
      do k = 1, km ! from bottom to top of resolved model layers
         ! kmod(1)     is 1  top model  layer
         ! kmod(km) is 65 top canopy layer (modified after mono adj.)
         kk = k ! kmod(i,k)
 
         sigmom_can(i, kk) = sigmom3(k) !
 
         ta_can3(kk) = ta3(k)   !          TA  (i, k) ! temperature [K]
         qv_can3(kk) = qv3(k)   ! Met_Data%QV  (i, k) ! spec. humidity
         prsl_can3(kk) = prsl3(k)   ! Met_Data%PRES(i, k) ! Pa
         prsi_can3(kk) = prsi3(k)   !
         dens_can3(kk) = dens3(k)   ! Met_Data%DENS(i, k) ! kg m-3
         ws_can3(kk) = ws3(k)   !                     ! m s-1
         dkt_can3(kk) = dkt3(k)   !          DKT (i, k) ! m2 s-2 atmos. thermal diffus.
         dku_can3(kk) = dku3(k)   !          DKU (i, k) ! m2 s-2 atmos. momentum diffus.
      end do
 
      do kc = 1, nkc       ! from top to bottom of canopy layers
         ! kk = 65 = kcan3(1) = km + 1
         ! kk = 66 = kcan3(2) = km + 2
         ! kk = 67 = kcan3(3) = km + 3
         kk = kc + km ! kcan3(i,kc)
 
!        zmom_can3 (i, kk) = zmom3  (km) ! full layer height [m]
         sigmom_can(i, kk) = sigmom3(km) !
 
         ta_can3(kk) = ta3(km)   !          TA  (i, k) ! temperature [K]
         qv_can3(kk) = qv3(km)   ! Met_Data%QV  (i, k) ! spec. humidity
         prsl_can3(kk) = prsl3(km)   ! Met_Data%PRES(i, k) ! Pa
         prsi_can3(kk) = prsi3(km)
         dens_can3(kk) = dens3(km)   ! Met_Data%DENS(i, k) ! kg m-3
         ws_can3(kk) = ws3(km)   !                     ! m s-1
         dkt_can3(kk) = dkt3(km)   !          DKT (i, k) ! m2 s-2 atmos. thermal diffus.
         dku_can3(kk) = dku3(km)   !          DKU (i, k) ! m2 s-2 atmos. momentum diffus.
      end do
! Lower interface at surface
      prsi_can3(   nkt+1) = prsi3(km+1)
      sigmom_can(i, nkt+1) = 1.0
 
!!!!! End non-canopy columns!!!!!
 
   ! Continuous forest canopy
   ELSE IF (frt_mask(i) > 0.) THEN
 
      hcan = cfch( i )
!!! Extract the canopy height (FCH)
 
! Generate initial canopy levels, as altitude above sea level
!
      do kc = 1, nkc
         zcan3(kc) = hcan * can_frac(kc)  ! Paul's hc is our hcan
                                          ! Paul's zcan is our zcan3
!!! Set the initial values of the heights of the inserted canopy layers to hc, 0.5 hc, and 0.2 hc
!!!
!!! NB. zcan3(1) is hc, top of canopy
!!!     zcan3(2) is 0.5 * hc
!!!     zcan3(3) is 0.2 * hc (bottom canopy level)
!
!        print*,'canopy_levs: ZCAN = ', i, kc, zcan3(kc)
      end do
 
! 1 = bottom (1st) model layer
! km= top model layer
      do k = 1, km     ! from bottom to top
         ii = km + 1 - k  ! from top to bottom of resolved model layers
! zl is height of layer center
! zmid3(1)  = zl(km) is top    model  layer height
! zmid3(km) = zl(1)  is bottom model  layer height
        ! Paul's zt is our zmid
        zmid3(ii) = zl(i,k) ! mid layer height [m]
!!! Heights of the original model layers for the canopy columns are extracted to the zmid array.
 
!       write(errmsg,*) 'canopy_levs: ZMID = ', i, II, zmid3(II)
 
        ! Paul's sigt2 is our sigmid2
        sigmid2(ii) = prsl(i,k)/ psfc(i)
 
      end do
      sigmid2(km+1) = 1.0
 
      do k = 1, km     ! from bottom to top
         ii = km + 1 - k  ! from top to bottom of resolved model layers
!!! Paul's zmom is our zmom
! zmom(1)     = ZFULL(km) is top    model  layer height
! zmom(km) = ZFULL(1)     is bottom model  layer height
         ! NB. zm(:,k) = zi(:,k+1)
         zmom3(ii) = zm(i,k) ! ZFULL(i,k)
!! Heights of the original model layers for the canopy columns are extracted to the zmom array.
 
         ! Create temperature & humidity array on reversed layer order for interpolation
         ta3(ii) = t1(i,k)               ! K
         qv3(ii) = q1(i,k,1)      ! 1=water vapor kg kg-1
         prsl3(ii) = prsl(i,k)  ! Pa  mean layer pressure
         dens3(ii) = dens(i,k)  ! kg m-3
         ws3(ii) = sqrt(u1(i,k)**2+v1(i,k)**2)
         dkt3(ii) = dkt(i,k)  ! m2 s-2
         dku3(ii) = dku(i,k)  ! m2 s-2
 
! From satmedmfvdifq.F:
! MOL  = zol(i)/zl(i,k)  !Monin-Obukhov Length in layer
! ZL is mid layer height [m]
         mol3(ii) = zol(i)/zl(i,k)  !Monin-Obukhov Length in layer
      end do
 
      do k = 1, km+1     ! from bottom to top
         ii = (km + 1) + 1 - k  ! from top to bottom of resolved model layers
 
         prsi3(ii) = prsi(i,k)  ! Pa  air pressure at model layer interfaces
! Paul's SIGM does not include surface layer lower interface (1.0) !!!
         sigmom3(ii) = prsi(i, k)/ psfc(i) ! PRES_FULL(i, k) / psfc(i)
 
      end do
 
!!! Find the resolved model level which lies above the top of the forest canopy,
!!! in each canopy column.  Usually the canopy is within the km or km-1
!!! level of the original model structure.
!
! The model level above the tallest canopy in grid
      kcan_top = 2         ! initialize to 2nd top model layer
      do l = km, 3, -1  ! from bottom to top model layer (going up)
         ! Mid-layer height m, zmid
         if (zmid3(l) > hcan) then ! Paul's zt is our zmid
            kcan_top = l - 1  ! level above the tallest canopy
            exit
         end if
      end do
! kcan_top = 62 or 63
!      print*,'canopy_levs: kcan_top = ', i, kcan_top
 
! MV2D_ILMO: Aggregated Inverse of Monin-Obukhov length
!  Setup of Monin-Obhukov Length similar to plumerise for upper limit:
! from satmedmfvdifq.F: ! MOL  = zol(i)/zl(i,k)  !Monin-Obukhov Length in layer
      safe_inv_mo_length(i) = ilmo(i)
      if (abs(ilmo(i)) > threshold) then
         safe_inv_mo_length(i) = sign(threshold, ilmo(i))
      end if
!
! Adjust the canopy levels: we don't want canopy levels to get closer than del (0.2m)
! to the model levels to prevent possible differencing errors in the diffusion.
! If zcan3 > zmid3 but is too close to zmid3, move zcan3 up by ddel.  If zcan3 < zmid3
! but is too close to zmid3, move zcan3 down by ddel.  The net result will be that the
! canopy levels are never closer than del from the original model levels.
      do k = kcan_top, km! from model layer above the canopy to bottom of model layer
         do kc = 1, nkc   ! from top to bottom of canopy
            if (abs(zmid3(k) - zcan3(kc)) < del) then
               ddel = max(0.0, del - abs(zcan3(kc) - zmid3(k)))
               zcan3(kc) = zcan3(kc) + sign(ddel, zcan3(kc) - zmid3(k))
 
!!! The reason why this section is necessary:  while it would be preferable for
!!! the canopy levels to stick with values of hc, 0.5 hc and 0.2 hc, somewhere in a
!!! large domain, there may be an overlap where one of these canopy levels is very
!!! close to or on top of an existing model level.  Which we dont want!
!!! What is done here, if the canopy levels come within "ddel" of an original model level,
!!! is to shift the canopy level in question a bit, to avoid overlaps.
            end if
         end do
      end do
 
!!! Starts the creation of the local array with the heights of the thermodynamic
!!! levels (layer midpoints) for the combined canopy + no canopy layers.
!!! Note that zmid_can at this point does not have these layers sorted in the
!!! correct order - the canopy layers have been tacked onto the bottom of the
!!! zmid_can array, but the values of zmid_can are not monotonically increasing with
!!! decreasing height index.
!
!  Set the initial values of the combined height array:
!
! Note that here, zmid_can is created, but the heights within each column have
! yet to be sorted to rearrange the layers in the correct order.
      do k = 1, km ! from top to bottom model layers
         zmid_can3(i,k)    = zmid3(k)
         ! Paul's zthrmcan is our zmid_can
!
      end do
 
! Add zcan3 additional thermo levels into zmid_can array for later sorting
      do kc = 1, nkc ! from top to bottom canopy layers
         zmid_can3(i,km+kc) = zcan3(kc)
      end do
 
!
!  Determine locations of canopy and resolved model levels within
!  the combined array for the canopy columns:
!
!!! This section sorts the zmid_can array to make sure that the new layers are
!!! all ordered so they monotonically increase with decreasing height.
 
! Top canopy layer height (km+1) is higher than the bottom model layer height (km)
      if (zmid_can3(i, km) < zmid_can3(i,km+1)) then
!
!  Non-trivial case:  the ancillary and original array levels intermingle.
!  Sort the combined height array to get the right order of the the heights:
!
!  zmid_can is the height locations of the combined array, which needs to be  sorted:
!  since there are only NKC levels in the canopy, and both zcan3 and z
!  decrease monotonically, only nkc+1 passes are needed to sort the combined array:
         do npass = 1, nkc+1
            flag_error = .false.
            do k = nkt, 2, -1
! Top    canopy layer height (nkt-2) is larger than the bottom model  layer height (nkt-3 = km)
! Middle canopy layer height (nkt-1) is larger than the top    canopy layer height (nkt-2)
! Bottom canopy layer height (nkt)   is larger than the middle canopy layer height (nkt-1)
               if (zmid_can3(i, k) > zmid_can3(i, k-1)) then
!  The combined array heights are out of order, sort them:
                  tmp = zmid_can3(i, k-1)
                  zmid_can3(i, k-1) = zmid_can3(i,k)
                  zmid_can3(i, k)   = tmp
                  flag_error = .true.
               end if
            end do
         end do
         if (flag_error) then
            write(errmsg,fmt='(*(a))') 'NKC+1 passes insufficient to sort canopy array ' // &
                                       'in canopy_levs.F90.  Scream and die.'
            errflg = 1
            return
         end if
      end if
 
!
!  Heights in zmid_can should now be monotonically decreasing.
 
!  Next, identify the locations of the vertical levels in the combined
!  array relative to the resolved model array and canopy array
!
!!! Now that the heights in zmid_can are in the right order, we can use them to
!!! identify the values of kcan and kmod:  the vertical locations of the canopy and
!!! original model layers in the augmented canopy layer code.
      do kc = 1, nkc         ! from top to bottom canopy layers
         do kk = nkt, 1, -1  ! from bottom to top of combined canopy and resolved model levels
            if (zmid_can3(i, kk) == zcan3(kc)) then
               kcan3(i,kc) = kk
               exit
            endif
         end do
      end do
 
! k=1        is top    model layer
! k=km       is bottom model layer
      do k = 1, km ! from top to bottom model layers
         do kk = k, nkt! from bottom to top of combined resolved plus canopy layers
 
! zmid_can3(1)    = zmid3(1)     is top    model  layer height
! ...
! zmid_can3(km)= zmid3(km) is bottom model  layer height
            if (zmid_can3(i, kk) == zmid3(k)) then
 
! kmod(1)       is 1      , top model  layer
! kmod(km-1) is km-1, 2nd model  layer
! kmod(km)   is          top canopy layer (modified after monotonic adj.)
               kmod(i,k) = kk
               exit
            endif
         end do
      end do
 
      if (local_dbg) then
      do kc = 1, nkc
         if (kcan3(i,kc) < 1) then
!           write(errmsg,fmt='(*(a,i0,a,i0))') 'canopy_levs: kcan undefined: kc=', kc, &
!                                             ' kcan3=', kcan3(i,kc)
            write(errmsg,*) 'canopy_levs: kcan undefined: ', kc, kcan3(i,kc)
            errflg = 1
            return
         end if
      end do
      do k = 1,km
         if (kmod(i,k) < 1) then
!           write(errmsg,fmt='(*(a,i0,a,i0))') 'canopy_levs: kmod undefined: k=', k, &
!                                             ' kmod=', kmod(i,k)
            write(errmsg,*) 'cannopy_levs: kmod undefined: ',k, kmod(i,k)
            errflg = 1
            return
         end if
      end do
      end if
 
 
!  Create the corresponding momentum height array
!
!  The original methodology adopted made use of the at2m array and the  thermodynamic heights determined above.
!  However, this methodology resulted in momentum levels which did not match the original model levels
!  above the region modified for canopy layers.  Here, the thermodynamic layers will be used to
!  (1) Determine whether the original model and canopy thermodynamic layers coincide, and if so,
!  (2) Use the existing model layer values for the momentum layers, while if not,
!  (3) Assign the new momentum layers as being 1/2 way between the canopy layers
! Note that these changes only exist inside the chemistry part of GEM-MACH and do not affect the model physics
!!!
!!! Create the momentum height (layer interface) array.  The original momentum layers are used above the canopy height.
!!! Below the canopy height, the "momentum" layers are assumed to be 1/2 way between the thermodynamic layers.
 
! Default case:  all added canopy thermodynamic layers are below the lowest resolved model thermodynamic layer
! kcan_top is either 2nd or 3rd (63 or 62) resolved model layer
      do k = 1, kcan_top - 1 ! from top model layer to model layer above the canopy
!          zmom(1) is top    model  layer height
! zmom(kcan_top-1) is model layer above the canopy < 234.061m
         zmom_can3(i,k) = zmom3(k) ! full layer height [m]
      end do
 
      ka(i) =km
      inner0: do k = kcan_top, km-1  ! from resolved model layer above the canopy to top model layer
!  Starting from the top, scan down through the original and combined mid layer heights, to see when
!  they first deviate from each other
         !Paul's zthrmcan is our zmid_can
         !Paul's zt is our zmid
         if (zmid_can3(i,k) == zmid3(k) .and. zmid_can3(i,k+1) == zmid3(k+1)) then
            ! Paul's zmom is our zmom
            ! Paul's zmomcan(nkt+1) is our zmom_can
            zmom_can3(i,k) = zmom3(k)  ! full layer height [m]
         else
            ka(i) = k
            exit inner0
         end if
      end do  inner0
 
! ka is 63 or 64
!      print*,'canopy_levs: ka = ', i, ka(i)
 
! ka  is the lower-most layer for which the combined layer zmom_can = zmom resolved model layer
      ! Paul's zmom is our zmom
      ! Paul's zmomcan is our zmom_can
      zmom_can3(i,ka(i)) = zmom3(ka(i))
      do k = ka(i)+1, nkt! from ka to bottom combined canopy and resolved layers
         zmom_can3(i,k) = (zmid_can3(i,k-1) + zmid_can3(i,k)) * 0.5
 
      end do
 
!########################################################################
 
! create original model arrays of z and sigma-t which include the surface, to
!  allow interpolation:
      !Paul's sigtcan is our sigmid_can
      sigmid_can(:,:) = 0.0
      do k = 1, km! from top to bottom of resolved model layers
 
! zmid3(1)  is top    model  layer height
! zmid3(km) is bottom model  layer height
         z2(k)    = zmid3(k)
 
!  Fill in the thermodynamic sigma levels (Pre-existing levels first):
! kmod(1)     is 1  top  (last) model layer
! kmod(km) is 64 bottom (1st)  model layer
         kk = kmod(i,k)
         sigmid_can(i, kk) = sigmid2(k)
 
      end do
      klower_can(:) = -999
      z2(km+1) = 0.0
 
!  fill in the remaining sigma levels by interpolating in z:
      do kc = 1, nkc   ! from top to bottom canopy layers
         do k2 = kcan_top, km+1  ! from resolved model layer above the canopy to top model layer
            if (zcan3(kc) > z2(k2) .and. zcan3(kc) <= z2(k2-1)) then
 
! Interpolate in sigma
               sigmid_can(i, kcan3(i,kc)) = sigmid2(k2-1)  +          &
                                 (sigmid2(k2) - sigmid2(k2-1)) /      &
                                   (   z2(k2) -    z2(k2-1)) *        &
                                   (zcan3(kc) -    z2(k2-1))
 
! Store grid locations for use in later interpolations
               klower_can(kc) = k2
            end if
 
         end do
 
         if (klower_can(kc) < 1) then
            errflg = 1
 
            write(errmsg,*) 'canopy_levs:  klower_can is unassigned at i, kc: ', i, kc
            write(errmsg,*) 'canopy_levs:  zcan3(kc): ',zcan3(kc)
!
            do kk = kcan_top, km+1
               write(errmsg,*) 'canopy_levs: kk z2(kk) which should bracket the above zcan3: ',kk, z2(kk)
            end do
            do kk = 1, km+1
               write(errmsg,*) 'canopy_levs:  kk z2(kk) full set of z2 values: ', kk, z2(kk)
            end do
            do kk = 1,nkc
               write(errmsg,*) 'canopy_levs:  kc zcan3(kc) hcan fr(kc) for full set of zcan3 values: ',kk, zcan3(kk), hcan, can_frac(kk)
            end do
 
            return
         end if
      end do
 
! NB. klower_can(1) is 64 or 65
!     klower_can(2) is 65 except for individual grid points near West coast
!     klower_can(3) is 65 uniformly
!
!
   if (local_dbg) then
!  Check on klower_can for NaN or out of bounds:
   do kc = 1,nkc
      if ((klower_can(kc) /= klower_can(kc)) .or. &
          (klower_can(kc) <= 0)              .or. &
          (klower_can(kc) > km+ 1) ) then
!        write(errmsg,*) 'canopy_levs: klower_can after creation NaN or <=0 or >km+1 : ', &
!                      kc, klower_can(kk)
         return
      end if
   end do
   end if
!
!  Create sigma coordinate  momentum levels:
!
!  As above, the existing momentum levels and the canopy values are used to create SIGM levels
!
!  (1) Determine whether the original model and canopy thermodynamic layers coincide, and if so,
!  (2) Use the existing model layer values for the momentum layers, while if not,
!  (3) Assign the new momentum layers as being 1/2 way between the canopy layers
! Note that these changes only exist inside the chemistry part of GEM-MACH and do not affect the
! model physics
 
 
!      Default case:  all added canopy half layers are
!      below the lowest resolved model half layer
      ka(i) = km
      inner2:   do k = 1, km-1
         if (sigmid_can(i, k) == sigmid2(k) .and. sigmid_can(i, k+1) == sigmid2(k+1) ) then
            sigmom_can(i, k) = sigmom3(k)
         else
            ka(i) = k
            exit inner2
         end if
      end do inner2
! ka  is the last layer for which sigmom_can= sigmom3(k)
      sigmom_can(i, ka(i)) = sigmom3(ka(i))
      do k = ka(i)+1,nkt
         sigmom_can(i, k) = (sigmid_can(i, k-1) + sigmid_can(i, k)) * 0.5
      end do
      sigmom_can(i, nkt+1) = 1.0
 
!  Next, do a sort of all of the variables in the original METV3D array into canopy.  Note that
!  the declaration of the met arrays for the new canopy subdomain has occurred earlier in the code.
!  Three-D variables are a bit more complicated, in that one must make decisions regarding
!  the values of the met variables in the canopy region.
!  The code which follows is based on chm_load_metvar.ftn90
!
!  First, carry over original model values for the matching layers
   do k = 1, km ! from bottom to top of resolved model layers
      ! kmod(1)     is 1  top model  layer
      ! kmod(km) is 65 top canopy layer (modified after mono adj.)
      kk = kmod(i,k)
      ta_can3(kk) = ta3(k)   !          TA  (i, k) ! temperature [K]
      qv_can3(kk) = qv3(k)   ! Met_Data%QV  (i, k) ! spec. humidity
      prsl_can3(kk) = prsl3(k)   ! Met_Data%PRES(i, k) ! Pa
      prsi_can3(kk) = prsi3(k)   !
      dens_can3(kk) = dens3(k)   ! Met_Data%DENS(i, k) ! kg m-3
      ws_can3(kk) = ws3(k)   !                     ! m s-1
      dkt_can3(kk) = dkt3(k)   !          DKT (i, k) ! m2 s-2 atmos. thermal diffus.
      dku_can3(kk) = dku3(k)   !          DKU (i, k) ! m2 s-2 atmos. momentum diffus.
 
   end do ! km
 
!----------------------------------------------------------------------------
!  Canopy region:  next, go through each variable to work out canopy values.
!
!  (1) Do those variables for which special canopy formulae will NOT be used:
      do kc = 1, nkc       ! from top to bottom of canopy layers
 
!  Each of the following 2 variables have a screen height (2m) value in the 2D met arrays
!          Temperature:         TA, T2M
!          Specific humidity:   Q,  Q2M
 
!   kcan3(1) = 65
!   kcan3(2) = 66
!   kcan3(3) = 67
         kk = kcan3(i,kc)
         if (klower_can(kc) <= km) then
!  Level is above first resolved model level
 
            k2 = klower_can(kc)
            zm2 = (zcan3(kc) - z2(k2-1)) / (z2(k2) - z2(k2-1))
!           zm2 = (zcan3(kc) - z2(k2-1)) / max(z2(k2) - z2(k2-1), epsilon)
 
            td = ( ta3(k2)  - ta3(k2-1)) * zm2
            hd = ( qv3(k2)  - qv3(k2-1)) * zm2
            ta_can3(kk) = ta3(k2-1) + td
            qv_can3(kk) = qv3(k2-1) + hd
 
         else
!  Level is below first resolved model level
 
            if (zcan3(kc) - z2(km+1) >= 2.0) then
            !  Level is below first resolved model level but above screen height
 
               zm2 = (zcan3(kc) - z2(km+1) - 2.0) / (z2(km) - z2(km+1) - 2.0)
!              zm2 = (zcan3(kc) - z2(km+1) - 2.0) / max(z2(km) - z2(km+1) - 2.0, epsilon)
 
               td = (ta3(km)  - t2m( i ) )  * zm2
               hd = (qv3(km)  - q2m( i ) )  * zm2
               ta_can3(kk)  = t2m( i ) + td
               qv_can3(kk)  = q2m( i ) + hd
 
            else
            ! Level in canopy is below screen height; assume constant values below screen height
 
               ta_can3(kk)  = t2m( i ) ! 2-m  temperature [K]
               qv_can3(kk)  = q2m( i ) ! 2-m  spec. humidity
            end if
 
         end if
 
!  Evaluate the air density in canopy columns using values determined above
!
! NB. PRSL  is air pressure on ZL, formerly ZH (mid-layers)
!     PRSI  is air pressure on ZI, formerly ZF (interfaces)
!     psfc   is surface air pressure psfc
 
! get pressure from sigma levels in Pa
         prsl_can3(kk) = sigmid_can(i, kk) * psfc(i)  ! ~zl mid-layers centers
         prsi_can3(kk) = sigmom_can(i, kk) * psfc(i)  ! ~zm/zi layers interfaces
 
! aqm_methods: dens: buffer(k) = stateIn % prl(c,r,l) / ( rdgas * stateIn % temp(c,r,l) )
         ! (1)       is top model  layer
         ! ...
         ! (km)   is 1hy model  layer
         ! (km+1) is top canopy layer
         ! (nkt)   is 1st canopy layer
         dens_can3(kk) = prsl_can3(kk) / ( rdgas * ta_can3(kk))  ! kg m-3
 
 
!  The following variables are assumed to have uniform values throughout the
!  lowest resolved model layer:
!
!    Cloud liquid water mass mixing ratio (QCPLUS)
!    Total cloud fraction (FTOT)
!    Stratospheric cloud fraction (FXP)
!    Convective cloud fraction (FDC)
!    Total liquid water flux (RNFLX)
!    Total solid water flux (SNOFLX)
!    Precipitation evaporation (FEVP)
!    Cloud to rain collection tendency (PPRO)
!  Search over the original model layers (k).  Note that the outer loop above this
!  one is over the canopy layers kc:  we are looking for the values to assign the
!  canopy layers in the combined canopy+resolved scale space.  For these variables,
!  the resolved scale values will be used, hence the aim is to determine the
!  resolved scale layer in which the canopy layer resides, and assign the
!  corresponding values to the locations of the canopy layers in the combined
!  canopy + resolved scale space (kk).
 
 
      end do ! kc = 1,nkc
! Surface layer lower interface
      prsi_can3(nkt+1) = prsi3(km+1)
 
   if (local_dbg) then
! Several checks for suspicious values:
      do kk = 1,nkt
         if ( ta_can3(kk) < 150.0) then
            write(errmsg,*) 'canopy_levs:  suspicious temperature detected in canopy_levs after creation (kk value): ',&
                        i, kk, ta_can3(kk)
            errflg = 1
 
            do kc = 1, nkc
               write(errmsg,*) 'canopy_levs: value of zcan(kc) z2(km+1) and difference  at this value of ic for kk: ',&
                            kc,' are: ',zcan3(kc),z2(km+1), zcan3(kc)-z2(km+1)
            end do
 
            do k = 1, nkt
               write(errmsg,*) 'canopy_levs: value of zmid_can for = ', i,' at k = ',k,' is: ',zmid_can3(i,k)
            end do
 
            do kc = 1,nkc
               write(errmsg,*) 'canopy_levs:  values of kcan zcan and original zcan for = ', i,' at kc = ',kc,' are: ',&
                           kcan3(i,kc), zcan3(kc), hcan * can_frac(kc)
            end do
 
            do k = 1,km
               write(errmsg,*) 'canopy_levs:  value of kmod and z for = ', i,' at k = ',k,' are: ',kmod(i,k), zmid3(k)
            end do
 
            do kc = 1,nkc
               write(errmsg,*) 'canopy_levs: value of klower_can at this grid point for kc: ',kc,' is: ',klower_can(kc)
            end do
 
            return
         end if
      end do
   end if
 
!  (2) For the last few variables, the value at the lowest resolved model layer and typical profiles for that variable
!  within the canopy will be used to create the canopy values:
      do kc = 1, nkc
         kk = kcan3(i,kc)
!  Ratio of lowest model level to canopy height:
!
         zr = (zmid3(km) - z2(km+1)) / hcan
!
!  Horizontal wind and KT profiles are from Raupach, Quarterly Journal
!  of the Royal Meteorological Society, vol 115, pp 609-632, 1989, examples
!  from page 626, equations (48) through (51).
!
!  Wind speed (equation 51), assumed to scale similarly in each horizontal dimension:
!
!   U(z) = ustar/karman * ln((z - d) / z0), where
!   k = 0.4
!   d = 0.75 hc
!   z0 = 0.07530 hc
! The next few lines calculate the average value of u(z), v(z), Raupach's eqn 51,
! at the first resolved level model height
! Paul's UE is our ustar,  surface friction velocity
         uh = ustar(i) * 3.0
         if (zr >= 1.0) then
         ! Paul's zt is our zmid (i.e. zmid(km) is zt(i,chm_nk))
         ! Paul's hc is our hcan
            uspr = ustar(i) / karman * &
                   log((zmid3(km) - z2(km+1) - 0.75 * hcan) / &
                   (0.07530 * hcan))
!                  log(max((zmid3(km) - z2(km+1) - 0.75 * hcan) / &
!                  (0.07530 * hcan), epsilon))
         else
            uspr = uh * exp(- 2.0 * ( 1.0 - zr))
         end if
!  wndr is the ratio of the wind to Raupach's average us(), eqn 51.
!  This is used to scale the wind speed with height values from eqn 51 to the current grid square
         ! Paul's WS(nk) is our spd1, wind speed at lowest model level m s-1
         wndr = spd1(i) / uspr
!        wndr = spd1(i) / max(uspr, epsilon)
!  Using Raupach's formulae for wind speed, multiplied by the above ratio, for the canopy layers:
!
         zr = (zcan3(kc) - z2(km+1)) / hcan
         if (zr >= 1.0) then
            uspr = log((zcan3(kc) - z2(km+1) - 0.75 * hcan) / &
                   (0.07530 * hcan)) * ustar(i)
         else
            uspr = uh * exp(- 2.0 * (1.0 - (zcan3(kc) - z2(km+1)) / hcan))
         end if
 
         ws_can3(kk) = wndr * uspr
!
!  Coefficients of diffusivity:
!  Find value of K at first model level from raupach's sigw and TL formulae (eqns 48, 49)
         zr = (zmid3(km) - z2(km+1)) / hcan
!  Gradient in stability under the canopy is reduced for higher stability conditions
!  in accord with Shaw, den Hartog and Neumann, BLM 45, 391-409, 1988, Fig 16.
         ! Paul's zl is our hol (as in satmedmfvdifq.F)
         hol = hcan * safe_inv_mo_length(i)
! Unstable:
         if(hol < -0.1) then
             a1 = 0.75
             b1 = 0.5
             c1 = 1.25
         end if
! Neutral:
         if(hol >= -0.1 .and. hol < 0.1) then
             a1 = 0.625
             b1 = 0.375
             c1 = 1.0
         end if
! Stable:
         if(hol >= 0.1 .and. hol < 0.9) then
             rat = 4.375 - 3.75 * hol
             a1 = 0.125 * rat + 0.125
             b1 = 0.125 * rat - 0.125
             c1 = 0.25 * rat
         end if
! Very stable (from extrapolation of Shaw et al's values at 0.1 and 0.5:
         ! Paul's MV3D_KT(nk) is our dkt3(km) m2 s-1 atmospheric heat diffusivity (thermal vertical diffusion coefficient)
         !  1st (bottom) model layer
         if(hol >= 0.9 .or. dkt3(km) <= min_kt) then
             a1 = 0.25
             b1 = 0.0
             c1 = 0.25
         end if
!  Raupach's originals:
!         if (zr >= 1.0) then
!            sigw = ustar(i) * 1.25
!         else
!            sigw = ustar(i) * ( 0.75 + 0.5 * cos(pi * (1.0 - (zmid3(km) - z2(km+1))/hcan) ) )
!         end if
!  Replace Raupach's originals with fit to Patton et al and Shaw et al 1988
         if(zr < 0.175) then
               sigw = ustar(i) * 0.25
         else
           if(zr < 1.25) then
               sigw = ustar(i) * ( a1 + b1 * cos(pi / 1.06818 * &
                      (1.25 - (zmid3(km) - z2(km+1)) / hcan)))
           else
               sigw = ustar(i) * c1
           end if
         end if
 
!        tl = hcan / max(ustar(i), epsilon)  * &
         tl = hcan / ustar(i)  * &
              (0.256 * ((zmid3(km) - z2(km+1) - 0.75 * hcan) / hcan) + &
               0.492 * exp(-(0.256 * ((zmid3(km) - z2(km+1)) / hcan) / 0.492)))
! ktr is the ratio of the resolved model diffusivity at the lowest resolved
! model level to that derived by Raupach's formula
!
         ktr =  dkt3(km) / (sigw * sigw * tl)
         kur =  dku3(km) / (sigw * sigw * tl)
!        ktr =  dkt3(km) / max(sigw * sigw * tl, epsilon)
!        kur =  dku3(km) / max(sigw * sigw * tl, epsilon)
 
!         print*, 'canopy_levs: KTR= ', i, ktr, dkt3(km), kk, kc
!
!  Use Raupach's formulae for diffusivity, multiplied by the above ratio, for the canopy layers:
!
         zr = (zcan3(kc) - z2(km+1)) / hcan
!  Gradient in stability under the canopy is reduced for higher stability conditions
!  in accord with Shaw, den Hartog and Neumann, BLM 45, 391-409, 1988, Fig 16.
!  Raupach's original:
!         if (zr >= 1.0) then
!            sigw = ustar(i) * 1.25
!         else
!            sigw = ustar(i) * ( 0.75 + 0.5 * cos(pi * (1.0 - (zcan3(kc) - z2(km+1))/hcan) ) )
!         end if
         if(zr < 0.175) then
               sigw = ustar(i) * 0.25
         else
           if(zr < 1.25) then
               sigw = ustar(i) * ( a1 + b1 * cos(pi / 1.06818 * &
                      (1.25 - (zcan3(kc) - z2(km+1))/hcan)))
           else
               sigw = ustar(i) * c1
           end if
         end if
!
!        tl = hcan / max(ustar(i), epsilon) *  &
         tl = hcan / ustar(i) *  &
              (0.256 * ( (zcan3(kc) - z2(km+1) - 0.75 * hcan) / hcan) + &
              (0.492 * exp(-(0.256 * (zcan3(kc) - z2(km+1)) / hcan) / 0.492) ) )
 
         dkt_can3(kk)  = (sigw * sigw * tl) * ktr
         dku_can3(kk)  = (sigw * sigw * tl) * kur
 
!        print*, 'canopy_levs: DKT_CAN= ', i, sigw, tl, dkt_can3(kk), kk, kc
      end do ! kc = 1,nkc
!
   if (local_dbg) then
      do kc = 1, nkc
         if (kcan3(i, kc) == 0) then
            write(errmsg,fmt='(*(a,i0,a,i0))') 'kcan zero inside canopy_levs at i=', i, ' kc=', kc
            errflg = 1
            return
         end if
      end do
   end if
!
      do k = 1, nkt! from top to bottom of combined layers
         ii = nkt + 1 - k  ! from bottom to top of combined layer
 
         ! Flip back meteo arrays on combined layers in same layer order as original model layer
         ! nkt  is top model layer         <= 1
         ! ...
         ! (4) is 1st (bottom) model layer   <= km
         ! (3) is 3rd (top) canopy layer     <= nkt-2
         ! (2) is 2nd canopy layer           <= nkt-1
         ! (1) is 1st (bottom) canopy layer  <=nkt
         zl_can(i,ii) = zmid_can3(i, k)
         zm_can(i,ii) = zmom_can3(i, k)
         prsl_can(i,ii) = prsl_can3(k)
 
         t1_can(i,ii) = ta_can3(k)
         qv_can(i,ii) = qv_can3(k)
         dens_can(i,ii) = dens_can3(k)
         ws_can(i,ii) = ws_can3(k)
         dkt_can(i,ii) = dkt_can3(k)
         dku_can(i,ii) = dku_can3(k)
 
      end do ! k = 1, nkt
 
! Pressure at layers interfaces
      do k = 1, nkt+1  ! from top to bottom of combined layers
         ii = (nkt+1) + 1 - k  ! from bottom to top of combined layer
 
! (km+1) (68=nkc+km +1) prsi3( 1)   Top    model layer upper interface  prsi_can3(1)
! (km)   (67=nkc+km   ) prsi3( 2)
! ...
! (2)  (5 =nkc    +2)   prsi3(km)   Bottom model  layer upper interface  prsi_can3(km)
!      (4 =nkc    +1)               Top    canopy layer upper interface  prsi_can3(km+1)
!      (3)                          Mid    canopy layer upper interface
!      (2)                          Bottom canopy layer upper interface prsi_can3(nkt)
! (1)  (1)              prsi3(km+1) Bottom model  layer LOWER interface prsi_can3(nkt+1)
!
         prsi_can(i,ii) = prsi_can3(k)
 
      end do ! k = 1, nkt+1
 
 
   END IF ! Continuous forest canopy: FRT_MASK == 1.
 
! ... have not finished Paul's code ...
 
 
   END DO  !I-index
 
   end subroutine canopy_levs_run
 
   end module canopy_levs_mod