module wetcells_module
   implicit none
   save
   private
   public compute_wetcells
contains
   subroutine compute_wetcells(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,         only: parameters
      use spaceparamsdef, only: spacepars
      use paramsconst

      IMPLICIT NONE

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

      integer                                 :: i
      integer                                 :: j
      integer                                 :: upwinddist
      integer,dimension(:,:),allocatable,save :: weteb

      if(.not.allocated(weteb)) allocate(weteb(s%nx+1,s%ny+1))

      ! 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,s%ny+1
         do i=1,s%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?
            !!
            s%hum(i,j)=max(.5d0*(s%hh(i,j)+s%hh(min(s%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:
      ! Wetting and drying criterion (only do momentum balance)
      do j=1,s%ny+1
         do i=1,s%nx+1
            if(s%hu(i,j)>par%eps .and. s%hum(i,j)>par%eps) then  ! Jaap and Pieter: If you want to compute correct advection term
               s%wetu(i,j)=1                                  ! then both s%hu and s%hum should be larger than par%eps. It is not
            else                                             ! necessarily true that if s%hu>par%eps also s%hum>par%eps.
               s%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,s%ny+1
         do i=1,s%nx+1
            ! Water depth in v-points do momentum equation: mean
            s%hvm(i,j)=max(.5d0*(s%hh(i,j)+s%hh(i,min(s%ny,j)+1)),par%eps)
         end do
      end do
      ! Wetting and drying criterion (only do momentum balance)
      do j=1,s%ny+1
         do i=1,s%nx+1
            if(s%hv(i,j)>par%eps .and. s%hvm(i,j)>par%eps) then
               s%wetv(i,j)=1
            else
               s%wetv(i,j)=0
            end if
         end do
      end do
      ! Jaap Wetting and drying criterion eta points
      do j=1,s%ny+1
         do i=1,s%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(s%hh(i,j)>par%eps) then
               s%wetz(i,j)=1
            else
               s%wetz(i,j)=0
            end if
         end do
      end do

      ! Wetting and drying for wave module
      where(s%hh+par%delta*s%H>par%eps .or. s%wetz==1)
         s%wete = 1
      elsewhere
         s%wete = 0
      endwhere
      ! also need to surround wet cells
      weteb = 0
      ! upwind distance in wave propagation routine
      if (par%scheme==SCHEME_UPWIND_1) then
         upwinddist = 1
      else
         upwinddist = 2
      endif
      forall(i=1:s%nx+1,j=1:s%ny+1,s%wete(i,j)==0)
         weteb(i,j) = max(   maxval(  s%wetz( max(i-upwinddist,1):min(i+upwinddist,s%nx+1) ,j )   ), &
         maxval(  s%wetz( i, max(j-upwinddist,1):min(j+upwinddist,s%ny+1) )   ), &
         maxval(  s%wetu( max(i-1,1):min(i,s%nx+1) ,j                     )   ), &
         maxval(  s%wetv( i, max(j-1,1):min(j,s%ny+1)                     )   ) &
         )
         !weteb(i,j) = maxval(s%wete(max(i-2,1):min(i+2,s%nx+1),max(j-2,1):min(j+2,s%ny+1))
         !weteb(i,j) = max(   maxval(  s%wete( max(i-2,1):min(i+2,s%nx+1) ,j                          )   ), &
         !                    maxval(  s%wete( i                          ,max(j-2,1):min(j+2,s%ny+1) )   ) &
         !            )
      endforall
      where(s%wete==0)
         s%wete=weteb
      endwhere

   end subroutine compute_wetcells
end module wetcells_module