subroutine dengra(icreep ,j ,nmmaxj ,nmmax ,kmax , &
& icx ,icy ,lstsci ,lsts ,lsal , &
& ltem ,lsed ,saleqs ,temeqs ,rhosol , &
& kcs ,kfu ,s0 ,dps ,hu , &
& thick ,sig ,guu ,gvu ,r0 , &
& dicuv ,dpu ,dpdksi ,dsdksi ,dtdksi , &
& dldksi ,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-----------------------------------------------------------------
!
! Function: Computes horizontal pressure gradient due
! to density gradients, along strictly
! horizontal planes, using a special limiter,
! to avoid 'artificial flow'. Recently the
! limiter was improved. The original minimum
! limiter leads to underestimation of the
! pressure gradient. The limiter is following
! Van Leer and preserves monotony.
! Horizontal diffusive fluxes of salt and
! temperature are stored in arrays for
! for transport equation.
! Also in case of sediment; but only in baroclinic
! density term (no correction for diffusive flux
! for sediment; i.e. water is assumed homogeneous
! HOWEVER array DLDKSI for this purpose is already
! RESERVED!)
!
! Method used: Reference : On the approximation of horizontal
! gradients in sigma co-ordinates for bathymetry
! with steep bottom slopes (G.S. Stelling and J.
! van Kester - International Journal for Methods
! in Fluids, Vol. 18 1994) and the discussion
! in "A comparison of two 3D shallow
! water models using sigma-coordinates and
! z-coordinates in the vertical direction."
! (M.D. Bijvelds, J. van Kester and G.S. Stelling -
! Proceedings 6-th International Conference on
! estuarine and coastal modelling, New Orleans,
! september 1999)
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
!
use globaldata
!
implicit none
!
type(globdat),target :: gdp
!
! The following list of pointer parameters is used to point inside the gdp structure
!
include 'pardef.igd'
real(fp) , pointer :: rhow
real(fp) , pointer :: ag
integer , pointer :: idensform
!
! Global variables
!
integer , intent(in) :: icreep ! Description and declaration in tricom.igs
integer , intent(in) :: icx !! Increment in the x-dir., if icx= nmax
!! then computation proceeds in the x-
!! dir. if icx=1 then computation pro-
!! ceeds in the y-dir.
integer :: icy !! Increment in the y-dir. (see icx)
integer :: j !! Begin pointer for arrays which have
!! been transformed into 1d arrays.
!! due to the shift in the 2nd (m-)
!! index, j = -2*nmax + 1
integer , intent(in) :: kmax ! Description and declaration in esm_alloc_int.f90
integer , intent(in) :: lsal ! Description and declaration in dimens.igs
integer , intent(in) :: lsed ! Description and declaration in esm_alloc_int.f90
integer , intent(in) :: lsts ! Description and declaration in dimens.igs
integer , intent(in) :: lstsci ! Description and declaration in esm_alloc_int.f90
integer , intent(in) :: ltem ! Description and declaration in dimens.igs
integer , intent(in) :: nmmax ! Description and declaration in dimens.igs
integer :: nmmaxj ! Description and declaration in dimens.igs
integer, dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kcs ! Description and declaration in esm_alloc_int.f90
integer, dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfu ! Description and declaration in esm_alloc_int.f90
real(fp) , intent(in) :: saleqs ! Description and declaration in tricom.igs
real(fp) , intent(in) :: temeqs ! Description and declaration in tricom.igs
real(prec), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: dps ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: dpu ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: guu ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gvu ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: hu ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: s0 ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax + 2) , intent(in) :: dicuv ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(out) :: dldksi ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: dpdksi ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: dsdksi ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: dtdksi ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), intent(in) :: r0 ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(kmax) , intent(in) :: sig ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(kmax) , intent(in) :: thick ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(lsed) , intent(in) :: rhosol ! Description and declaration in esm_alloc_real.f90
!
! Local variables
!
integer :: k
integer :: k1
integer :: k2
integer :: kf
integer :: kflux
integer :: kll
integer :: kpoint
integer :: krr
integer :: l ! Constituent number
integer :: lst
integer :: nm ! Index gridpoint
integer :: nmu ! Index neighbour (M+1)
integer, dimension(2*mxkmax + 1) :: kicol ! K-index concentration point left for flux kf
integer, dimension(2*mxkmax + 1) :: kicor ! K-index concentration point right for flux kf
real(fp) :: alph0
real(fp) :: cl
real(fp) :: cladt
real(fp) :: cr
real(fp) :: darea
real(fp) :: difco
real(fp) :: dpnm
real(fp) :: dpnmu
real(fp) :: drho
real(fp) :: dsal
real(fp) :: dtem
real(fp) :: dummy
real(fp) :: flux
real(fp) :: grad
real(fp) :: grad1
real(fp) :: grad2
real(fp) :: grmax
real(fp) :: grmin
real(fp) :: h0
real(fp) :: rhods
real(fp) :: rhodt
real(fp) :: sal
real(fp) :: sepnm
real(fp) :: sepnmu
real(fp) :: temp
real(fp) :: zubot
real(fp) :: zuk
real(fp), dimension(0:2*mxkmax + 1) :: point ! Merge arrays polal and polar
real(fp), dimension(0:mxkmax) :: polal ! Z-coordinate horizontal layers in nm
real(fp), dimension(0:mxkmax) :: polar ! Z-coordinate horizontal layers in nmu
real(fp), dimension(2*mxkmax + 1) :: poflu ! Z-coordinate gradient flux
real(fp), dimension(mxkmax) :: pocol ! Z-coordinate concentr. point nm ,k
real(fp), dimension(mxkmax) :: pocor ! Z-coordinate concentr. point nmu,k
!
!! executable statements -------------------------------------------------------
!
rhow => gdp%gdphysco%rhow
ag => gdp%gdphysco%ag
idensform => gdp%gdphysco%idensform
!
alph0 = 0.6980
!
! Initialise for all (nm, k)
!
dldksi = 0.0
dpdksi = 0.0
dsdksi = 0.0
dtdksi = 0.0
!
! pressure gradient; only for salinity and/or temperature and/or
! sediment and for anti creep
!
if (lsts==0 .or. icreep==0) goto 9999
!
nmu = icx
do nm = 1, nmmax
nmu = nmu + 1
!
! In the SED version the following line read:
! if (kfu(nm).eq.1 .and. kcs(nm)*kcs(nmu).ne.2) then
! The latter check was not introduced here on ground of
! backward compatibility for non-sed quantities.
!
if (kfu(nm)==1) then
sepnm = s0(nm)
sepnmu = s0(nmu)
dpnm = real(dps(nm) ,fp)
dpnmu = real(dps(nmu),fp)
!
! position horizontal interfaces left
!
polal(0) = sepnm
h0 = max(sepnm + dpnm, 0.01_fp)
do k = 1, kmax
polal(k) = (sig(k) - 0.5*thick(k))*h0 + sepnm
pocol(k) = 0.5*(polal(k - 1) + polal(k))
enddo
!
! position horizontal interfaces right
!
polar(0) = sepnmu
h0 = max(sepnmu + dpnmu, 0.01_fp)
do k = 1, kmax
polar(k) = (sig(k) - 0.5*thick(k))*h0 + sepnmu
pocor(k) = 0.5*(polar(k - 1) + polar(k))
enddo
!
! merge polal and polar
!
kll = 0
krr = 0
do k = 0, 2*kmax + 1
if (polal(kll)>polar(krr)) then
point(k) = polal(kll)
kll = kll + 1
if (kll>kmax) then
kpoint = k + 1
point(kpoint) = polar(krr)
goto 140
endif
else
point(k) = polar(krr)
krr = krr + 1
if (krr>kmax) then
kpoint = k + 1
point(kpoint) = polal(kll)
goto 140
endif
endif
enddo
kpoint = 2*kmax + 1
!
! position flux points
!
140 continue
kflux = kpoint
do k = 1, kflux
poflu(k) = 0.5*(point(k) + point(k - 1))
enddo
!
! k-index concentration points left and right for flux point
!
kll = 1
krr = 1
do kf = 1, kflux
kicol(kf) = 0
kicor(kf) = 0
do k = kll, kmax
if (poflu(kf)<=polal(k - 1) .and. poflu(kf)>=polal(k)) then
kicol(kf) = k
kll = k
exit
endif
enddo
do k = krr, kmax
if (poflu(kf)<=polar(k - 1) .and. poflu(kf)>=polar(k)) then
kicor(kf) = k
krr = k
exit
endif
enddo
enddo
!
! computation concentration flux
! using limiter
!
lst = max(lsal, ltem)
do l = 1, lsts
!
! computation of flux
!
do kf = 1, kflux
kll = kicol(kf)
krr = kicor(kf)
if (kll/=0 .and. krr/=0) then
!
! interpolation left point
!
if (pocor(krr)>pocol(kll)) then
k1 = 0
do k = kll - 1, 1, -1
if (pocol(k)>pocor(krr)) then
k1 = k
exit
endif
enddo
k2 = k1 + 1
if (k1<1) then
cl = r0(nm, k2, l)
else
cl = ((pocor(krr) - pocol(k2))/(pocol(k1) - pocol(k2)))&
& *r0(nm, k1, l) &
& + ((pocor(krr) - pocol(k1))/(pocol(k2) - pocol(k1)&
& ))*r0(nm, k2, l)
endif
else
k1 = kmax + 1
do k = kll + 1, kmax, 1
if (pocol(k)kmax) then
cl = r0(nm, k2, l)
else
cl = ((pocor(krr) - pocol(k2))/(pocol(k1) - pocol(k2)))&
& *r0(nm, k1, l) &
& + ((pocor(krr) - pocol(k1))/(pocol(k2) - pocol(k1)&
& ))*r0(nm, k2, l)
endif
endif
!
! interpolation right point
!
if (pocol(kll)>pocor(krr)) then
k1 = 0
do k = krr - 1, 1, -1
if (pocor(k)>pocol(kll)) then
k1 = k
exit
endif
enddo
k2 = k1 + 1
if (k1<1) then
cr = r0(nmu, k2, l)
else
cr = ((pocol(kll) - pocor(k2))/(pocor(k1) - pocor(k2)))&
& *r0(nmu, k1, l) &
& + ((pocol(kll) - pocor(k1))/(pocor(k2) - pocor(k1)&
& ))*r0(nmu, k2, l)
endif
else
k1 = kmax + 1
do k = krr + 1, kmax, 1
if (pocor(k)kmax) then
cr = r0(nmu, k2, l)
else
cr = ((pocol(kll) - pocor(k2))/(pocor(k1) - pocor(k2)))&
& *r0(nmu, k1, l) &
& + ((pocol(kll) - pocor(k1))/(pocor(k2) - pocor(k1)&
& ))*r0(nmu, k2, l)
endif
endif
grad1 = r0(nmu, krr, l) - cl
grad2 = cr - r0(nm, kll, l)
grmax = max(grad1, grad2)
grmin = min(grad1, grad2)
!
! new limiter following Van Leer
!
if (grmax>=0.0 .and. grmin<=0.0) then
grad = 0.0
else
grad = 2.*grad1*grad2/(grad1 + grad2)
endif
!
! flux
!
if (l==lsal) then
sal = 0.5*(r0(nm, kll, lsal) + r0(nmu, krr, lsal))
if (ltem==0) then
temp = temeqs
else
temp = 0.5*(r0(nm, kll, ltem) + r0(nmu, krr, ltem))
endif
select case(idensform)
case( dens_Eckart )
call dens_eck ( temp, sal, dummy, rhods, dummy )
case( dens_UNESCO )
call dens_unes( temp, sal, dummy, rhods, dummy )
end select
drho = rhods*grad
dsal = grad
dtem = 0.0
elseif (l==ltem) then
temp = 0.5*(r0(nm, kll, ltem) + r0(nmu, krr, ltem))
if (lsal==0) then
sal = saleqs
else
sal = 0.5*(r0(nm, kll, lsal) + r0(nmu, krr, lsal))
endif
select case(idensform)
case( dens_Eckart )
call dens_eck ( temp, sal, dummy, dummy, rhodt )
case( dens_UNESCO )
call dens_unes( temp, sal, dummy, dummy, rhodt )
end select
drho = rhodt*grad
dsal = 0.0
dtem = grad
else
drho = grad*(1.0 - rhow/rhosol(l - lst))
dsal = 0.0
dtem = 0.0
endif
darea = (point(kf - 1) - point(kf))*guu(nm)/gvu(nm)
difco = 0.5*(dicuv(nm, kll) + dicuv(nmu, krr))/0.7
!
! Bottom in U-point dpu(nm) may differ from MIN(dps(nm),dps(nmu))
! formula should be based on dpu
!
zubot = -dpu(nm)
if (zubot