!----- 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 !> <->[-[...]]