|
XBeach
|
00001 module params 00002 use typesandkinds, only: slen 00003 use mnemmodule, only: maxnamelen, numvars 00004 use xmpi_module 00005 use paramsconst 00006 implicit none 00007 save 00008 00009 private 00010 public parameters, params_inio, all_input 00011 #ifdef USEMPI 00012 public distribute_par 00013 #endif 00014 ! before using any routine from this module, take care of the value of params_inio 00015 ! .true. : values will be broadcasted to all processes, including output process 00016 ! otherwise: values will be broadcasted to compute processes only 00017 ! 00018 logical :: params_inio = .false. 00019 type parameters 00020 include 'paramsdecl.inc' 00021 end type parameters 00022 00023 contains 00024 00025 subroutine all_input(par) 00026 00027 use readkey_module, only: readkey_int, readkey_name, readkey_dbl, isSetParameter, parmapply, readkey_dblvec, readkey_inio 00028 use readkey_module, only: strippedline, setoldnames, setallowednames, readkey 00029 use xmpi_module 00030 use filefunctions, only: get_file_length, checkbcfilelength, check_file_exist, check_file_length 00031 use logging_module, only: writelog, report_file_read_error 00032 use constants, only: pi 00033 #ifdef USEMPI 00034 use xmpi_module 00035 #endif 00036 00037 implicit none 00038 00039 type(parameters), intent(inout) :: par 00040 00041 character(slen) :: testc,line 00042 character(slen) :: dummystring 00043 00044 integer :: filetype,mmax,nmax,ier,ic 00045 logical :: comment 00046 logical :: fe1,fe2,fe3 00047 00048 logical, parameter :: toall = .true. 00049 readkey_inio = toall 00050 ! 00051 call writelog('sl','','Reading input parameters: ') 00052 ! 00053 ! Check params.txt exists 00054 call check_file_exist('params.txt') 00055 ! 00056 ! 00057 ! Collections of default sets 00058 par%useXBeachGSettings = readkey_int ('params.txt','useXBeachGSettings', 0,0,1,strict=.true.,silent=.true.) 00059 ! par%useXBeachGSettings = 1 00060 ! 00061 ! 00062 ! Backward compatibility 00063 call check_instat_backward_compatibility(par) 00064 ! 00065 ! 00066 ! Physical processes 00067 call writelog('l','','--------------------------------') 00068 call writelog('l','','Physical processes: ') 00069 ! 00070 ! Wavemodel 00071 if (isSetParameter('params.txt','wavemodel')) then 00072 ! if defined in the params.txt 00073 call setallowednames('stationary', WAVEMODEL_STATIONARY, 'surfbeat', WAVEMODEL_SURFBEAT, 'nonh', WAVEMODEL_NONH) 00074 call setoldnames('0','1','2') 00075 if(par%useXBeachGSettings==0) then 00076 call parmapply('wavemodel',2, par%wavemodel, par%wavemodel_str) 00077 else 00078 call parmapply('wavemodel',3, par%wavemodel, par%wavemodel_str) 00079 endif 00080 elseif (par%wavemodel /= -123) then 00081 ! if already determined with backward compatibility (instat) 00082 if (par%wavemodel == 0) then 00083 call writelog('l','',' wavemodel =stationary') 00084 elseif (par%wavemodel == 1) then 00085 call writelog('l','',' wavemodel =surfbeat') 00086 else 00087 call writelog('l','',' wavemodel =nonh') 00088 endif 00089 elseif (par%wavemodel == -123) then 00090 ! no wavemodel determined -> error 00091 call writelog('lse','(a)','Error: XBeach cannot run without a hydrodynamic type') 00092 call writelog('lse','(a)',' Set wavemodel to stationary, surfbeat or nonh') 00093 call halt_program 00094 endif 00095 ! 00096 ! Rest of the processes 00097 par%cyclic = readkey_int ('params.txt','cyclic', 0, 0, 1,strict=.true.) 00098 if (par%wavemodel==WAVEMODEL_NONH) then 00099 par%swave = readkey_int ('params.txt','swave', 0, 0, 1,strict=.true.) 00100 else 00101 par%swave = readkey_int ('params.txt','swave', 1, 0, 1,strict=.true.) 00102 endif 00103 if(par%wavemodel==WAVEMODEL_SURFBEAT) then 00104 par%single_dir = readkey_int ('params.txt','single_dir', 0, 0, 1,strict=.true.) 00105 else 00106 par%single_dir = 0 00107 endif 00108 par%lwave = readkey_int ('params.txt','lwave', 1, 0, 1,strict=.true.) 00109 par%flow = readkey_int ('params.txt','flow', 1, 0, 1,strict=.true.) 00110 par%sedtrans = readkey_int ('params.txt','sedtrans', 1, 0, 1,strict=.true.) 00111 par%morphology = readkey_int ('params.txt','morphology', par%sedtrans, 0,1,strict=.true.) 00112 par%avalanching = readkey_int ('params.txt','avalanching', par%morphology,0,1,strict=.true.) 00113 par%gwflow = readkey_int ('params.txt','gwflow', par%useXBeachGSettings, 0, 1,strict=.true.) 00114 par%q3d = readkey_int ('params.txt','q3d', 0, 0, 1,silent=.true.,strict=.true.) 00115 par%swrunup = readkey_int ('params.txt','swrunup', 0, 0, 1,silent=.true.,strict=.true.) 00116 par%ships = readkey_int ('params.txt','ships', 0, 0, 1,strict=.true.) 00117 ! nship defined by shipfile 00118 if (par%ships == 0) par%nship = 0 00119 par%bchwiz = readkey_int ('params.txt','bchwiz', 0, 0, 1,silent=.true.,strict=.true.) 00120 par%vegetation = readkey_int ('params.txt','vegetation', 0, 0, 1,strict=.true.) 00121 par%setbathy = readkey_int ('params.txt','setbathy', 0, 0, 1,strict=.true.) 00122 par%viscosity = readkey_int ('params.txt','viscosity', 1, 0, 1,strict=.true.) 00123 par%advection = readkey_int ('params.txt','advection', 1, 0, 1,strict=.true.) 00124 par%wind = readkey_int ('params.txt','wind', 0, 0, 1,strict=.true.) 00125 par%rainfall = readkey_int ('params.txt','rainfall', 0, 0, 1,strict=.true.,silent=.true.) 00126 ! 00127 ! Grid parameters 00128 call writelog('l','','--------------------------------') 00129 call writelog('l','','Grid parameters: ') 00130 ! 00131 ! check gridform 00132 call setallowednames('xbeach',GRIDFORM_XBEACH,'delft3d',GRIDFORM_DELFT3D) 00133 call setoldnames('0','1') 00134 call parmapply('gridform',1,par%gridform,par%gridform_str) 00135 ! 00136 ! gridform switch 00137 if (par%gridform==GRIDFORM_XBEACH) then 00138 ! XBeach grid parameters 00139 par%xori = readkey_dbl('params.txt','xori', 0.d0, -1d9, 1d9) 00140 par%yori = readkey_dbl('params.txt','yori', 0.d0, -1d9, 1d9) 00141 par%alfa = readkey_dbl('params.txt','alfa', 0.d0, 0.d0, 360.d0) 00142 par%nx = readkey_int('params.txt','nx', 50, 2, 10000,required=.true.) 00143 par%ny = readkey_int('params.txt','ny', 2, 0, 10000,required=.true.) 00144 if(par%useXBeachGSettings==1) then 00145 if(par%ny>0) then 00146 call writelog('lws','(a)','XBeach-G settings cannot be applied in 2DH models') 00147 call writelog('lws','(a)','Set model to ny=0') 00148 call halt_program 00149 endif 00150 endif 00151 par%posdwn= readkey_dbl('params.txt','posdwn', 1.d0, -1.d0, 1.d0) 00152 if (par%setbathy .ne. 1 .and. par%hotstart .ne. 1) then 00153 par%depfile = readkey_name('params.txt','depfile',required=.true.) 00154 call check_file_exist(par%depfile) 00155 call check_file_length(par%depfile,par%nx+1,par%ny+1) 00156 else 00157 !par%depfile = readkey_name('params.txt','depfile') 00158 ! don't read this file 00159 endif 00160 par%vardx = readkey_int('params.txt','vardx', 0, 0, 1,strict=.true.) 00161 00162 if (par%vardx==0) then 00163 par%dx = readkey_dbl('params.txt','dx', -1.d0, 0.d0, 1d9,required=.true.) 00164 par%dy = readkey_dbl('params.txt','dy', -1.d0, 0.d0, 1d9,required=.true.) 00165 else 00166 par%dx = readkey_dbl('params.txt','dx', -1.d0, 0.d0, 1d9) ! bas: not required, but can be used in superfast 1D (smagorinsky and timestep) 00167 par%dy = readkey_dbl('params.txt','dy', -1.d0, 0.d0, 1d9) 00168 par%xfile = readkey_name('params.txt','xfile') 00169 call check_file_exist(par%xfile) 00170 call check_file_length(par%xfile,par%nx+1,par%ny+1) 00171 00172 par%yfile = readkey_name('params.txt','yfile') 00173 if (par%ny>0) then 00174 call check_file_exist(par%yfile) 00175 call check_file_length(par%yfile,par%nx+1,par%ny+1) 00176 end if 00177 00178 endif 00179 elseif (par%gridform==GRIDFORM_DELFT3D) then 00180 par%depfile = readkey_name('params.txt','depfile',required=.true.) 00181 call check_file_exist(par%depfile) 00182 par%dx = -1.d0 ! Why? 00183 par%dy = -1.d0 ! Why? 00184 par%xyfile = readkey_name('params.txt','xyfile',required=.true.) 00185 call check_file_exist(par%xyfile) 00186 ! read grid properties from xyfile 00187 if (xmaster) then 00188 open(31,file=par%xyfile,status='old',iostat=ier) 00189 ! skip comment text in file... 00190 comment=.true. 00191 do while (comment .eqv. .true.) 00192 read(31,'(a)',iostat=ier)line 00193 if (ier .ne. 0) then 00194 call report_file_read_error(par%xyfile) 00195 endif 00196 if (line(1:1)/='*') then 00197 comment=.false. 00198 endif 00199 enddo 00200 ! Check if grid coordinates are Cartesian 00201 ic=scan(line,'Cartesian') 00202 if (ic<=0) then 00203 call writelog('ewsl','','Delft3D grid is not Cartesian') 00204 call halt_program 00205 endif 00206 ! read grid dimensions 00207 read(31,*,iostat=ier) mmax,nmax 00208 ! catch new grid format that specifies missing value 00209 if (ier .ne. 0) then ! try reading the next line 00210 read(31,*,iostat=ier) mmax,nmax 00211 endif 00212 ! if still error, then XBeach cannot read this file 00213 if (ier .ne. 0) then 00214 call report_file_read_error(par%xyfile) 00215 endif 00216 close (31) 00217 endif 00218 #ifdef USEMPI 00219 call xmpi_bcast(mmax,toall) 00220 call xmpi_bcast(nmax,toall) 00221 #endif 00222 par%nx = mmax-1 00223 par%ny = nmax-1 00224 ! should we allow this input for gridform == 'Delft3D' 00225 par%xori = readkey_dbl('params.txt','xori', 0.d0, -1d9, 1d9) 00226 par%yori = readkey_dbl('params.txt','yori', 0.d0, -1d9, 1d9) 00227 par%alfa = readkey_dbl('params.txt','alfa', 0.d0, 0.d0, 360.d0) 00228 par%posdwn= readkey_dbl('params.txt','posdwn', 1.d0, -1.d0, 1.d0) 00229 endif 00230 ! Q3d grid 00231 par%nz = readkey_int ('params.txt','nz', 1, 1, 100) 00232 ! Wave directional grid 00233 if(par%swave==1) then 00234 par%thetamin = readkey_dbl ('params.txt','thetamin', -90.d0, -360.d0, 360.d0,required=.true.) 00235 par%thetamax = readkey_dbl ('params.txt','thetamax', 90.d0, -360.d0, 360.d0,required=.true.) 00236 par%thetanaut= readkey_int ('params.txt','thetanaut', 0, 0, 1,strict=.true.) 00237 if (par%single_dir==1) then 00238 call writelog('ls','','dtheta will automatically be computed from thetamin and thetamax for single_dir = 1') 00239 par%dtheta_s = readkey_dbl ('params.txt','dtheta_s', 10.d0, 0.1d0, 20.d0,required=.true.) 00240 else 00241 par%dtheta = readkey_dbl ('params.txt','dtheta', 10.d0, 0.1d0, 180.d0,required=.true.) 00242 endif 00243 endif 00244 ! 00245 ! Model time parameters 00246 call writelog('l','','--------------------------------') 00247 call writelog('l','','Model time parameters: ') 00248 if(par%useXBeachGSettings==0) then 00249 par%CFL = readkey_dbl ('params.txt','CFL', 0.7d0, 0.1d0, 0.9d0) 00250 else 00251 par%CFL = readkey_dbl ('params.txt','CFL', 0.6d0, 0.1d0, 0.9d0) 00252 endif 00253 par%dtset = readkey_dbl ('params.txt','dtset', 0.0d0, 0.001d0, 100.d0) 00254 par%tstop = readkey_dbl ('params.txt','tstop', 2000.d0, 1.d0, 1000000.d0,required=.true.) 00255 par%defuse = readkey_int ('params.txt','defuse', 1, 0, 1,strict=.true.,silent=.true.) 00256 if (par%wavemodel==WAVEMODEL_STATIONARY .or. par%wavemodel==WAVEMODEL_SURFBEAT) then 00257 par%maxdtfac = readkey_dbl ('params.txt','maxdtfac', 50.d0, 10.d0, 200.d0) 00258 else 00259 par%maxdtfac = readkey_dbl ('params.txt','maxdtfac', 500.d0, 100.d0, 1000.d0) 00260 endif 00261 ! 00262 ! Physical constants 00263 call writelog('l','','--------------------------------') 00264 call writelog('l','','Physical constants: ') 00265 par%rho = readkey_dbl ('params.txt','rho', 1025.0d0, 1000.0d0, 1040.0d0) 00266 par%g = readkey_dbl ('params.txt','g', 9.81d0, 9.7d0, 9.9d0) 00267 par%depthscale = readkey_dbl ('params.txt','depthscale',1.0d0, 1.0d0, 200.d0) 00268 ! 00269 ! Initial conditions 00270 call writelog('l','','--------------------------------') 00271 call writelog('l','','Initial conditions: ') 00272 par%zsinitfile = readkey_name('params.txt','zsinitfile') 00273 if (par%zsinitfile==' ') then 00274 ! do nothing 00275 else 00276 call check_file_exist(par%zsinitfile) 00277 if (par%gridform==GRIDFORM_XBEACH) then ! nx and ny not known in case of Delft3D 00278 call check_file_length(par%zsinitfile,par%nx+1,par%ny+1) 00279 endif 00280 endif 00281 par%hotstartflow = readkey_int ('params.txt','hotstartflow', 0, 0, 1,strict=.true.,silent=.true.) 00282 par%hotstart = readkey_int ('params.txt','hotstart', 0, 0, 1,strict=.true.,silent=.true.) 00283 if (par%hotstart==1) then 00284 par%hotstartfileno = readkey_int ('params.txt','hotstartfileno', 0, 0,999,required=.true.) 00285 endif 00286 ! 00287 ! Wave boundary condition parameters 00288 call writelog('l','','--------------------------------') 00289 call writelog('l','','Wave boundary condition parameters: ') 00290 ! 00291 ! New method: wbctype 00292 if (isSetParameter('params.txt','wbctype')) then 00293 call setallowednames('params', WBCTYPE_PARAMS, & 00294 'parametric', WBCTYPE_PARAMETRIC, & 00295 'swan', WBCTYPE_SWAN, & 00296 'vardens', WBCTYPE_VARDENS, & 00297 'off', WBCTYPE_OFF, & 00298 'jonstable', WBCTYPE_JONS_TABLE, & 00299 'reuse', WBCTYPE_REUSE, & 00300 'ts_1', WBCTYPE_TS_1, & 00301 'ts_2', WBCTYPE_TS_2, & 00302 'ts_nonh', WBCTYPE_TS_NONH) 00303 call setoldnames('0','1','2','3','4','5','6','7', '8', '9') 00304 call parmapply('wbctype', 1, par%wbctype, par%wbctype_str) 00305 elseif (par%wbctype /= -123) then 00306 ! if already determined with backward compatibility (instat) 00307 if (par%wbctype == WBCTYPE_PARAMS) then 00308 call writelog('l','',' wbctype =params') 00309 elseif (par%wbctype == WBCTYPE_PARAMETRIC) then 00310 call writelog('l','',' wbctype =parametric') 00311 elseif (par%wbctype == WBCTYPE_SWAN) then 00312 call writelog('l','',' wbctype =swan') 00313 elseif (par%wbctype == WBCTYPE_VARDENS) then 00314 call writelog('l','',' wbctype =vardens') 00315 elseif (par%wbctype == WBCTYPE_OFF) then 00316 call writelog('l','',' wbctype =off') 00317 elseif (par%wbctype == WBCTYPE_JONS_TABLE) then 00318 call writelog('l','',' wbctype =jonstable') 00319 elseif (par%wbctype == WBCTYPE_REUSE) then 00320 call writelog('l','',' wbctype =reuse') 00321 elseif (par%wbctype == WBCTYPE_TS_1) then 00322 call writelog('l','',' wbctype =ts_1') 00323 elseif (par%wbctype == WBCTYPE_TS_2) then 00324 call writelog('l','',' wbctype =ts_2') 00325 elseif (par%wbctype == WBCTYPE_TS_NONH) then 00326 call writelog('l','',' wbctype =ts_nonh') 00327 endif 00328 else 00329 ! no wbctype determined -> error 00330 call writelog('lse','(a)','Error: XBeach cannot run without a type of boundary condition') 00331 call writelog('lse','(a)',' Set wbctype') 00332 call halt_program 00333 endif 00334 ! 00335 if ( par%wbctype==WBCTYPE_PARAMETRIC .or. & 00336 par%wbctype==WBCTYPE_SWAN .or. & 00337 par%wbctype==WBCTYPE_VARDENS .or. & 00338 par%wbctype==WBCTYPE_JONS_TABLE & 00339 )then 00340 par%bcfile = readkey_name('params.txt','bcfile') 00341 call check_file_exist(par%bcfile) 00342 call checkbcfilelength(par%tstop,par%wbctype,par%bcfile,filetype) 00343 ! Only carried out on xmaster so: 00344 #ifdef USEMPI 00345 call xmpi_bcast(filetype,toall) 00346 #endif 00347 elseif (par%wbctype==WBCTYPE_REUSE .or. par%instat == INSTAT_REUSE) then 00348 ! See if this is reusing nonhydrostatic, or hydrostatic boundary conditions 00349 ! Note, check file length is done after recomputation of tstop due to morfacopt 00350 ! at the end of this subroutine. 00351 if (xmaster) then 00352 inquire(file='ebcflist.bcf',exist=fe1) 00353 inquire(file='qbcflist.bcf',exist=fe2) 00354 inquire(file='nhbcflist.bcf',exist=fe3) 00355 endif 00356 #ifdef USEMPI 00357 call xmpi_bcast(fe1,toall) 00358 call xmpi_bcast(fe2,toall) 00359 call xmpi_bcast(fe3,toall) 00360 #endif 00361 if (fe3 .and. .not. (fe1 .or. fe2)) then 00362 par%nonhspectrum = 1 00363 ! Check for file length is done later, after tstop is adjusted for morfac 00364 elseif (.not. fe3 .and. (fe1 .and. fe2)) then 00365 par%nonhspectrum = 0 00366 ! Check for file length is done later, after tstop is adjusted for morfac 00367 elseif (fe3 .and. (fe1 .or. fe2)) then 00368 call writelog('lswe','', & 00369 'If ''instat=reuse'' the model directory may not contain multiple boundary definition files.') 00370 call writelog('lswe','','Use either ebcflist.bcf/qbcflist.bcf, or nhbcflist.bcf') 00371 call halt_program 00372 elseif (.not. fe3 .and. .not. (fe1 .and. fe2)) then 00373 call writelog('lswe','', & 00374 'If ''instat=reuse'' the model directory may not contain sufficient boundary definition files.') 00375 if (.not. fe1) then 00376 call writelog('lswe','','Model currently missing ebcflist.bcf') 00377 elseif (.not. fe2) then 00378 call writelog('lswe','','Model currently missing qbcflist.bcf') 00379 endif 00380 call halt_program 00381 else 00382 call writelog('lswe','','If ''instat=reuse'' the model directory must contain boundary definition files.') 00383 call writelog('lswe','','Use either ebcflist.bcf/qbcflist.bcf, or nhbcflist.bcf') 00384 call halt_program 00385 endif 00386 else 00387 filetype=-1 00388 endif 00389 par%taper = readkey_dbl ('params.txt','taper', 100.d0, 0.0d0, 1000.d0) 00390 par%nmax = readkey_dbl ('params.txt','nmax', 0.8d0, 0.5d0, 1.d0) 00391 par%cyclicdiradjust = readkey_int ('params.txt','cyclicdiradjust', par%cyclic,0,1,strict=.true.,silent=.true.) 00392 if (par%wbctype == WBCTYPE_PARAMS.or. par%single_dir==1) then 00393 par%nonhspectrum = readkey_int ('params.txt','nonhspectrum', 0, 0, 1 ,strict=.true.) 00394 par%Hrms = readkey_dbl ('params.txt','Hrms', 1.d0, 0.d0, 10.d0) 00395 par%Tm01 = readkey_dbl ('params.txt','Tm01', 10.d0, 1.d0, 20.d0) 00396 par%Trep = readkey_dbl ('params.txt','Trep', par%Tm01, 1.d0, 20.d0) 00397 par%dir0 = readkey_dbl ('params.txt','dir0', 270.d0, -360.d0, 360.d0) 00398 par%m = readkey_int ('params.txt','m', 10, 2, 128) 00399 elseif (par%wbctype == WBCTYPE_TS_1 .or. par%wbctype == WBCTYPE_TS_2) then 00400 call check_file_exist('bc/gen.ezs') 00401 par%Tm01 = readkey_dbl ('params.txt','Tm01', 10.d0, 1.d0, 20.d0) 00402 par%Trep = readkey_dbl ('params.txt','Trep', par%Tm01, 1.d0, 20.d0) 00403 par%dir0 = readkey_dbl ('params.txt','dir0', 270.d0, -360.d0, 360.d0) 00404 par%m = readkey_int ('params.txt','m', 10, 2, 128) 00405 elseif (par%wbctype == WBCTYPE_TS_NONH) then 00406 call check_file_exist('boun_U.bcf') 00407 endif 00408 ! 00409 ! If Tlong is defined: than bichromatic waves 00410 if (isSetParameter('params.txt','Tlong')) then 00411 par%Tlong = readkey_dbl ('params.txt','Tlong', 80.d0, 20.d0, 300.d0) 00412 endif 00413 ! 00414 call setallowednames('neumann',LATERALWAVE_NEUMANN,'wavecrest',LATERALWAVE_WAVECREST,'cyclic',LATERALWAVE_CYCLIC) 00415 call setoldnames('0','1') 00416 call parmapply('lateralwave',1,par%lateralwave,par%lateralwave_str) 00417 par%bclwonly = readkey_int ('params.txt','bclwonly', 0,0,1,strict=.true.,silent=.true.) 00418 par%Sfold = readkey_int ('params.txt','Sfold', 0,0,1,strict=.true.,silent=.true.) 00419 ! 00420 ! 00421 ! Wave-spectrum boundary condition parameters 00422 if ( par%wbctype==WBCTYPE_PARAMETRIC .or. & 00423 par%wbctype==WBCTYPE_SWAN .or. & 00424 par%wbctype==WBCTYPE_VARDENS .or. & 00425 par%wbctype==WBCTYPE_JONS_TABLE & 00426 )then 00427 ! 00428 call writelog('l','','--------------------------------') 00429 call writelog('l','','Wave-spectrum boundary condition parameters: ') 00430 if (par%wavemodel==WAVEMODEL_NONH) then 00431 par%nonhspectrum= 1 00432 else 00433 par%nonhspectrum = 0 00434 endif 00435 par%random = readkey_int ('params.txt','random', 1, 0, 1 ,strict=.true.) 00436 par%fcutoff = readkey_dbl ('params.txt','fcutoff', 0.d0, 0.d0, 40.d0 ) 00437 par%nspr = readkey_int ('params.txt','nspr', 0, 0, 1 ,silent=.true.,strict=.true.) 00438 par%trepfac = readkey_dbl ('params.txt','trepfac', 0.01d0, 0.d0, 1.d0 ) 00439 if (par%nonhspectrum==1) then 00440 par%sprdthr = readkey_dbl ('params.txt','sprdthr', 0.00d0, 0.d0, 1.d0 ) 00441 else 00442 par%sprdthr = readkey_dbl ('params.txt','sprdthr', 0.08d0, 0.d0, 1.d0 ) 00443 endif 00444 par%correctHm0 = readkey_int ('params.txt','correctHm0', 1, 0, 1 ,silent=.true.,strict=.true.) 00445 par%Tm01switch = readkey_int ('params.txt','Tm01switch', 0, 0, 1 ,strict=.true.) 00446 ! 00447 if (filetype==0) then 00448 par%rt = readkey_dbl('params.txt','rt', min(3600.d0,par%tstop), 1200.d0, 7200.d0 ) 00449 par%dtbc = readkey_dbl('params.txt','dtbc', 1.0d0, 0.1d0, 2.0d0 ) 00450 endif 00451 ! 00452 if (par%wbctype==WBCTYPE_SWAN) then 00453 par%dthetaS_XB = readkey_dbl ('params.txt','dthetaS_XB', 0.0d0, -360.d0, 360.0d0,strict=.true. ) 00454 endif 00455 ! 00456 par%wbcversion = readkey_int ('params.txt','wbcversion', 3, 1, 3,strict=.true.,silent=.true.) ! wbcversion defaults to 3 00457 ! 00458 if (par%wbcversion>2) then 00459 par%nspectrumloc = readkey_int ('params.txt','nspectrumloc', 1, 1, par%ny+1 ) 00460 else 00461 par%nspectrumloc = 1 00462 endif 00463 par%wbcEvarreduce = readkey_dbl ('params.txt','wbcEvarreduce', 1.d0, 0.d0, 1.d0,strict=.true.,silent=.true. ) 00464 par%wbcQvarreduce = readkey_dbl ('params.txt','wbcQvarreduce', 1.d0, 0.d0, 1.d0,strict=.true.,silent=.true. ) 00465 endif 00466 ! 00467 ! Flow boundary condition parameters 00468 ! front 00469 call writelog('l','','--------------------------------') 00470 call writelog('l','','Flow boundary condition parameters: ') 00471 call setallowednames('abs_1d', FRONT_ABS_1D, & 00472 'abs_2d', FRONT_ABS_2D, & 00473 'wall', FRONT_WALL, & 00474 'wlevel', FRONT_WLEVEL, & 00475 'nonh_1d', FRONT_NONH_1D, & 00476 'waveflume', FRONT_WAVEFLUME) 00477 call setoldnames('0','1','2','3','4','5') 00478 if (par%nonhspectrum==1) then 00479 call parmapply('front',5,par%front,par%front_str) 00480 else 00481 if (par%ny>2) then 00482 call parmapply('front',2,par%front,par%front_str) 00483 else 00484 call parmapply('front',1,par%front,par%front_str) 00485 endif 00486 endif 00487 ! 00488 ! left and right 00489 call setallowednames('neumann', LR_NEUMANN, & 00490 'wall' , LR_WALL, & 00491 'no_advec', LR_NO_ADVEC, & 00492 'neumann_v', LR_NEUMANN_V, & 00493 'abs_1d', LR_ABS_1D) 00494 call setoldnames('0','1') 00495 call parmapply('left',1,par%left,par%left_str) 00496 00497 call parmapply('right',1,par%right,par%right_str) 00498 00499 ! back 00500 call setallowednames('wall', BACK_WALL, & 00501 'abs_1d', BACK_ABS_1D, & 00502 'abs_2d', BACK_ABS_2D, & 00503 'wlevel', BACK_WLEVEL) 00504 call setoldnames('0','1','2','3') 00505 if (par%ny>2) then 00506 call parmapply('back',3,par%back,par%back_str) 00507 else 00508 call parmapply('back',2,par%back,par%back_str) 00509 endif 00510 00511 ! Error for wlevel 00512 if ((par%front == FRONT_WLEVEL) .or. (par%back == BACK_WLEVEL)) then 00513 call writelog('lse','','Error: wlevel no longer supported. Change front and/or back boundary condition.') 00514 call halt_program 00515 endif 00516 00517 ! others 00518 par%ARC = readkey_int ('params.txt','ARC', 1, 0, 1 ,strict=.true.) 00519 par%order = readkey_dbl ('params.txt','order', 2.d0, 1.d0, 2.d0 ,strict=.true.) 00520 par%highcomp = readkey_int ('params.txt','highcomp', 0, 0, 1 ) 00521 par%freewave = readkey_int ('params.txt','freewave', 0, 0, 1 ,strict=.true.) 00522 par%epsi = readkey_dbl ('params.txt','epsi', -1.d0, -1.d0, 0.2d0 ) 00523 par%nc = readkey_int ('params.txt','nc', par%ny+1, 1, par%ny+1,strict=.true.,silent=.true.) 00524 00525 call setallowednames('instant', TIDETYPE_INSTANT, & 00526 'velocity', TIDETYPE_VELOCITY) 00527 call parmapply('tidetype',2,par%tidetype,par%tidetype_str) 00528 00529 ! 00530 ! 00531 ! Tide boundary conditions 00532 call writelog('l','','--------------------------------') 00533 call writelog('l','','Tide boundary conditions: ') 00534 par%tideloc = readkey_int ('params.txt','tideloc', 0, 0, 4) 00535 if (par%tideloc>0) then 00536 if (par%tideloc==2) then 00537 call setallowednames('land', PAULREVERE_LAND, & 00538 'sea', PAULREVERE_SEA) 00539 call setoldnames('0','1') 00540 call parmapply('paulrevere',1,par%paulrevere,par%paulrevere_str) 00541 endif 00542 par%zs0file = readkey_name('params.txt','zs0file') 00543 call check_file_exist(par%zs0file) 00544 else 00545 par%zs0 = readkey_dbl ('params.txt','zs0', 0.0d0, -5.d0, 5.d0) 00546 endif 00547 ! 00548 ! 00549 ! Discharge boundary conditions 00550 call writelog('l','','--------------------------------') 00551 call writelog('l','','Discharge boundary conditions: ') 00552 00553 par%disch_loc_file = readkey_name ('params.txt','disch_loc_file' ) 00554 par%disch_timeseries_file = readkey_name ('params.txt','disch_timeseries_file' ) 00555 00556 par%ndischarge = get_file_length(par%disch_loc_file ) 00557 par%ntdischarge = get_file_length(par%disch_timeseries_file ) 00558 00559 par%ndischarge = readkey_int ('params.txt', 'ndischarge', par%ndischarge, 0, 100) 00560 par%ntdischarge = readkey_int ('params.txt', 'ntdischarge', par%ntdischarge, 0, 100) 00561 00562 if (par%ndischarge>0) then 00563 call check_file_exist(par%disch_loc_file) 00564 00565 if (par%ntdischarge>0) then 00566 call check_file_exist(par%disch_timeseries_file) 00567 endif 00568 endif 00569 ! 00570 ! Wave breaking parameters 00571 par%beta = readkey_dbl ('params.txt','beta', 0.10d0, 0.05d0, 0.3d0) 00572 if (par%swave==1) then 00573 call writelog('l','','--------------------------------') 00574 call writelog('l','','Wave breaking parameters: ') 00575 call setallowednames('roelvink1', BREAK_ROELVINK1, & 00576 'baldock', BREAK_BALDOCK, & 00577 'roelvink2', BREAK_ROELVINK2, & 00578 'roelvink_daly', BREAK_ROELVINK_DALY, & 00579 'janssen', BREAK_JANSSEN) 00580 call setoldnames('1','2','3','4','5') 00581 if (par%wavemodel == WAVEMODEL_STATIONARY) then 00582 call parmapply('break',2,par%break,par%break_str) ! default: baldock 00583 par%gamma = readkey_dbl ('params.txt','gamma', 0.78d0, 0.4d0, 0.9d0) 00584 par%gammax = readkey_dbl ('params.txt','gammax', 0.6d0, .4d0, 5.d0) 00585 elseif (par%wavemodel == WAVEMODEL_SURFBEAT) then 00586 call parmapply('break',3,par%break,par%break_str) ! default: roelvink2 00587 par%gamma = readkey_dbl ('params.txt','gamma', 0.55d0, 0.4d0, 0.9d0) 00588 par%gammax = readkey_dbl ('params.txt','gammax', 2.d0, .4d0, 5.d0) 00589 else 00590 par%rollergammax = readkey_int ('params.txt','rollergammax', 1, 0, 1,silent=.true.,strict=.true.) 00591 endif 00592 if (par%break == BREAK_ROELVINK_DALY) then 00593 par%gamma2 = readkey_dbl ('params.txt','gamma2', 0.3d0, 0.0d0, 0.5d0) 00594 endif 00595 par%alpha = readkey_dbl ('params.txt','alpha', 1.0d0, 0.5d0, 2.0d0) 00596 par%n = readkey_dbl ('params.txt','n', 10.0d0, 5.0d0, 20.0d0) 00597 par%delta = readkey_dbl ('params.txt','delta', 0.0d0, 0.0d0, 1.0d0,strict=.true.) 00598 par%deltahmin = readkey_dbl ('params.txt','deltahmin', 0.0d0, 0.0d0, 1.0d0,strict=.true.) 00599 par%wavfriccoef = readkey_dbl ('params.txt','fw', 0.d0, 0d0, 1.0d0) 00600 ! try to read a wave friction file 00601 par%wavfricfile = readkey_name('params.txt','fwfile') 00602 if (par%wavfricfile .ne. ' ') then 00603 call check_file_exist(par%wavfricfile) 00604 if (par%gridform==GRIDFORM_XBEACH) then 00605 call check_file_length(par%wavfricfile,par%nx+1,par%ny+1) 00606 endif 00607 call writelog('lws','(a,a,a)','Warning: wave friction coefficient values from file ''',& 00608 trim(par%wavfricfile), & 00609 ''' will be used in computation') 00610 end if 00611 par%fwcutoff = readkey_dbl ('params.txt','fwcutoff', 1000.d0, 0d0, 1000.d0) 00612 par%breakerdelay = readkey_dbl ('params.txt','breakerdelay', 1.d0, 0.d0, 3.d0,strict=.true.) 00613 par%shoaldelay = readkey_int ('params.txt','shoaldelay', 0, 0, 1,silent=.true.,strict=.true.) 00614 if (par%shoaldelay==1) then 00615 par%facsd = readkey_dbl ('params.txt','facsd', 1.d0, 0d0, 2.0d0) 00616 endif 00617 if (par%swrunup==1) then 00618 par%facrun = readkey_dbl ('params.txt','facrun', 1.d0, 0d0, 2.0d0) 00619 endif 00620 ! 00621 ! Roller parameters 00622 call writelog('l','','--------------------------------') 00623 call writelog('l','','Roller parameters: ') 00624 par%roller = readkey_int ('params.txt','roller', 1, 0, 1,strict=.true.) 00625 par%rfb = readkey_int ('params.txt','rfb', 0, 0, 1,strict=.true.) 00626 ! 00627 ! Wave-current interaction parameters 00628 call writelog('l','','--------------------------------') 00629 call writelog('l','','Wave-current interaction parameters: ') 00630 par%wci = readkey_int ('params.txt','wci', 0, 0, 1,strict=.true.) 00631 par%hwci = readkey_dbl ('params.txt','hwci', 0.1d0, 0.001d0, 1.d0) 00632 par%hwcimax = readkey_dbl ('params.txt','hwcimax', 100.d0, 0.01d0, 100.d0) 00633 par%cats = readkey_dbl ('params.txt','cats', 4.d0, 1.d0, 50.d0) 00634 endif 00635 ! 00636 ! Flow parameters 00637 call writelog('l','','--------------------------------') 00638 call writelog('l','','Flow parameters: ') 00639 00640 call setallowednames('chezy', BEDFRICTION_CHEZY, & 00641 'cf', BEDFRICTION_CF, & 00642 'white-colebrook', BEDFRICTION_WHITE_COLEBROOK, & 00643 'manning', BEDFRICTION_MANNING, & 00644 'white-colebrook-grainsize', BEDFRICTION_WHITE_COLEBROOK_GRAINSIZE) 00645 if(par%useXBeachGSettings==0) then 00646 call parmapply('bedfriction',1,par%bedfriction,par%bedfriction_str) 00647 else 00648 call parmapply('bedfriction',5,par%bedfriction,par%bedfriction_str) 00649 endif 00650 00651 ! prevent using bed friction files without explicitly setting bedfriction type 00652 if (.not. isSetParameter('params.txt','bedfriction') .and. & 00653 isSetParameter('params.txt','bedfricfile')) then 00654 call writelog('lswe','','The use of keyword BEDFRICFILE without keyword BEDFRICTION is not allowed') 00655 call writelog('lswe','','Terminating simulation') 00656 call halt_program 00657 endif 00658 00659 if (par%bedfriction==BEDFRICTION_CHEZY .or. & 00660 par%bedfriction==BEDFRICTION_CF .or. & 00661 par%bedfriction==BEDFRICTION_MANNING .or. & 00662 par%bedfriction==BEDFRICTION_WHITE_COLEBROOK) then 00663 00664 ! try to read a bed friction file 00665 par%bedfricfile = readkey_name('params.txt','bedfricfile') 00666 if (par%bedfricfile .ne. ' ') then 00667 call check_file_exist(par%bedfricfile) 00668 if (par%gridform==GRIDFORM_XBEACH) then 00669 call check_file_length(par%bedfricfile,par%nx+1,par%ny+1) 00670 endif 00671 call writelog('lws','(a,a,a)','Warning: bed friction coefficient values from file ''',& 00672 trim(par%bedfricfile), & 00673 ''' will be used in computation') 00674 else 00675 if (par%bedfriction==BEDFRICTION_CHEZY) then 00676 par%bedfriccoef = readkey_dbl('params.txt','bedfriccoef', 55.d0, 20.d0, 100.d0) 00677 elseif (par%bedfriction==BEDFRICTION_CF) then 00678 par%bedfriccoef = readkey_dbl('params.txt','bedfriccoef', 3.d-3, 1.d-3, 0.1d0) 00679 elseif (par%bedfriction==BEDFRICTION_MANNING) then 00680 par%bedfriccoef = readkey_dbl('params.txt','bedfriccoef', 0.02d0, 0.01d0, 0.05d0) 00681 elseif (par%bedfriction==BEDFRICTION_WHITE_COLEBROOK) then 00682 par%bedfriccoef = readkey_dbl('params.txt','bedfriccoef', 0.01d0, 3.5d-5, 0.9d0) 00683 endif 00684 endif 00685 else 00686 ! other formulations do not require a bed friction coefficient 00687 par%bedfricfile = '' ! empty string so doesn't go searching for file called 'abc' 00688 par%bedfriccoef = -999.d0 00689 endif 00690 par%dynamrough = readkey_int ('params.txt','dynamicroughness ',0 , 0, 1,silent=.true.) 00691 par%droot = readkey_dbl('params.txt','droot', 0.5d0, 0.d0, 100d0) 00692 par%dstem = readkey_dbl('params.txt','dstem', 0.5d0, 0.d0, 100d0) 00693 par%maxcf = readkey_dbl ('params.txt','maxcf', 0.04d0, 0.0d0, 1.0d0) ! max cf, only used for Manning and White Colebrook (Chezy: 15) 00694 par%nuh = readkey_dbl ('params.txt','nuh', 0.1d0, 0.0d0, 1.0d0) 00695 par%nuhfac = readkey_dbl ('params.txt','nuhfac', 1.0d0, 0.0d0, 1.0d0) 00696 par%nuhv = readkey_dbl ('params.txt','nuhv', 1.d0, 1.d0, 20.d0,silent=.true.) 00697 par%smag = readkey_int ('params.txt','smag', 1, 0, 1,strict=.true.) 00698 00699 if (par%gwflow ==1) then 00700 if(par%useXBeachGSettings==1) then 00701 ! default on in XBeach-G! 00702 par%friction_infiltration = readkey_int('params.txt','friction_infiltration',1,0,1, strict=.true.) 00703 else 00704 par%friction_infiltration = readkey_int('params.txt','friction_infiltration',0,0,1, silent=.true.,strict=.true.) 00705 endif 00706 else 00707 par%friction_infiltration = 0 00708 endif 00709 par%friction_turbulence = readkey_int('params.txt','friction_turbulence' ,0,0,1, silent=.true.,strict=.true.) 00710 if (par%friction_turbulence==1) then 00711 par%gamma_turb = readkey_dbl ('params.txt','gamma_turb', 1.d0, 0.d0, 2.d0) 00712 endif 00713 00714 call setallowednames('none', CF_ACC_NONE, & 00715 'mccall', CF_ACC_MCCALL, & 00716 'nielsen',CF_ACC_NIELSEN) 00717 if(par%useXBeachGSettings==1) then 00718 ! Default should be MCCALL in XBeach-G 00719 call parmapply('friction_acceleration',2,par%friction_acceleration,par%friction_acceleration_str) 00720 else 00721 call parmapply('friction_acceleration',1,par%friction_acceleration,par%friction_acceleration_str,silent=.true.) 00722 endif 00723 00724 if(par%friction_acceleration==CF_ACC_MCCALL) then 00725 par%ci = readkey_dbl ('params.txt','ci',1.0d0, 0.5d0, 1.5d0) 00726 elseif (par%friction_acceleration==CF_ACC_NIELSEN) then 00727 par%phit = readkey_dbl ('params.txt','phit',25.0d0, 0.00d0, 90.0d0) 00728 endif 00729 00730 ! 00731 ! 00732 ! Coriolis force parameters 00733 call writelog('l','','--------------------------------') 00734 call writelog('l','','Coriolis force parameters: ') 00735 par%wearth = readkey_dbl ('params.txt','wearth', 1.d0/24.d0, 0.d0, 1.d0) 00736 par%lat = readkey_dbl ('params.txt','lat', 0.d0, -90.d0, 90.d0) 00737 ! 00738 ! 00739 ! Wind parameters 00740 if(par%wind==1) then 00741 call writelog('l','','--------------------------------') 00742 call writelog('l','','Wind parameters: ') 00743 par%rhoa = readkey_dbl ('params.txt','rhoa', 1.25d0, 1.0d0, 2.0d0) 00744 par%Cd = readkey_dbl ('params.txt','Cd', 0.002d0, 0.0001d0, 0.01d0) 00745 par%windfile = readkey_name('params.txt','windfile') 00746 if (par%windfile==' ') then 00747 par%windv = readkey_dbl ('params.txt','windv', 0.0d0, 0.0d0, 200.0d0) 00748 par%windth = readkey_dbl ('params.txt','windth', 270.0d0, -360.0d0, 360.0d0) 00749 else 00750 call check_file_exist(par%windfile) 00751 endif 00752 else 00753 par%windv = 0.d0 00754 par%windth = 0.d0 00755 par%windfile = ' ' 00756 endif 00757 ! 00758 ! 00759 ! Rainfall parameters 00760 if(par%rainfall==1) then 00761 call writelog('l','','--------------------------------') 00762 call writelog('l','','Rainfall parameters: ') 00763 par%rainfallratefile = readkey_name('params.txt','rainfallratefile') 00764 if (par%rainfallratefile==' ') then 00765 par%rainfallrate = readkey_dbl ('params.txt','rainfallrate', 0.0d0, 0.0d0, 2000.0d0) 00766 par%nrainfallrate = 1 00767 else 00768 call check_file_exist(par%rainfallratefile) 00769 par%nrainfallrate = readkey_int ('params.txt','nrainfallrate',1,1,1000) 00770 endif 00771 endif 00772 ! 00773 ! 00774 ! Groundwater parameters 00775 if (par%gwflow==1) then 00776 call writelog('l','','--------------------------------') 00777 call writelog('l','','Groundwater parameters: ') 00778 if(par%useXBeachGSettings==1) then 00779 par%kx = readkey_dbl ('params.txt','kx' , 0.01d0 , 0.001d0, 0.9d0) 00780 par%ky = readkey_dbl ('params.txt','ky' , par%kx , 0.001d0, 0.9d0) 00781 par%kz = readkey_dbl ('params.txt','kz' , par%kx , 0.001d0, 0.9d0) 00782 else 00783 par%kx = readkey_dbl ('params.txt','kx' , 0.0001d0 , 0.00001d0, 0.1d0) 00784 par%ky = readkey_dbl ('params.txt','ky' , par%kx , 0.00001d0, 0.1d0) 00785 par%kz = readkey_dbl ('params.txt','kz' , par%kx , 0.00001d0, 0.1d0) 00786 endif 00787 par%dwetlayer = readkey_dbl ('params.txt','dwetlayer' , 0.1d0 , 0.01d0 , 1.d0) 00788 par%aquiferbotfile = readkey_name('params.txt','aquiferbotfile') 00789 if (par%aquiferbotfile==' ') then 00790 !also read in groundwater.f90 which determines value 00791 par%aquiferbot = readkey_dbl('params.txt','aquiferbot',-10.d0,-100.d0,100.d0) 00792 else 00793 call check_file_exist(par%aquiferbotfile) 00794 endif 00795 par%gw0file = readkey_name('params.txt','gw0file') 00796 if (par%gw0file==' ') then 00797 par%gw0 = readkey_dbl('params.txt','gw0',0.d0,-5.d0,5.d0) 00798 else 00799 call check_file_exist(par%gw0file) 00800 endif 00801 par%gwnonh = readkey_int ('params.txt','gwnonh',par%useXBeachGSettings, 0, 1,strict=.true.) 00802 if (par%gwnonh==1) then 00803 if (par%ny>2) then 00804 par%gwfastsolve = readkey_int ('params.txt','gwfastsolve', 0, 0, 1,silent=.true.,strict=.true.) 00805 endif 00806 endif 00807 00808 ! Type of momentum equation 00809 call setallowednames('laminar', GWSCHEME_LAMINAR, & 00810 'turbulent', GWSCHEME_TURBULENT) 00811 call setoldnames('darcy','modflow') 00812 if(par%useXBeachGSettings==1) then 00813 call parmapply('gwscheme',2,par%gwscheme,par%gwscheme_str) 00814 else 00815 call parmapply('gwscheme',1,par%gwscheme,par%gwscheme_str) 00816 endif 00817 00818 if (par%gwscheme==GWSCHEME_TURBULENT) then 00819 par%gwReturb = readkey_dbl ('params.txt','gwReturb' , 100.d0 , 1.d0 , 600.d0) 00820 endif 00821 call setallowednames('parabolic', GWHEADMODEL_PARABOLIC, & 00822 'exponential', GWHEADMODEL_EXPONENTIAL) 00823 call parmapply('gwheadmodel',1,par%gwheadmodel,par%gwheadmodel_str) 00824 00825 par%gwhorinfil = readkey_int ('params.txt','gwhorinfil', 0, 0, 1,strict=.true.,silent=.true.) 00826 endif 00827 ! 00828 ! Non-hydrostatic correction parameters 00829 if (par%wavemodel==WAVEMODEL_NONH) then 00830 call writelog('l','','--------------------------------') 00831 call writelog('l','','Non-hydrostatic correction parameters: ') 00832 call setallowednames('sip', SOLVER_SIPP, & 00833 'tridiag', SOLVER_TRIDIAGG) 00834 call setoldnames('1','2') 00835 if (par%ny>2) then 00836 call parmapply('solver',1,par%solver,par%solver_str) ! default: sip 00837 else 00838 call parmapply('solver',2,par%solver,par%solver_str) ! default: tridiag 00839 endif 00840 if (par%solver==SOLVER_SIPP) then 00841 par%solver_maxit = readkey_int('params.txt','solver_maxit' ,30,1,1000) 00842 par%solver_acc = readkey_dbl('params.txt','solver_acc' ,0.005d0,0.00001d0,0.1d0) 00843 par%solver_urelax= readkey_dbl('params.txt','solver_urelax' ,0.92d0,0.5d0,0.99d0) 00844 endif 00845 par%kdmin = readkey_dbl('params.txt','kdmin' ,0.0d0,0.0d0,0.05d0) 00846 par%Topt = readkey_dbl('params.txt','Topt', 10.d0, 1.d0, 20.d0) 00847 par%nonhq3d = readkey_int('params.txt','nonhq3d' ,0,0,1,silent=.true.,strict=.true.) 00848 if (par%nonhq3d==1) then 00849 par%nhbreaker = readkey_int('params.txt','nhbreaker' ,1,0,1,strict=.true.) 00850 par%nhlay = readkey_dbl('params.txt','nhlay' ,0.33d0,0.d0,1.d0) 00851 else 00852 if(par%useXBeachGSettings==0) then 00853 par%nhbreaker = readkey_int('params.txt','nhbreaker' ,1,0,1,strict=.true.) 00854 else 00855 par%nhbreaker = readkey_int('params.txt','nhbreaker' ,1,0,3,strict=.true.) 00856 endif 00857 par%dispc = readkey_dbl('params.txt','dispc' ,-1.0d0,0.1d0,2.0d0) 00858 endif 00859 00860 if (par%nhbreaker==1) then 00861 par%breakviscfac = readkey_dbl('params.txt','breakviscfac',1.5d0, 1.d0, 3.d0) 00862 par%maxbrsteep = readkey_dbl('params.txt','maxbrsteep',0.4d0, 0.3d0, 0.8d0) 00863 par%reformsteep = readkey_dbl('params.txt','reformsteep',0.25d0*par%maxbrsteep,0.d0,0.95d0*par%maxbrsteep) 00864 elseif (par%nhbreaker==2) then 00865 par%breakvisclen = readkey_dbl('params.txt','breakvisclen',1.0d0, 0.75d0, 3.d0) 00866 par%maxbrsteep = readkey_dbl('params.txt','maxbrsteep',0.4d0, 0.3d0, 0.8d0) 00867 par%secbrsteep = readkey_dbl('params.txt','secbrsteep',0.5d0*par%maxbrsteep,0.d0,0.95d0*par%maxbrsteep) 00868 elseif (par%nhbreaker==3) then 00869 !par%avis = readkey_dbl ('params.txt','avis', 0.3d0, 0.0d0, 1.0d0) 00870 par%breakvisclen = readkey_dbl('params.txt','breakvisclen',1.0d0, 0.75d0, 3.d0) 00871 par%maxbrsteep = readkey_dbl('params.txt','maxbrsteep',0.6d0, 0.3d0, 0.8d0) 00872 par%secbrsteep = readkey_dbl('params.txt','secbrsteep',0.5d0*par%maxbrsteep,0.d0,0.95d0*par%maxbrsteep) 00873 endif 00874 endif 00875 ! 00876 ! 00877 ! Sediment transport parameters 00878 if (par%sedtrans==1) then 00879 call writelog('l','','--------------------------------') 00880 call writelog('l','','Sediment transport parameters: ') 00881 00882 call setallowednames('soulsby_vanrijn', FORM_SOULSBY_VANRIJN, & 00883 'vanthiel_vanrijn', FORM_VANTHIEL_VANRIJN, & 00884 'vanrijn1993', FORM_VANRIJN1993, & 00885 'nielsen2006', FORM_NIELSEN2006, & 00886 'mccall_vanrijn', FORM_MCCALL_VANRIJN, & 00887 'wilcock_crow', FORM_WILCOCK_CROW, & 00888 'engelund_fredsoe', FORM_ENGELUND_FREDSOE, & 00889 'mpm', FORM_MPM, & 00890 'wong_parker', FORM_WONG_PARKER, & 00891 'fl_vb', FORM_FL_VB, & 00892 'fredsoe_deigaard', FORM_FREDSOE_DEIGAARD) 00893 call setoldnames('1','2','3','4','5','6','7','8','9','10','11') 00894 if(par%useXBeachGSettings==0) then 00895 call parmapply('form',2,par%form, par%form_str) 00896 else 00897 call parmapply('form',5,par%form, par%form_str) 00898 endif 00899 00900 if (par%wavemodel==WAVEMODEL_STATIONARY .or. par%wavemodel==WAVEMODEL_SURFBEAT) then 00901 call setallowednames('ruessink_vanrijn', WAVEFORM_RUESSINK_VANRIJN, & 00902 'vanthiel', WAVEFORM_VANTHIEL) 00903 call setoldnames('1','2') 00904 call parmapply('waveform',2,par%waveform,par%waveform_str) 00905 par%sws = readkey_int ('params.txt','sws', 1, 0, 1,strict=.true.) 00906 par%lws = readkey_int ('params.txt','lws', 1, 0, 1,strict=.true.) 00907 par%BRfac = readkey_dbl ('params.txt','BRfac', 1.0d0, 0.d0, 1.d0) 00908 par%facua = readkey_dbl ('params.txt','facua ',0.10d0, 0.00d0, 1.0d0) 00909 par%facSk = readkey_dbl ('params.txt','facSk ',par%facua, 0.00d0, 1.0d0) 00910 par%facAs = readkey_dbl ('params.txt','facAs ',par%facua, 0.00d0, 1.0d0) 00911 if (par%waveform == WAVEFORM_VANTHIEL) then 00912 par%Tbfac = readkey_dbl ('params.txt','Tbfac ',1.0d0, 0.00d0, 1.0d0) 00913 endif 00914 call setallowednames('none', TURB_NONE, & 00915 'wave_averaged', TURB_WAVE_AVERAGED, & 00916 'bore_averaged', TURB_BORE_AVERAGED) 00917 call setoldnames('0','1','2') 00918 call parmapply('turb',3,par%turb,par%turb_str) 00919 if (par%waveform == WAVEFORM_RUESSINK_VANRIJN) then 00920 call parmapply('turb',2,par%turb,par%turb_str) 00921 endif 00922 call setallowednames('none', TURBADV_NONE, & 00923 'lagrangian', TURBADV_LAGRANGIAN, & 00924 'eulerian', TURBADV_EULERIAN) 00925 call parmapply('turbadv',1,par%turbadv,par%turbadv_str) 00926 00927 else 00928 par%sws = 0 00929 par%lws = 1 00930 par%BRfac = 0.d0 00931 par%facua = 0.d0 00932 par%facSk = 0.d0 00933 par%facAs = 0.d0 00934 par%turb = TURB_WAVE_AVERAGED ! note, needed in case we want to use lwt ... Better solution still waiting 00935 par%turbadv = TURBADV_LAGRANGIAN 00936 endif 00937 00938 ! Parameters for gravel-type equations 00939 if(par%form==FORM_NIELSEN2006 .or. & 00940 par%form==FORM_MCCALL_VANRIJN .or. & 00941 par%form==FORM_WILCOCK_CROW .or. & 00942 par%form==FORM_ENGELUND_FREDSOE .or. & 00943 par%form==FORM_MPM .or. & 00944 par%form==FORM_WONG_PARKER .or. & 00945 par%form==FORM_FL_VB .or. & 00946 par%form==FORM_FREDSOE_DEIGAARD) then 00947 00948 if(par%form==FORM_NIELSEN2006) then 00949 par%Ctrans = readkey_dbl ('params.txt','Ctrans',12d0, 0d0, 30d0) 00950 par%thetcr = readkey_dbl ('params.txt','thetcr ',0.05d0, 0.00d0, 1.0d0) 00951 par%phaselag = readkey_int ('params.txt','phaselag', 0, 0, 1) 00952 if (par%phaselag==1) then 00953 par%phit = readkey_dbl ('params.txt','phit ',25.00d0, 0.00d0, 90.0d0) 00954 endif 00955 call setallowednames('constant', SEDFRICFAC_CONSTANT, & 00956 'flowfricfac', SEDFRICFAC_FLOWFRIC, & 00957 'nielsen', SEDFRICFAC_NIELSEN, & 00958 'swart', SEDFRICFAC_SWART, & 00959 'wilson', SEDFRICFAC_WILSON) 00960 call setoldnames('0','1','2','3','4') 00961 call parmapply('sedfricfac',1,par%sedfricfac,par%sedfricfac_str) 00962 if (par%sedfricfac==SEDFRICFAC_SWART) then 00963 par%Arms = readkey_dbl ('params.txt','Arms', 10.d0, 0.d0, 100.d0) 00964 endif 00965 if (par%sedfricfac==SEDFRICFAC_CONSTANT) then 00966 par%fsed = readkey_dbl ('params.txt','fsed', 0.025d0, 0.005d0, 0.05d0) 00967 endif 00968 par%streaming = readkey_int ('params.txt','streaming', 0, 0, 1,silent=.true.,strict=.true.) 00969 else ! all other transport equations 00970 par%thetcr = readkey_dbl ('params.txt','thetcr ',-10.d0, 0.00d0, 1.0d0,silent=.true.) ! negative = self-compute 00971 par%sedfricfac = SEDFRICFAC_FLOWFRIC ! not actually used in the equations, all comes from taubx 00972 par%uprushfac = readkey_dbl ('params.txt','uprushfac', 1.d0, 0.d0, 3.d0) 00973 par%backwashfac = readkey_dbl ('params.txt','backwashfac', 1.d0, 0.d0, 3.d0) 00974 par%incldzdx = readkey_int('params.txt','incldzdx',0,0,1,strict=.true.,silent=.true.) 00975 if (par%gwflow==1) then 00976 par%inclrelweight = readkey_int('params.txt','inclrelweight',1,0,1,strict=.true.) 00977 else 00978 par%inclrelweight = 0 00979 endif 00980 endif 00981 00982 call setallowednames('none', SLOPECORR_NONE, & 00983 'nielsen', SLOPECORR_NIELSEN, & 00984 'fredsoe_deigaard', SLOPECORR_FREDSOE_DEIGAARD) 00985 call setoldnames('0','1','2') 00986 call parmapply('slopecorr',3,par%slopecorr,par%sedfricfac_str) 00987 00988 par%sus = readkey_int ('params.txt','sus ',0, 0, 0,strict=.true.) 00989 par%bed = readkey_int ('params.txt','bed ',1, 1, 1,strict=.true.) 00990 par%bulk = readkey_int ('params.txt','bulk ',0, 0, 0,strict=.true.) 00991 if(par%bulk==0) then 00992 par%facsl = 0.d0 00993 else 00994 par%facsl = readkey_dbl ('params.txt','facsl ',0.0d0, 0.d0, 1.6d0) 00995 endif 00996 par%bdslpeffmag = BDSLPEFFMAG_NONE 00997 par%bdslpeffini = BDSLPEFFINI_NONE 00998 par%bdslpeffdir = BDSLPEFFDIR_NONE 00999 par%smax = -1.d0 01000 else ! non-gravel transport equations 01001 par%sus = readkey_int ('params.txt','sus ',1, 0, 1,strict=.true.) 01002 par%bed = readkey_int ('params.txt','bed ',1, 0, 1,strict=.true.) 01003 par%bulk = readkey_int ('params.txt','bulk ',0, 0, 1,strict=.true.) 01004 par%facsl = readkey_dbl ('params.txt','facsl ', 0.15d0, 0.d0, 1.6d0) 01005 par%z0 = readkey_dbl ('params.txt','z0 ',0.006d0, 0.0001d0, 0.05d0) 01006 par%smax = readkey_dbl ('params.txt','smax', -1.d0, -1.d0, 3.d0) !changed 28/11 and back 10/2 01007 call setallowednames('none', BDSLPEFFMAG_NONE, & 01008 'roelvink_total', BDSLPEFFMAG_ROELV_TOTAL, & 01009 'roelvink_bed', BDSLPEFFMAG_ROELV_BED, & 01010 'soulsby_total', BDSLPEFFMAG_SOULS_TOTAL, & 01011 'soulsby_bed', BDSLPEFFMAG_SOULS_BED) 01012 call setoldnames('0','1','2','3','4') 01013 call parmapply('bdslpeffmag',2,par%bdslpeffmag) 01014 01015 call setallowednames('none', BDSLPEFFINI_NONE, & 01016 'total', BDSLPEFFINI_TOTAL, & 01017 'bed', BDSLPEFFINI_BED) 01018 call setoldnames('0','1','2') 01019 call parmapply('bdslpeffini',1,par%bdslpeffini) 01020 01021 call setallowednames('none', BDSLPEFFDIR_NONE, & 01022 'talmon', BDSLPEFFDIR_TALMON) 01023 call setoldnames('0','1') 01024 call parmapply('bdslpeffdir',1,par%bdslpeffdir) 01025 if (par%bdslpeffdir>0) then 01026 par%bdslpeffdirfac = readkey_dbl ('params.txt','bdslpeffdirfac', 1.d0, 0.d0, 2.d0) 01027 endif 01028 endif 01029 if(par%useXBeachGSettings==0) then 01030 par%reposeangle = readkey_dbl ('params.txt','reposeangle', 30.d0, 0.d0, 45.d0) 01031 else 01032 par%reposeangle = readkey_dbl ('params.txt','reposeangle', 35.d0, 20.d0, 60.d0) 01033 endif 01034 par%tsfac = readkey_dbl ('params.txt','tsfac', 0.1d0, 0.01d0, 1.d0) 01035 par%Tsmin = readkey_dbl ('params.txt','Tsmin ',0.5d0, 0.01d0, 10.d0) 01036 par%facDc = readkey_dbl ('params.txt','facDc ',1.0d0, 0.00d0, 1.0d0) 01037 par%jetfac = readkey_dbl ('params.txt','jetfac ',0.0d0, 0.00d0, 1.0d0,silent=.true.) 01038 par%lwt = readkey_int ('params.txt','lwt ',0, 0, 1,strict=.true.) 01039 if(par%lwt==1 .and. par%form==FORM_NIELSEN2006) then 01040 par%yturb = readkey_dbl ('params.txt','yturb', 1.d0, 0.d0, 1.d0) 01041 par%facthr = readkey_dbl ('params.txt','facthr', 1.0d0, 1.d0, 100.d0) 01042 endif 01043 par%betad = readkey_dbl ('params.txt','betad ',1.0d0, 0.00d0, 10.0d0) 01044 par%fallvelred = readkey_int ('params.txt','fallvelred', 0, 0, 1,strict=.true.) 01045 par%dilatancy = readkey_int ('params.txt','dilatancy', 0, 0, 1,strict=.true.) 01046 if (par%dilatancy==1) then 01047 par%rheeA = readkey_dbl ('params.txt','rheeA', 0.75d0, 0.75d0, 2.d0) ! Between 3/4 and 1/(1-n), see paper Van Rhee (2010) 01048 par%pormax = readkey_dbl ('params.txt','pormax', 0.5d0, 0.3d0, 0.6d0) 01049 endif 01050 par%bermslopetransport = readkey_int ('params.txt','bermslopetransport', 0, 0, 1,strict=.true.,silent=.true.) 01051 if (par%bermslopetransport==1) then 01052 par%bermslopebed = readkey_int ('params.txt','bermslopebed', 1, 0, 1,strict=.true.) 01053 par%bermslopesus = readkey_int ('params.txt','bermslopesus', 1, 0, 1,strict=.true.) 01054 par%bermslope = readkey_dbl ('params.txt','bermslope ',0.1d0, 0.00d0, 1.0d0) 01055 par%bermslopefac = readkey_dbl ('params.txt','bermslopefac ',15.0d0, 0.00d0, 30.0d0) 01056 if (par%wavemodel==WAVEMODEL_SURFBEAT) then 01057 par%bermslopegamma = readkey_dbl ('params.txt','bermslopegamma ',1.0d0, 0.75d0, 1.5d0) 01058 par%bermslopedepth = readkey_dbl ('params.txt','bermslopedepth ',0.0d0, 0.0d0, 0.5d0) 01059 else 01060 par%bermslopedepth = readkey_dbl ('params.txt','bermslopedepth ',1.0d0, 0.5d0, 1.5d0) 01061 endif 01062 else 01063 par%bermslopebed = 0 01064 par%bermslopesus = 0 01065 endif 01066 par%alfaD50 = readkey_dbl ('params.txt','alfaD50', 0.0d0, 0.0d0, 1.5d0, silent=.true.) 01067 endif 01068 ! 01069 ! Q3D sediment transport parameters 01070 ! -> also make this part active if van Rijn, 1993 is used as ceq solver 01071 if (par%q3d==1 .or. par%form==FORM_VANRIJN1993) then 01072 call writelog('l','','--------------------------------') 01073 call writelog('l','','Q3D sediment transport parameters: ') 01074 par%vonkar = readkey_dbl ('params.txt','vonkar', 0.4d0, 0.01d0, 1.d0) 01075 par%vicmol = readkey_dbl ('params.txt','vicmol', 0.000001d0, 0.d0, 0.001d0) 01076 par%kmax = readkey_int ('params.txt','kmax ', 100, 1, 1000) 01077 par%nz = par%kmax ! work-around; discuss with Dano later difference nz / kmax 01078 par%sigfac = readkey_dbl ('params.txt','sigfac ',1.3d0, 0.00d0, 10.d0) 01079 par%deltar = readkey_dbl ('params.txt','deltar ',0.025d0, 0.001d0, 1.0d0) 01080 par%rwave = readkey_dbl ('params.txt','rwave ',2.0d0, 0.1d0, 10.0d0) 01081 else ! workaround, needed for array allocation 01082 par%kmax = 1 01083 endif 01084 ! 01085 ! Bed composition parameters 01086 call writelog('l','','--------------------------------') 01087 call writelog('l','','Bed composition parameters: ') 01088 par%ngd = readkey_int ('params.txt','ngd', 1, 1, 20,strict=.true.) 01089 par%nd = readkey_int ('params.txt','nd ', 3, 3, 1000,strict = .true.) 01090 par%por = readkey_dbl ('params.txt','por', 0.4d0, 0.3d0, 0.5d0) 01091 if (par%dilatancy==1) then 01092 par%D15 = readkey_dblvec('params.txt','D15',par%ngd,size(par%D15),0.00015d0,0.00001d0,0.0008d0) ! Lodewijk 01093 endif 01094 if(par%useXBeachGSettings==1) then 01095 par%D50 = readkey_dblvec('params.txt','D50',par%ngd,size(par%D50),0.010d0,0.002d0,0.08d0) 01096 par%D90 = readkey_dblvec('params.txt','D90',par%ngd,size(par%D90),0.015d0,0.003d0,0.12d0) 01097 else 01098 par%D50 = readkey_dblvec('params.txt','D50',par%ngd,size(par%D50),0.0002d0,0.00005d0,0.0008d0) 01099 par%D90 = readkey_dblvec('params.txt','D90',par%ngd,size(par%D90),0.0003d0,0.00010d0,0.0015d0) 01100 endif 01101 01102 if (par%sedtrans==1) then 01103 par%rhos = readkey_dbl ('params.txt','rhos', 2650d0, 2400.d0, 2800.d0) 01104 par%dzg1 = readkey_dbl ('params.txt','dzg', 0.1d0, 0.01d0, 1.d0) 01105 par%dzg1 = readkey_dbl ('params.txt','dzg1', par%dzg1, 0.01d0, 1.d0) 01106 par%dzg2 = readkey_dbl ('params.txt','dzg2', par%dzg1, 0.01d0, 1.d0) 01107 par%dzg3 = readkey_dbl ('params.txt','dzg3', par%dzg1, 0.01d0, 1.d0) 01108 ! file,keyword,size read vector,size vector in par,default,min,max 01109 par%sedcal = readkey_dblvec('params.txt','sedcal',par%ngd,size(par%sedcal),1.d0,0.d0,2.d0) 01110 par%ucrcal = readkey_dblvec('params.txt','ucrcal',par%ngd,size(par%ucrcal),1.d0,0.d0,2.d0) 01111 endif 01112 ! 01113 ! Morphology parameters 01114 if (par%morphology==1) then 01115 call writelog('l','','--------------------------------') 01116 call writelog('l','','Morphology parameters: ') 01117 par%morfac = readkey_dbl ('params.txt','morfac', 1.0d0, 0.d0, 1000.d0) 01118 par%morfacopt= readkey_int ('params.txt','morfacopt', 1, 0, 1,strict=.true.) 01119 par%morstart = readkey_dbl ('params.txt','morstart',0.0d0, 0.d0, par%tstop) 01120 par%morstop = readkey_dbl ('params.txt','morstop', par%tstop, 0.d0, 10000000.d0) 01121 if(par%useXBeachGSettings==1) then 01122 par%wetslp = readkey_dbl ('params.txt','wetslp', tan(par%reposeangle/180.d0*4.d0*atan(1.d0)),0.1d0,1.d0) 01123 par%dryslp = readkey_dbl ('params.txt','dryslp', par%wetslp, 0.1d0, 2.d0) 01124 else 01125 par%wetslp = readkey_dbl ('params.txt','wetslp', 0.3d0, 0.1d0, 1.d0) 01126 par%dryslp = readkey_dbl ('params.txt','dryslp', 1.0d0, 0.1d0, 2.d0) 01127 endif 01128 if (par%ny == 0) then 01129 par%lsgrad = readkey_dbl ('params.txt','lsgrad', 0.0d0, -0.1d0, .1d0,silent=.true.) 01130 else 01131 par%lsgrad = 0.0d0 01132 endif 01133 par%hswitch = readkey_dbl ('params.txt','hswitch',0.1d0, 0.01d0, 1.0d0) 01134 par%dzmax = readkey_dbl ('params.txt','dzmax ',0.05d0, 0.00d0, 1.0d0) 01135 par%struct = readkey_int ('params.txt','struct ',0 , 0, 1,strict=.true.) 01136 if (par%struct==1 .and. par%hotstart .ne. 1) then 01137 par%ne_layer = readkey_name('params.txt','ne_layer') 01138 call check_file_exist(par%ne_layer) 01139 if (par%gridform==GRIDFORM_XBEACH) then 01140 call check_file_length(par%ne_layer,par%nx+1,par%ny+1) 01141 endif 01142 endif 01143 if (par%swrunup==1 .and. par%struct==0) then 01144 call writelog('lws','(a)', & 01145 'Warning: swrunup can only be used in combination with struct=1. swrunup will be turned off.') 01146 par%swrunup = 0; 01147 endif 01148 endif 01149 ! 01150 ! 01151 ! Output variables 01152 call writelog('l','','--------------------------------') 01153 call writelog('l','','Output variables: ') 01154 par%timings = readkey_int ('params.txt','timings', 1, 0, 1,strict=.true.) 01155 testc = readkey_name('params.txt','tunits') 01156 if (len(trim(testc)) .gt. 0) par%tunits = trim(testc) 01157 par%tstart = readkey_dbl ('params.txt','tstart', 0.d0, 0.d0,par%tstop,strict=.true.) 01158 par%tint = readkey_dbl ('params.txt','tint', 1.d0, .0d0, par%tstop-par%tstart,strict=.true.) ! Robert 01159 par%tsglobal = readkey_name('params.txt','tsglobal') 01160 if (par%tsglobal==' ') then 01161 par%tintg = readkey_dbl ('params.txt','tintg', par%tint, 0.d0, par%tstop-par%tstart,strict=.true.) ! Robert 01162 endif 01163 par%tspoints = readkey_name('params.txt','tspoints') 01164 if (par%tspoints==' ') then 01165 par%tintp = readkey_dbl ('params.txt','tintp', par%tint, 0.d0, par%tstop-par%tstart,strict=.true.) ! Robert 01166 endif 01167 par%tsmean = readkey_name('params.txt','tsmean') 01168 if (par%tsmean==' ') then 01169 par%tintm = readkey_dbl ('params.txt','tintm', par%tstop-par%tstart, 0.d0, par%tstop-par%tstart,strict=.true.) ! Robert 01170 endif 01171 01172 ! hotstart output 01173 par%writehotstart = readkey_int('params.txt','writehotstart',0,0,1,strict=.true.,silent=.true.) 01174 if (par%writehotstart==1) then 01175 par%tshotstart = readkey_name('params.txt','tshotstart',silent=.true.) 01176 if (par%tshotstart==' ') then 01177 par%tinth = readkey_dbl ('params.txt','tinth', par%tstop, 0.d0, par%tstop,strict=.true.) ! Robert 01178 endif 01179 endif 01180 01181 ! global output 01182 par%nglobalvar = readkey_int ('params.txt','nglobalvar', -1, -1, 20) 01183 call readglobalvars(par) 01184 01185 ! point output 01186 par%npoints = readkey_int ('params.txt','npoints', 0, 0, 50) 01187 par%nrugauge = readkey_int ('params.txt','nrugauge', 0, 0, 50) 01188 par%npointvar = readkey_int ('params.txt','npointvar', 0, 0, 50) 01189 call readpointvars(par) 01190 par%nrugdepth = readkey_int('params.txt','nrugdepth',1,1,10) 01191 par%rugdepth = readkey_dblvec('params.txt','rugdepth',par%nrugdepth,size(par%rugdepth),0.0d0,0.0d0,0.1d0) 01192 01193 ! mean output 01194 par%nmeanvar = readkey_int ('params.txt','nmeanvar' , 0, 0, 15) 01195 call readmeans(par) 01196 01197 call setallowednames('fortran', OUTPUTFORMAT_FORTRAN, & 01198 'netcdf ', OUTPUTFORMAT_NETCDF, & 01199 'debug ', OUTPUTFORMAT_DEBUG) 01200 #ifdef USENETCDF 01201 ! Default to NetCDF output in case of NetCDF-enabled executable 01202 call parmapply('outputformat',2,par%outputformat,par%outputformat_str, & 01203 required = .false.) ! wwvv-todo 01204 #else 01205 call parmapply('outputformat',1,par%outputformat,par%outputformat_str, & 01206 required = .false.) ! wwvv-todo 01207 #endif 01208 if(par%outputformat==OUTPUTFORMAT_NETCDF .or. & 01209 par%outputformat==OUTPUTFORMAT_DEBUG) then 01210 ! type of output precision 01211 call setallowednames('single', OUTPUTPRECISION_SINGLE, & 01212 'double', OUTPUTPRECISION_DOUBLE) 01213 call parmapply('outputprecision',2,par%outputprecision,par%outputprecision_str,required = .false.) 01214 ! get the nc output file name from the parameter file 01215 par%ncfilename = readkey_name('params.txt','ncfilename') 01216 if (len(trim(par%ncfilename)) .eq. 0) par%ncfilename = 'xboutput.nc' 01217 call writelog('ls','','netcdf output to:' // par%ncfilename) 01218 endif 01219 if(par%outputformat==OUTPUTFORMAT_NETCDF .and. par%useXBeachGSettings==0) then 01220 par%remdryoutput = readkey_int ('params.txt','remdryoutput', 1, 0, 1,strict=.true.) 01221 else 01222 par%remdryoutput = readkey_int ('params.txt','remdryoutput', 0, 0, 1,strict=.true.) 01223 endif 01224 ! 01225 ! 01226 ! Output projection 01227 call writelog('l','','--------------------------------') 01228 call writelog('l','','Output projection: ') 01229 testc = readkey_name('params.txt','projection') 01230 if (len(trim(testc)) .gt. 0) par%projection = testc 01231 par%rotate = readkey_int ('params.txt','rotate', 1, 0, 1,strict=.true.) 01232 ! 01233 ! 01234 ! Drifters parameters 01235 if (isSetParameter('params.txt','drifterfile')) then 01236 call writelog('l','','--------------------------------') 01237 call writelog('l','','Drifters parameters: ') 01238 par%drifterfile = readkey_name ('params.txt', 'drifterfile' ) 01239 call check_file_exist(par%drifterfile) 01240 par%ndrifter = get_file_length(par%drifterfile ) 01241 par%ndrifter = readkey_int ('params.txt', 'ndrifter', par%ndrifter, 0, 50 ) 01242 else 01243 par%ndrifter = 0 01244 endif 01245 ! 01246 ! 01247 ! Shipwaves parameters 01248 if (par%ships==1) then 01249 call writelog('l','','--------------------------------') 01250 call writelog('l','','Shipwaves parameters: ') 01251 par%shipfile = readkey_name ('params.txt', 'shipfile') 01252 call check_file_exist(par%shipfile) 01253 ! shipfile routine should set nship 01254 endif 01255 ! 01256 ! 01257 ! Vegetation parameters 01258 if (par%vegetation==1) then 01259 call writelog('l','','--------------------------------') 01260 call writelog('l','','Vegetation parameters: ') 01261 par%veggiefile = readkey_name ('params.txt', 'veggiefile' ) 01262 call check_file_exist(par%veggiefile) 01263 par%veggiemapfile = readkey_name ('params.txt', 'veggiemapfile' ) 01264 call check_file_exist(par%veggiemapfile) 01265 par%Trep = readkey_dbl ('params.txt','Trep', 1.d0, 0.01d0, 20.d0) 01266 par%vegnonlin = readkey_int ('params.txt', 'vegnonlin',0,0,1,silent=.true.) 01267 par%vegcanflo = readkey_int ('params.txt', 'vegcanflo',0,0,1,silent=.true.) 01268 par%veguntow = readkey_int ('params.txt', 'veguntow', 1,0,1,silent=.true.) 01269 par%porcanflow = readkey_int ('params.txt', 'porcanflow', 0,0,1) 01270 if (par%porcanflow ==1) then 01271 par%Kp = readkey_dbl ('params.txt','Kp', 0.0001d0, 0.d0, 1.d0) 01272 par%Cm = readkey_dbl ('params.txt','Cm', 1.0d0, 0.d0, 2.d0) 01273 endif 01274 ! veggiefile routine should set nveg 01275 endif 01276 ! 01277 ! 01278 ! Wave numerics parameters 01279 call writelog('l','','--------------------------------') 01280 call writelog('l','','Wave numerics parameters: ') 01281 01282 call setallowednames('upwind_1', SCHEME_UPWIND_1, & 01283 'lax_wendroff', SCHEME_LAX_WENDROFF, & 01284 'upwind_2', SCHEME_UPWIND_2, & 01285 'warmbeam', SCHEME_WARMBEAM) 01286 call setoldnames('1','2','3','4') 01287 call parmapply('scheme',4,par%scheme,par%scheme_str) 01288 01289 if (par%wavemodel == WAVEMODEL_STATIONARY .or. par%single_dir==1) then 01290 par%wavint = readkey_dbl ('params.txt','wavint', 60.d0, 1.d0, 3600.d0) 01291 par%maxerror = readkey_dbl ('params.txt','maxerror', 0.0001d0, 0.00001d0, 0.001d0) 01292 par%maxerror_angle = readkey_dbl ('params.txt','maxerror_angle', 1.d0, 0.1d0, 5.d0) 01293 par%maxerror_angle = par%maxerror_angle/180*par%px 01294 !if (par%single_dir==1) then 01295 ! par%maxerror = readkey_dbl ('params.txt','maxerror', 0.005d0, 0.00001d0, 0.001d0) 01296 !else 01297 ! par%maxerror = readkey_dbl ('params.txt','maxerror', 0.0005d0, 0.00001d0, 0.001d0) 01298 !endif 01299 par%maxiter = readkey_int ('params.txt','maxiter', 500, 2, 1000) 01300 endif 01301 ! only default to Snell's Law if in 1D and only one directional bin 01302 if (par%single_dir == 0) then 01303 if (par%ny==0 .and. nint(abs(par%thetamax-par%thetamin)/par%dtheta)<2) then 01304 par%snells = readkey_int ('params.txt','snells', 1, 0, 1,strict=.true.) 01305 else 01306 par%snells = readkey_int ('params.txt','snells', 0, 0, 1,strict=.true.) 01307 endif 01308 else 01309 par%snells = readkey_int ('params.txt','snells', 0, 0, 1,silent=.true.,strict=.true.) 01310 endif 01311 if(par%nonhspectrum == 0) then 01312 par%swkhmin = readkey_dbl ('params.txt','swkhmin', -0.01d0, -0.01d0, 0.35d0,silent=.true.) ! Robert for Daan 01313 endif 01314 par%oldhmin = readkey_int('params.txt','oldhmin' ,0,0,1,silent=.true.,strict=.true.) 01315 ! 01316 ! 01317 ! Flow numerics parameters 01318 call writelog('l','','--------------------------------') 01319 call writelog('l','','Flow numerics parameters: ') 01320 par%eps = readkey_dbl ('params.txt','eps', 0.005d0, 0.001d0, 0.1d0) 01321 par%eps_sd = readkey_dbl ('params.txt','eps_sd', 0.5d0, 0.000d0, 1.0d0) 01322 par%umin = readkey_dbl ('params.txt','umin', 0.0d0, 0.0d0, 0.2d0) 01323 par%hmin = readkey_dbl ('params.txt','hmin', 0.2d0, 0.001d0, 1.d0) 01324 if (par%wavemodel==WAVEMODEL_NONH) then 01325 par%secorder = readkey_int('params.txt','secorder' ,1,0,1,strict=.true.) 01326 else 01327 par%secorder = readkey_int('params.txt','secorder' ,0,0,1,strict=.true.) 01328 endif 01329 par%oldhu = readkey_int('params.txt','oldhu' ,0,0,1,silent=.true.,strict=.true.) 01330 ! 01331 ! 01332 ! Sediment transport numerics parameters 01333 if (par%sedtrans==1) then 01334 call writelog('l','','--------------------------------') 01335 call writelog('l','','Sediment transport numerics parameters: ') 01336 par%thetanum = readkey_dbl ('params.txt','thetanum', 1.d0, 0.5d0, 1.d0) 01337 par%sourcesink = readkey_int ('params.txt','sourcesink ',0, 0, 1,silent=.true.,strict=.true.) 01338 par%cmax = readkey_dbl ('params.txt','cmax', 0.1d0, 0.0d0, 1.d0) 01339 par%oldTsmin = readkey_int ('params.txt','oldTsmin ',0, 0, 1,silent=.true.,strict=.true.) 01340 if (par%oldTsmin==0) then 01341 par%dtlimTs = readkey_dbl ('params.txt','dtlimTs', 5.0d0, 0.0d0, 20.d0) 01342 endif 01343 endif 01344 ! 01345 ! 01346 ! Bed update numerics parameters 01347 if (par%morphology==1) then 01348 call writelog('l','','--------------------------------') 01349 call writelog('l','','Bed update numerics parameters: ') 01350 par%frac_dz = readkey_dbl ('params.txt','frac_dz', 0.7d0, 0.5d0, 0.98d0) 01351 par%nd_var = readkey_int ('params.txt','nd_var', 2, 1, par%nd,strict=.true.) 01352 par%split = readkey_dbl ('params.txt','split', 1.01d0, 1.005d0, 1.10d0) 01353 par%merge = readkey_dbl ('params.txt','merge', 0.01d0, 0.005d0, 0.10d0) 01354 par%fixedavaltime = readkey_int ('params.txt','fixedavaltime', 0, 0,1,strict=.true.,silent=.true.) 01355 if (par%fixedavaltime==1) then 01356 par%avaltime = readkey_dbl ('params.txt','avaltime', 10.0d0, 0.001d0, 1000.0d0) 01357 else 01358 par%nTrepavaltime = readkey_dbl ('params.txt','nTrepavaltime', 1.0d0, 0.001d0, 10.0d0) 01359 endif 01360 endif 01361 ! 01362 ! Prescribed bathy update 01363 if (par%setbathy==1) then 01364 par%nsetbathy = readkey_int ('params.txt','nsetbathy',1,1,1000) 01365 par%setbathyfile = readkey_name ('params.txt', 'setbathyfile',required=.true. ) 01366 call check_file_exist(par%setbathyfile) 01367 endif 01368 ! 01369 ! 01370 ! MPI parameters 01371 #ifdef USEMPI 01372 call writelog('l','','--------------------------------') 01373 call writelog('l','','MPI parameters: ') 01374 01375 call setallowednames('auto', MPIBOUNDARY_AUTO, & 01376 'x', MPIBOUNDARY_X, & 01377 'y', MPIBOUNDARY_Y, & 01378 'man', MPIBOUNDARY_MAN) 01379 call parmapply('mpiboundary',1,par%mpiboundary,par%mpiboundary_str) 01380 if (par%mpiboundary == MPIBOUNDARY_MAN) then 01381 par%mmpi= readkey_int('params.txt','mmpi',2,1,100) 01382 par%nmpi= readkey_int('params.txt','nmpi',4,1,100) 01383 endif 01384 #endif 01385 ! 01386 ! 01387 ! 01388 ! 01389 ! Finish 01390 call writelog('l','','--------------------------------') 01391 call writelog('sl','','Finished reading input parameters') 01392 call writelog('l','','--------------------------------') 01393 ! 01394 ! 01395 ! ------------------- Post-input processing ------------------------- 01396 ! 01397 ! 01398 ! Transport equations (gravel) that don't work in 2DH 01399 if(par%form==FORM_NIELSEN2006 .or. & 01400 par%form==FORM_MCCALL_VANRIJN .or. & 01401 par%form==FORM_WILCOCK_CROW .or. & 01402 par%form==FORM_ENGELUND_FREDSOE .or. & 01403 par%form==FORM_MPM .or. & 01404 par%form==FORM_WONG_PARKER .or. & 01405 par%form==FORM_FL_VB .or. & 01406 par%form==FORM_FREDSOE_DEIGAARD) then 01407 01408 ! These transport equations only work in 1D 01409 if (par%ny>0) then 01410 call writelog('lswe','(a)','Error: The following sediment transport equations cannot be used in 2DH simulations:') 01411 call writelog('lswe','(a)','nielsen2006') 01412 call writelog('lswe','(a)','mccall_vanrijn') 01413 call writelog('lswe','(a)','wilcock_crow') 01414 call writelog('lswe','(a)','engelund_fredsoe') 01415 call writelog('lswe','(a)','mpm') 01416 call writelog('lswe','(a)','wong_parker') 01417 call writelog('lswe','(a)','fl_vb') 01418 call writelog('lswe','(a)','fredsoe_deigaard') 01419 call writelog('lswe','(a)','Change to other transport equation or 1D model') 01420 call halt_program 01421 endif 01422 endif 01423 ! 01424 ! 01425 ! Fix input parameters for choosen depthscale 01426 if (par%depthscale .ne. 1.d0) then 01427 par%eps = par%eps/par%depthscale 01428 par%hmin = par%hmin/par%depthscale 01429 par%hswitch = par%hswitch/par%depthscale 01430 par%dzmax = par%dzmax/par%depthscale**1.5d0 01431 par%maxerror= par%maxerror/par%depthscale 01432 01433 call writelog('lws','(a)','Warning: input parameters eps, hmin, hswitch and dzmax are scaled with') 01434 call writelog('lws','(a)',' depthscale to:') 01435 call writelog('lws','(a,f0.4)','eps = ', par%eps) 01436 call writelog('lws','(a,f0.4)','hmin = ', par%hmin) 01437 call writelog('lws','(a,f0.4)','hswitch = ',par%hswitch) 01438 call writelog('lws','(a,f0.4)','dzmax = ', par%dzmax) 01439 call writelog('lws','(a,f0.4)','maxerror = ', par%maxerror) 01440 01441 endif 01442 ! 01443 ! 01444 ! Constants 01445 par%rhog8 = 1.0d0/8.0d0*par%rho*par%g 01446 ! 01447 ! 01448 if (par%posdwn<0.1d0) then 01449 par%posdwn=-1.d0 ! Backward compatibility, now posdwn = 0 also works in input (i.e. posdwn = false) 01450 endif 01451 ! 01452 ! 01453 ! Stop useless physical processes 01454 if (par%sedtrans==0 .and. par%morphology==1) then 01455 call writelog('lse','(a)','Error: Morphology cannot be computed without sediment transport.') 01456 call writelog('lse','(a)',' Set sedtrans=1 or morphology=0') 01457 call halt_program 01458 endif 01459 if (par%morphology==0 .and. par%avalanching==1) then 01460 call writelog('lsw','(a)','Warning: Avalanching cannot be computed without morphology.') 01461 call writelog('lsw','(a)',' Avalanching has been turned off') 01462 par%avalanching=0 01463 endif 01464 if (par%setbathy == 1 .and. par%morphology==1) then 01465 call writelog('lsw','(a)','Morphology and avalanching has been turned off for prescibed bed update (setbathy=1)') 01466 par%morphology=0 01467 par%avalanching=0 01468 endif 01469 ! 01470 ! Wavemodel and wbctypes: warnings 01471 ! swave and nonh 01472 if (par%swave == 1 .and. par%wavemodel == WAVEMODEL_NONH) then 01473 call writelog('lsw','(a)','Warning: swave with wavemodel = nonh not possible') 01474 par%swave = 0 01475 endif 01476 ! 01477 ! Wave model and wbctypes: errors (halt program) 01478 ! Bichromatic and not surfbeat 01479 if (par%wbctype == WBCTYPE_PARAMS .and. par%Tlong /= -123 .and. par%wavemodel /= WAVEMODEL_SURFBEAT) then 01480 call writelog('lwse','(a)','Error: bichromatic waves can only be used in surfbeat mode') 01481 call halt_program 01482 endif 01483 ! ts_nonh and not nonh 01484 if (par%wbctype == WBCTYPE_TS_NONH .and. par%wavemodel /= WAVEMODEL_NONH) then 01485 call writelog('lwse','(a)','Error: nonh time series can only be used in nonh mode') 01486 call halt_program 01487 endif 01488 ! ts_1 or ts_2 01489 if (par%wbctype == WBCTYPE_TS_1 .or. par%wbctype == WBCTYPE_TS_2) then 01490 if (par%wavemodel /= WAVEMODEL_SURFBEAT) then 01491 call writelog('lwse','(a)','Error: ts_1 and ts_2 time series can only be used in surfbeat mode') 01492 call halt_program 01493 endif 01494 endif 01495 ! errors for nonh = 1 and a wavemodel which is not wavemodel nonh 01496 if (par%nonh == 1 .and. par%wavemodel /= WAVEMODEL_NONH) then 01497 call writelog('lwse','(a)','Error: found both wavemodel and nonh=1. Define hydrodynamic mode with wavemodel') 01498 call halt_program 01499 endif 01500 ! errors for wavemodel = stationary and wbctype = reuse 01501 if (par%wbctype == WBCTYPE_REUSE .and. par%wavemodel == WAVEMODEL_STATIONARY) then 01502 call writelog('lwse','(a)','Error: wbctype = reuse cannot be used in stationary mode') 01503 call halt_program 01504 endif 01505 ! errors for wavemodel = stationary and wbctype = swan 01506 if (par%wbctype == WBCTYPE_SWAN .and. par%wavemodel == WAVEMODEL_STATIONARY) then 01507 call writelog('lwse','(a)','Error: wbctype = swan cannot be used in stationary mode') 01508 call halt_program 01509 endif 01510 ! errors for wavemodel = vardens and wbctype = swan 01511 if (par%wbctype == WBCTYPE_VARDENS .and. par%wavemodel == WAVEMODEL_STATIONARY) then 01512 call writelog('lwse','(a)','Error: wbctype = vardens cannot be used in stationary mode') 01513 call halt_program 01514 endif 01515 ! errors for wavemodel = vardens and wbctype = swan 01516 if (par%wbctype == WBCTYPE_PARAMETRIC .and. par%wavemodel == WAVEMODEL_STATIONARY) then 01517 call writelog('lwse','(a)','Error: wbctype = parametric cannot be used in stationary mode') 01518 call halt_program 01519 endif 01520 ! 01521 ! Cyclic boundary conditions 01522 #ifdef USEMPI 01523 if(par%cyclic==1) then 01524 call writelog('lsw','(a)','Warning: Cyclic boundary conditions only work in combination with MPI. Re-running this') 01525 call writelog('lsw','(a)',' simulation without MPI will not be possible') 01526 endif 01527 #else 01528 if(par%cyclic==1) then 01529 call writelog('lswe','(a)','Error: Cyclic boundary conditions only work in combination with MPI.') 01530 call writelog('lsw','(a)',' Choose an MPI-compatable version of XBeach for this option') 01531 call halt_program 01532 endif 01533 #endif 01534 ! 01535 ! 01536 ! Set taper to non-zero 01537 par%taper = max(par%taper,1.d-6) 01538 ! 01539 ! 01540 ! Compute Coriolis 01541 par%lat = par%lat*par%px/180.d0 01542 par%wearth = par%px*par%wearth/1800.d0 01543 ! 01544 ! 01545 ! Only allow Baldock or Janssen in stationary mode and Roelvink in non-stationary 01546 if (par%swave==1) then 01547 if (par%wavemodel == WAVEMODEL_STATIONARY) then 01548 if (par%break .ne. BREAK_BALDOCK .and. par%break .ne. BREAK_JANSSEN) then 01549 if(par%break == BREAK_ROELVINK_DALY) then 01550 call writelog('lwse','','Error: Roelvink-Daly formulations not implemented in stationary wave mode,') 01551 call writelog('lwse','',' use Baldock or Janssen formulation.') 01552 call halt_program 01553 else 01554 call writelog('lwse','','Error: Roelvink formulations not allowed in stationary,') 01555 call writelog('lwse','',' use Baldock or Janssen formulation.') 01556 call halt_program 01557 endif 01558 endif 01559 else 01560 if (par%break == BREAK_BALDOCK) then 01561 call writelog('lws','','Warning: Baldock formulation not allowed in non-stationary, use break=roelvink1') 01562 call writelog('lws','',' or roelvink2 or roelvink_daly formulation.') 01563 endif 01564 if (par%break == BREAK_JANSSEN) then 01565 call writelog('lws','','Warning: Janssen formulation not allowed in non-stationary, use break=roelvink1 ') 01566 call writelog('lws','',' or roelvink2 or roelvink_daly formulation.') 01567 endif 01568 endif 01569 endif 01570 ! 01571 ! 01572 ! Only allow bore-averaged turbulence in combination with vanthiel waveform 01573 if ((par%waveform .ne. WAVEFORM_VANTHIEL) .and. (par%turb .eq. TURB_BORE_AVERAGED)) then 01574 call writelog('lse','','Error: Cannot compute bore-averaged turbulence without vanthiel wave form.') 01575 call writelog('lse','',' Please set waveform=vanthiel in params.txt, or choose another') 01576 call writelog('lse','',' turbulence model') 01577 call halt_program 01578 endif 01579 ! 01580 ! 01581 ! Set smax to huge if default is specified, but allow for computations to be 01582 ! done with smax wwvv 01583 if (par%smax<0) par%smax=huge(0.d0)*1.0d-20 01584 ! 01585 ! 01586 ! Source-sink check 01587 if (par%morfac>1.d0) then 01588 if (par%sourcesink==1) then 01589 call writelog('lws','','Warning: Using source-sink terms for bed level change with morfac can lead to') 01590 call writelog('lws','',' loss of sediment mass conservation.') 01591 endif 01592 endif 01593 ! 01594 ! 01595 ! Check maximum grain size Soulsby-Van Rijn 01596 if ((par%form==FORM_SOULSBY_VANRIJN) .and. (maxval(par%D50(1:par%ngd))>= 0.05d0)) then 01597 call writelog('lews','','Error: Soulsby - Van Rijn cannot be used for sediment D50 > 0.05m') 01598 call halt_program 01599 endif 01600 ! 01601 ! 01602 ! Check on pormax value 01603 if (par%dilatancy==1) then 01604 if (par%pormax.le.par%por) then 01605 call writelog('lws','','Warning: Maximum porosity for dilatancy effect smaller than actual porosity.') 01606 call writelog('lws','',' Setting pormax equal to por, effectively switching off dilatancy.') 01607 par%pormax = par%por 01608 end if 01609 endif 01610 ! 01611 ! 01612 ! Bermslope only in 1D 01613 if(par%bermslopetransport==1 .and. par%ny>2) then 01614 call writelog('lws','(a)','Bermslope transport option should not be applied in 2DH models') 01615 call writelog('lws','(a)','Bermslope correction only applied in M-direction') 01616 !call halt_program 01617 endif 01618 01619 ! 01620 ! 01621 ! If using tide, epsi should be on 01622 if (par%tideloc>0) then 01623 if (par%epsi<-1.d0) then 01624 call writelog('lws','','Automatically computing epsi using offshore boundary conditions') 01625 ! par%epsi = 0.05d0 --> Jaap do this in boundary conditions to account for vary wave conditions during simulation 01626 endif 01627 endif 01628 ! 01629 ! is not using nonh with front = nonh_1d 01630 if (par%front==FRONT_NONH_1D) then 01631 if (par%wavemodel /= WAVEMODEL_NONH) then 01632 call writelog('lws','','Warning: not recommended to use front = nonh_1d without wavemodel = nonh') 01633 endif 01634 endif 01635 ! 01636 ! is using nonh without front = nonh_1d 01637 if (par%front/=FRONT_NONH_1D) then 01638 if (par%wavemodel == WAVEMODEL_NONH) then 01639 call writelog('lws','','Warning: automatically changing front to nonh_1d with wavemodel = nonh') 01640 par%front = FRONT_NONH_1D 01641 endif 01642 endif 01643 ! is using wavemodel = stat without front = abs_1d 01644 if (par%front/=FRONT_ABS_1D) then 01645 if (par%wavemodel == WAVEMODEL_STATIONARY ) then 01646 call writelog('lws','','Warning: automatically changing front to abs_1d with wavemodel = stationary') 01647 par%front = FRONT_ABS_1D 01648 endif 01649 endif 01650 ! 01651 ! If using nonh, secorder should always be on 01652 if (par%wavemodel==WAVEMODEL_NONH) then 01653 if (par%secorder==0) then 01654 call writelog('lws','','Warning: Automatically turning on 2nd order correction in flow for') 01655 call writelog('lws','',' non-hydrostatic module [secorder=1]') 01656 par%secorder = 1 01657 endif 01658 endif 01659 01660 ! 01661 ! If using nonh, then the solver type is set by the grid size 01662 if (par%wavemodel==WAVEMODEL_NONH) then 01663 if (par%ny>2 .and. par%solver==SOLVER_TRIDIAGG) then 01664 call writelog('lswe','','Tri-diagonal solver cannot be used if ny>2') 01665 call halt_program 01666 endif 01667 if (par%ny==0 .and. par%solver==SOLVER_SIPP) then 01668 call writelog('lse','','SIP solver cannot be used if ny==0') 01669 call halt_program 01670 endif 01671 endif 01672 ! 01673 ! 01674 ! If generating spectral time series for nonhydrostatic waves, you need at least wbcversion 3 01675 if (par%nonhspectrum==1 .and. par%wbcversion<3) then 01676 call writelog('lws','','Warning: Automatically changing to wbcversion 3 for') 01677 call writelog('lws','',' non-hydrostatic spectral boundary condition [nonhspectrum=1]') 01678 par%wbcversion=3 01679 endif 01680 if (par%wbcversion==1 .or. par%wbcversion==2) then 01681 call writelog('lwse','','************************** ERROR ******************************') 01682 call writelog('lwse','','wbcversion 1 and 2 are no longer supported, from v1.21 onwards') 01683 call writelog('lwse','','The current default wbcversion is 3') 01684 call writelog('lwse','','***************************************************************') 01685 call halt_program 01686 endif 01687 ! 01688 ! 01689 ! If using nhbreaker then maxbrsteep should be larger than reformsteep and reformsteep should be 01690 ! greater than zero 01691 if (par%nhbreaker==1) then 01692 if (par%reformsteep>0.95d0*par%maxbrsteep) then 01693 par%reformsteep=0.95d0*par%maxbrsteep 01694 call writelog('lws','(a)','Warning: Setting reformsteep to maximum of 95% of maxbrsteep') 01695 elseif (par%reformsteep<0.0d0) then 01696 par%reformsteep=0.0d0 01697 call writelog('lws','(a)','Warning: Setting reformsteep to minimum of zero') 01698 endif 01699 endif 01700 ! 01701 ! 01702 ! The number of layers in the bed should be at least 3 (par%nd>=3) 01703 if (par%nd<3) then 01704 call writelog('lwse','','The number of sediment layers in the bed (nd) may not be smaller than 3') 01705 call halt_program 01706 endif 01707 ! 01708 ! 01709 ! fix minimum runup depth 01710 if (par%nrugauge>0) then 01711 ! Fill up remaining part of the array with minimum value 01712 par%rugdepth(par%nrugdepth+1:) = (1.0d0 + epsilon(0.d0))*par%eps 01713 if (minval(par%rugdepth)<=par%eps) then 01714 where(par%rugdepth<=par%eps) 01715 par%rugdepth = (1.0d0 + epsilon(0.d0))*par%eps 01716 endwhere 01717 call writelog('lws','(a,f0.5,a)','Warning: Setting rugdepth to minimum value greater than eps (', & 01718 (1.0d0 + epsilon(0.d0))*par%eps,')') 01719 endif 01720 endif 01721 ! 01722 ! Give an error if you ask for netcdf output if you don't have a netcdf executable 01723 #ifndef USENETCDF 01724 if (par%outputformat .eq. OUTPUTFORMAT_NETCDF) then 01725 call writelog('lse', '', 'Error: You have asked for netcdf output [outputformat=netcdf] but this') 01726 call writelog('lse', '', ' executable is built without netcdf support. Use a netcdf enabled') 01727 call writelog('lse', '', ' executable or outputformat=fortran') 01728 call halt_program 01729 endif 01730 #endif 01731 ! 01732 ! Warm-beam scheme does not work for stationary wave model 01733 if (par%wavemodel==WAVEMODEL_STATIONARY) then 01734 if (par%scheme==SCHEME_LAX_WENDROFF) then 01735 par%scheme=SCHEME_UPWIND_2 01736 call writelog('lws','','Warning: Lax Wendroff [scheme=lax_wendroff] scheme is not supported.') 01737 call writelog('lws','',' Changed to 2nd order upwind [scheme=upwind_2]') 01738 elseif (par%scheme==SCHEME_WARMBEAM) then 01739 par%scheme=SCHEME_UPWIND_2 01740 call writelog('lws','','Warning: Warming and Beam [scheme=warmbeam] scheme is not supported in wave stationary mode.') 01741 call writelog('lws','',' Changed to 2nd order upwind [scheme=upwind_2]') 01742 endif 01743 elseif (par%wavemodel==WAVEMODEL_SURFBEAT .and. par%single_dir==1) then 01744 if (par%scheme==SCHEME_WARMBEAM) then 01745 call writelog('lws','','Warning: Warming and Beam [scheme=warmbeam] scheme is not supported in stationary solver') 01746 call writelog('lws','',' (wave direction) component of single_dir wave model. Wave directions will be solved') 01747 call writelog('lws','',' using 2nd order upwind scheme. Wave group propagation component will be soved using') 01748 call writelog('lws','',' Warming and Beam scheme') 01749 endif 01750 endif 01751 ! 01752 ! Lax-Wendroff not yet supported in curvilinear 01753 if (par%scheme==SCHEME_LAX_WENDROFF) then 01754 par%scheme=SCHEME_WARMBEAM 01755 call writelog('lws','','Warning: Lax Wendroff [scheme=lax_wendroff] scheme is not supported, changed') 01756 call writelog('lws','',' to Warming and Beam [scheme=SCHEME_WARMBEAM]') 01757 endif 01758 ! 01759 01760 ! 01761 ! Wave-current interaction with non-stationary waves still experimental 01762 if (par%wavemodel == WAVEMODEL_STATIONARY .and. par%wci==1) then 01763 call writelog('lws','','Warning: Wave-current interaction with non-stationary waves is still') 01764 call writelog('lws','',' experimental, continue with computation nevertheless') 01765 endif 01766 ! 01767 ! Check for setting Snells law and single_dir 01768 if (par%single_dir == 1 .and. par%snells==1) then 01769 call writelog('lse', '', 'The options ''single_dir = 1'' and ''snells = 1'' are not compatible') 01770 call writelog('lse', '', 'Terminating simulation') 01771 call halt_program 01772 endif 01773 ! 01774 ! 2D absorbing boundary limits to 1D absorbing boundary with 1D 01775 if (par%front==FRONT_ABS_2D .and. par%ny<3) then 01776 call writelog('lws','','Warning: 2D absorbing boundary condition [front=abs_2d] reduces to a') 01777 call writelog('lws','',' 1D absorbing boundary condition [front=abs_1d] in') 01778 call writelog('lws','',' 1D mode [ny=0]') 01779 par%front = FRONT_ABS_1D 01780 endif 01781 if (par%back==BACK_ABS_2D .and. par%ny<3) then 01782 call writelog('lws','','Warning: 2D absorbing boundary condition [back=abs_2d] reduces to a') 01783 call writelog('lws','',' 1D absorbing boundary condition [back=abs_1d] in') 01784 call writelog('lws','',' 1D mode [ny=0]') 01785 par%back = BACK_ABS_1D 01786 endif 01787 ! 01788 ! Mean output time fix 01789 if(par%tintm>(par%tstop-par%tstart) .and. par%nmeanvar>0) then 01790 call writelog('lws','','Warning: ''tintm'' is larger than output duration in the simulation.') 01791 call writelog('lws','',' Setting ''tintm'' = tstop-tstart = ',par%tstop-par%tstart,'s') 01792 par%tintm=par%tstop-par%tstart 01793 endif 01794 ! 01795 ! 01796 ! MPI domains 01797 #ifdef USEMPI 01798 if (par%swave==1) then 01799 if ((par%wavemodel == WAVEMODEL_STATIONARY .or. par%single_dir==1) .and. par%ny > 0) then 01800 ! We need to set to mpiboundary = x to solve the stationary wave model. 01801 ! However, this requires ny>3*xmpi_osize 01802 if (par%mpiboundary .ne. MPIBOUNDARY_MAN) then 01803 if(par%ny<=2*xmpi_size) then 01804 call writelog('ewsl','','This simulation cannot be run in current MPI mode:') 01805 call writelog('ewsl','','The stationary wave solver requires MPI subdivision by "x" (split ny).') 01806 call writelog('ewsl','','The number of subdomains selected to run the model is ',xmpi_size,'.') 01807 call writelog('ewsl','','The total number of grid cells in y (ny) is ',par%ny,'.') 01808 call writelog('ewsl','(a,f0.2,a)','The number of cells per domain is ',dble(par%ny)/xmpi_size, & 01809 ', which is less than the minimum value of 3') 01810 call writelog('ewsl','','If you really (!) know what you''re doing, use "mpiboundary = man" and deal') 01811 call writelog('ewsl','','with any unsatisfactory results') 01812 call halt_program 01813 else 01814 par%mpiboundary=MPIBOUNDARY_X 01815 par%mpiboundary_str='x' 01816 call writelog('wsl','','Changing mpiboundary to "x" for stationary wave model') 01817 endif 01818 else 01819 call writelog('wsl','','Warning: the stationary wave model only works with "mpiboundary=x"!') 01820 endif 01821 endif 01822 endif 01823 #endif 01824 ! 01825 ! 01826 ! fix tint 01827 par%tint = min(par%tintg,par%tintp,par%tintm,par%tinth) 01828 ! 01829 ! 01830 ! All input time frames converted to XBeach hydrodynamic time 01831 if (par%morfacopt==1) then 01832 par%tstart = par%tstart / max(par%morfac,1.d0) 01833 par%tint = par%tint / max(par%morfac,1.d0) 01834 par%tintg = par%tintg / max(par%morfac,1.d0) 01835 par%tintp = par%tintp / max(par%morfac,1.d0) 01836 par%tintm = par%tintm / max(par%morfac,1.d0) 01837 par%tinth = par%tinth / max(par%morfac,1.d0) 01838 par%wavint = par%wavint / max(par%morfac,1.d0) 01839 par%tstop = par%tstop / max(par%morfac,1.d0) 01840 par%morstart= par%morstart / max(par%morfac,1.d0) 01841 par%morstop = par%morstop / max(par%morfac,1.d0) 01842 par%rt = par%rt / max(par%morfac,1.d0) 01843 endif 01844 ! 01845 ! 01846 ! Check bc file length in case of instat = reuse. In this case time is defined by 01847 ! morfacopt, which is not known at earlier stage 01848 if (par%wbctype == WBCTYPE_REUSE) then 01849 if (par%wavemodel==WAVEMODEL_NONH) then 01850 dummystring='nhbcflist.bcf' 01851 call checkbcfilelength(par%tstop,par%wbctype,dummystring,filetype,nonh=.true.) 01852 else 01853 dummystring='ebcflist.bcf' 01854 call checkbcfilelength(par%tstop,par%wbctype,dummystring,filetype) 01855 dummystring='qbcflist.bcf' 01856 call checkbcfilelength(par%tstop,par%wbctype,dummystring,filetype) 01857 endif 01858 endif 01859 ! 01860 ! 01861 ! Check for unknown parameters 01862 if (xmaster) call readkey('params.txt','checkparams',dummystring) 01863 ! 01864 ! 01865 ! Distribute over MPI processes 01866 #ifdef USEMPI 01867 call distribute_par(par) 01868 #endif 01869 ! 01870 ! 01871 ! Write settings to file 01872 !include 'parameters.inc' 01873 !call outputparameters(par) 01874 01875 end subroutine all_input 01876 01877 01878 #ifdef USEMPI 01879 subroutine distribute_par(par) 01880 use mpi 01881 use xmpi_module 01882 implicit none 01883 type(parameters) :: par 01884 integer :: ierror,i,npoints 01885 01886 ! We're sending these over by hand, because intel fortran + vs2008 breaks things... 01887 integer, allocatable :: pointtypes(:) ! [-] Point types (0 = point, 1=rugauge) 01888 double precision, allocatable :: xpointsw(:) ! world x-coordinate of output points 01889 double precision, allocatable :: ypointsw(:) ! world y-coordinate of output points 01890 01891 logical, parameter :: toall = .true. 01892 integer :: parlen 01893 01894 01895 ! 01896 ! distribute parameters 01897 01898 ! This distributes all of the properties of par, including pointers. These point to memory adresses on the master 01899 ! We need to reset these on the non masters 01900 01901 call xmpi_bcast(par%swave,toall) 01902 parlen = int(sizeof(par)) 01903 01904 if (toall) then 01905 call MPI_Bcast(par,parlen,MPI_BYTE,xmpi_imaster,xmpi_ocomm,ierror) 01906 else 01907 call MPI_Bcast(par,parlen,MPI_BYTE,xmpi_master,xmpi_comm,ierror) 01908 endif 01909 01910 ! Ok now for the manual stuff to circumvent a bug in the intel compiler, which doesn't allow to send over arrays in derived types 01911 ! The only way to do it on all 3 compilers (gfortran, CVF, ifort) is with pointers. 01912 ! First let's store the number of variables, we need this to reserve some memory on all nodes 01913 01914 do i=1,size(par%globalvars) 01915 call xmpi_bcast(par%globalvars(i),toall) 01916 enddo 01917 01918 do i=1,size(par%pointvars) 01919 call xmpi_bcast(par%pointvars(i),toall) 01920 enddo 01921 01922 do i=1,size(par%meanvars) 01923 call xmpi_bcast(par%meanvars(i),toall) 01924 enddo 01925 01926 if (xmaster) npoints = size(par%pointtypes) 01927 ! send it over 01928 call xmpi_bcast(npoints,toall) 01929 01930 ! now on all nodes allocate a array outside the par structure 01931 allocate(pointtypes(npoints)) 01932 01933 ! Par is only filled on the master, so use that one and put it in the seperate array 01934 if (xmaster) pointtypes = par%pointtypes 01935 01936 ! Now for another ugly step, we can't broadcast the whole array but have to do it per variable. 01937 call xmpi_bcast(pointtypes,toall) 01938 01939 ! so now everyone has the pointtypes, let's put it back in par 01940 ! first dereference the old one, otherwise we get a nasty error on ifort.... 01941 if (.not. xmaster) par%pointtypes => NULL() 01942 01943 ! now we can allocate the memory again 01944 if (.not. xmaster) allocate(par%pointtypes(npoints)) 01945 01946 ! and store the values from the local array 01947 if (.not. xmaster) par%pointtypes = pointtypes 01948 01949 ! and clean up the local one. 01950 deallocate(pointtypes) 01951 01952 if (xmaster) npoints = size(par%xpointsw) 01953 ! send it over 01954 call xmpi_bcast(npoints,toall) 01955 ! now on all nodes allocate a array outside the par structure 01956 allocate(xpointsw(npoints)) 01957 ! Par is only filled on the master, so use that one and put it in the seperate array 01958 if (xmaster) xpointsw = par%xpointsw 01959 ! Now for another ugly step, we can't broadcast the whole array but have to do it per variable. 01960 call xmpi_bcast(xpointsw,toall) 01961 ! so now everyone has the xpointsw, let's put it back in par 01962 ! first dereference the old one, otherwise we get a nasty error on ifort.... 01963 if (.not. xmaster) par%xpointsw => NULL() 01964 ! now we can allocate the memory again 01965 if (.not. xmaster) allocate(par%xpointsw(npoints)) 01966 ! and store the values from the local array 01967 if (.not. xmaster) par%xpointsw = xpointsw 01968 ! and clean up the local one. 01969 deallocate(xpointsw) 01970 01971 if (xmaster) npoints = size(par%ypointsw) 01972 ! send it over 01973 call xmpi_bcast(npoints,toall) 01974 ! now on all nodes allocate a array outside the par structure 01975 allocate(ypointsw(npoints)) 01976 ! Par is only filled on the master, so use that one and put it in the seperate array 01977 if (xmaster) ypointsw = par%ypointsw 01978 ! Now for another ugly step, we can't broadcast the whole array but have to do it per variable. 01979 call xmpi_bcast(ypointsw,toall) 01980 ! so now everyone has the ypointsw, let's put it back in par 01981 ! first dereference the old one, otherwise we get a nasty error on ifort.... 01982 if (.not. xmaster) par%ypointsw => NULL() 01983 ! now we can allocate the memory again 01984 if (.not. xmaster) allocate(par%ypointsw(npoints)) 01985 ! and store the values from the local array 01986 if (.not. xmaster) par%ypointsw = ypointsw 01987 ! and clean up the local one. 01988 deallocate(ypointsw) 01989 01990 01991 01992 01993 return 01994 01995 ! so, the following code is NOT used anymore. I left this here 01996 ! maybe method above does not work everywhere. wwvv 01997 01998 ! For efficiency, this subroutine should use MPI_Pack and 01999 ! MPI_Unpack. However, since this subroutine is only called a 02000 ! few times, a more simple approach is used. 02001 ! 02002 02003 !call xmpi_bcast(par%px) 02004 !call xmpi_bcast(par%Hrms) 02005 !.... 02006 02007 end subroutine distribute_par 02008 #endif 02009 02010 ! 02011 ! Some extra functions to make reading the output variables possible 02012 ! 02013 subroutine readglobalvars(par) 02014 use logging_module 02015 use mnemmodule, only: mnemonics 02016 implicit none 02017 type(parameters), intent(inout) :: par 02018 integer :: i 02019 02020 if (xmaster) then 02021 if (par%nglobalvar == -1) then 02022 par%globalvars(1:21) = (/'H ', 'zs ', 'zs0 ', 'zb ', 'hh ', 'u ', 'v ', 'ue ',& 02023 've ', 'urms ', 'Fx ', 'Fy ', 'ccg ', 'ceqsg', 'ceqbg', 'Susg ',& 02024 'Svsg ', 'E ', 'R ', 'D ', 'DR ' /) 02025 par%nglobalvar = 21 02026 elseif (par%nglobalvar == 999) then ! Output all 02027 par%nglobalvar = 0 02028 do i=1,numvars 02029 ! Don't output variables that don't work (misalignment) 02030 if (mnemonics(i) .eq. 'xyzs01') then 02031 cycle 02032 elseif (mnemonics(i) .eq. 'xyzs02') then 02033 cycle 02034 elseif (mnemonics(i) .eq. 'xyzs03') then 02035 cycle 02036 elseif (mnemonics(i) .eq. 'xyzs04') then 02037 cycle 02038 elseif (mnemonics(i) .eq. 'tideinpz') then 02039 cycle 02040 !elseif (mnemonics(i) .eq. 'umean') then 02041 ! cycle 02042 !elseif (mnemonics(i) .eq. 'vmean') then 02043 ! cycle 02044 elseif (mnemonics(i) .eq. 'gw0back') then 02045 cycle 02046 elseif (mnemonics(i) .eq. 'zi') then 02047 cycle 02048 elseif (mnemonics(i) .eq. 'wi') then 02049 cycle 02050 elseif (mnemonics(i) .eq. 'tideinpz') then 02051 cycle 02052 end if 02053 par%globalvars(par%nglobalvar+1) = mnemonics(i) ! list of all s% variables 02054 par%nglobalvar = par%nglobalvar + 1 02055 end do 02056 else 02057 ! User specified output 02058 ! Find nglobalvar keyword in params.txt 02059 call readOutputStrings(par,'global') 02060 end if ! globalvar 02061 end if ! xmaster 02062 end subroutine readglobalvars 02063 02064 ! 02065 ! FB: 02066 ! Now for a long one, reading the point and rugauges output options 02067 ! Just moved this from varoutput, so it can be used in ncoutput also 02068 ! Keeping as much as possible local to the subroutine... 02069 ! It can be reduced a bit by combining points and rugauges 02070 ! also instead of nvarpoints it can use a ragged array: 02071 ! For example: http://coding.derkeiler.com/Archive/Fortran/comp.lang.fortran/2004-05/0774.html 02072 ! 02073 ! The outputformat for points is xworld, yworld, nvars, var1#var2# 02074 ! Split for nrugauges and npoints 02075 ! This is a bit inconvenient because you have different sets of points 02076 ! I think it's more logical to get value per variable than per point.... 02077 ! So I read the data as follows: 02078 ! make a collection of all points 02079 ! store the types per point 02080 ! store all found variables in a combined list (not per point) 02081 ! So for each point all variables are stored 02082 ! 02083 02084 ! TODO after some discussion with Robert we will change this to: 02085 ! Rugauges don't need any output variables other then zs, x, y, t 02086 ! Pointvars need a different input specification: 02087 ! npointvars=3 02088 ! zs 02089 ! H 02090 ! u 02091 ! npoints=2 02092 ! 3.0 2.0 02093 ! 10.0 3.1 02094 ! nrugauges=1 02095 ! 4.1 3.2 02096 ! Using the old notation should give a warning and an explanation what will be output. 02097 02098 ! Tasks: 02099 ! Create tests: FB+RMC 02100 ! Implement npointvars, stop using vars in #: RMC 02101 ! Give error if # present in points: RMC 02102 ! Use fixed outputvars for rugauges: RMC 02103 ! Fix matlab read routines (toolbox + zelt): FB 02104 ! Reimplement rugauges in ncoutput: FB 02105 ! Update documentation: FB, check RMC 02106 02107 02108 subroutine readpointvars(par) 02109 use logging_module, only: writelog 02110 use mnemmodule, only: chartoindex, mnem_xz, mnem_yz, mnem_yz, mnem_zs 02111 use readkey_module, only: isSetParameter 02112 implicit none 02113 type(parameters), intent(inout) :: par 02114 02115 double precision, dimension(:),allocatable :: xpointsw,ypointsw 02116 integer, dimension(:), allocatable :: pointtypes 02117 integer :: i,j 02118 logical :: xzfound, yzfound, zsfound 02119 02120 ! These "targets" must be allocated by all processes 02121 allocate(pointtypes(par%npoints+par%nrugauge)) 02122 allocate(xpointsw(par%npoints+par%nrugauge)) 02123 allocate(ypointsw(par%npoints+par%nrugauge)) 02124 if (xmaster) then 02125 ! set the point types to the indicator 02126 pointtypes(1:par%npoints)=0 02127 pointtypes(par%npoints+1:par%npoints+par%nrugauge)=1 02128 par%pointvars='' 02129 if (par%npoints>0) then 02130 call readPointPosition(par,'point',xpointsw,ypointsw) 02131 ! Is the keyword npointvar specified? 02132 if (isSetParameter('params.txt','npointvar',bcast=.false.)) then 02133 ! Look for keyword npointvar in params.txt 02134 call readOutputStrings(par,'point') 02135 else 02136 if(par%nrugauge<=0) then 02137 ! This branch of the else will change to a halt_program statement in later versions (written on 13 January 2011) 02138 call writelog('lswe','','Point output must be specified using keyword ''npointvar''') 02139 call writelog('lswe','','Stopping simulation') 02140 call halt_program 02141 else 02142 call writelog('lsw','','All point output will contain same data as rugauge output (x,y,zs).') 02143 call writelog('lsw','','Other point output must be specified using keyword ''npointvar''') 02144 endif 02145 endif ! isSetParameter 02146 endif ! par%npoints>0 02147 02148 if (par%nrugauge>0) then 02149 call readPointPosition(par,'rugauge',xpointsw,ypointsw) 02150 ! Make sure xz, yz and zs are in pointvars (default) 02151 xzfound = .false. 02152 yzfound = .false. 02153 zsfound = .false. 02154 do i=1,par%npointvar 02155 if (par%pointvars(i) .eq. 'xz') xzfound = .true. 02156 if (par%pointvars(i) .eq. 'yz') yzfound = .true. 02157 if (par%pointvars(i) .eq. 'zs') zsfound = .true. 02158 end do 02159 if (.not.xzfound) then 02160 par%pointvars(par%npointvar+1)='xz' 02161 par%npointvar = par%npointvar + 1 02162 end if 02163 if (.not.yzfound) then 02164 par%pointvars(par%npointvar+1)='yz' 02165 par%npointvar = par%npointvar + 1 02166 end if 02167 if (.not.zsfound) then 02168 par%pointvars(par%npointvar+1)='zs' 02169 par%npointvar = par%npointvar + 1 02170 end if 02171 end if 02172 02173 ! Now check all point and runup gauge names 02174 if(par%nrugauge+par%npoints>0) then 02175 do i=1,par%nrugauge+par%npoints-1 02176 do j=i+1,par%nrugauge+par%npoints 02177 if(par%stationid(i)==par%stationid(j)) then 02178 call writelog('lswe','','Duplicate names used for point station ID:') 02179 call writelog('lswe','',par%stationid(i)) 02180 call writelog('lswe','','Stopping simulation') 02181 call halt_program 02182 endif 02183 enddo 02184 enddo 02185 endif 02186 02187 ! Set pointers correct 02188 allocate(par%pointtypes(size(pointtypes))) 02189 allocate(par%xpointsw(size(xpointsw))) 02190 allocate(par%ypointsw(size(ypointsw))) 02191 par%pointtypes=pointtypes 02192 par%xpointsw= xpointsw 02193 par%ypointsw=ypointsw 02194 endif ! xmaster 02195 ! 02196 end subroutine readpointvars 02197 02198 02199 subroutine readmeans(par) 02200 use logging_module 02201 use mnemmodule 02202 02203 implicit none 02204 02205 type(parameters),intent(inout) :: par 02206 02207 if (par%nmeanvar>0) then 02208 call readOutputStrings(par,'mean') 02209 endif 02210 end subroutine readmeans 02211 02212 subroutine readOutputStrings(par,readtype) 02213 use logging_module, only: writelog, report_file_read_error 02214 use mnemmodule, only: chartoindex 02215 use readkey_module 02216 implicit none 02217 type(parameters), intent(inout) :: par 02218 character(*), intent(in) :: readtype 02219 02220 character(slen) :: okline,errline 02221 character(slen) :: line,keyword,keyread 02222 character(slen+slen) :: tempout 02223 integer :: i,imax,id,ic,index,ier 02224 character(maxnamelen),dimension(numvars) :: tempnames 02225 02226 imax = -123 02227 select case (readtype) 02228 case ('global') 02229 imax = par%nglobalvar 02230 keyword = 'nglobalvar' 02231 okline = 'nglobalvar: Will generate global output for variable: ' 02232 errline = ' Unknown global output variable: ''' 02233 case ('mean') 02234 imax = par%nmeanvar 02235 keyword = 'nmeanvar' 02236 okline = 'nmeanvar: Will generate mean, min, max and variance output for variable: ' 02237 errline = ' Unknown mean output variable: ''' 02238 case ('point') 02239 imax = par%npointvar 02240 keyword = 'npointvar' 02241 okline = 'npointvar: Will generate point output for variable: ' 02242 errline = ' Unknown global output variable: ''' 02243 case ('rugauge') 02244 imax = 0 02245 keyword = '' 02246 okline = '' 02247 errline = '' 02248 case default 02249 write(*,*)'Programming error calling readOutputStrings' 02250 write(*,*)'Unknown calling type ''',trim(readtype),'''' 02251 call halt_program 02252 end select 02253 tempnames='' 02254 02255 if (xmaster) then 02256 id=0 02257 open(10,file='params.txt') ! (this is done by xmaster only) 02258 do while (id == 0) 02259 read(10,'(a)',iostat=ier)line 02260 if (ier .ne. 0) then 02261 tempout = 'params.txt (looking for '//trim(keyword)//')' 02262 call report_file_read_error(tempout) 02263 endif 02264 ic=scan(line,'=') 02265 if (ic>0) then 02266 keyread=adjustl(line(1:ic-1)) 02267 if (keyread == keyword) id=1 02268 endif 02269 enddo 02270 ! Read through the variables lines, 02271 do i=1,imax 02272 read(10,'(a)',iostat=ier)line 02273 if (ier .ne. 0) then 02274 tempout = 'params.txt (reading '//trim(keyword)//')' 02275 call report_file_read_error(tempout) 02276 endif 02277 line = line 02278 ! Check if this is a valid variable name 02279 index = chartoindex(line) 02280 if (index/=-1) then 02281 tempnames(i)=trim(line) ! wwvv use trim() to avoid compiler warning 02282 call writelog('ls','',trim(okline),trim(tempnames(i))) 02283 else 02284 call writelog('sle','',trim(errline),trim(line),'''') 02285 call halt_program 02286 endif 02287 end do 02288 close(10) 02289 endif 02290 02291 ! only useful information on xmaster, but distributed later by distribute_pars 02292 select case (readtype) 02293 case ('global') 02294 par%globalvars=tempnames 02295 case ('mean') 02296 par%meanvars=tempnames 02297 case ('point') 02298 par%pointvars=tempnames 02299 call writelog('ls','','Order of point output variables stored in ''pointvars.idx''') 02300 case ('rugauge') 02301 ! no variables to store 02302 end select 02303 02304 end subroutine readOutputStrings 02305 02306 subroutine readPointPosition(par,readtype,xpoints,ypoints) 02307 use logging_module, only: writelog 02308 use mnemmodule 02309 use readkey_module, only: strippedline 02310 implicit none 02311 type(parameters), intent(inout) :: par 02312 character(*), intent(in) :: readtype 02313 double precision, dimension(:),intent(inout) :: xpoints,ypoints 02314 02315 character(slen) :: line,keyword,keyread,varline 02316 character(maxnamelen) :: varstr 02317 character(slen) :: fullline,errmes1,errmes2,okaymes 02318 integer :: i,imax,id,ic,imark,imarkold,imin,nvar,ivar,index,j,ier,ier2 02319 logical :: varfound,readerror 02320 02321 imin = 0 02322 imax = 0 02323 select case (readtype) 02324 case('point') 02325 imin = 0 02326 imax = par%npoints 02327 keyword = 'npoints' 02328 errmes1 = ' points ' 02329 errmes2 = 'State point output variables using the "npointvar" keyword' 02330 okaymes = ' Output point ' 02331 case('rugauge') 02332 imin = par%npoints 02333 imax = par%nrugauge 02334 keyword = 'nrugauge' 02335 errmes1 = ' runup gauge ' 02336 errmes2 = 'Runup gauge automatically returns t,xz,yz and zs only' 02337 okaymes = ' Output runup gauge ' 02338 case default 02339 write(*,*)'Programming error calling readOutputStrings' 02340 write(*,*)'Unknown calling type ''',trim(readtype),'''' 02341 call halt_program 02342 end select 02343 02344 if (xmaster) then 02345 id=0 02346 open(10,file='params.txt') ! (this is done by xmaster only) 02347 do while (id == 0) 02348 read(10,'(a)')line 02349 ic=scan(line,'=') 02350 if (ic>0) then 02351 keyread=adjustl(line(1:ic-1)) 02352 if (keyread == keyword) id=1 02353 endif 02354 enddo 02355 ! Read positions 02356 do i=1,imax 02357 read(10,'(a)')fullline 02358 ! remove tab characters 02359 fullline = strippedline(fullline) 02360 ! Check params.txt has old (unsupported) method of defining points 02361 ic=scan(fullline,'#') 02362 if (ic .ne. 0) then ! This branch of the if/else will change to a halt_program statement in later versions 02363 ! (written on 13 January 2011) 02364 call writelog('lswe','','Error in definition of point output.') 02365 call writelog('lswe','','Use of #var1#var2#etc. is no longer valid') 02366 call writelog('lswe','','Stopping simulation') 02367 call halt_program 02368 else ! not old method of setting point output 02369 read(fullline,*,iostat=ier)xpoints(i+imin),ypoints(i+imin),par%stationid(i+imin) 02370 ! error checking 02371 if(ier==0) then 02372 ! all fine 02373 readerror=.false. 02374 elseif(ier==-1) then 02375 ! line is too short, probably missing name of station 02376 ! try reading just the coordinates 02377 read(fullline,*,iostat=ier2)xpoints(i+imin),ypoints(i+imin) 02378 if (ier2==0) then 02379 ! coordinates okay, just name missing 02380 select case (readtype) 02381 case('point') 02382 write(par%stationid(i+imin),'("point",i0.3)') i 02383 case('rugauge') 02384 write(par%stationid(i+imin),'("rugau",i0.3)') i 02385 end select 02386 readerror=.false. 02387 else 02388 readerror=.true. 02389 endif 02390 else 02391 readerror=.true. 02392 endif 02393 if(readerror) then 02394 ! Error reading point/rugauge input. Stop 02395 select case (readtype) 02396 case('point') 02397 call writelog('lswe','','Error reading output point location/name in the following line in params.txt:') 02398 case('rugauge') 02399 call writelog('lswe','','Error reading runup gauge location/name in the following line in params.txt:') 02400 end select 02401 call writelog('lswe','',trim(fullline)) 02402 call writelog('lswe','','Stopping simulation') 02403 call halt_program 02404 else 02405 call writelog('ls','(a,a,a,f0.2,a,f0.2)',& 02406 trim(okaymes),trim(par%stationid(i+imin)),' xpoint: ',& 02407 xpoints(i+imin),' ypoint: ',ypoints(i+imin)) 02408 endif 02409 endif ! old method of point output 02410 enddo 02411 close(10) 02412 endif 02413 02414 end subroutine readPointPosition 02415 02416 subroutine check_instat_backward_compatibility(par) 02417 use logging_module 02418 use mnemmodule 02419 use readkey_module 02420 02421 implicit none 02422 02423 type(parameters), intent(inout) :: par 02424 02425 if (isSetParameter('params.txt','nonh') .or. isSetParameter('params.txt','instat')) then 02426 call writelog('l','','--------------------------------') 02427 call writelog('l','','Backward compatibility:') 02428 endif 02429 ! 02430 ! Part 1: nonh still needs to work 02431 if (isSetParameter('params.txt','nonh') .or. par%useXBeachGSettings==1) then 02432 call writelog('ws','(a,a,a)','Warning: Specification of nonh using parameter ''wavemodel''') 02433 if (par%useXBeachGSettings==1) then 02434 par%nonh = readkey_int ('params.txt','nonh', 1, 0, 1) 02435 else 02436 par%nonh = readkey_int ('params.txt','nonh', 0, 0, 1) 02437 endif 02438 if (par%nonh==1) then 02439 par%wavemodel = WAVEMODEL_NONH 02440 par%wavemodel_str = 'nonh' 02441 endif 02442 else 02443 endif 02444 ! 02445 ! Part 2:instat still needs to work and defines both wavemodel and wbctypef 02446 if (isSetParameter('params.txt','instat')) then 02447 ! 02448 ! Write warning message 02449 call writelog('ws','(a,a,a)','Warning: Specification of instat using parameter ''wbctype''') 02450 ! 02451 ! Read instat anyway 02452 call setallowednames('stat', INSTAT_STAT, & 02453 'bichrom', INSTAT_BICHROM, & 02454 'ts_1', INSTAT_TS_1, & 02455 'ts_2', INSTAT_TS_2, & 02456 'jons', INSTAT_JONS, & 02457 'swan', INSTAT_SWAN, & 02458 'vardens', INSTAT_VARDENS, & 02459 'reuse', INSTAT_REUSE, & 02460 'ts_nonh', INSTAT_TS_NONH, & 02461 'off', INSTAT_OFF, & 02462 'stat_table', INSTAT_STAT_TABLE, & 02463 'jons_table', INSTAT_JONS_TABLE) 02464 call setoldnames('0','1','2','3','4','5','6','7','8','9','40','41') 02465 call parmapply('instat',2,par%instat,par%instat_str) 02466 ! 02467 ! Conditions without spectra (e.g. stat, ts_1, etc.) 02468 ! 1) INSTAT = STAT 02469 if (par%instat==INSTAT_STAT) then 02470 if (par%nonh==1) then 02471 par%wavemodel = WAVEMODEL_NONH 02472 par%wavemodel_str = 'nonh' 02473 else 02474 par%wavemodel = WAVEMODEL_STATIONARY 02475 par%wavemodel_str = 'stationary' 02476 endif 02477 par%wbctype = WBCTYPE_PARAMS 02478 par%wbctype_str = 'params' 02479 endif 02480 ! 2) INSTAT = TS_1 02481 if (par%instat==INSTAT_TS_1) then 02482 par%wavemodel = WAVEMODEL_SURFBEAT 02483 par%wavemodel_str = 'surfbeat' 02484 par%wbctype = WBCTYPE_TS_1 02485 par%wbctype_str = 'ts_1' 02486 endif 02487 ! 3) INSTAT = TS_2 02488 if (par%instat==INSTAT_TS_2) then 02489 par%wavemodel = WAVEMODEL_SURFBEAT 02490 par%wavemodel_str = 'surfbeat' 02491 par%wbctype = WBCTYPE_TS_2 02492 par%wbctype_str = 'ts_2' 02493 endif 02494 ! 4) INSTAT = OFF 02495 if (par%instat==INSTAT_OFF .and. par%wavemodel == WAVEMODEL_NONH) then 02496 par%wavemodel = WAVEMODEL_NONH 02497 par%wavemodel_str = 'nonh' 02498 par%wbctype = WBCTYPE_OFF 02499 par%wbctype_str = 'off' 02500 elseif (par%instat==INSTAT_OFF) then 02501 par%wavemodel = WAVEMODEL_SURFBEAT 02502 par%wavemodel_str = 'surfbeat' 02503 par%wbctype = WBCTYPE_OFF 02504 par%wbctype_str = 'off' 02505 endif 02506 ! 5) INSTAT = BICHROM 02507 if (par%instat==INSTAT_BICHROM) then 02508 par%wavemodel = WAVEMODEL_SURFBEAT 02509 par%wavemodel_str = 'surfbeat' 02510 par%wbctype = WBCTYPE_PARAMS 02511 par%wbctype_str = 'params' 02512 ! But we should define the Tlong 02513 if (isSetParameter('params.txt','Tlong')) then 02514 ! we will read it in later 02515 else 02516 par%Tlong = 80 02517 endif 02518 endif 02519 ! 6) INSTAT = STAT_TABLE 02520 if (par%instat==INSTAT_STAT_TABLE) then 02521 par%wavemodel = WAVEMODEL_STATIONARY 02522 par%wavemodel_str = 'stationary' 02523 par%wbctype = WBCTYPE_JONS_TABLE 02524 par%wbctype_str = 'jons_table' 02525 endif 02526 ! 7) INSTAT = ts_nonh 02527 if (par%instat==INSTAT_TS_NONH) then 02528 par%wavemodel = WAVEMODEL_NONH 02529 par%wavemodel_str = 'nonh' 02530 par%wbctype = WBCTYPE_TS_NONH 02531 par%wbctype_str = 'ts_nonh' 02532 endif 02533 ! 02534 ! Conditions with spectra 02535 ! 1) INSTAT_JONS 02536 if (par%instat==INSTAT_JONS .and. par%wavemodel == WAVEMODEL_NONH) then 02537 par%wavemodel = WAVEMODEL_NONH 02538 par%wavemodel_str = 'nonh' 02539 par%wbctype = WBCTYPE_PARAMETRIC 02540 par%wbctype_str = 'parametric' 02541 elseif (par%instat==INSTAT_JONS) then 02542 par%wavemodel = WAVEMODEL_SURFBEAT 02543 par%wavemodel_str = 'surfbeat' 02544 par%wbctype = WBCTYPE_PARAMETRIC 02545 par%wbctype_str = 'parametric' 02546 endif 02547 ! 02548 ! 2) INSTAT_JONS_TABLE 02549 if (par%instat==INSTAT_JONS_TABLE .and. par%wavemodel == WAVEMODEL_NONH) then 02550 par%wavemodel = WAVEMODEL_NONH 02551 par%wavemodel_str = 'nonh' 02552 par%wbctype = WBCTYPE_JONS_TABLE 02553 par%wbctype_str = 'jons_table' 02554 elseif (par%instat==INSTAT_JONS_TABLE) then 02555 par%wavemodel = WAVEMODEL_SURFBEAT 02556 par%wavemodel_str = 'surfbeat' 02557 par%wbctype = WBCTYPE_JONS_TABLE 02558 par%wbctype_str = 'jons_table' 02559 endif 02560 ! 02561 ! 3) INSTAT = SWAN 02562 if (par%instat==INSTAT_SWAN .and. par%wavemodel == WAVEMODEL_NONH) then 02563 par%wavemodel = WAVEMODEL_NONH 02564 par%wavemodel_str = 'nonh' 02565 par%wbctype = WBCTYPE_SWAN 02566 par%wbctype_str = 'swan' 02567 elseif (par%instat==INSTAT_SWAN) then 02568 par%wavemodel = WAVEMODEL_SURFBEAT 02569 par%wavemodel_str = 'surfbeat' 02570 par%wbctype = WBCTYPE_SWAN 02571 par%wbctype_str = 'swan' 02572 endif 02573 ! 02574 ! 4) INSTAT = VARDENS 02575 if (par%instat==INSTAT_VARDENS .and. par%wavemodel == WAVEMODEL_NONH) then 02576 par%wavemodel = WAVEMODEL_NONH 02577 par%wavemodel_str = 'nonh' 02578 par%wbctype = WBCTYPE_VARDENS 02579 par%wbctype_str = 'vardens' 02580 elseif (par%instat==INSTAT_VARDENS) then 02581 par%wavemodel = WAVEMODEL_SURFBEAT 02582 par%wavemodel_str = 'surfbeat' 02583 par%wbctype = WBCTYPE_VARDENS 02584 par%wbctype_str = 'vardens' 02585 endif 02586 ! 02587 ! Other: re-use 02588 if (par%instat==INSTAT_REUSE .and. par%wavemodel == WAVEMODEL_NONH) then 02589 par%wavemodel = WAVEMODEL_NONH 02590 par%wavemodel_str = 'nonh' 02591 par%wbctype = WBCTYPE_REUSE 02592 par%wbctype_str = 'reuse' 02593 elseif (par%instat==INSTAT_REUSE) then 02594 par%wavemodel = WAVEMODEL_SURFBEAT 02595 par%wavemodel_str = 'surfbeat' 02596 par%wbctype = WBCTYPE_REUSE 02597 par%wbctype_str = 'reuse' 02598 endif 02599 endif 02600 end subroutine 02601 02602 end module params
1.7.6.1
