subroutine cofrbc(j ,nmmaxj ,norow ,icx ,icy , &
& urf ,hu ,dpu ,umean ,vmean , &
& crbc ,ctbf ,ctbl ,ctif ,ctil , &
& ctrf ,ctrl ,stbf ,stbl ,stif , &
& stil ,cgdghf ,cgdghl ,zetabf ,zetabl , &
& zetaif ,zetail ,zmeanf ,zmeanl ,umeanf , &
& umeanl ,circ2d ,gvu ,wsu ,irocol , &
& timsec ,hdt ,gdp )
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2011-2014.
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation version 3.
!
! This program 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 General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id$
! $HeadURL$
!!--description-----------------------------------------------------------------
! NONE
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
use mathconsts
use globaldata
!
implicit none
!
type(globdat),target :: gdp
!
! The following list of pointer parameters is used to point inside the gdp structure
!
integer , pointer :: nmskf
integer , pointer :: nmskl
integer , pointer :: mmskf
integer , pointer :: mmskl
real(fp) , pointer :: rhow
real(fp) , pointer :: ag
integer , pointer :: iro
!
! Global variables
!
integer :: icx
integer :: icy
integer :: j
integer :: nmmaxj ! Description and declaration in dimens.igs
integer :: norow ! Description and declaration in esm_alloc_int.f90
integer, dimension(5, norow) :: irocol ! Description and declaration in esm_alloc_int.f90
real(fp) :: hdt ! Description and declaration in esm_alloc_real.f90
real(fp) :: timsec ! Description and declaration in inttim.igs
real(fp) :: urf
real(fp), dimension(norow) :: cgdghf ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: cgdghl ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: ctbf ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: ctbl ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: ctif ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: ctil ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: ctrf ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: ctrl ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: stbf ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: stbl ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: stif ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: stil ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: umeanf ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: umeanl ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: zetabf ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: zetabl ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: zetaif ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: zetail ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: zmeanf ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(norow) :: zmeanl ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(4, norow) :: circ2d ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(12, norow) :: crbc ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: dpu ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: gvu ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: hu ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: umean ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: vmean ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: wsu ! Description and declaration in esm_alloc_real.f90
!
! Local variables
!
integer :: ddb
integer :: ibf
integer :: ibl
integer :: ic
integer :: icxy
integer :: mf
integer :: ml
integer :: n
integer :: nmf
integer :: nmfu
integer :: nmfuu
integer :: nml
integer :: nmld
integer :: nmldd
real(fp):: alfas
real(fp):: ctr
real(fp):: ctr0
real(fp):: depth
real(fp):: dt
real(fp):: f
real(fp):: fdenom
real(fp):: rmux
real(fp):: sepu
real(fp):: tapf
real(fp):: tau
real(fp):: temp1
real(fp):: temp2
real(fp):: temp3
real(fp):: temp4
real(fp):: urfa
real(fp):: zbfic
real(fp):: zblic
!
!! executable statements -------------------------------------------------------
!
nmskf => gdp%gdbcdat%nmskf
nmskl => gdp%gdbcdat%nmskl
mmskf => gdp%gdbcdat%mmskf
mmskl => gdp%gdbcdat%mmskl
rhow => gdp%gdphysco%rhow
ag => gdp%gdphysco%ag
iro => gdp%gdphysco%iro
!
ctr0 = cos(8.0*pi/18.0)
dt = 2.0*hdt
!
! First smoothen the angles of the reflected outgoing waves
!
if (timsec<4*hdt) then
do ic = 1, norow
ctrf(ic) = 1.0
ctrl(ic) = 1.0
enddo
endif
if (norow<2) then
temp1 = 1.0
temp3 = 1.0
else
temp1 = 0.5*(ctrf(1) + ctrf(2))
temp3 = 0.5*(ctrl(1) + ctrl(2))
endif
do ic = 2, norow - 1
temp2 = 0.25*(ctrf(ic - 1) + 2.0*ctrf(ic) + ctrf(ic + 1))
temp4 = 0.25*(ctrl(ic - 1) + 2.0*ctrl(ic) + ctrl(ic + 1))
ctrf(ic - 1) = temp1
ctrl(ic - 1) = temp3
temp1 = temp2
temp3 = temp4
enddo
if (norow>1) then
ctrf(norow) = 0.5*(ctrf(norow - 1) + ctrf(norow))
ctrl(norow) = 0.5*(ctrl(norow - 1) + ctrl(norow))
ctrf(norow - 1) = temp1
ctrl(norow - 1) = temp3
endif
!
icxy = max(icx, icy)
ddb = gdp%d%ddbound
!
do ic = 1, norow
n = irocol(1, ic)
mf = irocol(2, ic) - 1
ml = irocol(3, ic)
ibf = irocol(4, ic)
ibl = irocol(5, ic)
nmf = (n + ddb)*icy + (mf + ddb)*icx - icxy
nmfu = nmf + icx
nmfuu = nmfu + icx
nml = (n + ddb)*icy + (ml + ddb)*icx - icxy
nmld = nml - icx
nmldd = nmld - icx
depth = hu(nmf)
if ((ibf==6 .or. ibf==2).and.depth>1.0e-6) then
!
! WEAKLY REFLECTIVE BOUNDARY CONDITION AT LEFT BOUNDARY
!
sepu = hu(nmf) - dpu(nmf)
urfa = urf
!
! reduction of locked incoming wave due to mask effects
! (nearly) equidistant grid assumed
!
zbfic = zetabf(ic)
!
! left BC
!
if (icy==1) then
if (mmskf>0) then
zbfic = 0.0
else
if (nmskf>0 .and. ic0 .and. ic>nmskl) then
tau = real(ic - nmskl, fp)/real(icx - nmskl, fp)
zbfic = zbfic*tapf(tau)
endif
endif
endif
!
! bottom BC
!
if (icx==1) then
if (nmskf>0) then
zbfic = 0.0
else
if (mmskf>0 .and. ic0 .and. ic>mmskl) then
tau = real(ic - mmskl, fp)/real(icy - mmskl, fp)
zbfic = zbfic*tapf(tau)
endif
endif
endif
!
alfas = sqrt(ag/depth)
fdenom = umean(nmf) - alfas*(sepu - zetaif(ic)*(1.0 - ctif(ic)) &
& - zbfic*(1.0 - cgdghf(ic)*ctbf(ic)))
if (abs(fdenom)>1.0E-4) then
f = (0.25*(3.0*vmean(nmfu) - vmean(nmfuu) + 3.0*vmean(nmfu - icy) &
& - vmean(nmfuu - icy)) - alfas*stif(ic)*zetaif(ic) &
& - alfas*cgdghf(ic)*stbf(ic)*zbfic)/fdenom
if (f*f<=1.0) then
ctr = (1.0 - f*f)/(1.0 + f*f)
else
ctr = 1.0
urfa = 0.0
endif
ctr = max(ctr0, ctr)
ctrf(ic) = ctrf(ic)*(1.0 - urfa) + urfa*ctr
endif
!
! In the case of tidal computations the mean waterlevel at
! the sea boundary can be prescribed through circ2d
!
if (ibf==2) then
zmeanf(ic) = circ2d(1, ic)
umeanf(ic) = 0.0
endif
crbc(1, ic) = 0.5*alfas*ctrf(ic)
crbc(2, ic) = 0.5*alfas*ctrf(ic)
crbc(3, ic) = umeanf(ic) + alfas*(zmeanf(ic)*ctrf(ic) &
& + zetaif(ic)*(ctrf(ic) + ctif(ic)) &
& + zbfic*(ctrf(ic) + cgdghf(ic)*ctbf(ic)))
!
! LAYER VELOCITIES ( VELOCITY PROFILE )
!
rmux = sqrt(ag*depth)*dt/gvu(nmf)
crbc(4, ic) = ctrf(ic) + rmux
crbc(5, ic) = ctrf(ic) - rmux
crbc(6, ic) = 2.0*zetaif(ic)*(ctrf(ic) + ctif(ic)) &
& + 2.0*zbfic*(ctrf(ic) + cgdghf(ic)*ctbf(ic)) &
& + 2.0*zmeanf(ic)*ctrf(ic) &
& - (2.0*(umean(nmf) - umeanf(ic)) + dt*wsu(nmf) &
& /rhow/depth)/alfas
endif
depth = hu(nml)
if ((ibl==6 .or. ibl==2) .and. depth>1.0e-6) then
!
! WEAKLY REFLECTIVE BOUNDARY CONDITION AT RIGHT BOUNDARY
!
sepu = hu(nml) - dpu(nml)
urfa = urf
!
! reduction of locked incoming wave due to mask strips
! (nearly) equidistant grid assumed
!
zblic = zetabl(ic)
!
! right BC
!
if (icy==1) then
if (mmskl>0) then
zblic = 0.0
else
if (nmskf>0 .and. ic0 .and. ic>nmskl) then
tau = real(ic - nmskl, fp)/real(icx - nmskl, fp)
zblic = zblic*tapf(tau)
endif
endif
endif
!
! top BC
!
if (icx==1) then
if (nmskl>0) then
zblic = 0.0
else
if (mmskf>0 .and. ic0 .and. ic>mmskl) then
tau = real(ic - mmskl, fp)/real(icy - mmskl, fp)
zblic = zblic*tapf(tau)
endif
endif
endif
!
alfas = sqrt(ag/depth)
fdenom = umean(nml) + alfas*(sepu - zetail(ic)*(1.0 - ctil(ic)) &
& - zblic*(1.0 - cgdghl(ic)*ctbl(ic)))
if (abs(fdenom)>1.0E-4) then
f = (0.25*(3.0*vmean(nml) - vmean(nmld) + 3.0*vmean(nml - icy) &
& - vmean(nmld - icy)) - alfas*stil(ic)*zetail(ic) &
& - alfas*cgdghl(ic)*stbl(ic)*zblic)/fdenom
if (f*f<=1.0) then
ctr = (1.0 - f*f)/(1.0 + f*f)
else
ctr = 1.0
urfa = 0.0
endif
ctr = max(ctr0, ctr)
ctrl(ic) = ctrl(ic)*(1.0 - urfa) + urfa*ctr
endif
!
! In the case of tidal computations the mean waterlevel at
! the sea boundary can be prescribed through circ2d
!
if (ibl==2) then
zmeanl(ic) = circ2d(2, ic)
umeanl(ic) = 0.0
endif
!
crbc(7, ic) = -0.5*alfas*ctrl(ic)
crbc(8, ic) = -0.5*alfas*ctrl(ic)
crbc(9, ic) = umeanl(ic) - alfas*(zmeanl(ic)*ctrl(ic) &
& + zetail(ic)*(ctrl(ic) + ctil(ic)) &
& + zblic*(ctrl(ic) + cgdghl(ic)*ctbl(ic)))
!
rmux = sqrt(ag*depth)*dt/gvu(nml)
crbc(10, ic) = ctrl(ic) - rmux
crbc(11, ic) = ctrl(ic) + rmux
crbc(12, ic) = 2.0*zetail(ic)*(ctrl(ic) + ctil(ic)) &
& + 2.0*zblic*(ctrl(ic) + cgdghl(ic)*ctbl(ic)) &
& + 2.0*zmeanl(ic)*ctrl(ic) &
& + (2.0*(umean(nml) - umeanl(ic)) + dt*wsu(nml) &
& /rhow/depth)/alfas
endif
enddo
end subroutine cofrbc