module params use mnemmodule use xmpi_module type parameters ! These parameters are constants, variables read from params.txt, or are scalars derived directly from read input ! ! Type name initialize ! [unit] description ! Physical processes integer*4 :: swave = -123 ! [-] Include short waves (1), exclude short waves (0) integer*4 :: lwave = -123 ! [-] Include short wave forcing on NLSW equations and boundary conditions (1), or exclude (0) integer*4 :: flow = -123 ! [-] Include flow calculation (1), or exclude (0) integer*4 :: sedtrans = -123 ! [-] Include sediment transport (1) or exclude (0) integer*4 :: morphology = -123 ! [-] Include morphology (1) or exclude (0) integer*4 :: nonh = -123 ! [-] Non-hydrostatic pressure option: 0 = NSWE, 1 = NSW + non-hydrostatic pressure compensation Stelling & Zijlema, 2003 integer*4 :: gwflow = -123 ! [-] Turn on (1) or off (0) groundwater flow module integer*4 :: q3d = -123 ! [-] Turn on (1) or off (0) quasi-3D sediment transport module ! Grid parameters character(256):: depfile = 'abc' ! [-] Name of the input bathymetry file real*8 :: posdwn = -123 ! [-] Bathymetry is specified positive down (1) or positive up (-1) integer*4 :: nx = -123 ! [-] Number of computiation cell corners in x-direction integer*4 :: ny = -123 ! [-] Number of computiation cell corners in y-direction real*8 :: alfa = -123 ! [deg] Angle of x-axis from East integer*4 :: vardx = -123 ! [-] Switch for variable grid spacing: 1 = irregular spacing, 0 = regular grid spacing real*8 :: dx = -123 ! [m] Regular grid spacing in x-direction real*8 :: dy = -123 ! [m] Regular grid spacing in y-direction character(256) :: xfile = 'abc' ! [name] Name of the file containing x-coordinates of the calculation grid character(256) :: yfile = 'abc' ! [name] Name of the file containing y-coordinates of the calculation grid real*8 :: xori = -123 ! [m] X-coordinate of origin of axis real*8 :: yori = -123 ! [m] Y-coordinate of origin of axis real*8 :: thetamin = -123 ! [deg] Lower directional limit (angle w.r.t computational x-axis) real*8 :: thetamax = -123 ! [deg] Higher directional limit (angle w.r.t computational x-axis) real*8 :: dtheta = -123 ! [deg] Directional resolution integer*4 :: thetanaut = -123 ! [-] Thetamin,thetamax in cartesian (0) or nautical convention (1) ! Model time real*8 :: tstop = -123 ! [s] Stop time of simulation, in morphological time real*8 :: CFL = -123 ! [-] Maximum Courant-Friedrichs-Lewy number character*80 :: tunits = 's' ! Units can be defined in udunits format (seconds since 1970-01-01 00:00:00.00 +1:00) ! Physical constants real*8 :: g = -123 ! [ms^-2] Gravitational acceleration real*8 :: rho = -123 ! [kgm^-3] Density of water ! Initial conditions character(256):: zsinitfile = 'abc' ! [name] Name of inital condition file zs ! Wave boundary condition parameters character(24) :: instat = 'abc' ! [-] Wave boundary condtion type real*8 :: taper = -123 ! [s] Spin-up time of wave boundary conditions, in morphological time real*8 :: Hrms = -123 ! [m] Hrms wave height for instat = 0,1,2,3 real*8 :: Tm01 = -123 ! [s] Old name for Trep real*8 :: Trep = -123 ! [s] Representative wave period for instat = 0,1,2,3 real*8 :: Tlong = -123 ! [s] Wave group period for case instat = 1 real*8 :: dir0 = -123 ! [deg] Mean wave direction (Nautical convention) for instat = 0,1,2,3 integer*4 :: m = -123 ! [-] Power in cos^m directional distribution for instat = 0,1,2,3 character(24) :: lateralwave = 'neumann' ! [name] Switch for lateral boundary at left, 'neumann' = E Neumann, 'wavefront' = along wave front character(24) :: leftwave = 'abc' ! [-] old name for lateralwave character(24) :: rightwave = 'abc' ! [-] old name for lateralwave ! Wave-spectrum boundary condition parameters character(256):: bcfile = 'abc' ! Note, will replace current lookup in boundary conditions [name] Name of spectrum file integer*4 :: random = -123 ! [-] Random seed on (1) or off (0) for instat = 4,5,6 boundary conditions real*8 :: fcutoff = -123 ! [Hz] Low-freq cutoff frequency for instat = 4,5,6 boundary conditions integer*4 :: nspr = -123 ! [-] nspr = 1 bin all wave components for generation of qin (instat 4,5,6) in one direction, nspr = 0 regular long wave spreadin real*8 :: trepfac = -123 ! [-] Compute mean wave period over energy band: par%trepfac*maxval(Sf) for instat 4,5,6; converges to Tm01 for trepfac = 0.0 and real*8 :: sprdthr = -123 ! [-] Threshold ratio to maxval of S above which spec dens are read in (default 0.08*maxval) integer*4 :: oldwbc = -123 ! [-] (1) Use old version wave boundary conditions for instat 4,5,6 real*8 :: rt = -123 ! [s] Duration of wave spectrum at offshore boundary, in morphological time real*8 :: dtbc = -123 ! [s] Timestep used to describe time series of wave energy and long wave flux at offshore boundary (not affected by morfac) real*8 :: dthetaS_XB = -123 ! [deg] If SWAN input is not in nautical degrees, dthetaS_XB is the angle from SWAN x-axis to XBeach x-axis in cathesian degrees ! Flow boundary condition parameters character(24) :: front = 'abc' ! [-] Switch for seaward flow boundary: 0 = radiating boundary(Ad), 1 = Van Dongeren, 1997 character(24) :: left = 'abc' ! [name] Switch for lateral boundary at ny+1, 'neumann' = vv computed from NSWE, 'wall' = reflective wall; vv=0 character(24) :: right = 'abc' ! [-] Switch for lateral boundary at right, 0 = vv computed from NSWE, 1 = reflective wall; vv=0 character(24) :: back = 'abc' ! [-] Switch for boundary at bay side, 0 = radiating boundary (Ad), 1 = reflective boundary; uu=0 integer*4 :: ARC = -123 ! [-] Switch for active reflection compensation at seaward boundary: 0 = reflective, 1 = weakly (non) reflective real*4 :: order = -123 ! [-] 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*4 :: carspan = -123 ! [-] Switch for Carrier-Greenspan test 0 = use cg (default); 1 = use sqrt(gh) in instat = 3 for c&g tests real*8 :: epsi = -123 ! [-] Ratio of mean current to time varying current through offshore boundary ! Tide boundary conditions real*8 :: zs0 = -123 ! [m] Inital water level character(256):: zs0file = 'abc' ! Note, will replace lookup in readtide [name] Name of tide boundary condition series integer*4 :: tideloc = -123 ! [-] Number of corner points on which a tide time series is specified integer*4 :: paulrevere = -123 ! [-] Specifies tide on sea and land (0) or two sea points (1) if tideloc = 2 ! if tideloc =>2, then this indicates where the time series are to be ! applied. Input for tidal information to xbeach options (3): ! 1. one tidal record --> specify tidal record everywhere ! 2. two tidal records --> Need to specify keyword 'paulrevere' ! paulrevere==0 implies to apply one tidal record to ! both sea corners and one tidal record to both land corners ! paulrevere==1 implies to apply the first tidal record ! (column 2 in zs0input.dat) to the (x=1,y=1) sea corner and ! the second tidal record (third column) to the (x=1,y=N) sea corner ! 3. four tidal records --> Need to list tidal records in ! zs0input.dat in order of: ! (x=1,y=1) ! (x=1,y=N) ! (x=N,y=N) ! (x=N,y=1) ! NOTE: clockwise from (1,1) corner ! Discharge boundary conditions character(256):: disch_loc_file = 'abc' ! Note: will replace lookup in boundary conditions [name] Name of discharge locations file character(256):: disch_timeseries_file = 'abc' ! Note: will replace lookup in boundary conditions [name] Name of discharge timeseries file ! Wave breaking parameters character(24) :: break = 'abc' ! [-] Type of breaker formulation (1=roelvink, 2=baldock, 3=roelvink adapted, 4=roelvink on/off breaking) real*8 :: gamma = -123 ! [-] Breaker parameter in Baldock or Roelvink formulation real*8 :: gamma2 = -123 ! [-] End of breaking parameter in break = 4 formulation real*8 :: alpha = -123 ! [-] Wave dissipation coefficient in Roelvink formulation real*8 :: n = -123 ! [-] Power in Roelvink dissipation model real*8 :: gammax = -123 ! [-] Maximum ratio wave height to water depth real*8 :: delta = -123 ! [-] Fraction of wave height to add to water depth real*8 :: fw = -123 ! [-] Bed friction factor integer*4 :: breakerdelay = -123 ! [-] Turn on (1) or off (0) breaker delay model ! Roller parameters integer*4 :: roller = -123 ! [-] Turn on (1) or off(0) roller model real*8 :: beta = -123 ! [-] Breaker slope coefficient in roller model integer*4 :: rfb = -123 ! [-] If rfb = 1 then maximum wave surface slope is feeded back in roller energy balance; else rfb = par%Beta ! Wave-current interaction parameters integer*4 :: wci = -123 ! [-] Turns on (1) or off (0) wave-current interaction real*8 :: hwci = -123 ! [m] Minimum depth until which wave-current interaction is used real*8 :: cats = -123 ! [Trep] Current averaging time scale for wci, in terms of mean wave periods ! Flow parameters real*8 :: C = -123 ! [m^0.5s^-1] Chezy coefficient real*8 :: cf = -123 ! [-] Friction coefficient flow real*8 :: nuh = -123 ! [m^2s^-1] Horizontal background viscosity real*8 :: nuhfac = -123 ! [-] Viscosity switch for roller induced turbulent horizontal viscosity real*8 :: nuhv = -123 ! [-] Longshore viscosity enhancement factor, following Svendsen (?) integer*4 :: smag = -123 ! [-] (=1) Use smagorinsky subgrid model for viscocity ! Coriolis force parameters real*8 :: wearth = -123 ! [hour^-1] Angular velocity of earth calculated as: 1/rotation_time (in hours), later changed in calculation code to rad/s real*8 :: lat = -123 ! [deg] Latitude at model location for computing Coriolis ! Wind parameters real*8 :: rhoa = -123 ! [kgm^-3] Air density real*8 :: Cd = -123 ! [-] Wind drag coefficient real*8 :: windv = -123 ! [ms^-1] Wind velocity, in case of stationary wind real*8 :: windth = -123 ! [deg] Nautical wind direction, in case of stationary wind character(256) :: windfile = 'abc' ! Note, will replace lookup in readwind [name] Name of file with non-stationary wind data ! Groundwater parameters real*8 :: kx = -123 ! [ms^-1] Darcy-flow permeability coefficient in x-direction [m/s] real*8 :: ky = -123 ! [ms^-1] Darcy-flow permeability coefficient in y-direction [m/s] real*8 :: kz = -123 ! [ms^-1] Darcy-flow permeability coefficient in z-direction [m/s] real*8 :: dwetlayer = -123 ! [m] Thickness of the top soil layer interacting more freely with the surface water real*8 :: aquiferbot = -123 ! Note, will replace lookup in groundwater module [m] Level of uniform aquifer bottom character(256):: aquiferbotfile = 'abc' ! Note, will replace lookup in groundwater module [name] Name of the aquifer bottom file real*8 :: gw0 = -123 ! Note, will replace lookup in groundwater module [m] Level initial groundwater level character(256):: gw0file = 'abc' ! Note, will replace lookup in groundwater module [name] Name of initial groundwater level file ! Q3D sediment transport parameters real*8 :: vonkar = -123 ! von Karman constant real*8 :: vicmol = -123 ! molecular viscosity integer*4 :: kmax = -123 ! [-] Number of sigma layers in Quasi-3D model; kmax = 1 (default) is without vertical structure of flow and suspensions real*8 :: sigfac = -123 ! [-] dsig scales with log(sigfac). Default = 1.3 ! Non-hydrostatic correction parameters integer*4 :: solver_maxit = -123 ! [-] Maximum number of iterations in the linear SIP solver real*8 :: solver_acc = -123 ! [?] accuracy with respect to the right-hand side used ! in the following termination criterion: ! ||b-Ax || < acc*||b|| real*8 :: solver_urelax = -123 ! [?] Underrelaxation parameter integer*4 :: solver = -123 ! [-] Solver used to solve the linear system, 1=SIP, 2=TRIDIAG (only for 1d) real*8 :: kdmin = -123 ! Minimum value of kd ( pi/dx > minkd ) real*8 :: dispc = -123 ! Coefficient in front of the vertical pressure gradient, Default = 1. real*8 :: Topt = -123 ! Absolute period to optimize coefficient ! Bed composition parameters real*8 :: rhos = -123 ! [kgm^-3] Solid sediment density (no pores) integer*4 :: ngd = -123 ! [-] Number of sediment classes integer*4 :: nd = -123 ! [-] Number of computational layers in the bed real*8 :: dzg1 = -123 ! [m] Thickness of top sediment class layers real*8 :: dzg2 = -123 ! [m] Nominal thickness of variable sediment class layer real*8 :: dzg3 = -123 ! [m] Thickness of bottom sediment class layers real*8 :: por = -123 ! [-] Porosity real*8,dimension(99) :: D50 = -123 ! Note, will replace lookup in initialize [m] D50 grain size per grain type real*8,dimension(99) :: D90 = -123 ! Note, will replace lookup in initialize [m] D90 grain size per grain type real*8,dimension(99) :: sedcal = -123 ! Note, will replace lookup in initialize [-] Sediment transport calibration coefficient per grain type real*8,dimension(99) :: ucrcal = -123 ! Note, will replace lookup in initialize [-] Critical velocity calibration coefficient per grain type ! Sediment transport parameters character(24) :: waveform = 'abc' ! [-] Option for waveshape model: 1 = Ruessink & Van Rijn, 2 = Van Thiel de Vries, 2009 character(24) :: form = 'abc' ! [-] Equilibrium sed. conc. formulation: 1 = Soulsby van Rijn, 1997, 2 = Van Rijn 2008 with modifications by Van Thiel integer*4 :: sws = -123 ! [-] 1 = short wave & roller stirring and undertow, 0 = no short wave & roller stirring and undertow integer*4 :: lws = -123 ! [-] 1 = long wave stirring, 0 = no long wave stirring real*8 :: BRfac = -123 ! [-] Calibration factor surface slope real*8 :: facsl = -123 ! [-] Factor bedslope effect real*8 :: z0 = -123 ! [m] Zero flow velocity level in Soulsby van Rijn (1997) sed.conc. expression real*8 :: smax = -123 ! [-] Being tested: maximum Shields parameter for ceq Diane Foster real*8 :: tsfac = -123 ! [-] Coefficient determining Ts = tsfac * h/ws in sediment source term real*8 :: facua = -123 ! [-] Calibration factor time averaged flows due to wave asymmetry character(24) :: turb = 'abc' ! [-] Equlibrium sediment concentration is computed as function of: ! none = no turbulence, ! wave_averaged = wave averaged turbulence ! bore_averaged = maximum turbulence real*8 :: Tbfac = -123 ! [-] Calibration factor for bore interval Tbore: Tbore = Tbfac*Tbore real*8 :: Tsmin = -123 ! [s] Minimum adaptation time scale in advection diffusion equation sediment integer*4 :: lwt = -123 ! [-] Switch 0/1 long wave turbulence model real*8 :: betad = -123 ! [-] Dissipation parameter long wave breaking turbulence character(256) :: swtable = 'abc' ! [-] Name of intra short wave assymetry and skewness table integer*4 :: sus = -123 ! [-] Calibration factor for suspensions transports [0..1] integer*4 :: bed = -123 ! [-] Calibration factor for bed transports [0..1] ! Morphology parameters real*8 :: morfac = -123 ! [-] Morphological acceleration factor integer*4 :: morfacopt = -123 ! [-] Option indicating whether times should be adjusted (1) or not(0) for morfac real*8 :: morstart = -123 ! [s] Start time morphology, in morphological time real*8 :: wetslp = -123 ! [-] Critical avalanching slope under water (dz/dx and dz/dy) real*8 :: dryslp = -123 ! [-] Critical avalanching slope above water (dz/dx and dz/dy) real*8 :: hswitch = -123 ! [m] Water depth at which is switched from wetslp to dryslp real*8 :: dzmax = -123 ! [m/s/m] Maximum bedlevel change due to avalanching integer*4 :: struct = -123 ! [-] Switch for hard structures character(256):: ne_layer = 'abc' ! [name] Name of file containing depth of hard structure ! Output variables integer*4 :: timings = -123 ! [-] Switch to turn on (1) or off (0) progress output to screen real*8 :: tstart = -123 ! [s] Start time of output, in morphological time real*8 :: tint = -123 ! [s] Interval time of global output (replaced by tintg) real*8 :: tintg = -123 ! [s] Interval time of global output real*8 :: tintp = -123 ! [s] Interval time of point and runup gauge output real*8 :: tintc = -123 ! [s] Interval time of cross section output real*8 :: tintm = -123 ! [s] Interval time of mean,var,max,min output character(256):: tsglobal = 'abc' ! [name] Name of file containing timings of global output character(256):: tspoints = 'abc' ! [name] Name of file containing timings of point output character(256):: tscross = 'abc' ! [name] Name of file containing timings of cross section output character(256):: tsmean = 'abc' ! [name] Name of file containing timings of mean, max, min and var output integer*4 :: nglobalvar = -123 ! [-] Number of global output variables (as specified by user) character(len=maxnamelen), dimension(numvars) :: globalvars = 'abc' ! [-] Mnems of global output variables, not per se the same sice as nglobalvar (invalid variables, defaults) integer*4 :: nmeanvar = -123 ! [-] Number of mean,min,max,var output variables integer*4 :: npoints = -123 ! [-] Number of output point locations character(len=maxnamelen), dimension(numvars) :: pointvars = 'abc' ! [-] Mnems of point output variables (by variables) integer, dimension(:), pointer :: pointtypes => NULL() ! [-] Point types (0 = point, 1=rugauge) real*8 ,dimension(:), pointer :: xpointsw => NULL() ! world x-coordinate of output points real*8 ,dimension(:), pointer :: ypointsw => NULL() ! world y-coordinate of output points character(len=maxnamelen), dimension(numvars) :: meansvars = 'abc' ! [-] Mnems of mean output variables (by variables) integer*4 :: nrugauge = -123 ! [-] Number of output runup gauge locations integer*4 :: ncross = -123 ! [-] Number of output cross sections character(256):: outputformat = 'debug' ! [name] Choice of output file format: 'netcdf', 'fortran', or 'debug' ! Drifters parameters integer*4 :: ndrifter = -123 ! Note: will replace lookup in drifters module [-] Number of drifers character(256) :: drifterfile = 'abc' ! Note: will replace lookup in drifters module [name] Name of drifter data file ! Wave numerics parameters character(256):: scheme = 'abc' ! [-] Use first-order upwind (upwind_1), second order upwind (upwind_2) or Lax-Wendroff (lax_wendroff) ! for wave propagation real*8 :: wavint = -123 ! [s] Interval between wave module calls (only in stationary wave mode) real*8 :: maxerror = -123 ! [m] Maximum wave height error in wave stationary iteration integer*4 :: maxiter = -123 ! [-] Maximum number of iterations in wave stationary ! Flow numerics parameters real*8 :: eps = -123 ! [m] Threshold water depth above which cells are considered wet real*8 :: umin = -123 ! [m/s] Threshold velocity for upwind velocity detection and for vmag2 in eq. sediment concentration real*8 :: hmin = -123 ! [m] Threshold water depth above which Stokes drift is included integer*4 :: secorder = -123 ! [-] Use second order corrections to advection/non-linear terms based on mcCormack scheme ! Sediment transport numerics parameters real*8 :: thetanum = -123 ! [-] Coefficient determining whether upwind (1) or central scheme (0.5) is used. integer*4 :: sourcesink = -123 ! [-] In suspended transport use source-sink terms to calculate bed level change (1) or sus transport gradients (0) ! Bed update numerics parameters real*8 :: frac_dz = -123 ! [-] Relative thickness to split time step for bed updating integer*4 :: nd_var = -123 ! [-] Index of layer with variable thickness real*8 :: split = -123 ! [-] Split threshold for variable sediment layer (ratio to nominal thickness) real*8 :: merge = -123 ! [-] Merge threshold for variable sediment layer (ratio to nominal thickness) ! MPI parameters character(4) :: mpiboundary = 'abc' ! Fix mpi boundaries along y-lines ('y'), x-lines ('x'), or find shortest boundary ('auto') ! Constants, not read in params.txt real*8 :: px = -123 ! [-] Pi complex(kind(0.0d0)) :: compi = -123 ! [-] Imaginary unit real*8 :: rhog8 = -123 ! [Nm^-3] 1/8*rho*g ! Variables, not read in params.txt real*8 :: dt = -123 ! [s] Computational time step, in hydrodynamic time real*8 :: t = -123 ! [s] Computational time, in hydrodynamic time real*8 :: tnext = -123 ! [s] Next time point for output or wave stationary calculation, in hydrodynamic time real*8 :: Emean = -123 ! Note: can be made local [Jm^-2] Mean wave energy at boundary real*8 :: w = -123 ! Note: can be made local [ms^-1] Fall velocity sediment integer*4 :: listline = -123 ! Note: can be made local [-] Keeps track of the record line in bcf-files in case of instat 4,5,6 real*8 :: fc = -123 ! Note: can be made local [s^-1] Coriolis forcing parameter integer*4 :: tidelen = -123 ! [-] Length of input tidal time series integer*4 :: windlen = -123 ! [-] Length of input wind time series real*8 :: Llong = -123 ! Note: can be made local [m] Alongshore wave group length for case instat=1 real*8 :: zs01 = -123 ! [m] Initial water level first sea boundary real*8 :: zs02 = -123 ! [m] Initial water level second sea boundary real*8 :: zs03 = -123 ! [m] Initial water level first land boundary real*8 :: zs04 = -123 ! [m] Initial water level second land boundary real*8,dimension(2) :: xyzs01 = -123 ! global xy coordinates of corner (x=1,y=1) real*8,dimension(2) :: xyzs02 = -123 ! global xy coordinates of corner (x=1,y=N) real*8,dimension(2) :: xyzs03 = -123 ! global xy coordinates of corner (x=N,y=N) real*8,dimension(2) :: xyzs04 = -123 ! global xy coordinates of corner (x=N,y=1) real*8 :: waverr = -123 ! [-] max. absolute wave height difference between time steps end type parameters contains subroutine all_input(par) use readkey_module use xmpi_module use filefunctions use logging_module implicit none type(parameters) :: par character(128) :: testc character(len=256) :: line, keyword character(24),dimension(:),allocatable :: allowednames,oldnames integer :: filetype integer :: i, ic, nvars ! integers for reading output variables call writelog('sl','','Reading input parameters: ') ! Physical processes call writelog('l','','--------------------------------') call writelog('l','','Physical processes: ') par%swave = readkey_int ('params.txt','swave', 1, 0, 1) par%lwave = readkey_int ('params.txt','lwave', 1, 0, 1) par%flow = readkey_int ('params.txt','flow', 1, 0, 1) par%sedtrans = readkey_int ('params.txt','sedtrans', 1, 0, 1) par%morphology = readkey_int ('params.txt','morphology', 1, 0, 1) par%nonh = readkey_int ('params.txt','nonh', 0, 0, 1) par%gwflow = readkey_int ('params.txt','gwflow', 0, 0, 1) par%q3d = readkey_int ('params.txt','q3d', 0, 0, 1) ! Grid parameters call writelog('l','','--------------------------------') call writelog('l','','Grid parameters: ') par%xori = readkey_dbl('params.txt','xori', 0.d0, -1d9, 1d9) par%yori = readkey_dbl('params.txt','yori', 0.d0, -1d9, 1d9) par%alfa = readkey_dbl('params.txt','alfa', 0.d0, 0.d0, 360.d0) par%nx = readkey_int('params.txt','nx', 50, 2, 10000,required=.true.) par%ny = readkey_int('params.txt','ny', 2, 0, 10000,required=.true.) par%posdwn= readkey_dbl('params.txt','posdwn', 1.d0, -1.d0, 1.d0) par%depfile = readkey_name('params.txt','depfile',required=.true.) call check_file_exist(par%depfile) call check_file_length(par%depfile,par%nx+1,par%ny+1) par%vardx = readkey_int('params.txt','vardx', 0, 0, 1) if (par%vardx==0) then par%dx = readkey_dbl('params.txt','dx', -1.d0, 0.d0, 1d9,required=.true.) par%dy = readkey_dbl('params.txt','dy', -1.d0, 0.d0, 1d9,required=.true.) else par%xfile = readkey_name('params.txt','xfile') call check_file_exist(par%xfile) call check_file_length(par%xfile,par%nx+1,par%ny+1) par%yfile = readkey_name('params.txt','yfile') call check_file_exist(par%yfile) call check_file_length(par%yfile,par%nx+1,par%ny+1) endif par%thetamin = readkey_dbl ('params.txt','thetamin', -90.d0, -180.d0, 180.d0,required=(par%swave==1)) par%thetamax = readkey_dbl ('params.txt','thetamax', 90.d0, -180.d0, 180.d0,required=(par%swave==1)) par%dtheta = readkey_dbl ('params.txt','dtheta', 10.d0, 0.1d0, 20.d0,required=(par%swave==1)) par%thetanaut= readkey_int ('params.txt','thetanaut', 0, 0, 1) ! ! ! Model time parameters call writelog('l','','--------------------------------') call writelog('l','','Model time parameters: ') par%CFL = readkey_dbl ('params.txt','CFL', 0.7d0, 0.1d0, 0.9d0) par%tstop = readkey_dbl ('params.txt','tstop', 2000.d0, 1.d0, 1000000.d0,required=.true.) ! ! ! Physical constants call writelog('l','','--------------------------------') call writelog('l','','Physical constants: ') par%rho = readkey_dbl ('params.txt','rho', 1025.0d0, 1000.0d0, 1040.0d0) par%g = readkey_dbl ('params.txt','g', 9.81d0, 9.7d0, 9.9d0) ! ! ! ! ! Initial conditions call writelog('l','','--------------------------------') call writelog('l','','Initial conditions: ') par%zsinitfile = readkey_name('params.txt','zsinitfile') if (par%zsinitfile==' ') then ! do nothing else call check_file_exist(par%zsinitfile) call check_file_length(par%zsinitfile,par%nx+1,par%ny+1) endif ! ! ! Wave boundary condition parameters call writelog('l','','--------------------------------') call writelog('l','','Wave boundary condition parameters: ') allocate(allowednames(12),oldnames(12)) allowednames=(/'stat ','bichrom ','ts_1 ','ts_2 ','jons ','swan ', & 'vardens ','reuse ','nonh ','off ','stat_table ','jons_table '/) oldnames=(/'0 ','1 ','2 ','3 ','4 ','5 ','6 ','7 ','8 ','9 ','40','41'/) ! function = file key default n allowed n old allowed allowed names old allowed names par%instat = readkey_str('params.txt','instat','bichrom',12 ,12 ,allowednames ,oldnames,required=(par%swave==1)) deallocate(allowednames,oldnames) if ( trim(par%instat)=='jons' .or. & trim(par%instat)=='swan' .or. & trim(par%instat)=='vardens'.or. & trim(par%instat)=='stat_table' .or. & trim(par%instat)=='jons_table' & )then par%bcfile = readkey_name('params.txt','bcfile') call check_file_exist(par%bcfile) call checkbcfilelength(par%tstop,par%instat,par%bcfile,filetype) ! Only carried out on xmaster so: #ifdef USEMPI call xmpi_bcast(filetype) #endif else filetype=-1 endif par%taper = readkey_dbl ('params.txt','taper', 100.d0, 0.0d0, 1000.d0) if (trim(par%instat) == 'stat') then par%Hrms = readkey_dbl ('params.txt','Hrms', 1.d0, 0.d0, 10.d0) par%Tm01 = readkey_dbl ('params.txt','Tm01', 10.d0, 1.d0, 20.d0) par%Trep = readkey_dbl ('params.txt','Trep', par%Tm01, 1.d0, 20.d0) par%dir0 = readkey_dbl ('params.txt','dir0', 270.d0, 180.d0, 360.d0) par%m = readkey_int ('params.txt','m', 10, 2, 128) elseif (trim(par%instat) == 'bichrom') then par%Hrms = readkey_dbl ('params.txt','Hrms', 1.d0, 0.d0, 10.d0) par%Tm01 = readkey_dbl ('params.txt','Tm01', 10.d0, 1.d0, 20.d0) par%Trep = readkey_dbl ('params.txt','Trep', par%Tm01, 1.d0, 20.d0) par%Tlong = readkey_dbl ('params.txt','Tlong', 80.d0, 20.d0, 300.d0) par%dir0 = readkey_dbl ('params.txt','dir0', 270.d0, 180.d0, 360.d0) par%m = readkey_int ('params.txt','m', 10, 2, 128) elseif (trim(par%instat) == 'ts_1' .or. trim(par%instat) == 'ts_2') then par%Hrms = readkey_dbl ('params.txt','Hrms', 1.d0, 0.d0, 10.d0) par%Tm01 = readkey_dbl ('params.txt','Tm01', 10.d0, 1.d0, 20.d0) par%Trep = readkey_dbl ('params.txt','Trep', par%Tm01, 1.d0, 20.d0) par%dir0 = readkey_dbl ('params.txt','dir0', 270.d0, 180.d0, 360.d0) par%m = readkey_int ('params.txt','m', 10, 2, 128) endif allocate(allowednames(2),oldnames(2)) allowednames=(/'neumann ','wavecrest'/) oldnames=(/'0','1'/) par%leftwave = readkey_str('params.txt','leftwave','neumann',2,2,allowednames,oldnames) par%rightwave = readkey_str('params.txt','rightwave','neumann',2,2,allowednames,oldnames) deallocate(allowednames,oldnames) ! ! ! Wave-spectrum boundary condition parameters if (trim(par%instat) == 'jons' .or. & trim(par%instat) == 'swan' .or. & trim(par%instat) == 'vardens'.or. & trim(par%instat) == 'jons_table' & )then call writelog('l','','--------------------------------') call writelog('l','','Wave-spectrum boundary condition parameters: ') par%random = readkey_int ('params.txt','random', 1, 0, 1) par%fcutoff = readkey_dbl ('params.txt','fcutoff', 0.d0, 0.d0, 40.d0) par%nspr = readkey_int ('params.txt','nspr', 0, 0, 1) par%trepfac = readkey_dbl ('params.txt','trepfac', 0.01d0, 0.d0, 1.d0) par%sprdthr = readkey_dbl ('params.txt','sprdthr', 0.08d0, 0.d0, 1.d0) par%oldwbc = readkey_int ('params.txt','oldwbc', 0, 0, 1) if (filetype==0) then par%rt = readkey_dbl('params.txt','rt' , 3600.d0, 1200.d0,7200.d0) par%dtbc = readkey_dbl('params.txt','dtbc', 0.5d0,0.1d0,2.0d0) endif if (trim(par%instat)=='swan') then par%dthetaS_XB = readkey_dbl ('params.txt','dthetaS_XB', 0.0d0, -360.0d0, 360.0d0) endif endif ! ! ! Flow boundary condition parameters ! front call writelog('l','','--------------------------------') call writelog('l','','Flow boundary condition parameters: ') allocate(allowednames(5),oldnames(5)) allowednames=(/'abs_1d ','abs_2d ','wall ','wlevel ','nonh_1d'/) oldnames=(/'0','1','2','3','4'/) par%front = readkey_str('params.txt','front','abs_2d',5,5,allowednames,oldnames) deallocate(allowednames,oldnames) ! left and right allocate(allowednames(2),oldnames(2)) allowednames=(/'neumann','wall '/) oldnames=(/'0','1'/) par%left = readkey_str('params.txt','left','neumann',2,2,allowednames,oldnames) par%right = readkey_str('params.txt','right','neumann',2,2,allowednames,oldnames) deallocate(allowednames,oldnames) ! Back allocate(allowednames(4),oldnames(4)) allowednames=(/'wall ', 'abs_1d ','abs_2d ','wlevel '/) oldnames=(/'0','1','2','3'/) par%back = readkey_str('params.txt','back','abs_2d',4,4,allowednames,oldnames) deallocate(allowednames,oldnames) ! Others par%ARC = readkey_int ('params.txt','ARC', 1, 0, 1) par%order = readkey_dbl ('params.txt','order', 2.d0, 1.d0, 2.d0) par%carspan = readkey_int ('params.txt','carspan', 0, 0, 1) par%epsi = readkey_dbl ('params.txt','epsi', 0.d0, 0.d0, 0.2d0) ! ! ! Tide boundary conditions call writelog('l','','--------------------------------') call writelog('l','','Tide boundary conditions: ') par%tideloc = readkey_int ('params.txt','tideloc', 0, 0, 4) if (par%tideloc>0) then if (par%tideloc==2) then par%paulrevere = readkey_int ('params.txt','paulrevere', 0, 0, 1) endif par%zs0file = readkey_name('params.txt','zs0file') call check_file_exist(par%zs0file) else par%zs0 = readkey_dbl ('params.txt','zs0', 0.0d0, -5.d0, 5.d0) endif ! ! ! Discharge boundary conditions call writelog('l','','--------------------------------') call writelog('l','','Discharge boundary conditions: ') par%disch_loc_file = readkey_name('params.txt','disch_loc_file') if (par%disch_loc_file==' ') then ! do nothing else call check_file_exist(par%disch_loc_file) endif par%disch_timeseries_file = readkey_name('params.txt','disch_timeseries_file') if (par%disch_timeseries_file==' ') then ! do nothing else call check_file_exist(par%disch_timeseries_file) endif ! ! ! Wave breaking parameters call writelog('l','','--------------------------------') call writelog('l','','Wave breaking parameters: ') allocate(allowednames(4),oldnames(4)) allowednames=(/'roelvink1 ','baldock ','roelvink2 ','roelvink_daly'/) oldnames =(/'1','2','3','4'/) par%break = readkey_str('params.txt','break','roelvink2',4,4,allowednames,oldnames) deallocate(allowednames,oldnames) par%gamma = readkey_dbl ('params.txt','gamma', 0.55d0, 0.4d0, 0.9d0) !changed 28/11 if (trim(par%break)=='roelvink_daly') par%gamma2 = readkey_dbl ('params.txt','gamma2', 0.3d0, 0.0d0, 0.5d0) par%alpha = readkey_dbl ('params.txt','alpha', 1.0d0, 0.5d0, 2.0d0) par%n = readkey_dbl ('params.txt','n', 10.0d0, 5.0d0, 20.0d0) !changed 28/11 par%gammax = readkey_dbl ('params.txt','gammax', 2.d0, .4d0, 5.d0) !changed 28/11 par%delta = readkey_dbl ('params.txt','delta', 0.0d0, 0.0d0, 1.0d0) par%fw = readkey_dbl ('params.txt','fw', 0.d0, 0d0, 1.0d0) par%breakerdelay = readkey_int ('params.txt','breakerdelay', 1, 0, 1) ! ! ! Roller parameters call writelog('l','','--------------------------------') call writelog('l','','Roller parameters: ') par%roller = readkey_int ('params.txt','roller', 1, 0, 1) par%beta = readkey_dbl ('params.txt','beta', 0.15d0, 0.05d0, 0.3d0) par%rfb = readkey_int ('params.txt','rfb', 1, 0, 1) ! ! Wave-current interaction parameters call writelog('l','','--------------------------------') call writelog('l','','Wave-current interaction parameters: ') par%wci = readkey_int ('params.txt','wci', 0, 0, 1) par%hwci = readkey_dbl ('params.txt','hwci', 0.1d0, 0.001d0, 1.d0) par%cats = readkey_dbl ('params.txt','cats', 20.d0, 1.d0, 50.d0) ! ! ! Flow parameters call writelog('l','','--------------------------------') call writelog('l','','Flow parameters: ') par%cf = readkey_dbl ('params.txt','cf', 3.d-3, 0.d0, 0.1d0) par%C = readkey_dbl ('params.txt','C', sqrt(par%g/par%cf), 20.d0, 100.d0) par%nuh = readkey_dbl ('params.txt','nuh', 0.15d0, 0.0d0, 1.0d0) par%nuhfac = readkey_dbl ('params.txt','nuhfac', 0.0d0, 0.0d0, 1.0d0) par%nuhv = readkey_dbl ('params.txt','nuhv', 1.d0, 1.d0, 20.d0) par%smag = readkey_int ('params.txt','smag', 1, 0, 1) ! ! ! Coriolis force parameters call writelog('l','','--------------------------------') call writelog('l','','Coriolis force parameters: ') par%wearth = readkey_dbl ('params.txt','wearth', 1.d0/24.d0, 0.d0, 1.d0) par%lat = readkey_dbl ('params.txt','lat', 0.d0, -90.d0, 90.d0) ! ! ! Wind parameters call writelog('l','','--------------------------------') call writelog('l','','Wind parameters: ') par%rhoa = readkey_dbl ('params.txt','rhoa', 1.25d0, 1.0d0, 2.0d0) par%Cd = readkey_dbl ('params.txt','Cd', 0.002d0, 0.0001d0, 0.01d0) par%windfile = readkey_name('params.txt','windfile') if (par%windfile==' ') then par%windv = readkey_dbl ('params.txt','windv', 0.0d0, 0.0d0, 200.0d0) par%windth = readkey_dbl ('params.txt','windth', 270.0d0, -360.0d0, 360.0d0) else call check_file_exist(par%windfile) endif ! ! ! Groundwater parameters if (par%gwflow==1) then call writelog('l','','--------------------------------') call writelog('l','','Groundwater parameters: ') par%kx = readkey_dbl ('params.txt','kx' , 0.0001d0 , 0.00001d0, 0.01d0) par%ky = readkey_dbl ('params.txt','ky' , par%kx , 0.00001d0, 0.01d0) par%kz = readkey_dbl ('params.txt','kz' , par%kx , 0.00001d0, 0.01d0) par%dwetlayer = readkey_dbl ('params.txt','dwetlayer' , 0.2d0 , 0.01d0 , 1.d0) par%aquiferbotfile = readkey_name('params.txt','aquiferbotfile') if (par%aquiferbotfile==' ') then par%aquiferbot = readkey_dbl('params.txt','aquiferbot',-10.d0,-100.d0,100.d0) !also read in groundwater.f90 which determines value else call check_file_exist(par%aquiferbotfile) endif par%gw0file = readkey_name('params.txt','gw0file') if (par%gw0file==' ') then par%gw0 = readkey_dbl('params.txt','gw0',0.d0,-5.d0,5.d0) else call check_file_exist(par%gw0file) endif endif ! ! ! Q3D sediment transport parameters if (par%q3d==1) then call writelog('l','','--------------------------------') call writelog('l','','Q3D sediment transport parameters: ') par%vonkar = readkey_dbl ('params.txt','vonkar', 0.4d0, 0.01d0, 1.d0) par%vicmol = readkey_dbl ('params.txt','vicmol', 0.000001d0, 0.d0, 0.001d0) par%kmax = readkey_int ('params.txt','kmax ', 1, 1, 1000) par%sigfac = readkey_dbl ('params.txt','sigfac ',1.3d0, 0.00d0, 10.d0) endif ! ! ! Non-hydrostatic correction parameters if (par%nonh==1) then call writelog('l','','--------------------------------') call writelog('l','','Non-hydrostatic correction parameters: ') par%solver_maxit = readkey_int('params.txt','solver_maxit' ,30,1,1000) par%solver_acc = readkey_dbl('params.txt','solver_acc' ,0.005d0,0.00001d0,0.1d0) par%solver_urelax= readkey_dbl('params.txt','solver_urelax' ,0.92d0,0.5d0,0.99d0) par%solver = readkey_int('params.txt','solver' ,1,0,2) par%kdmin = readkey_dbl('params.txt','kdmin' ,0.0d0,0.0d0,0.05d0) par%dispc = readkey_dbl('params.txt','dispc' ,1.0d0,0.1d0,2.0d0) par%Topt = readkey_dbl('params.txt','Topt', 10.d0, 1.d0, 20.d0) endif ! ! ! Bed composition parameters par%ngd = readkey_int ('params.txt','ngd', 1, 1, 20) par%nd = readkey_int ('params.txt','nd ', 3, 3, 1000) if (par%sedtrans==1) then call writelog('l','','--------------------------------') call writelog('l','','Bed composition parameters: ') par%rhos = readkey_dbl ('params.txt','rhos', 2650d0, 2400.d0, 2800.d0) par%dzg1 = readkey_dbl ('params.txt','dzg', 0.1d0, 0.01d0, 1.d0) par%dzg1 = readkey_dbl ('params.txt','dzg1', par%dzg1, 0.01d0, 1.d0) par%dzg2 = readkey_dbl ('params.txt','dzg2', par%dzg1, 0.01d0, 1.d0) par%dzg3 = readkey_dbl ('params.txt','dzg3', par%dzg1, 0.01d0, 1.d0) par%por = readkey_dbl ('params.txt','por', 0.4d0, 0.3d0, 0.5d0) testc = readkey_name('params.txt','D50') if (testc==' ') then par%D50(1:par%ngd)=0.0002d0 ! Default call writelog('l','','Setting D50 to default value 0.0002') else read(testc,*) par%D50(1:par%ngd) endif testc = readkey_name('params.txt','D90') if (testc==' ') then par%D90(1:par%ngd)=0.0003d0 ! Default call writelog('l','','Setting D90 to default value 0.0003') else read(testc,*) par%D50(1:par%ngd) endif testc = readkey_name('params.txt','sedcal') if (testc==' ') then par%sedcal(1:par%ngd)=1.d0 ! Default call writelog('l','','Setting sedcal to default value 1.0') else read(testc,*) par%sedcal(1:par%ngd) endif testc = readkey_name('params.txt','ucrcal') if (testc==' ') then par%ucrcal(1:par%ngd)=1.d0 ! Default call writelog('l','','Setting ucrcal to default value 1.0') else read(testc,*) par%ucrcal(1:par%ngd) endif endif ! ! ! Sediment transport parameters if (par%sedtrans==1) then call writelog('l','','--------------------------------') call writelog('l','','Sediment transport parameters: ') allocate(allowednames(2),oldnames(2)) allowednames=(/'soulsby_vanrijn ','vanthiel_vanrijn'/) oldnames=(/'1','2'/) par%form = readkey_str('params.txt','form','soulsby_vanrijn',2,2,allowednames,oldnames) deallocate(allowednames,oldnames) allocate(allowednames(2),oldnames(2)) allowednames=(/'ruessink_vanrijn','vanthiel '/) oldnames=(/'1','2'/) par%waveform = readkey_str('params.txt','waveform','vanthiel',2,2,allowednames,oldnames) deallocate(allowednames,oldnames) par%sws = readkey_int ('params.txt','sws', 1, 0, 1) par%lws = readkey_int ('params.txt','lws', 1, 0, 1) par%BRfac = readkey_dbl ('params.txt','BRfac', 1.0d0, 0.d0, 1.d0) par%facsl = readkey_dbl ('params.txt','facsl ',0.00d0, 0.00d0, 1.6d0) par%z0 = readkey_dbl ('params.txt','z0 ',0.006d0, 0.0001d0, 0.05d0) par%smax = readkey_dbl ('params.txt','smax', -1.d0, -1.d0, 3.d0) !changed 28/11 and back 10/2 par%tsfac = readkey_dbl ('params.txt','tsfac', 0.1d0, 0.01d0, 1.d0) par%facua = readkey_dbl ('params.txt','facua ',0.00d0, 0.00d0, 1.0d0) allocate(allowednames(3),oldnames(3)) allowednames=(/'none ','wave_averaged','bore_averaged'/) oldnames=(/'0','1','2'/) par%turb = readkey_str('params.txt','turb','bore_averaged',3,3,allowednames,oldnames) deallocate(allowednames,oldnames) par%Tbfac = readkey_dbl ('params.txt','Tbfac ',1.0d0, 0.00d0, 1.0d0) par%Tsmin = readkey_dbl ('params.txt','Tsmin ',0.2d0, 0.01d0, 10.d0) par%lwt = readkey_int ('params.txt','lwt ',0, 0, 1) par%betad = readkey_dbl ('params.txt','betad ',1.0d0, 0.00d0, 10.0d0) if (trim(par%form)=='vanthiel_vanrijn') then par%swtable = readkey_name('params.txt','swtable',required=.true.) call check_file_exist(par%swtable) call check_file_length(par%swtable,18,33,40) endif par%sus = readkey_int ('params.txt','sus ',1, 0, 1) par%bed = readkey_int ('params.txt','bed ',1, 0, 1) endif ! ! ! Morphology parameters if (par%morphology==1) then call writelog('l','','--------------------------------') call writelog('l','','Morphology parameters: ') par%morfac = readkey_dbl ('params.txt','morfac', 0.0d0, 0.d0, 1000.d0) par%morfacopt= readkey_int ('params.txt','morfacopt', 1, 0, 1) par%morstart = readkey_dbl ('params.txt','morstart',120.d0, 0.d0, 10000.d0) par%wetslp = readkey_dbl ('params.txt','wetslp', 0.3d0, 0.1d0, 1.d0) par%dryslp = readkey_dbl ('params.txt','dryslp', 1.0d0, 0.1d0, 2.d0) par%hswitch = readkey_dbl ('params.txt','hswitch',0.1d0, 0.01d0, 1.0d0) par%dzmax = readkey_dbl ('params.txt','dzmax ',0.05d0, 0.00d0, 1.0d0) par%struct = readkey_int ('params.txt','struct ',0 , 0, 1) if (par%struct==1) then par%ne_layer = readkey_name('params.txt','ne_layer') call check_file_exist(par%ne_layer) call check_file_length(par%ne_layer,par%nx+1,par%ny+1) endif endif ! ! ! Output variables call writelog('l','','--------------------------------') call writelog('l','','Output variables: ') par%timings = readkey_int ('params.txt','timings', 1, 0, 1) testc = readkey_name('params.txt','tunits') if (len(trim(testc)) .gt. 0) par%tunits = trim(testc) par%tstart = readkey_dbl ('params.txt','tstart', 1.d0, 0.d0,1000000.d0) par%tint = readkey_dbl ('params.txt','tint', 1.d0, .01d0, 100000.d0) ! Robert par%tsglobal = readkey_name('params.txt','tsglobal') if (par%tsglobal==' ') then par%tintg = readkey_dbl ('params.txt','tintg', par%tint, .01d0, 100000.d0) ! Robert endif par%tspoints = readkey_name('params.txt','tspoints') if (par%tspoints==' ') then par%tintp = readkey_dbl ('params.txt','tintp', par%tint, .01d0, 100000.d0) ! Robert endif par%tscross = readkey_name('params.txt','tscross') if (par%tscross==' ') then par%tintc = readkey_dbl ('params.txt','tintc', par%tint, .01d0, 100000.d0) ! Robert endif par%tsmean = readkey_name('params.txt','tsmean') if (par%tsmean==' ') then par%tintm = readkey_dbl ('params.txt','tintm', par%tstop-par%tstart, 1.d0, par%tstop-par%tstart) ! Robert endif !!!!! GLOBAL OUTPUT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! par%nglobalvar = readkey_int ('params.txt','nglobalvar', -1, -1, 20) call readglobalvars(par) par%npoints = readkey_int ('params.txt','npoints', 0, 0, 50) par%nrugauge = readkey_int ('params.txt','nrugauge', 0, 0, 50) ! update the pointvariables !call readpointvars(par, par%xpointsw, par%ypointsw, par%pointtypes, par%pointvars) ! Robert: to deal with MPI some changes here call readpointvars(par) ! par%nmeanvar = readkey_int ('params.txt','nmeanvar' , 0, 0, 15) call readmeans(par) par%ncross = readkey_int ('params.txt','ncross', 0, 0, 50) allocate(allowednames(3)) allocate(oldnames(0)) allowednames = (/'fortran', 'netcdf ', 'debug '/) par%outputformat= readkey_str ('params.txt','outputformat','fortran',3, 0, allowednames ,oldnames,required=.false.) deallocate(allowednames) deallocate(oldnames) ! ! ! Drifters parameters call writelog('l','','--------------------------------') call writelog('l','','Drifters parameters: ') par%ndrifter = readkey_int ('params.txt','ndrifter', 0, 0, 50) if (par%ndrifter>0) then par%drifterfile = readkey_name('params.txt','drifterfile') call check_file_exist(par%drifterfile) endif ! ! ! Wave numerics parameters call writelog('l','','--------------------------------') call writelog('l','','Wave numerics parameters: ') allocate(allowednames(3)) allocate(oldnames(3)) allowednames = (/'upwind_1 ', 'lax_wendroff', 'upwind_2 '/) oldnames = (/'1', '2', '3'/) par%scheme= readkey_str ('params.txt','scheme','lax_wendroff',3, 3, allowednames ,oldnames) deallocate(allowednames) deallocate(oldnames) if (trim(par%instat) == 'stat' .or. trim(par%instat) == 'stat_table') then par%wavint = readkey_dbl ('params.txt','wavint', 60.d0, 1.d0, 3600.d0) par%maxerror = readkey_dbl ('params.txt','maxerror', 0.00005d0, 0.00001d0, 0.001d0) par%maxiter = readkey_int ('params.txt','maxiter', 500, 2, 1000) endif ! ! ! Flow numerics parameters call writelog('l','','--------------------------------') call writelog('l','','Flow numerics parameters: ') par%eps = readkey_dbl ('params.txt','eps', 0.005d0, 0.001d0, 0.1d0) par%umin = readkey_dbl ('params.txt','umin', 0.0d0, 0.0d0, 0.2d0) par%hmin = readkey_dbl ('params.txt','hmin', 0.05d0, 0.001d0, 1.d0) par%secorder = readkey_int('params.txt','secorder' ,0,0,1) ! ! ! Sediment transport numerics parameters if (par%sedtrans==1) then call writelog('l','','--------------------------------') call writelog('l','','Sediment transport numerics parameters: ') par%thetanum = readkey_dbl ('params.txt','thetanum', 1.d0, 0.5d0, 1.d0) par%sourcesink = readkey_int ('params.txt','sourcesink ',0, 0, 1) endif ! ! ! Bed update numerics parameters if (par%morphology==1) then call writelog('l','','--------------------------------') call writelog('l','','Bed update numerics parameters: ') par%frac_dz = readkey_dbl ('params.txt','frac_dz', 0.7d0, 0.5d0, 0.98d0) par%nd_var = readkey_int ('params.txt','nd_var', 2, 2, par%nd) par%split = readkey_dbl ('params.txt','split', 1.01d0, 1.005d0, 1.10d0) par%merge = readkey_dbl ('params.txt','merge', 0.01d0, 0.005d0, 0.10d0) endif ! ! ! MPI parameters call writelog('l','','--------------------------------') call writelog('l','','MPI parameters: ') par%mpiboundary = readkey_name('params.txt','mpiboundary') if (par%mpiboundary==' ') then par%mpiboundary='auto' ! Default endif ! ! ! Finish call writelog('l','','--------------------------------') call writelog('sl','','Finished reading input parameters') call writelog('l','','--------------------------------') ! ! ! ------------------- Post-input processing ------------------------- ! ! ! Constants par%px = 4.d0*atan(1.d0) par%compi = (0.0d0,1.0d0) par%rhog8 = 1.0d0/8.0d0*par%rho*par%g ! ! ! Grid orientation if (par%alfa.lt.0) then par%alfa = 360.d0+par%alfa endif par%alfa = par%alfa*atan(1.0d0)/45.d0 ! All input converted directly to Cartesian XBeach grid direction if (par%posdwn<0.1d0) then par%posdwn=-1.d0 ! Backward compatibility, now posdwn = 0 also works in input (i.e. posdwn = false) endif ! ! ! Stop useless physical processes if (par%sedtrans==0 .and. par%morphology==1) then call writelog('lse','(a)','Error: morphology cannot be computed without sediment transport') call writelog('lse','(a)','Set sedtrans=1 or morphology=0') call writelog('lse','(a)','Stopping calculation') call halt_program endif ! ! ! Set taper to non-zero par%taper = max(par%taper,1.d-6) ! ! ! Set tide par%zs01=par%zs0 ! ! ! Compute Coriolis par%lat = par%lat*par%px/180.d0 par%wearth = par%px*par%wearth/1800.d0 par%fc = 2.d0*par%wearth*sin(par%lat) ! ! ! Only allow Baldock in stationary mode and Roelvink in non-stationary if (trim(par%instat) == 'stat' .or. trim(par%instat) == 'stat_table') then if (trim(par%break) .ne. 'baldock') then call writelog('ls','','Warning: Roelvink formulations not allowed in stationary calculation, use Baldock formulation.') endif else if (trim(par%break)=='baldock') then call writelog('ls','','Warning: Baldock formulation not allowed in non-stationary calculation, use a Roelvink formulation.') endif endif ! ! ! Convert cf from C par%cf = par%g/par%C**2 ! ! ! Set smax to huge if default is specified if (par%smax<0) par%smax=huge(0.d0) ! ! ! fix tint par%tint = min(par%tintg,par%tintp,par%tintm,par%tintc) ! Robert ! ! ! Source-sink check if (par%morfac>1.d0) then if (par%sourcesink==1) then call writelog('ls','','Warning: Using source-sink terms for bed level change with morfac can lead to') call writelog('ls','','loss of sediment mass conservation.') endif endif ! ! ! If using tide, epsi should be on if (par%tideloc>0) then if (par%epsi NULL() ! now we can allocate the memory again if (.not. xmaster) allocate(par%pointtypes(npoints)) ! and store the values from the local array if (.not. xmaster) par%pointtypes = pointtypes ! and clean up the local one. deallocate(pointtypes) if (xmaster) npoints = size(par%xpointsw) ! send it over call xmpi_bcast(npoints) ! now on all nodes allocate a array outside the par structure allocate(xpointsw(npoints)) ! Par is only filled on the master, so use that one and put it in the seperate array if (xmaster) xpointsw = par%xpointsw ! Now for another ugly step, we can't broadcast the whole array but have to do it per variable. call xmpi_bcast(xpointsw) ! so now everyone has the xpointsw, let's put it back in par ! first dereference the old one, otherwise we get a nasty error on ifort.... if (.not. xmaster) par%xpointsw => NULL() ! now we can allocate the memory again if (.not. xmaster) allocate(par%xpointsw(npoints)) ! and store the values from the local array if (.not. xmaster) par%xpointsw = xpointsw ! and clean up the local one. deallocate(xpointsw) if (xmaster) npoints = size(par%ypointsw) ! send it over call xmpi_bcast(npoints) ! now on all nodes allocate a array outside the par structure allocate(ypointsw(npoints)) ! Par is only filled on the master, so use that one and put it in the seperate array if (xmaster) ypointsw = par%ypointsw ! Now for another ugly step, we can't broadcast the whole array but have to do it per variable. call xmpi_bcast(ypointsw) ! so now everyone has the ypointsw, let's put it back in par ! first dereference the old one, otherwise we get a nasty error on ifort.... if (.not. xmaster) par%ypointsw => NULL() ! now we can allocate the memory again if (.not. xmaster) allocate(par%ypointsw(npoints)) ! and store the values from the local array if (.not. xmaster) par%ypointsw = ypointsw ! and clean up the local one. deallocate(ypointsw) return ! so, the following code is NOT used anymore. I left this here ! maybe method above does not work everywhere. wwvv ! For efficiency, this subroutine should use MPI_Pack and ! MPI_Unpack. However, since this subroutine is only called a ! few times, a more simple approach is used. ! !call xmpi_bcast(par%px) !call xmpi_bcast(par%Hrms) !.... end subroutine distribute_par #endif ! ! printparams for debugging only ! subroutine printparams(par,str) #ifdef USEMPI use xmpi_module #endif type (parameters) :: par character(*), intent(in) :: str integer :: f integer :: id id = 0 #ifdef USEMPI id = xmpi_rank #endif f = 70+id if (id .gt. 0) then return endif ! write(f,*) 'printparams ', id, ' ', str ! write(f,*) 'printpar ',id,' ','px:',par%px ! .... end subroutine printparams ! ! Some extra functions to make reading the output variables possible ! subroutine add_outmnem(mnem, vars) use logging_module implicit none character(len=maxnamelen), intent(in) :: mnem character(len=maxnamelen), dimension(:), allocatable, intent(inout) :: vars character(len=maxnamelen), dimension(:), allocatable :: temp integer :: i, nvars i = chartoindex(mnem) if (i .lt. 1) then if(xmaster) then call writelog('els','','Error: cannot locate variable "',trim(mnem),'". Program terminating') call halt_program endif return endif ! make a copy of the old vars allocate(temp(size(vars))) ! now copy the values temp =vars ! create room for an extra variable deallocate(vars) allocate(vars(size(temp)+1)) ! copy back vars(1:size(temp)) = temp ! and add the new one vars(size(vars)) = mnem ! clean up temp deallocate(temp) if(xmaster) then call writelog('ls','','Will generate global output for variable "',trim(mnem),'"') endif return end subroutine add_outmnem subroutine readglobalvars(par) use logging_module use mnemmodule implicit none type(parameters), intent(inout) :: par character(len=maxnamelen), dimension(:), allocatable :: globalvars character(len=256) :: line, keyword integer :: id, ic, i if (xmaster) then if (par%nglobalvar == -1) then allocate(globalvars(21)) globalvars = (/'H ', 'zs ', 'zs0 ', 'zb ', 'hh ', 'u ', 'v ', 'ue ', 've ', 'urms ', 'Fx ', & 'Fy ', 'ccg ', 'ceqsg', 'ceqbg', 'Susg ', 'Svsg ', 'E ', 'R ', 'D ', 'DR ' /) elseif (par%nglobalvar == 999) then ! Output all allocate(globalvars(size(mnemonics))) globalvars = mnemonics else ! User specified output ! This is important because there might be less usable variables than specified. ! Look for keyword nglobalvar in params.txt id=0 ! we allocate to 0, we'll let it grow inside add_mnem.... allocate(globalvars(0)) open(10,file='params.txt') do while (id == 0) read(10,'(a)')line ic=scan(line,'=') if (ic>0) then keyword=adjustl(line(1:ic-1)) if (keyword == 'nglobalvar') id=1 endif enddo ! Read through the variables lines, do i=1,par%nglobalvar read(10,'(a)')line line = trim(line) ! store the mnemonic in globalvars call add_outmnem(line, globalvars) end do close(10) end if ! globalvar par%globalvars(1:par%nglobalvar) = globalvars end if ! xmaster end subroutine ! ! FB: ! Now for a long one, reading the point and rugauges output options ! Just moved this from varoutput, so it can be used in ncoutput also ! Keeping as much as possible local to the subroutine... ! It can be reduced a bit by combining points and rugauges ! also instead of nvarpoints it can use a ragged array: ! For example: http://coding.derkeiler.com/Archive/Fortran/comp.lang.fortran/2004-05/0774.html ! ! The outputformat for points is xworld, yworld, nvars, var1#var2# ! Split for nrugauges and npoints ! This is a bit inconvenient because you have different sets of points ! I think it's more logical to get value per variable than per point.... ! So I read the data as follows: ! make a collection of all points ! store the types per point ! store all found variables in a combined list (not per point) ! So for each point all variables are stored ! subroutine readpointvars(par) use logging_module use mnemmodule implicit none type(parameters), intent(inout) :: par real*8,dimension(:),allocatable :: xpointsw,ypointsw integer, dimension(:), allocatable :: pointtypes ! [-] Point types (0 = point, 1=rugauge) character(len=maxnamelen), dimension(:), allocatable :: pointvars integer*4,dimension(:),allocatable :: nvarpoint ! vector with number of output variable per output point character(len=maxnamelen), dimension(:),allocatable :: temparray character(80) :: line, keyword character(len=maxnamelen) :: var logical :: varfound integer :: i, ic, icold, id, ii, index, nvars, ivar ! These must be allocated by all processes allocate(pointtypes(par%npoints+par%nrugauge)) allocate(xpointsw(par%npoints+par%nrugauge)) allocate(ypointsw(par%npoints+par%nrugauge)) nvars = 0 ! the total number of variables ! The rest can all be done by xmaster and then distributed by distribute_par if (xmaster) then allocate(nvarpoint(par%npoints+par%nrugauge)) allocate(temparray(99)) ! set the point types to the indicator pointtypes(1:par%npoints)=0 pointtypes(par%npoints+1:par%npoints+par%nrugauge)=1 temparray='' varfound = .false. if (par%npoints>0) then id=0 ! Look for keyword npoints in params.txt open(10,file='params.txt') do while (id == 0) read(10,'(a)')line ic=scan(line,'=') if (ic>0) then keyword=adjustl(line(1:ic-1)) if (keyword == 'npoints') id=1 endif enddo do i=1,par%npoints call writelog('ls','(a,i0)',' Output point ',i) read(10,*) xpointsw(i),ypointsw(i),nvarpoint(i),line call writelog('ls','(a,f0.2,a,f0.2)',' xpoint: ',xpointsw(i),' ypoint: ',ypointsw(i)) icold=0 do ii =1,nvarpoint(i) ic=scan(line(icold+1:80),'#') ic=ic+icold var=line(icold+1:ic-1) index = chartoindex(var) if (index/=-1) then ! ! see if we already found this variable.... ! varfound = .false. do ivar=1,nvars if (trim(temparray(ivar)) == trim(var)) then varfound = .true. end if end do ! we have a new variable, store it if (varfound .eqv. .false.) then nvars = nvars + 1 temparray(nvars)=trim(var) end if call writelog('sl','',' Output type: ''',trim(var),'''') else call writelog('sle','',' Unknown point output type: ''',trim(var),'''') call halt_program endif icold=ic enddo enddo close(10) endif ! par%npoints>0 if (par%nrugauge>0) then id=0 ! Look for keyword nrugauge in params.txt open(10,file='params.txt') do while (id == 0) read(10,'(a)')line ic=scan(line,'=') if (ic>0) then keyword=adjustl(line(1:ic-1)) if (keyword == 'nrugauge') id=1 endif enddo do i=1+par%npoints,par%nrugauge+par%npoints read(10,*) xpointsw(i),ypointsw(i),nvarpoint(i),line call writelog('ls','(a,i0)',' Output runup gauge ',i-par%npoints) call writelog('ls','(a,f0.2,a,f0.2)',' xpoint: ',xpointsw(i),' ypoint: ',ypointsw(i)) icold=0 do ii =1,nvarpoint(i) ic=scan(line(icold+1:80),'#') ic=ic+icold var=line(icold+1:ic-1) index = chartoindex(var) if (index/=-1) then ! ! see if we already found this variable.... ! varfound = .false. do ivar=1,nvars if (trim(temparray(ivar)) == trim(var)) then varfound = .true. end if end do ! we have a new variable, store it if (varfound .eqv. .false.) then nvars = nvars + 1 temparray(nvars)=trim(var) end if call writelog('sl','',' Output type: ''',trim(var),'''') else call writelog('sle','',' Unknown point output type: ''',trim(var),'''') call halt_program endif icold=ic enddo enddo close(10) endif ! copy values to the pointvars array allocate(pointvars(nvars)) pointvars(:)=temparray(1:nvars) deallocate(temparray) deallocate(nvarpoint) allocate(par%pointtypes(size(pointtypes))) allocate(par%xpointsw(size(xpointsw))) allocate(par%ypointsw(size(ypointsw))) par%pointvars(1:nvars)=pointvars par%pointtypes=pointtypes par%xpointsw= xpointsw par%ypointsw=ypointsw endif ! xmaster ! end subroutine readpointvars subroutine readmeans(par) use logging_module use mnemmodule implicit none type(parameters),intent(inout) :: par integer :: id,ic,index,i character(128) :: line,keyword if (par%nmeanvar>0) then id=0 ! Look for keyword nmeanvar in params.txt if(xmaster) then open(10,file='params.txt') do while (id == 0) read(10,'(a)')line ic=scan(line,'=') if (ic>0) then keyword=adjustl(line(1:ic-1)) if (keyword == 'nmeanvar') id=1 endif enddo ! Read through the variables lines do i=1,par%nmeanvar read(10,'(a)')line index = chartoindex(trim(line)) if (index/=-1) then par%meansvars(i)=trim(line) call writelog('ls','(a)',' Will generate mean, min, max and variance output for variable: '//trim(mnemonics(index))) else call writelog('sle','',' Unknown point output type: ''',trim(line),'''') call halt_program endif end do close(10) endif endif end subroutine readmeans end module params