!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
! 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 turbulence_module
   private    ! set private default
   ! set these public to users of groundwater module
   public waveturb
  
contains
   
 subroutine waveturb(s,par)
    use params
    use spaceparams
    use xmpi_module

    IMPLICIT NONE

    type(spacepars),target                   :: s
    type(parameters)                         :: par

    integer                                  :: i
    integer                                  :: j
    real*8, save                             :: twothird
    real*8,dimension(:,:),allocatable,save   :: ksource, kturbu,kturbv,Sturbu,Sturbv,dzsdt_cr,ML,dcfin,dcf

    include 's.ind'
    include 's.inp'

    if (.not. allocated(kturbu)) then
       allocate(ksource (nx+1,ny+1))
       allocate(kturbu (nx+1,ny+1))
       allocate(kturbv (nx+1,ny+1))
       allocate(Sturbu (nx+1,ny+1))
       allocate(Sturbv (nx+1,ny+1))
       allocate(dzsdt_cr (nx+1,ny+1))
       allocate(ML (nx+1,ny+1))
       allocate(dcfin (nx+1,ny+1))
       allocate(dcf (nx+1,ny+1))
       twothird=2.d0/3.d0
       Sturbu = 0.d0
       Sturbv = 0.d0
    endif
    ! use lagrangian velocities
    kturbu       = 0.0d0  !Jaap
    kturbv       = 0.0d0  !Jaap		
    ! Update roller thickness
    ! This is Jaap's code:
         !dzsdt_cr=par%beta*c
         !rolthick=rolthick+par%dt*(abs(dzsdt)-dzsdt_cr)
         !rolthick=max(rolthick,0.d0)
    ! This is what Robert changed
    dzsdt_cr = par%secbrsteep*sqrt(par%g*hh)
    rolthick = rolthick+par%dt*(max(dzsdt,0.d0)-dzsdt_cr)
    rolthick = min(max(rolthick,0.d0),hh)*wetz
    
    ! Jaap compute sources and sinks
    ! long wave source
    ksource = 0.d0
    if (par%lwt == 1) then
       ksource = par%yturb*par%g*rolthick*par%secbrsteep*sqrt(par%g*hh)
    endif
    if (trim(par%turbadv) == 'none') then  ! Robert: To Do: is this correct summation?
       kturb = ((DR/par%rho+ksource)/par%betad)**twothird          ! See Battjes, 1975 / 1985
    else
       ksource = ksource + DR/par%rho
     
    
       ! Jaap do long wave turb approach for both short waves and long waves
       !
       ! Turbulence in uu-points
       !
       do j=1,ny+1
         do i=1,nx
            if(wetu(i,j)==1) then
               if(uu(i,j)>0.d0) then
                 kturbu(i,j)=par%thetanum*kturb(i,j)+(1.d0-par%thetanum)*kturb(min(i+1,nx),j)
               elseif(uu(i,j)<0.d0) then
                 kturbu(i,j)=par%thetanum*kturb(i+1,j)+(1.d0-par%thetanum)*kturb(max(i,2),j)
               else
                 kturbu(i,j)=0.5d0*(kturb(i,j)+kturb(i+1,j))
               endif
            endif
         enddo
       enddo
       !do j=1,ny+1
       !  do i=1,nx
       !    if(uu(i,j)>=0.d0) then
       !      kturbu(i,j)=kturb(i,j)
       !    else
       !      kturbu(i,j)=kturb(i+1,j)
       !    endif
       !  enddo
       !enddo
       kturbu(nx+1,:) = kturb(nx+1,:) !Robert
    
       ! wwvv fix this in parallel case
#ifdef USEMPI
   call xmpi_shift(kturbu,'m:')
#endif
       !
       ! Turbulence in vv-points
       !
       do j=1,ny
         do i=1,nx+1
            if(wetv(i,j)==1) then
               if(vv(i,j)>0) then
                  kturbv(i,j)=par%thetanum*kturb(i,j)+(1.d0-par%thetanum)*kturb(i,min(j+1,ny))
               else if(vv(i,j)<0) then
                  kturbv(i,j)=par%thetanum*kturb(i,j+1)+(1.d0-par%thetanum)*kturb(i,max(j,2))
               else
                  kturbv(i,j)=0.5d0*(kturbv(i,j)+kturbv(i,j+1))
               endif
            endif
         enddo
       enddo
       kturbv(:,ny+1) = kturb(:,ny+1) !Robert
       !
       ! Turbulence advection in X and Y direction
       !
       if (trim(par%turbadv) == 'lagrangian') then
          Sturbu=kturbu*uu*hu*wetu
          Sturbv=kturbv*vv*hv*wetv
       elseif (trim(par%turbadv) == 'eulerian') then
          Sturbu=kturbu*ueu*hu*wetu
          Sturbv=kturbv*vev*hv*wetv
       endif
       !
       ! Update turbulence
       !
       if (ny>0) then
         do j=2,ny+1
           do i=2,nx+1
             if (wetz(i,j)==1) then
                kturb(i,j) = hold(i,j)*kturb(i,j)-par%dt*(       &
	                             (Sturbu(i,j)*dnu(i,j)-Sturbu(i-1,j)*dnu(i-1,j)+&
                                 Sturbv(i,j)*dsv(i,j)-Sturbv(i,j-1)*dsv(i,j-1))*dsdnzi(i,j)-&
                                (ksource(i,j)-par%betad*kturb(i,j)**1.5d0))
                kturb(i,j)=max(kturb(i,j),0.0d0)
             else
                kturb(i,j) = 0.d0
             endif  
            enddo
         enddo
       else
         j=1
         do i=2,nx+1
           if(wetz(i,j)==1) then 
              kturb(i,j) = hold(i,j)*kturb(i,j)-par%dt*(       &
	                           (Sturbu(i,j)-Sturbu(i-1,j))/dsz(i,j)-&
                              (ksource(i,j)-par%betad*kturb(i,j)**1.5d0))
              kturb(i,j)=max(kturb(i,j),0.0d0)
           else
              kturb(i,j) = 0.d0
           endif
         enddo
       endif
       
       kturb = kturb/max(hh,0.01d0)
       
       ! Jaap only required for advection mode?
       kturb(1,:)=kturb(2,:)
       kturb(nx+1,:)=kturb(nx+1-1,:)
       if (ny>0) then
         kturb(:,1)=kturb(:,2)
         kturb(:,ny+1)=kturb(:,ny+1-1)
       endif
       
       ! wwvv fix the first and last rows and columns of kturb in parallel case
#ifdef USEMPI
       call xmpi_shift(kturb,'1:')
       call xmpi_shift(kturb,'m:')
       call xmpi_shift(kturb,':1')   ! Maybe not necessary as zb calculated from 1:ny+1, but just in case...
       call xmpi_shift(kturb,':n')   ! Dito
#endif
    
    endif ! turbadv == 'none'
    !
    where(kturb>0.d0)
       ML = (rolthick*par%Trep/2*par%g*hh)   ! sqrt(area_of _roller) Van Thiel de Vries PhD
       ML = min(ML,hh);
       ! exponential decay turbulence over depth
       dcfin = exp(min(100.d0,hh/max(ML,0.01d0)))
       dcf = min(1.d0,1.d0/(dcfin-1.d0))
       kb = kturb*dcf
    elsewhere
       kb = 0.d0
    endwhere

   end subroutine waveturb
   
end module turbulence_module