subroutine wave_bc(s,par,it)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
! 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 waveparams
use spaceparams

IMPLICIT NONE

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

integer                                     :: it,nt
integer                                     :: i,ii,new
integer, save                               :: old
integer                                     :: j
integer                                     :: itheta
integer                                     :: instat,E_idx
real*8                                      :: E1,ei,dum
real*8, save                                :: dtbcfile,rt,bcendtime
real*8                                      :: em,tshifted,h0t0,tnew
real*8,dimension(:)     ,allocatable,save   :: e01       ! [J/m2/rad] directional distribution of wave energy at boundary
real*8,dimension(:)     ,allocatable,save   :: fac1,fac2
real*8,dimension(:)     ,allocatable,save   :: tE,dataE,databi
real*8,dimension(:,:)   ,allocatable,save   :: ht
real*8,dimension(:)     ,allocatable,save   :: q1,q2,q
real*8,dimension(:,:)   ,allocatable,save   :: ee1, ee2
character*1                                 :: bl
character*80                                :: fname, ebcfname,qbcfname
real*8                                      :: degrad,E0
real*8,dimension(:),allocatable             :: dist,factor
logical                                     :: startbcf


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

if (.not. allocated(fac1)) then
   allocate(ht      (2,ny+1))
   allocate(fac1    (nx))
   allocate(fac2    (nx))
   allocate(e01(1:ntheta))
   allocate(dist(1:ntheta))
   allocate(factor(1:ntheta))
endif
!
!  GENERATE AND READ-IN WAVE BOUNDARY CONDITIONS
!
! added for bound long wave comp Ad 28 march 2006
dtheta = par%dtheta*par%px/180;
startbcf=.false.
if(par%t==par%dt) then
    write(*,*)'Setting up boundary conditions'
    startbcf=.true.                     ! trigger read from bcf for instat 3,4,5,7
    bcendtime=huge(0.0)                 ! initial assumption for instat 3,4,5,7
    if (par%instat==2) then
       open( unit=7, file='bc\gen.ezs')
5      read(7,'(a)')bl
       if(bl.eq.'*') goto 5
       read(7,*)nt
       allocate(dataE  (nt))
       allocate(tE     (nt))
       do i=1,nt
          read(7,*) tE(i),dum,dataE(i)
       end do
       close(7)
       par%Emean=sum(dataE)/nt
    elseif (par%instat==3) then
       open( unit=7, file='bc\gen.ezs')
6      read(7,'(a)')bl
       if(bl.eq.'*') goto 6
       read(7,*)nt
       allocate(dataE  (nt))
       allocate(databi (nt))
       allocate(tE     (nt))
       do i=1,nt
          read(7,*) tE(i),databi(i),dataE(i)
       end do
       close(7)
       par%Emean=sum(dataE)/nt
    elseif (par%instat==4) then
       call makebcf(par,s,wp)
    elseif (par%instat==5) then
       call makebcf(par,s,wp)
    elseif (par%instat==6) then
       call makebcf(par,s,wp)
    elseif (par%instat==7) then
       par%listline=1
    endif
    !
    ! Directional distribution
    !
    if(par%instat==0 .or. par%instat==1 .or. par%instat==2 .or. par%instat==3)then
        degrad=par%px/180.
        dist=(cos(theta-theta0))**par%m;
        do i=1,ntheta
            if(abs(theta(i)-theta0)>par%px/2.) then
                dist(i)=0.;
            end if
        end do
        E0=par%rhog8*par%Hrms**2.;
    
        ! energy density distribution

        if (sum(dist)>0.d0) then 
           factor = (dist/sum(dist))/dtheta
        else
           factor=0.d0
        endif
        e01    = factor*E0;
        e01    = max(e01,0.0d0);

        if (abs(theta0)<1.e-9) theta0=1.e-9
        par%Llong=par%Tlong*cg(1,1)/sin(theta0)

    endif
	write(*,*)'Boundary conditions complete, starting computation'
end if

!write(*,*)'par%t =',par%t,'bcendtime = ',bcendtime
if (par%t .ge. bcendtime) then  ! Recalculate bcf-file 
    if (par%instat==4) then
        close(71)
        close(72)
        call makebcf(par,s,wp)
        startbcf=.true.
    elseif (par%instat==5) then 
        close(71)
        close(72)
        call makebcf(par,s,wp)
        startbcf=.true.
    elseif (par%instat==6) then 
        close(71)
        close(72)
        call makebcf(par,s,wp)
        startbcf=.true.
    elseif (par%instat==7) then
        close(71)
        close(72)
        startbcf=.true.
        if (par%t <= (par%tstop-par%dt)) then
            par%listline=par%listline+1
        end if
    end if
end if  

!
! COMPUTE WAVE BOUNDARY CONDITIONS CURRENT TIMESTEP
!
! instat = 0 => stationary wave conditions
! instat = 1 => wave energy fluctuations associated with Tlong in params.txt; bound long wave from LHS 1962
! instat = 2 => wave energy time series from Gen.ezs; bound long wave from LHS 1962
! instat = 3 => wave energy time series and (bound) long waves from Gen.ezs
! instat = 4 => directional wave energy time series for Jonswap spectrum from user specified file; bound long wave from van Dongeren, 19??
! instat = 5 => directional wave energy time series from SWAN 2D spectrum file; bound long wave from van Dongeren, 19??    
! instat = 6 => directional wave energy time series from spectrum file; bound long wave from van Dongeren, 19??
! instat = 7 => as instat = 4/5/6; reading from previously computed wave boundary condition file.
if (par%instat==0) then
   do j=1,ny+1
       ee(1,j,:)=e01*min(par%t/par%taper,1.)
       bi(1) = 0.d0
       ui(1,j) = 0.d0
   end do
elseif (par%instat==1) then
   do j=1,ny+1
       ee(1,j,:)=e01*0.5d0*(1.d0+cos(2.d0*par%px*(par%t/par%Tlong-sin(theta0)*y(1,j)/par%Llong))) *min(par%t/par%taper,1.d0);
       em = (sum(0.5d0*e01))*dtheta *min(par%t/par%taper,1.d0)
       ei =  sum(ee(1,j,1:ntheta))*dtheta
       bi(1) = -(2.d0*cg(1,j)/c(1,j)-0.5d0)*(em-ei)/(cg(1,j)**2.d0-par%g*hh(1,j))/par%rho
	   ht=s%zs0(1:2,:)-zb(1:2,:)
       ui(1,j) = cg(1,j)*bi(1)/ht(1,j)*cos(theta0)
   end do
elseif (par%instat==2) then
   do j=1,ny+1
      if (abs(theta0)<1e-3) then
         call linear_interp(tE,dataE,nt,par%t,E1,E_idx)
      else
         tshifted=max(par%t-(y(1,j)-y(1,1))*sin(theta0)/cg(1,1),0.)
         call linear_interp(tE,dataE,nt,tshifted,E1,E_idx) 
      endif
      ee(1,j,:)=e01*E1/max(par%Emean,0.000001d0)*min(par%t/par%taper,1.)
      em = par%Emean *min(par%t/par%taper,1.)
      ei = sum(ee(1,j,:))*dtheta
      bi(1) = -(2.*cg(1,j)/c(1,j)-0.5)*(em-ei)/(cg(1,j)**2.-par%g*hh(1,j))/par%rho
      ht=s%zs0(1:2,:)-zb(1:2,:)
      ui(1,j) = cg(1,j)*bi(1)/ht(1,j)*cos(theta0)
   end do
elseif (par%instat==3) then
    ht=s%zs0(1:2,:)-zb(1:2,:)
    do j=1,ny+1
        if (abs(theta0)<1e-3) then
           call linear_interp(tE,dataE,nt,par%t,E1,E_idx) 
           call linear_interp(tE,databi,nt,par%t,bi,E_idx)
        else
           tshifted=max(par%t-(y(1,j)-y(1,1))*sin(theta0)/cg(1,1),0.)
           call linear_interp(tE,dataE,nt,tshifted,E1,E_idx) 
           call linear_interp(tE,databi,nt,tshifted,bi(1),E_idx)
        endif
        ee(1,j,:)=e01*E1/max(par%Emean,0.000001d0)*min(par%t/par%taper,1.);
        ui(1,j) = cg(1,j)*bi(1)/ht(1,j)*cos(theta0)*min(par%t/par%taper,1.)
		if (par%carspan==1) then
        ui(1,j) = sqrt(par%g/ht(1,j))*bi(1)! Carrier and Greenspan
		end if

    end do
elseif ((par%instat==4).or.(par%instat==5) .or. (par%instat==6) .or. (par%instat==7)) then  
    ! open file if first time
    if (startbcf) then
        open(53,file='ebcflist.bcf',form='formatted',position='rewind')
        open(54,file='qbcflist.bcf',form='formatted',position='rewind')
!        write(*,*)'par%listline = ',par%listline
        do i=1,par%listline
            read(53,*)bcendtime,rt,dtbcfile,par%Trep,s%theta0,ebcfname
            read(54,*)bcendtime,rt,dtbcfile,par%Trep,s%theta0,qbcfname
        end do
        close(53)
        close(54)
        ! Robert and Jaap : Initialize for new wave conditions
        par%Trep = par%Trep
        par%omega = 2.d0*par%px/par%Trep            
        do itheta=1,s%ntheta
            s%sigt(:,:,itheta) = 2.*par%px/par%Trep;
        end do
        s%sigm = sum(s%sigt,3)/s%ntheta;
        call dispersion(par,s);     
        ! End initialize
        open(71,file=ebcfname,status='old',form='unformatted')
        open(72,file=qbcfname,status='old',form='unformatted')
        if (.not. allocated(q1) ) then
            allocate(q1(ny+1),q2(ny+1),q(ny+1))
            allocate(ee1(ny+1,ntheta),ee2(ny+1,ntheta))
        end if
!        write(*,*)'e1'
        read(71)ee1       ! Earlier in time
!        write(*,*)'e2'
        read(71)ee2       ! Later in time
!        write(*,*)'q1'
        read(72)q1        ! Earlier in time
!        write(*,*)'q2'
        read(72)q2        ! Later in time
        old=floor((par%t/dtbcfile)+1.)
    end if 

    new=floor((par%t/dtbcfile)+1.)

    ! Check for next level in boundary condition file
    if (new/=old) then
        ! Check for how many bcfile steps are jumped
        ! Difference = new-old, but as the files jump 1 anyway,
        ! only do loop for difference-1 times
        do i=1,(new-old)-1
!            write(*,*)i,'qi'
            read(72)q2
!            write(*,*)i,'ei'
            read(71)ee2
        end do
        ! Continue. Switch q2 to q1 and read new q2 
        q1=q2
        ee1=ee2
!        write(*,*)'read q'
        read(72)q2
!        write(*,*)'read e'
        read(71)ee2
        old=new
    end if
    ht=s%zs0(1:2,:)-zb(1:2,:)
    tnew = new*dtbcfile
    ee(1,:,:) = (dtbcfile-(tnew-par%t))/dtbcfile*ee2 + & !Jaap
                (tnew-par%t)/dtbcfile*ee1
    q = (dtbcfile-(tnew-par%t))/dtbcfile*q2 + &          !Jaap
        (tnew-par%t)/dtbcfile*q1
    ui(1,:) = q/ht(1,:)*min(par%t/par%taper,1.)
    ee(1,:,:)=ee(1,:,:)*min(par%t/par%taper,1.)
else
   write(*,*)' instat = ',par%instat, ' invalid option'
   stop
endif

if (par%t>0.) then
    !
    ! Lateral boundary at y=0;
    !
    if(par%instat>0) then
        fac1=(y(2:nx+1,2)-y(2:nx+1,1))*abs(tan(thetamean(2:nx+1,2)))/(x(2:nx+1,2)-x(1:nx,2));
    else
        fac1=0;
    end if
    fac1=min(fac1,1.);
    fac1=max(fac1,0.);
    fac2=1.-fac1;
    do itheta=1,ntheta 

          ee(2:nx+1,1,itheta)=ee(1:nx+1-1,2,itheta)*fac1+ee(2:nx+1,2,itheta)*fac2

    end do
    !
    ! lateral; boundary at y=ny*dy
    !
    if(par%instat>0) then
        fac1=(y(2:nx+1,ny+1)-y(2:nx+1,ny))*abs(tan(thetamean(2:nx+1,ny)))/(x(2:nx+1,ny)-x(1:nx,ny));
    else
        fac1=0.
    end if
    fac1=min(fac1,1.);
    fac1=max(fac1,0.);
    fac2=1.-fac1;
     do itheta=1,ntheta

           ee(2:nx+1,ny+1,itheta)= &
           ee(1:nx+1-1,ny+1-1,itheta)*fac1+ee(2:nx+1,ny+1-1,itheta)*fac2;

     end do
end if

end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! FLOW BOUNDARY CONDITIONS
!
subroutine flow_bc(s,par,it)
use params
use spaceparams

IMPLICIT NONE

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

integer                                     :: it
integer                                     :: i,ig
integer                                     :: j,jj,indt
real*8                                      :: omega,aux1,qxr,alphanew,vert
real*8 , dimension(:)    ,allocatable       :: yzs0,szs0
real*8 , dimension(:,:)  ,allocatable,save  :: zs0old
real*8 , dimension(:,:)  ,allocatable,save  :: ht,beta,betanp1
real*8 , dimension(:)    ,allocatable,save  :: bn,alpha2
real*8 , dimension(:,:)  ,allocatable,save  :: dhdx,dhdy,dvdx,dvdy
real*8 , dimension(:,:)  ,allocatable,save  :: dbetadx,dbetady,dvudy
real*8 , dimension(:)    ,allocatable,save  :: inv_ht

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

if (.not. allocated(ht)) then
   allocate(ht      (2,ny+1))
   allocate(dhdx    (2,ny+1))
   allocate(dhdy    (2,ny+1))
   allocate(dvdx    (2,ny+1))
   allocate(dvdy    (2,ny+1))
   allocate(dvudy   (1,ny+1))
   allocate(inv_ht  (ny+1))
   allocate(dbetadx (2,ny+1))
   allocate(dbetady (2,ny+1))
   allocate(beta    (2,ny+1))
   allocate(bn      (ny+1))
   allocate(alpha2  (ny+1))
   allocate(betanp1 (1,ny+1))
   allocate(zs0old(nx+1,ny+1))
endif

allocate(szs0(1:2)) 
allocate(yzs0(1:2)) 
!
! Sea boundary at x=0;
!
!
! UPDATE TIDE AND SURGE
! 
! Need to interpolate input tidal signal to xbeach par%t to 
! compute proper tide contribution 
!write(*,*)  'made it into flow_bc and about to interp tide(t)'

if (par%tideloc>0) then

  call LINEAR_INTERP(s%tideinpt, s%tideinpz(:,1), par%tidelen, par%t, par%zs01, indt)

  if(par%tideloc.eq.1) then 
  par%zs02=par%zs01

  end if

  if(par%tideloc.eq.2 .and. par%paulrevere.eq.0) then
    call LINEAR_INTERP(s%tideinpt, s%tideinpz(:,2), par%tidelen, par%t, par%zs03, indt)
    par%zs02=par%zs01
    par%zs04=par%zs03
  endif

  if(par%tideloc.eq.2 .and. par%paulrevere.eq.1) then
    call LINEAR_INTERP(s%tideinpt, s%tideinpz(:,2), par%tidelen, par%t, par%zs02, indt)
    par%zs03=0.d0
    par%zs04=0.d0
  endif

  if(par%tideloc.eq.4) then
    call LINEAR_INTERP(s%tideinpt, s%tideinpz(:,2), par%tidelen, par%t, par%zs02, indt)
    call LINEAR_INTERP(s%tideinpt, s%tideinpz(:,3), par%tidelen, par%t, par%zs03, indt)
    call LINEAR_INTERP(s%tideinpt, s%tideinpz(:,4), par%tidelen, par%t, par%zs04, indt)
  endif

  if(par%tideloc.eq.1) s%zs0 = par%zs01

  if(par%tideloc.eq.2 .and. par%paulrevere.eq.1) then
    yzs0(1)=s%yz(1)
    yzs0(2)=s%yz(s%ny+1)
    szs0(1)=par%zs01
    szs0(2)=par%zs02
    do i = 1,s%ny+1
        call LINEAR_INTERP(yzs0, szs0, 2, s%yz(i), s%zs0(1,i), indt)
    enddo
    do j = 1,s%ny+1 
        do i = 1,s%nx+1
            s%zs0(i,j) = s%zs0(1,j)
        enddo
    enddo   
  endif

  if(par%tideloc.eq.2 .and. par%paulrevere.eq.0) then
    yzs0(1)=s%xz(1)
    yzs0(2)=s%xz(s%nx+1)
    szs0(1)=par%zs01
    szs0(2)=par%zs04
    s%zs0(1,:)=par%zs01
    s%zs0(s%nx+1,:)=par%zs03
    do j = 1,s%ny+1 
        do i = 1,s%nx+1
            if (s%zb(i,j).gt.s%zs0(1,j)+par%eps) goto 302
            s%zs0(i,j) = s%zs0(1,j)
        enddo
302     if (i.lt.s%nx+1) then
            do ig=i+1,s%nx+1
                s%zs0(ig,j) = s%zs0(s%nx+1,j)   
            enddo
        else
            do i = 1,s%nx+1
                call LINEAR_INTERP(yzs0, szs0, 2, s%xz(i), s%zs0(i,j), indt)
            enddo
        endif
    enddo
  endif

  if(par%tideloc.eq.4) then
    yzs0(1)=s%yz(1)
    yzs0(2)=s%yz(s%ny+1)
    szs0(1)=par%zs01
    szs0(2)=par%zs02
    do i = 1,s%ny+1
        call LINEAR_INTERP(yzs0, szs0, 2, s%yz(i), s%zs0(1,i), indt)
    enddo
    yzs0(1)=s%yz(1)
    yzs0(2)=s%yz(s%ny+1)
    szs0(1)=par%zs04
    szs0(2)=par%zs03
    do i = 1,s%ny+1
        call LINEAR_INTERP(yzs0, szs0, 2, s%yz(i), s%zs0(s%nx+1,i), indt)
    enddo

    do j = 1,s%ny+1 
        do i = 1,s%nx+1
            if (s%zb(i,j).gt.s%zs0(1,j)+par%eps) goto 303
            s%zs0(i,j) = s%zs0(1,j)
        enddo
        
303     if (i.lt.s%nx+1) then
            do ig=i+1,s%nx+1
                s%zs0(ig,j) = s%zs0(s%nx+1,j)   
            enddo
        else
            do i = 1,s%nx+1
                call LINEAR_INTERP(yzs0, szs0, 2, s%xz(i), s%zs0(i,j), indt)
            enddo
        endif
    enddo
  endif

endif
!
! UPDATE (LONG) WAVES
!
if (par%instat<8)then
if (par%front==0) then ! Ad's radiating boundary
   uu(1,:)=2.d0*ui(1,:)-(sqrt(par%g/hh(1,:))*(zs(2,:)-s%zs0(2,:)))
   vv(1,:)=vv(2,:);
elseif (par%front==1) then ! Van Dongeren (1997), weakly reflective boundary condition
   ht(1:2,:)=s%zs0(1:2,:)-zb(1:2,:)
   beta=uu(1:2,:)-2.*dsqrt(par%g*hum(1:2,:)) !cjaap : replace hh with hum

   do j=2,ny
        ! compute gradients in u-points....
        dvdy(1,j)=(vu(1,min(j,ny)+1)-vu(1,max(j,2)-1))/(yz(min(j,ny)+1)-yz(max(j,2)-1))
        dhdx(1,j)=(ht(2,j)-ht(1,j))/(xz(2)-xz(1))
        dbetadx(1,j)=(beta(2,j)-beta(1,j))/(xu(2)-xu(1))
        dbetady(1,j)=(beta(1,min(j,ny)+1)-beta(1,max(j,2)-1))/(yz(min(j,ny)+1)-yz(max(j,2)-1))

        inv_ht(j) = 1.d0/hum(1,j)                                  !Jaap replaced hh with hum

        bn(j)=-(uu(1,j)-dsqrt(par%g*hum(1,j)))*dbetadx(1,j) &      !Jaap replaced hh with hum 
              -vu(1,j)*dbetady(1,j)& !Ap vu
              +dsqrt(par%g*hum(1,j))*dvdy(1,j)&                    !Jaap replaced hh with hum 
              +Fx(1,j)*inv_ht(j)/par%rho-par%g/par%C**2.d0&        !Ap
              *sqrt(uu(1,j)**2+vu(1,j)**2)*uu(1,j)/hum(1,j)&       !Jaap replaced hh with hum
              +par%g*dhdx(1,j)
   end do

   do j=2,ny
      betanp1(1,j) = beta(1,j)+ bn(j)*par%dt
      alpha2(j)=-theta0
      alphanew = 0.d0
      s%umean(1,j) = (par%epsi*uu(1,j)+(1-par%epsi)*s%umean(1,j))  
      do jj=1,50
         !---------- Lower order bound. cond. ---
         qxr = dcos(alpha2(j))/(dcos(alpha2(j))+1.d0)&
           *(0.5*(ht(1,j)+ht(2,j))*(betanp1(1,j)-s%umean(1,j)+2.d0*DSQRT(par%g*0.5*(ht(1,j)+ht(2,j))))&  !Jaap replaced ht(1,j) with 0.5*(ht(1,j)+ht(2,j))
           -(ui(1,j)*hum(1,j))*(dcos(theta0)-1.d0)/dcos(theta0))   !Jaap replaced hh with hu

         !vert = velocity of the reflected wave = total-specified
         vert = vu(1,j)-ui(1,j)*tan(theta0)
         alphanew = datan(vert*hh(1,j)/(qxr+1.d-16))         
         if (alphanew .gt. (par%px*0.5d0)) alphanew=alphanew-par%px
         if (alphanew .le. (-par%px*0.5d0)) alphanew=alphanew+par%px
         if(dabs(alphanew-alpha2(j)).lt.0.001) goto 1000
         alpha2(j) = alphanew 
      end do
1000  continue
      if (par%ARC==0) then
         uu(1,j) = (par%order-1)*ui(1,j)
         zs(1,j) = zs(2,j)
      else
         uu(1,j) = (par%order-1.)*ui(1,j)+2.*qxr/(ht(1,j)+ht(2,j)) + s%umean(1,j) !Jaap: replaced ht(1,j) with 0.5*(ht(1,j)+ht(2,j))
         !with a taylor expansion to get to the zs point at index 1 from uu(1) and uu(2)
	     zs(1,j) = 1.5*((betanp1(1,j)-uu(1,j))**2/4./par%g+.5*(zb(1,j)+zb(2,j)))- &
	               0.5*((beta(2,j)-uu(2,j))**2/4./par%g+.5*(zb(2,j)+zb(3,j)))    
	  end if
   end do
endif
!
! Radiating boundary at x=nx*dx
!
if (par%back==1) then !solve water level at nx+1 from continuity and set uu(nx+1,:)=0
    uu(nx,:) = 0.d0;
    zs(nx+1,:) = zs(nx,:);
    !zs(nx+1,2:ny) = zs(nx+1,2:ny) + par%dt*hh(nx,2:ny)*uu(nx,2:ny)/(xu(nx+1)-xu(nx)) -par%dt*(hv(nx+1,2:ny+1)*vv(nx+1,2:ny+1)-hv(nx+1,1:ny)*vv(nx+1,1:ny))/(yv(2:ny+1)-yv(1:ny))
elseif (par%back==0) then
    ! uu(nx+1,:)=sqrt(par%g/hh(nx+1,:))*(zs(nx+1,:)-max(zb(nx+1,:),s%zs0(nx+1,:))) ! cjaap: make sure if the last cell is dry no radiating flow is computed... 
    !uu(nx+1,:)=sqrt(par%g/(s%zs0(nx+1,:)-zb(nx+1,:)))*(zs(nx+1,:)-max(zb(nx+1,:),s%zs0(nx+1,:)))
    s%umean(2,:) = par%epsi*uu(nx,:)+(1-par%epsi)*s%umean(2,:)    !Ap
    zs(nx+1,:)=max(s%zs0(nx+1,:),s%zb(nx+1,:))+(uu(nx,:)-s%umean(2,:))*sqrt(max((s%zs0(nx+1,:)-zb(nx+1,:)),par%eps)/par%g)    !Ap
elseif (par%back==2) then
    ht(1:2,:)=s%zs0(s%nx:s%nx+1,:)-zb(s%nx:s%nx+1,:)
    beta=uu(s%nx-1:s%nx,:)+2.*dsqrt(par%g*hum(s%nx-1:s%nx,:)) !cjaap : replace hh with hum

   do j=2,ny
      if (wetu(s%nx,j)==1) then   ! Robert: dry back boundary points
        ! compute gradients in u-points....
        dvdy(2,j)=(vu(s%nx,min(j,ny)+1)-vu(s%nx,max(j,2)-1))/(yz(min(j,ny)+1)-yz(max(j,2)-1))
        dhdx(2,j)=(ht(2,j)-ht(1,j))/(xz(s%nx+1)-xz(s%nx))
        dbetadx(2,j)=(beta(2,j)-beta(1,j))/(xu(s%nx)-xu(s%nx-1))
        dbetady(2,j)=(beta(1,min(j,ny)+1)-beta(1,max(j,2)-1))/(yz(min(j,ny)+1)-yz(max(j,2)-1))

        inv_ht(j) = 1.d0/hum(s%nx,j)                                           !Jaap replaced hh with hum

        bn(j)=-(uu(s%nx,j)+dsqrt(par%g*hum(s%nx,j)))*dbetadx(2,j) &            !Ap says plus  !Jaap replaced hh with hum 
              -vu(s%nx,j)*dbetady(2,j)& !Ap vu
              -dsqrt(par%g*hum(s%nx,j))*dvdy(2,j)&                             !Jaap replaced hh with hum 
              +Fx(s%nx,j)*inv_ht(j)/par%rho-par%g/par%C**2.d0&                 !Ap
              *sqrt(uu(s%nx,j)**2+vu(s%nx,j)**2)*uu(s%nx,j)/hum(s%nx,j)&       !Jaap replaced hh with hum
              +par%g*dhdx(2,j)
      endif   ! Robert: dry back boundary points
   enddo

   do j=2,ny
     if (wetu(s%nx,j)==1) then                                                   ! Robert: dry back boundary points
      betanp1(1,j) = beta(2,j)+ bn(j)*par%dt                                   !Ap toch?
      alpha2(j)= theta0
      alphanew = 0.d0
      s%umean(2,j) = (par%epsi*uu(s%nx,j)+(1-par%epsi)*s%umean(2,j))           !Ap 
      do jj=1,50
         !---------- Lower order bound. cond. ---
         qxr = dcos(alpha2(j))/(dcos(alpha2(j))+1.d0)&  
           *(0.5*(ht(1,j)+ht(2,j))*(betanp1(1,j)-s%umean(2,j)-2.d0*DSQRT(par%g*0.5*(ht(1,j)+ht(2,j))))) 

         !vert = velocity of the reflected wave = total-specified
         vert = vu(s%nx,j)
         alphanew = datan(vert*hh(s%nx+1,j)/(qxr+1.d-16))                      !Ap
         if (alphanew .gt. (par%px*0.5d0)) alphanew=alphanew-par%px
         if (alphanew .le. (-par%px*0.5d0)) alphanew=alphanew+par%px
         if(dabs(alphanew-alpha2(j)).lt.0.001) goto 2000
         alpha2(j) = alphanew 
      end do
2000  continue
	  uu(s%nx,j) = 2.*qxr/(ht(1,j)+ht(2,j)) + s%umean(2,j)                       !Jaap: replaced ht(1,j) with 0.5*(ht(1,j)+ht(2,j))
	  ! Ap replaced zs with extrapolation.
      zs(s%nx+1,j) = 1.5*((betanp1(1,j)-uu(s%nx,j))**2.d0/4.d0/par%g+.5*(zb(s%nx,j)+zb(s%nx+1,j)))-&
				     0.5*((beta(1,j)-uu(s%nx-1,j))**2.d0/4.d0/par%g+.5*(zb(s%nx-1,j)+zb(s%nx,j)))
   
     endif   ! Robert: dry back boundary points
   enddo
endif
endif 

end subroutine