!----- 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: modules.f90 43424 2015-12-04 17:30:45Z kernkam $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/modules.f90 $
!> @file modules.f90
!! Modules with global data.
!! call default_*() routines upon program startup and when loading a new MDU.
!! call only reset_*() routines when reinitialising an active flow model.
! NOTE: this document is automatically parsed
! CONVENTION:
! please use the following variable notation so the parser will pickup variables for dynamic exchange
! {=optional}
! typename, {allocatable, }target :: name{(:)} !< {(altname)} [units] description {JSON}
! NOTE: only one variable definition per line, the variable should not continue on the next line.
!
! The JSON part can contain additional key-value pairs in JSON format, e.g.:
! !< [m] waterlevel at previous timestep {"state":true,"slice":"1:nodtot","standard_name":"sea_surface_height"}
!
! For state variables values the following JSON key-value pairs are required:
! "standard_name" is the netcdf standard name for the variable, e.g. "sea_surface_height"
module m_physcoef
implicit none
double precision :: ag !< 10.0 ! 9.81 ! (m/s2)
double precision :: sag !< sqrt(ag)
double precision :: vonkar !< von Karman constant ()
double precision :: frcuni !< uniform friction coeff 2D
double precision :: frcuni1D !< uniform friction coeff 1D
double precision :: frcuni1D2D !< uniform friction coeff 1D2D
integer :: ifrctypuni !< 0=chezy, 1=manning, 2=white colebrook D3D, 3=white colebrook Waqua (now only 2D)
double precision :: frcunilin !<60. ! 6 ! 66 ! uniform friction coeff
double precision :: umodlin !< linear friction umod , ifrctyp 4,5,6
double precision :: wall_ks !< vertical wall nIKURADSE ROUGHNESSs (m)
double precision :: wall_z0 !< z0 for vertical walls, ~= Ks/30 (m)
!! z0 for bottom follows from ifrctyp==3 and z0=frcuni
double precision :: z0 !< z0
double precision :: ee !< natural e
double precision :: ee9 !< natural e/c9of1
double precision :: vicouv !< constant horizontal eddy viscosity (m2/s) mom
double precision :: dicouv !< constant horizontal eddy diffusivity (m2/s) sal, sed
double precision :: Elder !< add Elder viscosity
double precision :: Smagorinsky !< add Smagorinsky Cs coefficient, vic = vic + (Cs*dx)**2 * S
double precision :: viuchk !< if < 0.5 then eddy viscosity cell peclet check viu Sets ALL (scalar) variables in this module to their default values.
subroutine default_physcoef()
ag = 9.81d0 ! 10.0 ! (m/s2)
sag = sqrt(ag)
vonkar = 0.41d0 ! von Karman constant ()
ee = exp(1d0) ! natural e ()
ee9 = 9d0*ee !
frcuni = 0.023d0 ! 60. ! 6 ! 66 ! uniform friction coeff
frcuni1D = 0.023d0 ! 60. ! 6 ! 66 ! uniform friction coeff
frcuni1D2D = 0.023d0 ! 60. ! 6 ! 66 ! uniform friction coeff
ifrctypuni = 1 ! 0=chezy, 1=manning, 2=white colebrook (D3D), 3=white colebrook (WAQUA)
frcunilin = 0d0 !
umodlin = 1.d0 ! linear friction umod , ifrctyp 4,5,6
wall_ks = 0.0d0 ! vertical wall nIKURADSE ROUGHNESSs (m)
vicouv = 1d0 ! constant horizontal eddy viscosity (m2/s) mom
dicouv = 1d0 ! constant horizontal eddy diffusivity (m2/s) sal, sed
Elder = 0d0 ! add Elder viscosity
Smagorinsky = 0d0 ! add Smagorinsky Cs coefficient, vic = vic + (Cs*dx)**2 * S
viuchk = 0.24 ! if < 0.5 then eddy viscosity cell check viu Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, only call reset_waves() instead.
subroutine default_waves()
rouwav = ''
gammax = 1.0d0 !< Maximum wave height/water depth ratio
jatpwav = TPWAVDEFAULT !< TPWAV, TPWAVSMOOTH, TPWAVRELATIVE
call reset_waves()
end subroutine default_waves
!> Resets only waves variables intended for a restart of flow simulation.
!! Upon loading of new model/MDU, call default_waves() instead.
subroutine reset_waves()
extfor_wave_initialized = .false. !< is set to .true. when the "external forcing"-part that must be initialized for WAVE during running (instead of during initialization) has actually been initialized
end subroutine reset_waves
end module m_waves
module m_xbeach_data
use m_xbeach_typesandkinds
!==================================================================================================================================
! XBeach related variables
!==================================================================================================================================
!! Hydrodynamics arrays, allocatables
double precision, allocatable :: s1initial(:) !< initial waterlevel (m ) for Riemann boundary condition
double precision, allocatable :: ee0(:,:) !< wave energy at begin of timestep
double precision, allocatable :: ee1(:,:) !< wave energy at end of timestep
double precision, allocatable :: cwav(:) !< phase speed (m/s)
double precision, allocatable :: cgwav(:) !< wave group velocity (m/s)
double precision, allocatable :: kwav(:) !< wavenumber k (rad/m)
double precision, allocatable :: nwav(:) !< cg/c (-)
double precision, allocatable :: ctheta(:,:) !< propagation speed in theta space
double precision, allocatable :: sigmwav(:) !< wave frequency (rad/s)
double precision, allocatable :: sigt(:,:)
double precision, allocatable :: ee1sum(:) !< sum over directional bins
double precision, allocatable :: horadvec(:,:) !< horizontal advection
double precision, allocatable :: thetaadvec(:,:) !< directional advection
double precision, allocatable :: rrthetaadvec(:,:) !< directional advection roller
double precision, allocatable :: rrhoradvec(:,:) !< horz advection roller
double precision, allocatable :: rr(:,:) !< directional advection roller
double precision, allocatable :: csx(:)
double precision, allocatable :: snx(:)
double precision, allocatable :: H(:) !< golfhoogte, onafh van instat
double precision, allocatable :: E(:) !< bulk wave energy in nodes
double precision, allocatable :: DR(:) !< Bulk roller dissipation
double precision, allocatable :: R(:) !< Bulk roller energy
double precision, allocatable :: thet(:,:) !< centre angle dir bin in each node
double precision, allocatable :: costh(:,:)
double precision, allocatable :: sinth(:,:)
double precision, allocatable :: Sxx(:) !< Radiation stresses
double precision, allocatable :: Syy(:)
double precision, allocatable :: Sxy(:)
double precision, allocatable :: Fx(:) !< wave forces, on links
double precision, allocatable :: Fy(:)
double precision, allocatable :: urms(:) !< contribution for wave induced bed shear stress
double precision, allocatable :: ust(:) !< Stokes drift east
double precision, allocatable :: vst(:) !< Stokes drift north
double precision, allocatable :: thetamean(:) !< mean wave angle
double precision, allocatable :: Qb(:) !< Wave breaking proportion
double precision, allocatable :: D(:) !< Wave breaking dissipation
double precision, allocatable :: BR(:) !< Roller surface slope, also important for morph
double precision, allocatable :: uin(:) !< xcomponent incoming long wave induced velocity
double precision, allocatable :: vin(:) !< ycomponent incoming long wave induced velocity
double precision, allocatable :: bi(:) !< long wave component bichromatic bc
double precision, allocatable :: Ltemp(:)
double precision, allocatable :: L1(:)
double precision, allocatable :: e01(:)
double precision, allocatable :: tE(:), dataE(:), databi(:)
double precision, allocatable :: L0(:),khdisp(:),hdisp(:)
double precision :: newstatbc !< stationary bc generated
double precision :: xref0, yref0 !< reference coordinates phase shift bc
integer :: randomseed
integer, allocatable, dimension(:) :: kbndu2kbndw !< mapping velocity bc pts to wave bc pts
integer, allocatable, dimension(:) :: kbndw2kbndu !< mapping velocity bc pts to wave bc pts
integer, allocatable, dimension(:) :: kbndz2kbndw !< mapping water level bc pts to wave bc pts
!! Model parameters
!! 1. DFLOW specific
integer :: limtypw = 4 !< 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC, 21=central voor wave action transport
double precision :: dtmaxwav !< subtimestepping for xbeach wave driver
!! 2. XBeach specific
! Type name initialize ! [unit] (advanced/deprecated) description
! [Section] Physical processes
integer :: swave = -123 ! [-] Include short waves (1), exclude short waves (0)
integer :: lwave = -123 ! [-] Include short wave forcing on NLSW equations and boundary conditions (1), or exclude (0)
! [Section] Wave boundary condition parameters
character(slen) :: instat = 'abc' ! [-] Wave boundary condtion type
double precision :: taper = -123 ! [s] Spin-up time of wave boundary conditions, in morphological time
double precision :: Hrms = -123 ! [m] Hrms wave height for instat = 0,1,2,3
double precision :: Tm01 = -123 ! [s] (deprecated) Old name for Trep
double precision :: Trep = -123 ! [s] Representative wave period for instat = 0,1,2,3
double precision :: Tlong = -123 ! [s] Wave group period for case instat = 1
double precision :: dir0 = -123 ! [deg] Mean wave direction (Nautical convention) for instat = 0,1,2,3
double precision :: nwavmax = -123 ! [-] (advanced) maximum ratio of cg/c fro computing long wave boundary conditions
integer :: m = -123 ! [-] Power in cos^m directional distribution for instat = 0,1,2,3
! [Section] Wave-spectrum boundary condition parameters
character(slen) :: bcfile = 'abc' ! [-] Name of spectrum file
integer :: random = -123 ! [-] (advanced) Random seed on (1) or off (0) for instat = 4,5,6 boundary conditions
double precision :: fcutoff = -123 ! [Hz] (advanced) Low-freq cutoff frequency for instat = 4,5,6 boundary conditions
integer :: nspr = -123 ! [-] (advanced) nspr = 1 long wave direction forced into centres of short wave bins, nspr = 0 regular long wave spreadin
double precision :: trepfac = -123 ! [-] (advanced) Compute mean wave period over energy band: trepfac*maxval(Sf) for instat 4,5,6; converges to Tm01 for trepfac = 0.0 and
double precision :: sprdthr = -123 ! [-] (advanced) Threshold ratio to maxval of S above which spec dens are read in (default 0.08*maxval)
integer :: correctHm0 = -123 ! [-] (advanced) Turn off or on Hm0 correction
integer :: Tm01switch = -123 ! [-] (advanced) Turn off or on Tm01 or Tm-10 switch
double precision :: rt = -123 ! [s] Duration of wave spectrum at offshore boundary, in morphological time
double precision :: dtbc = -123 ! [s] (advanced) Timestep used to describe time series of wave energy and long wave flux at offshore boundary (not affected by morfac)
double precision :: dthetaS_XB = -123 ! [deg] (advanced) The (counter-clockwise) angle in the degrees needed to rotate from the x-axis in SWAN to the x-axis pointing East
integer :: nspectrumloc = -123 ! [-] (advanced) Number of input spectrum locations
integer :: oldnyq = -123 ! [-] (advanced) Turn off or on old nyquist switch
! [Section] Flow boundary condition parameters
integer :: order = -123 ! [-] (advanced) Switch for order of wave steering, 1 = first order wave steering (short wave energy only), 2 = second oder wave steering (bound long wave corresponding to short wave forcing is added)
integer :: freewave = -123 ! [-] (advanced) Switch for free wave propagation 0 = use cg (default); 1 = use sqrt(gh) in instat = 3
double precision :: epsi = -123 ! [-] (advanced) Ratio of mean current to time varying current through offshore boundary
character(slen) :: tidetype = 'abc' ! [-] (advanced) Switch for offfshore boundary, velocity boundary or instant water level boundary (default)
integer :: ARC = -123 ! [-] (advanced) Switch for active reflection compensation at seaward boundary: 0 = reflective, 1 = weakly (non) reflective
double precision :: hminlw = -123 ! [-] minimum depth for wave forcing in flow momentum equation RHS
! [Section] Wave breaking parameters
character(slen) :: break = 'abc' ! [-] Type of breaker formulation (1=roelvink, 2=baldock, 3=roelvink adapted, 4=roelvink on/off breaking)
double precision :: gamma = -123 ! [-] Breaker parameter in Baldock or Roelvink formulation
double precision :: gamma2 = -123 ! [-] End of breaking parameter in break = 4 formulation
double precision :: alpha = -123 ! [-] (advanced) Wave dissipation coefficient in Roelvink formulation
double precision :: nroelvink = -123 ! [-] (advanced) Power in Roelvink dissipation model
double precision :: gammax = -123 ! [-] (advanced) Maximum ratio wave height to water depth
double precision :: deltaH = -123 ! [-] (advanced) Fraction of wave height to add to water depth
double precision :: fw = -123 ! [-] (advanced) Bed friction factor
double precision :: fwcutoff = -123 ! [-] Depth greater than which the bed friction factor is NOT applied
integer :: breakerdelay = -123 ! [-] (advanced) Turn on (1) or off (0) breaker delay model
! [Section] Roller parameters
integer :: roller = -123 ! [-] (advanced) Turn on (1) or off(0) roller model
double precision :: beta = -123 ! [-] (advanced) Breaker slope coefficient in roller model
! [Section] Wave-current interaction parameters
integer :: wci = -123 ! [-] Turns on (1) or off (0) wave-current interaction
double precision :: hwci = -123 ! [m] (advanced) Minimum depth until which wave-current interaction is used
double precision :: cats = -123 ! [Trep] (advanced) Current averaging time scale for wci, in terms of mean wave periods
! [Section] Wave numerics parameters
double precision :: wavint = -123 ! [s] Interval between wave module calls (only in stationary wave mode)
double precision :: maxerror = -123 ! [m] (advanced) Maximum wave height error in wave stationary iteration
integer :: maxiter = -123 ! [-] (advanced) Maximum number of iterations in wave stationary
!
!! spectral_wave_bc_module in XBeach (waveparamsnew.f90)
!type filenames ! Place to store multiple file names
! character(slen) :: fname ! file name of boundary condition file
! integer :: listline ! read position in FILELIST files
! logical :: repeat = .false. ! indicate to reuse this file every rtbc cycle
!endtype filenames
!type(filenames),dimension(:),allocatable,save :: bcfiles ! input wave spectrum files
end module m_xbeach_data
module m_xbeach_avgoutput
double precision, allocatable :: H_mean(:), H_var(:), H_min(:), H_max(:), H_varcross(:), H_varsquare(:) !< Sign wave height
double precision, allocatable :: urms_mean(:), urms_var(:), urms_max(:), urms_min(:), urms_varcross(:), urms_varsquare(:) !< orbital velocity
double precision, allocatable :: ust_mean(:), ust_var(:), ust_max(:), ust_min(:), ust_varcross(:), ust_varsquare(:) !< orbital velocity
double precision, allocatable :: vst_mean(:), vst_var(:), vst_max(:), vst_min(:), vst_varcross(:), vst_varsquare(:) !< orbital velocity
double precision, allocatable :: Fx_mean(:),Fx_var(:), Fx_min(:), Fx_max(:) , Fx_varcross(:), Fx_varsquare(:) !< x-comp wave forcing
double precision, allocatable :: Fy_mean(:),Fy_var(:), Fy_min(:), Fy_max(:), Fy_varcross(:), Fy_varsquare(:) !< y-comp wave forcing
double precision, allocatable :: E_mean(:),E_var(:), E_min(:), E_max(:), E_varcross(:), E_varsquare(:) !< Bulk wave energy
double precision, allocatable :: R_mean(:),R_var(:), R_min(:), R_max(:), R_varcross(:), R_varsquare(:) !< Bulk roller energy
double precision, allocatable :: D_mean(:),D_var(:), D_min(:), D_max(:), D_varcross(:), D_varsquare(:) !< Bulk wave dissipation
double precision, allocatable :: DR_mean(:),DR_var(:), DR_min(:), DR_max(:), DR_varcross(:), DR_varsquare(:) !< Bulk roller dissipation
double precision, allocatable :: s1_mean(:),s1_var(:), s1_min(:), s1_max (:), s1_varcross(:), s1_varsquare(:) !< Water level
double precision, allocatable :: u_mean(:),u_var(:), u_min(:), u_max(:), u_varcross(:), u_varsquare(:) !< velocity
double precision, allocatable :: v_mean(:),v_var(:),v_min(:), v_max(:), v_varcross(:), v_varsquare(:) !< velocity
double precision, allocatable :: thetamean_mean(:),thetamean_var(:), thetamean_min(:), thetamean_max(:), &
thetamean_varcross(:), thetamean_varsquare(:), &
thetamean_sin(:), thetamean_cos(:)
double precision, allocatable :: cwav_mean(:),cwav_var(:), cwav_min(:), cwav_max(:), cwav_varcross(:), cwav_varsquare(:)
double precision, allocatable :: cgwav_mean(:),cgwav_var(:), cgwav_min(:), cgwav_max(:), cgwav_varcross(:), cgwav_varsquare(:)
double precision, allocatable :: sigmwav_mean(:),sigmwav_var(:), sigmwav_min(:), sigmwav_max(:), sigmwav_varcross(:), sigmwav_varsquare(:)
end module m_xbeach_avgoutput
!>
!! Trachytope module containing FM data for trachytopes.
!! Note that arrays are dimensioned on the number of net links
!! This was done so that roughness characteristics can be specified
!! not related to the location of open boundaries, which are not
!! yet available at the time the model is constructed.
!! Besides that certain sediment transport formulae require a
!! Cf at the flow node, which can only accurately be deterimined if the
!! values at net links are known.
!!
module m_trachy
use trachytopes_data_module
use properties
implicit none
type(TREE_DATA), pointer :: trtdef_ptr !< Property tree structure containing the trachytopes definitions from MDU [Trachytopes] chapter.
!
!type(trachy_type) :: trachy_nl !< Structure containing trachytope definitions on net links
type(trachy_type) :: trachy_fl !< Structure containing trachytope definitions on flow links
!
logical :: linit !< Logical for initial step of trachytopes computation (not used)
logical :: waqol !< Logical for waq-online coupling in trachytopes computation
logical :: lfdxx !< Logical for sediment diameters of trachytopes computation (not used yet)
logical :: spatial_bedform !< Logical for inclusion of spatially varying dune properties in trachytopes computation (not used yet)
logical :: update_umag !< Logical for updating cell-centred velocity magnitude in trachytopes computation
!
double precision, allocatable :: rhosol(:) !< Density of sediment (lsedtot)
double precision, allocatable :: sig(:) !< sigma layer notation as in Delft3D in trachytopes computation
double precision, allocatable :: umag(:) !< velocity magnitude in trachytopes computation (ndx)
double precision, allocatable :: bedformD50(:) !< 50-th percentile of the sediment considered for bedform in trachytopes computation (ndx)
double precision, allocatable :: bedformD90(:) !< 90-th percentile of the sediment considered for bedform in trachytopes computation (ndx)
double precision, allocatable :: rksr(:) !< roughness due to ripples in trachytopes computation (cf. Van Rijn 20..) (ndx)
double precision, allocatable :: rksmr(:) !< roughness due to mega-ripples in trachytopes computation (cf. Van Rijn 20..) (ndx)
double precision, allocatable :: rksd(:) !< roughness due to dunes in trachytopes computation (cf. Van Rijn 20..)(ndx)
double precision, allocatable :: dxx(:,:) !< sediment percentiles in trachytopes computation (cf. Van Rijn 20..) (ndx,nxx)
double precision, allocatable :: z0rou(:) !< z0rou in trachytopes computation (numl)
double precision, allocatable :: hu_trt(:) !< water depth on net links in trachytopes computation (numl)
double precision, allocatable :: dx_trt(:) !< length of net links in trachytopes computation (numl)
!
integer :: kmaxtrt !< number of sigma layers as in Delft3D in trachytopes computation
integer :: lsedtot !< number of sediment fractions in trachytopes computation
integer :: i50 !< index of sediment percentile in trachytopes computation
integer :: i90 !< index of sediment percentile in trachytopes computation
integer :: itimtt !< update frequency (multiple of dt_max, default = 1 --> every dt_max timestep)
integer :: nxx !< dimension of sediment percentiles in trachytopes computation
integer, allocatable :: kcu_trt(:) !< temporary array for kcu on net-links instead of flow-links in trachytopes computation (numl)
!
character(4) :: rouflo !< roughness method as described by Delft3D in trachytopes computation
!
contains
!> Sets ALL (scalar) variables in this module to their default values.
!! Make sure to call this at least once for each newly loaded model.
subroutine default_trachy()
call tree_destroy(trtdef_ptr)
linit = .false. !< Logical for initial step of trachytopes computation (not used)
waqol = .false. !< Logical for waq-online coupling in trachytopes computation
lfdxx = .false. !< Logical for sediment diameters of trachytopes computation (not used yet)
spatial_bedform = .false. !< Logical for inclusion of spatially varying dune properties in trachytopes computation (not used yet)
lsedtot = 1 !< number of sediment fractions in trachytopes computation
i50 = 1 !< index of sediment percentile in trachytopes computation
i90 = 2 !< index of sediment percentile in trachytopes computation
itimtt = 1 !< update frequency (multiple of dt_max, default = 1 --> every dt_max timestep)
nxx = 2 !< dimension of sediment percentiles in trachytopes computation
end subroutine default_trachy
end module m_trachy
module m_sediment
use precision, only: fp
use m_rdstm, only: stmtype
use morphology_data_module, only: sedtra_type
use m_waves
implicit none
!-------------------------------------------------- new sediment transport and morphology
type mortmpdummy
logical :: have_waves
real(fp), dimension(:), pointer :: hrms !< Temporary dummy variable for wave hrms
real(fp), dimension(:), pointer :: tp !< Temporary dummy variable for wave tp
real(fp), dimension(:), pointer :: teta !< Temporary dummy variable for wave teta
real(fp), dimension(:), pointer :: rlabda !< Temporary dummy variable for wave rlabda
real(fp), dimension(:), pointer :: uorb !< Temporary dummy variable for wave uorb
real(fp), dimension(:), pointer :: ubot !< Temporary dummy variable for wave ubot
real(fp), dimension(:), pointer :: rksr !< Temporary dummy variable for bedform rksr
real(fp), dimension(:), pointer :: rhowat !< Temporary dummy variable for flow rhowat
real(fp), dimension(:,:), pointer :: seddif !< vertical sediment diffusivity
real(fp), dimension(:,:), pointer :: caksrho !< equilibrium density in kg/m3 at reference height
real(fp), dimension(:,:), pointer :: sed !< sediment concentration
real(fp), dimension(:,:), pointer :: ws !< sediment settling velocity
real(fp), dimension(:), pointer :: depchg
end type mortmpdummy
!
logical :: stm_included !< Include sediment transport (and optionally morphology)
type(stmtype), target :: stmpar !< All relevant parameters for sediment-transport-morphology module.
type(sedtra_type) :: sedtra !< All sediment-transport-morphology fields.
type(mortmpdummy), target :: mtd !< Dummy quantities not yet available in D-Flow FM
!-------------------------------------------------- old sediment transport and morphology
integer :: jased !< Include sediment, 1=Krone, 2=Soulsby van Rijn 2007, 3=Bert's morphology module
double precision :: vismol !< molecular viscosity (m2/s)
integer :: mxgr !< nr of grainsizes
integer :: mxgrKrone !< mx grainsize index nr that followsKrone. Rest follows v.Rijn
double precision, allocatable :: D50(:) !< mean sand diameter (m) ! used only if Ws ==0
double precision, allocatable :: D90(:) !< 90percentile sand diameter (m) ! not in Krone Partheniades
double precision, allocatable :: rhosed(:) !< rho of sediment (kg/m3)
double precision, allocatable :: rhodelta(:) !< relative density diff (rhosed-rhomean)/rhomean ( )
double precision, allocatable :: dstar(:) !< dimensionless particle diameter( )
double precision, allocatable :: dstar03(:) !< dimensionless particle diameter( ) **-0.3d0
double precision, allocatable :: Ws(:) !< Fall velocity (m/s) ( used only if D50=0)
double precision, allocatable :: erosionpar(:) !< Pickup erosion parameter ( kg/(m2s) ) Krone
double precision, allocatable :: Ustcre2(:) !< ustar critic erosion **2 ( m2/s2)
double precision, allocatable :: sqsgd50(:) !< sqrt( ((s-1)gD50) ) (m/s)
double precision, allocatable :: Accr(:) ! save time
double precision, allocatable :: Awcr(:) ! save time, see below
double precision, allocatable :: Bwcr(:) ! save time, see below
double precision, allocatable :: D50ca(:), D50cb(:), D50wa(:), D50wb(:), D50wc(:) !< SvR definitions + user defined for < 0.000062 (m)
double precision, allocatable :: Uniformerodablethickness(:) !< Uniform erodable thickness per fraction (m)
double precision, allocatable :: sedini(:) !< uniform initial sedcon (kg/m3)
double precision :: rhobulkrhosed = 1650d0/2650d0 !< rho of bulk sand / rho of sedimentmaterial
double precision :: sedmax !< user defined max sediment concentration (kg/m3)
double precision :: dmorfac ! morphological acceleration factor() , 0.0 = no bottom update, 1.0 = realtime, 10.0 = ten times faster
double precision :: tmorfspinup ! time period without morfac
double precision :: alfabed=1d0 ! calibration par bed load
double precision :: alfasus=1d0 ! calibration par suspended load
double precision :: crefcav=20d0 ! cref / caverage in Engelund Hansen wse = ws*crefcav
integer :: jamorf ! 0 or 1 do morf
double precision, allocatable :: sed (:,:) !< sediment concentraton kg/m3 (mxgr,ndkx)
double precision, allocatable :: sedi (:,:) !< sediment concentraton increment, kg/m3 only needed for jaceneqtr == 2
double precision, allocatable :: sdupq (:,:) !< sediment flux kg/s
double precision, allocatable :: blinc (:) !< bottom level increment (m)
double precision, allocatable :: grainlay(:,:) !< spatial erodable grain layer thickness for each grain size fraction (m)
integer :: jagrainlayerthicknessspecified = 0 !< specified non-uniformly yes or no
integer :: isusandorbed = 2 !< Supended and or Bedload: 1= S, 2=S+B
integer :: jaceneqtr = 2 !< equilibrium transport in cell centre=1, in net nodes=2
integer :: jgrtek = 1 !< grainsize fraction nr to plot
integer :: numintverticaleinstein = 10 !< number of vertical intervals in einstein integrals
double precision, allocatable :: taucx (:) !< cell centre tau current in x dir N/m2
double precision, allocatable :: taucy (:) !< cell centre tau in y dir N/m2
double precision, allocatable :: tauu (:) !< link tau in link dir N/m2, to cx,cy through Perot
contains
subroutine default_sediment()
use m_physcoef
implicit none
mxgr = 0
mxgrKrone = 0
wavenikuradse = 0.01d0
z0wav = wavenikuradse / 30d0
sedmax = 30d0
dmorfac = 1d0
tmorfspinup = 0d0
alfabed = 1d0
alfasus = 1d0
jamorf = 0
end subroutine default_sediment
subroutine allocgrains() ! for all fractions:
use MessageHandling
use m_physcoef
implicit none
integer :: m
double precision :: Taucre
if (allocated (D50) ) then
deallocate (D50, rhosed, erosionpar, Ustcre2, Ws, sedini, Uniformerodablethickness, &
D50ca, D50cb, D50wa, D50wb, D50wc, Bwcr )
endif
if (mxgr == 0) return
m = mxgr
allocate (D50(m), rhosed(m), erosionpar(m), Ustcre2(m), Ws(m), sedini(m), Uniformerodablethickness(m), &
D50ca(m), D50cb(m), D50wa(m), D50wb(m), D50wc(m), Bwcr(m) )
D50 = 0.2d-3 ! 1d-3
rhosed = 2650.0
erosionpar = 1d-4 ! krone
Taucre = 0.3d0
Ustcre2 = Taucre/rhomean ! krone, i.e. taucre = 0.3
ws = 3d-4
sedini = 0d0
Uniformerodablethickness = 1d0
D50ca = 0.19d0
D50cb = 0.1d0
D50wa = 0.24d0
D50wb = 0.66d0
D50wc = 0.33d0
Bwcr = 0.33d0
end subroutine allocgrains
end module m_sediment
module m_grw
integer :: jagrw=0 !< include ground water
double precision, allocatable :: sgrw0(:) !< ground water level start
double precision, allocatable :: sgrw1(:) !< ground water level end of timestep
double precision, allocatable :: pgrw (:) !< pressure and plotting of sgrw
double precision, allocatable :: h_unsat(:) !< initial height unsaturated zone
double precision :: Conductivity = 0d-4 !< non dimensionless K conductivity saturated (m/s), Q = K*A*i (m3/s)
double precision :: Unsatfac = 1.0d0 !< reduction factor for conductivity in unsaturated zone
double precision :: h_aquiferuni = 20d0 !< uniform height of carrying layer
double precision :: h_unsatini = 0.2 !< initial level groundwater is bedlevel - h_unsatini
double precision :: h_capillair = 0.5 !< Capillary rising height (m)
double precision :: h_transfer = 0.1d0 !< uniform thickness (numerical) transfer zone grw <-> openw
double precision :: porosgrw = 0.35d0 !< porosity of soil = Vair / (Vsoil+Vair) , or,
!< porosity of soil = (Rhoparticle - Rhobulk) / Rhoparticle
!< e.g.
end module m_grw
module m_ship
integer :: nshiptxy = 0, iniship !< nr of ships / initialised 0,1
integer, allocatable :: kship(:) !< index array
double precision, target, allocatable :: xyship(:) !< new position or velocity provided by module
double precision, target, allocatable :: shx(:) !< [m] current position {"shape": ["nshiptxy"]}
double precision, target, allocatable :: shy(:) !< [m] current position {"shape": ["nshiptxy"]}
double precision, target, allocatable :: shi(:) !< [m] current position {"shape": ["nshiptxy"]}
double precision, target, allocatable :: shu(:), shv(:), sho(:) !< current velocity
double precision, target, allocatable :: zsp(:) !< [m] ship depth at flownodes {"shape": ["nshiptxy"]}
double precision, allocatable :: adve0(:) !< saved adve array at flowlinks
double precision, allocatable, target :: shL(:) !< [m] ship size L/2, B/2, D ! for now, fixed max nr =2 {"shape": [2]}
double precision, allocatable, target :: shB(:) !< [m] ship size L/2, B/2, D ! for now, fixed max nr =2 {"shape": [2]}
double precision, allocatable, target :: shd(:) !< [m] ship size L/2, B/2, D ! for now, fixed max nr =2 {"shape": [2]}
double precision :: epsi=1d0
double precision :: fx2(2) =0d0 , fy2(2) =0d0 , fm2(2) =0d0 !< pressure force in global coordinate sys (interacting with flow)
double precision :: fricx (2)=0d0, fricy (2)=0d0, fricm (2)=0d0 !< frictiuon force in global coordinate sys (interacting with flow)
double precision :: fricxe(2)=0d0, fricye (2)=0d0, fricme(2)=0d0 !< frictiuon force in global coordinate sys (interacting with flow) explicit
double precision :: fricxi(2)=0d0, fricyi (2)=0d0, fricmi(2)=0d0 !< frictiuon force in global coordinate sys (interacting with flow) implicit
double precision :: stuwx (2)=0d0, stuwy (2)=0d0, stuwm (2)=0d0 !< thrust force in global coordinate sys (interacting with flow)
double precision :: fextx (2)=0d0, fexty (2)=0d0, fextm (2)=0d0 !< external force in global coordinate sys ( not on flow)
double precision, allocatable, target :: stuw(:) !< [N] actual thrust force in ship dir {"shape": [2]}
double precision, allocatable, target :: fstuw(:) !< [-] thrust setting 0-1 {"shape": [2]}
double precision, allocatable, target :: stuwmx(:) !< [N] max thrust {"shape": [2]}
double precision, allocatable, target :: roer(:) !< [degree] actual rudder angle {"shape": [2]}
double precision, allocatable, target :: froer(:) !< [degree] actual rudder setting 0-1 {"shape": [2]}
double precision, allocatable, target :: roermx(:) !< [degree] max rudder angle {"shape": [2]}
double precision :: powermx(2) , speedmx(2) !< mx engine power (Hp on input, then Watts), max ship velocity (Knts on input, then m/s)
double precision :: deadw(2) , deadwi (2), checkdw(2) !< inertia (x,y), moment
double precision :: xmx, xmn, ymx, ymn !< minmax of shipping domain
double precision :: Trelax = 4d0, depmin = 18d0 !< relax period pressureforces (s), ships no deeper than depmi
double precision :: Cfskin = 0.0015d0 !< skin friction coefficient tau/rho=Cfskin*Udif**2
double precision :: alfahull = 0d0 !< 0d0 = pressure forcing just hydrostatic, 1.0 = plus correction previous step
double precision :: vicuship = 1d0 !< increase background eddy viscosity under ship
integer :: icontroltyp(2) = 3 !< 1 = prescribed t,x,y and flow blocakage sluides,
!< 2 = prescribed t,x,y, ship
!< 3 = prescribed t,u,v, ship
!< 4 = keycontrolled ship
!< 5 = keycontrolled ship plus gyring
integer :: japressurehull = 1 !< apply pressure field on hull yes/no
integer :: japrop = 1 !< apply propellor forces
integer :: jafric = 1 !< apply friction forces
integer :: jashfricimpl = 1 !< frcition forces on ship implicit yes/no
integer :: ithull, ithullmx
end module m_ship
module m_wind
implicit none
double precision, allocatable, target :: wx(:) !< [m/s] wind x velocity (m/s) at u point {"shape": ["lnx"]}
double precision, allocatable, target :: wy(:) !< [m/s] wind y velocity (m/s) at u point {"shape": ["lnx"]}
double precision, allocatable, target :: ec_pwxwy_x(:) !< Temporary array, for comparing EC-module to Meteo1.
double precision, allocatable, target :: ec_pwxwy_y(:) !< Temporary array, for comparing EC-module to Meteo1.
double precision, allocatable, target :: patm(:) !< atmospheric pressure user specified in (N/m2), internally reworked to (m2/s2)
!! so that it can be merged with tidep later and difpatm/dx = m/s2, saves 1 array , using mode = 'add'
double precision, allocatable, target :: rain(:) !< [mm/day] rain at xz,yz {"shape": ["ndx"]}
double precision, allocatable :: pwxwythc(:) !< airpressure_windx_windy_humidity_airtemperature_cloudiness
double precision, allocatable, target :: tair(:) !< air temperature (degC)
double precision, allocatable, target :: rhum(:) !< air relative humidity (%)
double precision, allocatable, target :: clou(:) !< air cloudiness (%)
double precision, allocatable, target :: qrad(:) !< solar radiation (W/m2)
double precision, allocatable :: heatsrc (:) !< resulting 2D or 3D heat source per cell (Km3/s)
double precision, allocatable :: heatsrc0(:) !< resulting 2D or 3D heat source per cell, only set at timeuser (Km3/s)
double precision, allocatable :: salsrc (:) !< salinity source per cell (pptm3/s)
double precision, allocatable :: cdwcof(:) !< wind stress cd coefficient () , only if jatemp ==5
integer , allocatable :: kcw (:) !< mask array
integer :: jawind !< use wind yes or no
integer :: japatm !< use patm yes or no
integer :: jarain !< use rain yes or no
integer :: jatair !< use air temperature yes or no
integer :: jarhum !< use relative humidity yes or no
integer :: jaclou !< use cloudiness yes or no
integer :: jasol !< use 1 = use solrad, 2 = use cloudiness
integer :: jaqin !< use qin , sum of all in fluxes
double precision :: windxav, windyav !< average wind for plotting
double precision :: windsp
double precision :: winddir !< deg from north sailor
double precision :: wsx
double precision :: wsy
double precision :: rhoair !< (kg/m3)
double precision :: PavBnd !< average ambient pressure (N/m2) for correction on open boundaries
double precision :: PavIni !< average ambient pressure (N/m2) for initial waterlevel correction
double precision :: paver !< Average ambient pressure (N/m2)
double precision :: patmfac !< 100 if Mbar, 1 if Pascal
double precision :: cdb(3) !< breakpoints cd function cd coefficient
double precision :: wdb(3) !< breakpoints cd function windspeed
integer :: ICdtyp !< 1=Const; 2=Smith&Banke; 3=Smith&Banks, 3 breakpoints; 4 = Whang 2004 !SPvdP: comment out-of-date
contains
!> Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, only call reset_wind() instead.
subroutine default_wind()
use m_physcoef, only : rhomean
windsp = 0
winddir = 90d0 !< deg from north sailor
rhoair = 1.2d0
paver = 101325.0
Pavini = 0d0
PavBnd = 0d0 !< default: no pressure correction on open boundaries.
!< choose ambient pressure on boundaries equal to overall standard ambient pressure
patmfac = 1d0 !< 100 if Mbar, 1 if Pascal
cdb(1) = 0.00063d0 !< first wind breakpoint
wdb(1) = 0
cdb(2) = 0.00723d0 !< second wind breakpoint
wdb(2) = 100
cdb(3) = 0.003d0 !< third wind breakpoint
wdb(3) = 30
icdtyp = 2
windxav = 0d0
windyav = 0d0
! Rain+qin not reset every re-init, only upon new MDU load, because rain can be
! enabled by user in MDU (for BMI use, even without rain in external forcings file)
jarain = 0 !< use rain yes or no
jaqin = 0 !< use qin , sum of all in fluxes
! Remaining of variables is handled in reset_wind()
call reset_wind()
end subroutine default_wind
!> Resets only wind variables intended for a restart of flow simulation.
!! Upon loading of new model/MDU, call default_wind() instead.
subroutine reset_wind()
jawind = 0 !< use wind yes or no
japatm = 0 !< use patm yes or no
end subroutine reset_wind
end module m_wind
module m_heatfluxes
implicit none
double precision :: albedo ! reflection coefficient of water () at average incidence angle of 60 deg,
! (albedo is .025 at angle 0 deg, 0.13 at angle 70 deg)
double precision :: em ! Emissivity ()
double precision :: cpa ! Specific heat air [J/kg/K]
double precision :: rcpa !
double precision :: cpw ! Specific heat water [J/kg/K]
double precision :: rcpi ! m3K/J
double precision :: stf ! Stefan's constant =5.67e-8 [W/m^2/K^4]
double precision :: emstf ! Em*Stf [W/m^2/K^4]
double precision :: tkelvn ! Absolute zero
double precision :: QSUNav ! Solar influx (W/m2)
double precision :: QEVAav ! Evaporative heat loss (W/m2)
double precision :: QCONav ! Convective heat loss (W/m2)
double precision :: QLongav ! Long wave back radiation (W/m2)
double precision :: Qfreeav ! Free conv + evap heat loss (W/m2)
double precision :: Qfrconav ! Free convection heat loss (W/m2)
double precision :: Qfrevaav ! Free evaporation heat loss (W/m2)
double precision :: sarea ! Only for excess temp model jatem=3, lake area
double precision :: fwind ! Only for excess temp model jatem=3, wind factor
integer :: jamapheatflux !< write heatfluxes to map
integer :: jaRichardsononoutput !< write Richardson nr to his
double precision, allocatable :: Qsunmap(:)
double precision, allocatable :: Qevamap(:)
double precision, allocatable :: Qconmap(:)
double precision, allocatable :: Qlongmap(:)
double precision, allocatable :: Qfrevamap(:)
double precision, allocatable :: Qfrconmap(:)
double precision, allocatable :: Qtotmap(:)
double precision, allocatable :: Rich(:)
contains
subroutine default_heatfluxes()
use m_physcoef, only : rhomean
use m_wind , only : rhoair
!< Heat flux model comnstants
albedo = 0.06d0 !< reflection coefficient of water () at average incidence angle of 60 deg,
!< (albedo is .025 at angle 0 deg, 0.13 at angle 70 deg)
em = 0.985d0 !< Emissivity ()
cpa = 1004d0 !< Specific heat air [J/kg/K]
rcpa = rhoair*cpa !
cpw = 3986d0 !< Specific heat water [J/kg/K]
rcpi = 1d0/(rhomean*cpw) !< m3K/J
stf = 5.67d-8 !< Stefan's constant =5.67e-8 [W/m^2/K^4]
emstf = em*stf
tkelvn = 273.15d0 !< Absolute 0
jamapheatflux = 0
jaRichardsononoutput = 0
end subroutine default_heatfluxes
end module m_heatfluxes
! todo: MERGE THIS WITH UNSTRUC_BOUNDARIES
module m_bnd !< boundary-type module
implicit none
integer, parameter :: NAMLEN = 128
type bndtype
character(len=NAMLEN) :: name !< boundary-type name
integer :: N !< number of boundary points
double precision, dimension(:), allocatable :: x !< inner node x-coordinates
double precision, dimension(:), allocatable :: y !< inner node y-coordinates
double precision, dimension(:), allocatable :: sigma !< sigma-values
double precision, dimension(:), allocatable :: z !< boundary condition values
double precision, dimension(:,:), allocatable :: xy2 !< outer-node (x,y)-coordinates
integer, dimension(:), allocatable :: kd !< boundary points
integer, dimension(:,:), allocatable :: k !< index array, see e.g. kbnd
double precision, dimension(:), allocatable :: tht !< Thatcher-Harleman outflow time
double precision, dimension(:), allocatable :: thz !< Thatcher-Harleman concentration
end type bndtype
contains
!> deallocate boundary-type
subroutine dealloc_bnd(bnd)
implicit none
type(bndtype), intent(inout) :: bnd !< boundary data
if ( allocated(bnd%x) ) deallocate(bnd%x)
if ( allocated(bnd%y) ) deallocate(bnd%y)
if ( allocated(bnd%sigma) ) deallocate(bnd%sigma )
if ( allocated(bnd%z) ) deallocate(bnd%z)
if ( allocated(bnd%xy2) ) deallocate(bnd%xy2)
if ( allocated(bnd%kd) ) deallocate(bnd%kd)
if ( allocated(bnd%k) ) deallocate(bnd%k)
return
end subroutine dealloc_bnd
!> (re)allocate boundary-type
subroutine alloc_bnd(N, kmx, bnd)
implicit none
integer, intent(in) :: N !< number of boundary points
integer, intent(in) :: kmx !< maximum number of layers
type(bndtype), intent(inout) :: bnd !< boundary data
call dealloc_bnd(bnd)
if ( N.gt.0 ) then
allocate(bnd%x(N))
allocate(bnd%y(N))
allocate(bnd%xy2(2,N))
allocate(bnd%kd(N))
allocate(bnd%k(5,N))
bnd%x = 0d0
bnd%y = 0d0
bnd%xy2 = 0d0
bnd%kd = 0
bnd%k = 0
if ( kmx.gt.0 ) then
allocate(bnd%sigma(kmx*N))
allocate(bnd%z(kmx*N))
bnd%sigma = 0d0
bnd%z = 0d0
else
allocate(bnd%z(N))
bnd%z = 0d0
end if
end if
return
end subroutine alloc_bnd
!> deallocate bndtr
subroutine dealloc_bndarr(bndarr)
implicit none
type(bndtype), dimension(:), allocatable, intent(inout) :: bndarr
integer :: i, numbnd
if ( allocated(bndarr) ) then
numbnd = ubound(bndarr,1)
do i=1,numbnd
call dealloc_bnd(bndarr(i))
end do
deallocate(bndarr)
end if
return
end subroutine dealloc_bndarr
end module m_bnd
! unstruc.f90
module m_flowexternalforcings
use m_wind
use m_bnd
implicit none
logical :: success !< want je wil maar liever succes
integer :: jatimespace !< doen ja/nee 1/0
integer :: mhis !< unit nr external forcings history *.exthis
integer, save :: kx, filetype, mext
character(len=256) :: qid
character(len=1) :: operand
integer :: numbnp !< total nr of open boundary cells for network extension
! For postprocessing, each boundary polyline is named as a open boundary section.
! For each section, all underlying netlinks are administered, along with the name of the pli file.
integer, parameter :: IBNDTP_UNKNOWN = -1 !< Uninitialized boundary type/fill value
integer, parameter :: IBNDTP_CLOSED = 0 !< Closed boundary (not used yet)
integer, parameter :: IBNDTP_ZETA = 1 !< Water level boundary
integer, parameter :: IBNDTP_U = 2 !< Velocity boundary (not detailed in q/u/normal/.. yet)
integer, parameter :: IBNDTP_1D2D = 3 !< Special 1D2D boundary
integer :: nopenbndsect !< Nr. of open boundary sections.
integer, allocatable :: openbndtype(:) !< (nopenbndsect) Type of this boundary section (one of IBNDTP_ZETA, etc...)
integer, allocatable :: nopenbndlin(:) !< (nopenbndsect) Nr. of links/exchanges per section.
integer, allocatable :: openbndlin(:) !< (sum(nopenbndlin(1:nopenbndsect)) Net link nrs for each open boundary.
character(len=256), allocatable :: openbndname(:) !< (nopenbndsec)
character(len=256), allocatable :: openbndfile(:) !< (nopenbndsec)
double precision :: transformcoef(25) !< Transform coefficients a+b*x
integer , allocatable :: kez (:) !< temp (numl) edge oriented z lev
integer , allocatable :: keu (:) !< temp (numl) edge oriented u vel
integer , allocatable :: kes (:) !< temp (numl) edge oriented s sal
integer , allocatable :: ketm (:) !< temp (numl) edge oriented tm sem
integer , allocatable :: kesd (:) !< temp (numl) edge oriented sd sed
integer , allocatable :: ket (:) !< temp (numl) edge oriented t tang. vel.
integer , allocatable :: keuxy(:) !< temp (numl) edge oriented uxuy vel.
integer , allocatable :: ken (:) !< temp (numl) edge oriented n normal vel.
integer , allocatable, target :: ke1d2d(:) !< temp (numl) edge oriented 1d2d bnd
integer , allocatable :: keg (:) !< temp (numl) edge oriented g gate
integer , allocatable :: ked (:) !< temp (numl) edge oriented d cdam
integer , allocatable :: kegen(:) !< temp (numl) edge oriented general structure
integer , allocatable :: kegs (:) !< temp (numl) edge oriented general structure new style
integer , allocatable :: kep (:) !< temp (numl) edge oriented p pump
integer , allocatable :: kew (:) !< temp (numl) edge oriented w waves
integer , allocatable :: ketr (:,:) !< temp (numl) edge oriented tracer
integer, allocatable :: itpez(:) !< temp (numl) edge oriented,
!! 1,*=boundary typ, see type indicator kbndz(4,*) below
integer, allocatable :: itpeu(:) !< hulp (numl) edge oriented,
!! 1,*=boundary typ, see type indicator kbndz(4,*) below
integer, allocatable :: itpenz(:) !< temp (numl) edge oriented,
!! 1,*=boundary number (nopenbndsect), see type indicator kbndz(5,*) below
integer, allocatable :: itpenu(:) !< temp (numl) edge oriented,
!! 1,*=boundary number (nopenbndsect), see type indicator kbndu(5,*) below
double precision, allocatable :: threttim(:,:) !< (NUMCONST,nopenbndsect) Thatcher-Harleman return times
character(len=256), allocatable :: thrtq(:) !< temp array for Thatcher-Harleman return time readout, stores constituents
double precision, allocatable :: thrtt(:) !< temp array for Thatcher-Harleman return time readout, stores return times
integer, allocatable :: thrtn(:) !< temp array for Thatcher-Harleman return time readout, stores cell indices (first one)
integer :: nzbnd !< number of waterlevel boundary segments
integer :: nbndz !< waterlevel boundary points dimension
double precision, allocatable :: xbndz(:) !< waterlevel boundary points xcor
double precision, allocatable :: ybndz(:) !< waterlevel boundary points ycor
double precision, allocatable, target :: zbndz(:) !< waterlevel boundary points function
double precision, allocatable :: zbndz0(:) !< waterlevel boundary points function
double precision, allocatable :: xy2bndz(:,:) !< waterlevel boundary 'external tolerance point'
integer , allocatable :: kdz (:) !< waterlevel boundary points temp array
integer , allocatable :: kbndz(:,:) !< waterlevel boundary points index array
!! 1,* = index in s1 boundary point
!! 2,* = index in s1 first point on the inside
!! 3,* = index in u1 of their connecting link (always positive to the inside)
!! 4,* = type indicator :
!! 1 = waterlevel boundary
!! 2 = waterlevel neumann
!! 3 = velocity normal ingoing component
!! 4 = velocity flux boundary
!! 5 = velocity Riemann boundary
!! 6 = waterlevel outflow
!! 5,* = member of boundary number somuch of this type
double precision, allocatable :: zkbndz(:,:) !< only for jaceneqtr == 2 : left and right vertical netnode zk levels
integer :: nbndu !< velocity boundary points dimension
double precision, allocatable :: xbndu(:) !< velocity boundary points xcor
double precision, allocatable :: ybndu(:) !< velocity boundary points ycor
double precision, allocatable, target :: zbndu(:) !< velocity boundary points function
double precision, allocatable :: zbndu_store(:) !< store of velocity boundary conditions
double precision, allocatable :: xy2bndu(:,:) !< velocity boundary 'external tolerance point'
integer , allocatable :: kdu (:) !< velocity boundary points temp array
integer , allocatable :: kbndu(:,:) !< velocity boundary points index array, see lines above
integer , allocatable :: L1qbnd(:) !< first nbndu point in discharge bnd nqbnd
integer , allocatable :: L2qbnd(:) !< second nbndu point in discharge bnd nqbnd
double precision, allocatable :: at_all(:) !< "at" for all qbnd's, dim(nqbnd)
double precision, allocatable :: at_sum(:) !< "at" for all qbnd's, summed over all domains, dim(nqbnd)
double precision, allocatable :: wwssav_all(:,:) !< "wwav" and "ssav" for all qnbd's, dim(2,nqbnd)
double precision, allocatable :: wwssav_sum(:,:) !< "wwav" and "ssav" for all qnbd's, summed over all domains, dim(2,nqbnd)
integer :: japartqbnd !< one or more of the discharge boundaries is partitioned (1) or not (0)
double precision, allocatable :: huqbnd(:) !< hu used in normalised Manning discharge boundary condition, based on average water-level
integer :: nqbnd !<
double precision :: qbndhutrs = 0.1d0 !< only discharge bnd here if hu>qbndhutrs
double precision, allocatable :: zkbndu(:,:) !< only for jaceneqtr == 2 : left and right vertical netnode zk levels
integer :: nbnds !< salinity boundary points dimension in 1D and 2D
double precision, allocatable :: xbnds(:) !< salinity boundary points xcor
double precision, allocatable :: ybnds(:) !< salinity boundary points ycor
double precision, allocatable, target :: sigmabnds(:) !< salinity boundary points sigma coordinates (for now: dim = (nbnds*kmx) )
double precision, allocatable, target :: zbnds(:) !< salinity boundary points function
double precision, allocatable :: xy2bnds(:,:) !< salinity boundary 'external tolerance point'
integer , allocatable :: kds (:) !< satinity boundary points temp array
integer , allocatable :: kbnds(:,:) !< salinity boundary points index array, see lines above
double precision, allocatable :: thtbnds(:) !< salinity Thatcher-Harleman outflow times (dim = (nbnds))
double precision, allocatable :: thzbnds(:) !< salinity Thatcher-Harleman outflow concentrations (dim = (nbnds*kmx))
integer :: nbndw !< wave boundary points dimension
double precision, allocatable :: xbndw(:) !< wave boundary points xcor
double precision, allocatable :: ybndw(:) !< wave boundary points ycor
double precision, allocatable :: zbndw(:,:) !< wave boundary points function
double precision, allocatable :: xy2bndw(:,:) !< wave boundary 'external tolerance point'
integer , allocatable :: kdw (:) !< wave boundary points temp array
integer , allocatable :: kbndw(:,:) !< wave boundary points index array, see lines above
integer :: nbndtm !< temperature boundary points dimension
double precision, allocatable :: xbndtm(:) !< temperature boundary points xcor
double precision, allocatable :: ybndtm(:) !< temperature boundary points ycor
double precision, allocatable, target :: sigmabndTM(:) !< temperature boundary points sigma coordinates
double precision, allocatable, target :: zbndtm(:) !< temperature boundary points function
double precision, allocatable :: xy2bndtm(:,:) !< temperature external tolerance point'
integer , allocatable :: kdtm (:) !< temperature boundary points temp array
integer , allocatable :: kbndtm(:,:) !< temperature boundary points index array, see lines above
double precision, allocatable :: thtbndtm(:) !< temperature Thatcher-Harleman outflow times
double precision, allocatable :: thzbndtm(:) !< temperature Thatcher-Harleman outflow concentrations
integer :: nbndsd !< sediment boundary points dimension
double precision, allocatable :: xbndsd(:) !< sediment boundary points xcor
double precision, allocatable :: ybndsd(:) !< sediment boundary points ycor
double precision, allocatable, target :: zbndsd(:) !< sediment boundary points function
double precision, allocatable :: xy2bndsd(:,:) !< sediment boundary 'external tolerance point'
integer , allocatable :: kdsd (:) !< sediment boundary points temp array
integer , allocatable :: kbndsd(:,:) !< sediment boundary points index array, see lines above
double precision, allocatable :: thtbndsd(:) !< sediment Thatcher-Harleman outflow times
double precision, allocatable :: thzbndsd(:) !< sediment Thatcher-Harleman outflow concentrations
integer, allocatable :: nbndtr(:) !< tracer boundary points dimension
integer :: nbndtr_all !< all tracer boundary points dimension (max(nbndtr))
integer :: numtracers !< number of tracers with boundary conditions
integer, parameter :: NAMTRACLEN = 128
character(len=NAMTRACLEN), allocatable :: trnames(:)!< tracer names (boundary conditions only, used for look-up)
type(bndtype), allocatable, target :: bndtr(:)
integer :: nbndt !< tang.velocity boundary points dimension
double precision, allocatable :: xbndt(:) !< tang.velocity boundary points xcor
double precision, allocatable :: ybndt(:) !< tang.velocity boundary points ycor
double precision, allocatable, target :: zbndt(:) !< tang.velocity boundary points function
double precision, allocatable :: xy2bndt(:,:) !< tang.velocity boundary 'external tolerance point'
integer , allocatable :: kdt (:) !< tang.velocity boundary points temp array
integer , allocatable :: kbndt(:,:) !< tang.velocity boundary points index array, see lines above
integer :: nbnduxy !< uxuyadvectionvelocity boundary points dimension
double precision, allocatable :: xbnduxy(:) !< uxuyadvectionvelocity boundary points xcor
double precision, allocatable :: ybnduxy(:) !< uxuyadvectionvelocity boundary points ycor
double precision, allocatable, target :: sigmabnduxy(:) !< uxuyadvectionvelocity boundary points sigma coordinates (for now: dim = (nbnds*kmx) )
double precision, allocatable, target :: zbnduxy(:) !< uxuyadvectionvelocity boundary points function
double precision, allocatable :: xy2bnduxy(:,:) !< uxuyadvectionvelocity boundary 'external tolerance point'
integer , allocatable :: kduxy (:) !< uxuyadvectionvelocity boundary points temp array
integer , allocatable :: kbnduxy(:,:) !< uxuyadvectionvelocity boundary points index array, see lines above
integer :: nbndn !< norm.velocity boundary points dimension
double precision, allocatable :: xbndn(:) !< norm.velocity boundary points xcor
double precision, allocatable :: ybndn(:) !< norm.velocity boundary points ycor
double precision, allocatable, target :: zbndn(:) !< norm.velocity boundary points function
double precision, allocatable :: xy2bndn(:,:) !< norm.velocity boundary 'external tolerance point'
integer , allocatable :: kdn (:) !< norm.velocity boundary points temp array
integer , allocatable :: kbndn(:,:) !< norm.velocity boundary points index array, see lines above
integer :: ngate !< gates links dimension, to specify gate lower edge level
double precision, allocatable :: xgate(:) !< gates links xcor = xz(k1)
double precision, allocatable :: ygate(:) !< gates links ycor
double precision, allocatable, target :: zgate(:) !< gates lower_edge_level value
double precision, allocatable :: xy2gate(:,:) !< gates links second point xcor = xz(k2)
integer , allocatable :: kgate(:,:) !< gates links index array, see lines above
integer , allocatable, target :: kdg (:) !< helper for multiple_uni
integer , allocatable :: L1gatesg(:) !< first ngate point in gate signal ngatesg
integer , allocatable :: L2gatesg(:) !< second ngate point in gate signal ngatesg
integer :: ngatesg !< nr of gate signals specified
character(len=128), allocatable, target :: gate_ids(:)
integer :: ncdam !< nr of controllable dam points
double precision, allocatable :: xcdam(:) !< dam nodes xcor = xz(k1)
double precision, allocatable :: ycdam(:) !< dam nodes ycor
double precision, allocatable, target :: zcdam(:) !< dam nodes zvalue {"shape": ["ncdam"]}
double precision, allocatable :: xy2cdam(:,:) !< cdams links second point xcor = xz(k2)
integer , allocatable :: kcdam(:,:) !< cdams links index array, see lines above
integer , allocatable, target :: kdd(:) !< helper for multiple_uni_damlevel
integer , allocatable :: L1cdamsg(:) !< first ncdam point in cdam signal ncdamsg
integer , allocatable :: L2cdamsg(:) !< second ncdam point in cdam signal ncdamsg
integer :: ncdamsg !< nr of cdam signals specified
character(len=128), allocatable, target :: cdam_ids(:)
integer :: ncgen !< nr of controllable generalstr points
double precision, allocatable :: xcgen(:) !< generalstr nodes xcor = xz(k1)
double precision, allocatable :: ycgen(:) !< generalstr nodes ycor
double precision, allocatable, target :: zcgen(:) !< generalstr nodes zvalue (kx=3)
double precision, allocatable :: xy2cgen(:,:) !< cgen links second point xcor = xz(k2)
double precision, allocatable :: Fusav(:,:) !< only needed if gatedoorheight > 0 , dim = ncgen
double precision, allocatable :: Rusav(:,:) !< only needed if gatedoorheight > 0
double precision, allocatable :: Ausav(:,:) !< only needed if gatedoorheight > 0
integer , allocatable :: kcgen(:,:) !< cgen links index array, see lines above
!! 1,* = index in s1 point "left" of genstru
!! 2,* = index in s1 point "right" of genstru
!! 3,* = index in u1 of their connecting link (may point from #2 -> #1 if flow link is in opposite direction through the genstru polyline)
!! 4,* = pointer to general structure signal nr n
integer , allocatable, target :: kdgen(:) !< helper for multiple_uni_damlevel
integer , allocatable :: L1cgensg(:) !< first ncdam point in cdam signal ncdamsg
integer , allocatable :: L2cgensg(:) !< second ncdam point in cdam signal ncdamsg
integer :: ncgensg !< nr of cdam signals specified
integer, parameter :: ICGENTP_WEIR = 1 !< general structure type: a weir
integer, parameter :: ICGENTP_GATE = 2 !< general structure type: a gate
integer, parameter :: ICGENTP_GENSTRU = 3 !< general structure type: a true general structure
character(len=128), allocatable, target :: cgen_ids(:)
integer, allocatable :: cgen_type(:) !< (1:ngensg) The type for each general structure, one of ICGENTP_WEIR|GATE|GENSTRU
integer, allocatable :: cgen2str(:) !< (1:ngensg) Mapping from overall ngensg index to underlying structure index in either 1:nweirgen, 1:ngategen, or 1:ngenstru (inverse from *2cgen arrays below)
! The user may specify different 'gate'/'weir'/'generalstructure',
! and all are translated into a general structure (in computations
! and external forcings). To distinguish them, maintain counters for each.
! This should hold: ncgensg = nweirgen + ngategen + ngenstru
integer :: nweirgen !< nr of weirs in the generalstructure set
integer :: ngategen !< nr of gates in the generalstructure set
integer :: ngenstru !< nr of real general structures in the generalstructure set
integer , allocatable, target :: weir2cgen(:) !< (1:nweirgen) Mapping from weir number to underlying generalstructure number
integer , allocatable, target :: gate2cgen(:) !< (1:ngategen) Mapping from gate number to underlying generalstructure number
integer , allocatable, target :: genstru2cgen(:) !< (1:ngenstru) Mapping from true general structure number to underlying generalstructure number
integer :: npump !< nr of pump links
double precision, allocatable :: xpump(:) !< pump nodes xcor = xz(k1)
double precision, allocatable :: ypump(:) !< pump nodes ycor
double precision, allocatable, target :: qpump(:) !< pump discharge m3/s
double precision, allocatable :: xy2pump(:,:) !< pump links second point xcor = xz(k2)
integer , allocatable :: kpump(:,:) !< pump links index array, see lines above
integer , allocatable, target :: kdp(:) !< helper for multiple_uni_pump
integer , allocatable :: L1pumpsg(:) !< first npump point in pump signal npumpsg
integer , allocatable :: L2pumpsg(:) !< second npump point in pump signal npumpsg
integer :: npumpsg !< nr of pump signals specified
character(len=128), allocatable, target :: pump_ids(:)
integer :: nbndqh !< q-h boundary points dimension
double precision, allocatable :: xbndqh(:) !< q-h boundary points xcor
double precision, allocatable :: ybndqh(:) !< q-h boundary points ycor
double precision, allocatable :: zbndqh(:) !< q-h boundary points function
integer , allocatable :: kdqh (:) !< q-h boundary points temp array
integer , allocatable :: kbndqh(:,:) !< q-h boundary points index array
!! 1,* = index in s1 boundary point
!! 2,* = index in s1 first point on the inside
!! 3,* = index in u1 of their connecting link (always positive to the inside)
!! 4,* = type indicator :
!! 1 = waterlevel boundary
!! 2 = waterlevel neumann
!! 3 = velocity normal ingoing component
!! 4 = velocity flux boundary
!! 5 = velocity Riemann boundary
!! 6 = waterlevel outflow
!! 7 = q-h boundary
integer :: nqhbnd !< number of qh boundaries
integer , allocatable :: L1qhbnd(:) !< first nbndz point in discharge bnd nqbnd
integer , allocatable :: L2qhbnd(:) !< second nbndz point in discharge bnd nqbnd
double precision, allocatable, target :: qhbndz(:) !< temporary array for storing boundary values per qh boundary segment
double precision, allocatable :: atqh_all(:) !< temporary array for computing discharge through qh boundary per domain
double precision, allocatable :: atqh_sum(:) !< temporary array for computing total discharge through qh boundary
double precision :: qhrelax = 1d-2 !< relaxation factor for h signal
integer :: nwbnd !< number of wave-energy boundaries
character(len=255), dimension(:), allocatable :: fnamwbnd !< polyline filenames associated with wave-energy boundary
integer :: numsrc !< nr of point sources/sinks
integer, allocatable :: ksrc(:,:) !< index array, 1=nodenr sink, 2 =kbsin , 3=ktsin, 4 = nodenr source, 5 =kbsor , 6=ktsor
double precision, allocatable :: qsrc(:) !< cell influx (m3/s) if negative: outflux
double precision, allocatable :: sasrc(:) !< salinity (ppt) if ksrc 3,4 == 0, else delta salinity
double precision, allocatable :: tmsrc(:) !< temperature (degC) if ksrc 3,4 == 0, else delta temperature
double precision, allocatable :: arsrc(:) !< pipe cross sectional area (m2). If 0, no net momentum
double precision, allocatable :: cssrc(:,:) !< (1:2,numsrc) cosine discharge dir pipe on start side (1) and end side (2) of pipe.
double precision, allocatable :: snsrc(:,:) !< (1:2,numsrc) sine discharge dir pipe on start side (1) and end side (2) of pipe.
double precision, allocatable :: zsrc (:,:) !< vertical level (m) bot
double precision, allocatable :: zsrc2(:,:) !< vertical level (m) top (optional)
double precision, allocatable :: srsn (:,:) !< 6,numsrc, to be reduced
integer, allocatable :: jamess(:) !< issue message mess for from or to point, 0, 1, 2
integer, allocatable, target :: kdss (:) !< helper for multiple_uni_discharge_salinity_temperature
double precision, allocatable, target :: qstss(:) !< array to catch multiple_uni_discharge_salinity_temperature
character(len=255), allocatable :: srcname(:) !< sources/sinks name (numsrc)
integer, allocatable :: ksrcwaq(:) !< index array, starting point in qsrcwaq
double precision, allocatable :: qsrcwaq (:) !< Cumulative qsrc within current waq-timestep
contains
!> Sets ALL (scalar) variables in this module to their default values.
!! For external forcings it is equivalent with default_flowexternalforcings().
subroutine reset_flowexternalforcings()
call default_flowexternalforcings()
end subroutine reset_flowexternalforcings
!> Resets external forcing variables intended for a restart of flow simulation.
!! For external forcings it is equivalent with reset_flowexternalforcings().
subroutine default_flowexternalforcings()
jatimespace = 0 ! doen ja/nee 1/0
mhis = 0 ! unit nr external forcings history *.exthis
numbnp = 0 ! total nr of open boundary cells for network extension
nopenbndsect = 0 ! Nr. of open boundary sections.
nbndz = 0 ! waterlevel boundary points dimension
nbndu = 0 ! velocity boundary points dimension
nbndqh = 0 ! q-h boundary points dimension
nbnds = 0 ! salinity boundary points dimension
nbndtm = 0 ! temperature boundary points dimension
nbndsd = 0 ! sediment boundary points dimension
nbndw = 0 ! JRE: wave boundary points dimension
nbndt = 0 ! tang.velocity boundary points dimension
nbnduxy = 0 ! uxuy adv vel bnd
nbndn = 0 ! norm.velocity boundary points dimension
ngate = 0 ! gates links dimension, to specify gate lower edge level
ngatesg = 0 ! nr of gate control signals
ncdam = 0 ! controllable dams nodes dimension, to specify local bottom level
ncdamsg = 0 ! nr of controllable dam signals
ncgen = 0 ! general structure nodes dimension, to apply gen struc defs
ncgensg = 0 ! nr of general structure signals
ncgen = 0 ! general structure nodes dimension, to apply gen struc defs
ncgensg = 0 ! nr of general structure signals
nweirgen = 0 ! nr of weirs in the generalstructure set
ngategen = 0 ! nr of gates in the generalstructure set
ngenstru = 0 ! nr of real general structures in the generalstructure set
npump = 0 ! npump dimension
npumpsg = 0 ! nr of pump signals
nqbnd = 0 ! nr of q bnd's
numsrc = 0
end subroutine default_flowexternalforcings
end module m_flowexternalforcings
module unstruc_channel_flow
use m_network
implicit none
type(NetworkType) :: network
contains
!> Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, only call reset_*() instead.
subroutine default_channel_flow()
call dealloc(network)
end subroutine default_channel_flow
end module unstruc_channel_flow
module m_structures
use properties
implicit none
type(tree_data), pointer, public :: strs_ptr !< A property list with all input structure specifications of the current model. Not the actual structure set.
integer :: jaoldstr !< tmp backwards comp: we cannot mix structures from EXT and from structure-input files. Use one or the other.
! Structure Parameters ! Nabi
double precision, dimension(:), allocatable :: cgendisch ! Array for structure discharge parameter (for output to his file)
double precision, dimension(:), allocatable :: pumpdisch ! Array for pump discharge parameter (for output to his file)
double precision, dimension(:), allocatable :: gatedisch ! Array for gate discharge parameter (for output to his file)
double precision, dimension(:), allocatable :: gategendisch ! Array for gate(new) discharge parameter (for output to his file)
double precision, dimension(:), allocatable :: gategenflowh ! Array for gate(new) flow through height parameter (for output to his file)
double precision, dimension(:), allocatable :: cdamdisch ! Array for dam discharge parameter (for output to his file)
double precision, dimension(:), allocatable :: weirdisch ! Array for weir discharge parameter (for output to his file)
double precision, dimension(:), allocatable :: pumpcapac ! Array for pump capacity parameter (for output to his file)
integer :: jahiscgen !< Write structure parameters to his file, 0: n0, 1: yes
integer :: jahispump !< Write pump parameters to his file, 0: n0, 1: yes
integer :: jahisgate !< Write gate parameters to his file, 0: n0, 1: yes
integer :: jahiscdam !< Write dam parameters to his file, 0: n0, 1: yes
integer :: jahisweir !< Write weir parameters to his file, 0: n0, 1: yes
integer, parameter :: IOPENDIR_FROMLEFT = -1 !< Gate door opens/closes from left side.
integer, parameter :: IOPENDIR_FROMRIGHT = 1 !< Gate door opens/closes from right side.
integer, parameter :: IOPENDIR_SYMMETRIC = 0 !< Gate door opens/closes symmetrically (from center).
type tgate !< Gate structure type, before it gets evaluated as a general structure.
!double precision :: sill_level !< Not used: stored in zcgen(1,igen)
!double precision :: lower_edge_level !< Not used: stored in zcgen(2,igen)
!double precision :: opening_width !< Not used: stored in zcgen(3,igen)
double precision :: door_height !< Height of the door, used for 'double-barrier' overflow. Time-INDEPENDENT.
double precision :: sill_width !< Width of the sill, may be larger than the opening width, such that in open part we have weir flow and in closed part we have gate flow. Time-INDEPENDENT.
integer :: opening_direction !< Direction from which the gate opens/closes, IOPENDIR_FROMLEFT|FROMRIGHT|SYMMETRIC.
end type tgate
! TIDAL TURBINES: Insert allocatable of type structure_turbines here
type(tgate), allocatable :: gates(:)
contains
subroutine init_structure_hisvalues()
use m_flowexternalforcings , only: npumpsg, ncgensg, ngatesg, ncdamsg, ngategen, ngenstru, nweirgen
jahiscgen = 1
jahispump = 1
jahisgate = 1
jahiscdam = 1
jahisweir = 1
if( jahispump > 0 .and. npumpsg > 0) then
if( allocated( pumpdisch ) ) deallocate( pumpdisch, pumpcapac )
allocate( pumpdisch(npumpsg) ) ; pumpdisch = 0d0
allocate( pumpcapac(npumpsg) ) ; pumpcapac = 0d0
endif
if( jahiscgen > 0 .and. ncgensg+ngenstru > 0 ) then
if( allocated( cgendisch ) ) deallocate( cgendisch )
allocate( cgendisch(ncgensg+ngenstru) ) ; cgendisch = 0d0
endif
if( jahisgate > 0 ) then
if( ngatesg > 0 ) then
if( allocated( gatedisch ) ) deallocate( gatedisch )
allocate( gatedisch(ngatesg) ) ; gatedisch = 0d0
endif
if( ngategen > 0 ) then
if( allocated( gategendisch ) ) deallocate( gategendisch, gategenflowh )
allocate( gategendisch(ngategen) ) ; gategendisch = 0d0
allocate( gategenflowh(ngategen) ) ; gategenflowh = 0d0
endif
endif
if( jahiscdam > 0 .and. ncdamsg > 0) then
if( allocated( cdamdisch) ) deallocate( cdamdisch )
allocate( cdamdisch(ncdamsg) ) ; cdamdisch = 0d0
endif
if( jahisweir > 0 .and. nweirgen > 0) then
if( allocated( weirdisch) ) deallocate( weirdisch )
allocate( weirdisch(nweirgen) ) ; weirdisch = 0d0
endif
! TIDAL TURBINES: Insert init_turbines here
end subroutine init_structure_hisvalues
!> Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, only call reset_structures() instead.
subroutine default_structures()
call tree_destroy(strs_ptr)
call reset_structures()
! TIDAL TURBINES: Insert calls to deallocate_turbines and init_turbines here
end subroutine default_structures
!> Resets only structures variables intended for a restart of an existing flow simulation (same MDU).
!! Upon loading of new model/MDU, call default_structures() instead.
subroutine reset_structures()
if (allocated(gates)) deallocate(gates)
end subroutine reset_structures
end module m_structures
module m_turbulence
implicit none
! Coefficients of k-e model:
double precision :: cmukep
double precision :: sqcmukep
double precision :: sqcmukepi
double precision :: cewall
double precision :: c1e
double precision :: c2e
double precision :: sigdif
double precision :: sigtke
double precision :: sigeps
double precision :: sigsal, sigsali
double precision :: sigtem, sigtemi
double precision :: sigsed, sigsedi
double precision :: sigrho
! c1e = c2e-vonkar**2/(sigeps*sqrt(cmukep))
double precision :: c1t
double precision :: c2t
double precision :: c3tsta
double precision :: c3tuns
double precision :: coefn2
double precision :: skmy
double precision :: a1ph
double precision :: a2
double precision :: b1
double precision :: b2
double precision :: c1
double precision :: e1
double precision :: e2
!double precision :: e2 = e1-1+vonkar**2*b1**0.6667/skmy
double precision :: ghmin
double precision :: ghmax
!quasi-equilibrium coefficients:
!double precision :: csb1 = a2*(1-6*a1ph/b1)
!double precision :: csb2 = 3*a2*(6*a1ph+b2)
!double precision :: csu1 = a1*(1-3*c1-6*a1ph/b1)
!double precision :: csu2 = 3*a1ph*a2*((b2-3*a2)*(1-6*a1ph/b1)-3*c1*(6*a1ph+b2))
!double precision :: csu3 = 3*a2*(6*a1ph+b2)
!double precision :: csu4 = 9*a1*a2
integer, parameter :: kmxx = 2000 ! max dim of nr of vertical layers
integer, parameter :: mg = 4 ! max dim of nr of sediment fractions
double precision :: dijdij (0:kmxx) ! dudz(k)**2+dvdz(k)**2 vertical shear squared
double precision :: buoflu ( kmxx)
double precision :: bruva ( kmxx)
double precision :: ak (0:kmxx) ! local arrays, (0:
double precision :: bk (0:kmxx)
double precision :: ck (0:kmxx)
double precision :: dk (0:kmxx)
double precision :: ek (0:kmxx)
double precision :: dz (0:kmxx)
double precision :: ucxref ( kmxx) ! and for reference/plotting:
double precision :: ucm ( kmxx) ! and for reference/plotting:
double precision :: dijdijref(0:kmxx)
double precision :: tkin1ref(0:kmxx)
double precision :: teps1ref(0:kmxx)
double precision :: vicwref(0:kmxx)
double precision :: hcref( kmxx) ! mid-layer heigths
double precision :: hwref(0:kmxx) ! layer interface height, 0=bed
double precision :: epstke = 1d-32 ! D3D: - 7, dpm: -32
double precision :: epseps = 1d-32 ! D3D: - 7, dpm: -32
double precision :: epsd = 1d-32 ! D3D: - 7, dpm: -32
double precision, allocatable :: turkin0 (:) ! k old (m2/s2) , at layer interface at u these will become global, rename to : turkinwu0
double precision, allocatable :: turkin1 (:) ! k new , at layer interface at u
double precision, allocatable :: tureps0 (:) ! eps old (1/s) , at layer interface at u
double precision, allocatable :: tureps1 (:) ! eps new , at layer interface at u
double precision, allocatable :: vicwwu (:) ! vertical eddy viscosity (m2/s) at layer interface at u point
double precision, allocatable :: vicwws (:) ! vertical eddy viscosity (m2/s) at layer interface at s point
real , allocatable :: tkepro (:) ! vertical production t
real , allocatable :: tkedis (:) ! vertical dissipation
double precision, allocatable :: rho (:) ! density at cell centres (kg/m3)
double precision, allocatable :: rho0 (:) ! density at cell centres (kg/m3), previous step
double precision, allocatable :: dpbdx0 (:) ! previous step baroclinic pressure gradient, at u points
double precision, allocatable :: rhou (:) ! density at flow links (kg/m3)
double precision, allocatable :: sigdifi (:) ! inverse prandtl schmidt nrs
double precision, allocatable :: turkinepsws (:,:) ! k and eps,1,2 at layer interface at c , horizontal transport of k and eps
double precision, allocatable :: tqcu(:) ! sum of q*turkinws at layer interface at cupw , horizontal transport of k and eps
double precision, allocatable :: eqcu(:) ! sum of q*turepsws at layer interface at cupw , horizontal transport of k and eps
double precision, allocatable :: sqcu(:) ! sum of q at layer interface at cupw , horizontal transport of k and eps
double precision, allocatable :: tttu(:), ttqc(:), tttc(:) ! test12
integer , allocatable :: ln0(:,:) ! links in transport trimmed to minimum of ktop,ktop0 for z-layers
contains
!> Sets ALL (scalar) variables in this module to their default values.
subroutine default_turbulence()
use m_physcoef
! Coefficients of k-e model:
sigdif = 1d0
sigtke = 1.0d0
sigeps = 1.3d0
sigrho = 0.7d0 ! bouyancy
sigsal = 0.7d0 ; sigsali = 1d0/sigsal
sigtem = 0.7d0 ; sigtemi = 1d0/sigtem
sigsed = 1.0d0 ; sigsedi = 1d0/sigsed
cmukep = 0.09d0
sqcmukep = sqrt(cmukep)
sqcmukepi = 1.d0 / sqcmukep
cewall = cmukep**0.75d0/vonkar ! 0.4769d0 ! 0.09**0.75/0.41 ! /vonkar
c2e = 1.92d0
c1e = 1.44d0
c1e = c2e-vonkar**2/(sigeps*sqcmukep)
c1t = 1d0-c1e
c2t = 1d0-c2e
c3tsta = 1d0
c3tuns = 1d0-c1e
coefn2 = - ag / (sigrho*rhomean)
skmy = 1.96d0
a1ph = 0.92d0
a2 = 0.74d0
b1 = 16.6d0
b2 = 10.1d0
c1 = 0.08d0
e1 = 1.80d0
e2 = 1.33d0
!e2 = e1-1+vonkar**2*b1**0.6667/skmy
ghmin =-0.280d0
ghmax = 0.0233d0
end subroutine default_turbulence
end module m_turbulence
module m_flowparameters
use m_sediment, only: jased
! use m_d3ddimens
implicit none
integer :: jatransportmodule = 1 !< use transport module (1) or subroutine (0)
integer :: itstep !< time step 0=no, 1 =step_explicit, 2=step_reduce, 3=step_jacobi, 4: explicit
integer :: iadvec !< adv type, 0=no, 1 = Wenneker vol, qu-udzt array, 2=1, function,
!< 3 =Perot in uit, 4 =Perot in, explicit
!< 5 =Perot in uit, 6 =Perot in, piaczek teta
!< 7 =Perot in uit, 8 =Perot in, piaczek fully implicit
!< 9 =Perot in uit, 10=Perot in, piaczek fully implicit, weir + factor
!< 11=Perot in uit, 12=Perot in, piaczek fully implicit, gate + factor
!< 20=Energy conserving compact, piaczek fully implicit, weir
integer :: iadvec1D !< same, now for 1D links
integer :: lincontin !< 0 = no, 1 = yes linear continuity
integer :: iPerot !< Perot weigthing type of cell center velocities ucx, ucy
!! in vectoren:
!! 0 : uc*sum(w) = sum (u W)
!! 1 : uc*A = sum(u dxa W)
!! 2 : uc*A*hs = sum(u dxa W hu ), ie waterdepth dependent
!! 2 : uc*V = sum(q dxa ), ie waterdepth dependent
!! 3 : uc*A*humx = sum(u dxa W hu ), humx = max(hu)
!! 4 : uc*A*humx = sum(u dxa W hu ), humx = max(hu)
integer :: icorio !< Coriolis weigthing
!! 1 : v = tang. comp of total velocity at u-point
!! 2 : v = idem, depth weighted
!! 3 : v = idem, depumin weighted
!! 4 : v = idem, continuity weighted ucnxq, ucnyq = default so far
!! 5 : 4, scaling down coriolis force below trshcorio
double precision :: trshcorio !< below this depth coriolis force scaled down linearly to 0
integer :: jatidep !< use tide potential forcing yes no
double precision :: doodsonstart, doodsonstop , doodsoneps
integer :: jatrt !< Trtrou = #Y# --> 1 , Trtrou = #N# --> 0 (Delf3D style input)
integer :: jasal !< Include salinity set in mdf
integer :: jasalsrc !< Need to use salinity source array or not
integer :: jatem !< Temperature model (0=no, 5=heatfluxmodel)
integer :: jaorgsethu = 1 !< 1 = before everything, as original
integer :: jarhoxu !< rho effects in momentum, 0=no, 1=in horizontal adv, 2=+ in vertical adv, 3 = + in pressure term
integer :: jawave !< Include wave model nr, 0=no, 1=fetchlimited hurdle stive + swart, 4: XBeach wave driver
integer :: jaavgwavquant = 0 !< 0=no time avg'd output for jawave==4; 1=time avg'd output for jawave==4. Avg'ing period = ti_avgwav
integer :: ihorvic !< 0=no visc, 1=do visc
integer :: jacreep !< Include anti-creep calculation, (0=no, 1=yes)
integer :: jasecflow !< 0: no, 1 : yes
integer :: jasecf !< variable same as jasecflow for kmx<2 , 0 for kmax >=2
integer :: javiusp !< if 1 spatially varying horizontal eddy viscosity
integer :: jadiusp !< if 1 spatially varying horizontal eddy viscosity
integer :: jaCdwusp !< if 1 spatially varying windstress coefficient
integer :: javiuplus3D = 1 !< add vertical eddy viscosity to horizontal eddy viscosity (1 = yes, 0 = no)
integer :: jafrculin !< use linear friction yes/no
integer :: iuvfield !< intialise this velocityfield: 0 = no
!! 1:u=y**2, 2:idem, 60 deg, 3:rotation, 4=lin, 5=lin 60 deg
integer :: istresstyp !< 1 : full stress tensor, semi link oriented horvic2
!! 2 : full stress tensor, fully link oriented dvxc = ok and fast
!! 3 : 2, volume weighted
!! 4 : full node oriented
!! 5 : 4, volume weighted
integer :: irov !< 0 : free slip
!! 1 : partial slip
!! 2 : no slip
!! 3 : glass (mind you, in D3DFLOW 3 means free slip)
integer :: ibedlevmode !< 1 : See BLMODE_DFM
!! 2 : See BLMODE_D3D
integer, parameter :: BLMODE_DFM = 1 !< ibedlevmode value: Compute bed levels solely by ibedlevtyp, i.e., derived from velocity points (or direct bl tiles)
integer, parameter :: BLMODE_D3D = 2 !< ibedlevmode value: Compute bed levels as D3D, i.e., derived from corner point depths. Currently always deepest (== DPSOPT=MAX).
integer :: ibedlevtyp !< 1 : Bottom levels at waterlevel cells (=flow nodes), like tiles xz, yz, bl , bob = max(bl left, bl right)
!! 2 : Bottom levels at velocity points (=flow links), xu, yu, blu, bob = blu, bl = lowest connected link
!! 3 : Bottom levels at velocity points (=flow links), using mean network levels xk, yk, zk bl = lowest connected link
!! 4 : Bottom levels at velocity points (=flow links), using min network levels xk, yk, zk bl = lowest connected link
integer :: ibedlevtyp1D !< 1 : same, 1D, 1 = tiles, xz(flow)=zk(net), bob(1,2) = max(zkr,zkl) , 3=mean netnode based
integer :: izbndpos !< 0 : waterlevel boundary location as in D3DFLOW, 1=on network boundary, 2=on specified boundary polyline
double precision :: blmeanbelow !< : if not -999d0, below this level the cell centre bedlevel is the mean of surrouding netnodes
double precision :: blminabove !< : if not -999d0, above this level the cell centre bedlevel is the min of surrouding netnodes
integer :: nonlin !< 1 : non-linear continuity eq
integer :: nonlin2D !< 1 : non-linear continuity eq for 2D points, only for jaconveyance2D >= 3 .and.
integer :: jaembed1D !< 1 : use embedded 1d channels, (first run 'Operations: Embed 1D channels'), 2=idem, non-linear cont.
integer :: numembed !< : nr of embedded (1D2D) cells
integer, allocatable :: kembed(:) !< : Embedded 2D flowcell nrs in ndx
double precision, allocatable :: dbl21(:) !< ! bl2D - bl1D
integer :: iproftypuni !< 1 : circle, 2 : rectan
double precision :: slotw2D !< minimum slotwidth 2D
double precision :: slotw1D !< minimum slotwidth 1D
double precision :: ditchh = 0.3 !< ditch depth
double precision :: ditchwfrac = 0d0 !< ditch width percentage
integer :: jaconveyance2D !< 1 : yes, 0 : no
integer :: nums1it !< : nr of non-linear continuity iterations
double precision :: dnums1it !< : total nr of non-linear continuity iterations
integer :: isimplefixedweirs !< 1=only links stored, 0=complete crossection paths stored
! Secondary Flow
! integer :: jasftesting !< (secondary flow testing: 0: no just compute fm velocitie, 1: prescribe exact ucx,ucy point values 2: prescribe exact edge velocities (FM reconstructs ucx,ucy) 3: prescribe exact ucx,ucy cell-averaged values
! integer :: jasfinttype !< 1: perot reconstruction 2: standard gauss
! integer :: jasftesttype !< -1: u_theta = 0.1*200/r; u_r = 0 | k: k=0..3, u_x = x^k, u_y = 0
double precision :: Uniformhu !< Uniformhu for arjen's membranes
double precision :: bedslope !< bottom inclination testcases
double precision :: Slopedrop2D !< Apply losses for 'rain from the roof', only if local bottom slope > Slopedrop2D, only for Slopedrop2D > 0.0
double precision :: drop3D !< Apply losses in or 3D if downwind z below bob + 2/3 hu
integer :: jaslopedr !< Apply this yes/no 0/1
double precision :: cflmx !< max Courant nr ()
double precision :: cflw !< wave velocity fraction, total courant vel = u + cflw*wavevelocity
double precision :: teta0 !< 1.00d0 ! .52 ! uniform teta in horizontal (),
integer :: ivariableteta !< 0=fully implicit, 1=teta constant, 2=variable teta
!! (set teta=1.0) (set teta=0.51->0.99) (set teta<0)
integer :: jacstbnd !< Delft-3D type cell-centered velocities at boundaries (ucx, ucy)
integer :: limtypsa !< 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC, 21=central voor stof transport
integer :: limtypTM !< 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC, 21=central voor stof transport
integer :: limtypsed !< 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC, 21=central voor stof transport
integer :: limtyphu !< 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC, 21=central voor hu
integer :: limtypmom !< 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC, 21=central voor momentum transport
integer :: ifixedweirscheme !< 0 = no, 1 = compact stencil, 2 = whole tile lifted, full subgrid weir + factor
integer :: ifxedweirfrictscheme !< 0 = friction based on hu, 1 = friction based on subgrid weirfriction scheme
double precision :: fixedweircontraction !< flow width = flow width*fixedweircontraction
double precision :: fixedweirtopwidth !< , e.g. 4.00 (m)
double precision :: fixedweirtopfrictcoef !< if .ne. dmiss, use this friction coefficient on top width
double precision :: fixedweirtalud !< , e.g. 0.25 ( ) for 1 to 4 talud
double precision :: waquaweirthetaw=0.6d0 !< , e.g. 0.6
double precision :: sini !< uniform initial waterlevel (m), (uniform bottom level = zkuni)
double precision :: uini !< uniform initial velociy (m/s),
double precision :: salini !< uniform initial sal (ppt)
double precision :: deltasalinity=-999d0 !< uniform initial sal (ppt)
double precision :: Sal0abovezlev !< sal=0 above lev= zlev (m)
double precision :: temini !< uniform initial temp (degC)
double precision :: spirini !< uniform initial spirint (m/s)
double precision :: zbnd !< for now only, uniform waterlevel on boundary
double precision :: zkdropstep !< Amount of bottomlevel to be added with dropland (m)
double precision :: sdropstep !< Amount of water to be added with dropwater (m)
double precision :: eps4 !< min au in poshchk
double precision :: eps6 !<
double precision :: eps8 !< implicit diffusion
double precision :: eps10 !<
double precision :: epshsdif=1d-2 !< hs < epshsdif: no vertical diffusion if hs < epshsdif
double precision :: s01max !< water level threshold (m) between s0 and s1 in validation routine
double precision :: u01max !< velocity threshold (m/s) between u0 and u1 in validation routine
! See also m_flowtimes::dtminbreak
! parameters controlling flooding/drying/solving
double precision :: epshu !< minimum waterdepth for setting hu>0
double precision :: epshs !< minimum waterdepth for setting cfu
double precision :: epswav !< minimum waterdepth for wave calculations
double precision :: epswavxbeach=0.2d0 !< thrhehshhhohlhdh waterdepth for wave calculations
double precision :: chkhuexpl !< only for step_explicit: check computed flux beneath this waterdepth
double precision :: chkadvd !< check advection for 'drying' below this (upwind) waterdepth
integer :: jposhchk !< check for positive waterdepth; 0 = no
!! 1 = 0.7*dts, just redo
!! 2 = 1.0*dts, close all links
!! 3 = 0.7*dts, close all links
!! 4 = 1.0*dts, reduce au
!! 5 = 0.7*dts, reduce au
integer :: jsolpos !< in iterative solver force solution above bottom level
integer :: Icgsolver !< 'Solver type , 1 = sobekGS_OMP, 2 = sobekGS_OMPthreadsafe, 3 = sobekGS, 4 = sobekGS + Saadilud, 5 = parallel/global Saad, 6 = parallel/Petsc, 7 = parallel/GS '
integer :: ipre !< Preconditioner, 0=rowscaling, 1=GS, 2=trial
integer :: mdump ! dump file unit nr
double precision :: hwetbed !< for case wetbed
integer :: javau !< vert. adv. u1 : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL
integer :: javakeps !< vert. adv. keps : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL
integer :: javasal !< vert. adv. sa1 : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL, switched to 3 for neg. strat.
integer :: javatem !< vert. adv. tem1 : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL, switched to 3 for neg. strat.
integer :: javased !< vert. adv. suspended sediment concentrations : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL, 5=switched to 3 for neg stratif., 6=higher-order upwind/explicit
integer :: javatest !< vert. adv. keps : test, 0 = no
integer :: jahazlayer !< vertical treatment of horizontal advection in z layers 1=org, 2=sigma, 3=node volumes
integer :: jabaroctimeint !< time integration baroclini pressure, 1 = Euler, abs() = 2; rho at n+1/2, 3: AdamsB
integer :: jalts = 0 ! local time-stepping (1) or not (0)
integer :: jaanalytic !< analytic solution available in black sideview => do not also show computed surface in black
integer :: jaustarint !< 1=integral bed layer velocity, 0=velocity at half bed layer
integer :: jaqaisq1 = 0 !< 1 : qa = q1, 0 : qa = au*u1
integer :: inisal2D !< 1 if specified through meteo module
integer :: initem2D !< 1 if specified through meteo module
double precision :: cffacver = 0d0 !< switch to low order at high cf in constituent transport vertical, 1d0=yes, 0d0 = no
double precision :: toplayminthick !< minimum top layer thickness (m)
double precision :: botlayminthick !< minimum bot layer thickness (m)
double precision :: s1incinterval !< incremental interval. if > 0 write incremental
integer :: jbasqbnddownwindhs !< 0 : original hu on qbnd, 1 = downwind hs on qbnd
integer :: maxitverticalforestersal !< 100, max iterations vertical forester
integer :: maxitverticalforestertem !< 100, max iterations vertical forester
double precision :: salmax !< filter if sal > maxsal
! written to his file yes or no
integer :: jahisbal !< Write mass balance/volume totals to his file, 0: no, 1: yes
! written to map file yes or no
integer :: jamaps0 !< previous step water levels to map file, 0: no, 1: yes
integer :: jamaps1 !< water levels to map file, 0: no, 1: yes
integer :: jamapu1 !< velocities to map file, 0: no, 1: yes
integer :: jamapu0 !< previous step velocities to map file, 0: no, 1: yes
integer :: jamapucvec !< velocity vectors to map file, 0: no, 1: yes
integer :: jamapww1 !< upward velocity on flow link to map file, 0: no, 1: yes
integer :: jamapnumlimdt !< num limdt to map file, 0: no, 1: yes
integer :: jamaptaucurrent !< shear stress to map file, 0: no, 1: yes
integer :: jamapchezy !< chezy to map file, 0: no, 1: yes
integer :: jamapsal !< salinity to map file, 0: no, 1: yes
integer :: jamaptem !< temperature to map file, 0: no, 1: yes
integer :: jamapconst !< constituents to map file, 0: no, 1: yes
integer :: jamapsed !< sediment fractions to map file, 0: no, 1: yes
integer :: jamaptur !< k, eps and vicww to map file, 0: no, 1: yes
integer :: jamaptrachy !< trachytope roughnesses to map file, 0: no, 1: yes
integer :: jamapwind !< wind velocities to map file, 0: no, 1: yes
integer :: jamapviu !< horizontal viscosity to map file, 0: no, 1: yes
integer :: jamapdiu !< horizontal diffusity to map file, 0: no, 1: yes
integer :: jamaprho !< flow density to map file, 0: no, 1: yes
integer :: jamapq1 !< flow flux to map file, 0: no, 1: yes
integer :: jamapspir !< spiral flow to map file, 0: no, 1: yes
integer :: jafullgridoutput !< 0:compact, 1:full time-varying grid data
integer :: jaeulervel !< 0:GLM, 1:Euler velocities
! parameters for parms solver
integer, parameter :: NPARMS_INT=2 !< for parms solver, number of integer parameters
integer, parameter :: IPARMS_ILUTYPE=1
integer, parameter :: IPARMS_NLEVEL=2
character(len=128), dimension(NPARMS_INT), parameter :: iparmsnam= [ character(len=128):: 'ilutype', 'nlevel' ]
integer, dimension(NPARMS_INT) :: iparms
integer, parameter :: NPARMS_DBL=1 !< for parms solver, number of double precision parameters
integer, parameter :: IPARMS_DTOL=1
character(len=128), dimension(NPARMS_DBL), parameter :: dparmsnam= [ character(len=128):: 'dtol' ]
double precision, dimension(NPARMS_DBL) :: dparms
contains
!> Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, only call reset_flowparameters() instead.
subroutine default_flowparameters()
jatransportmodule = 1 ! transportmethod 1 (module) or 0 hk (subroutine)
itstep = 2 ! time step 0=only transport, 1=transport + velocity update, 2=full implicit step_reduce
iadvec = 33 ! adv type, 0=no, 1= Wenneker vol, qu-udzt array, 2=1, function, 3=Perot in uit, 4=Perot in, 5=3,piaczek
iadvec1D = 33 ! same, now for 1D links
lincontin= 0 ! 0 = no, 1 = yes linear continuity
iPerot = 1 ! Perot weigthing type of cell center velocities ucx, ucy
! in vectoren:
! 0 : uc*sum(w) = sum (u W)
! 1 : uc*A = sum(u dxa W)
! 2 : uc*A*hs = sum(u dxa W hu ), ie waterdepth dependent
! 2 : uc*V = sum(q dxa ), ie waterdepth dependent
! 3 : uc*A*humx = sum(u dxa W hu ), humx = max(hu)
! 4 : uc*A*humx = sum(u dxa W hu ), humx = max(hu)
! 5 : uc*Vc = sum(u dxa W hu ), Vc = dxa W hu based volume in cell
! 6 : as 5, also for Coriolis
icorio = 5 ! Coriolis weigthing
! 1 : v = tang. comp of total velocity at u-point
! 2 : v = idem, depth weighted
! 3 : v = idem, depumin weighted
! 4 : v = idem, continuity weighted ucnxq, ucnyq = default so far
! 5 : 4, scaling down coriolis force below trshcorio
trshcorio = 1.0 ! below this depth coriolis force scaled down linearly to 0
jatidep = 1 ! use tide potential forcing yes no
! DOODSONSTART = 55.565D0 ; DOODSONSTOP = 375.575D0 ; Doodsoneps = .00D0 ! standaard triwaq alle 484 cmp
DOODSONSTART = 55.565D0 ; DOODSONSTOP = 375.575D0 ; Doodsoneps = .03D0 ! standaard triwaq 60 cmp
! DOODSONSTART = 57.555D0 ; DOODSONSTOP = 275.555D0 ; Doodsoneps = .03D0 ! Delft3d
jasecflow = 0 ! include secondary flow (0=no, 1=yes)
jasecf = 0
! Anti-creep
jacreep = 0 ! Include anti-creep calculation, (0=no, 1=yes)
jasal = 0 ! Include salinity (autoset by flow_initexternalforcings())
jatem = 0 ! Temperature model
jarhoxu = 0 ! rho effects in momentum, 0=no, 1=in horizontal adv, 2=+ in vertical adv, 3 = + in pressure term
jased = 0 ! Include sediment
jatrt = 0 !< Include alluvial and vegetation roughness (trachytopes)
jawave = 0 ! Include wave model nr
jafrculin = 0 !< do not use linear friction
javiusp = 0 !< spatially varying eddyviscosity yes/no 1/0
jadiusp = 0 !< spatially varying eddydiffusivity yes/no 1/0
ihorvic = 0 !< 0=no visc, 1=do visc
iuvfield = 0 ! intialise this velocityfield: 0 = no
! 1:u=y**2, 2:idem, 60 deg, 3:rotation, 4=lin, 5=lin 60 deg
istresstyp = 3 ! 1 : full stress tensor, semi link oriented horvic2
! 2 : full stress tensor, fully link oriented dvxc = ok and fast
! 3 : 2, volume weighted
! 4 : full node oriented
! 5 : 4, volume weighted
irov = 0 ! 0 : free slip
! 1 : partial slip
! 2 : no slip
! 3 : glass (mind you, in D3DFLOW 3 means free slip)
ibedlevmode = BLMODE_DFM !< Default: Compute bed levels solely by ibedlevtyp, i.e., derived from velocity points (or direct bl tiles).
ibedlevtyp = 3 ! 1 : Bottom levels at waterlevel cells (=flow nodes), like tiles xz, yz, bl , bob = max(bl left, bl right)
! 2 : Bottom levels at velocity points (=flow links), xu, yu, blu, bob = blu, bl = lowest connected link
! 3 : Bottom levels at velocity points (=flow links), using mean network levels xk, yk, zk bl = lowest connected link
! 4 : Bottom levels at velocity points (=flow links), using min network levels xk, yk, zk bl = lowest connected link
! 5 : Bottom levels at velocity points (=flow links), using max network levels xk, yk, zk bl = lowest connected link
blmeanbelow = -999d0
blminabove = -999d0
ibedlevtyp1D = 3 !< 1 : same, 1D, 1 = tiles, xz(flow)=zk(net), bob(1,2) = max(zkr,zkl) , 3=mean netnode based
izbndpos = 0 !< 0 : waterlevel boundary location as in D3DFLOW, 1=on network boundary, 2=on specified boundary polyline
nonlin = 0 !< 1 : non-linear continuity eq, now governed by iproftyp in volsur, 0 for rectan, else 1
nonlin2D = 0 !< 1 : non-linear continuity eq in 2D, sets nonlin
jaembed1D = 0 !< 1 : use embedded 1d channels
iproftypuni = 3 ! 1 : circle, 2 : rectan R=A/P , 3 = rectan, R=H, 4 = 3,nonlin
slotw2D = 0d-3
slotw1D = 1d-3
jafullgridoutput = 0 !output grid in compact manner or full manner
jaeulervel = 0 !GLM velocities
jaconveyance2D = 1 !
bedslope = 0d0 ! bottom inclination testcases
Slopedrop2D = 0d0 ! Apply droplosses only if local bottom slope > Slopedrop2D, negative = no droplosses
drop3D = -999d0 ! Apply droplosses in 3D yes or no 1 or 0
jaslopedr = 0 ! Flag only
jacstbnd = 0
cflmx = 0.7d0 ! max Courant nr ()
cflw = 0.1d0 ! wave velocity fraction, total courant vel = u + cflw*wavevelocity
teta0 = 0.55d0 ! 1.00d0 ! .52 ! uniform teta in horizontal (),
ivariableteta = 0 ! 0=fully implicit, 1=teta constant, 2=variable teta
! (set teta=1.0) (set teta=0.51->0.99) (set teta<0)
limtypsa = 4 ! 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC voor scalar transport SALINITY
limtypTM = 4 ! 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC voor scalar transport TEMPERATURE
limtypsed = 4 ! 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC voor scalar transport SEDIMENT
limtyphu = 0 ! 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC voor hu WATERHEIGHT AT U POINT
limtypmom = 4 ! 0=no, 1=minmod, 2=vanleer, 3=koren 4=MC voor MOMENTUM transport
ifixedweirscheme = 6 !< 0 = no special treatment, setbobs only, 1 = compact stencil, 2 = whole tile lifted, full subgrid weir + factor
fixedweircontraction = 1d0 !< flow width = flow width*fixedweircontraction
fixedweirtopwidth = 3d0 !< e.g. 4.00 (m)
fixedweirtopfrictcoef = -999d0 !< if .ne. dmiss, use this friction coefficient on top width
fixedweirtalud = 0.25d0 !< e.g. 0.25 ( ) for 1 to 4 talud
isimplefixedweirs = 1
ifxedweirfrictscheme = 1 !< 0 = friction based on hu, 1 = friction based on subgrid weirfriction scheme
sini = 0d0 ! uniform initial waterlevel (m), (uniform bottom level = zkuni)
uini = 0 ! uniform initial velocity (m/s)
salini = 0d0 !< uniform initial sal (ppt)
temini = 6d0 !< uniform initial temp (degC)
spirini = 0d0 !< uniform initial spirint (m/s)
Sal0abovezlev = -999d0 !< if not dmiss, only seta salini below this level
zkdropstep = 1d0 !< Amount of bottomlevel to be added with dropland (m)
sdropstep = 1d0 !< Amount of water to be added with dropwater (m)
uniformhu = -999d0 !< Uniformhu
zbnd = 2d0 ! for now only, uniform waterlevel on boundary
eps4 = 1d-4 ! min au in poshchk
eps6 = 1d-6 !
eps8 = 1d-8 ! implicit diffusion
eps10 = 1d-10 !
s01max = 0d0 ! max. water level change: off
u01max = 0d0 ! max. velocity change: off
! See also: m_flowtimes::dtminbreak
! parameters controlling flooding/drying/solving
epshu = 1d-4 ! minimum waterdepth for setting hu>0
epshs = .2d0*epshu ! minimum waterdepth for setting cfu
epswav = 1d-2 ! minimum waterdepth for wave calculations
chkhuexpl = 0.1d0 ! only for step_explicit: check computed flux beneath this waterdepth
chkadvd = 0.1d0 ! check advection for 'drying' below this (upwind) waterdepth
jposhchk = 2 ! check for positive waterdepth; 0 = no
! 1 = 0.7*dts, just redo
! 2 = 1.0*dts, close all links
! 3 = 0.7*dts, close all links
! 4 = 1.0*dts, reduce au
! 5 = 0.7*dts, reduce au
jsolpos = 0 ! in iterative solver force solution above bottom level
Icgsolver = 4 ! Icgsolver = 1 ! 1 = GS_OMP, 2 = GS_OMPthreadsafe, 3 = GS, 4 = Saadilud
ipre = 0 ! preconditioner, 0=rowscaling, 1=GS, 2=trial
hwetbed = 0.2d0 ! for case wetbed
javau = 3 !< vert. adv. u1 : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL
javakeps = 3 !< vert. adv. keps : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL
javasal = 6 !< vert. adv. sa1 : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL, 5=switched to 3 for neg stratif, 6=hoexplicit.
javatem = 6 !< vert. adv. tem1 : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL, 5=switched to 3 for neg stratif.
javased = 6 !< vert. adv. suspended sediment concentrations : 0=No, 1=UpwexpL, 2=Centralexpl, 3=UpwimpL, 4=CentraLimpL, 5=switched to 3 for neg stratif., 6=higher-order upwind/explicit
jahazlayer = 0 !<
jabaroctimeint = -2 !< time integration baroclini pressure, 1 = explicit, abs() = 2; adams bashford , 3 = ab3, 5 = adv rho
jaanalytic = 0 !< analytic solution available in black sideview => do not also show computed surface in black
jaustarint = 1 !< 1=integral bed layer velocity, 0=velocity at half bed layer
inisal2D = 0 !< 1 if specified through meteo module
initem2D = 0 !< 1 if specified through meteo module
toplayminthick = 0d0 ! minimum top layer thickness (m)
s1incinterval=0 !< incremental interval. if > 0 write incremental
jbasqbnddownwindhs = 0 !< 0 : original hu on qbnd, 1 = downwind hs on qbnd
maxitverticalforestersal = 100 !< 100, max iterations vertical forester
maxitverticalforestertem = 0 !< 100, max iterations vertical forester
salmax = 0d0 !< filter if sal > maxsal
! Remaining of variables is handled in reset_flowparameters()
! call reset_flowparameters()
iparms = 0 ! parms-default
dparms = 0d0 ! parms-default
jahisbal = 1
jamaps0 = 1
jamaps1 = 1
jamapu0 = 1
jamapu1 = 1
jamapucvec = 1
jamapww1 = 1
jamapnumlimdt = 1
jamaptaucurrent = 1
jamapchezy = 1
jamapsal = 1
jamaptem = 1
jamapconst = 1
jamapsed = 1
jamaptur = 1
jamaptrachy = 1
jamapwind = 1
jamapviu = 1
jamapdiu = 1
jamaprho = 1
jamapq1 = 1
jamapspir = 1
call reset_flowparameters()
end subroutine default_flowparameters
!> Resets only flowparameters variables intended for a restart of an existing flow simulation (same MDU).
!! Upon loading of new model/MDU, call default_flowparameters() instead.
subroutine reset_flowparameters()
end subroutine reset_flowparameters
end module m_flowparameters
module m_statistics
implicit none
double precision :: avedif !< for now only, cum dif with analytic sol
double precision :: sqadif !< for now only, cum dif with analytic sol
double precision :: rmsdif !< for now only, cum dif with analytic sol
double precision :: dmxdif !< for now only, cum dif with analytic sol
integer :: numdif
double precision :: cumavedif !< for now only, cum dif with analytic sol
double precision :: cumrmsdif !< for now only, cum dif with analytic sol
double precision :: cumdmxdif !< for now only, cum dif with analytic sol
integer :: numcum
contains
subroutine reset_statistics()
avedif = 0d0 ! for now only, cum dif with analytic sol
sqadif = 0d0
rmsdif = 0d0 ! for now only, cum dif with analytic sol
dmxdif = 0d0 ! for now only, cum dif with analytic sol
numdif = 0
cumavedif = 0d0 ! for now only, cum dif with analytic sol
cumrmsdif = 0d0 ! for now only, cum dif with analytic sol
cumdmxdif = 0d0 ! for now only, cum dif with analytic sol
numcum = 0
end subroutine reset_statistics
end module m_statistics
module m_missing
implicit none
double precision :: dmiss = -999d0 !
double precision :: xymis = -999d0 !
double precision :: dxymis = -999d0
integer :: LMOD, KMOD ! TBV READDY
integer :: jins = 1
end module m_missing
module m_flow ! flow arrays-999
use m_flowparameters
use m_flowexternalforcings
use m_physcoef
use m_turbulence
use m_grw
! use m_fixedweirs
use m_heatfluxes
use m_alloc
implicit none
! 3D parameters
integer :: kmx !< nr of 3d layers, increasing in positive upward direction
!! if kmx==0 then 2D code. if kmx==1 then 3D code
integer :: kmx1 !< kmx + 1, for dimensioning arrays that used be (0:kmax)
integer :: kmxd !< dim of kmx, >= 1
integer :: ndkx !< dim of 3d flow nodes (internal + boundary)
integer :: ndkx1 !< dim of 3d flow horizontal interfaces (internal + boundary), (0:kmx)
integer :: lnkx !< dim of 3d flow links (internal + boundary)
integer :: numvertdis !< number of vertical layer distributions
integer :: mxlayz !< max nr of z layers in flow domain
integer :: mxlays !< max nr of sigma layers in flow domain
integer :: kplot !< layer nr to be plotted
integer :: nplot !< vertical profile to be plotted at node nr
integer :: layertype !< 1= all sigma, 2 = all z, 3 = left sigma, 4 = left z
integer :: numtopsig = 0 !< number of top layers in sigma
double precision :: Tsigma = 100 !< relaxation period density controlled sigma
integer, parameter :: LAYTP_SIGMA = 1
integer, parameter :: LAYTP_Z = 2
integer, parameter :: LAYTP_LEFTSIGMA = 3
integer, parameter :: LAYTP_LEFTZ = 4
integer :: iStrchType = 0
integer, parameter :: STRCH_USER = 1
integer, parameter :: STRCH_EXPONENT = 2
integer :: iturbulencemodel !< 0=no, 1 = constant, 2 = algebraic, 3 = k-eps
integer :: ieps !< bottom boundary type eps. eqation, 1=dpmorg, 2 = dpmsandpit, 3=D3D, 4=Dirichlethdzb
double precision :: sigmagrowthfactor ! Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, only call reset_flow() instead.
subroutine default_flow()
! 3D parameters
! kmx = 0 ! nr of 3d layers, increasing in positive upward direction
! if kmx==0 then 2D code. if kmx==1 then 3D code
! kmx1 = kmx+1 ! kmx + 1, for dimensioning arrays that used be (0:kmax)
! kmxd = 0 ! dim of kmx, >= 1
! ndkx = 0 ! dim of 3d flow nodes (internal + boundary)
! ndkx1 = 1 ! dim of 3d flow horizontal interfaces (internal + boundary), (0:kmx)
! lnkx = 0 ! dim of 3d flow links (internal + boundary)
mxlayz = 1 ! max nr of z layers in flow domain
mxlays = 1 ! max nr of sigma layers in flow domain
kplot = 1 ! layer nr to be plotted
nplot = 1 ! vertical profile to be plotted at node nr
layertype = 1 !< 1= all sigma, 2 = all z, 3 = left sigma, 4 = left z
iturbulencemodel = 3 !< 0=no, 1 = constant, 2 = algebraic, 3 = k-eps, 4 = k-tau
ieps = 2 !< bottom boundary type eps. eqation, 1=dpmorg, 2 = dpmsandpit, 3=D3D, 4=Dirichlethdzb
sigmagrowthfactor = 1d0 ! Resets only flow variables intended for a restart of flow simulation.
!! Upon loading of new model/MDU, call default_flow() instead.
subroutine reset_flow()
! node related
! basis
vol0tot = 0 ! total volume start of timestep (m3)
vol1tot = 0 ! total volume end of timestep (m3)
vol1ini = -1d0 ! total volume initially (m3)
qinbnd = 0 ! total inflow boundaries (m3/s) Actual values
qoutbnd = 0 ! total outflow boundaries (m3/s)
qincel = 0 ! total inflow cells (m3/s) Actual values
qoutcel = 0 ! total outflow cells (m3/s)
vinbnd = 0 ! total volume in boundaries (m3) Actual values
voutbnd = 0 ! total volume out boundaries (m3)
vincel = 0 ! total volume in cells (m3) Actual values
voutcel = 0 ! total volume out cells (m3)
volerr = 0 ! vol1tot - vol0tot - vinbnd + voutbnd - vincel + voutcel (m3)
vinbndcum = 0 ! total inflow boundaries (m3) Cumulative values
voutbndcum = 0 ! total outflow boundaries (m3)
vincelcum = 0 ! total volume in cells (m3) Actual values
voutcelcum = 0 ! total volume out cells (m3)
volerrcum = 0 ! (m3)
dvolbot = 0d0 ! (m3)
! extra
qinrain = 0 ! total inflow rain (m3/s)
qouteva = 0 ! total outflow evaporation (m3/s)
qinlat = 0 ! total inflow diffuse laterals (m3/s)
qoutlat = 0 ! total outflow diffuse laterals (m3/s)
qinsoil = 0 ! total inflow groundwater (m3/s)
qoutsoil = 0 ! total outflow groundwater (m3/s)
qinsrc = 0 ! total inflow local point sources (m3/s)
qoutsrc = 0 ! total outflow local pount sources (m3/s)
vinrain = 0 ! total volume in rain (m3)
vouteva = 0 ! total volume out evaporation (m3)
vinlat = 0 ! total volume in diffuse laterals (m3)
voutlat = 0 ! total volume out diffuse laterals (m3)
vinsoil = 0 ! total volume in groundwater (m3)
voutsoil = 0 ! total volume out groundwater (m3)
vinsrc = 0 ! total volume in local point sources (m3)
voutsrc = 0 ! total volume out local pount sources (m3)
vinraincum = 0 ! total inflow rain (m3)
voutevacum = 0 ! total outflow evaporation (m3)
vinlatcum = 0 ! total inflow diffuse laterals (m3)
voutlatcum = 0 ! total outflow diffuse laterals (m3)
vinsoilcum = 0 ! total inflow groundwater (m3)
voutsoilcum = 0 ! total outflow groundwater (m3)
vinsrccum = 0 ! total inflow local point sources (m3)
voutsrccum = 0 ! total outflow local pount sources (m3)
a0tot = 0 ! total wet surface area start of timestep (m2)
a1tot = 0 ! total wet surface area end of timestep (m2)
a1ini = 0 ! total wet surface area initially (m2)
ek1tot = 0 ! volume averaged kin energy (m2/s2) end of timestep
ep1tot = 0 ! volume averaged pot energy (m2/s2) end of timestep
ep1rela = 0 ! time relaxe ep1tot
hsaver = 0 ! average waterdepth (m), vol/are
epsmaxvol = 1d-9 ! eps vol diff (m3) ! both not used now
difmaxlev = 0 ! max lev diff (m3)
epsmaxlev = 1d-8 ! eps lev diff (m3)
debugon = .false. ! texts yes or no
validateon = .true. ! validate flow state yes or no
nodneg = 0 ! node nr with negative hs
kkcflmx = 0 ! 2D node nr with max courant
itsol = 0 ! act nr. of iterations in solve
nochkadv = 0 ! nr of chkadvd checks
numnodneg = 0 ! nr of posh checks
nrimptran = 0 ! nr of implicit transport points
ndmin = 0 ! node nr where min znod is found in viewing area
ndmax = 0 ! node nr where max znod is found in viewing area
Lnmin = 0 ! link nr where min zlin is found in viewing area
Lnmax = 0 ! link nr where max zlin is found in viewing area
sam0tot = 0d0 !< Total mass start of timestep (m3ppt)
sam1tot = 0d0 !< Total mass end of timestep (m3ppt)
samerr = 0d0 !< vol1tot - vol0tot - vinbnd + voutbnd - vincel + voutcel (m3)
voltot(:) = 0d0
voltotname(IDX_VOLTOT ) = 'total_volume'
voltotname(IDX_STOR ) = 'storage'
voltotname(IDX_VOLERR ) = 'volume_error'
voltotname(IDX_BNDIN ) = 'boundaries_in'
voltotname(IDX_BNDOUT ) = 'boundaries_out'
voltotname(IDX_BNDTOT ) = 'boundaries_total'
voltotname(IDX_EXCHIN ) = 'exchange_with_1D_in'
voltotname(IDX_EXCHOUT) = 'exchange_with_1D_out'
voltotname(IDX_EXCHTOT) = 'exchange_with_1D_total'
voltotname(IDX_PARTIP ) = 'precipitation'
voltotname(IDX_SOUR ) = 'source_sink'
end subroutine reset_flow
end module m_flow
module m_profiles
! profile related :
type tprof !< this is a profile type
integer :: ityp !< 1 = circle, 2=rectan1dlumped, 3=rectan2d, 9=yzlumped, 10=yzconveyance
integer :: frctp !< friction type chezy manning etc
double precision :: frccf !< friction coefficient
double precision :: width !< max width
double precision :: height !< max height
real, allocatable :: xx(:), yy(:) !< y z point coordinates
real, allocatable :: y (:), z (:) !< y and z arrays
end type tprof
integer :: nprofdefs !< nr of unique profile definitions
type(tprof), allocatable :: profiles1D(:) !< these are the profiles
integer :: nproflocs !< nr of profile locations, always <= nprofdefs
integer :: maxproflocnr !< highest referred profnr in profloc.xyz
integer :: minproflocnr !< lowest referred profnr in profloc.xyz
double precision, allocatable :: xpr(:), ypr(:), zpr(:) !< profile locations, x,y,z
integer, allocatable :: npr(:) !< at these locations, reference to profdefs
end module m_profiles
!> in m_flowgeom: nd and ln apply to waterlevel nodes and links
!! in m_netw : nod and lin apply to 'grid' or 'net' nodes and links
module m_flowgeom
use m_profiles
use grid_dimens_module
use m_d3ddimens
use m_flowparameters, only: jawave
implicit none
! node (s) related : dim=ndx
type tnode !< node administration
integer :: lnx !< max nr of links attached to this node
integer, allocatable :: ln (:) !< linknrs attached to this node, >0: to this flownode, <0: from this flownode
integer, allocatable :: nod(:) !< Mapping to net nodes
double precision, allocatable :: x (:) !< for now, this is only for quick/aligned plotting, the corners of a cell
double precision, allocatable :: y (:) !< for now, this is only for quick/aligned plotting, the corners of a cell
end type tnode
double precision :: bamin !< minimum 2D cell area
double precision :: bamin1D !< minimum cell area 1d nodes
double precision :: dxmin !< minimum link length (m)
double precision :: wu1DUNI !< uniform 1D profile width
double precision :: hh1DUNI !< uniform 1D profile height
type (griddimtype) :: griddim
type (gd_dimens), target :: gddimens
type (gd_dimens), pointer :: gddimens_ptr
! Flow node numbering:
! 1:ndx2D, ndx2D+1:ndxi, ndxi+1:ndx1Db, ndx1Db:ndx
! ^ 2D int ^ 1D int ^ 1D bnd ^ 2D bnd ^ total
integer :: ndx2D !< nr of 2d FLOW CELLS = NUMP
integer :: ndxi !< nr of internal flowcells (internal = 2D + 1D )
integer :: ndx1Db !< nr of flow nodes incl. 1D bnds (internal 2D+1D + 1D bnd)
integer :: ndx !< nr of flow nodes (internal + boundary)
type (tnode), allocatable :: nd(:) !< (ndx) node administration
integer, allocatable :: kcs(:) !< node code permanent
integer, allocatable, target :: kfs(:) !< [-] node code flooding {"shape": ["ndx"]}
integer, allocatable, target :: kfst0(:) !< [-] node code flooding {"shape": ["ndx"]}
double precision, allocatable, target :: ba (:) !< [m2] bottom area, if < 0 use table in node type {"shape": ["ndx"]}
double precision, allocatable :: bai(:) !< inv bottom area (m2), if < 0 use table in node type
double precision, allocatable, target :: bl(:) !< [m] bottom level (m) (positive upward) {"shape": ["ndx"]}
double precision, allocatable :: xz (:) !< waterlevel point / cell centre x (m)
double precision, allocatable :: yz (:) !< waterlevel point / cell centre y (m)
double precision, allocatable :: aif(:) !< cell based skewness ai factor sqrt(1+(dz/dy)**2) = abed/asurface
!< so that cfu=g(Au/conveyance)**2 = g*aif*(Au/convflat)**2
!< convflat is flat-bottom conveyance
double precision, allocatable :: aifu(:) !< bed skewness at u point (Lnx)
double precision, allocatable :: bz(:) !< averaged bed level at cell center (Ndx)
! link (u) related : dim = lnx
! Flow link numbering:
! 1:lnx1d, lnx1d+1:lnxi, lnxi+1:lnx1Db, lnx1Db+1:lnx
! ^ 1D int ^ 2D int ^ 1D bnd ^ 2D bnd ^ total
integer :: lnx1D !< nr of 1D flow links (so first 1D, next 2D, next boundaries)
integer :: lnxi !< nr of flow links (internal, 1D+2D )
integer :: lnx1Db !< nr of flow links including 1D bnds (internal, 1D+2D, boundary: only 1D. 2D bnd behind it)
integer :: lnx !< nr of flow links (internal + boundary). First we have 1D links, next 2D links, next boundary links (first 1D, then 2D)
integer, allocatable, target :: ln (:,:) !< [-] link (2,*) node administration, 1=nd1, 2=nd2 linker en rechter celnr {"shape": [2, "lnkx"]}
integer, allocatable, target :: lncn (:,:) !< [-] link (2,*) corner administration, 1=nod1, 2=nod2 linker en rechter netnr {"shape": [2, "lnkx"]}
integer, allocatable :: lnk (:,:) !< link (2,*) 3D cel administration, 1=ndk1, 2=ndk2 linker en rechter kcelnr
integer, allocatable :: kcu (:) !< link code, 1=1D link, 2=2D link, -1= bc 1D, -2=bc 2D, 3=2D parall wall, 4=1D2Dlink, 5=Pump
integer, allocatable, target :: iadv(:) !< [-] type of advection for this link {"shape": ["lnx"]}
double precision, allocatable :: teta (:) !< link teta (m)
integer, allocatable :: klnup (:,:) !< link upwind cell pointer if q> 0 use (1:3,L), else (4:6,L)
double precision, allocatable :: dx (:) !< link length (m)
double precision, allocatable :: dxi (:) !< inverse dx
double precision, allocatable :: wu (:) !< link initial width (m), if < 0 pointer to convtab
double precision, allocatable :: wui (:) !< inverse link initial width (m), if < 0 pointer to convtab
real, allocatable :: prof1D (:,:) !< dim = (3,lnx1D) 1= 1D prof width, 2=1D profile height, 3=proftyp, or: if 1,2< 0, pointers to prof 1,2, then 3=alfa1
double precision, allocatable, target :: bob (:,:) !< [m] left and right inside lowerside tube (binnenkant onderkant buis) HEIGHT values (m) (positive upward) {"shape": [2, "lnx"]}
integer, allocatable :: ibot (:) !< local ibedlevtype for setting min or max network depths (temporary, result goes to bobs)
double precision, allocatable :: acl ( :) !< left dx fraction, alfacl
double precision, allocatable :: acn (:,:) !< 2,L left and right wu fraction
double precision, allocatable :: xu (:) !< velocity point x (m)
double precision, allocatable :: yu (:) !< velocity point y (m)
double precision, allocatable :: blu (:) !< velocity point bottom level positive up (m)
double precision, allocatable :: csu (:) !< cosine comp of u0, u1
double precision, allocatable :: snu (:) !< sine comp of u0, u1
double precision, allocatable :: wcl (:,:) !< link weights (2,lnx) for center scalar , 1,L for k1, 2,L for k2
!double precision, allocatable :: wcxyl (:,:) !< link weights (4,lnx) for center velocities 1=1-x, 2=1-y, 3=2-x, 4=2y
!double precision, allocatable :: wcnxyl(:,:) !< link weights (4,lnx) for corner velocities 1=3-x, 2=3-y, 3=4-x, 4=4y
double precision, allocatable :: wcx1(:) !< link weights (lnx) for corner velocities k3
double precision, allocatable :: wcy1(:) !< link weights (lnx) for corner velocities k3
double precision, allocatable :: wcx2(:) !< link weights (lnx) for corner velocities k4
double precision, allocatable :: wcy2(:) !< link weights (lnx) for corner velocities k4
double precision, allocatable :: wcnx3(:) !< link weights (lnx) for corner velocities k3
double precision, allocatable :: wcny3(:) !< link weights (lnx) for corner velocities k3
double precision, allocatable :: wcnx4(:) !< link weights (lnx) for corner velocities k4
double precision, allocatable :: wcny4(:) !< link weights (lnx) for corner velocities k4
double precision, allocatable :: slnup (:,:) !< link upwind cell weight, if q> 0 use (1:3,L), else (4:6,L)
integer, allocatable :: ln2lne(:) !< flowlink to netlink nr dim = lnx
integer, allocatable :: lne2ln(:) !< netlink to flowlink nr dim = numL
! cell corner related, the links attached to a cell corner
type tcorn !< corner administration
integer :: lnx !< max nr of links attached to this corner
integer, allocatable :: ln (:) !< linknrs attached to this corner
integer :: nwx !< nr of walls attached
integer, allocatable :: nw(:) !< wallnrs attached to this corner
end type tcorn !< corner administration
type(tcorn) , allocatable :: cn (:) !< cell cornerpoints, (in counting order of nod)
double precision, allocatable :: ucnx(:) !< cell corner velocity, global x-dir (m/s)
double precision, allocatable :: ucny(:) !< cell corner velocity, global y-dir (m/s) (in m_flowgeom...)
double precision, allocatable, target :: vort(:) !< [s-1] vorticity at netnodes {"shape": ["ndx"], "comment": "Currently not available, is nowhere allocated nor filled."}
! fixed wall related, may be expanded to closed internal walls later for now, dim=(7,*)
integer :: mxwalls !< max nr of walls
double precision, allocatable :: walls(:,:) !< 1,* : inside waterlevel point (node)
!! 2,* : first cornerpoint
!! 3,* : second cornerpoint
!! 4,* : flow link 1 attached to first cornerpoint
!! 5,* : flow link 2 attached to second cornerpoint
!! 6,* : stress contribution to link 1
!! 7,* : stress contribution to link 1
integer :: nwcnx !< max nr of cornerpoints to which walls are attached
integer, allocatable :: nwalcnw(:,:) !< pointer to those walls, 1 = left wall, 2 =right wall
! closed wall corner (netnode) related
integer :: nrcnw !< nr of cn points attached to 2 closed walls
integer , allocatable :: kcnw (:) !< closed corner point nr k, reference to net nodes
real , allocatable :: cscnw (:) !< closed corner alignment cos (1:nrcnw)
real , allocatable :: sncnw (:) !< closed corner alignment sin (1:nrcnw)
real , allocatable :: sfcnw (:) !< closed corner partial slip sf = u*/u (based on walls average)
! branch related :
type tbranch !< this is a branch type
integer :: nx !< with nx links and nx + 1 nodes in it
integer, allocatable :: ln (:) !< successive flow linknrs
end type tbranch
integer :: mxflowbr !< max nr of flow branches
type(tbranch), allocatable :: flowbr(:) !< this is a list of flow branches
integer, allocatable :: Lbnd1D(:) !< for prof1D, boundary links refer to regular attached 1D links
! 1D endnode related
integer :: mx1Dend !< nr of 1D endnodes
integer, allocatable :: n1Dend(:) !< node nrs of 1D endnodes
! netnode/flownode related, dim = mxban
double precision, allocatable :: banf (:) !< horizontal netnode/flownode area (m2)
double precision, allocatable :: ban (:) !< horizontal netnode area (m2)
integer , allocatable :: nban (:,:) !< base area pointers to banf, 1,* = netnode number, 2,* = flow node number
integer :: mxban !< max dim of ban
! useful parameters :
double precision :: rrtol !< relative cellsize factor in search tolerance ()
double precision, allocatable :: xyen(:,:) !< temp boundary opposite point (end of EdgeNormal) (replaces ebtol tolerance)
integer :: jasnelucxy !< method for cell center velocities
!! 1 = slow node oriented, 2 = link csu etc, 3= fast linkweight
integer :: jarenumber !< renumberFlowNodes
integer :: jaFlowNetChanged !< To enforce various net(link)-related init routines after renumbering
! JRE Stuff related to setting up wave directional grid
integer :: ntheta !< Number of wave direction bins
double precision :: thetamax !< upper limit wave directional sector
double precision :: thetamin !< lower limit wave directional sector
double precision :: thetanaut !< nautical convention or not
double precision :: dtheta !< directional resolution
double precision :: theta0 !< mean theta-grid direction
double precision, allocatable :: thetabin(:) !< bin-means of theta-grid
! Debug parameter
integer :: jabatch = 0
contains
!> Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, call reset_flowgeom() instead.
subroutine default_flowgeom()
bamin = 1d-6 ! 1d0 ! minimum 2D cell area
bamin1D = 0d-2 ! minimum cell area 1d nodes
dxmin = 1D-3 ! minimum link length (m)
wu1DUNI = 2d0 ! Uniform 1D profile width
hh1DUNI = 3d0 ! Uniform 1D profile height
! useful parameters :
rrtol = 3d0 ! relative cellsize factor in search tolerance ()
jasnelucxy = 3 ! method for cell center velocities
! 1 = slow node oriented, 2 = link csu etc, 3= fast linkweight
jarenumber = 1
! Remaining of variables is handled in reset_flowgeom()
call reset_flowgeom()
end subroutine default_flowgeom
!> Resets only flow geometry variables intended for a restart of flow simulation.
!! Upon loading of new model/MDU, use default_flowgeom() instead.
subroutine reset_flowgeom()
! node (s) related : dim=ndx
ndx2D = 0 ! nr of 2d FLOW CELLS = NUMP
ndxi = 0 ! max nr of internal flowcells (internal = 2D + 1D )
ndx1Db = 0 ! nr of flow nodes incl. 1D bnds (internal 2D+1D + 1D bnd)
ndx = 0 ! nr of flow nodes (internal + boundary)
! link (u) related : dim = lnx
lnx1D = 0 ! nr of 1D flow links
lnxi = 0 ! nr of flow links (internal )
lnx1Db = 0 ! nr of flow links incl. 1D bnds (internal 1D+2D + 1D bnd)
lnx = 0 ! nr of flow links (internal + boundary)
! useful parameters :
jaFlowNetChanged = 1 ! To enforce various net(link)-related init routines after renumbering
if (jawave .eq. 4) then ! reinitialize wave directional grid
ntheta = 0
end if
end subroutine reset_flowgeom
end module m_flowgeom
!> this module contains the real flow times, only to be managed by setting times in module m_usertimes
module m_flowtimes
implicit none
character (len=8) :: refdat !< Reference date (e.g., '20090101'). All times (tstart_user, tend_user, etc.) are w.r.t. to this date.
integer :: julrefdat !< will be set by calling settimespacerefdat
double precision :: Tzone !< Data Sources in GMT are interrogated with time in minutes since refdat-Tzone*60
double precision :: Timjan !< time in hours of refdat relative to Januari 1 of the same year
double precision :: dt_user !< User specified time step (s) for external forcing update.
double precision :: dt_nodal !< User specified time step (s) for nodal factors update.
double precision :: dt_max !< Computational timestep limit by user.
double precision :: dt_init !< dt of first timestep, if not specified, use dt_max, if that also not specified, use 1 s
integer :: ja_timestep_auto !< Use CFL-based dt (with dt_max as upper bound)
double precision :: tstart_user !< User specified time start (s) w.r.t. refdat
double precision :: tstop_user !< User specified time stop (s) w.r.t. refdat
double precision :: time_user !< Next time of external forcings update (steps increment by dt_user).
double precision :: dts !< internal computational timestep (s)
double precision :: dti !< inverse computational timestep (1/s)
double precision :: dtprev !< previous computational timestep (s) (1s is a bit like sobek)
double precision :: dtmin !< dt < dtmin : surely crashed
double precision :: dtminbreak !< smallest allowed timestep (in s), checked on a sliding average of several timesteps in validation routine.
double precision, allocatable :: tvalswindow(:) !< (NUMDTWINDOWSIZE) Time1 values in a moving time window to compute sliding average dt
integer , parameter :: NUMDTWINDOWSIZE = 100 !< Number of time steps to include in the sliding average, don't set this too optimistic to avoid too fast simulation breaks.
integer :: idtwindow_start !< Current start index in tvalswindow(:) array. This array is filled in a cyclic order, with never more than NUMDTWINDOWSIZE time values.
double precision :: time0 !< current julian (s) of s0
double precision :: time1 !< current julian (s) of s1 ! and of course, time1 = time0 + dt
double precision :: tim1bnd !< last time boundary signals were given
double precision :: dnt_user !< counter for nr of user steps ( )
double precision :: dnt !< number of timesteps ( )
double precision :: fhr !< Factor sec hrs
double precision :: fday !< Factor sec day
double precision :: ti_map !< map interval (s)
double precision :: ti_maps !< Start of map output period (as assigned in mdu-file) (s)
double precision :: ti_mape !< End of map output period (as assigned in mdu-file) (s)
double precision :: ti_his !< history interval (s)
double precision :: ti_hiss !< Start of his output period (as assigned in mdu-file) (s)
double precision :: ti_hise !< End of his output period (as assigned in mdu-file) (s)
double precision :: ti_wav !< averaging interval spatial wave quantities (s) !! JRE
double precision :: ti_wavs !< averaging interval spatial wave quantities !! JRE
double precision :: ti_wave !< averaging interval spatial wave quantities !! JRE
double precision :: ti_xls !< history interval (s) xls
double precision :: ti_rst !< restart interval (s)
double precision :: ti_rsts !< Start of restart output period (as assigned in mdu-file) (s)
double precision :: ti_rste !< End of restart output period (as assigned in mdu-file) (s)
double precision :: ti_waq !< Interval between output in delwaq files (s).
double precision :: ti_stat !< Interval between simulation statistics output (s).
double precision :: ti_timings !< (parallel) timings output interval
double precision :: ti_split !< Time interval for time splitting: time after which new his/map file will be created (e.g. montly), see also the unit below.
!! Default is 0 to have no time-splitting of output files.
character(len=1) :: ti_split_unit !< Unit for time splitting interval: Y: years, M: months, D: days, h:hours, m: minutes, s: seconds.
double precision, allocatable :: ti_mpt(:) !< times for writing map-files (s), possibly non-equidistant in time
double precision, allocatable :: ti_mpt_rel(:) !< times for writing map-files (s) relative to current time, possibly non-equidistant in time
double precision :: tmini !< Initial time for updating map/his/rst
double precision :: time_choice !< Time consisting the next time_user / time_map
double precision :: time_out !< Next time for output in the most general sense (map, his, etc.)
double precision :: time_map !< Map output interval
double precision :: time_wav !< Time-avg'd output interval xb JRE
double precision :: time_his !< Next time for his output
double precision :: time_xls !< Next time for his output
double precision :: time_rst !< Next time for restart output
double precision :: time_waq !< Next time for delwaq output
double precision :: time_stat !< Next time for simulation statistics output
double precision :: time_timings !< Next time for timings output
double precision :: time_split !< Next time for a new time-split output file.
double precision :: time_split0 !< Start time for the current time-split output file.
integer :: it_map !< Nr of snapshots presently in map file
integer :: it_wav !< Nr of snapshots presently in time-avg'd wave output file JRE
integer :: it_map_tec !< Nr of snapshots presently in map file, Tecplot format
integer :: it_his !< Nr of snapshots presently in his file
integer :: it_inc !< Nr of lines presently in inc file
integer :: it_rst !< Nr of snapshots presently in rst file
integer :: it_waq !< Nr of snapshots presently in delwaq files.
integer :: it_stat !< Nr of simulation statistics presently in log file.
! for performance timings
logical :: debugtimeon !< timing yes or no
double precision :: cpusteps(3) !< cputime timesteps (1)=start,(2)=stop,(3)=total
double precision :: cpuumod (3) !< cputime set-umod (1)=start,(2)=stop,(3)=total
double precision :: cpusol (3) !< cputime conj-grad (1)=start,(2)=stop,(3)=total
double precision :: cpufuru (3) !< cputime conj-grad (1)=start,(2)=stop,(3)=total
double precision :: cpuall (3) !< cputime steps + plots 1,2,3 idem
double precision :: cpuinistep(3) !< cputime inistep 1,2,3 idem
double precision :: cpuext (3) !< cputime externalforcings (1)=start,(2)=stop,(3)=total
double precision :: cpuextbnd (3) !< cputime externalforcingsonbnd (1)=start,(2)=stop,(3)=total
double precision :: cpuwri !< cputime writing (s)
double precision :: dsetb !< number of setbacks ()
character(len=20) :: rundat0 !< start and end date (wallclock) of computer run
character(len=20) :: rundat2 !< start and end date (wallclock) of computer run format = _yymmddhhmmss
character(len=20) :: restartdatetime = ' '!< desired time to be taken from restart map files
double precision :: tlfsmo = 0d0 !< fourier bnd smoothing times
double precision :: alfsmo = 1d0 !< fourier bnd smoothing weight factor
double precision :: Tspinupturblogprof = 0d0 !< From Tstart to Tstart+Tspinupturblogprof, Turbulent profiles based on log profiles
!< 0d0 = No
contains
!> Sets ALL (scalar) variables in this module to their default values.
!! For a reinit prior to flow computation, call reset_flowtimes() instead.
subroutine default_flowtimes()
refdat = '20010101' !< Reference date (e.g., '20090101'). All times (tstart_user, tend_user, etc.) are w.r.t. to this date.
Tzone = 0d0
dt_user = 300d0 !< User specified time step (s) for external forcing update.
dt_nodal = 0d0 !< User specified time step (s) for nodal factors update.
! dt_nodal = 21600d0 !< User specified time step (s) for nodal factors update.
dt_max = 30d0 !< Computational timestep limit by user.
dtmin = 1d-4 !< dt < dtmin : surely crashed
dtminbreak = 0d0 !< smallest allowed timestep, otherwise break: off
dt_init = 1d0
ja_timestep_auto = 1 !< Use CFL-based dt (with dt_max as upper bound)
tstart_user = 0d0 !< User specified time start (s) w.r.t. refdat
tstop_user = 100*24*3600 !< User specified time stop (s) w.r.t. refdat
time_user = tstart_user !< Next time of external forcings update (steps increment by dt_user).
dnt_user = 0 !< counter for nr of user steps ( )
dnt = 0 !< number of timesteps ( )
fhr = 1d0/3600d0 !< Factor sec hrs
fday = 1d0/(3600d0*24d0) !< Factor sec day
ti_map = 1200 !< map interval (s)
ti_wav = 1200 !< wave avg'ing interval (s), 20 minutes okay default JRE
ti_maps = 1200 !< start interval (s)
ti_mape = 1200 !< end interval (s)
ti_his = 120 !< history interval (s)
ti_xls = 0d0 !< history interval (s) xls
ti_rst = 24*3600 !< restart interval (s)
ti_waq = 0d0 !< delwaq interval (s) (Default: off)
ti_stat = -600d0 !< simulation statistics interval (s) (Default: off, will later default to dt_user), <0: use wc-time
ti_timings = 0d0 !< timings output interval
ti_split = 0d0 !< Time interval for time splitting of output files.
ti_split_unit= 'X' !< Unit for time partitioning interval
tmini = -1d9 !< initial time for updating the 4 above
! Remaining of variables is handled in reset_flowtimes()
call reset_flowtimes()
end subroutine default_flowtimes
!> Resets only flow times variables intended for a restart of flow simulation.
!! Upon loading of new model/MDU, use default_flowtimes() instead.
subroutine reset_flowtimes()
dtprev = dt_init !< previous computational timestep (s) (1s is a bit like sobek)
dts = dt_init !< internal computational timestep (s)
time0 = 0d0 !< current julian (s) of s0
time1 = 0d0 !< current julian (s) of s1 ! and of course, time1 = time0 + dt
tim1bnd = -9d9 !< last time bnd signals were given
time_map = tstart_user !< next time for map output
time_wav = tstart_user !< same, xb JRE
time_his = tstart_user !< next time for his output
time_xls = tstart_user !< next time for his output
time_rst = tstart_user !< next time for restart output
time_waq = tstart_user !< next time for waq output
time_stat = tstart_user !< next time for simulation statistics output
time_timings = tstart_user !< next time for timing output
time_split = tstart_user !< next time for a new time-split output file
time_split0 = time_split !< Start time for the current time-split output file.
if (dtminbreak > 0d0) then
if (.not. allocated(tvalswindow)) then
allocate(tvalswindow(NUMDTWINDOWSIZE))
end if
idtwindow_start = 1 ! Start fresh, with first time0 on pos #1.
tvalswindow(idtwindow_start) = tstart_user
end if
it_map = 0 !< Nr of snapshots presently in map file
it_wav = 0 !< Nr of snapshots presently in time-avg'd file JRE
it_map_tec = 0 !< Nr of snapshots presently in map file
it_his = 0 !< Nr of snapshots presently in his file
it_inc = 0 !< Nr of lines presently in inc file
it_rst = 0 !< Nr of snapshots presently in rst file
it_waq = 0 !< Nr of snapshots presently in waq couple files
it_stat = 0 !< Nr of simulation statistics presently in log file.
! for performance timings
debugtimeon = .false. !< timing yes or no
cpusteps (3) = 0 !< cputime timesteps (1)=start,(2)=stop,(3)=total
cpufuru (3) = 0 !< cputime furu (1)=start,(2)=stop,(3)=total
cpusol (3) = 0 !< cputime conj-grad (1)=start,(2)=stop,(3)=total
cpuall (3) = 0 !< cputime steps + plots 1,2,3 idem
cpuinistep(3) = 0 !< cputime steps + plots 1,2,3 idem
cpuext (3) = 0 !< external forcing
cpuextbnd (3) = 0 !< external forcing bnd
cpuwri = 0 !< cputime writing (s)
dsetb = 0 !< number of setbacks ()
alfsmo = 1d0 !<
end subroutine reset_flowtimes
end module m_flowtimes
module m_jacobi ! arrays needed for solving jacobi
implicit none
integer :: ndxjac = 0 ! nr. of nodes already allocated for jacobi should be ndx
integer :: lnxjac = 0 ! nr. of links already allocated for jacobi should be lnx
integer :: itmxjac = 200 ! max nr. of iterations in solve-jacobi
double precision :: epsjac = 1d-10 ! epsilon waterlevels jacobi method (maximum)
double precision, allocatable :: tetaf (:) ! teta(L)*au(L)*fu(L)
double precision, allocatable :: db (:) ! solving related, right-hand side
double precision, allocatable :: bbi (:) ! solving related, diagonal
end module m_jacobi
module m_sferic
implicit none
integer :: jsferic = 0 ! xy pair is in : 0=cart, 1=sferic coordinates
integer :: jsfertek= 0 ! drawn in 0=cart, 1=stereografisch
integer :: jasferdistance = 0 ! 1 = local dx and dy based upon greatcircledistances , 0=original form
integer :: jglobe = 0 ! if (jsferic==1) do we need extra tests for 360-0 transgression
double precision :: pi ! pi
double precision :: twopi ! 2pi
double precision :: dg2rd ! degrees to radians
double precision :: rd2dg ! and vice versa
double precision :: ra = 6378137d0 ! earth radius (m)
double precision :: omega ! earth angular velocity (rad/s)
double precision :: fcorio ! 2omegasinfi
double precision :: anglat = 0d0 ! 26.0 ! dubai 52.5 ! angle of latitude
double precision :: anglon = 0d0 ! 26.0 ! dubai 52.5 ! angle of latitude
double precision :: dy2dg ! from dy in m to lat in degrees
double precision :: csphi ! cosphi of latest requested
double precision, parameter :: dtol_pole = 1d-4 ! pole tolerance in degrees
end module m_sferic
! unstruc.f90
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! solve_guus.f90
module m_reduce
implicit none
!
!
! Derived Type definitions
!
type list
integer l
integer, allocatable :: j(:)
end type list
type listb
integer l
integer, allocatable :: j(:)
logical, allocatable :: b(:)
end type listb
type listc
integer l
integer, allocatable :: j(:)
integer, allocatable :: a(:)
end type listc
type listd
integer l
integer, allocatable :: j(:)
integer, allocatable :: a(:)
integer, allocatable :: f(:)
end type listd
!
! Local variables
!
integer :: maxdge = 6 ! 500
integer :: noactive = 0
integer :: nodtot = 0 ! (nodtot=ndx)
integer :: lintot = 0 ! (lintot=lnx)
integer :: nocg = 0
integer :: nocg0 = 0
integer :: nogauss = 0
integer :: nogauss0 = 0
integer :: noexpl = 0
integer :: nowet = 0
integer :: ijstck = 0
integer :: npmax = 0 ! ccc arrays size in gauss_elimination
integer, allocatable :: noel(:)
integer, allocatable :: noel0(:)
integer, allocatable :: nbrstk(:)
integer, allocatable :: nodstk(:)
integer, allocatable :: nodbr2(:)
integer, allocatable :: lv2(:) ! old linev(2,L), linev(5,L) and (6,L) are now ln(1,L) and ln(2,L)
integer, allocatable :: jagauss(:)
type (listb), allocatable :: ij(:)
type (listc), allocatable :: ia(:)
type (listd), allocatable :: row(:)
double precision :: epscg = 1d-14 ! epsilon waterlevels cg method (maximum)
double precision :: epsdiff = 1d-3 ! tolerance in (outer) Schwarz iterations (for Schwarz solver)
integer :: maxmatvecs = 1000 ! maximum number of matrix-vector multiplications in Saad solver
double precision, allocatable :: bbr(:), bbl(:) ! not left !
double precision, allocatable :: ccr(:), ccrsav(:)
double precision, allocatable :: ddr(:)
double precision, allocatable :: d0 (:)
double precision, allocatable :: zkr(:)
double precision, allocatable :: pk (:)
double precision, allocatable :: apk(:)
double precision, allocatable :: rk (:)
double precision, allocatable :: ccc(:) !< work array in gauss_elimination
integer :: nbr
integer :: nodl
integer :: nodr
integer :: ndn
integer :: mindgr
integer :: nocgiter
double precision, allocatable, dimension(:) :: s1_ghost ! for testsolver
end module m_reduce
! solve_guus.f90
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! rest.f90
MODULE M_BITMAP
implicit none
INTEGER, ALLOCATABLE, SAVE :: IPIX(:)
INTEGER, SAVE :: MXP, NXP
double precision, SAVE :: XB(4), YB(4), XP(4), YP(4)
END MODULE M_BITMAP
MODULE M_ARCINFO
implicit none
double precision, ALLOCATABLE :: D(:,:)
INTEGER :: MCa = 0, NCa
double precision :: X0=0, Y0=0, DXa=1d0, DYa=1d0, RMIS=-999d0
END MODULE M_ARCINFO
!> Main sample set
MODULE M_SAMPLES
implicit none
double precision, ALLOCATABLE :: XS(:), YS(:), ZS(:)
INTEGER, ALLOCATABLE :: IPSAM(:) !< permutation array (increasing x-coordinate)
integer, parameter :: IPSTAT_OK=0 !< permutation array is OK
integer, parameter :: IPSTAT_NOTOK=1 !< permutation array is out of date
integer :: IPSTAT=IPSTAT_NOTOK !< permutation array status
INTEGER :: NS =0, NSMAX
integer :: MXSAM=0, MYSAM=0 !< structured block sizes (.gt.0), or not structured (0)
double precision :: xsammin, ysammin, xsammax, ysammax !< bounding box corner coordinates
END MODULE M_SAMPLES
!> Backup of main sample set
!! @see savesam()
MODULE M_SAMPLES2
implicit none
double precision, ALLOCATABLE :: XS2(:), YS2(:), ZS2(:)
INTEGER :: NS2, MXSAM2=0, MYSAM2=0
END MODULE M_SAMPLES2
!> Alternate sample set.
!! @see SWAPSAMPLES()
MODULE M_SAMPLES3
implicit none
double precision, ALLOCATABLE :: XS3(:), YS3(:), ZS3(:)
INTEGER :: NS3
END MODULE M_SAMPLES3
MODULE M_TRIANGLE
implicit none
double precision, ALLOCATABLE :: XCENT(:), YCENT(:)
INTEGER, ALLOCATABLE :: INDX(:,:)
INTEGER, ALLOCATABLE :: EDGEINDX(:,:)
INTEGER, ALLOCATABLE :: TRIEDGE(:,:)
INTEGER :: NUMTRI
INTEGER :: NUMTRIINPOLYGON
INTEGER :: NUMEDGE
INTEGER, PARAMETER :: ITYPE = 2 ! 1 = ORIGINAL FORTRAN ROUTINE, 2 = NEW C ROUTINE
integer :: jagetwf = 0 ! if 1, also assemble weightfactors and indices in:
INTEGER, ALLOCATABLE :: indxx(:,:) ! to be dimensioned by yourselves 3,*
double precision, ALLOCATABLE :: wfxx (:,:)
double precision :: TRIANGLEMINANGLE = 5d0 ! MINIMUM ANGLE IN CREATED TRIANGLES IF MINANGLE > MAXANGLE: NO CHECK
double precision :: TRIANGLEMAXANGLE = 150 ! MAXIMUM ANGLE IN CREATED TRIANGLES
double precision :: TRIANGLESIZEFAC = 1.0 ! TRIANGLE SIZEFACTOR, SIZE INSIDE VS AVERAGE SIZE ON POLYGON BORDER
TYPE T_NODI
INTEGER :: NUMTRIS ! total number of TRIANGLES ATtached to this node
INTEGER, allocatable :: TRINRS(:) ! numbers of ATTACHED TRIANGLES
END TYPE T_NODI
TYPE (T_NODI), DIMENSION(:), ALLOCATABLE :: NODE !
!integer, dimension(:,:), allocatable :: trinods ! triangle nodes, dim(3,numtri)
integer, dimension(:,:), allocatable :: LNtri ! triangles connected to edges, dim(2,numedges)
integer :: IDENT ! identifier
integer, dimension(:), allocatable :: imask ! mask array for triangles
END MODULE M_TRIANGLE
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! net.f90
module m_dimens
implicit none
integer :: MMAX_old = 3, NMAX_old = 3
integer :: KMAX, LMAX, KNX, MXB
end module m_dimens
module m_landboundary
implicit none
double precision, allocatable :: XLAN (:), YLAN(:), ZLAN(:)
integer, allocatable :: NCLAN(:)
integer :: MXLAN, MAXLAN
! SPvdP: segments
integer :: MXLAN_loc ! actual MXLAN
integer :: Nlanseg ! number of land boundary segments
integer, allocatable, dimension(:,:) :: lanseg_startend ! segment start and end indices, dim(2,Nlanseg)
integer, allocatable, dimension(:) :: lanseg_map ! node to land boundary segment mapping, dim(numk)
integer :: jleft, jright !< outer land boundary segments in projection
double precision :: rLleft, rLright !< fractional location of the projected outer nodes (min and max) on the land boundary segment
double precision :: DCLOSE_bound = 5d0 !< close-to-landboundary tolerance, netbound only, measured in number of meshwidths
double precision :: DCLOSE_whole = 1d0 !< close-to-landboundary tolerance, whole network, measured in number of meshwidths
double precision :: DCLOSE = 1d0 ! close-to-landboundary tolerance, measured in number of meshwidths
logical :: Ladd_land = .true. ! add land boundary between land boundary segments that are close to each other
end module m_landboundary
Module m_oldz
implicit none
double precision :: OZ = 999
end module m_oldz
module m_polygon
implicit none
double precision, allocatable :: XPL (:), YPL (:), ZPL (:), XPH(:), YPH(:), ZPH(:), DZL(:), DZR(:)
integer :: NPL, NPH, MAXPOL, MP, MPS, jakol45 = 0
character(len=64), allocatable :: nampli(:) ! Names of polylines, set in reapol,
! not shifted/updated during editpol.
end module m_polygon
module m_makenet
implicit none
integer :: NTYP = 0, NRX = 3, NRY = 3
double precision :: ANGLE = 0, SIZE = 50, THICK = 8, HSIZE = 50
double precision :: X0 = 0, Y0 = 0, Z0 = 0, DX0 = 10, DY0 = 10
end module m_makenet
module m_mergenet
implicit none
integer :: NUMM = 2, JXYZ = 1
end module m_mergenet
module m_seastate
implicit none
double precision :: TWOPI, WAVLEN, WAVKX, WAVOM
double precision :: WAVCEL = 5 , WAVPER = 15, WAVAMP = .1
end module m_seastate
module m_boat
implicit none
double precision, allocatable :: XBOAT (:), YBOAT(:), ZBOAT(:) ! INLEZEN &
double precision, allocatable :: XBOOT (:), YBOOT(:), ZBOOT(:) ! AFBEELDEN LIJNEN
integer, allocatable :: NCBOAT(:) !
integer :: MXBOAT, MAXBOAT, NCLBOAT
double precision :: BLEN = 25, BHEIGHT = 4, BWIDTH = 5, BHPMAX = 1000, BHPPERC = 0.5d0
integer :: KKB(10) = 0
end module m_boat
module M_afmeting
implicit none
double precision :: RLENGTH, RWIDTH, RTHICK, RDIAM, RLMIN
integer :: JVAST, MC, NC, K0, LFAC
end module M_afmeting
MODULE M_ISOSCALEUNIT
implicit none
CHARACTER (LEN=10) :: UNIT(2) = ' '
CHARACTER (LEN=20) :: PARAMTEX(2) = ' '
END MODULE M_ISOSCALEUNIT
MODULE M_RAAITEK
implicit none
double precision :: ZMINrai=-999, ZMAXrai=-999
INTEGER :: JTEXTFLOW = 1
real :: xs1m,ys1m,xs2m,ys2m
double precision :: xw1m,yw1m,xw2m,yw2m ! screen and world corners in module
END MODULE M_RAAITEK
module m_grid
implicit none
integer :: MET = 1 ! drawing method
integer :: mc = 0 , nc = 0 ! actuele dimensies
integer :: MMAX = 0 , NMAX = 0 , MNMAX = 0 ! array dimensies
integer :: mch = 0 , nch = 0 ! actuele dimensies van hulprooster
double precision, allocatable :: xc (:,:), yc (:,:), zc (:,:)
double precision, allocatable :: xch(:,:), ych(:,:), zch(:,:) ! save/restore
integer, allocatable :: ijc(:,:), ijyes(:,:)
end module m_grid
!> Regular grid generation settings. All orthogonalisation settings are in
!! module m_orthosettings.
MODULE M_GRIDSETTINGS
implicit none
integer :: MFAC = 2000 !< M-refinement factor for regular grid generation.
integer :: NFAC = 40 !< N-refinement factor for regular grid generation.
integer :: ITSMO = 10 !< Nr. of inner iterations in regular grid smoothing.
integer :: ITSMA !< Not in use, old rgfgrid
integer :: JADEPDESIGN = 0
integer :: MDESIGN
double precision :: BFAC=1d0, CSMO = 0.5d0, RFAC
double precision :: SRM,SRN,DEPSLO,FSMA, ALINEN, ALINEM
INTEGER :: KEEPSTARTDIR = 1
double precision :: BAAS2 = 0.5d0, FACMIR = 1.2d0
double precision :: SPLFAC, SPLFAC2
INTEGER :: JDEMO = 0
! Pillar grid settings
double precision :: pil_rad = 0d0 !< pillar radius
double precision :: pil_x = 0d0 !< pillar center point x-coordinate
double precision :: pil_y = 0d0 !< pillar center point y-coordinate
double precision :: pil_grow = 1d0 !< pillar grid growth factor *not used*
END MODULE M_GRIDSETTINGS
!> Orthogonalisation settings, both for regular grids and unstructured nets.
module m_orthosettings
implicit none
integer :: ITATP = 25 !< Nr. of outer iterations in grid/net orthogonalisation.
integer :: ITBND = 1 !< Nr. of boundary iterations in grid/net orthogonalisation. (within ITATP)
integer :: ITIN = 25 !< Nr. of inner iterations in grid/net orthogonalisation. (within ITBND)
!! Also used within transfinite regular grid generation.
integer :: JAPROJECT = 1 !< Project nodes back to boundary (2: yes, all, 1:yes, net bounds only, 0:no)
double precision :: ATPF = 0.95d0 !< Factor (0.<=ATPF<=1.) between grid smoothing and grid ortho resp.
double precision :: ATPF_B = 1d0 !< minimum ATPF on the boundary
double precision :: circumormasscenter = 1d0 !< 1.0 = circumcentre, 0.0 = masscentre, 1.0 -> 0.0 : weighted
double precision :: smoothorarea = 1d0 !< Factor between smoother (1d0) and area-homogenizer (0d0)
integer :: adapt_method = 1 !< Mesh-adaptation method; 0: Winslow, 1: arc-length, 2: harmonic map
double precision :: adapt_beta = 0.0d0 !< Mesh-refinement factor; between 0d0 and 1d0
integer :: adapt_niter_u = 0 !< number of smoothing iterations of `solution` u in adaptation
integer :: adapt_niter_G = 4 !< number of smoothing iterations of monitor matrix G in adaptation
double precision :: ortho_pure = 0.5d0 !< curvi-linear-like (0d0) or pure (1d0) orthogonalisation
end module m_orthosettings
MODULE M_INTERPOLATIONSETTINGS
implicit none
integer, parameter :: INTP_INTP = 1
integer, parameter :: INTP_AVG = 2
INTEGER :: INTERPOLATIONTYPE ! 1 = TRIANGULATION/BILINEAR INTERPOLATION 2= CELL AVERAGING
INTEGER :: JTEKINTERPOLATIONPROCESS ! TEKEN INTERPOLATION PROCESS YES/NO 1/0
INTEGER :: IAV ! AVERAGING METHOD, 1 = SIMPLE AVERAGING, 2 = CLOSEST POINT, 3 = MAX, 4 = MIN, 5 = INVERSE WEIGHTED DISTANCE, 6 = MINABS, 7 = KDTREE
INTEGER :: NUMMIN ! MINIMUM NR OF POINTS NEEDED INSIDE CELL TO HANDLE CELL
DOUBLE PRECISION, parameter :: RCEL_DEFAULT = 1.01d0 ! we need a default
DOUBLE PRECISION :: RCEL ! RELATIVE SEARCH CELL SIZE, DEFAULT 1D0 = ACTUAL CELL SIZE, 2D0=TWICE AS LARGE
INTEGER :: Interpolate_to ! 1=bathy, 2=zk, 3=s1, 4=Zc
contains
!> set default interpolation settings
subroutine default_interpolationsettings()
implicit none
INTERPOLATIONTYPE = INTP_INTP ! 1 = TRIANGULATION/BILINEAR INTERPOLATION 2= CELL AVERAGING
JTEKINTERPOLATIONPROCESS = 0 ! TEKEN INTERPOLATION PROCESS YES/NO 1/0
IAV = 1 ! AVERAGING METHOD, 1 = SIMPLE AVERAGING, 2 = CLOSEST POINT, 3 = MAX, 4 = MIN, 5 = INVERSE WEIGHTED DISTANCE, 6 = MINABS, 7 = KDTREE
NUMMIN = 1 ! MINIMUM NR OF POINTS NEEDED INSIDE CELL TO HANDLE CELL
RCEL = RCEL_DEFAULT ! RELATIVE SEARCH CELL SIZE, DEFAULT 1D0 = ACTUAL CELL SIZE, 2D0=TWICE AS LARGE
Interpolate_to = 2 ! 1=bathy, 2=zk, 3=s1, 4=Zc
return
end subroutine default_interpolationsettings
END MODULE M_INTERPOLATIONSETTINGS
MODULE M_XYTEXTS
implicit none
integer, parameter :: maxtxt = 2000
double precision :: xtxt(maxtxt), ytxt(maxtxt),heitxt(maxtxt)
integer :: symtxt(maxtxt), coltxt(maxtxt)
integer :: ntxt
character(len=120) :: xytexts(maxtxt)
END MODULE M_XYTEXTS
module m_equatorial
double precision :: x,fr,Ue0,k,h,g,L,utyp,period,ap,am,om,Zp,ufac
integer :: nmode=0, nfreq=1
integer :: ndxfreeL, ndxforced, ndtfreeL, ndtforced
double precision :: app = 0.0d0 , amm = 0d0, ztyp = 0.20d0
double precision :: ampliforced, amplitotal, amplifreeL, cflfreeL, cfLforced, TfreeL, Tforce, amplicomp
end module m_equatorial
!> inverse-map smoother in orthogonalisenet
module m_inverse_map
type tops !< operator type
double precision, allocatable, dimension(:, :) :: Az !< cell-center coefficient matrix; node-to-cell
double precision, allocatable, dimension(:, :) :: Gxi, Geta !< netcell-gradient coefficient matrices; node-to-link
double precision, allocatable, dimension(:) :: Divxi, Diveta !< netnode-gradient coefficient matrices; link-to-node
double precision, allocatable, dimension(:) :: Jxi, Jeta !< netnode-gradient coefficient matrices; node-to-node
double precision, allocatable, dimension(:) :: ww2 !< weights in Laplacian smoother
end type
type tadm !< administration type (per node)
integer :: Ncell !< number of netcells connected to the center node
integer, allocatable, dimension(:) :: icell !< netcells connected to node k0
integer :: nmk
integer :: nmk2
integer, allocatable, dimension(:) :: kk2 ! local node administration
integer, allocatable, dimension(:,:) :: kkc ! position of netcell nodes in kk2 array
end type
type ttop !< topology array type (for unique topologies)
integer, allocatable, dimension(:) :: nmk, nmk2 !< stored number of links and nodes in stencil resp.
double precision, allocatable, dimension(:,:) :: xi, eta !< stored node coordinates (xi, eta)
end type
!---------------------------
! administration
!---------------------------
integer :: nmkx !< maximum number of links connected to center node k0
integer, save :: nmkx2=4 !< maximum number of nodes in stencil
integer, allocatable, dimension(:) :: nmk2 !< number of nodes in stencil
integer, allocatable, dimension(:,:) :: kk2 !< node administration; first row, i.e. kk2(1,:), points to center nodes
!---------------------------
! unique topologies
!---------------------------
integer :: numtopo !< number of unique topologies
integer, allocatable, dimension(:) :: ktopo !< index in topology array
!--------------------------------------------
! parameters
!--------------------------------------------
integer, parameter :: M=6 !< maximum number of nodes per netcell
!---------------------------
end module m_inverse_map
!> global variables for spline2curvi
module m_spline2curvi
integer, parameter :: Nsubmax = 10 !< maximum number of subintervals of grid layers, each having their own exponential grow factor
type tspline !< center spline type that contains information derived from cross splines
! cross spline data
integer :: id !< 0: center spline, 1: true cross spline, 2: unassociated bounding spline, 3: artificial cross spline, -i, where i>0: associated bounding spline, with centerspline -i
integer :: ncs !< number of cross splines
double precision :: length !< spline path length
double precision :: hmax !< maximum grid height
integer, allocatable, dimension(:) :: ics !< cross spline numbers
Logical, allocatable, dimension(:) :: Lorient !< cross spline is left to right (.true.) or not (.false.) w.r.t. center spline
double precision, allocatable, dimension(:) :: t !< center spline coordinates of cross splines
double precision, allocatable, dimension(:) :: cosphi !< cosine of crossing angle
double precision, allocatable, dimension(:,:) :: hL !< left-hand side grid heights at cross spline locations for each grid layer subinterval, hL(1,:) being the height of the first subinterval, etc.
double precision, allocatable, dimension(:,:) :: hR !< right-hand side grid heights at cross spline locations for each grid layer subinterval, hR(1,:) being the height of the first subinterval, etc.
integer, allocatable, dimension(:) :: NsubL !< number of subintervals of grid layers at cross spline locations at the left-hand side of the spline, each having their own exponential grow factor
integer, allocatable, dimension(:) :: NsubR !< number of subintervals of grid layers at cross spline locations at the right-hand side of the spline, each having their own exponential grow factor
! grid data
integer :: mfac !< number of grid intervals on the spline
integer, dimension(Nsubmax) :: nfacL !< number of grid layers in each subinterval at the left-hand side of the spline *not used yet*
integer, dimension(Nsubmax) :: nfacR !< number of grid layers in each subinterval at the right-hand side of the spline *not used yet*
integer :: iL !< index in the whole gridline array of the first grid point on the left-hand side of the spline
integer :: iR !< index in the whole gridline array of the first grid point on the right-hand side of the spline
end type
type(tspline), allocatable, dimension(:) :: splineprops !< spline properties
double precision, allocatable, dimension(:) :: xg1, yg1 !< coordinates of the first gridline
double precision, allocatable, dimension(:) :: sg1 !< center spline coordinates of the first gridline
integer :: jacirc = 0 !< circularly connected grid (1) or not (0) (for netbound2curvi), note: valid for one gridline only
integer :: jaoutside = 1 ! grow the grid outside the prescribed grid height
double precision :: daspect = 0.1d0 ! aspect ratio
double precision :: dgrow = 1.1d0 ! grow factor of aspect ratio
double precision :: dheight0 = 1.0d1 ! grid layer height
double precision :: maxaspect = 1.0d0 ! maximum cell aspect ratio *inoperative*
double precision :: dwidth = 0.5d3 ! average mesh width on center spline
double precision :: dtolLR = 1d-4 ! on-top-of-each-other tolerance *IMPORTANT*
double precision :: dtolcos = 0.95d0 ! minimum allowed absolute value of crossing-angle cosine
integer :: NFACUNIMAX = 5 ! maximum number of layers in the uniform part
integer :: jaCheckFrontCollision = 0 ! check for collisions with other parts of the front (1) or not (0)
double precision :: dunigridsize = 0d0 ! uniform grid size (netboundary to grid only)
integer :: jacurv = 1 ! curvature adapted grid spacing (1) or not (0)
end module m_spline2curvi
module m_samples_refine ! used in refinecellsandfaces2 and in sample paths
integer :: NDIM=5 !< sample vector dimension
double precision, allocatable, dimension(:,:,:) :: zss !< sample data [zs, direction vector x-component, direction vector y-component, refinement criterion, ridge distance], dim(NDIM,MXSAM,MYSAM)
integer :: Nsamplesmooth = 0 !< number of sample smoothing iterations
integer :: Nsamplesmooth_last = -1 !< last number of sample smoothing iterations
integer :: MAXLEVEL = 10 !< maximum number of refinement levels
double precision :: threshold = 1d2 !< typical obstacle height in grid refinement
double precision :: thresholdmin = 1d0 !< minimum obstacle height grid refinement
double precision :: hmin = 5d0 !< minimum cell size
integer :: jadirectional = 0 !< directional refinement (1) or not (0)
integer, parameter :: iHesstat_OK = 0 !< sample Hessians up-to-date
integer, parameter :: iHesstat_DIRTY = 1 !< sample Hessians out-of-date
integer :: iHesstat = 0 !< sample Hessians up-to-date (0) or not (1)
integer, parameter :: ITYPE_RIDGE = 1 !< critetrion based on ridge-detection
integer, parameter :: ITYPE_WAVECOURANT = 2 !< critetrion based on wave Courant number
integer :: irefinetype = ITYPE_WAVECOURANT !< refinement criterion type
integer :: jaconnect = 1 !< connect hanging nodes (1) or not (0)
double precision :: Dt_maxcour = 0d0 !< maximum time-step in courant grid
double precision :: dminsampledist = 0d0 !< minimum sample distance
integer :: jaoutsidecell = 1 !< take samples outside cell into account (1) or not (0)
end module m_samples_refine
module m_kml_parameters
implicit none
integer :: kml_janet !< Whether or not (1/0) to export flat view of 2D+1D grid (faster)
integer :: kml_jadepth !< Whether or not (1/0) to export bathymetry view of grid cells (nicer).
integer :: kml_jadepth3d !< Whether or not (1/0) to export bathymetry view in 3D.
double precision :: kml_altfact !< Altitude exaggeration factor: altitude differences are multiplied by this.
integer :: kml_jaoffsetzk !< Whether or not (1/0) to offset all altitudes with deepest zk-value.
double precision :: kml_useroffset !< Additional user offset for altitude values.
double precision :: kml_dmiss !< Dummy altitude to replace missing zk values.
double precision :: kml_zmin, kml_zmax !< Min/max values used for color scaling.
contains
!> This subroutine should be called during program initialization.
subroutine default_kml_parameters()
kml_janet = 1 !< Whether or not (1/0) to export flat view of 2D+1D grid (faster)
kml_jadepth = 0 !< Whether or not (1/0) to export bathymetry view of grid cells (nicer).
kml_jadepth3d = 0 !< Whether or not (1/0) to export bathymetry view in 3D.
kml_altfact = 5 !< Altitude exaggeration factor: altitude differences are multiplied by this.
kml_jaoffsetzk = 1 !< Whether or not (1/0) to offset all altitudes with deepest zk-value.
kml_useroffset = 0d0 !< Additional user offset for altitude values.
kml_dmiss = 99d0 !< Dummy altitude to replace missing zk values.
end subroutine default_kml_parameters
end module m_kml_parameters
module m_timer
implicit none
integer, parameter :: jatimer = 1 !< time parallel solver (1) or not (0)
integer, parameter :: NUMT=21 !< number of timings
double precision, dimension(3,NUMT) :: t !< wall-clock timings, (1,:): start, (2,:): end, (3,:): sum
double precision, dimension(3,NUMT) :: tcpu !< CPU timings, (1,:): start, (2,:): end, (3,:): sum
integer, dimension(NUMT) :: itstat !< timer status, 0: not timing (stopped), 1: timing (started)
integer :: numtsteps !< number of time steps
integer :: numcgits !< number of CG iterations
integer, parameter :: IREDUCE = 1
integer, parameter :: IMAKEMAT = 2
integer, parameter :: IPACKMAT = 3
integer, parameter :: IGAUSSEL = 4
integer, parameter :: ICG = 5
integer, parameter :: IGAUSSSU = 6
integer, parameter :: IMPICOMM = 7
integer, parameter :: ITOTALSOL = 8
integer, parameter :: ITIMESTEP = 9
integer, parameter :: ITOTAL = 10
integer, parameter :: IPOSTMPI = 11
integer, parameter :: IUPDSALL = 12
integer, parameter :: IUPDU = 13
integer, parameter :: IMPIREDUCE = 14
integer, parameter :: IPARTINIT = 15
integer, parameter :: IOUTPUT = 16
integer, parameter :: IOUTPUTMPI = 17
integer, parameter :: ITRANSPORT = 18
integer, parameter :: IXBEACH = 19
integer, parameter :: IAXPY = 20
integer, parameter :: IDEBUG = 21
character(len=10), dimension(numt), parameter :: tnams= [character(len=10) :: &
'reduce', &
'make_mat.', &
'pack', &
'Gauss_elem', &
'CG', &
'Gauss_subs', &
'MPI_comm', &
'total_sol', &
'timestep', &
'total', &
'MPI_post', &
'MPI_sall', &
'MPI_u', &
'MPI_reduce', &
'part_init', &
'output', &
'MPI_output', &
'transport', &
'XBeach', &
'Axpy', &
'debug']
contains
!> initialize timers
subroutine initimer()
implicit none
t = 0d0
tcpu = 0d0
itstat = 0
numtsteps = 0
numcgits = 0
end subroutine initimer
!> start the timer
subroutine starttimer(itvar)
implicit none
integer, intent(in) :: itvar !< timer number
double precision :: tloc
integer :: i
! check status
if ( itstat(itvar).ne.0 ) then
write (6,'("starttimer: status error for timer ", I0)') itvar
end if
call klok(tloc)
t(1,itvar) = tloc
call cpu_time(tcpu(1,itvar))
! set status
itstat(itvar) = 1
return
end subroutine starttimer
!> stop the timer
subroutine stoptimer(itvar)
use MessageHandling
implicit none
integer, intent(in) :: itvar !< timer number
double precision :: tloc
double precision, parameter :: dtol=1d-3 !< timer tolerance
integer :: i
! check status
if ( itstat(itvar).ne.1 ) then
write (6,'("stoptimer: status error for timer ", I0)') itvar
else
call klok(tloc)
t(2,itvar) = tloc
t(3,itvar) = t(3,itvar) + tloc - t(1,itvar)
call cpu_time(tcpu(2,itvar))
tcpu(3,itvar) = tcpu(3,itvar) + tcpu(2,itvar) - tcpu(1,itvar)
! if ( (tcpu(3,itvar)-t(3,itvar)).gt.dtol ) then ! should not happen
! write(6,*) tnams(itvar)
! write(6,*) 't=', t(:,itvar)
! write(6,*) 'tcpu=', tcpu(:,itvar)
! call mess(LEVEL_ERROR, 'stoptimer: timer error')
! end if
end if
! set status
itstat(itvar) = 0
return
end subroutine stoptimer
!> get timer value
double precision function gettimer(itype, itvar)
implicit none
integer, intent(in) :: itype !< timer type, cpu-time (0) or wall-clock time (other)
integer, intent(in) :: itvar !< timer number
if ( itype.eq.0 ) then
gettimer = tcpu(3,itvar)
else
gettimer = t(3,itvar)
end if
return
end function gettimer
end module m_timer
!> process command line option
module m_commandline_option
implicit none
integer, parameter :: MAXOPTLEN = 255 !< maximum option string length
integer, parameter :: MAXKEYLEN = 255 !< maximum key string lenght
integer, parameter :: MAXKEYS = 16 !< maximum number of key-value pairs
integer, parameter :: IMISS = -999999 !< unset value
! input files from command line
integer, parameter :: lenfile = 255 !< maximum filename length
integer, parameter :: maxnumfiles = 10 !< maximum number of files
integer :: numfiles !< number of files to be loaded
character(len=lenfile), dimension(maxnumfiles) :: inputfiles !< files to be loaded
integer :: iarg_autostart !< autostart/autstartstop or not set (-1)
contains
!> read, from a string, command line option with key-value pair(s) of the form
!> <->[-[...]]