!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Copyright (C) 2011 UNESCO-IHE, WL|Delft Hydraulics and Delft University !
! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski,            !
! Jaap van Thiel de Vries, Robert McCall                                  !
!                                                                         !
! d.roelvink@unesco-ihe.org                                               !
! UNESCO-IHE Institute for Water Education                                !
! P.O. Box 3015                                                           !
! 2601 DA Delft                                                           !
! The Netherlands                                                         !
!                                                                         !
! This library is free software; you can redistribute it and/or           !
! modify it under the terms of the GNU Lesser General Public              !
! License as published by the Free Software Foundation; either            !
! version 2.1 of the License, or (at your option) any later version.      !
!                                                                         !
! This library is distributed in the hope that it will be useful,         !
! but WITHOUT ANY WARRANTY; without even the implied warranty of          !
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU        !
! Lesser General Public License for more details.                         !
!                                                                         !
! You should have received a copy of the GNU Lesser General Public        !
! License along with this library; if not, write to the Free Software     !
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307     !
! USA                                                                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module canopyflow_module

    subroutine canflow(s,par,veg)
        use params
        use spaceparams
        use readkey_module
        use xmpi_module
        use filefunctions
        use interp
        use logging_module

        type(parameters)                            :: par
        type(spacepars)                             :: s
        type(veggie), dimension(:), pointer         :: veg

        ! local variables
        integer                                     :: i,j,ind,m
        real*8                                      :: Cdterm

        use typesandkinds
        implicit none
        
        ! Initialization

        h_can =0.1d0 ! Ryan: define canopy parameters
        lam_f =0.5d0
        lam_p =0.2d0
        Cd_can =2.0d0
        CM_can =2.0d0

        Ld_can=2.*h_can*(1-lam_p)/(Cd_can*lam_f) ! Ryan: Canopy drag lengthscale (Lowe et al. 2005) Eq.(20)
        M=1+CM_can*lam_p/(1-lam_p) ! Ryan: added mass term for canopy momentum equation
   
        ! Ryan: In canopy momentum balance (Lowe et al. 2005)
        do j=j1,max(ny,1)
            if (xmpi_isbot) then
                imax = nx-1
            else
                imax=nx
            endif
            do i=2,imax
                if(wetu(i,j)==1) then
                    ! Option 1 (canopy drag computed from velocity at previous time step)
                    u_can(i,j)=u_can(i,j)-par%dt/M*(par%g*dzsdx(i,j) &
                    + ABS(u_can(i,j))*u_can(i,j)/Ld_can) ! Ryan: Add canopy dissipation term
                    ! Option 2 (canopy drag computed by staggering velocity across time steps)
                    !D_can(i,j)=M/par%dt+ABS(u_can(i,j))/Ld_can ! Ryan: Eq. (A5) in Lowe et al. (2007)
                    !u_can(i,j)=1/D_can(i,j)*(M/par%dt*u_can(i,j)-par%g*dzsdx(i,j)) ! Ryan: Eq. (A4) in Lowe et al. (2007)
                else
                    u_can(i,j)=0.0d0
                end if
            end do
        end do
        ! Lateral boundary conditions for u_can
        if (ny>0) then
            if (xmpi_isleft) then !Dano/Robert only on outer boundary
                u_can(1:nx+1,1)=u_can(1:nx+1,2) ! RJ: can also be done after continuity but more appropriate here
            endif
            ! Lateral boundary at y=ny*dy
            if (xmpi_isright) then !Dano/Robert only at outer boundary
                u_can(1:nx+1,ny+1)=u_can(1:nx+1,ny) ! RJ: can also be done after continuity but more appropriate here
            endif
        endif
    
        Fcanu   = abs(u_can(i,j))*u_can(i,j)*h_can/(Ld_can*hu(i,j))
        Fcanv   = 0.d0
    
        s%Fcanu = Fcnu
        s%Fcanv = Fcnv
    
end module canopyflow_module

! FLOW TIMESTEP:
!if(wetu(i,j)==1) then
!413 uu(i,j)=uu(i,j)-par%dt*(ududx(i,j)+vdudy(i,j)-viscu(i,j) & !Ap,Robert,Jaap
!414 + par%g*dzsdx(i,j) &
!415 + taubx(i,j)/(par%rho*hu(i,j)) & ! Dano: changed hum to hu NOT cf volume approach
!416 - par%lwave*Fx(i,j)/(par%rho*max(hum(i,j),par%hmin)) &
!417 - fc*vu(i,j) &
!418 - par%rhoa*par%Cd*windsu(i,j)**2/(par%rho*hum(i,j)) &
!419 + ABS(u_can(i,j))*u_can(i,j)*h_can/(Ld_can*hu(i,j))) ! Ryan: Add canopy dissipation term
!420 else
!421 uu(i,j)=0.0d0
!422 end if