!----- AGPL --------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2015.
!
! This file is part of Delft3D (D-Flow Flexible Mesh component).
!
! Delft3D is free software: you can redistribute it and/or modify
! it under the terms of the GNU Affero General Public License as
! published by the Free Software Foundation version 3.
!
! Delft3D 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 Affero General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License
! along with Delft3D. 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",
! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting
! Deltares, and remain the property of Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: transport.f90 43363 2015-12-01 15:14:53Z kleczek $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/transport.f90 $
!> This subroutine transports an array of scalars.
!> In light of future vectorization, the aim is to:
!> -use as few module variables as possible,
!> -do not use if statements within do-loops.
!> updates sedl with
!> d(h sedl)/dt = -div (q1 sedl) + div (diag(NU) grad sedl) + h ( -sink sedl + source )
!> solves for each column {k | 1<=k aaj(k) sedj(k-1) + bbj(k) sedj(k) + ccj(k) sedj(k+1) = ddj(k)
subroutine update_constituents(jarhoonly)
use m_flowgeom, only: Ndx, Ndxi, Lnxi, Lnx, ln, nd ! static mesh information
use m_flow, only: Ndkx, Lnkx, u1, q1, au, qw, zws, sq, sqi, vol1, kbot, ktop, Lbot, Ltop, kmxn, kmxL, kmx, viu, vicwws, plotlin, jalts
use m_flowtimes, only: dts, ja_timestep_auto
use m_turbulence, only: sigdifi
use m_physcoef, only: dicoww, vicouv, difmolsal
use m_transport
use m_flowparameters, only: limtypsa, limtyptm, limtypsed
use m_alloc
use m_partitioninfo
use m_timer
use unstruc_messages
implicit none
integer :: jarhoonly
integer :: ierror
integer :: limtyp !< limiter type (>0), or first-order upwind (0)
double precision :: dvoli
double precision :: dt, dts_store
integer :: k, L, j, numconst_store
integer :: istep
integer :: numstepssync
if ( NUMCONST.eq.0 ) return ! nothing to do
ierror = 1
limtyp = max(limtypsa, limtyptm, limtypsed)
if (jarhoonly == 1) then
call fill_rho() ; numconst_store = numconst
else
call fill_constituents()
endif
call get_dtmax()
if ( jalts.eq.1 ) then ! local time-stepping
call get_ndeltasteps()
else
nsubsteps = 1
ndeltasteps = 1
numnonglobal = 0
end if
! store dts
dts_store = dts
! set dts to smallest timestep
dts = dts/nsubsteps
if ( jampi.ne.0 ) then
! determine at which sub timesteps to update
if ( limtyp.eq.0 ) then
numstepssync = max(numlay_cellbased-1,1) ! one level used in first-order upwind order reconstruction
else
numstepssync = max(int(numlay_cellbased/2),1) ! two levels used in higher-order upwind reconstruction
end if
end if
jaupdate = 1
fluxhor = 0d0 ! not necessary
sumhorflux = 0d0
do istep=0,nsubsteps-1
if ( kmx.gt.0 ) then
fluxver = 0d0
end if
! BEGIN DEBUG
! difsedu = 0d0
! difsedw = 0d0
! END DEBUG
! BEGIN DEBUG
! call comp_sq(Ndkx, Lnkx, kbot, ktop, Lbot, Ltop, q1, qw, sq)
! END DEBUG
! determine which fluxes need to be updated
if ( nsubsteps.gt.1 ) then
call get_jaupdatehorflux(nsubsteps, limtyp, jaupdate,jaupdatehorflux)
end if
! compute horizontal fluxes, explicit part
call comp_fluxhor3D(NUMCONST, limtyp, Ndkx, Lnkx, u1, q1, au, sqi, vol1, kbot, Lbot, Ltop, kmxn, kmxL, constituents, difsedu, sigdifi, viu, vicouv, nsubsteps, jaupdate, jaupdatehorflux, ndeltasteps, fluxhor)
call starttimer(IDEBUG)
call comp_sumhorflux(NUMCONST, kmx, Lnkx, Ndkx, Lbot, Ltop, fluxhor, sumhorflux)
call stoptimer(IDEBUG)
! determine which cells need to be updated
if ( nsubsteps.gt.1 ) then
call get_jaupdate(istep,nsubsteps,Ndxi,Ndx,ndeltasteps,jaupdate)
end if
if ( kmx.lt.1 ) then ! 2D, call to 3D as well for now
call solve_2D(NUMCONST, Ndkx, Lnkx, vol1, kbot, ktop, Lbot, Ltop, sumhorflux, fluxver, const_sour, const_sink, nsubsteps, jaupdate, ndeltasteps, constituents, rhs)
else
call comp_fluxver( NUMCONST, limtyp, thetavert, Ndkx, kmx, zws, qw, kbot, ktop, constituents, nsubsteps, jaupdate, ndeltasteps, fluxver)
call solve_vertical(NUMCONST, ISED1, ISEDN, limtyp, thetavert, Ndkx, Lnkx, kmx, &
zws, qw, vol1, kbot, ktop, Lbot, Ltop, &
sumhorflux, fluxver, const_sour, const_sink, &
difsedw, sigdifi, vicwws, nsubsteps, jaupdate, ndeltasteps, constituents, &
a, b, c, d, e, sol, rhs)
end if
if ( jampi.gt.0 ) then
! communicate every numstepssync'th or last subtimestep
if ( mod(istep+1,numstepssync).eq.0 .or. istep+1.eq.nsubsteps ) then
if ( jatimer.eq.1 ) call starttimer(IUPDSALL)
if ( kmx.lt.1 ) then ! 2D
call update_ghosts(ITYPE_Sall, NUMCONST, Ndx, constituents, ierror)
else ! 3D
call update_ghosts(ITYPE_Sall3D, NUMCONST, Ndkx, constituents, ierror)
end if
if ( jatimer.eq.1 ) call stoptimer(IUPDSALL)
end if
end if
end do
!! communicate
! if ( jampi.gt.0 ) then
! if ( jatimer.eq.1 ) call starttimer(IUPDSALL)
! if ( kmx.lt.1 ) then ! 2D
! call update_ghosts(ITYPE_Sall, NUMCONST, Ndx, constituents, ierror)
! else ! 3D
! call update_ghosts(ITYPE_Sall3D, NUMCONST, Ndkx, constituents, ierror)
! end if
! if ( jatimer.eq.1 ) call stoptimer(IUPDSALL)
! end if
if (jarhoonly == 1) then
call extract_rho() ; numconst = numconst_store
else
call extract_constituents()
endif
ierror = 0
1234 continue
! restore dts
dts = dts_store
return
end subroutine update_constituents
!> compute horizontal transport fluxes at flowlink
subroutine comp_fluxhor3D(NUMCONST, limtyp, Ndkx, Lnkx, u1, q1, au, sqi, vol1, kbot, Lbot, Ltop, kmxn, kmxL, sed, difsed, sigdifi, viu, vicouv, nsubsteps, jaupdate, jaupdatehorflux, ndeltasteps, flux)
use m_flowgeom, only: Ndx, Lnx, Lnxi, ln, nd, klnup, slnup, dxi, acl ! static mesh information
use m_flowtimes, only: dts, dnt
use m_flowparameters, only: cflmx
use m_flow, only: jadiusp, diusp, dicouv, jacreep, dsalL, dtemL, hu, epshu
use m_transport, only: ISALT, ITEMP
use m_missing
implicit none
integer, intent(in) :: NUMCONST !< number of transported quantities
integer, intent(in) :: limtyp !< limiter type
integer, intent(in) :: Ndkx !< total number of flownodes (dynamically changing)
integer, intent(in) :: Lnkx !< total number of flowlinks (dynamically changing)
double precision, dimension(Lnkx), intent(in) :: u1 !< flow-field face-normal velocities
double precision, dimension(Lnkx), intent(in) :: q1 !< flow-field discharges
double precision, dimension(Lnkx), intent(in) :: au !< wet area of flowlinks, note: q1=au*u1
double precision, dimension(Ndkx), intent(in) :: sqi !< total outward-fluxes at flownodes
double precision, dimension(Ndkx), intent(in) :: vol1 !< volumes
integer, dimension(Ndx), intent(in) :: kbot !< flow-node based layer administration
integer, dimension(Lnx), intent(in) :: Lbot !< flow-link based layer administration
integer, dimension(Lnx), intent(in) :: Ltop !< flow-link based layer administration
integer, dimension(Ndx), intent(in) :: kmxn !< flow-link based layer administration
integer, dimension(Lnx), intent(in) :: kmxL !< flow-link based layer administration
double precision, dimension(NUMCONST,Ndkx), intent(in) :: sed !< transported quantities
double precision, dimension(NUMCONST), intent(in) :: difsed !< scalar-specific diffusion coefficent (dicouv)
real, dimension(Lnkx), intent(in) :: viu !< spatially varying horizontal eddy viscosity, NOTE: real, not double
double precision, intent(in) :: vicouv !< uniform horizontal eddy viscosity
double precision, dimension(NUMCONST), intent(in) :: sigdifi !< 1/(Prandtl number) for heat, 1/(Schmidt number) for mass
integer, intent(in) :: nsubsteps !< number of substeps
integer, dimension(Ndx), intent(in) :: jaupdate !< cell updated (1) or not (0)
integer, dimension(Lnx), intent(in) :: jaupdatehorflux !< update horizontal flux (1) or not (0)
integer, dimension(Ndx), intent(in) :: ndeltasteps !< number of substeps between updates
double precision, dimension(NUMCONST,Lnkx), intent(inout) :: flux !< adds horizontal advection and diffusion fluxes
double precision :: sl1L, sl2L, sl3L, sl1R, sl2R, sl3R
double precision :: cf, sedkuL, sedkuR, ds1L, ds2L, ds1R, ds2R
double precision :: fluxL, fluxR
double precision :: sedL, sedR
double precision :: sumdiffusion
double precision :: fluxfac, fluxfacMaxL, fluxfacMaxR
double precision :: dfac1, dfac2
double precision :: difcoeff, QL, QR, diuspL
double precision :: dt_loc
integer :: j, iswitchL, iswitchR, jahigherL, jahigherR
integer :: k1, k2, LL, L, Lb, Lt, laydif, jaL, jaR
integer :: kk1L, kk2L, kk1R, kk2R, k1L, k2L, k1R, k2R
double precision, external :: dlimiter
continue
dt_loc = dts
!$OMP PARALLEL DO &
!$OMP PRIVATE(LL,L,Lb,Lt,kk1L, kk2L, kk1R, kk2R, k1L, k2L, k1R, k2R, iswitchL, iswitchR, sl1L, sl2L, sl3L, sl1R, sl2R, sl3R) &
!$OMP PRIVATE(cf, k1, k2, laydif, j, sedL, sedR, sedkuL, sedkuR, ds1L, ds2L, ds1R, ds2R, jaL, jaR, QL, QR)
! advection
do LL=1,Lnx
if ( nsubsteps.gt.1 ) then
if ( jaupdatehorflux(LL).eq.0 ) then
cycle
else
dt_loc = dts * min(ndeltasteps(ln(1,LL)),ndeltasteps(ln(2,LL)))
end if
end if
Lb = Lbot(LL)
Lt = Ltop(LL)
! if (limtyp .ne. 0) then
! get the 2D flownodes in the stencil
kk1L = klnup(1,LL)
iswitchL = 1-min(max(kk1L,0),1) ! 1 if kk1L<0, 0 otherwise
kk2L = (1-iswitchL)*klnup(2,LL) + iswitchL*kk1L ! make kk2L safe for when it is not intented to be used
kk1R = klnup(4,LL)
iswitchR = 1-min(max(kk1R,0),1) ! 1 if kk1R<0, 0 otherwise
kk2R = (1-iswitchR)*klnup(5,LL) + iswitchR*kk1R ! make kk2R safe for when it is not intented to be used
! get the weights in the stencil
sl1L = (dble(1-iswitchL)*slnup(1,LL) + dble(iswitchL)*1d0)
sl2L = dble(1-iswitchL)*slnup(2,LL)
sl3L = slnup(3,LL)
sl1R = (dble(1-iswitchR)*slnup(4,LL) + dble(iswitchR)*1d0)
sl2R = dble(1-iswitchR)*slnup(5,LL)
sl3R = slnup(6,LL)
! make cell indices safe
! kk1L = max(iabs(kk1L),1)
! kk2L = max(iabs(kk2L),1)
!
! kk1R = max(iabs(kk1R),1)
! kk2R = max(iabs(kk2R),1)
! make cell indices safe
kk1L = iabs(kk1L)
kk2L = iabs(kk2L)
kk1R = iabs(kk1R)
kk2R = iabs(kk2R)
! endif
! loop over vertical flowlinks
do L=Lb,Lt
! get left and right neighboring flownodes
k1 = ln(1,L)
k2 = ln(2,L)
! if (limtyp > 0) then
! compute Courant number
cf = dt_loc*abs(u1(L))*dxi(LL)
if ( cf.gt.cflmx ) then
! write(6,*) cf
continue
end if
laydif = L-Lb
! reconstuct from the left and from the right
jaL = 0
if ( kk1L.ne.0 ) then
k1L = kbot(kk1L) + laydif + kmxn(kk1L) - kmxL(LL)
if ( kk2L.ne.0 ) then
k2L = kbot(kk2L) + laydif + kmxn(kk2L) - kmxL(LL)
jaL = min(max(k1L-kbot(kk1L)+1,0),1)*min(max(k2L-kbot(kk2L)+1,0),1)*k1L
end if
end if
jaR = 0
if ( kk1R.ne.0 ) then
k1R = kbot(kk1R) + laydif + kmxn(kk1R) - kmxL(LL)
if ( kk2R.ne.0 ) then
k2R = kbot(kk2R) + laydif + kmxn(kk2R) - kmxL(LL)
jaR = min(max(k1R-kbot(kk1R)+1,0),1)*min(max(k2R-kbot(kk2R)+1,0),1)*k1R
end if
end if
! endif
! BEGIN DEBUG
! if ( dnt.eq.5 .and. ( ( k1.eq.1736 .and. q1(L).gt.0d0 ) .or. ( k2.eq.1736 .and. q1(L).lt.0d0 ) ) ) then
! continue
! end if
! END DEBUG
QL = max(q1(L),0d0)
QR = min(q1(L),0d0)
do j=1,NUMCONST
sedL = sed(j,k1)
sedR = sed(j,k2)
! if (limtyp > 0) then
if ( kk1L.ne.0 .and. q1(L).gt.0d0 .and. jaL > 0) then
sedkuL = sed(j,k1L)*sl1L + sed(j,k2L)*sl2L
ds2L = sed(j,k2) - sed(j,k1)
ds1L = (sed(j,k1) - sedkuL)*sl3L
sedL = sedL + acl(LL) *max(0d0,1d0-cf) * dlimiter(ds1L,ds2L,limtyp)
end if
if ( kk1R.ne.0 .and. q1(L).lt.0d0 .and. jaR > 0) then
sedkuR = sed(j,k1R)*sl1R + sed(j,k2R)*sl2R
ds2R = sed(j,k1) - sed(j,k2)
ds1R = (sed(j,k2) - sedkuR)*sl3R
sedR = sedR + (1d0-acl(LL))*max(0d0,1d0-cf) * dlimiter(ds1R,ds2R,limtyp)
end if
! endif
flux(j,L) = QL*sedL + QR*sedR
end do
end do
end do
!$OMP END PARALLEL DO
! BEGIN DEBUG
! return
! END DEBUG
! diffusion
sumdiffusion = sum(difsed(1:NUMCONST))
if ( sumdiffusion.gt.0d0 .and. dicouv >= 0d0) then
!$OMP PARALLEL DO &
!$OMP PRIVATE(LL,dfac1,dfac2,Lb,Lt,L,k1,k2,fluxfacMaxL,fluxfacMaxR,j,difcoeff,fluxfac )
do LL=1,Lnx
if ( nsubsteps.gt.1 ) then
if ( jaupdatehorflux(LL).eq.0 ) then
cycle
else
dt_loc = dts * min(ndeltasteps(ln(1,LL)),ndeltasteps(ln(2,LL)))
end if
end if
!monotinicity criterion, safe for triangles, quad and pentagons, but not for hexahedrons
!dfac1 = 0.2d0
!dfac2 = 0.2d0
dfac1 = 1d0/dble(nd(ln(1,LL))%lnx)
dfac2 = 1d0/dble(nd(ln(2,LL))%lnx)
if (jadiusp == 1) then
diuspl = diusp(LL)
else
diuspL = dicouv
endif
Lb = Lbot(LL)
Lt = Ltop(LL)
do L=Lb,Lt
k1 = ln(1,L)
k2 = ln(2,L)
fluxfacMaxL = dfac1*vol1(k1)/dt_loc - sqi(k1)
fluxfacMaxR = dfac2*vol1(k2)/dt_loc - sqi(k2)
do j=1,NUMCONST
difcoeff = sigdifi(j)*viu(L) + difsed(j) + diuspL ! without smagorinsky, viu is 0 ,
! difsed only contains molecular value,
! so then you only get user specified value
fluxfac = difcoeff*au(L)*dxi(LL)
fluxfac = min(fluxfac, fluxfacMaxL, fluxfacMaxR) ! zie Borsboom sobek note
fluxfac = max(fluxfac, 0d0)
if (jacreep == 0) then
flux(j,L) = flux(j,L) - fluxfac*(sed(j,k2) - sed(j,k1))
else
if (j == ISALT) then
!if (dsalL(L) > 0d0 ) then
! dsalL(L) = max(0d0, min(dsalL(L), sed(j,k2) - sed(j,k1) ) )
!else if (dsalL(L) < 0d0 ) then
! dsalL(L) = min(0d0, max(dsalL(L), sed(j,k2) - sed(j,k1) ) )
!endif
flux(j,L) = flux(j,L) - fluxfac*dsalL(L)
else if (j == Itemp) then
!if (dtemL(L) > 0 ) then
! dtemL(L) = max(0d0, min(dtemL(L), sed(j,k2) - sed(j,k1) ) )
!else if (dtemL(L) < 0 ) then
! dtemL(L) = min(0d0, max(dtemL(L), sed(j,k2) - sed(j,k1) ) )
!endif
flux(j,L) = flux(j,L) - fluxfac*dtemL(L)
endif
endif
end do
end do
end do
!$OMP END PARALLEL DO
end if
return
end subroutine comp_fluxhor3D
!> compute vertical fluxes
subroutine comp_fluxver(NUMCONST, limtyp, thetavert, Ndkx, kmx, zws, qw, kbot, ktop, sed, nsubsteps, jaupdate, ndeltasteps, flux)
use m_flowgeom, only: Ndx, ba, kfs ! static mesh information
use m_flowtimes, only: dts
use m_flowparameters, only: cflmx
use m_flow, only : hs, epshs, s1, epshsdif, cffacver ! do not use m_flow, please put this in the argument list
implicit none
integer, intent(in) :: NUMCONST !< number of transported quantities
integer, intent(in) :: limtyp !< limiter type
double precision, dimension(NUMCONST), intent(in) :: thetavert !< compute fluxes (<1) or not (1)
integer, intent(in) :: Ndkx !< total number of flownodes (dynamically changing)
integer, intent(in) :: kmx !< maximum number of layers (dynamically changing)
double precision, dimension(Ndkx), intent(in) :: zws !< vertical coordinate of layers at interface/center locations
double precision, dimension(Ndkx), intent(in) :: qw !< flow-field vertical discharges
integer, dimension(Ndx), intent(in) :: kbot !< flow-node based layer administration
integer, dimension(Ndx), intent(in) :: ktop !< flow-node based layer administration
double precision, dimension(NUMCONST,Ndkx), intent(in) :: sed !< transported quantities
integer, intent(in) :: nsubsteps !< number of substeps
integer, dimension(Ndx), intent(in) :: jaupdate !< update cell (1) or not (0)
integer, dimension(Ndx), intent(in) :: ndeltasteps !< number of substeps between updates
double precision, dimension(NUMCONST,Ndkx), intent(inout) :: flux !< adds vertical advection fluxes
double precision, dimension(2049) :: dz
double precision :: sedL, sedR, ds1L, ds2L, ds1R, ds2R, sl3L, sl3R, cf, dum
double precision :: dt_loc
integer :: kk, k, kb, kt, kL, kR, kLL, kRR
integer :: j
double precision, external :: dlimiter
double precision, parameter :: DTOL = 1d-8
if ( sum(1d0-thetavert(1:NUMCONST)).lt.DTOL ) goto 1234 ! nothing to do
dt_loc = dts
!$OMP PARALLEL DO &
!$OMP PRIVATE(kk,kb,kt,dz,k,cf,kL,kR,j,sedL,sedR,kLL,kRR,sl3L,sl3R,ds1L,ds1R,ds2L,ds2R )
do kk=1,Ndx
if (kfs(kk) == 0) cycle
if ( nsubsteps.gt.1 ) then
if (jaupdate(kk).eq.0 ) then
cycle
else
dt_loc = dts * ndeltasteps(kk)
end if
end if
kb = kbot(kk)
kt = ktop(kk)
dz(1) = max(dtol, zws(kb)-zws(kb-1) )
! dz(2:kt-kb+1) = max(dtol, 0.5d0*(zws(kb+1:kt)-zws(kb-1:kt-1)) ) ! thickness between cell centers org
dz(2:kt-kb+1) = max(dtol, 0.5d0*(zws(kb+1:kt)-zws(kb-1:kt-2)) ) ! thickness between cell centers
dz(kt-kb+2) = max(dtol, zws(kt)-zws(kt-1) )
! dz = max(dz,dtol) ! fix for zero-thickness layer
do k=kb,kt-1
! cf = dt_loc*qw(k)/(ba(kk)*dz(k-kb+2))
cf = cffacver*dt_loc*abs(qw(k))/(ba(kk)*dz(k-kb+2))
if ( cf.gt.cflmx ) then
continue
end if
kL = max(k,kb)
kR = min(k+1,kt)
do j=1,NUMCONST
if ( thetavert(j).eq.1d0 ) cycle
sedL = sed(j,kL)
sedR = sed(j,kR)
if ( thetavert(j).gt.0d0 ) then ! semi-explicit, use central scheme
flux(j,k) = flux(j,k) + qw(k)*0.5d0*(sedL+sedR)
else ! fully explicit
! if (limtyp.ne.0 ) then
if ( k.gt.kb-1 .and. qw(k).gt.0d0 ) then
kLL = max(k-1,kb)
sL3L = dz(k-kb+2)/dz(k-kb+1)
! if ( abs(sL3L-1d0).gt.1d-4 ) then
! continue
! end if
ds2L = sedR-sedL
ds1L = (sedL-sed(j,kLL))*sl3L
sedL = sedL + 0.5d0 * max(0d0,1d0-cf) * dlimiter(ds1L,ds2L,limtyp)
end if
if ( k.lt.kt .and. qw(k).lt.0d0 .and. s1(kk)-zws(kb-1) > epshsdif ) then
kRR = min(k+2,kt)
sL3R = dz(k-kb+2)/dz(k-kb+3)
! if ( abs(sL3R-1d0).gt.1d-4 ) then
! continue
! end if
ds2R = sedL-sedR
ds1R = (sedR-sed(j,kRR))*sl3R
sedR = sedR + 0.5d0 * max(0d0,1d0-cf) * dlimiter(ds1R,ds2R,limtyp)
end if
! end if
flux(j,k) = flux(j,k) + max(qw(k),0d0)*sedL + min(qw(k),0d0)*sedR
end if
end do
end do
end do
!$OMP END PARALLEL DO
1234 continue
return
end subroutine comp_fluxver
!> compose right-hand side
subroutine make_rhs(NUMCONST, thetavert, Ndkx, Lnkx, kmx, vol1, kbot, ktop, Lbot, Ltop, sumhorflux, fluxver, source, sed, nsubsteps, jaupdate, ndeltasteps, rhs)
use m_flowgeom, only: Ndxi, Ndx, Lnx, Ln, ba ! static mesh information
use m_flowtimes, only: dts
implicit none
integer, intent(in) :: NUMCONST !< number of transported quantities
double precision, dimension(NUMCONST), intent(in) :: thetavert !< vertical advection explicit (0) or implicit (1)
integer, intent(in) :: Ndkx !< total number of flownodes (dynamically changing)
integer, intent(in) :: Lnkx !< total number of flowlinks (dynamically changing)
integer, intent(in) :: kmx !< maximum number of layers
! double precision, dimension(Ndkx), intent(in) :: sq !< flux balance (inward positive)
double precision, dimension(Ndkx), intent(in) :: vol1 !< volumes
integer, dimension(Ndx), intent(in) :: kbot !< flow-node based layer administration
integer, dimension(Ndx), intent(in) :: ktop !< flow-node based layer administration
integer, dimension(Lnx), intent(in) :: Lbot !< flow-link based layer administration
integer, dimension(Lnx), intent(in) :: Ltop !< flow-link based layer administration
double precision, dimension(NUMCONST,Ndkx), intent(inout) :: sumhorflux !< sum of horizontal fluxes
double precision, dimension(NUMCONST,Ndkx), intent(in) :: fluxver !< vertical fluxes
double precision, dimension(NUMCONST,Ndkx), intent(in) :: source !< sources
double precision, dimension(NUMCONST,Ndkx), intent(inout) :: sed !< transported quantities
integer, intent(in) :: nsubsteps !< number of substeps
integer, dimension(Ndx), intent(in) :: jaupdate !< update cell (1) or not (0)
integer, dimension(Ndx), intent(in) :: ndeltasteps !< number of substeps between updates
double precision, dimension(NUMCONST,Ndkx), intent(out) :: rhs ! right-hand side, dim(NUMCONST,Ndkx)
double precision :: dvoli
double precision :: dt_loc
integer :: LL, L, Lb, Lt
integer :: kk, k, kb, kt
integer :: k1, k2, j
double precision, parameter :: dtol=1d-8
dt_loc = dts
! rhs = 0d0
!
!! add horizontal fluxes to right-hand side
! do LL=1,Lnx
! Lb = Lbot(LL)
! Lt = Ltop(LL)
! do L=Lb,Lt
!! get neighboring flownodes
! k1 = ln(1,L)
! k2 = ln(2,L)
! do j=1,NUMCONST
! rhs(j,k1) = rhs(j,k1) - fluxhor(j,L)
! rhs(j,k2) = rhs(j,k2) + fluxhor(j,L)
! end do
! end do
! end do
if ( kmx.gt.0 ) then
! add vertical fluxes, sources, storage term and time derivative to right-hand side
!$OMP PARALLEL DO &
!$OMP PRIVATE(kk,kb,kt,k,dvoli,j )
do kk=1,Ndxi
if ( jaupdate(kk).eq.0 ) then
cycle
else
dt_loc = dts * ndeltasteps(kk)
end if
kb = kbot(kk)
kt = ktop(kk)
do k=kb,kt
dvoli = 1d0/max(vol1(k),dtol)
do j=1,NUMCONST
! rhs(j,k) = ((rhs(j,k) - (1d0-thetavert(j))*(fluxver(j,k) - fluxver(j,k-1)) - sed(j,k)*sq(k)) * dvoli + source(j,k))*dts + sed(j,k)
rhs(j,k) = ((sumhorflux(j,k)/ndeltasteps(kk) - (1d0-thetavert(j))*(fluxver(j,k) - fluxver(j,k-1))) * dvoli + source(j,k))*dt_loc + sed(j,k)
sumhorflux(j,k) = 0d0
! BEGIN DEBUG
! rhs(j,k) = source(j,k)*dts + sed(j,k)
! END DEBUG
end do
end do
end do
!$OMP END PARALLEL DO
else
! add time derivative
if ( nsubsteps.eq.1 ) then
!$OMP PARALLEL DO &
!$OMP PRIVATE(k,j,dvoli )
do k=1,Ndxi
dvoli = 1d0/max(vol1(k),dtol)
do j=1,NUMCONST
rhs(j,k) = (sumhorflux(j,k) * dvoli + source(j,k)) * dts + sed(j,k)
sumhorflux(j,k) = 0d0
end do
end do
!$OMP END PARALLEL DO
else
!$OMP PARALLEL DO &
!$OMP PRIVATE(k,j,dvoli )
do k=1,Ndxi
if ( jaupdate(k).eq.0 ) then
cycle
else
dt_loc = dts * ndeltasteps(k)
end if
dvoli = 1d0/max(vol1(k),dtol)
do j=1,NUMCONST
rhs(j,k) = (sumhorflux(j,k)/ndeltasteps(k) * dvoli + source(j,k)) * dt_loc + sed(j,k)
sumhorflux(j,k) = 0d0
end do
end do
!$OMP END PARALLEL DO
end if
end if
end subroutine make_rhs
!> compose right-hand side
subroutine solve_2D(NUMCONST, Ndkx, Lnkx, vol1, kbot, ktop, Lbot, Ltop, sumhorflux, fluxver, source, sink, nsubsteps, jaupdate, ndeltasteps, sed, rhs)
use m_flowgeom, only: Ndxi, Ndx, Lnx, Ln, ba ! static mesh information
use m_flowtimes, only: dts
implicit none
integer, intent(in) :: NUMCONST !< number of transported quantities
integer, intent(in) :: Ndkx !< total number of flownodes (dynamically changing)
integer, intent(in) :: Lnkx !< total number of flowlinks (dynamically changing)
! double precision, dimension(Ndkx), intent(in) :: sq !< flux balance (inward positive)
double precision, dimension(Ndkx), intent(in) :: vol1 !< volumes
integer, dimension(Ndkx), intent(in) :: kbot !< flow-node based layer administration
integer, dimension(Ndkx), intent(in) :: ktop !< flow-node based layer administration
integer, dimension(Lnkx), intent(in) :: Lbot !< flow-link based layer administration
integer, dimension(Lnkx), intent(in) :: Ltop !< flow-link based layer administration
double precision, dimension(NUMCONST,Ndkx), intent(in) :: sumhorflux !< sum of horizontal fluxes
double precision, dimension(NUMCONST,Ndkx), intent(in) :: fluxver !< vertical fluxes
double precision, dimension(NUMCONST,Ndkx), intent(in) :: source !< sources
double precision, dimension(NUMCONST,Ndkx), intent(in) :: sink !< linearized sinks
integer, intent(in) :: nsubsteps !< number of substeps
integer, dimension(Ndx), intent(in) :: jaupdate !< update cell (1) or not (0)
integer, dimension(Ndx), intent(in) :: ndeltasteps !< number of substeps between updates
double precision, dimension(NUMCONST,Ndkx), intent(inout) :: sed !< transported quantities
double precision, dimension(NUMCONST,Ndkx) :: rhs ! work array: right-hand side, dim(NUMCONST,Ndkx)
double precision, dimension(NUMCONST) :: thetavert
double precision :: dt_loc
integer :: j, k
thetavert = 0d0
dt_loc = dts
call make_rhs(NUMCONST, thetavert, Ndkx, Lnkx, 0, vol1, kbot, ktop, Lbot, Ltop, sumhorflux, fluxver, source, sed, nsubsteps, jaupdate, ndeltasteps, rhs)
!$OMP PARALLEL DO &
!$OMP PRIVATE(k,j)
do k=1,Ndxi
if ( nsubsteps.gt.1 ) then
if ( jaupdate(k).eq.0 ) then
cycle
else
dt_loc = dts*ndeltasteps(k)
end if
end if
do j=1,NUMCONST
sed(j,k) = rhs(j,k) / (1d0 + dt_loc*sink(j,k))
end do
end do
!$OMP END PARALLEL DO
return
end subroutine solve_2D
!> solve equations implicitly in vertical direction
subroutine solve_vertical(NUMCONST, ISED1, ISEDN, limtyp, thetavert, Ndkx, Lnkx, kmx, &
zws, qw, vol1, kbot, ktop, Lbot, Ltop, &
sumhorflux, fluxver, source, sink, &
difsed, sigdifi, vicwws, &
nsubsteps, jaupdate, ndeltasteps, sed, &
a, b, c, d, e, sol, rhs)
use m_flowgeom, only: Ndxi, Ndx, Lnx, Ln, ba, kfs, bl ! static mesh information
use m_flowtimes, only: dts
use m_flow, only: epshsdif, s1, kmxn ! do not use m_flow, please put this in the argument list
use m_sediment, only: mtd
implicit none
integer, intent(in) :: NUMCONST !< number of transported quantities
integer, intent(in) :: ISED1 !< index of first sediment fraction in constituents array
integer, intent(in) :: ISEDN !< index of last sediment fraction in constituents array
integer, intent(in) :: limtyp !< limiter type
double precision, dimension(NUMCONST), intent(in) :: thetavert !< vertical advection explicit (0) or implicit (1)
integer, intent(in) :: Ndkx !< total number of flownodes (dynamically changing)
integer, intent(in) :: Lnkx !< total number of flowlinks (dynamically changing)
integer, intent(in) :: kmx !< maximum number of layers
double precision, dimension(Ndkx), intent(in) :: zws !< vertical coordinate of layers at interface/center locations
double precision, dimension(Ndkx), intent(in) :: qw !< flow-field vertical discharges
! double precision, dimension(Ndkx), intent(in) :: sq !< flux balance (inward positive)
double precision, dimension(Ndkx), intent(in) :: vol1 !< volumes
integer, dimension(Ndx), intent(in) :: kbot !< flow-node based layer administration
integer, dimension(Ndx), intent(in) :: ktop !< flow-node based layer administration
integer, dimension(Lnx), intent(in) :: Lbot !< flow-link based layer administration
integer, dimension(Lnx), intent(in) :: Ltop !< flow-link based layer administration
double precision, dimension(NUMCONST,Ndkx), intent(in) :: sumhorflux !< sum of horizontal fluxes
double precision, dimension(NUMCONST,Ndkx), intent(in) :: fluxver !< vertical fluxes
double precision, dimension(NUMCONST,Ndkx), intent(in) :: source !< sources
double precision, dimension(NUMCONST,Ndkx), intent(in) :: sink !< sinks
double precision, dimension(NUMCONST), intent(in) :: difsed !< scalar-specific diffusion coefficent (dicoww+difmod)
double precision, dimension(Ndkx), intent(in) :: vicwws !< vertical eddy viscosity, NOTE: real, not double
double precision, dimension(NUMCONST), intent(in) :: sigdifi !< 1/(Prandtl number) for heat, 1/(Schmidt number) for mass
integer, intent(in) :: nsubsteps !< number of substeps
integer, dimension(Ndx), intent(in) :: jaupdate !< update cell (1) or not (0)
integer, dimension(Ndx), intent(in) :: ndeltasteps !< number of substeps
double precision, dimension(NUMCONST,Ndkx), intent(inout) :: sed !< transported quantities
double precision, dimension(kmx,NUMCONST) :: a,b,c,d ! work array: aj(i,j)*sed(j,k-1) + bj(i,j)*sed(j,k) + c(i,j)*sed(j,k+1) = d(i), i=k-kb+1
double precision, dimension(kmx) :: sol, e ! work array: solution and dummy array in tridag, respectively
double precision, dimension(NUMCONST,Ndkx) :: rhs ! work array: right-hand side, dim(NUMCONST,Ndkx)
double precision :: dvoli, fluxfac, dvol1i, dvol2i
double precision :: dum, rhstot, dtbazi
integer :: LL, L, Lb, Lt
integer :: kk, k, kb, kt, ktx
integer :: k1, k2, i, j, n
double precision :: dt_loc
double precision, parameter :: dtol=1d-8
dt_loc = dts
rhs = 0d0
call make_rhs(NUMCONST, thetavert, Ndkx, Lnkx, kmx, vol1, kbot, ktop, Lbot, Ltop, sumhorflux, fluxver, source, sed, nsubsteps, jaupdate, ndeltasteps, rhs)
! construct and solve system
!$OMP PARALLEL DO &
!$OMP PRIVATE(kk,kb,ktx,kt,a,b,c,sol,j,d,k,n,dvol1i,dvol2i,fluxfac,e,dtbazi)
do kk=1,Ndxi
if ( nsubsteps.gt.1 ) then
if ( jaupdate(kk).eq.0 ) then
cycle
else
dt_loc = dts * ndeltasteps(kk)
end if
end if
kb = kbot(kk)
kt = ktop(kk)
if (kfs(kk) == 0) cycle
ktx = kb + kmxn(kk) - 1
a = 0d0
c = 0d0
sol = 0d0
! add linearized sinks to diagonal
do k=kb,kt
n = k-kb+1 ! layer number
do j=1,NUMCONST
b(n,j) = 1d0 + dt_loc*sink(j,k)
end do
end do
do j=1,NUMCONST
d(1,j) = rhs(j,kb)
do k=kb,kt ! assume zero-fluxes at boundary and top
n = k-kb+1 ! layer number
d(n,j) = rhs(j,k)
end do
end do
! if ( s1(kk)-bl(kk) > epshsdif ) then
if ( s1(kk)-zws(kb-1) > epshsdif ) then
do k=kb,kt-1 ! assume zero-fluxes at boundary and top
n = k-kb+1 ! layer number
dvol1i = 1d0/max(vol1(k),dtol) ! dtol: safety
dvol2i = 1d0/max(vol1(k+1),dtol) ! dtol: safety
dtbazi = dt_loc*ba(kk)/ ( 0.5d0*(zws(k+1)-zws(k-1)) )
do j=1,NUMCONST
! ! diffusion
if ( j.lt.ISED1 .or. j.gt.ISEDN ) then ! non-sediment
fluxfac = (sigdifi(j)*vicwws(k) + difsed(j))*dtbazi
else
fluxfac = (mtd%seddif(j-ISED1+1,k))*dtbazi
end if
! BEGIN DEBUG
! fluxfac = dt_loc * (difsed(j)) *ba(kk) / ( 0.5d0*(zws(k+1) - zws(k-1)) ) ! m3
! END DEBUG
b(n,j) = b(n,j) + fluxfac*dvol1i
c(n,j) = c(n,j) - fluxfac*dvol1i
b(n+1,j) = b(n+1,j) + fluxfac*dvol2i
a(n+1,j) = a(n+1,j) - fluxfac*dvol2i
! advection
if ( thetavert(j).gt.0d0 ) then ! semi-implicit, use central scheme
! BEGIN DEBUG
! if ( .false. .and. thetavert(j).gt.0d0 ) then ! semi-implicit, use central scheme
! END DEBUG
fluxfac = qw(k)*0.5d0*thetavert(j)*dt_loc
a(n+1,j) = a(n+1,j) - fluxfac*dvol2i
b(n+1,j) = b(n+1,j) - fluxfac*dvol2i
b(n,j) = b(n,j) + fluxfac*dvol1i
c(n,j) = c(n,j) + fluxfac*dvol1i
end if
end do
end do
endif
! solve system(s)
do j=1,NUMCONST
! if ( kk.eq.2 .and. j.eq.1 ) then
! do k=kb,kt
! n = k-kb+1
! write(6,*) n, a(n,j), b(n,j), c(n,j), d(n,j)
! end do
! end if
call tridag(a(1,j), b(1,j), c(1,j), d(1,j), e, sol, kt-kb+1)
sed(j,kb:kt) = sol(1:kt-kb+1)
sed(j,kt+1:ktx) = sed(j,kt)
! BEGIN DEBUG
! do k=kb,kt
! if ( j.eq.1 .and. ( sed(j,k).gt.30.0001 .or. sed(j,k).lt.-0.0001 ) ) then
! continue
! write(6,*) 'kk=', kk, 'lay=', k-kb+1
! write(6,*) 'rhs=', rhs(j,k)
! write(6,*) 'sed=', sed(j,k)
! call qnerror(' ', ' ', ' ')
! end if
! end do
! END DEBUG
end do
end do
!$OMP END PARALLEL DO
return
end subroutine solve_vertical
!> limiter function
double precision function dlimiter(d1,d2,limtyp)
implicit none
double precision, intent(in) :: d1, d2 !< left and right slopes
integer , intent(in) :: limtyp !< first order upwind (0) or MC (>0)
double precision :: r
double precision, parameter :: dtol=1d-16
double precision, parameter :: TWO=2.0d0
dlimiter = 0d0
if (limtyp == 0) return
if ( d1*d2.lt.dtol ) return
r = d1/d2 ! d1/d2
! Monotinized Central
dlimiter = dble(min(limtyp,1)) * max(0d0, min(TWO*r,TWO,0.5d0*(1d0+r)) ) * d2
end function dlimiter
double precision function dlimitercentral(d1,d2,limtyp) ! as dlimiter, now for central gradient instead of slope
implicit none
double precision, intent(in) :: d1, d2 !< left and right slopes
integer , intent(in) :: limtyp !< first order upwind (0) or MC (>0)
double precision :: r
double precision, parameter :: dtol=1d-16
dlimitercentral = 0d0
if (limtyp == 0) return
if ( d1*d2.lt.dtol ) return
r = d1/d2 ! d1/d2
r = 2d0*r - 1d0
dlimitercentral = d2 * max(0d0, min(2d0*r,0.5d0*(1d0+r),2d0) ) ! Monotonized Central
end function dlimitercentral
!> initialize transport, set the enumerators
subroutine ini_transport()
use m_transport
use m_flowparameters
use m_sediment
use m_physcoef
implicit none
character(len=8) :: str
integer :: i, itrace, ised
NUMCONST = 0
ISALT = 0
ITEMP = 0
ISED1 = 0
ISEDN = 0
ISPIR = 0
ITRA1 = 0
ITRAN = 0
if ( jasal.ne.0 ) then
NUMCONST = NUMCONST+1
ISALT = NUMCONST
end if
if ( jatem.ne.0 ) then
NUMCONST = NUMCONST+1
ITEMP = NUMCONST
end if
if ( jased.ne.0 ) then
NUMCONST = NUMCONST+1
ISED1 = NUMCONST
NUMCONST = NUMCONST+mxgr-1
ISEDN = NUMCONST
end if
if ( jasecf > 0 ) then
NUMCONST = NUMCONST+1
ISPIR = NUMCONST
endif
call alloc_transport()
if ( ISALT.gt.0 ) then
if ( javasal.ne.6 ) then
thetavert(ISALT) = tetav ! Central implicit
else
thetavert(ISALT) = 0d0 ! Ho explicit
end if
const_names(ISALT) = 'salt'
end if
if ( ITEMP.gt.0 ) then
if ( javatem.ne.6 ) then
thetavert(ITEMP) = tetav ! 0.55d0
else
thetavert(ITEMP) = 0d0
end if
const_names(ITEMP) = 'temperature'
end if
if ( ISED1.gt.0 ) then
if ( javased.ne.6 ) then
thetavert(ISED1:ISEDN) = tetav
else
thetavert(ISED1:ISEDN) = 0d0
end if
do i=ISED1,ISEDN
ised = i-ISED1+1
write(str,"(I0)") ised
const_names(i) = 'sediment_'//trim(str)
end do
end if
if ( jasecf > 0 ) then
const_names(ISPIR) = 'secondary_flow_intensity'
end if
if ( ITRA1.gt.0 ) then
do i=ITRA1,ITRAN
itrace = i-ITRA1+1
write(str,"(I0)") itrace
const_names(i) = 'tracer_'//trim(str)
end do
end if
! iconst_cur = min(NUMCONST,ITRA1)
iconst_cur = min(NUMCONST,1)
! local timestepping
time_dtmax = -1d0 ! cfl-numbers not evaluated
nsubsteps = 1
ndeltasteps = 1
jaupdatehorflux = 1
numnonglobal = 0
return
end subroutine ini_transport
!> allocate transport arrays
subroutine alloc_transport()
use m_flowgeom, only: Ndx, Lnx
use m_flow, only: Lnkx, Ndkx, kmx, sigdifi
use m_transport
use m_alloc
use m_meteo, only: numtracers
implicit none
! allocate and initialize fluxes
call realloc(fluxhor, (/ NUMCONST, Lnkx /), keepExisting=.true., fill=0d0)
call realloc(fluxver, (/ NUMCONST, Ndkx /), keepExisting=.true., fill=0d0)
call realloc(difsedu, NUMCONST, keepExisting=.true., fill=0d0)
call realloc(difsedw, NUMCONST, keepExisting=.true., fill=0d0)
call realloc(sigdifi, NUMCONST, keepExisting=.true., fill=0d0)
call realloc(constituents, (/ NUMCONST, Ndkx /), keepExisting=.true., fill=0d0)
call realloc(const_sour, (/ NUMCONST, Ndkx /), keepExisting=.true., fill=0d0)
call realloc(const_sink, (/ NUMCONST, Ndkx /), keepExisting=.true., fill=0d0)
call realloc(thetavert, NUMCONST, keepExisting=.true., fill=0d0)
call realloc(const_names, NUMCONST, keepExisting=.true., fill='')
call realloc(id_const, (/ 2, NUMCONST /), keepExisting=.true., fill = 0)
call realloc(sumhorflux, (/ NUMCONST, Ndkx /), keepExisting=.false., fill=0d0)
call realloc(ndeltasteps, Ndx, keepExisting=.false., fill=1)
call realloc(jaupdate, Ndx, keepExisting=.false., fill=1)
call realloc(jaupdatehorflux, Lnx, keepExisting=.false., fill=1)
call realloc(dtmax, Ndx, keepExisting=.false., fill=0d0)
! work arrays
if ( allocated(rhs) ) deallocate(rhs)
allocate(rhs(NUMCONST,Ndkx))
if ( kmx.gt.0 ) then ! 3D
if ( allocated(a) ) deallocate(a,b,c,d,e,sol)
allocate(a(kmx,NUMCONST),b(kmx,NUMCONST),c(kmx,NUMCONST),d(kmx,NUMCONST),e(kmx),sol(kmx))
a = 0d0
b = 0d0
c = 0d0
d = 0d0
e = 0d0
end if
! tracer boundary condition
call realloc(itrac2const, numtracers, keepExisting=.true., fill=0)
return
end subroutine alloc_transport
!> fill constituent array
subroutine fill_constituents()
use m_transport
use m_flowgeom
use m_flow
use m_sediment
use m_transport
use m_sferic
use m_flowtimes , only : dnt
implicit none
double precision :: dvoli
integer :: i, iconst, j, kk, k, kb, kt
double precision, parameter :: dtol=1d-8
double precision :: spir_ce, spir_be, spir_e, alength_a, time_a, alpha, fcoriocof
const_sour = 0d0
const_sink = 0d0
do k=1,Ndkx
if ( ISALT.ne.0 ) then
constituents(ISALT,k) = sa1(k)
end if
if ( ITEMP.ne.0 ) then
constituents(ITEMP,k) = tem1(k)
end if
if( JASECF > 0 ) then
constituents(ISPIR,k) = spirint(k)
endif
if ( ISED1.ne.0 ) then
do i=1,mxgr
iconst = ISED1+i-1
constituents(iconst,k) = sed(i,k)
end do
end if
end do
difsedu = 0d0 ; difsedw = 0d0 ; sigdifi = 0d0
! diffusion coefficients
if ( ISALT.ne.0 ) then
if (dicouv .ge. 0d0) then
difsedu(ISALT) = difmolsal
endif
if (dicoww .ge. 0d0) then
difsedw(ISALT) = dicoww + difmolsal
sigdifi(ISALT) = 1d0/sigsal
endif
end if
if ( ITEMP.ne.0 ) then
if (dicouv .ge. 0d0) then
difsedu(ITEMP) = difmoltem
endif
if (dicoww .ge. 0d0) then
difsedw(ITEMP) = dicoww + difmoltem
sigdifi(ITEMP) = 1d0/sigtem
endif
end if
if( jasecf > 0 ) then
difsedu(ISPIR) = 0d0
difsedw(ISPIR) = 0d0 !dicoww + difmoltem
sigdifi(ISPIR) = 0d0 !/sigspi
endif
if ( ISED1.ne.0 ) then
do i=1,mxgr
iconst = ISED1+i-1
if (dicouv .ge. 0d0) difsedu(iconst) = 0d0
if (dicoww .ge. 0d0) difsedw(iconst) = dicoww
sigdifi(iconst) = 1d0/sigsed
end do
end if
if ( ITRA1.gt.0 ) then
do i=ITRA1,ITRAN
difsedu(i) = difmoltr
difsedw(i) = dicoww + difmoltr
sigdifi(i) = 1d0
end do
end if
! sources
do kk=1,Ndx
call getkbotktop(kk,kb,kt)
do k=kb,kt
dvoli = 1d0/max(vol1(k),dtol)
! salt
if ( ISALT.gt.0 .and. jasalsrc.gt.0 ) then
const_sour(ISALT,k) = salsrc(k)*dvoli
end if
! temperature
if ( ITEMP.gt.0 ) then
const_sour(ITEMP,k) = heatsrc(k)*dvoli
end if
! terms due to non-conservative formulation
do j=1,NUMCONST
const_sour(j,k) = const_sour(j,k) - constituents(j,k) * sq(k) * dvoli
end do
end do
! Note: from now on, only _add_ to sources
! spiral flow source term
if ( jasecf > 0 .and. kmx.eq.0 ) then
if( spirucm(kk) < 1d-3 .or. hs(kk) < epshu ) then
const_sour(ISPIR,kk) = 0d0
const_sink(ISPIR,kk) = 0d0
else
fcoriocof = fcorio ; if( icorio > 0 .and. jsferic == 1 ) fcoriocof = fcoris(kk)
alpha = sqrt( ag ) / vonkar / czssf(kk)
spir_ce = fcorio * hs(kk) * 0.5d0
spir_be = hs(kk) * spircrv(kk) * spirucm(kk)
spir_e = spir_be - spir_ce
alength_a = ( 1.0d0 - 2.0d0 * alpha ) * hs(kk) / ( 2.0d0 * vonkar * vonkar * alpha )
time_a = alength_a / spirucm(kk)
!const_sour(ISPIR,kk) = ( spir_e - spirint(kk) ) / time_a !* dvoli ! S=(I_eq - I)/Ta
const_sour(ISPIR,kk) = spir_e / time_a
const_sink(ISPIR,kk) = 1d0 / time_a
endif
end if
! sediment (2D sources, also in 3D)
if ( stm_included ) then ! 2D only
if ( ISED1.gt.0 .and. kk.le.Ndxi ) then
do i=1,mxgr
iconst = i+ISED1-1
const_sour(iconst,kb) = const_sour(iconst,kb) + sedtra%sourse(kk,i)
const_sink(iconst,kb) = const_sink(iconst,kb) + sedtra%sinkse(kk,i)
end do
end if
end if
end do
return
end subroutine fill_constituents
subroutine fill_rho()
use m_transport
use m_flowgeom
use m_flow
use m_sediment
use m_transport
use m_sferic
use m_flowtimes , only : dnt
implicit none
integer :: kk, k, kb, kt
double precision :: dvoli, dtol=1d-8
do k=1,Ndkx
constituents(1,k) = rho(k)
enddo
! sources
do kk=1,Ndx
call getkbotktop(kk,kb,kt)
do k=kb,kt
dvoli = 1d0/max(vol1(k),dtol)
const_sour(1,k) = - rho(k) * sq(k) * dvoli
const_sink(1,k) = 0d0
end do
enddo
return
end subroutine fill_rho
!> extract constituent array
subroutine extract_constituents()
use m_transport
use m_flow
use m_sediment
use m_transport
implicit none
integer :: i, iconst, k
do k=1,Ndkx
if ( ISALT.ne.0 ) then
sa1(k) = constituents(ISALT,k)
end if
if ( ITEMP.ne.0 ) then
tem1(k) = constituents(ITEMP,k)
end if
if( jasecf > 0 ) then
spirint(k) = constituents(ISPIR,k)
endif
if ( ISED1.ne.0 ) then
do i=1,mxgr
iconst = ISED1+i-1
sed(i,k) = constituents(iconst,k)
end do
end if
end do
return
end subroutine extract_constituents
subroutine extract_rho()
use m_transport
use m_flow
use m_sediment
use m_transport
implicit none
integer :: k
do k=1,Ndkx
rho(k) = constituents(1,k)
enddo
return
end subroutine extract_rho
!> add tracer to constituents, or get constituents number if tracer already exists
subroutine add_tracer(tracer_name, iconst)
use m_transport
use unstruc_messages
use m_meteo, only: numtracers, trnames
implicit none
character(len=*), intent(in) :: tracer_name ! tracer name, or '' for default name
integer, intent(out) :: iconst ! constituent number
character(len=8) :: str
integer :: ierror, i, itrac
integer, external :: findname
ierror = 1
! check if tracer already exist, based on tracer name
iconst = 0
if ( ITRA1.gt.0 .and. trim(tracer_name).ne.'') then
iconst = findname(NUMCONST, const_names, tracer_name)
if ( iconst.ge.ITRA1 .and. iconst.le.ITRAN ) then ! existing tracer found
call mess(LEVEL_INFO, 'add_tracer: tracer ' // trim(tracer_name) // ' already exists. No update required.')
goto 1234
end if
end if
! apend tracer
if ( ITRA1.eq.0) then
NUMCONST = NUMCONST+1
ITRA1 = NUMCONST
ITRAN = NUMCONST
else
! check if tracers are at the back
if ( iTRAN.ne.NUMCONST ) then
call mess(LEVEL_ERROR, 'add_tracer: tracer(s) not at the back of the constituents array')
goto 1234
end if
NUMCONST = NUMCONST+1
ITRAN = NUMCONST
end if
! output tracer number
iconst = ITRAN
! reallocate arrays
call alloc_transport()
! set name
if ( trim(tracer_name).ne.'' ) then
const_names(ITRAN) = trim(tracer_name)
else
write(str,"(I0)") ITRAN-ITRA1+1
const_names(ITRAN) = 'tracer_'//trim(str)
end if
if ( numtracers.gt.0 ) then ! number of tracers
! generate tracer to constituent
itrac = findname(numtracers,trnames,trim(tracer_name))
if ( itrac.gt.0 ) then
itrac2const(itrac) = iconst
end if
end if
ierror = 0
1234 continue
return
end subroutine add_tracer
subroutine comp_sq(Ndkx, Lnkx, kbot, ktop, Lbot, Ltop, q1, qw, sq)
use m_flowgeom, only: Ndx, Lnx, ln
implicit none
integer, intent(in) :: Ndkx !< total number of flownodes (dynamically changing)
integer, intent(in) :: Lnkx !< total number of flowlinks (dynamically changing)
integer, dimension(Ndx), intent(in) :: kbot !< flow-node based layer administration
integer, dimension(Ndx), intent(in) :: ktop !< flow-node based layer administration
integer, dimension(Lnx), intent(in) :: Lbot !< flow-link based layer administration
integer, dimension(Lnx), intent(in) :: Ltop !< flow-link based layer administration
double precision, dimension(Lnkx), intent(in) :: q1 !< flow-field discharges
double precision, dimension(Ndkx), intent(in) :: qw !< flow-field vertical discharges
double precision, dimension(Ndkx), intent(out) :: sq !< flux balance (inward positive)
double precision :: dum, sumsq
integer :: k1, k2, k, kk, L, LL
sq = 0d0
! write(6,*) q1(1714)-q1(1735)
! write(6,*) qw(1736)-qw(1735)
! vertical
do kk=1,Ndx
do k = kbot(kk),ktop(kk)-1
sq(k) = sq(k) - qw(k)
sq(k+1) = sq(k+1) + qw(k)
end do
if ( abs(qw(ktop(kk))).gt.1d-8 ) then
continue
end if
if ( abs(qw(kbot(kk)-1)).gt.1d-8 ) then
continue
end if
end do
! write(6,"('after vertical: ',E)") sq(1736)
! horizontal
do LL=1,Lnx
do L=Lbot(LL),Ltop(LL)
k1 = ln(1,L)
k2 = ln(2,L)
sq(k1) = sq(k1) - q1(L)
sq(k2) = sq(k2) + q1(L)
end do
end do
! write(6,"('after horizontal: ',E)") sq(1736)
sumsq = 0d0
do kk=1,Ndx
do k=kbot(kk),ktop(kk)
sumsq = sumsq + sq(k)
end do
end do
! write(6,*) sumsq
return
end subroutine comp_sq
subroutine drop_tracer(xp, yp, dval)
use m_transport
use m_flowgeom
use m_flow, only: kmx
use m_polygon
implicit none
double precision :: xp, yp !< point coordinates
double precision :: dval !< value
integer, dimension(:), allocatable :: icelllist
integer :: Ncells
integer :: i, k, kk, kb, kt
integer :: N, in
integer :: ja
integer :: iconst
! allocate
allocate(icelllist(Ndx))
icelllist = 0
! add a tracer if no current tracer is selected (in visualization)
if ( ITRA1.eq.0 .or. iconst_cur.eq.0 .or. iconst_cur.lt.ITRA1 ) then ! note: tracers always at the back
call add_tracer('', iconst)
! set current tracer (for visualization)
iconst_cur = iconst
else
! select current tracer
iconst = iconst_cur
end if
! find active flow nodes
Ncells = 0
if ( NPL.le.2 ) then ! no (usable) polygon
call in_flowcell(xp,yp,kk)
if ( kk.gt.0 ) then
Ncells = Ncells + 1
icelllist(1) = kk
end if
else
in = -1
do kk=1,Ndx
N = size(nd(kk)%x)
call dbpinpol(xz(kk), yz(kk), in)
if ( in.eq.1 ) then
Ncells = Ncells+1
icelllist(Ncells) = kk
end if
end do
end if
! fill active flow nodes
do i=1,Ncells
kk = icelllist(i)
call getkbotktop(kk,kb,kt)
do k=kb,kb+max(kmx,1)-1
constituents(iconst,k) = constituents(iconst,k) + dval
end do
end do
! plot
call tekflowstuff(ja)
! deallocate
if ( allocated(icelllist) ) deallocate(icelllist)
return
end subroutine drop_tracer
!> find index of string in array of strings
integer function findname(N, snames, sname)
implicit none
integer, intent(in) :: N
character(len=*), dimension(N), intent(in) :: snames
character(len=*), intent(in) :: sname
integer :: i
findname = 0
do i=1,N
if ( trim(sname).eq.trim(snames(i)) ) then
findname = i
return
end if
end do
return
end function findname
!> Convert qid (from .ext file) to tracer name (split in generic qidname and specific tracer name).
!! If the input qid is not tracer, then the same qid is returned (and no tracer name)
subroutine get_tracername(qid, trname, qidname)
use m_transport, only: DEFTRACER
implicit none
character(len=*), intent(in) :: qid !< Original quantityid, e.g., 'tracerbndfluor'.
character(len=*), intent(out) :: trname !< The trimmed tracer name, e.g., 'fluor'.
character(len=*), intent(out) :: qidname !< The base quantity name for further use in external forcing, e.g., 'tracerbnd'.
trname = ''
qidname = qid
if ( qid(1:9).eq.'tracerbnd' ) then
qidname = qid(1:9)
if ( len_trim(qid).gt.9 ) then
trname = trim(qid(10:))
else
trname = trim(DEFTRACER)
end if
else if (qid(1:13).eq.'initialtracer' ) then
qidname = qid(1:13)
if ( len_trim(qid).gt.13 ) then
trname = trim(qid(14:))
else
trname = trim(DEFTRACER)
end if
end if
return
end subroutine get_tracername
!> apply tracer boundary conditions
subroutine apply_tracer_bc()
use m_transport
use m_meteo
use m_flowgeom, only: ln
use m_flow, only: kmxd, q1
implicit none
character (len=NAMTRACLEN) :: tracnam
integer :: itrac, iconst
integer :: k, kk, ki, kb
integer :: L, LL, Lb, Lt
! loop over the tracer boundary conditions
do itrac=1,numtracers
iconst = itrac2const(itrac)
do k=1,nbndtr(itrac)
LL = bndtr(itrac)%k(3,k)
call getLbotLtop(LL,Lb,Lt)
do L = Lb,Lt
kb = ln(1,L)
ki = ln(2,L)
if ( q1(L).gt.0 ) then ! inflow
kk = kmxd*(k-1)+L-Lb+1
constituents(iconst,kb) = bndtr(itrac)%z(kk)
else ! outflow
constituents(iconst,kb) = constituents(iconst,ki)
end if
end do
end do
end do
return
end subroutine apply_tracer_bc
! begin DEBUG
subroutine dum_dostep() ! dit heb ik hier gezet omdat deze code veel lijkt op unstruc code, ik raakte in de war
use m_flowgeom
use m_flow
use m_flowparameters
use m_flowtimes
use m_physcoef, only : difmolsal
implicit none
double precision, dimension(:), allocatable :: sa_temp, DsedDt
limtypsed = 4 ! 4
call sethu(0)
call setau()
call setucxucyucxuucyu()
call setkbotktop(0)
call setdtorg()
! allocate
allocate(sa_temp(Ndkx), DsedDt(Ndkx))
!call update_sed(1, limtypsa, (/ difmolsal /), dts, sa1)
time1 = time0 + dts
time0 = time1
! deallocate
if ( allocated(sa_temp) ) deallocate(sa_temp)
if ( allocated(DsedDt) ) deallocate(DsedDt)
return
end subroutine dum_dostep
subroutine dum_makeflowfield()
use m_flowgeom
use m_flow
use m_flowtimes
implicit none
double precision :: x0 = 0.5d0
double precision :: y0 = 0.5d0
double precision :: ux, uy
integer :: kk, k, kb, kt, k1, k2, LL, L, Lb, Lt
squ = 0d0
do LL=1,Lnx
ux = -(yu(LL)-y0)
uy = xu(LL)-x0
!ux = 1d0
!uy = 0d0
call getLbotLtop(LL,Lb,Lt)
Lb = Lbot(LL)
do L=Lb,Lt
u1(L) = csu(LL)*ux + snu(LL)*uy
end do
!do L = Lt+1, Lb + kmxL(LL) - 1 ! copy top inactive part of column == utop
! u1(L) = u1(Lt)
!enddo
end do
if ( kmx.eq.0 ) then
do L=1,Lnx
q1(L) = au(L)*u1(L)
k1 = ln(1,L)
k2 = ln(2,L)
squ(k1) = squ(k1) + max(q1(L),0d0)
sqi(k1) = sqi(k1) - min(q1(L),0d0)
squ(k2) = squ(k2) - min(q1(L),0d0)
sqi(k2) = sqi(k2) + max(q1(L),0d0)
end do
else
do LL=1,Lnx
q1(LL) = 0d0
u1(LL) = 0d0
au(LL) = 0d0
call getlbotltop(LL,Lb,Lt)
do L=Lb,Lt
q1(L) = au(L)*u1(L)
q1(LL) = q1(LL) + q1(L) ! depth integrated result
au(LL) = au(LL) + au(L)
k1 = ln(1,L)
k2 = ln(2,L)
squ(k1) = squ(k1) + max(q1(L),0d0)
sqi(k1) = sqi(k1) - min(q1(L),0d0)
squ(k2) = squ(k2) - min(q1(L),0d0)
sqi(k2) = sqi(k2) + max(q1(L),0d0)
end do
if (au(LL) > 0d0 ) then ! depth averaged velocity
u1(LL) = q1(LL) / au(LL)
endif
end do
if (ja_timestep_auto >= 3) then ! 2D timestep
squ2d = squ ! needs to be larger than eps10 for setdtorg
if (ja_timestep_auto >= 4) then
squ2d = squ2d + sqi
endif
squ2d(1:Ndx) = 1d0
endif
end if
if ( kmx.gt.0 ) then
! close flow-field, assume zero flux through water surface
do kk=1,Ndxi
call getkbotktop(kk,kb,kt)
qw(kt) = 0d0
do k=kt,kb,-1
qw(k-1) = qw(k) + squ(k) - sqi(k)
end do
end do
end if
sq = sqi - squ !< note: inward fluxes positive
return
end subroutine dum_makeflowfield
subroutine dum_makesal()
use m_flowgeom
use m_flow
use m_flowparameters
implicit none
double precision :: x0 = 0.25d0
double precision :: y0 = 0.25d0
double precision :: sigma0 = 0.25d0
double precision :: R= 0.125d0
double precision :: xx, yy, zz
double precision :: R2, sa, sigma
integer :: kk, k, kb, kt
jasal = 1
if ( allocated(sa1) ) deallocate(sa1)
allocate(sa1(Ndkx))
if ( allocated(sa0) ) deallocate(sa0)
allocate(sa0(Ndkx))
R2 = R**2
sa1 = 0d0
do kk=1,Ndx
xx = xz(kk)-x0
yy = yz(kk)-y0
sa = (0.5d0+sign(0.5d0,R2-(xx**2+yy**2)))
! sa = xx
if ( kmx.eq.0 ) then
sa1(kk) = sa
else
call getkbotktop(kk,kb,kt)
do k=kb,kt
sigma = (0.5d0*(zws(k-1)+zws(k))-zws(kb-1))/(zws(kt)-zws(kb))
zz = sigma-sigma0
sa1(k) = (0.5d0+sign(0.5d0,R2-(xx**2+yy**2+zz**2)))
end do
end if
end do
if ( kmx.gt.0 ) then
do kk=1,Ndx
sa1(kk) = 0d0
call getkbotktop(kk,kb,kt)
do k=kb,kt
sa1(kk) = sa1(kk) + vol1(k)*sa1(k)
if ( sa1(kk).ne.0d0 .and. k.eq.kb ) then
continue
end if
end do
sa1(kk) = sa1(kk)/vol1(kk)
end do
end if
return
end subroutine dum_makesal
! end DEBUG
!> get maximum timestep for water columns (see setdtorg)
subroutine get_dtmax()
use m_flowgeom, only: Ndx, bl, nd, Dx
use m_flow, only: s1, epshu, squ, sqi, vol1, kmx, u1, hu
use m_flowparameters, only: eps10, cflmx
use m_flowtimes, only: time1
use m_timer
use m_transport
use m_partitioninfo
implicit none
integer :: kk, k, kb, kt
integer :: L, LL
integer :: ierror
double precision, parameter :: dtmax_default = 1d4
dtmin_transp = huge(1d0)
kk_dtmin = 0
if ( kmx.eq.0 ) then
do k=1,Ndx
dtmax(k) = dtmax_default
if ( s1(k)-bl(k).gt.epshu ) then
if ( squ(k).gt.eps10 ) then
dtmax(k) = min(dtmax(k),cflmx*vol1(k)/squ(k))
end if
! BEGIN DEBUG
! do LL=1,nd(k)%lnx
! L = iabs(nd(k)%ln(LL))
! if ( hu(L).gt.0d0 .and. u1(L).gt.0d0 ) then
! dtmax(k) = min(dtmax(k),cflmx*Dx(L)/u1(L))
! end if
! end do
! END DEBUG
if ( jampi.eq.1 ) then
! do not include ghost cells
if ( idomain(k).ne.my_rank ) cycle
end if
if ( dtmax(k).lt.dtmin_transp ) then
dtmin_transp = dtmax(k)
kk_dtmin = k
end if
end if
end do
else
do kk=1,Ndx
dtmax(kk) = dtmax_default
if ( s1(kk)-bl(kk).gt.epshu ) then
call getkbotktop(kk,kb,kt)
do k=kb,kt
if ( squ(k).gt.eps10 .or. sqi(k).gt.eps10 ) then
dtmax(kk) = min(dtmax(kk),vol1(k)/max(squ(k),sqi(k)))
end if
end do
dtmax(kk) = cflmx*dtmax(kk)
if ( jampi.eq.1 ) then
! do not include ghost cells
if ( idomain(kk).ne.my_rank ) cycle
end if
if ( dtmax(kk).lt.dtmin_transp ) then
dtmin_transp = dtmax(kk)
kk_dtmin = kk
end if
end if
end do
end if
time_dtmax = time1
if ( jampi.eq.1 ) then
! update dtmax
call update_ghosts(ITYPE_Sall, 1, Ndx, dtmax, ierror)
! globally reduce maximum time-step
if ( jatimer.eq.1 ) call starttimer(IMPIREDUCE)
call reduce_double_min(dtmin_transp)
if ( jatimer.eq.1 ) call stoptimer(IMPIREDUCE)
end if
return
end subroutine get_dtmax
!> get number of subtimesteps and delta subtimesteps
subroutine get_ndeltasteps()
use m_flowgeom, only: Ndxi, Lnxi, Lnx, ln
use m_flowtimes, only: dts
use m_transport
implicit none
double precision :: dt, dtmin
double precision :: logtwo
integer :: kk, LL
double precision, external :: get_dt
numnonglobal = 0
! get smallest and largest time steps
dtmin = dtmin_transp
if ( dtmin.ge.dts ) then
nsubsteps = 1
ndeltasteps = 1
else
logtwo = log(2d0)
nsubsteps = max(1,2**int(log(dts/dtmin)/logtwo+0.9999d0))
dtmin = dts/nsubsteps
! get number of substeps
do kk=1,Ndxi
dt = dtmax(kk)
if ( dt.lt.dts ) then
ndeltasteps(kk) = min(2**int(log(dt/dtmin)/logtwo),nsubsteps)
numnonglobal = numnonglobal+1
else
ndeltasteps(kk) = nsubsteps
end if
end do
! fictitious boundary cells
do LL=Lnxi+1,Lnx
ndeltasteps(ln(1,LL)) = ndeltasteps(ln(2,LL))
end do
! if ( nsubsteps.gt.1 ) then
! write(6,*) dtmin
! end if
end if
return
end subroutine get_ndeltasteps
!> sum horizontal fluxes
subroutine comp_sumhorflux(NUMCONST, kmx, Lnkx, Ndkx, Lbot, Ltop, fluxhor, sumhorflux)
use m_flowgeom, only: Lnx, Ln, nd, Ndx ! static mesh information
implicit none
integer, intent(in) :: NUMCONST !< number of transported quantities
integer, intent(in) :: kmx !< number of layers
integer, intent(in) :: Ndkx !< total number of flownodes (dynamically changing)
integer, intent(in) :: Lnkx !< total number of flowlinks (dynamically changing)
integer, dimension(Lnx), intent(in) :: Lbot !< flow-link based layer administration
integer, dimension(Lnx), intent(in) :: Ltop !< flow-link based layer administration
double precision, dimension(NUMCONST,Lnkx), intent(in) :: fluxhor !< horizontal advection fluxes
double precision, dimension(NUMCONST,Ndkx), intent(inout) :: sumhorflux ! sum of horizontal fluxes, dim(NUMCONST,Ndkx)
integer :: LL, L, Lb, Lt
integer :: j, k1, k2
integer :: k
if ( kmx.lt.1 ) then
! add horizontal fluxes to right-hand side
do L=1,Lnx
! get neighboring flownodes
k1 = ln(1,L)
k2 = ln(2,L)
do j=1,NUMCONST
sumhorflux(j,k1) = sumhorflux(j,k1) - fluxhor(j,L)
sumhorflux(j,k2) = sumhorflux(j,k2) + fluxhor(j,L)
end do
end do
else
! add horizontal fluxes to right-hand side
do LL=1,Lnx
Lb = Lbot(LL)
Lt = Ltop(LL)
do L=Lb,Lt
! get neighboring flownodes
k1 = ln(1,L)
k2 = ln(2,L)
do j=1,NUMCONST
sumhorflux(j,k1) = sumhorflux(j,k1) - fluxhor(j,L)
sumhorflux(j,k2) = sumhorflux(j,k2) + fluxhor(j,L)
end do
end do
end do
end if
return
end subroutine comp_sumhorflux
!> determine if the cells have to be updated (1) or not (0)
subroutine get_jaupdate(istep,nsubsteps,Ndxi,Ndx,ndeltasteps,jaupdate)
implicit none
integer, intent(in) :: istep !< substep number
integer, intent(in) :: nsubsteps !< number of substeps
integer, intent(in) :: Ndxi !< number of cells, excluding virtual boundary cells
integer, intent(in) :: Ndx !< number of cells, including virtual boundary cells
integer, dimension(Ndx), intent(in) :: ndeltasteps !< number of substeps between updates
integer, dimension(Ndx), intent(out) :: jaupdate !< update cell (1) or not (0)
integer :: kk
integer :: num
jaupdate = 0
num = 0
do kk=1,Ndxi
if ( mod(istep+1, ndeltasteps(kk)).eq.0 ) then
! if ( int((istep+1)/ndeltasteps(kk))*ndeltasteps(kk).eq.istep+1 ) then
jaupdate(kk) = 1
num = num+1
end if
end do
! BEGIN DEBUG
! jaupdate = 1
! END DEBUG
! if ( istep.lt.nsubsteps ) then
! write(6,"(I0,':',I0, ' ', $)") istep+1, num
! end if
return
end subroutine get_jaupdate
!> determine if the horizontal fluxes have to be updated (1) or not (0) from cell-based mask
subroutine get_jaupdatehorflux(nsubsteps, limtyp, jaupdate,jaupdatehorflux)
use m_flowgeom, only: Ndx, Lnx, ln, klnup
implicit none
integer, intent(in) :: nsubsteps !< number of substeps
integer, intent(in) :: limtyp !< limited higher-order upwind (>0) or first-order upwind (0)
integer, dimension(Ndx), intent(in) :: jaupdate !< cell updated (1) or not (0)
integer, dimension(Lnx), intent(out) :: jaupdatehorflux !< update horizontal flux (1) or not (0)
integer :: kk, k1, k2, LL
integer :: kk1L, kk2L, iswitchL
integer :: kk1R, kk2R, iswitchR
if ( nsubsteps.eq.1 ) then
jaupdatehorflux = 1
else
jaupdatehorflux = 0
if ( limtyp.eq.0 ) then
do LL=1,Lnx
k1 = ln(1,LL)
k2 = ln(2,LL)
if ( jaupdate(k1).eq.1 .or. jaupdate(k2).eq.1 ) then
jaupdatehorflux(LL) = 1
end if
end do
else
do LL=1,Lnx
k1 = ln(1,LL)
k2 = ln(2,LL)
if ( jaupdate(k1).eq.1 .or. jaupdate(k2).eq.1 ) then
jaupdatehorflux(LL) = 1
cycle
end if
kk1L = klnup(1,LL)
if ( kk1L.ne.0 ) then
if ( jaupdate(iabs(kk1L)).eq.1 ) then
jaupdatehorflux(LL) = 1
cycle
end if
if ( kk1L.gt.0 ) then
kk2L = klnup(2,LL)
if ( jaupdate(iabs(kk2L)).eq.1 ) then
jaupdatehorflux(LL) = 1
cycle
end if
end if
end if
kk1R = klnup(4,LL)
if ( kk1R.ne.0 ) then
if ( jaupdate(iabs(kk1R)).eq.1 ) then
jaupdatehorflux(LL) = 1
cycle
end if
if ( kk1R.gt.0 ) then
kk2R = klnup(5,LL)
if ( jaupdate(iabs(kk2R)).eq.1 ) then
jaupdatehorflux(LL) = 1
end if
end if
end if
end do
end if
end if
return
end subroutine get_jaupdatehorflux