!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
! Copyright (C) 2007 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 bedroughness_module



implicit none

real*8,dimension(:,:),allocatable,save,private  :: kru,krv,kru50,krv50,kru90,krv90
real*8,dimension(:,:),allocatable,save,private  :: urms_upd,u2_upd
real*8,dimension(:,:),allocatable,save,private  :: facbl,blphi,infilb,Ubed,Ventilation
real*8,dimension(:,:),allocatable,save,private  :: ueuf,uevf,vevf,veuf
real*8,dimension(:,:),allocatable,save,private  :: ueuold,uevold,vevold,veuold
real*8,dimension(:,:),allocatable,save,private  :: dudtsmooth,dvdtsmooth
real*8,dimension(:,:),allocatable,save,private  :: shieldsu,shieldsv
real*8,dimension(2),save,private                :: dtold
real*8,save,private                             :: delta,rhogdelta
real*8,private                                  :: Fd,Fi,Fl,Wg
! parameters for Van Gent model
!real*8,parameter,private                        :: cd = 0.018d0
!real*8,parameter,private                        :: cm = 0.08d0
!real*8,parameter,private                        :: cl = 0.075d0
!real*8,parameter,private                        :: k1 = 0.66d0
!real*8,parameter,private                        :: k2 = 0.90d0
! parameters for Okayasu
!real*8,parameter,private                        :: cd = 0.40d0
!real*8,parameter,private                        :: cm = 1.50d0
!real*8,parameter,private                        :: cl = 0.20d0
!real*8,parameter,private                        :: k1 = 0.52d0
!real*8,parameter,private                        :: k2 = 0.79d0
! parameters borrowed from Wu and Chou
real*8,parameter,private                        :: cd = 0.36d0
real*8,parameter,private                        :: cm = 1.00d0
real*8,parameter,private                        :: cl = 0.36d0
real*8,parameter,private                        :: k1 = 0.52d0
real*8,parameter,private                        :: k2 = 0.79d0


public bedroughness_init
public bedroughness_update
private update_urms
private infiltration_boundary_layer_effect
private acceleration_boundary_layer_effect
private acceleration_boundary_layer_effect_Nielsen
private turbulence_boundary_layer_effect

contains

subroutine bedroughness_init(s,par)
   use params
   use spaceparams
   
   IMPLICIT NONE

   type(parameters),intent(in)                 :: par
   type(spacepars)                             :: s
   integer                                     :: i,j

   include 's.ind'
   include 's.inp'
      
   select case (trim(par%bedfriction))
      
   case('chezy')
      cfu = par%g/s%bedfriccoef**2
      cfv = par%g/s%bedfriccoef**2
   case('cf')
      cfu = s%bedfriccoef
      cfv = s%bedfriccoef
   case ('manning')
      where(wetu==1)
         cfu = par%g*s%bedfriccoef**2/hu**(1.d0/3)
         cfu = min(cfu,0.1d0)
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1)
         cfv = par%g*s%bedfriccoef**2/hv**(1.d0/3)
         cfv= min(cfv,0.1d0)
      elsewhere
         cfv = 0.1d0
      endwhere
   case ('white-colebrook-grainsize')
      allocate (kru(nx+1,ny+1))
      allocate (krv(nx+1,ny+1))
      ! compute initial D90 at cell interface
      if (par%ngd==1) then
         kru = 3.d0*D90top
         krv = 3.d0*D90top
      else
         ! Rougness height in u-points
         do j=1,ny+1
            do i=1,nx
               kru(i,j) = 1.5d0*(D90top(i,j)+D90top(i+1,j))
            enddo
         enddo
         kru(nx+1,:) = kru(nx,:)
         ! Roughness height in v-points
         if (par%ny==0) then
            ! v-point is central on top of z-point
            krv = 3*D90top
         else             
            do j=1,ny
               do i=1,nx+1
                  krv(i,j) = 1.5d0*(D90top(i,j)+D90top(i,j+1))
               enddo
            enddo
            krv(:,ny+1) = krv(:,ny)
         endif
      endif
      where(wetu==1) 
         cfu = par%g/(18*log10(12*max(hu,kru)/kru))**2
         cfu = min(cfu,0.1d0)
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1) 
         cfv = par%g/(18*log10(12*max(hv,krv)/krv))**2
         cfv = min(cfv,0.1d0)
      elsewhere
         cfv = 0.1d0
      endwhere
   case ('white-colebrook')
      allocate (kru(nx+1,ny+1))
      allocate (krv(nx+1,ny+1))
      kru = s%bedfriccoef
      krv = s%bedfriccoef
      where(wetu==1) 
         cfu = par%g/(18*log10(12*max(hu,kru)/kru))**2
         cfu = min(cfu,0.1d0)
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1) 
         cfv = par%g/(18*log10(12*max(hv,krv)/krv))**2
         cfv = min(cfv,0.1d0)
      elsewhere
         cfv = 0.1d0
      endwhere
   case ('wilson')
      delta = (par%rhos-par%rho)/par%rho
      cfu = 0.5d0*0.114d0*(0.1d0/(par%g*delta*10.d0**2))**0.4d0
   case('wilson-hughes')
      allocate (kru(nx+1,ny+1))
      allocate (krv(nx+1,ny+1))
      allocate (shieldsu(nx+1,ny+1))
      allocate (shieldsv(nx+1,ny+1))
      delta = (par%rhos-par%rho)/par%rho
      ! compute initial D90 at cell interface
      if (par%ngd==1) then
         kru = 3.d0*D90top
         krv = 3.d0*D90top
      else
         ! Rougness height in u-points
         do j=1,ny+1
            do i=1,nx
               kru(i,j) = 1.5d0*(D90top(i,j)+D90top(i+1,j))
            enddo
         enddo
         kru(nx+1,:) = kru(nx,:)
         ! Roughness height in v-points
         if (par%ny==0) then
            ! v-point is central on top of z-point
            krv = 3*D90top
         else             
            do j=1,ny
               do i=1,nx+1
                  krv(i,j) = 1.5d0*(D90top(i,j)+D90top(i,j+1))
               enddo
            enddo
            krv(:,ny+1) = krv(:,ny)
         endif
      endif
   case('puleo-holland')
      delta = (par%rhos-par%rho)/par%rho
      rhogdelta = par%rho*par%g*delta
      allocate (kru90(nx+1,ny+1))
      allocate (krv90(nx+1,ny+1))
      allocate (kru50(nx+1,ny+1))
      allocate (krv50(nx+1,ny+1))
      allocate (shieldsu(nx+1,ny+1))
      allocate (shieldsv(nx+1,ny+1))
      ! compute initial D90 at cell interface
      if (par%ngd==1) then
         kru90 = D90top
         krv90 = D90top
         kru50 = D50top
         krv50 = D50top
      else
         ! Rougness height in u-points
         do j=1,ny+1
            do i=1,nx
               kru90(i,j) = 0.5d0*(D90top(i,j)+D90top(i+1,j))
               kru50(i,j) = 0.5d0*(D50top(i,j)+D50top(i+1,j))
            enddo
         enddo
         kru90(nx+1,:) = kru90(nx,:)
         kru50(nx+1,:) = kru50(nx,:)
         ! Roughness height in v-points
         if (par%ny==0) then
            ! v-point is central on top of z-point
            krv90 = D90top
            krv50 = D50top
         else             
            do j=1,ny
               do i=1,nx+1
                  krv90(i,j) = 0.5d0*(D90top(i,j)+D90top(i,j+1))
                  krv50(i,j) = 0.5d0*(D50top(i,j)+D50top(i,j+1))
               enddo
            enddo
            krv90(:,ny+1) = krv90(:,ny)
            krv50(:,ny+1) = krv50(:,ny)
         endif
      endif
      where(wetu==1) 
         cfu = 1.d0/(2.5d0*max(log(10*hu/kru90),1.d-12))**2  ! note from f = 2/etc. to cf = 1/etc. with cf = 1/2*f
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1) 
         cfv = 1.d0/(2.5d0*max(log(10*hv/krv90),1.d-12))**2  ! note from f = 2/etc. to cf = 1/etc. with cf = 1/2*f
      elsewhere
         cfu = 0.1d0
      endwhere
   case('swart')
      cfu = 0.5d0*exp(5.213d0*(2.5d0/0.1d0*D50top)**(0.194)-5.977)
   case ('vangent-okayasu')
      delta = (par%rhos-par%rho)/par%rho
      allocate (shieldsu(nx+1,ny+1))
      allocate (shieldsv(nx+1,ny+1))
      cfu = 0.1d0
      cfv = 0.1d0
   end select



end subroutine bedroughness_init


subroutine bedroughness_update(s,par)

   use params
   use spaceparams
   
   IMPLICIT NONE

   type(parameters),intent(in)                 :: par
   type(spacepars)                             :: s
   integer                                     :: i,j
   real*8                                      :: ubed,dudtbed,hda,klocal
   
   include 's.ind'
   include 's.inp'
   
   select case (trim(par%bedfriction))
      
   case('chezy')
      ! do nothing, this is constant bed roughness, already set in initialise
   case('cf')
      ! do nothing, this is constant bed roughness, already set in initialise
   case ('manning')
      ! C = H**(1/6)/n
      ! cf = g/C**2 = g/(hu**(1/6)/n)**2 = g*n**2/hu**(1/3)
      where(wetu==1) 
         cfu = par%g*s%bedfriccoef**2/hu**(1.d0/3)
         cfu = min(cfu,0.1d0)
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1) 
         cfv = par%g*s%bedfriccoef**2/hv**(1.d0/3)
         cfv = min(cfv,0.1d0)
      elsewhere
         cfv = 0.1d0
      endwhere
   case ('white-colebrook-grainsize')
      if (par%ngd>1) then
         ! Rougness height in u-points
         do j=1,ny+1
            do i=1,nx
               if (wetu(i,j)==1) then
                  kru(i,j) = 1.5d0*(D90top(i,j)+D90top(i+1,j))
               endif
            enddo
         enddo
         kru(nx+1,:) = kru(nx,:)
         ! Roughness height in v-points
         if (par%ny==0) then
            ! v-point is central on top of z-point
            where (wetv==1)
               krv = 3*D90top
            endwhere
         else             
            do j=1,ny
               do i=1,nx+1
                  if (wetv(i,j)==1) then
                     krv(i,j) = 1.5d0*(D90top(i,j)+D90top(i,j+1))
                  endif
               enddo
            enddo
            krv(:,ny+1) = krv(:,ny)
         endif
      endif
      where(wetu==1) 
         cfu = par%g/(18*log10(12*max(hu/kru,2.d0)))**2
         cfu = min(cfu,0.1d0)
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1) 
         cfv = par%g/(18*log10(12*max(hv/krv,2.d0)))**2
         cfv = min(cfv,0.1d0)
      elsewhere
         cfv = 0.1d0
      endwhere
   case ('white-colebrook')
      where(wetu==1) 
         cfu = par%g/(18*log10(12*max(hu/kru,2.d0)))**2
         cfu = min(cfu,0.1d0)
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1) 
         cfv = par%g/(18*log10(12*max(hv/krv,2.d0)))**2
         cfv = min(cfv,0.1d0)
      elsewhere
         cfv = 0.1d0
      endwhere
   case ('wilson')
      call update_urms(s,par)
      where(wetu==1) 
         cfu = 0.5d0*0.114d0*(max(urms_upd,0.1d0)/(par%g*delta*par%Trep**2))**0.4d0
      elsewhere
         cfu = 0.1d0
      endwhere
   case ('wilson-hughes')
      if (par%ngd>1) then
         ! Rougness height in u-points
         do j=1,ny+1
            do i=1,nx
               if (wetu(i,j)==1) then
                  kru(i,j) = 1.5d0*(D90top(i,j)+D90top(i+1,j))
               endif
            enddo
         enddo
         kru(nx+1,:) = kru(nx,:)
         ! Roughness height in v-points
         if (par%ny==0) then
            ! v-point is central on top of z-point
            where (wetv==1)
               krv = 3*D90top
            endwhere
         else             
            do j=1,ny
               do i=1,nx+1
                  if (wetv(i,j)==1) then
                     krv(i,j) = 1.5d0*(D90top(i,j)+D90top(i,j+1))
                  endif
               enddo
            enddo
            krv(:,ny+1) = krv(:,ny)
         endif
      endif
      where(wetu==1)
         ! first estimate f in eq 13 of Hughes 1995. Note that taubx = 1/8*rho*f*U^2 in Wilson, 
         ! therefore cf = 1/8f
         cfu = 1.d0/(2.5d0*log(30*max(hu/kru,1.d0)))**2
         ! now estimate shields using eq 19 of Hughes 1995. Note, assuming approximation of shoreline
         ! velocity according to eq 20
         shieldsu = cfu*vmageu**2/par%g/D50top/delta
         ! finally, update cf
         cfu = 1.d0/(2.5d0*log(0.5d0*max(hu/shieldsu/D50top,1.d0)))**2
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1)
         ! first estimate f in eq 13 of Hughes 1995. Note that taubx = 1/8*rho*f*U^2 in Wilson, 
         ! therefore cf = 1/8f
         cfv = 1.d0/(2.5d0*log(30*max(hv/krv,1.d0)))**2
         ! now estimate shields using eq 19 of Hughes 1995. Note, assuming approximation of shoreline
         ! velocity according to eq 20
         shieldsv = cfu*vmagev**2/par%g/D50top/delta
         ! finally, update cf
         cfu = 1.d0/(2.5d0*log(0.5d0*max(hv/shieldsv/D50top,1.d0)))**2
      elsewhere
         cfu = 0.1d0
      endwhere
      where(wetv==1)
         ! first estimate f in eq 13 of Hughes 1995. Note that taubx = 1/8*rho*f*U^2 in Wilson, 
         ! therefore cf = 1/8f
         cfv = 1.d0/(2.5d0*log(30*max(hv/krv,1.d0)))**2
         ! now estimate shields using eq 19 of Hughes 1995. Note, assuming approximation of shoreline
         ! velocity according to eq 20
         shieldsv = cfv*vmagev**2/par%g/D50top/delta
         ! finally, update cf
         cfv = 1.d0/(2.5d0*log(0.5d0*max(hv/shieldsv/D50top,1.d0)))**2
      elsewhere
         cfv = 0.1d0
      endwhere   
   case ('puleo-holland')
      ! Update sediment properties in u and v points
      if (par%ngd>1) then
         ! Rougness height in u-points
         do j=1,ny+1
            do i=1,nx
               if (wetu(i,j)==1) then
                  kru90(i,j) = 0.5d0*(D90top(i,j)+D90top(i+1,j))
                  kru50(i,j) = 0.5d0*(D50top(i,j)+D50top(i+1,j))
               endif
            enddo
         enddo
         kru90(nx+1,:) = kru90(nx,:)
         kru50(nx+1,:) = kru50(nx,:)
         ! Roughness height in v-points
         if (par%ny==0) then
            ! v-point is central on top of z-point
            where (wetv==1)
               krv90 = D90top
               krv50 = D50top
            endwhere
         else             
            do j=1,ny
               do i=1,nx+1
                  if (wetv(i,j)==1) then
                     krv90(i,j) = 0.5d0*(D90top(i,j)+D90top(i,j+1))
                     krv50(i,j) = 0.5d0*(D50top(i,j)+D50top(i,j+1))
                  endif
               enddo
            enddo
            krv90(:,ny+1) = krv90(:,ny)
            krv50(:,ny+1) = krv50(:,ny)
         endif
      endif
      ! check Shields values in u and v points
      where (wetu==1)
         shieldsu = abs(taubx)/(rhogdelta*kru50)
      endwhere
      where (wetv==1)
         shieldsv = abs(tauby)/(rhogdelta*krv50)
      endwhere
      ! compute bed friction coefficient based on shields value
      where (wetu==1 .and. shieldsu > 0.8d0)  ! sediment-laden flow
         cfu = max(   1.d0/(2.5d0*log(30    *max(2.d0,hu/(3*kru90)          )))**2, &
                      1.d0/(2.5d0*log(53.2d0*max(2.d0,hu/(10*shieldsu*kru50))))**2  &
                   )
         !
         !
         !cfu = max(1.d0/(2.5d0*max(log(5.32d0*hu*max(rhogdelta/abs(taubx),10*1/hu)),1.d-12))**2, &   ! sediment-laden
         !          1.d0/(2.5d0*log(10*max(hu,kru90)/kru90))**2)                          ! clear water
      elsewhere (wetu==1 .and. shieldsu <= 0.8d0)  ! clear water flow
        cfu = 1.d0/(2.5d0*log(30    *max(2.d0,hu/(3*kru90)          )))**2   
       !cfu = 1.d0/(2.5d0*log(10*max(hu,kru90)/kru90))**2
      elsewhere
         cfu = 0.1d0
      endwhere
      where (wetv==1 .and. shieldsv > 0.8d0)  ! sediment-laden flow
         cfv = max(   1.d0/(2.5d0*log(30    *max(2.d0,hv/(3*krv90)          )))**2, &
                      1.d0/(2.5d0*log(53.2d0*max(2.d0,hv/(10*shieldsv*krv50))))**2  &
                   )
         !cfv = max(1.d0/(2.5d0*max(log(5.32d0*hv*max(rhogdelta/abs(tauby),10*1/hv)),1.d-12))**2, &   ! sediment-laden
         !          1.d0/(2.5d0*log(10*max(hv,krv90)/krv90))**2)                          ! clear water
      elsewhere (wetv==1 .and. shieldsv <= 0.8d0)  ! clear water flow
         cfv = 1.d0/(2.5d0*log(30    *max(2.d0,hv/(3*krv90)          )))**2 
         !cfv = 1.d0/(2.5d0*log(10*max(hv,krv90)/krv90))**2
      elsewhere
         cfv = 0.1d0
      endwhere
   case('swart')
      call update_urms(s,par)
      where(wetu==1) 
         cfu = 0.5d0*exp(5.213d0*(2.5d0/urms_upd*D50top)**(0.194)-5.977)
      elsewhere
         cfu = 0.1d0
      endwhere
   case ('vangent-okayasu')
      ! compute all components of Shields 
      ! Shields = (Fd + Fi)/(W-Fl) = tbx/(delta*rho*g*D50) -> tbx = Shields*delta*rho*g*D50
      ! tbx = cf*rho*u^2 -> cf = tbx/(rho*u^2)
      ! cf = (Shields*delta*g*D50)/(u^2)
      !
      do j=1,ny+1
         do i=1,nx
            if (wetu(i,j)==1) then
               !
               ! http://en.wikipedia.org/wiki/Law_of_the_wall
               hda = max(0.37d0*hh(i,j),par%eps)
               klocal = min(max(3*D90top(i,j),5*shieldsu(i,j)*D50top(i,j)),hh(i,j))
               !ubed = 0.41d0/log(hda/min(par%eps,klocal/30.d0))*uu(i,j)
               ubed = uu(i,j)/(2.5d0*log(11.0d0*max(hh(i,j)/klocal,1.d0)))
               ! 
               dudtbed = dudt(i,j)*abs(ubed/uu(i,j))
               Fd = 0.5d0*par%rho*cd*k2*D50top(i,j)**2*ubed*abs(ubed)
               Fi = 0*par%rho*cm*k1*D50top(i,j)**3*dudtbed
               Fl = 0.5d0*par%rho*cl*k2*D50top(i,j)**2*ubed**2
               Wg = (par%rhos-par%rho)*par%g*k1*D50top(i,j)**3
               shieldsu(i,j) = abs(Fd+Fi)/max(Wg-Fl,1.d-8)
               cfu(i,j) = min(shieldsu(i,j)*delta*par%g*D50top(i,j)/uu(i,j)**2,0.1d0)
               !cfu(i,j) = cfu(i,j)*0.5d0 ! try this out to see how implicit infilration may affect bed roughness results
            else
               cfu(i,j) = 0.1d0
            endif
         enddo
      enddo
                  
   end select
   
   ! Additional bed friction terms due to acceleration
   taubx_add = 0.d0
   tauby_add = 0.d0
   
   ! Below: Nielsen not preferred
   if (par%friction_nielsen==1) then
      call acceleration_boundary_layer_effect_Nielsen(s,par)
   else
      ! turbulence and acceleration add to shear stress, but do
      ! not necessarily also need enhanced cf from infiltration
      if (par%friction_turbulence==1) then
         call turbulence_boundary_layer_effect(s,par)   
      endif
      
      if (par%friction_acceleration==1) then
         call acceleration_boundary_layer_effect(s,par)
      endif
      
      ! Infiltration enhances normal bed shear stress through 
      ! increase of cfu and cfv
      if (par%friction_infiltration==1) then
         call infiltration_boundary_layer_effect(s,par)
      endif
   
      cfu = min(cfu,1.d0)
      cfv = min(cfv,1.d0)
   endif
   
   
   
end subroutine bedroughness_update 
   
subroutine update_urms(s,par)

   use params
   use spaceparams
   
   IMPLICIT NONE

   type(parameters),intent(in)                 :: par
   type(spacepars)                             :: s
   real*8                                      :: Tref
   
   include 's.ind'
   include 's.inp'
   
   if(.not. allocated(urms_upd)) then 
      allocate(urms_upd(nx+1,ny+1))
      allocate(u2_upd(nx+1,ny+1))
      urms_upd = 0.d0
      if (par%swave==1) then 
         u2_upd = sqrt(uu**2 + urms**2)
      else
         u2_upd = uu
      endif
   endif
      
   Tref = 15*par%Trep
   
   if (par%swave==1) then
      where(wetu==1)
         u2_upd = (1.d0-par%dt/Tref)*u2_upd + (par%dt/Tref)*(uu**2 + urms**2)
      endwhere
   else
      where(wetu==1)
         u2_upd = (1.d0-par%dt/Tref)*u2_upd + (par%dt/Tref)*uu**2 
      endwhere
   endif
   
   urms_upd = sqrt(u2_upd)
   
end subroutine update_urms

subroutine infiltration_boundary_layer_effect(s,par)

   use params
   use spaceparams
   
   IMPLICIT NONE

   type(parameters),intent(in)                 :: par
   type(spacepars)                             :: s
   integer                                     :: i,j
   real*8,parameter                            :: epsVentilation = 1.d0
   real*8,parameter                            :: maxEnhancement1 = 3.d0 ! Pers. Corr. Conley 2014
   real*8,parameter                            :: maxEnhancement2 = 0.1d0 ! Pers. Corr. Conley 2014
   
   include 's.ind'
   include 's.inp'
   
   if (.not.allocated(facbl)) then
      allocate (facbl(nx+1,ny+1))
      allocate (blphi(nx+1,ny+1))
      allocate (infilb(nx+1,ny+1))
      allocate (Ubed(nx+1,ny+1))
      allocate (Ventilation(nx+1,ny+1))
   endif
   
   
   ! Include infiltration effect on boundary layer thinning
   ! Taken from eq 13 Butt, Russell, Turner 2001, based on
   ! Conley and Inman 1994
   ! In C&I derivation f = 2*cf and w is positive up, and b=0.9-2.0 so modified to:
   ! Tau_w/Tau_0 = C/(2*cf) * -infil/abs(u)
   ! where C is a constant of ~2
   !
   ! Tau ~ cf
   ! cf_w/cf_0 = Tau_w/Tau_0

   ! U-points
   do j=1,ny+1
      do i=1,nx
         if(wetu(i,j)==1) then
            infilb(i,j) = 0.5d0*(infil(i,j)+infil(i+1,j))
         else
            infilb(i,j) = 0.d0
         endif
      enddo
   enddo
   infilb(nx+1,:) = infilb(nx,:)
   where(wetu==1)
      !Ubed = sqrt(cfu)*abs(vmageu)
      Ubed = abs(uu)
      Ventilation  = -infilb/max(Ubed,1.d-6)
      Ventilation  = max(min(Ventilation,epsVentilation),-epsVentilation)
      blphi = 0.9d0/2/cfu*Ventilation
      ! need to account for devision by zero
      where(infilb>=0.d0)
         blphi = min(blphi,-1d-4)
         facbl = blphi/(exp(blphi)-1.d0)
         facbl = min(facbl,maxEnhancement1) 
      elsewhere
         blphi = max(blphi,1d-4)
         facbl = blphi/(exp(blphi)-1.d0)
         facbl = max(facbl,maxEnhancement2)
      endwhere
      cfu = cfu*facbl
      cfu = min(cfu,1.d0)
   endwhere
   ! V-points
   if(par%ny>0) then
      do j=1,ny
         do i=1,nx+1
            if(wetv(i,j)==1) then
               infilb(i,j) = 0.5d0*(infil(i,j)+infil(i,j+1))
            else
               infilb(i,j) = 0.d0
            endif
         enddo
      enddo
      infilb(:,ny+1) = infil(:,ny)
   else
      ! v-points centered on cell centre
      infilb = infil
   endif
   where(wetv==1)
      !Ubed = sqrt(cfv)*abs(vmagev)
      Ubed = abs(vv)
      Ventilation  = -infilb/max(Ubed,1.d-6)
      Ventilation  = max(min(Ventilation,epsVentilation),-epsVentilation)
      blphi = 0.9d0/2/cfv*Ventilation
      ! need to account for devision by zero
      where(infilb>=0.d0)
         blphi = min(blphi,-1d-4)
         facbl = blphi/(exp(blphi)-1.d0)
         facbl = min(facbl,maxEnhancement1)
      elsewhere
         blphi = max(blphi,1d-4)
         facbl = blphi/(exp(blphi)-1.d0)
         facbl = max(facbl,maxEnhancement2)
      endwhere
      cfv = cfv*facbl
      cfv = min(cfv,1.d0)
   endwhere

end subroutine infiltration_boundary_layer_effect

subroutine acceleration_boundary_layer_effect(s,par)

   use params
   use spaceparams
   
   IMPLICIT NONE

   type(parameters),intent(in)                 :: par
   type(spacepars)                             :: s
   
!   real*8,parameter                            :: k1 = 0.5236d0 ! pi/6
   real*8                                      :: Tsmooth,factime
   real*8,dimension(:,:),allocatable,save      :: Fi
!   real*8,dimension(:,:),allocatable,save      :: dynpresupd
!   integer                                     :: i,j
!   integer,dimension(:,:),allocatable,save     :: wetzp1,wetu_ududx
  
   include 's.ind'
   include 's.inp'
   
   if(.not.allocated(dudtsmooth)) then
      allocate (dudtsmooth(nx+1,ny+1))
      allocate (dvdtsmooth(nx+1,ny+1))
      allocate (ueuold    (nx+1,ny+1))
      allocate (uevold    (nx+1,ny+1))
      allocate (vevold    (nx+1,ny+1))
      allocate (veuold    (nx+1,ny+1))
      allocate (ueuf      (nx+1,ny+1))
      allocate (uevf      (nx+1,ny+1))
      allocate (vevf      (nx+1,ny+1))
      allocate (veuf      (nx+1,ny+1))
      allocate (Fi        (nx+1,ny+1))
!      allocate (wetzp1    (nx+1,ny+1))
!      allocate (wetu_ududx(nx+1,ny+1))
      ueuold = 0.d0
      uevold = 0.d0
      vevold = 0.d0
      veuold = 0.d0
      ueuf   = 0.d0
      uevf   = 0.d0
      vevf   = 0.d0
      veuf   = 0.d0
      dtold = par%dt
!      allocate(dynpresupd(1:nx+1,1:ny+1))
!      dynpresupd = 0.d0
   endif
   
   ! Problem with derivation of DUdt in U and V points
   ! ueu and veu are not known until end of flow timestep, so
   ! both are based on previous time step. The temporal gradient
   ! can only be found by going back one more time step.
   ! We need to keep track of par%dt so the correct time step for
   ! the gradient is used:
   !
   !    i-2            i-1           i (now)          i+1 (future)
   !         dtold(1)       dtold(2)          par%dt   
   !   ueuold          ueuf         unknown
   !
   ! Note that this assumes that the acceleration term does not change significantly
   ! over 2 timesteps and that the location mask wetu can still be used,
   ! even though it is not the wetu mask of 2 timesteps ago.
   Tsmooth = par%Trep/20
   !Tsmooth = min(par%Trep/100,0.0333d0)
   factime = min(dtold(1)/Tsmooth,1.d0)
   !factime = 1.d0
   !wetzp1(1:nx,:) = wetz(2:nx+1,:)
   !wetzp1(nx+1,:) = wetz(nx,:)
   !wetu_ududx = min(wetz,wetzp1)
   where(wetu==1)
      ueuf = (1-factime)*ueuf + factime*ueu
      dudtsmooth = (ueuf-ueuold)/dtold(1)
      !dudtsmooth = (ueuf-ueuold)/dtold(1)+ududx*wetu_ududx
      where (dudtsmooth>100*par%g)
         dudtsmooth = 100*par%g
      elsewhere (dudtsmooth<-100*par%g)
         dudtsmooth = -100*par%g
      endwhere
      !veuf = (1-factime)*veuf + factime*veu
      !dudtsmooth = sqrt(((ueuf-ueuold)/dtold(1))**2 + &
      !                  ((veuf-veuold)/dtold(1))**2)
   elsewhere
      ueuf = 0.d0
      veuf = 0.d0
      dudtsmooth = 0.d0
   endwhere
   where(wetv==1)
      vevf = (1-factime)*vevf + factime*vev
      dvdtsmooth = (vevf-vevold)/dtold(1)
      !uevf = (1-factime)*uevf + factime*uev
      !dvdtsmooth = sqrt(((vevf-vevold)/dtold(1))**2 + &
      !                  ((uevf-uevold)/dtold(1))**2)
      where (dvdtsmooth>100*par%g)
         dvdtsmooth = 100*par%g
      elsewhere (dvdtsmooth<-100*par%g)
         dvdtsmooth = -100*par%g
      endwhere
   elsewhere
      vevf = 0.d0
      uevf = 0.d0
      dvdtsmooth = 0.d0
   endwhere
   ! now set the old variable to the i-1 variables
   ueuold = ueuf
   uevold = uevf
   vevold = vevf
   veuold = veuf
   dtold(1) = dtold(2)
   dtold(2) = par%dt 
   !
   !
   ! Van Gent: Fi = rho*cm*k1*D^3*du/dt. [N] This is force on one stone, not per m^2
   ! Assume 1/D^2 stones per m^2
   ! Robert :  Fi = rho*cm*k1*D*du/dt. [N/m^2]
   !
   ! Fupd = F0 + Fi
   ! F ~ cf, therefore
   ! cfupd = cf0 + cfi
   ! 
   ! cfi = Fi / (par%rho*u^2).
   
   !factime = min(4*par%dt/par%Trep,1.d0)
   !dynpresupd = (1-factime) * dynpresupd + factime * pres
   !do j=1,s%ny+1
   !   do i=1,s%nx
   !      if(wetz(i,j)==1 .and. wetz(i+1,j)==1) then
   !         dudtsmooth(i,j) = -(par%g*(zs(i+1,j)-zs(i,j))+(dynpresupd(i+1,j)-dynpresupd(i,j)))/dsu(i,j)
   !      else
   !         dudtsmooth(i,j) = 0.d0
   !      endif
   !   enddo
   !enddo
   !where (dudtsmooth>100*par%g)
   !   dudtsmooth = 100*par%g
   !elsewhere (dudtsmooth<-100*par%g)
   !   dudtsmooth = -100*par%g
   !endwhere
      
   where(wetu==1)
      taubx_add = taubx_add + par%ci*par%rho*min(D50top,hu)*dudtsmooth
      !taubx_add = taubx_add + par%cm*par%rho*k1*min(D50top,hu*par%por)*dudtsmooth
   endwhere
   where(wetv==1)
      tauby_add = tauby_add + par%ci*par%rho*min(D50top,hv)*dvdtsmooth 
   endwhere
   
end subroutine acceleration_boundary_layer_effect

subroutine acceleration_boundary_layer_effect_Nielsen(s,par)

   use params
   use spaceparams
   
   IMPLICIT NONE

   type(parameters),intent(in)                 :: par
   type(spacepars)                             :: s
   real*8                                      :: omegap,Tsmooth,factime
   real*8,save                                 :: phirad
   
   include 's.ind'
   include 's.inp'
   
   if(.not.allocated(dudtsmooth)) then
      allocate (dudtsmooth(nx+1,ny+1))
      allocate (dvdtsmooth(nx+1,ny+1))
      allocate (ueuold    (nx+1,ny+1))
      allocate (uevold    (nx+1,ny+1))
      allocate (vevold    (nx+1,ny+1))
      allocate (veuold    (nx+1,ny+1))
      allocate (ueuf      (nx+1,ny+1))
      allocate (uevf      (nx+1,ny+1))
      allocate (vevf      (nx+1,ny+1))
      allocate (veuf      (nx+1,ny+1))
      ueuold = 0.d0
      uevold = 0.d0
      vevold = 0.d0
      veuold = 0.d0
      ueuf   = 0.d0
      uevf   = 0.d0
      vevf   = 0.d0
      veuf   = 0.d0
      dtold = par%dt
      phirad  = par%phit/180*par%px
   endif
   
   omegap = 2*par%px/par%Trep
   ! Problem with derivation of DUdt in U and V points
   ! ueu and veu are not known until end of flow timestep, so
   ! both are based on previous time step. The temporal gradient
   ! can only be found by going back one more time step.
   ! We need to keep track of par%dt so the correct time step for
   ! the gradient is used:
   !
   !    i-2            i-1           i (now)          i+1 (future)
   !         dtold(1)       dtold(2)          par%dt   
   !   ueuold          ueuf         unknown
   !
   ! Note that this assumes that the acceleration term does not change significantly
   ! over 2 timesteps and that the location mask wetu can still be used,
   ! even though it is not the wetu mask of 2 timesteps ago.
   Tsmooth = par%Trep/20
   factime = min(dtold(1)/Tsmooth,1.d0)
   where(wetu==1)
      ueuf = (1-factime)*ueuf + factime*ueu
      veuf = (1-factime)*veuf + factime*veu
      dudtsmooth = sqrt(((ueuf-ueuold)/dtold(1))**2 + &
                        ((veuf-veuold)/dtold(1))**2)
   elsewhere
      ueuf = 0.d0
      veuf = 0.d0
      dudtsmooth = 0.d0
   endwhere
   where(wetv==1)
      vevf = (1-factime)*vevf + factime*vev
      uevf = (1-factime)*uevf + factime*uev
      dvdtsmooth = sqrt(((vevf-vevold)/dtold(1))**2 + &
                        ((uevf-uevold)/dtold(1))**2)
   elsewhere
      vevf = 0.d0
      uevf = 0.d0
      dvdtsmooth = 0.d0
   endwhere
   ! now set the old variable to the i-1 variables
   ueuold = ueuf
   uevold = uevf
   vevold = vevf
   veuold = veuf
   dtold(1) = dtold(2)
   dtold(2) = par%dt
   !
   ! The acceleration term is written for ustar, but can be rewritten
   ! to affect cf:
   ! Nielsen, including acceleration effects:
   ! ustar = sqrt(cf) * (cos(phi)*ue + 1/w*sin(phi)*dudt)
   ! should be rewritten in form of:
   ! ustar = sqrt(cfw) * ue
   ! 
   ! sqrt(cfw) = sqrt(cf)*(cos(phi)*ue + 1/w*sin(phi)*dudt) / ue
   ! 
   where(wetu==1)
      cfu = cfu*(cos(phirad)*ueu + &
               sin(phirad)/omegap*dudtsmooth)**2 / &
            max(ueu**2,1.d-12)
   endwhere
   where(wetv==1)
      cfv = cfv*(cos(phirad)*vev + &
               sin(phirad)/omegap*dvdtsmooth)**2 / &
            max(vev**2,1.d-12)
   endwhere
   

end subroutine acceleration_boundary_layer_effect_Nielsen

subroutine turbulence_boundary_layer_effect(s,par)

   use params
   use spaceparams
   
   IMPLICIT NONE

   type(parameters),intent(in)                 :: par
   type(spacepars)                             :: s
   real*8,dimension(:,:),allocatable,save      :: kbl
   integer                                     :: i,j
   
   include 's.ind'
   include 's.inp'
   
   ! Following Reniers et al (2004)
   ! ubed_turb = sqrt(ubed^2+gamma*kb) :: gamma ~ 1
   ! u- = u without turbulence
   ! u+ = u with turbulence
   ! 
   ! New way: compute taubx_add and tauby_add
   ! taubx_total = cf*rho*(ubed^2+gamma*kb)
   ! taubx = cf*rho*(ubed^2)
   ! taubx_add = cf*rho*(gamma*kb)
   
   
   ! Old way: increase cfu and cfv
   !!!! u+ = sqrt(u-^2 + kb)
   !!!! u+^2 = u-^2 + kb
   !!!! taubx- ~ u-^2, taubx+ ~ u+^2
   !!!! factor = taubx+/taubx- = u+^2/u-^2 = (u-^2 + kb)/u-^2 = 1 + kb/u-^2
   
    if (.not.allocated(kbl)) then
      allocate (kbl(nx+1,ny+1))
    endif
    
    ! U-points
   do j=1,ny+1
      do i=1,nx
         if(wetu(i,j)==1) then
            kbl(i,j) = 0.5d0*(kb(i,j)+kb(i+1,j))
         else
            kbl(i,j) = 0.d0
         endif
      enddo
   enddo
   kbl(nx+1,:) = kbl(nx,:)
   where(wetu==1 .and. kbl>0.d0)
      ! New way:
      taubx_add = taubx_add + cfu*par%rho*kbl*sign(1.d0,ue)
      
      ! Old way:
      !cfu = cfu*(1.d0+(kbl/max(uu**2,1.d-6)))
      !cfu = cfu*(1.d0+(2.d0**(1.25d0)*sqrt(kbl))/max(abs(uu),1.d-6)+(sqrt(2.d0)*kbl)/max(uu**2,1.d-6))
   endwhere
   ! V-points
   if(par%ny>0) then
      do j=1,ny
         do i=1,nx+1
            if(wetv(i,j)==1) then
               kbl(i,j) = 0.5d0*(kb(i,j)+kb(i,j+1))
            else
               kbl(i,j) = 0.d0
            endif
         enddo
      enddo
      kbl(:,ny+1) = kbl(:,ny)
   else
      ! v-points centered on cell centre
      kbl = kb
   endif
   where(wetv==1 .and. kbl>0.d0)
       ! New way:
      tauby_add = tauby_add + cfv*par%rho*kbl*sign(1.d0,ve)
      
      ! Old way:
      !cfv = cfv*(1.d0+(kbl/max(vv**2,1.d-6)))
      !cfv = cfv*(1.d0+(2.d0**(1.25d0)*sqrt(kbl))/max(abs(vv),1.d-6)+(sqrt(2.d0)*kbl)/max(vv**2,1.d-6))
   endwhere
   
   


end subroutine turbulence_boundary_layer_effect


end module bedroughness_module