module flow_timestep_module
contains
  subroutine flow(s,par)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! 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                                                                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    use params
    use spaceparams
    use xmpi_module
    use boundaryconditions

#ifndef USEMPI
    use flow_secondorder_module
    use nonh_module
#endif


    IMPLICIT NONE

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

    integer                                 :: i
    integer                                 :: j,j1,jp1
    real*8,dimension(:,:),allocatable,save  :: vv_old                   !Velocity at previous timestep
    real*8,dimension(:,:),allocatable,save  :: uu_old                   !Velocity at previous timestep
    real*8,dimension(:,:),allocatable,save  :: zs_old
    real*8,dimension(:,:),allocatable,save  :: vsu,usu,vsv,usv,veu,uev
    real*8,dimension(:,:),allocatable,save  :: ududx,vdvdy,udvdx,vdudy
    real*8,dimension(:,:),allocatable,save  :: viscu,viscv
    real*8,dimension(:,:),allocatable,save  :: us,vs
    real*8                                  :: nuh1,nuh2
    real*8                                  :: dudx1,dudx2,dudy1,dudy2
    real*8                                  :: dvdy1,dvdy2,dvdx1,dvdx2  !Jaap
	real*8                                  :: dalfa                    !difference in grid angles
	real*8                                  :: uin,vin                  !u resp v-velocity corrected for grid angle change
	real*8                                  :: qin                      !specific discharge entering cell
	real*8                                  :: dzsdnavg                 !alongshore water level slope
	real*8,save                             :: fc
	real*8,dimension(:,:),allocatable,save  :: sinthm,costhm

    integer                                 :: imax,jmax,jmin

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

    if (.not. allocated(vsu) ) then
       allocate (   vsu(nx+1,ny+1))
       allocate (   usu(nx+1,ny+1))
       allocate (   vsv(nx+1,ny+1))
       allocate (   usv(nx+1,ny+1))
       allocate (   veu(nx+1,ny+1))
       allocate (   uev(nx+1,ny+1))
       allocate ( ududx(nx+1,ny+1))
       allocate ( vdvdy(nx+1,ny+1))
       allocate ( udvdx(nx+1,ny+1))
       allocate ( vdudy(nx+1,ny+1))
       allocate ( viscu(nx+1,ny+1))
       allocate ( viscv(nx+1,ny+1))
       allocate (    us(nx+1,ny+1))
       allocate (    vs(nx+1,ny+1))
       allocate (sinthm(nx+1,ny+1))
       allocate (costhm(nx+1,ny+1))

       if (par%secorder == 1) then
          allocate(vv_old(nx+1,ny+1)); vv_old = vv
          allocate(uu_old(nx+1,ny+1)); uu_old = uu
          allocate(zs_old(nx+1,ny+1)); zs_old = zs
       endif

       vu      =0.d0
       vsu     =0.d0
       usu     =0.d0
       vsv     =0.d0
       usv     =0.d0
       uv      =0.d0
       veu     =0.d0
       uev     =0.d0
       ueu     =0.d0
       vev     =0.d0
       ududx   =0.d0
       vdvdy   =0.d0
       udvdx   =0.d0
       vdudy   =0.d0
       viscu   =0.d0
       viscv   =0.d0
       us      =0.d0
       vs      =0.d0
       hum     =0.d0
       hvm     =0.d0
       u       =0.d0
       v       =0.d0
       ue      =0.d0
       ve      =0.d0
       fc      =2.d0*par%wearth*sin(par%lat)
    endif
       
    ! update bedfriction coefficient cf
    if (trim(par%bedfriction)=='white-colebrook') then
       cf = par%g/(18.d0*log10(4.*hh/min(hh,D90top)))**2 ! cf = g/C^2 where C = 18*log(4*hh/D90)
    endif

    ! Super fast 1D
    if (ny==0) then
      j1 = 1
    else
      j1 = 2
    endif
    
    ! Add vertical discharges
    call discharge_boundary_v(s,par)
    
    !
    ! zs=zs*wetz
    ! Water level slopes
    do j=1,ny+1
       do i=2,nx
          dzsdx(i,j)=(zs(i+1,j)+ph(i+1,j)-zs(i,j)-ph(i,j))/dsu(i,j)
       end do
    end do
    !    do j=2,ny
    do j=1,ny ! Dano need to get correct slope on boundary y=0
       do i=1,nx+1
          dzsdy(i,j)=(zs(i,j+1)+ph(i,j+1)-zs(i,j)-ph(i,j))/dnv(i,j)
       end do
    end do
    
    ! wwvv in the next lines
    !  hu(i,j) is more or less a function of hh(i+1,j)
    !  In the parallel case, this needs some action because
    !  for processes not at the bottom, the last row of
    !  hu (hu(nx+1,:)) has to be taken from the neighbour below
    !  hu is used later on in this subroutine, so we have to insert
    !  an mpi call.
    !  The same for hum
    ! Of course, no action is necessary if hu(nx+1:) is never used...
    do j=1,ny+1
       do i=1,nx+1 !Ap
          ! Water depth in u-points do momentum equation: mean
		  !!
		  !! ARBJ: mean water depth or weighted water depth? How to deal with this in curvi-linear?
		  !!
          hum(i,j)=max(.5d0*(hh(i,j)+hh(min(nx,i)+1,j)),par%eps) 
       end do
    end do
    ! wwvv here the mpi code to communicate a row of hu
    ! we send to the neighbour above and receive from the neighbour
    ! below:
#ifdef USEMPI
    call xmpi_shift(hum,'m:')
#endif
    ! Wetting and drying criterion (only do momentum balance)
    do j=1,ny+1
       do i=1,nx+1
          if(hu(i,j)>par%eps .and. hum(i,j)>par%eps) then  ! Jaap and Pieter: If you want to compute correct advection term
             wetu(i,j)=1                                  ! then both hu and hum should be larger than par%eps. It is not
          else                                             ! necessarily true that if hu>par%eps also hum>par%eps.
             wetu(i,j)=0
          end if
       end do
    end do
    ! wwvv about the same for hv, only in the left-right direction
    ! hv(i,j) is more or less a function of hh(i,j+1)
    ! so in the parallel case, hv(:,ny+1) has to be collected
    ! from the right neighbour
    ! the same for hvm
    do j=1,ny+1
       do i=1,nx+1
          ! Water depth in v-points do momentum equation: mean
          hvm(i,j)=max(.5d0*(hh(i,j)+hh(i,min(ny,j)+1)),par%eps)
       end do
    end do
    ! send to the left, read from the right
#ifdef USEMPI
    call xmpi_shift(hvm,':n')
#endif
    ! Wetting and drying criterion (only do momentum balance)
    do j=1,ny+1
       do i=1,nx+1
          if(hv(i,j)>par%eps .and. hvm(i,j)>par%eps) then
             wetv(i,j)=1
          else
             wetv(i,j)=0
          end if
       end do
    end do

    ! Jaap Wetting and drying criterion eta points
    do j=1,ny+1
       do i=1,nx+1
          !A eta point is wet if any of the surrounding velocity points is wet...
          !            wetz(i,j) = min(1,wetu(max(i,2)-1,j)+wetu(i,j)+wetv(i,j)+wetv(i,max(j,2)-1))
          if(hh(i,j)>par%eps) then
             wetz(i,j)=1
          else
             wetz(i,j)=0
          end if
       end do
    end do

    !
    ! X-direction
    !
    do j=j1,max(ny,1)
       do i=2,nx
          ududx(i,j)        = 0.d0
          qin               = .5d0*(qx(i,j)+qx(i-1,j))
          if (qin>0) then
             dalfa          = alfau(i,j)-alfau(i-1,j)
             uin            = uu(i-1,j)*cos(dalfa) + vu(i-1,j)*sin(dalfa)
             ududx(i,j)     = ududx(i,j) + qin*(uu(i,j)-uin)*dnz(i,j)/hum(i,j)*dsdnui(i,j)
          endif
          qin               = -.5d0*(qx(i,j)+qx(i+1,j))
          if (qin>0) then
             dalfa          = alfau(i,j)-alfau(i+1,j)
             uin            = uu(i+1,j)*cos(dalfa) + vu(i+1,j)*sin(dalfa)
             ududx(i,j)     = ududx(i,j) + qin*(uu(i,j)-uin)*dnz(i+1,j)/hum(i,j)*dsdnui(i,j)
          endif
       end do
    end do
    ! wwvv fix border rows and columns of ududx
#ifdef USEMPI
    call xmpi_shift(ududx,'1:')
    call xmpi_shift(ududx,'m:')
    call xmpi_shift(ududx,':1')
    call xmpi_shift(ududx,':n')
#endif
    do j=2,ny
       do i=1,nx
          vdudy(i,j)        = 0.d0
          qin               = .5d0*(qy(i,j-1)+qy(i+1,j-1))
          if (qin>0) then
             dalfa          = alfau(i,j)-alfau(i,j-1)
             uin            = uu(i,j-1)*cos(dalfa) + vu(i,j-1)*sin(dalfa)
             vdudy(i,j)     = vdudy(i,j) + qin*(uu(i,j)-uin)*dsc(i,j-1)/hum(i,j)*dsdnui(i,j)
          endif
          qin               = -.5d0*(qy(i,j)+qy(i+1,j))
          if (qin>0) then
             dalfa          = alfau(i,j)-alfau(i,j+1)
             uin            = uu(i,j+1)*cos(dalfa) + vu(i,j+1)*sin(dalfa)
             vdudy(i,j)     = vdudy(i,j) + qin*(uu(i,j)-uin)*dsc(i,j)/hum(i,j)*dsdnui(i,j)
          endif
       end do
    end do
    ! wwvv fix border rows and columns of vdudy
#ifdef USEMPI
    ! call xmpi_shift(vdudy,'1:') ! not necessary wwvv
    call xmpi_shift(vdudy,'m:')
    call xmpi_shift(vdudy,':1')
    call xmpi_shift(vdudy,':n')
#endif
    !

    ! Jaap: Slightly changes approach; 1) background viscosity is user defined or obtained from Smagorinsky, 2) nuh = max(nuh,roller induced viscosity)
    if (par%smag == 1) then
       !Use smagorinsky subgrid model
       call visc_smagorinsky(s,par)
    else
       nuh = par%nuh
    endif

    do j=j1,max(ny,1)
       do i=2,nx
          nuh(i,j) = max(nuh(i,j),par%nuhfac*hh(i,j)*(DR(i,j)/par%rho)**(1.0d0/3.0d0)) ! Ad: change to max
       end do
    end do

    do j=j1,max(ny,1)
       do i=2,nx
       !write(*,*)i,j,2
          dudx1 = nuh(i+1,j)*hh(i+1,j)*(uu(i+1,j)-uu(i,j))/dsz(i+1,j)         
          dudx2 = nuh(i,j)  *hh(i  ,j)*(uu(i,j)-uu(i-1,j))/dsz(i,j)           
          viscu(i,j) = (1.0d0/hum(i,j))*( 2*(dudx1-dudx2)/(dsz(i,j)+dsz(i+1,j)) )
       end do
    end do

    if (par%smag == 1) then
       !
       !   For non constant eddy viscosity the stress terms read:
       !
       !   d           d           d      d         d
       !   -- [ 2* mu -- (U) ]  + --[ mu --(U) +mu --(V) ]
       !   dx         dx          dy     dy        dx
       !
       !                d         d
       !   Only when   -- [mu] = --[mu] = 0 we have (for incompressible flow)
       !               dx        dy
       !
       !         2            2
       !        d            d
       !    mu --2 [U] + mu --2 [U]
       !       dx           dy
       !
       viscu = 2.0d0*viscu
    endif

    do j=2,ny
       do i=2,nx
          !Nuh is defined at eta points, interpolate from four surrounding points
          nuh1  = .25d0*(nuh(i,j)+nuh(i+1,j)+nuh(i+1,j+1)+nuh(i,j+1))
          nuh2  = .25d0*(nuh(i,j)+nuh(i+1,j)+nuh(i+1,j-1)+nuh(i,j-1))

          dudy1 = nuh1 *.5d0*(hvm(i,j  )+hvm(i+1,j  ))*(uu(i,j+1)-uu(i,j))/dnc(i,j)
          dudy2 = nuh2 *.5d0*(hvm(i,j-1)+hvm(i+1,j-1))*(uu(i,j)-uu(i,j-1))/dnc(i,j-1)
          viscu(i,j) = viscu(i,j) + (1.0d0/hum(i,j))*( 2.0d0*(dudy1-dudy2)/(dnc(i,j)+dnc(i,j-1)) )*wetu(i,j+1)*wetu(i,j-1)
       end do
    end do

    if (par%smag == 1) then
       do j=2,ny
          do i=2,nx
             !Nuh is defined at eta points, interpolate from four surrounding points
             nuh1  = .25d0*(nuh(i,j)+nuh(i+1,j)+nuh(i+1,j+1)+nuh(i,j+1))
             nuh2  = .25d0*(nuh(i,j)+nuh(i+1,j)+nuh(i+1,j-1)+nuh(i,j-1))

             dvdx1 = nuh1*.5d0*(hvm(i,j  )+hvm(i+1,j  ))*(vv(i+1,j  )-vv(i,j  ))/dsc(i,j)
             dvdx2 = nuh2*.5d0*(hvm(i,j-1)+hvm(i+1,j-1))*(vv(i+1,j-1)-vv(i,j-1))/dsc(i,j-1)
             viscu(i,j) = viscu(i,j) + (1.d0/hum(i,j))*(dvdx1-dvdx2)/dnz(i,j)*real(wetv(i+1,j) &
                  * wetv(i,j)*wetv(i+1,j-1)*wetv(i,j-1),8)
          enddo
       enddo
    endif !smag ==1 and ny>0
    !
    ! Bed friction term
	where (wetu==1)
	   taubx=cf*par%rho*ueu*sqrt((1.16d0*urms)**2+vmageu**2) !Ruessink et al, 2001
	elsewhere
	   taubx = 0.d0
	endwhere
    !
    ! Explicit Euler step momentum u-direction
    !
    do j=j1,max(ny,1)
       ! do i=2,nx-1   ! wwvv uu(nx,:) is never computed in this subroutine, is that ok?
       if (xmpi_isbot) then
          imax = nx-1
       else
          imax=nx
       endif
       do i=2,imax ! wwvv with this modification, parallel and serial version
          ! give the same results. If this modification is not ok, then
          ! we have a problem
          if(wetu(i,j)==1) then
             uu(i,j)=uu(i,j)-par%dt*(ududx(i,j)+vdudy(i,j)-viscu(i,j) & !Ap,Robert,Jaap
                  + par%g*dzsdx(i,j) &
                  + taubx(i,j)/(par%rho*hu(i,j)) &  ! Dano: changed hum to hu NOT cf volume approach
                  - par%lwave*Fx(i,j)/(par%rho*max(hum(i,j),par%hmin)) &
                  - fc*vu(i,j) &
                  - par%rhoa*par%Cd*windsu(i,j)**2/(par%rho*hum(i,j)))
          else
             uu(i,j)=0.0d0
          end if
       end do
    end do
   	! Lateral boundary conditions for uu
   	if (ny>0) then
	   if (xmpi_isleft) then !Dano/Robert only on outer boundary
          uu(1:nx+1,1)=uu(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
          uu(1:nx+1,ny+1)=uu(1:nx+1,ny) ! RJ: can also be done after continuity but more appropriate here
       endif
    endif
#ifdef USEMPI
    ! wwvv qx is used later on, also the first row, fix the first row
    ! of uu first
    call xmpi_shift(uu,'1:')
    call xmpi_shift(uu,'m:')
    call xmpi_shift(uu,':1')
    call xmpi_shift(uu,':n')
#endif
    !
    ! Y-direction
    !
	! Robert: Complicated
	!         In overall model (not mpi-subdomains) vv(:,ny+1) is never calculated and has dummy values
	!         In overall model vv(:,1) and vv(:,ny) are used as boundary conditions
	!                          vv(:,1) and vv(:,ny) = 0 for wall boundary conditions
	!                          vv(:,1) and vv(:,ny) = copied from internal in Neumann conditions
	!                          vv(:,1) and vv(:,ny) = calculated without advection terms in no_advec conditions
	!                          vv(:,1) and vv(:,ny) = calculated with advection terms in free conditions 
	!         In overall model (or in non-MPI) vv is calculated over j=2,ny-1
	!
	!         In mpi subdomain vv(:,1) and vv(:,ny+1) are communicated and no boundary conditions are given.
	!         In mpi subdomain vv(:,ny+1) is NOT a dummy and can be used to calculate advection terms
	!         In mpi subdomain vv is calculated over j=2,ny

	! Advection term vdvdy
	! Robert: in MPI every model subdomain has vv(:,ny+1) from other domain
	!         at overall model boundaries vv(:,ny+1) assumed to be vv(:,ny) and vv(:,0) = vv(:,1)
	!         so calculate vdvdy where we can
    if (xmpi_isright .and. ny>0) then  ! no such condition needed for _isleft, because vdvdy(:,1) not needed in mpi subdomain
       jmax = ny-1
    elseif (xmpi_isright .and. ny==0) then
	   jmax = 1
	else
       jmax = ny
    endif
	vdvdy=0.d0
	! calculate true vdvdy up to ny in central domains and up to ny-1 on isright
	do j=2,jmax
       do i=2,nx
          qin               = .5d0*(qy(i,j)+qy(i,j-1))
          if (qin>0) then
             dalfa          = alfav(i,j)-alfav(i,j-1)
             vin            = vv(i,j-1)*cos(dalfa) - uv(i,j-1)*sin(dalfa)
             vdvdy(i,j)     = vdvdy(i,j) + qin*(vv(i,j)-vin)*dsz(i,j)/hvm(i,j)*dsdnvi(i,j)
          endif
          qin               = -.5d0*(qy(i,j)+qy(i,j+1))
          if (qin>0) then
             dalfa          = alfav(i,j)-alfav(i,j+1)
             vin            = vv(i,j+1)*cos(dalfa) - uv(i,j+1)*sin(dalfa)
             vdvdy(i,j)     = vdvdy(i,j) + qin*(vv(i,j)-vin)*dsz(i,j+1)/hvm(i,j)*dsdnvi(i,j)
          endif
       enddo
    enddo
	if (ny>0) then 
       ! Global boundary conditions for vdvdy(:,1) and vdvdy(:,ny), global vdvdy(:,ny+1) not needed anywhere
       if (xmpi_isleft) then
          ! (vv(:,1)-vv(:,0))/dy == 0 so only second part of the vdvdy equation:
          do i=2,nx
		     qin            = -.5d0*(qy(i,1)+qy(i,2))
             if (qin>0) then
                dalfa       = alfav(i,1)-alfav(i,2)
                vin         = vv(i,2)*cos(dalfa) - uv(i,2)*sin(dalfa)
                vdvdy(i,1)  = vdvdy(i,1) + qin*(vv(i,1)-vin)*dsz(i,2)/hvm(i,1)*dsdnvi(i,1)
             endif
		  enddo
       endif
	   if (xmpi_isright) then
	      ! (vv(:,ny+1)-vv(:,ny))/dy == 0 so only first part of the vdvdy equation:
		  do i=2,nx
		     qin            = .5d0*(qy(i,ny)+qy(i,ny-1))
		     if (qin>0) then
                dalfa       = alfav(i,ny)-alfav(i,ny-1)
                vin         = vv(i,ny-1)*cos(dalfa) - uv(i,ny-1)*sin(dalfa)
                vdvdy(i,ny) = vdvdy(i,ny) + qin*(vv(i,ny)-vin)*dsz(i,ny)/hvm(i,ny)*dsdnvi(i,ny)
             endif
          enddo
       endif
    endif

    udvdx=0.d0
    if (ny>0) then
	   ! Robert: udvdx not usually needed at j = 1
	   do j=1,ny !1,ny instead of 2,ny
          do i=2,nx
             qin            = .5d0*(qx(i-1,j)+qx(i-1,j+1))
             if (qin>0) then
                dalfa       = alfav(i,j)-alfav(i-1,j)
                vin         = vv(i-1,j)*cos(dalfa) - uv(i-1,j)*sin(dalfa)
                udvdx(i,j)  = udvdx(i,j) + qin*(vv(i,j)-vin)*dnc(i-1,j)/hvm(i,j)*dsdnvi(i,j)
             endif
             qin            = -.5d0*(qx(i,j)+qx(i,j+1))
             if (qin>0) then
                dalfa       = alfav(i,j)-alfav(i+1,j)
                vin         = vv(i+1,j)*cos(dalfa) - uv(i+1,j)*sin(dalfa)
                udvdx(i,j)  = udvdx(i,j) + qin*(vv(i,j)-vin)*dnc(i,j)/hvm(i,j)*dsdnvi(i,j)
             endif
          end do
       end do
    else
	   do i=2,nx
	      qin               = qx(i-1,1)
          if (qin>0) then
             dalfa          = alfav(i,1)-alfav(i-1,1)
             vin            = vv(i-1,1)*cos(dalfa) - uv(i-1,1)*sin(dalfa)
             udvdx(i,1)     = udvdx(i,1) + qin*(vv(i,1)-vin)*dnc(i-1,1)/hvm(i,1)*dsdnvi(i,1)
          endif
          qin               = -qx(i,1)
          if (qin>0) then
             dalfa          = alfav(i,1)-alfav(i+1,1)
             vin            = vv(i+1,1)*cos(dalfa) - uv(i+1,1)*sin(dalfa)
             udvdx(i,1)     = udvdx(i,1) + qin*(vv(i,1)-vin)*dnc(i,1)/hvm(i,1)*dsdnvi(i,1)
          endif
       enddo
    endif
	!
	viscv =0.d0	
    do j=2,ny
       do i=2,nx
          dvdy1 = nuh(i,j+1)*hh(i,j+1)*(vv(i,j+1)-vv(i,j))/dnz(i,j+1)
          dvdy2 = nuh(i,j)  *hh(i,j  )*(vv(i,j)-vv(i,j-1))/dnz(i,j)
          viscv(i,j) = (1.0d0/hvm(i,j))* 2*(dvdy1-dvdy2)/(dnz(i,j)+dnz(i,j+1))*wetv(i,j+1)*wetv(i,j-1)
       end do
    end do
	! Robert: global boundary at (:,1) edge
	if (ny>0) then
       if (xmpi_isleft) then
          viscv(:,1) = viscv(:,2)
       endif
       if (xmpi_isright) then
          viscv(:,ny) = viscv(:,ny-1)
       endif
    endif
	!
	! Viscosity
	!
    if (par%smag == 1) then
       viscv = 2.0d0*viscv
    endif

    nuh = par%nuhv*nuh !Robert en Ap: increase nuh interaction in d2v/dx2
    do j=1,max(ny,1)
       jp1 = min(j+1,ny+1)
       do i=2,nx
          !Nuh is defined at eta points, interpolate from four surrounding points
          nuh1  = .25d0*(nuh(i,j)+nuh(i+1,j)+nuh(i+1,jp1)+nuh(i,jp1))
          nuh2  = .25d0*(nuh(i,j)+nuh(i-1,j)+nuh(i-1,jp1)+nuh(i,jp1))

          dvdx1 = nuh1*.5d0*(hum(i  ,j)+hum(i  ,jp1))*(vv(i+1,j)-vv(i,j))/dsc(i,j)
          dvdx2 = nuh2*.5d0*(hum(i-1,j)+hum(i-1,jp1))*(vv(i,j)-vv(i-1,j))/dsc(i-1,j)
          viscv(i,j) = viscv(i,j) + (1.0d0/hvm(i,j))*( 2*(dvdx1-dvdx2)/(dsc(i-1,j)+dsc(i,j)) )*wetv(i+1,j)*wetv(i-1,j)
       end do
    end do
    !
    if (par%smag == 1) then
       do j=1,max(ny,1)
          jp1 = min(j+1,ny+1)
          do i=2,nx
             !Nuh is defined at eta points, interpolate from four surrounding points
             nuh1  = .25d0*(nuh(i,j)+nuh(i+1,j)+nuh(i+1,jp1)+nuh(i,jp1))
             nuh2  = .25d0*(nuh(i,j)+nuh(i-1,j)+nuh(i-1,jp1)+nuh(i,jp1))

             dudy1 = nuh1 *.5d0*(hum(i  ,j)+hum(i  ,jp1))*(uu(i,jp1  )-uu(i,j  ))/dnc(i,j)
             dudy2 = nuh2 *.5d0*(hum(i-1,j)+hum(i-1,jp1))*(uu(i-1,jp1)-uu(i-1,j))/dnc(i-1,j)
             
             viscv(i,j) = viscv(i,j) + (1.d0/hvm(i,j))*(dudy1-dudy2)/dsz(i,j)  &
                  * real(wetu(i,jp1)*wetu(i,j)*wetu(i-1,jp1)*wetv(i-1,j),8)
          enddo
       enddo
    endif
	!
    ! Bed friction term
    !
	where (wetv==1)
	   tauby=cf*par%rho*vev*sqrt((1.16d0*urms)**2+vmagev**2) !Ruessink et al, 2001
	elsewhere
	   tauby = 0.d0
	endwhere
    !
    ! Explicit Euler step momentum v-direction
    !
	if (ny==0) then
	   jmin = 1
	else
	   jmin = 2
	   if (ny==2) then
	      jmax = 2 ! Robert: very special case of ny=2 and xmpi_isright would otherwise lead to no calculation of vv with Neumann boundaries
	   endif
	endif
    !
    do j=jmin,jmax
       do i=2,nx !jaap instead of nx+1
          if(wetv(i,j)==1) then
             ! Robert: ensure taubx always has the same sign as uu (always decelerates)
             ! Dano: I don't agree 
             vv(i,j)=vv(i,j)-par%dt*(udvdx(i,j)+vdvdy(i,j)-viscv(i,j)& !Ap,Robert,Jaap
                  + par%g*dzsdy(i,j)&
                  + tauby(i,j)/(par%rho*hv(i,j)) &  ! Dano: hv instead of hvm, NOT cf volume approach
                  - par%lwave*Fy(i,j)/(par%rho*max(hvm(i,j),par%hmin)) &
                  + fc*uv(i,j) &
                  - par%rhoa*par%Cd*windnv(i,j)**2/(par%rho*hvm(i,j)))
          else
             vv(i,j)=0.0d0
          end if
       end do
    end do
    ! Communicate vv at internal boundaries
#ifdef USEMPI
    call xmpi_shift(vv,'1:')
    call xmpi_shift(vv,'m:')
    call xmpi_shift(vv,':1')
    call xmpi_shift(vv,':n')
#endif
    ! Robert: Boundary conditions along the global boundaries
	! Function flow_lat_bc located in boundaryconditions.F90
	! function call takes care of 1D vs 2D models and boundary condition types
	if (ny>0) then
       if (xmpi_isleft) then
	      vv(:,1)=flow_lat_bc(s,par,par%right,1,2,udvdx(:,1),vdvdy(:,1),viscv(:,1))
	   endif
       if (xmpi_isright) then
	      vv(:,ny)=flow_lat_bc(s,par,par%left,ny,ny-1,udvdx(:,ny),vdvdy(:,ny),viscv(:,ny))
	   endif
	endif

#ifndef USEMPI
    if (par%nonh==1) then
       !Do explicit predictor step with pressure
       call nonh_explicit(s,par,nuh)
    end if

    if (par%secorder==1) then
       !Call second order correction to the advection
       call flow_secondorder_advUV(s,par,uu_old,vv_old)
    end if
#endif

! Robert: include again when 2nd order returned
#ifdef USEMPI
    call xmpi_shift(uu,':1')
    call xmpi_shift(uu,':n')
    call xmpi_shift(uu,'1:')
    call xmpi_shift(uu,'m:')
    call xmpi_shift(vv,':1')
    call xmpi_shift(vv,':n')
    call xmpi_shift(vv,'1:')
    call xmpi_shift(vv,'m:')
#endif

    ! Pieter and Jaap: update hu en hv for continuity
    do j=1,ny+1
       do i=1,nx+1 !Ap
          ! Water depth in u-points do continuity equation: upwind
          if (uu(i,j)>par%umin) then
             if (par%oldhu == 1) then 
			    hu(i,j)=hh(i,j)
			 else
                hu(i,j)=zs(i,j)-max(zb(i,j),zb(min(nx,i)+1,j))
			 endif
          elseif (uu(i,j)<-par%umin) then
		     if (par%oldhu == 1) then 
                hu(i,j)=hh(min(nx,i)+1,j)
			 else
                hu(i,j)=zs(min(nx,i)+1,j)-max(zb(i,j),zb(min(nx,i)+1,j))
		     endif
          else
             hu(i,j)=max(max(zs(i,j),zs(min(nx,i)+1,j))-max(zb(i,j),zb(min(nx,i)+1,j)),par%eps)
          end if
       end do
    end do

    hu = max(hu,0.d0)

    do j=1,ny+1
       do i=1,nx+1
          ! Water depth in v-points do continuity equation: upwind
          if (vv(i,j)>par%umin) then
		     if (par%oldhu == 1) then 
                hv(i,j)=hh(i,j)
			 else
                hv(i,j)=zs(i,j)-max(zb(i,j),zb(i,min(ny,j)+1))
			 endif
          elseif (vv(i,j)<-par%umin) then
		     if (par%oldhu == 1) then 
                hv(i,j)=hh(i,min(ny,j)+1)
			 else
                hv(i,j)=zs(i,min(ny,j)+1)-max(zb(i,j),zb(i,min(ny,j)+1))
			 endif
          else
             hv(i,j)=max(max(zs(i,j),zs(i,min(ny,j)+1))-max(zb(i,j),zb(i,min(ny,j)+1)),par%eps)
          end if
       end do
    end do
    hv = max(hv,0.d0)

#ifdef USEMPI
    call xmpi_shift(hu ,'m:')
    call xmpi_shift(hv ,':n')
#endif


#ifndef USEMPI
    if (par%nonh==1) then
       ! do non-hydrostatic pressure compensation to solve short waves
       call nonh_cor(s,par)
    end if
#endif

    ! Flux in u-point
    qx=uu*hu
    ! Flux in v-points
    ! first column of qy is used later, and it is defined in the loop above
    ! no communication  necessary at this point
    qy=vv*hv
    !
    ! Add horizontal discharges
	!
    call discharge_boundary_h(s,par)
    !
    ! Update water level using continuity eq.
    !
    if (xmpi_isright) then
       jmax=ny
    else
       jmax=ny+1
    endif
    if (xmpi_isbot) then
       imax=nx
    else
       imax=nx+1
    endif
    if (ny>0) then
      do j=2,jmax
        do i=2,imax
          dzsdt(i,j) = (-1.d0)*( qx(i,j)*dnu(i,j)-qx(i-1,j)*dnu(i-1,j)  &
                               + qy(i,j)*dsv(i,j)-qy(i,j-1)*dsv(i,j-1) )*dsdnzi(i,j) &
                        - gww(i,j)
        end do
      end do
      zs(2:nx,2:ny) = zs(2:nx,2:ny)+dzsdt(2:nx,2:ny)*par%dt !Jaap nx instead of nx+1
    else
       j=1
       do i=2,imax
          dzsdt(i,j) = (-1.d0)*( qx(i,j)*dnu(i,j)-qx(i-1,j)*dnu(i-1,j) )*dsdnzi(i,j) &
                        - gww(i,j)
       end do
       zs(2:nx,1) = zs(2:nx,1)+dzsdt(2:nx,1)*par%dt !Jaap nx instead of nx+1
    endif !ny>0
    ! call discharge_boundary(s,par)
    !

#ifndef USEMPI
    if (par%secorder == 1) then
       !Second order correction
       call flow_secondorder_con(s,par,zs_old)
    endif
#endif
    !
    ! Lateral boundary conditions
    !
    if (ny>0) then
      ! RJ: Neumann water levels in case of right = 1 or right = 0
      do i=1,nx+1
         ! Jaap multiply with wetz(i,ny+1)*wetz(i,1) here to prevent prsssure grdaient over land
         dzsdnavg=wetz(i,ny+1)*wetz(i,1)*(zs0(i,ny+1)-zs0(i,1))/(ndist(i,ny+1)-ndist(i,1))
         ! Lateral boundary at y=0
         if (xmpi_isleft) then
           zs(i,1)=max(zs(i,2) - dzsdnavg*dnv(i,1),zb(i,1))
         endif
		 ! Lateral boundary at y=ny+1
         if (xmpi_isright) then
            zs(i,ny+1)=max(zs(i,ny) + dzsdnavg*dnv(i,ny),zb(i,ny+1))
         endif
       enddo
    endif !ny>0


    ! wwvv zs, uu, vv have to be communicated now, because they are used later on
#ifdef USEMPI
    call xmpi_shift(zs,':1')
    call xmpi_shift(zs,':n')
    call xmpi_shift(zs,'1:')
    call xmpi_shift(zs,'m:')
    call xmpi_shift(dzsdt,':1')  ! wwvv dzsdt maybe not necessary because first and last columns are not used
    call xmpi_shift(dzsdt,':n')
    call xmpi_shift(dzsdt,'1:')
    call xmpi_shift(dzsdt,'m:')
    call xmpi_shift(qx,':1')  ! wwvv qx maybe not necessary because first and last columns are not used
    call xmpi_shift(qx,':n')
    call xmpi_shift(qx,'1:')
    call xmpi_shift(qx,'m:')
    call xmpi_shift(qy,':1')
    call xmpi_shift(qy,':n')
    call xmpi_shift(qy,'1:')! wwvv qy maybe not necessary because first and last rows are not used
    call xmpi_shift(qy,'m:')
#endif

    if (par%secorder == 1) then
       vv_old = vv
       uu_old = uu
       zs_old = zs
    endif

    ! offshore boundary
    !
    ! U and V in cell centre; do output and sediment stirring
    !
    u(2:nx,:)=0.5d0*(uu(1:nx-1,:)+uu(2:nx,:))
    if(xmpi_istop) then
       u(1,:)=uu(1,:)
    endif
    if(xmpi_isbot) then
       u(nx+1,:)=u(nx,:)
    endif
    
#ifdef USEMPI
    call xmpi_shift(u,'1:')
    call xmpi_shift(u,'m:')
    call xmpi_shift(u,':1')
    call xmpi_shift(u,':n')
#endif

    if (ny>0) then
        v(:,2:ny)=0.5d0*(vv(:,1:ny-1)+vv(:,2:ny))
        if(xmpi_isleft) then
           v(:,1)=vv(:,1)
        endif
        if(xmpi_isright) then
           v(:,ny+1)=v(:,ny)        ! bas: need this for calculation of ee in wci routine
        endif
        !Ap
        v(nx+1,:)=v(nx,:)
    else ! Dano
       v=vv
    endif !ny>0

#ifdef USEMPI
    call xmpi_shift(v,':1')
    call xmpi_shift(v,':n')
    call xmpi_shift(v,'1:')
    call xmpi_shift(v,'m:')
#endif
    ! Robert + Jaap: compute derivatives of u and v
    !
    
    sinthm = sin(thetamean-alfaz)
    costhm = cos(thetamean-alfaz)
    
    ! V-velocities at u-points
    if (ny>0) then
      vu(1:nx,2:ny)= 0.25d0*(vv(1:nx,1:ny-1)+vv(1:nx,2:ny)+ &
                             vv(2:nx+1,1:ny-1)+vv(2:nx+1,2:ny))
      ! how about boundaries?
      if(xmpi_isleft) then
        vu(:,1) = vu(:,2)
      endif
      if(xmpi_isright) then
        vu(:,ny+1) = vu(:,ny)
      endif
    else 
      vu(1:nx,1)= 0.5d0*(vv(1:nx,1)+vv(2:nx+1,1))
    endif !ny>0
    ! wwvv fill in vu(:1) and vu(:ny+1) for non-left and non-right processes
    !  and vu(nx+1,:)
#ifdef USEMPI
    call xmpi_shift(vu,':1')
    call xmpi_shift(vu,':n')
    call xmpi_shift(vu,'m:')
    ! Jaap: whu not vu,'1:' --> seesm to be necessary to compute vmagu that is used to compute Su
#endif
    vu=vu*wetu
    ! V-stokes velocities at U point
    if (ny>0) then
      vsu(1:nx,2:ny)=0.5d0*(ust(1:nx,2:ny)*sinthm(1:nx,2:ny)+ &
                            ust(2:nx+1,2:ny)*sinthm(2:nx+1,2:ny))
      if(xmpi_isleft) then
        vsu(:,1)=vsu(:,2)
      endif
      if(xmpi_isright) then
        vsu(:,ny+1) = vsu(:,ny)
      endif
    else
      vsu(1:nx,1)=0.5d0*(ust(1:nx,1)*sinthm(1:nx,1)+ &
                         ust(2:nx+1,1)*sinthm(2:nx+1,1))
    endif !ny>0
    ! wwvv same for vsu
#ifdef USEMPI
    call xmpi_shift(vsu,':1')
    call xmpi_shift(vsu,':n')
    call xmpi_shift(vsu,'m:')
#endif
    vsu = vsu*wetu
    ! U-stokes velocities at U point
    if (ny>0) then
      usu(1:nx,2:ny)=0.5d0*(ust(1:nx,2:ny)*costhm(1:nx,2:ny)+ &
                            ust(2:nx+1,2:ny)*costhm(2:nx+1,2:ny))
      if(xmpi_isleft) then
        usu(:,1)=usu(:,2)
      endif
      if(xmpi_isright) then
        usu(:,ny+1)=usu(:,ny)
      endif
    else
      usu(1:nx,1)=0.5d0*(ust(1:nx,1)*costhm(1:nx,1)+ &
                         ust(2:nx+1,1)*costhm(2:nx+1,1))
    endif !ny>0
    ! wwvv same for usu   
#ifdef USEMPI
    call xmpi_shift(usu,':1')
    call xmpi_shift(usu,':n')
    call xmpi_shift(usu,'m:')
#endif
    usu=usu*wetu

    ! V-euler velocities at u-point
    veu = vu - vsu
    ! U-euler velocties at u-point
    ueu = uu - usu
    ! Velocity magnitude at u-points
    ! vmagu=sqrt(uu**2+vu**2)
    ! Eulerian velocity magnitude at u-points
    vmageu=sqrt(ueu**2+veu**2)

    ! U-velocities at v-points
    if (ny>0) then
      uv(2:nx,1:ny)= .25d0*(uu(1:nx-1,1:ny)+uu(2:nx,1:ny)+ &
                            uu(1:nx-1,2:ny+1)+uu(2:nx,2:ny+1))
      ! boundaries?
      ! wwvv and what about uv(:,1) ?
      if(xmpi_isright) then
        uv(:,ny+1) = uv(:,ny)
      endif
    else
      uv(2:nx,1)= .5d0*(uu(1:nx-1,1)+uu(2:nx,1))
    endif !ny>0
    ! wwvv fix uv(:,ny+1) for non-right processes
    ! uv(1,:) and uv(nx+1,:) need to be filled in for
    ! non-bot or top processes
#ifdef USEMPI
    call xmpi_shift(uv,':n')
    call xmpi_shift(uv,':1')
    call xmpi_shift(uv,'1:')
    call xmpi_shift(uv,'m:')
#endif
    uv=uv*wetv
    ! V-stokes velocities at V point
    if (ny>0) then
      vsv(2:nx,1:ny)=0.5d0*(ust(2:nx,1:ny)*sinthm(2:nx,1:ny)+&
                            ust(2:nx,2:ny+1)*sinthm(2:nx,2:ny+1))
      if(xmpi_isleft) then
        vsv(:,1) = vsv(:,2)
      endif
      if(xmpi_isright) then
        vsv(:,ny+1) = vsv(:,ny)
      endif
    else
      vsv(2:nx,1)= ust(2:nx,1)*sinthm(2:nx,1)
    endif !ny>0
    ! wwvv fix vsv(:,1) and vsv(:,ny+1) and vsv(1,:) and vsv(nx+1,:)
#ifdef USEMPI
    call xmpi_shift(vsv,':n')
    call xmpi_shift(vsv,':1')
    call xmpi_shift(vsv,'1:')
    call xmpi_shift(vsv,'m:')
#endif

    vsv=vsv*wetv
    ! U-stokes velocities at V point
    if (ny>0) then
      usv(2:nx,1:ny)=0.5d0*(ust(2:nx,1:ny)*costhm(2:nx,1:ny)+&
                            ust(2:nx,2:ny+1)*costhm(2:nx,2:ny+1))
      if(xmpi_isleft) then
        usv(:,1) = usv(:,2)
      endif
      if(xmpi_isleft) then
        usv(:,ny+1) = usv(:,ny)
      endif
    else
      usv(2:nx,1)=ust(2:nx,1)*costhm(2:nx,1)
    endif !ny>0
    ! wwvv fix usv(:,1) and usv(:,ny+1) and usv(1,:) and usv(nx+1,:)
#ifdef USEMPI
    call xmpi_shift(usv,':n')
    call xmpi_shift(usv,':1')
    call xmpi_shift(usv,'1:')
    call xmpi_shift(usv,'m:')
#endif
    usv=usv*wetv

    ! V-euler velocities at V-point
    vev = vv - vsv
    ! U-euler velocties at V-point
    uev = uv - usv
    ! Velocity magnitude at v-points
    ! vmagv=sqrt(uv**2+vv**2)
    ! Eulerian velocity magnitude at v-points
    vmagev=sqrt(uev**2+vev**2)


    ! Ue and Ve in cell centre; do output and sediment stirring
    ue(2:nx,:)=0.5d0*(ueu(1:nx-1,:)+ueu(2:nx,:))
    ue(1,:)=ueu(1,:)
    ! wwvv ue(nx+1,:) ?
#ifdef USEMPI
    call xmpi_shift(ue,':1')
    call xmpi_shift(ue,':n')
    call xmpi_shift(ue,'1:')
    call xmpi_shift(ue,'m:')
#endif

    if (ny>0) then
      ve(:,2:ny)=0.5d0*(vev(:,1:ny-1)+vev(:,2:ny)) !Jaap ny+1
      ve(:,1)=vev(:,1)
    else
      ve(:,1) = vev(:,1)
    endif !ny>0
    ! wwvv vev(nx+1,:) ?
#ifdef USEMPI
    call xmpi_shift(ve,':1')
    call xmpi_shift(ve,':n')
    call xmpi_shift(ve,'1:')
    call xmpi_shift(ve,'m:')
#endif
    !
    hold =hh    ! wwvv ?  hold is never else used
    !
    hh=max(zs-zb,par%eps)

    maxzs=max(zs,maxzs)
    minzs=min(zs,minzs)

  end subroutine flow


subroutine visc_smagorinsky(s,par)
  use params
  use spaceparams
  use xmpi_module

  IMPLICIT NONE

! DATE               AUTHOR               CHANGES
!
! December 2010       Pieter Bart Smit     New Subroutine
! March    2010       Pieter Bart Smit     Changed formulation to standard smag. model

!-------------------------------------------------------------------------------
!                             DECLARATIONS
!-------------------------------------------------------------------------------

!--------------------------        PURPOSE         ----------------------------
!
!  Calculates the turbulent viscocity coefficient nuh according to the smagorinsky
!  subgrid model.
!
!--------------------------        METHOD          ----------------------------
!
! The turbulent viscocity is given as:
!
! nuh = C^2*dx*dy*Tau
!
! Tau =2^(1/2) *  [ (du/dx)^2 + (dv/dy)^2 + 1/2 * (du/dy + dv/dx)^2 ] ^ (1/2)
!
! Where
!
! dx,dy : grid size
! C     : Constant ~0.15 (set by par%nuh)
! Tau   : Measure for the magnitude of the turbulent stresses
!
!--------------------------     ARGUMENTS          ----------------------------

  type(spacepars),target                   ,intent(inout) :: s
  type(parameters)                         ,intent(in)    :: par

!--------------------------     LOCAL VARIABLES    ----------------------------

  real*8                                                  :: dudx !U Velocity gradient in x-dir
  real*8                                                  :: dudy !U Velocity gradient in y-dir
  real*8                                                  :: dvdx !V Velocity gradient in x-dir
  real*8                                                  :: dvdy !V Velocity gradient in y-dir
  real*8                                                  :: Tau  !Measure for magnitude viscous stresses
  real*8                                                  :: l    !Local gridcell area
  integer                                                 :: i    !Index variable
  integer                                                 :: j    !Index variable

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

!-------------------------------------------------------------------------------
!                             IMPLEMENTATION
!-------------------------------------------------------------------------------

  !MPI WARNING -> Check loop indices
  if (ny>2) then
    do j=2,ny
      do i=2,nx
        dudx = (uu(i,j)-uu(i-1,j))/dsz(i,j)
        dudy = .5d0*(uu(i,j+1) - uu(i,j-1) + uu(i-1,j+1) - uu(i-1,j-1))/(dnv(i,j)+dnv(i,j-1))
        dvdx = .5d0*(vv(i+1,j) - vv(i-1,j) + vv(i+1,j-1) - vv(i-1,j-1))/(dsu(i,j)+dsu(i-1,j))
        dvdy = (vv(i,j)-vv(i,j-1))/dnz(i,j)
        Tau  = sqrt(2.0d0 * dudx**2+2.0d0 * dvdy**2 + (dvdx+dudy)**2)
        l    = 1.d0/dsdnzi(i,j)
        nuh(i,j) = par%nuh**2 * l * Tau * real(wetu(i,j)*wetu(i-1,j)*wetv(i,j)*wetv(i,j-1),kind=8)
      enddo
    enddo
  
  else
  
    j = max(ny,1)
    
    do i=2,nx
        dudx = (uu(i,j)-uu(i-1,j))/dsz(i,j)
        dvdx = (vv(i+1,j) - vv(i-1,j) )/(dsu(i,j)+dsu(i-1,j))
        Tau  = sqrt(2.0d0 * dudx**2 + dvdx**2)
        
        if (par%dy > -1.d0) then
          l = dsz(i,j)*par%dy
        else
          l = dsz(i,j)**2
        endif
        
        nuh(i,j) = par%nuh**2 * l * Tau * real(wetu(i,j)*wetu(i-1,j),kind=8)
    enddo
  
  endif !ny>2

  if (ny>0) then
    if (xmpi_isleft)    nuh(:,1)      = nuh(:,2)
    if (xmpi_isright)   nuh(:,ny+1)   = nuh(:,ny)
  endif
  
  if (xmpi_istop)       nuh(1,:)      = nuh(2,:)
  if (xmpi_isbot)       nuh(nx+1,:)   = nuh(nx,:)

#ifdef USEMPI
  call xmpi_shift(nuh,'1:')
  call xmpi_shift(nuh,'m:')
  call xmpi_shift(nuh,':1')
  call xmpi_shift(nuh,':n')
#endif  
end subroutine visc_smagorinsky

end module flow_timestep_module