|
XBeach
|
00001 ! Everything in this module should be handled by xmaster, so only call using xmaster 00002 module spectral_wave_bc_module 00003 use typesandkinds, only: slen 00004 use constants, only: pi, compi 00005 use paramsconst 00006 implicit none 00007 save 00008 00009 private 00010 public lastwaveelevation, n_index_loc, bcfiles, bccount, nspectra, reuseall, spectrumendtime 00011 public spectral_wave_bc 00012 00013 type spectrum ! These are related to input spectra 00014 real*8,dimension(:,:),pointer :: S ! 2D variance density spectrum 00015 real*8,dimension(:),pointer :: f,ang ! 1D frequency and direction vectors for S 00016 real*8,dimension(:),pointer :: Sf ! S integrated over all directions 00017 real*8,dimension(:),pointer :: Sd ! S integrated over all frequencies 00018 integer :: nf,nang ! number of frequencies and angles 00019 real*8 :: df,dang ! frequency and angle step size 00020 real*8 :: hm0,fp,dir0,scoeff ! imposed significant wave height, peak frequency, 00021 ! main wave angle and spreading coefficient 00022 real*8 :: trep,dirm ! representative period and mean wave direction 00023 endtype spectrum 00024 type shortspectrum 00025 real*8,dimension(:),pointer :: Sf ! S integrated over all directions 00026 ! real*8,dimension(:),pointer :: Srep ! S at frequency/direction locations of wave train components 00027 ! real*8 :: dangrep ! representative angle step size for Srep 00028 endtype shortspectrum 00029 type waveparamsnew ! These are related to the generated bc series 00030 real*8 :: h0 ! average water depth on offshore boundary 00031 integer :: K ! number of wave train components 00032 real*8 :: rtbc,dtbc ! duration and time step to be written to boundary condition file 00033 real*8 :: rtin,dtin ! duration and time step for the internal time axis, based on Fourier 00034 ! limitations and the number of wave train components (wp%K) 00035 real*8,dimension(:),pointer :: tin ! internal time axis 00036 real*8,dimension(:),pointer :: taperf,taperw ! internal taper function for flow and waves 00037 integer :: tslen ! internal time axis length (is always even, not odd) 00038 integer :: tslenbc ! time axis length for boundary condition file 00039 logical :: dtchanged ! quick check to see if dtin == dtbc (useful for interpolation purpose) 00040 real*8,dimension(:),pointer :: fgen,thetagen,phigen,kgen,wgen ! frequency, angle, phase, 00041 ! wave number, radian frequency 00042 ! of wavetrain components for boundary signal 00043 real*8 :: dfgen ! frequency grid size in the generated components 00044 type(shortspectrum),dimension(:),pointer :: vargen ! This is the variance for each wave train at each spectrum location 00045 ! is stored 00046 type(shortspectrum),dimension(:),pointer :: vargenq ! This is the variance for each wave train at each spectrum location 00047 ! is stored, which is not scaled and is used for the generation of bound waves 00048 !real*8,dimension(:),pointer :: danggen ! Representative integration angle for Srep to Sf, per offshore grid point 00049 real*8,dimension(:,:),pointer :: A ! Amplitude, per wave train component, per offshore grid point A(ny+1,K) 00050 real*8,dimension(:,:),pointer :: Sfinterp ! S integrated over all directions at frequency locations of fgen, 00051 ! per offshore grid point Sfinterp(ny+1,K) 00052 real*8,dimension(:,:),pointer :: Sfinterpq ! S integrated over all directions at frequency locations of fgen, 00053 ! per offshore grid point Sfinterpq(ny+1,K), uncorrected for generation of bound waves 00054 real*8,dimension(:),pointer :: Hm0interp ! Hm0 per offshore point, based on intergration of Sfinterp. 00055 ! Used to scale spectrum 00056 ! final time series 00057 integer, dimension(:),pointer :: Findex ! Index of wave train component locations on frequency/Fourier axis 00058 integer, dimension(:),pointer :: WDindex ! Index of wave train component locations on wave directional bin axis 00059 integer, dimension(:),pointer :: PRindex ! Index of wave train components to be phase-resolved (rather than 00060 ! using energy balance) 00061 double complex,dimension(:,:),pointer :: CompFn ! Fourier components of the wave trains 00062 character(slen) :: Efilename,qfilename,nhfilename,Esfilename 00063 real*8,dimension(:,:),pointer :: zsits ! time series of total surface elevation for nonhspectrum==1 00064 real*8,dimension(:,:),pointer :: uits,vits ! time series of depth-averaged horizontal velocity nonhspectrum==1 00065 real*8,dimension(:,:),pointer :: duits,dvits! time series of depth-averaged horizontal velocity nonhspectrum==1 00066 00067 real*8,dimension(:,:),pointer :: wits ! time series of depth-averaged vertical velocity for nonhspectrum==1 ?? 00068 endtype waveparamsnew 00069 type filenames ! Place to store multiple file names 00070 character(slen) :: fname ! file name of boundary condition file 00071 integer :: listline ! read position in FILELIST files 00072 logical :: reuse = .false. ! indicate to reuse this file every rtbc cycle 00073 endtype filenames 00074 00075 ! These are for administration purposes and are initialized in initialize.F90 00076 integer,dimension(:),allocatable :: n_index_loc ! y-index locations of all input spectra, set in init spectrum 00077 integer :: nspectra ! number of input spectrs, set in init spectrum 00078 type(filenames),dimension(:),allocatable :: bcfiles ! input wave spectrum files 00079 logical :: reuseall ! switch to reuse all of the wave boundary conditions 00080 integer,save :: bccount ! number of times boundary conditions have been generated, set in init spectrum 00081 integer,save :: fidelist,fidqlist,fidnhlist,fideslist ! file identifiers for ebcflist and qbcflist 00082 real*8 :: spectrumendtime ! end time of boundary condition written to administration file 00083 real*8,dimension(:,:),allocatable :: lastwaveelevation ! wave height at the end of the last spectrum 00084 integer :: ind_end_taper ! index of where the taper function equals rtbc 00085 ! These parameters control a lot how the spectra are handled. They could be put in params.txt, 00086 ! but most users will want to keep these at their default values anyway 00087 integer,parameter,private :: nfint = 801 ! size of standard 2D spectrum in frequency dimension 00088 integer,parameter,private :: naint = 401 ! size of standard 2D spectrum in angular dimension 00089 integer,parameter,private :: Kmin = 200 ! minimum number of wave train components 00090 real*8,parameter,private :: wdmax = 5.d0 ! maximum depth*reliable angular wave frequency that can be resolved by 00091 ! nonhydrostatic wave model. All frequencies above this are removed 00092 ! from nonhspectrum generation 00093 ! Constants, cannot be modified 00094 real*8,parameter,private :: par_pi = 4.d0*atan(1.d0) 00095 00096 contains 00097 00098 subroutine spectral_wave_bc(s,par,curline) 00099 use params, only: parameters 00100 use spaceparamsdef, only: spacepars 00101 use logging_module, only: writelog 00102 use filefunctions, only: create_new_fid 00103 use xmpi_module 00104 00105 implicit none 00106 type(spacepars),intent(inout) :: s 00107 type(parameters),intent(inout) :: par 00108 integer,intent(out) :: curline 00109 type(spectrum),dimension(:),allocatable :: specin,specinterp 00110 type(spectrum) :: combspec 00111 type(waveparamsnew) :: wp ! Most will be deallocated, but some variables useful to keep? 00112 integer :: iloc 00113 integer :: iostat 00114 real*8 :: spectrumendtimeold,fmax 00115 real*8,dimension(:),allocatable :: spectrumendtimearray 00116 real*8,save :: rtbc_local,dtbc_local 00117 real*8,save :: maindir_local 00118 00119 spectrumendtimeold = -123 00120 00121 !fidqlist = -123 00122 !fidnhlist = -123 00123 !fidelist = -123 00124 ! Update the number of boundary conditions generated by this module 00125 bccount = bccount+1 00126 00127 ! Arnold: check if this is the last computational timestep. If so, don't recalculate 00128 ! boundary conditions 00129 if (par%t>(par%tstop-par%dt)) then 00130 return 00131 end if 00132 00133 ! wwvv copied call to set_bcfilenames to here: (see later) 00134 call set_bcfilenames(wp) 00135 if (xmaster) then 00136 if (reuseall) then 00137 ! Wave boundary conditions have already been computed 00138 ! Filename, rtbc, dtbc, etc. still kept in memory. 00139 00140 ! The end time for the boundary conditions calculated in the last call 00141 ! is equal to the old end time, plus the duration of the boundary conditions 00142 ! created at the last call 00143 spectrumendtimeold = spectrumendtime ! Robert: need this in case new run with 00144 ! instat='reuse' does not have exact same 00145 ! time steps as original simulation 00146 spectrumendtime = spectrumendtime + rtbc_local 00147 00148 else 00149 call writelog('l','','--------------------------------') 00150 call writelog('ls','','Calculating spectral wave boundary conditions ') 00151 call writelog('l','','--------------------------------') 00152 00153 ! allocate temporary input storage 00154 if (.not. allocated(specin)) then 00155 allocate(specin(nspectra)) 00156 allocate(specinterp(nspectra)) 00157 allocate(spectrumendtimearray(nspectra)) 00158 endif 00159 00160 ! calculate the average offshore water depth, which is used in various 00161 ! computation in this module 00162 wp%h0 = calculate_average_water_depth(s) 00163 00164 ! Read through input spectra files 00165 fmax = 1.d0 ! assume 1Hz as maximum frequency. Increase in loop below if needed. 00166 spectrumendtimearray = spectrumendtime 00167 do iloc = 1,nspectra 00168 call writelog('sl','(a,i0)','Reading spectrum at location ',iloc) 00169 ! read input files until endtime exceeds current time. 00170 ! necessary in case of external time management 00171 00172 do while (par%t .ge. spectrumendtimearray(iloc)) 00173 call read_spectrum_input(par,wp,bcfiles(iloc),specin(iloc)) 00174 00175 ! update end time 00176 if (iloc==1) then 00177 spectrumendtimeold = spectrumendtime 00178 spectrumendtime = spectrumendtime + wp%rtbc 00179 endif 00180 spectrumendtimearray(iloc) = spectrumendtimearray(iloc)+wp%rtbc 00181 end do 00182 fmax = max(fmax,maxval(specin(iloc)%f)) 00183 enddo 00184 do iloc = 1,nspectra 00185 call writelog('sl','(a,i0)','Interpreting spectrum at location ',iloc) 00186 ! Interpolate input 2D spectrum to standard 2D spectrum 00187 call interpolate_spectrum(specin(iloc),specinterp(iloc),par,fmax) 00188 call writelog('sl','','Values calculated from interpolated spectrum:') 00189 call writelog('sl','(a,f0.2,a)','Hm0 = ',specinterp(iloc)%hm0,' m') 00190 call writelog('sl','(a,f0.2,a)','Trep = ',specinterp(iloc)%trep,' s') 00191 call writelog('sl','(a,f0.2,a)','Mean dir = ',specinterp(iloc)%dirm,' degN') 00192 enddo 00193 00194 ! Determine whether all the spectra are to be reused, which implies that the global reuseall should be 00195 ! set to true (no further computations required in future calls) 00196 call set_reuseall 00197 00198 ! Determine the file names of the output boundary condition time series 00199 ! wwvv copied call to set_bcfilenames to start of subroutine. 00200 ! to make sure that filenames are set, regardless the value of reuseall 00201 call set_bcfilenames(wp) 00202 00203 ! calculate the mean combined spectra (used for combined Trep, determination of wave components, etc.) 00204 ! now still uses simple averaging, but could be improved to use weighting for distance etc. 00205 call generate_combined_spectrum(specinterp,combspec,par%px,par%trepfac,par%Tm01switch) 00206 par%Trep = combspec%trep 00207 call writelog('sl','(a,f0.2,a)','Overall Trep from all spectra calculated: ',par%Trep,' s') 00208 00209 if (par%single_dir==1) then 00210 call set_stationary_spectrum (s,par,wp,combspec) 00211 endif 00212 00213 ! Wave trains that are used by XBeach. The number of wave trains, their frequencies and directions 00214 ! are based on the combined spectra of all the locations to ensure all wave conditions are 00215 ! represented in the XBeach model 00216 ! call generate_wavetrain_components(combspec,wp,par) 00217 ! max statement below needed for compiler in case of 1D simulation without cyclic boundaries 00218 call generate_wavetrain_components(combspec,wp,par,s%xz(1,1),s%yz(1,1),s%xz(1,max(s%ny-2,1)),s%yz(1,max(s%ny-2,1))) 00219 00220 ! We can now apply a correction to the wave train components if necessary. This section can be 00221 ! improved later 00222 if (nspectra==1 .and. specin(1)%scoeff>1000.d0) then 00223 ! this can be used both for Jonswap and vardens input 00224 wp%thetagen = specin(1)%dir0 00225 endif 00226 00227 ! Set up time axis, including the time axis for output to boundary condition files and an 00228 ! internal time axis, which may differ in length to the output time axis 00229 call generate_wave_time_axis(par,wp) 00230 00231 ! Determine the variance for each wave train component, at every spectrum location point 00232 call generate_wave_train_variance(wp,specinterp) 00233 00234 ! Determine the amplitude of each wave train component, at every point along the 00235 ! offshore boundary 00236 call generate_wave_train_properties_per_offshore_point(wp,s,par%nonhspectrum,par%swkhmin) 00237 00238 ! Generate Fourier components for all wave train component, at every point along 00239 ! the offshore boundary 00240 call generate_wave_train_Fourier(wp,s,par) 00241 00242 ! time series of short wave energy or surface elevation 00243 if (par%nonhspectrum==0) then 00244 ! Distribute all wave train components among the wave direction bins. Also rearrage 00245 ! the randomly drawn wave directions to match the centres of the wave bins if the 00246 ! user-defined nspr is set on. 00247 call distribute_wave_train_directions(wp,s,par%px,par%nspr,.false.) 00248 00249 ! if we want to send some low-frequency swell waves into the model in the NLSWE then 00250 ! separate here into a mix of wave action balance and NLSWE components 00251 if (par%swkhmin>0.d0) then 00252 ! recalculate Trep 00253 call tpDcalc(sum(wp%Sfinterp,DIM=1)/(s%ny+1)*(1-wp%PRindex),wp%fgen,par%Trep,par%trepfac,par%Tm01switch) 00254 call writelog('sl','','Trep recomputed to account only for components in wave action balance.') 00255 call writelog('sl','(a,f0.2,a)','New Trep in wave action balance: ',par%Trep,' s') 00256 ! 00257 ! Calculate the wave energy envelope per offshore grid point and write to output file 00258 call generate_ebcf(wp,s,par) ! note, this subroutine will account for only non-phase-resolved components 00259 ! Generate time series of surface elevation and horizontal velocity, only for phase-resolved components 00260 call generate_swts(wp,s,par) 00261 else 00262 ! Only do wave action balance stuff 00263 ! 00264 ! Calculate the wave energy envelope per offshore grid point and write to output file 00265 call generate_ebcf(wp,s,par) 00266 endif ! swkhmin>0.d0 00267 else 00268 ! If user set nspr on, then force all short wave components to head along computational 00269 ! x-axis. 00270 call distribute_wave_train_directions(wp,s,par%px,par%nspr,.true.) 00271 ! Generate time series of surface elevation and horizontal velocity 00272 call generate_swts(wp,s,par) 00273 endif 00274 00275 ! Calculate the bound long wave from the wave train components and write to output file 00276 if ((par%nonhspectrum==0) .or. (par%nonhspectrum==1 .and. par%order>1)) then 00277 call generate_qbcf(wp,s,par) 00278 endif 00279 00280 ! Calculate second order bound components. 00281 if (par%nonhspectrum==1 .and. par%highcomp==1 .and. par%order>1) then 00282 call generate_secondorder(par,s, wp) 00283 endif 00284 00285 00286 ! Write non-hydrostatic time series of combined short and long waves if necessary 00287 if (par%nonhspectrum==1) then 00288 call generate_nhtimeseries_file(wp,par) 00289 endif 00290 ! 00291 ! 00292 ! Save key variables to memory, in case of 'reuseall' 00293 rtbc_local = wp%rtbc 00294 dtbc_local = wp%dtbc 00295 maindir_local = combspec%dir0 00296 ! 00297 ! 00298 ! Deallocate a lot of memory 00299 deallocate(specin,specinterp,spectrumendtimearray) 00300 deallocate(wp%tin,wp%taperf,wp%taperw) 00301 deallocate(wp%fgen,wp%thetagen,wp%phigen,wp%kgen,wp%wgen) 00302 deallocate(wp%vargen) 00303 deallocate(wp%vargenq) 00304 deallocate(wp%Sfinterp) 00305 deallocate(wp%Sfinterpq) 00306 deallocate(wp%Hm0interp) 00307 deallocate(wp%A) 00308 deallocate(wp%Findex) 00309 deallocate(wp%CompFn) 00310 deallocate(wp%PRindex) 00311 if (par%nonhspectrum==0) then 00312 deallocate(wp%WDindex) 00313 if (par%swkhmin>0.d0) then 00314 deallocate(wp%zsits) 00315 deallocate(wp%uits) 00316 endif 00317 else 00318 deallocate(wp%zsits) 00319 deallocate(wp%uits) 00320 endif 00321 ! Send message to screen and log 00322 call writelog('l','','--------------------------------') 00323 call writelog('ls','','Spectral wave boundary conditions complete ') 00324 call writelog('l','','--------------------------------') 00325 endif ! reuseall 00326 00327 ! Collect new file identifiers for administration list files 00328 call generate_admin_files(par,wp,rtbc_local,dtbc_local,maindir_local,spectrumendtimeold,spectrumendtime) 00329 00330 ! Set output of the correct line 00331 curline = bccount 00332 endif 00333 00334 end subroutine spectral_wave_bc 00335 00336 ! -------------------------------------------------------------- 00337 ! ---------------- Read input spectra files -------------------- 00338 ! -------------------------------------------------------------- 00339 subroutine read_spectrum_input(par,wp,fn,specin) 00340 use params, only: parameters 00341 use filefunctions, only: create_new_fid 00342 use logging_module, only: report_file_read_error 00343 use paramsconst 00344 00345 implicit none 00346 ! Interface 00347 type(parameters),intent(in) :: par 00348 type(waveparamsnew),intent(inout) :: wp 00349 type(filenames),intent(inout) :: fn 00350 type(spectrum),intent(inout) :: specin 00351 ! internal 00352 integer :: fid 00353 character(8) :: testline 00354 character(slen) :: readfile 00355 logical :: filelist 00356 integer :: i,ier 00357 00358 ! check the first line of the boundary condition file for FILELIST keyword 00359 fid = create_new_fid() 00360 open(fid,file=fn%fname,status='old',form='formatted') 00361 read(fid,*,iostat=ier)testline 00362 if (ier .ne. 0) then 00363 call report_file_read_error(fn%fname) 00364 endif 00365 if (trim(testline)=='FILELIST') then 00366 filelist = .true. 00367 ! move listline off its default position of zero to the first row 00368 if (fn%listline==0) then 00369 fn%listline = 1 00370 endif 00371 fn%reuse = .false. 00372 else 00373 filelist = .false. 00374 if (par%wbctype /= WBCTYPE_JONS_TABLE) then 00375 fn%reuse = .true. 00376 else 00377 fn%reuse = .false. 00378 endif 00379 endif 00380 close(fid) 00381 ! If file has list of spectra, read through the lines to find the correct 00382 ! spectrum file name and rtbc and dtbc. Else the filename and rtbc, dtbc are in 00383 ! params.txt 00384 if (filelist) then 00385 fid = create_new_fid() 00386 open(fid,file=fn%fname,status='old',form='formatted') 00387 do i=1,fn%listline 00388 read(fid,*)testline ! old stuff, not needed anymore 00389 enddo 00390 read(fid,*,iostat=ier)wp%rtbc,wp%dtbc,readfile ! new boundary condition 00391 if (ier .ne. 0) then 00392 call report_file_read_error(fn%fname) 00393 endif 00394 ! we have to adjust this to morphological time, as done in params.txt 00395 if (par%morfacopt==1) then 00396 wp%rtbc = wp%rtbc/max(par%morfac,1.d0) 00397 endif 00398 fn%listline = fn%listline + 1 ! move one on from the last time we opened this file 00399 close(fid) 00400 else 00401 wp%rtbc = par%rt ! already set to morphological time 00402 wp%dtbc = par%dtbc 00403 readfile = fn%fname 00404 endif 00405 00406 ! based on the value of instat, we need to read either Jonswap, Swan or vardens files 00407 ! note: jons_table is also handeled by read_jonswap_file subroutine 00408 !select case (par%instat(1:4)) 00409 select case(par%wbctype) 00410 case (WBCTYPE_PARAMETRIC, WBCTYPE_JONS_TABLE) 00411 ! wp type sent in to receive rtbc and dtbc from jons_table file 00412 ! fn%listline sent in to find correct row in jons_table file 00413 ! pfff.. 00414 call read_jonswap_file(par,wp,readfile,fn%listline,specin) 00415 !case ('swan') 00416 case (WBCTYPE_SWAN) 00417 call read_swan_file(par,readfile,specin) 00418 !case ('vard') 00419 case (WBCTYPE_VARDENS) 00420 call read_vardens_file(par,readfile,specin) 00421 endselect 00422 00423 end subroutine read_spectrum_input 00424 00425 ! -------------------------------------------------------------- 00426 ! ------------------- Read JONSWAP files ----------------------- 00427 ! -------------------------------------------------------------- 00428 subroutine read_jonswap_file(par,wp,readfile,listline,specin) 00429 00430 use readkey_module, only: readkey_int, readkey_dbl, readkey_dblvec, isSetParameter, readkey_intvec, readkey 00431 use params, only: parameters 00432 use logging_module, only: writelog, report_file_read_error 00433 use filefunctions, only: create_new_fid 00434 use wave_functions_module, only: iteratedispersion 00435 use xmpi_module 00436 use paramsconst 00437 00438 IMPLICIT NONE 00439 00440 ! Input / output variables 00441 type(parameters), INTENT(IN) :: par 00442 type(waveparamsnew),intent(inout) :: wp 00443 character(len=*), intent(IN) :: readfile 00444 integer, intent(INOUT) :: listline 00445 type(spectrum),intent(inout) :: specin 00446 00447 ! Internal variables 00448 integer :: i,ii,ier,ip,ind 00449 integer :: nmodal 00450 integer :: fid 00451 integer :: forcepartition = -123 00452 real*8,dimension(:),allocatable :: x, y, Dd, tempdir 00453 real*8,dimension(:),allocatable :: Hm0,fp,gam,mainang,scoeff 00454 integer,dimension(:),allocatable :: tma 00455 real*8 :: dfj, fnyq 00456 real*8 :: Tp 00457 real*8 :: L,L0,sigmatma,k,n 00458 character(len=80) :: dummystring 00459 type(spectrum),dimension(:),allocatable :: multinomalspec 00460 00461 ! First part: read JONSWAP parameter data 00462 00463 ! Check whether spectrum characteristics or table should be used 00464 if (par%wbctype /= WBCTYPE_JONS_TABLE) then 00465 ! Use spectrum characteristics 00466 call writelog('sl','','waveparams: Reading from: ',trim(readfile),' ...') 00467 ! 00468 ! First read if the spectrum is multinodal, and how many peaks there should be 00469 ! 00470 nmodal = readkey_int (readfile, 'nmodal', 1, 1, 4, bcast=.false.) 00471 if (nmodal<1) then 00472 call writelog('lswe','','Error: number of spectral partions may not be less than 1 in ',trim(readfile)) 00473 call halt_program 00474 endif 00475 ! 00476 ! Allocate space for all spectral parameters 00477 ! 00478 allocate(Hm0(nmodal)) 00479 allocate(fp(nmodal)) 00480 allocate(gam(nmodal)) 00481 allocate(mainang(nmodal)) 00482 allocate(scoeff(nmodal)) 00483 allocate(tma(nmodal)) 00484 ! 00485 ! Read the spectral parameters for all spectrum components 00486 ! 00487 ! Wave height (required) 00488 Hm0 = readkey_dblvec(readfile, 'Hm0',nmodal,nmodal, 0.0d0, 0.0d0, 5.0d0, bcast=.false.,required=.true. ) 00489 ! 00490 ! Wave period (required) 00491 ! allow both Tp and fp specification to bring in line with params.txt 00492 if (isSetParameter(readfile,'Tp',bcast=.false.) .and. .not. isSetParameter(readfile,'fp',bcast=.false.)) then 00493 fp = 1.d0/readkey_dblvec(readfile, 'Tp',nmodal,nmodal, 12.5d0, 2.5d0, 20.0d0, bcast=.false.) 00494 elseif (isSetParameter(readfile,'fp',bcast=.false.) .and. .not. isSetParameter(readfile,'Tp',bcast=.false.)) then 00495 fp = readkey_dblvec(readfile, 'fp',nmodal,nmodal, 0.08d0,0.0625d0, 0.4d0, bcast=.false.) 00496 elseif (.not. isSetParameter(readfile,'fp',bcast=.false.) .and. .not. isSetParameter(readfile,'Tp',bcast=.false.)) then 00497 call writelog('lswe','','Error: missing required value for parameter ''Tp'' or ''fp'' in ',trim(readfile)) 00498 call halt_program 00499 else 00500 fp = 1.d0/readkey_dblvec(readfile, 'Tp',nmodal,nmodal, 12.5d0, 2.5d0, 20.0d0, bcast=.false.) 00501 call writelog('lsw','','Warning: selecting to read peak period (Tp) instead of frequency (fp) in ',trim(readfile)) 00502 endif 00503 ! 00504 ! Wave spreading in frequency domain (peakedness) 00505 ! 00506 gam = readkey_dblvec(readfile, 'gammajsp',nmodal,nmodal, 3.3d0, 1.0d0, 5.0d0, bcast=.false.) 00507 ! 00508 ! Wave spreading in directional domain 00509 ! 00510 scoeff = readkey_dblvec(readfile, 's',nmodal,nmodal, 10.0d0, 1.0d0, 1000.0d0, bcast=.false.) 00511 ! 00512 ! TMA 00513 ! 00514 tma = readkey_intvec(readfile, 'tma',nmodal,nmodal, 0, 0, 1, bcast=.false.) 00515 ! 00516 ! Main wave direction 00517 ! allow both mainang and dir0 specification to bring in line with params.txt 00518 if (isSetParameter(readfile,'mainang',bcast=.false.) .and. & 00519 .not. isSetParameter(readfile,'dir0',bcast=.false.)) then 00520 mainang = readkey_dblvec(readfile, 'mainang',nmodal,nmodal, & 00521 270.0d0, 0.0d0, 360.0d0, bcast=.false.) 00522 elseif (isSetParameter(readfile,'dir0',bcast=.false.) .and. & 00523 .not. isSetParameter(readfile,'mainang',bcast=.false.)) then 00524 mainang = readkey_dblvec(readfile, 'dir0',nmodal,nmodal, & 00525 270.0d0, 0.0d0, 360.0d0, bcast=.false.) 00526 elseif (.not. isSetParameter(readfile,'dir0',bcast=.false.) .and. & 00527 .not. isSetParameter(readfile,'mainang',bcast=.false.)) then 00528 mainang = 270.d0 00529 else 00530 mainang = readkey_dblvec(readfile, 'mainang',nmodal,nmodal, 270.0d0, 0.0d0, 360.0d0, bcast=.false.) 00531 call writelog('lsw','','Warning: selecting to read ''mainang'' instead of ''dir0'' in ',trim(readfile)) 00532 endif 00533 ! 00534 ! Nyquist parameters used only in this subroutine 00535 ! are not read individually for each spectrum partition 00536 fnyq = readkey_dbl (readfile, 'fnyq',max(0.3d0,3.d0*maxval(fp)), 0.2d0, 1.0d0, bcast=.false. ) 00537 dfj = readkey_dbl (readfile, 'dfj', fnyq/200, fnyq/1000, fnyq/20, bcast=.false. ) 00538 ! 00539 ! Finally check if XBeach should accept even the most stupid partioning (sets error level to warning 00540 ! level when computing partition overlap 00541 if (nmodal>1) then 00542 forcepartition = readkey_int (readfile, 'forcepartition', 0, 0, 1, bcast=.false.) 00543 endif 00544 ! check for other strange values in this file 00545 call readkey(readfile,'checkparams',dummystring) 00546 else 00547 nmodal = 1 00548 allocate(Hm0(nmodal)) 00549 allocate(fp(nmodal)) 00550 allocate(gam(nmodal)) 00551 allocate(mainang(nmodal)) 00552 allocate(scoeff(nmodal)) 00553 allocate(tma(nmodal)) 00554 ! Use spectrum table 00555 fid = create_new_fid() 00556 call writelog('sl','','waveparams: Reading from table: ',trim(readfile),' ...') 00557 open(fid,file=readfile,status='old',form='formatted') 00558 ! read junk up to the correct line in the file 00559 do i=1,listline 00560 read(fid,*,iostat=ier)dummystring 00561 if (ier .ne. 0) then 00562 call report_file_read_error(readfile) 00563 endif 00564 enddo 00565 read(fid,*,iostat=ier)Hm0(1),Tp,mainang(1),gam(1),scoeff(1),wp%rtbc,wp%dtbc 00566 if (ier .ne. 0) then 00567 call writelog('lswe','','Error reading file ',trim(readfile)) 00568 close(fid) 00569 call halt_program 00570 endif 00571 ! move the line pointer in the file 00572 listline = listline+1 00573 ! convert to morphological time 00574 if (par%morfacopt==1) then 00575 wp%rtbc = wp%rtbc/max(par%morfac,1.d0) 00576 endif 00577 fp(1)=1.d0/Tp 00578 fnyq=3.d0*fp(1) 00579 dfj=fp(1)/50 00580 tma(1) = 0 00581 close(fid) 00582 endif 00583 ! 00584 ! Second part: generate 2D spectrum from input parameters 00585 ! 00586 ! Define number of frequency bins by defining an array of the necessary length 00587 ! using the Nyquist frequency and frequency step size 00588 specin%nf = ceiling((fnyq-dfj)/dfj) 00589 specin%df = dfj 00590 ! 00591 ! Define array with actual eqidistant frequency bins 00592 allocate(specin%f(specin%nf)) 00593 do i=1,specin%nf 00594 specin%f(i)=i*dfj 00595 enddo 00596 ! 00597 ! we need a normalised frequency and variance vector for JONSWAP generation 00598 allocate(x(size(specin%f))) 00599 allocate(y(size(specin%f))) 00600 ! 00601 ! Define 200 directions relative to main angle running from 0 to 2*pi 00602 ! we also need a temporary vector for direction distribution 00603 specin%nang = naint 00604 allocate(tempdir(specin%nang)) 00605 allocate(Dd(specin%nang)) 00606 allocate(specin%ang(specin%nang)) 00607 specin%dang = 2*par%px/(naint-1) 00608 do i=1,specin%nang 00609 specin%ang(i)=(i-1)*specin%dang 00610 enddo 00611 ! 00612 ! Generate density spectrum for each spectrum partition 00613 ! 00614 allocate(multinomalspec(nmodal)) 00615 ! 00616 do ip=1,nmodal 00617 ! relative frequenct vector 00618 x=specin%f/fp(ip) 00619 ! Calculate unscaled and non-directional JONSWAP spectrum using 00620 ! peak-enhancement factor and pre-determined frequency bins 00621 call jonswapgk(x,gam(ip),y) 00622 ! Determine scaled and non-directional JONSWAP spectrum using the JONSWAP 00623 ! characteristics 00624 if (tma(ip)==1) then 00625 do ii=1,specin%nf 00626 L0 = par%g*(1/specin%f(ii))**2/2/par%px ! deep water wave length 00627 L = iteratedispersion(L0,L0,par%px,wp%h0) 00628 if (L<0.d0) then 00629 call writelog('lsw','','No dispersion convergence found for wave train ',i, & 00630 ' in boundary condition generation') 00631 L = -L 00632 endif 00633 k = 2*par%px/L 00634 ! h = wp%h0 00635 n = 0.5d0*(1+k*wp%h0*((1-tanh(k*wp%h0)**2)/(tanh(k*wp%h0)))) 00636 sigmatma = ((1/(2.d0*n))*tanh(k*wp%h0)**2) 00637 y(ii) = y(ii) * sigmatma 00638 enddo 00639 endif 00640 y=(Hm0(ip)/(4.d0*sqrt(sum(y)*dfj)))**2*y 00641 ! Convert main angle from degrees to radians and from nautical convention to 00642 ! internal grid 00643 mainang(ip)=(1.5d0*par%px)-mainang(ip)*par%px/180 00644 ! Make sure the main angle is defined between 0 and 2*pi 00645 do while (mainang(ip)>2*par%px .or. mainang(ip)<0.d0) !Robert en Ap 00646 if (mainang(ip)>2*par%px) then 00647 mainang(ip)=mainang(ip)-2*par%px 00648 elseif (mainang(ip)<0.d0) then 00649 mainang(ip)=mainang(ip)+2*par%px 00650 endif 00651 enddo 00652 ! Convert 200 directions relative to main angle to directions relative to 00653 ! internal grid 00654 ! Bas: apparently division by 2 for cosine law happens already here 00655 tempdir = (specin%ang-mainang(ip))/2 00656 ! Make sure all directions around the main angle are defined between 0 and 2*pi 00657 do while (any(tempdir>2*par%px) .or. any(tempdir<0.d0)) 00658 where (tempdir>2*par%px) 00659 tempdir=tempdir-2*par%px 00660 elsewhere (tempdir<0.d0) 00661 tempdir=tempdir+2*par%px 00662 endwhere 00663 enddo 00664 ! Calculate directional spreading based on cosine law 00665 Dd = dcos(tempdir)**(2*nint(scoeff(ip))) ! Robert: apparently nint is needed here, else MATH error 00666 ! Scale directional spreading to have a surface of unity by dividing by its 00667 ! own surface 00668 Dd = Dd / (sum(Dd)*specin%dang) 00669 ! Define two-dimensional variance density spectrum array and distribute 00670 ! variance density for each frequency over directional bins 00671 allocate(multinomalspec(ip)%S(specin%nf,specin%nang)) 00672 do i=1,specin%nang 00673 do ii=1,specin%nf 00674 multinomalspec(ip)%S(ii,i)=y(ii)*Dd(i) 00675 end do 00676 end do 00677 enddo ! 1,nmodal 00678 do ip=1,nmodal 00679 write(*,*)'Hm0 = ',4*sqrt(sum(multinomalspec(ip)%S)*specin%dang*specin%df) 00680 enddo 00681 ! 00682 ! Combine spectrum partitions so that the total spectrum is correct 00683 ! 00684 if (nmodal==1) then 00685 ! Set all the useful parameters and arrays 00686 specin%Hm0 = Hm0(1) 00687 specin%fp = fp(1) 00688 specin%dir0 = mainang(1) 00689 specin%scoeff = scoeff(1) 00690 allocate(specin%S(specin%nf,specin%nang)) 00691 specin%S = multinomalspec(1)%S 00692 else 00693 specin%Hm0 = sqrt(sum(Hm0**2)) ! total wave height 00694 ind = maxval(maxloc(Hm0)) ! spectrum that is largest 00695 specin%fp = fp(ind) ! not really used in further calculation 00696 specin%dir0 = mainang(ind) ! again not really used in further calculation 00697 ! if all scoeff>1000 then all waves should be in the same direction exactly 00698 if (all(scoeff>=1024.d0) .and. all(mainang==mainang(ind))) then 00699 specin%scoeff = 1024.d0 00700 else 00701 specin%scoeff = min(scoeff(ind),999.d0) 00702 endif 00703 ! 00704 ! Now set total spectrum by summation of all subspectra, no checks 00705 allocate(specin%S(specin%nf,specin%nang)) 00706 specin%S = 0.d0 00707 do ip=1,nmodal 00708 specin%S = specin%S+multinomalspec(ip)%S 00709 enddo 00710 !! Now we have to loop over all partitioned spectra. Where two or more spectra 00711 !! overlap, only the largest is counted, and all others are set to zero. Afterwards 00712 !! all spectra are scaled so that the total energy is maintained. Since the scaling 00713 !! affects where spectra overlap, this loop is repeated until only minor changes 00714 !! in the scaling occur. Warnings or errors are given if the spectrum is scaled too 00715 !! much, i.e. spectra overlap too much. 00716 !! 00717 !! Allocate space 00718 !allocate(scaledspec(nmodal)) ! used to store scaled density spectra 00719 !do ip=1,nmodal 00720 ! allocate(scaledspec(ip)%S(specin%nf,specin%nang)) 00721 ! scaledspec(ip)%S = multinomalspec(ip)%S 00722 !enddo 00723 !allocate(scalefac1(nmodal)) ! this is the scaling factor required to maintain Hm0 00724 ! ! including that parts of the spectrum are set to zero 00725 !scalefac1 = 1.d0 00726 !allocate(scalefac2(nmodal)) ! this is the scaling factor required to maintain Hm0 00727 ! ! in the previous iteration 00728 !scalefac2 = 1.d0 00729 !allocate(avgscale(nmodal)) 00730 !avgscale = 1.d0 00731 !! these are convergence criteria 00732 !newconv = 0.d0 00733 !oldconv = huge(0.d0) 00734 !cont = .true. 00735 !allocate(tempmax(nmodal)) ! used to store maximum value in f,theta space 00736 !allocate(oldvariance(nmodal)) ! used to store the sum of variance in the original partitions 00737 !do ip=1,nmodal 00738 ! oldvariance(ip) = sum(multinomalspec(ip)%S) 00739 !enddo 00740 !allocate(newvariance(nmodal)) ! used to store the sum of variance in the new partitions 00741 !! 00742 !! Start convergence loop 00743 !do while (cont) 00744 ! avgscale = (avgscale+scalefac1)/2 00745 ! ! First scale the spectra 00746 ! do ip=1,nmodal 00747 ! ! scale the spectrum using less than half the additional scale factor, this to 00748 ! ! ensure the method does not become unstable 00749 ! ! scaledspec(ip)%S = multinomalspec(ip)%S*(0.51d0+scalefac1(ip)*0.49d0) 00750 ! ! scaledspec(ip)%S = multinomalspec(ip)%S*(scalefac1(ip)+scalefac2(ip))/2 00751 ! scaledspec(ip)%S = multinomalspec(ip)%S*avgscale(ip) 00752 ! 00753 ! enddo 00754 ! do i=1,specin%nang 00755 ! do ii=1,specin%nf 00756 ! ! vector of variance densities at this point in f,theta space 00757 ! do ip=1,nmodal 00758 ! tempmax(ip) = scaledspec(ip)%S(ii,i) 00759 ! end do 00760 ! ! All spectra other than the one that is greatest in this point 00761 ! ! is set to zero. In case multiple spectra are equal largest, 00762 ! ! the first spectrum in the list is chosen (minval(maxloc)) 00763 ! ind = minval(maxloc(tempmax)) 00764 ! do ip=1,nmodal 00765 ! if (ip/=ind) then 00766 ! scaledspec(ip)%S(ii,i) = 0.d0 00767 ! endif 00768 ! enddo 00769 ! enddo 00770 ! end do 00771 ! ! 00772 ! ! Now rescale adjusted partition spectra so that they match the incident Hm0 00773 ! scalefac2 = scalefac1 ! keep previous results 00774 ! do ip=1,nmodal 00775 ! newvariance(ip) = sum(scaledspec(ip)%S) 00776 ! if (newvariance(ip)>0.01d0*oldvariance(ip)) then 00777 ! scalefac1(ip) = oldvariance(ip)/newvariance(ip) ! want to maximise to a factor 2 00778 ! ! else can generate rediculous results 00779 ! else 00780 ! scalefac1(ip) = 0.d0 ! completely remove this spectrum 00781 ! endif 00782 ! enddo 00783 ! ! 00784 ! ! check convergence criteria (if error is increasing, we have passed best fit) 00785 ! newconv = maxval(abs(scalefac2-scalefac1)/scalefac1) 00786 ! if (newconv<0.0001d0 .or. abs(newconv-oldconv)<0.0001d0) then 00787 ! cont = .false. 00788 ! endif 00789 ! oldconv = newconv 00790 !end do 00791 !! Ensure full energy conservation by scaling now according to full scalefac, not just the 00792 !! half used in the iteration loop 00793 !do ip=1,nmodal 00794 ! scaledspec(ip)%S = multinomalspec(ip)%S*scalefac1(ip) 00795 !enddo 00796 !! 00797 !! Now check the total scaling that has taken place across the spectrum, also accounting 00798 !! for the fact that some parts of the spectra are set to zero. This differs from the 00799 !! scalefac1 and scalefac2 vectors in the loop that only adjust the part of the spectrum 00800 !! that is not zero. 00801 !do ip=1,nmodal 00802 ! write(*,*)'Hm0m = ',4*sqrt(sum(scaledspec(ip)%S)*specin%dang*specin%df) 00803 !enddo 00804 !do ip=1,nmodal 00805 ! indvec = maxloc(scaledspec(ip)%S) 00806 ! scalefac1(ip) = scaledspec(ip)%S(indvec(1),indvec(2)) / & 00807 ! multinomalspec(ip)%S(indvec(1),indvec(2))-1.d0 00808 !enddo 00809 !! 00810 !! Warning and/or error criteria here if spectra overlap each other too much 00811 !do ip=1,nmodal 00812 ! if(scalefac1(ip)>0.5d0) then 00813 ! if (forcepartition==1) then 00814 ! call writelog('lsw','(a,f0.0,a,i0,a,a,a)', & 00815 ! 'Warning: ',scalefac1(ip)*100,'% of energy in spectrum partition ''',ip, & 00816 ! ''' in ', trim(readfile),' is overlapped by other partitions') 00817 ! call writelog('lsw','',' Check spectral partitioning in ',trim(readfile)) 00818 ! else 00819 ! call writelog('lswe','(a,f0.0,a,i0,a,a,a)', & 00820 ! 'Error: ',scalefac1(ip)*100,'% of energy in spectrum partition ''',ip, & 00821 ! ''' in ', trim(readfile),' is overlapped by other partitions') 00822 ! call writelog('lswe','','This spectrum overlaps too much with another spectrum partition.', & 00823 ! ' Check spectral partitioning in ',trim(readfile)) 00824 ! call writelog('lswe','','If the partitioning should be carried out regardles of energy loss, ', & 00825 ! 'set ''forcepartition = 1'' in ',trim(readfile)) 00826 ! call halt_program 00827 ! endif 00828 ! elseif (scalefac1(ip)>0.2d0 .and. scalefac1(ip)<=0.5d0) then 00829 ! call writelog('lsw','(a,f0.0,a,i0,a,a,a)', & 00830 ! 'Warning: ',scalefac1(ip)*100,'% of energy in spectrum partition ''',ip, & 00831 ! ''' in ', trim(readfile),' is overlapped by other partitions') 00832 ! call writelog('lsw','',' Check spectral partitioning in ',trim(readfile)) 00833 ! elseif (scalefac1(ip)<0.d0) then 00834 ! call writelog('lsw','(a,i0,a,a,a)','Warning: spectrum partition ''',ip,''' in ',trim(readfile), & 00835 ! ' has been removed') 00836 ! call writelog('lsw','','This spectrum is entirely overlapped by another spectrum partition.',& 00837 ! ' Check spectral partitioning in ',trim(readfile)) 00838 ! endif 00839 !enddo 00840 !! 00841 !! Now set total spectrum 00842 !allocate(specin%S(specin%nf,specin%nang)) 00843 !specin%S = 0.d0 00844 !do ip=1,nmodal 00845 ! specin%S = specin%S+scaledspec(ip)%S 00846 !enddo 00847 ! 00848 !deallocate(scaledspec) 00849 !deallocate(tempmax) 00850 !deallocate(scalefac1,scalefac2) 00851 !deallocate(newvariance,oldvariance) 00852 endif 00853 00854 ! We need frequency spectrum to ensure Sf remains correct between interpoloation 00855 ! routines 00856 allocate(specin%Sf(specin%nf)) 00857 specin%Sf=0.d0 00858 do i=1,specin%nf 00859 do ii=1,specin%nang 00860 if (ii==1) then 00861 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(2)-specin%ang(1)) 00862 elseif (ii==specin%nang) then 00863 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii)-specin%ang(ii-1)) 00864 else 00865 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii+1)-specin%ang(ii-1))/2 00866 endif 00867 enddo 00868 enddo 00869 00870 deallocate (Dd) 00871 deallocate(Hm0,fp,gam,mainang,scoeff,tma) 00872 deallocate (x,y) 00873 deallocate(tempdir) 00874 ! TODO dereference pointers first.... 00875 do ip=1,nmodal 00876 deallocate(multinomalspec(ip)%S) 00877 end do 00878 deallocate(multinomalspec) 00879 00880 00881 end subroutine read_jonswap_file 00882 00883 00884 00885 ! ----------------------------------------------------------- 00886 ! --------- JONSWAP unscaled JONSWAP spectrum -------------- 00887 ! -------------(used by read_jonswap_files)------------------ 00888 subroutine jonswapgk(x,gam,y) 00889 00890 IMPLICIT NONE 00891 ! Required input: - x : nondimensional frequency, divided by the peak frequency 00892 ! - gam : peak enhancement factor, optional parameter (DEFAULT 3.3) 00893 ! - y is output : nondimensional relative spectral density, equal to one at the peak 00894 00895 real*8, INTENT(IN) :: gam 00896 real*8,dimension(:), INTENT(IN) :: x 00897 real*8,dimension(:), INTENT(INOUT) :: y 00898 00899 ! Internal variables 00900 real*8,dimension(size(x)) :: xa, sigma, fac1, fac2, fac3, temp 00901 00902 xa=abs(x) 00903 00904 where (xa==0) 00905 xa=1e-20 00906 end where 00907 00908 sigma=xa 00909 00910 where (sigma<1.) 00911 sigma=0.07 00912 end where 00913 00914 where (sigma>=1.) 00915 sigma=0.09 00916 end where 00917 00918 temp=0*xa+1 00919 00920 fac1=xa**(-5) 00921 fac2=exp(-1.25*(xa**(-4))) 00922 fac3=(gam*temp)**(exp(-((xa-1)**2)/(2*(sigma**2)))) 00923 00924 y=fac1*fac2*fac3 00925 y=y/maxval(y) 00926 00927 return 00928 00929 end subroutine jonswapgk 00930 00931 00932 ! -------------------------------------------------------------- 00933 ! ----------------------Read SWAN files ------------------------ 00934 ! -------------------------------------------------------------- 00935 subroutine read_swan_file(par,readfile,specin) 00936 00937 use params, only: parameters 00938 use logging_module, only: writelog, report_file_read_error 00939 use filefunctions, only: create_new_fid 00940 use xmpi_module, only: Halt_Program 00941 use math_tools, only: flipv,flipa 00942 00943 IMPLICIT NONE 00944 00945 ! Input / output variables 00946 type(parameters), INTENT(IN) :: par 00947 character(len=*), intent(IN) :: readfile 00948 type(spectrum),intent(inout) :: specin 00949 00950 ! Internal variables 00951 character(6) :: rtext 00952 real*8 :: factor,exc 00953 integer :: fid,switch 00954 integer :: i,ii,ier,ier2,ier3 00955 logical :: flipped 00956 00957 flipped =.false. 00958 switch = 0 00959 00960 call writelog('sl','','Reading from SWAN file: ',trim(readfile),' ...') 00961 fid = create_new_fid() 00962 open(fid,file=readfile,form='formatted',status='old') 00963 00964 ! Read file until RFREQ or AFREQ is found 00965 do while (switch==0) 00966 read(fid,'(a)',iostat=ier)rtext 00967 if (ier .ne. 0) then 00968 call report_file_read_error(readfile) 00969 endif 00970 if (rtext == 'RFREQ ') then 00971 switch = 1 00972 elseif (rtext == 'AFREQ ') then 00973 switch = 2 00974 end if 00975 end do 00976 00977 ! Read nfreq and f 00978 ! Note f is not monotonically increasing in most simulations 00979 read(fid,*,iostat=ier)specin%nf 00980 if (ier .ne. 0) then 00981 call report_file_read_error(readfile) 00982 endif 00983 allocate(specin%f(specin%nf)) 00984 do i=1,specin%nf 00985 read(fid,*,iostat=ier)specin%f(i) 00986 if (ier .ne. 0) then 00987 call report_file_read_error(readfile) 00988 endif 00989 end do 00990 00991 ! Convert to absolute frequencies: 00992 ! STILL TO BE DONE 00993 if (switch == 1) then 00994 specin%f = specin%f 00995 else 00996 specin%f = specin%f 00997 end if 00998 00999 ! Read CDIR or NDIR 01000 read(fid,'(a)',iostat=ier)rtext 01001 if (ier .ne. 0) then 01002 call report_file_read_error(readfile) 01003 endif 01004 if (rtext == 'NDIR ') then 01005 switch = 1 01006 elseif (rtext == 'CDIR ') then 01007 switch = 2 01008 else 01009 call writelog('ewls','', 'SWAN directional bins keyword not found') 01010 call halt_program 01011 endif 01012 01013 ! Read ndir, theta 01014 read(fid,*,iostat=ier)specin%nang 01015 if (ier .ne. 0) then 01016 call report_file_read_error(readfile) 01017 endif 01018 allocate(specin%ang(specin%nang)) 01019 do i=1,specin%nang 01020 read(fid,*,iostat=ier)specin%ang(i) 01021 if (ier .ne. 0) then 01022 call report_file_read_error(readfile) 01023 endif 01024 end do 01025 01026 ! Convert angles to cartesian degrees relative to East 01027 if (switch == 1) then 01028 ! nautical to cartesian East 01029 specin%ang = 270.d0-specin%ang 01030 else 01031 ! cartesian to cartesian East 01032 specin%ang = specin%ang-par%dthetaS_XB 01033 ! dthetaS_XB is the (counter-clockwise) angle in the degrees to rotate from the x-axis in SWAN to the 01034 ! x-axis pointing East 01035 end if 01036 01037 ! Ensure angles are increasing instead of decreasing 01038 if (specin%ang(2)<specin%ang(1)) then 01039 call flipv(specin%ang,size(specin%ang)) 01040 flipped=.true. 01041 end if 01042 01043 ! convert to radians 01044 specin%ang=specin%ang*par%px/180 01045 specin%dang=specin%ang(2)-specin%ang(1) 01046 01047 ! Skip Quant, next line, read VaDens or EnDens 01048 read(fid,'(a)',iostat=ier)rtext 01049 read(fid,'(a)',iostat=ier2)rtext 01050 read(fid,'(a)',iostat=ier3)rtext 01051 if (ier+ier2+ier3 .ne. 0) then 01052 call report_file_read_error(readfile) 01053 endif 01054 if (rtext == 'VaDens') then 01055 switch = 1 01056 elseif (rtext == 'EnDens') then 01057 switch = 2 01058 else 01059 call writelog('slwe','', 'SWAN VaDens/EnDens keyword not found') 01060 call halt_program 01061 end if 01062 read(fid,'(a)',iostat=ier)rtext 01063 read(fid,*,iostat=ier2)exc 01064 if (ier+ier2 .ne. 0) then 01065 call report_file_read_error(readfile) 01066 endif 01067 01068 i=0 01069 ! Find FACTOR keyword 01070 do while (i==0) 01071 read(fid,'(a)',iostat=ier)rtext 01072 if (ier .ne. 0) then 01073 call report_file_read_error(readfile) 01074 endif 01075 if (rtext == 'FACTOR') then 01076 i=1 01077 elseif (rtext == 'ZERO ') then 01078 call writelog('lswe','','Zero energy density input for this point') 01079 call halt_program 01080 elseif (rtext == 'NODATA') then 01081 call writelog('lwse','','SWAN file has no data for this point') 01082 call halt_program 01083 end if 01084 end do 01085 read(fid,*,iostat=ier)factor 01086 if (ier .ne. 0) then 01087 call report_file_read_error(readfile) 01088 endif 01089 01090 ! Read 2D S array 01091 allocate(specin%S(specin%nf,specin%nang)) 01092 do i=1,specin%nf 01093 read(fid,*,iostat=ier)(specin%S(i,ii),ii=1,specin%nang) 01094 if (ier .ne. 0) then 01095 call report_file_read_error(readfile) 01096 endif 01097 end do 01098 01099 ! Finished reading file 01100 close(fid) 01101 01102 ! Replace exception value 01103 where (specin%S == exc) 01104 specin%S = 0.d0 01105 endwhere 01106 01107 ! If angles were decreasing, flip S_array as also dir is flipped 01108 if (flipped) then 01109 call flipa(specin%S,specin%nf,specin%nang,2) 01110 end if 01111 01112 ! multiply by SWAN output factor 01113 specin%S=specin%S*factor 01114 01115 ! Convert from energy density to variance density 01116 if (switch == 2) then 01117 specin%S=specin%S/(par%rho*par%g) 01118 end if 01119 01120 ! Convert to m2/Hz/rad 01121 specin%S=specin%S*180/par%px 01122 01123 ! We need a value for spreading. The assumption is that it is less than 1000 01124 ! This way, wp%fgen will not be set to just one angle. 01125 specin%scoeff = -1.d0 01126 ! We need to know if hm0 was set explicitly, not the case for Swan files 01127 specin%hm0 = -1.d0 01128 01129 ! We need frequency spectrum to ensure Sf remains correct between interpoloation 01130 ! routines 01131 allocate(specin%Sf(specin%nf)) 01132 specin%Sf=0.d0 01133 do i=1,specin%nf 01134 do ii=1,specin%nang 01135 if (ii==1) then 01136 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(2)-specin%ang(1)) 01137 elseif (ii==specin%nang) then 01138 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii)-specin%ang(ii-1)) 01139 else 01140 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii+1)-specin%ang(ii-1))/2 01141 endif 01142 enddo 01143 enddo 01144 01145 end subroutine read_swan_file 01146 01147 01148 ! -------------------------------------------------------------- 01149 ! -----------------Read variance density files ----------------- 01150 ! -------------------------------------------------------------- 01151 subroutine read_vardens_file(par,readfile,specin) 01152 01153 use params, only: parameters 01154 use logging_module, only: writelog, report_file_read_error 01155 use filefunctions, only: create_new_fid 01156 use xmpi_module, only: Halt_Program 01157 use math_tools, only: flipv,flipa 01158 01159 IMPLICIT NONE 01160 01161 ! Input / output variables 01162 type(parameters), INTENT(IN) :: par 01163 character(len=*), intent(IN) :: readfile 01164 type(spectrum),intent(inout) :: specin 01165 01166 ! Internal variables 01167 integer :: fid,i,ii,nnz,ier 01168 01169 ! Open file to start read 01170 call writelog('sl','','Reading from vardens file: ',trim(readfile),' ...') 01171 fid = create_new_fid() 01172 open(fid,file=readfile,form='formatted',status='old') 01173 01174 ! Read number of frequencies and frequency vector 01175 read(fid,*,iostat=ier)specin%nf 01176 if (ier .ne. 0) then 01177 call report_file_read_error(readfile) 01178 endif 01179 allocate(specin%f(specin%nf)) 01180 do i=1,specin%nf 01181 read(fid,*,iostat=ier)specin%f(i) 01182 if (ier .ne. 0) then 01183 call report_file_read_error(readfile) 01184 endif 01185 end do 01186 01187 ! Read number of angles and angles vector 01188 read(fid,*,iostat=ier)specin%nang 01189 if (ier .ne. 0) then 01190 call report_file_read_error(readfile) 01191 endif 01192 allocate(specin%ang(specin%nang)) 01193 do i=1,specin%nang 01194 read(fid,*,iostat=ier)specin%ang(i) 01195 if (ier .ne. 0) then 01196 call report_file_read_error(readfile) 01197 endif 01198 end do 01199 01200 ! Convert from degrees to rad 01201 specin%ang=specin%ang*par%px/180 01202 specin%dang=specin%ang(2)-specin%ang(1) 01203 01204 ! Read 2D S array 01205 allocate(specin%S(specin%nf,specin%nang)) 01206 do i=1,specin%nf 01207 read(fid,*,iostat=ier)(specin%S(i,ii),ii=1,specin%nang) 01208 if (ier .ne. 0) then 01209 call report_file_read_error(readfile) 01210 endif 01211 end do 01212 01213 ! Finished reading file 01214 close(fid) 01215 01216 ! Convert to m2/Hz/rad 01217 specin%S=specin%S*180/par%px 01218 01219 ! We need a value for spreading. The assumption is that it is less than 1000 01220 ! if more than one direction column has variance. If only one column has variance 01221 ! we assume the model requires long crested waves, and wp%thetagen should be exactly 01222 ! equal to the input direction. We then also need to find the dominant angle. 01223 ! First flatten spectrum to direction only, then determine if all but one are zero, 01224 ! use the direction of the non-zero as the main wave direction, set spreading to 01225 ! a high value (greater than 1000). Else set spreading to a negative value. 01226 allocate(specin%Sd(specin%nang)) 01227 specin%Sd = sum(specin%S,DIM = 1) 01228 nnz = count(specin%Sd>0.d0) 01229 if (nnz == 1) then 01230 specin%scoeff = 1024.d0 01231 specin%dir0 = specin%ang(minval(maxloc(specin%Sd))) 01232 else 01233 specin%scoeff = -1.d0 01234 endif 01235 01236 ! We need frequency spectrum to ensure Sf remains correct between interpoloation 01237 ! routines 01238 allocate(specin%Sf(specin%nf)) 01239 specin%Sf=0.d0 01240 do i=1,specin%nf 01241 do ii=1,specin%nang 01242 if (ii==1) then 01243 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(2)-specin%ang(1)) 01244 elseif (ii==specin%nang) then 01245 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii)-specin%ang(ii-1)) 01246 else 01247 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii+1)-specin%ang(ii-1))/2 01248 endif 01249 enddo 01250 enddo 01251 01252 01253 ! We need to know if hm0 was set explicitly, not the case for vardens files 01254 specin%hm0 = -1.d0 01255 01256 end subroutine read_vardens_file 01257 01258 01259 ! -------------------------------------------------------------- 01260 ! ------------- Interpolate to standard spectrum --------------- 01261 ! -------------------------------------------------------------- 01262 subroutine interpolate_spectrum(specin,specinterp,par,fmax) 01263 01264 use interp, only: linear_interp, interp_in_cyclic_function, interp_using_trapez_rule 01265 use params, only: parameters 01266 01267 IMPLICIT NONE 01268 01269 ! input/output 01270 type(spectrum),intent(in) :: specin 01271 type(spectrum),intent(inout) :: specinterp 01272 type(parameters),intent(in) :: par 01273 real*8, intent(in) :: fmax 01274 ! internal 01275 integer :: i,j,dummy 01276 real*8 :: m0,df,dang 01277 ! real*8,dimension(naint) :: Sd 01278 real*8 :: hm0pre,hm0post,Sfnow,factor,xcycle 01279 real*8, dimension(:,:),allocatable :: Stemp 01280 01281 ! allocate size of f,ang,Sf and S arrays in specinterp 01282 allocate(specinterp%f(nfint)) 01283 allocate(specinterp%ang(naint)) 01284 allocate(specinterp%Sf(nfint)) 01285 allocate(specinterp%Sd(naint)) 01286 allocate(specinterp%S(nfint,naint)) 01287 allocate(Stemp(specin%nf,naint)) 01288 01289 ! fill f and ang arrays 01290 specinterp%nf = nfint 01291 specinterp%df = fmax/(nfint-1) ! this usually gives a range of 0-1Hz, 01292 ! with 1/800Hz resolution 01293 do i=1,nfint 01294 specinterp%f(i)=(i-1)*specinterp%df 01295 enddo 01296 specinterp%nang = naint 01297 specinterp%dang = 2*par%px/(naint-1) ! this is exactly the same as in the JONSWAP construction 01298 do i=1,specinterp%nang 01299 specinterp%ang(i)=(i-1)*specinterp%dang 01300 enddo 01301 01302 ! If hm0 was set explicitly, then use that, else calculate hm0 01303 if (specin%hm0>0.d0) then 01304 hm0pre = specin%hm0 01305 else 01306 ! pre-interpolation hm0 value (can be on a non-monotonic f,ang grid) 01307 m0 = 0 01308 do j=1,specin%nang 01309 do i=1,specin%nf 01310 df = specin%f(max(2,i))-specin%f(max(1,i-1)) 01311 dang = specin%ang(max(2,j))-specin%ang(max(1,j-1)) 01312 m0 = m0 + specin%S(i,j)*df*dang 01313 enddo 01314 enddo 01315 hm0pre = 4*sqrt(m0) 01316 endif 01317 01318 ! interpolation (no extrapolation) of input 2D spectrum to standard 2D spectrum 01319 ! do j=1,naint 01320 ! do i=1,nfint 01321 ! call linear_interp_2d(specin%f,specin%nf,specin%ang,specin%nang,specin%S, & 01322 ! specinterp%f(i),specinterp%ang(j),specinterp%S(i,j),'interp',0.d0) 01323 ! enddo 01324 ! enddo 01325 xcycle=2.d0*par%px 01326 do i=1,specin%nf 01327 call interp_in_cyclic_function(specin%ang,specin%S(i,:),specin%nang,xcycle,specinterp%ang,naint,Stemp(i,:)) 01328 enddo 01329 do j=1,naint 01330 do i=1,nfint 01331 if (specinterp%f(i)>specin%f(specin%nf).or.specinterp%f(i)<specin%f(1) ) then 01332 specinterp%S(i,j)=0.d0 01333 else 01334 call linear_interp(specin%f,Stemp(:,j),specin%nf,specinterp%f(i),specinterp%S(i,j),dummy) 01335 endif 01336 end do 01337 enddo 01338 01339 deallocate(Stemp) 01340 ! hm0 post (is always on a monotonic f,ang grid 01341 m0 = sum(specinterp%S)*specinterp%df*specinterp%dang 01342 hm0post = 4*sqrt(m0) 01343 01344 ! correct the wave energy to ensure hm0post == hm0pre 01345 ! specinterp%S = (hm0pre/hm0post)**2*specinterp%S 01346 ! This is done later using Sf 01347 01348 ! calculate 1D spectrum, summed over directions 01349 specinterp%Sf = sum(specinterp%S, DIM = 2)*specinterp%dang 01350 01351 ! correct the wave energy by setting Sfpost == Sfpre 01352 do i=1,nfint 01353 if (specinterp%f(i)>=minval(specin%f) .and. specinterp%f(i)<=maxval(specin%f)) then 01354 call LINEAR_INTERP(specin%f,specin%Sf,specin%nf,specinterp%f(i),Sfnow,dummy) 01355 if (specinterp%Sf(i)>0.d0 .and. Sfnow>0.d0) then 01356 factor = Sfnow/specinterp%Sf(i) 01357 specinterp%Sf(i) = specinterp%Sf(i)*factor 01358 specinterp%S(i,:) = specinterp%S(i,:)*factor 01359 elseif (specinterp%Sf(i)==0.d0 .and. Sfnow>=0.d0) then 01360 specinterp%Sf(i) = Sfnow 01361 dummy = maxval(maxloc(specinterp%S(i,:))) 01362 specinterp%S(i,dummy) = Sfnow/specinterp%dang 01363 else 01364 specinterp%Sf(i) = 0.d0 01365 specinterp%S(i,:) = 0.d0 01366 endif 01367 endif 01368 enddo 01369 01370 ! calculate other wave statistics from interpolated spectrum 01371 m0 = sum(specinterp%Sf)*specinterp%df 01372 specinterp%hm0 = 4*sqrt(m0) 01373 specinterp%Sd = sum(specinterp%S, DIM = 1)*specinterp%df 01374 i = maxval(maxloc(specinterp%Sd)) 01375 specinterp%dir0 = 270.d0 - specinterp%ang(i)*180/par%px ! converted back into nautical degrees 01376 i = maxval(maxloc(specinterp%Sf)) 01377 specinterp%fp = specinterp%f(i) 01378 call tpDcalc(specinterp%Sf,specinterp%f,specinterp%trep,par%trepfac,par%Tm01switch) 01379 specinterp%dirm = 270.d0-180.d0/par%px*atan2( sum(sin(specinterp%ang)*specinterp%Sd)/sum(specinterp%Sd),& 01380 sum(cos(specinterp%ang)*specinterp%Sd)/sum(specinterp%Sd) ) 01381 01382 specinterp%dirm = mod(specinterp%dirm,360.d0) 01383 01384 end subroutine interpolate_spectrum 01385 01386 ! -------------------------------------------------------------- 01387 ! ----------- Small subroutine to determine if the ------------ 01388 ! ---------- global reuseall should be true or false ----------- 01389 subroutine set_reuseall 01390 01391 implicit none 01392 ! internal 01393 integer :: i 01394 01395 ! set default to reuse all boundary conditions 01396 reuseall = .true. 01397 ! find any spectrum that cannot be reused, then change global 01398 ! reuseall to false 01399 do i=1,nspectra 01400 if (.not. bcfiles(i)%reuse) then 01401 reuseall = .false. 01402 endif 01403 enddo 01404 01405 end subroutine set_reuseall 01406 01407 ! -------------------------------------------------------------- 01408 ! ----------- Small subroutine to set the filenames ------------ 01409 ! ----------- of the boundary condition output files ----------- 01410 subroutine set_bcfilenames(wp) 01411 01412 implicit none 01413 ! input/output 01414 type(waveparamsnew),intent(inout) :: wp 01415 ! internal 01416 integer :: i1,i2,i3,i4,i5 01417 01418 if (reuseall) then 01419 wp%Efilename = 'E_reuse.bcf' 01420 wp%Esfilename = 'Es_reuse.bcf' 01421 wp%qfilename = 'q_reuse.bcf' 01422 wp%nhfilename = 'nh_reuse.bcf' 01423 else 01424 i1=floor(real(bccount)/10000) 01425 i2=floor(real(bccount-i1*10000)/1000) 01426 i3=floor(real(bccount-i1*10000-i2*1000)/100) 01427 i4=floor(real(bccount-i1*10000-i2*1000-i3*100)/10) 01428 i5=bccount-i1*10000-i2*1000-i3*100-i4*10 01429 wp%Efilename='E_series'//char(48+i1)//char(48+i2)//char(48+i3)//char(48+i4)//char(48+i5)//'.bcf' 01430 wp%Esfilename='Es_series'//char(48+i1)//char(48+i2)//char(48+i3)//char(48+i4)//char(48+i5)//'.bcf' 01431 wp%qfilename='q_series'//char(48+i1)//char(48+i2)//char(48+i3)//char(48+i4)//char(48+i5)//'.bcf' 01432 wp%nhfilename='nh_series'//char(48+i1)//char(48+i2)//char(48+i3)//char(48+i4)//char(48+i5)//'.bcf' 01433 endif 01434 01435 end subroutine set_bcfilenames 01436 01437 ! -------------------------------------------------------------- 01438 ! ----------- Merge all separate spectra into one -------------- 01439 ! -------------- average spectrum for other use ---------------- 01440 subroutine generate_combined_spectrum(specinterp,combspec,px,trepfac,Tm01switch) 01441 01442 implicit none 01443 type(spectrum),dimension(nspectra),intent(in) :: specinterp 01444 type(spectrum),intent(inout) :: combspec 01445 real*8, intent(in) :: px,trepfac 01446 integer,intent(in) :: Tm01switch 01447 integer :: iloc 01448 ! real*8,dimension(naint) :: Sd 01449 real*8,dimension(3) :: peakSd,peakang 01450 real*8,dimension(:),allocatable :: tempSf 01451 01452 allocate(combspec%f(nfint)) 01453 allocate(combspec%Sf(nfint)) 01454 allocate(combspec%Sd(naint)) 01455 allocate(combspec%ang(naint)) 01456 allocate(combspec%S(nfint,naint)) 01457 combspec%f=specinterp(1)%f 01458 combspec%nf=nfint 01459 combspec%df=specinterp(1)%df 01460 combspec%ang=specinterp(1)%ang 01461 combspec%nang=naint 01462 combspec%dang=specinterp(1)%dang 01463 combspec%trep = 0.d0 01464 combspec%S=0.d0 01465 combspec%Sf=0.d0 01466 01467 do iloc = 1,nspectra 01468 !combspec%trep = combspec%trep+specinterp(iloc)%trep/nspectra 01469 combspec%S = combspec%S +specinterp(iloc)%S /nspectra 01470 combspec%Sf = combspec%Sf +specinterp(iloc)%Sf /nspectra 01471 call tpDcalc(combspec%Sf,combspec%f,combspec%trep,trepfac,Tm01switch) 01472 enddo 01473 ! 01474 ! Calculate peak wave angle 01475 ! 01476 ! frequency integrated variance array 01477 combspec%Sd = sum(combspec%S,DIM=1)*combspec%df 01478 ! peak location of f-int array233333333 01479 iloc = maxval(maxloc(combspec%Sd)) 01480 ! pick two neighbouring directional bins, including effect of closing circle at 01481 ! 0 and 2pi rad 01482 if (iloc>1 .and. iloc<naint) then 01483 peakSd = (/combspec%Sd(iloc-1),combspec%Sd(iloc),combspec%Sd(iloc+1)/) 01484 peakang = (/combspec%ang(iloc-1),combspec%ang(iloc),combspec%ang(iloc+1)/) 01485 elseif (iloc==1) then 01486 peakSd = (/combspec%Sd(naint),combspec%Sd(1),combspec%Sd(2)/) 01487 peakang = (/combspec%ang(naint)-2*px,combspec%ang(1),combspec%ang(2)/) 01488 elseif (iloc==naint) then 01489 peakSd = (/combspec%Sd(naint-1),combspec%Sd(naint),combspec%Sd(1)/) 01490 peakang = (/combspec%ang(naint-1),combspec%ang(naint),combspec%ang(1)+2*px/) 01491 endif 01492 ! dir0 calculated as mean over peak and two neighbouring cells 01493 combspec%dir0 = sum(peakSd*peakang)/sum(peakSd) 01494 ! return to 0<=dir0<=2*px 01495 if (combspec%dir0>2*px) then 01496 combspec%dir0 = combspec%dir0-2*px 01497 elseif (combspec%dir0<0) then 01498 combspec%dir0 = combspec%dir0+2*px 01499 endif 01500 ! 01501 ! Compute peak wave frequency 01502 ! 01503 allocate(tempSf(size(combspec%Sf))) 01504 tempSf=0*combspec%Sf 01505 ! find frequency range of 95% wave energy around peak 01506 where (combspec%Sf>0.95d0*maxval(combspec%Sf)) 01507 tempSf=combspec%Sf 01508 end where 01509 ! smoothed peak is weighted mean frequency of 95% range 01510 combspec%fp = sum(tempSf*combspec%f)/sum(tempSf) 01511 end subroutine generate_combined_spectrum 01512 01513 01514 ! -------------------------------------------------------------- 01515 ! -------------- Calculate average water depth ----------------- 01516 ! ------------------ on offshore boundary ---------------------- 01517 pure function calculate_average_water_depth(s) result(h) 01518 01519 use spaceparamsdef, only: spacepars 01520 01521 implicit none 01522 01523 ! input 01524 type(spacepars),intent(in) :: s 01525 ! output 01526 real*8 :: h 01527 ! internal 01528 01529 h = sum((s%zs0(1,:)-s%zb(1,:))*s%dnc(1,:))/sum(s%dnc(1,:)) 01530 01531 end function calculate_average_water_depth 01532 01533 01534 ! -------------------------------------------------------------- 01535 ! ----------- Choose wave train components based on ------------ 01536 ! --------------------- combined spectrum ---------------------- 01537 subroutine generate_wavetrain_components(combspec,wp,par,x0,y0,xe,ye) 01538 01539 use logging_module, only: writelog 01540 use params, only: parameters 01541 use interp, only: linear_interp 01542 use wave_functions_module, only: iteratedispersion 01543 use math_tools, only: random 01544 01545 implicit none 01546 ! input/output 01547 type(spectrum),intent(in) :: combspec 01548 type(waveparamsnew),intent(inout) :: wp 01549 type(parameters),intent(in) :: par 01550 ! internal 01551 integer :: i,ii 01552 integer :: ind1,ind2,dummy,seed,ni 01553 real*8,dimension(:),allocatable :: randnums,pdflocal,cdflocal 01554 real*8 :: L0,L,kmax,fmax 01555 real*8, intent(in) :: x0,y0,xe,ye 01556 real*8 :: dx,dy,L_grid,alfa_grid,alfa_adj,alfa_wave,L_proj,L_wave,nr,L_proj_adj 01557 logical :: cont 01558 01559 01560 ! If we are running non-hydrostatic boundary conditions, we want to remove 01561 ! very high frequency components, as these will not be computed well anyway 01562 if (par%nonhspectrum==1) then 01563 kmax = wdmax/wp%h0 01564 fmax = par%g*kmax*tanh(kmax*wp%h0)/2/par%px 01565 else 01566 ! this is really already taken into account 01567 ! by specifying par%sprdthr 01568 fmax = 2*maxval(combspec%f) 01569 endif 01570 01571 ! Determine frequencies around peak frequency of one-dimensional 01572 ! non-directional variance density spectrum, based on factor sprdthr, which 01573 ! should be included in the determination of the wave boundary conditions 01574 call frange(par,combspec%f,combspec%Sf,fmax,ind1,ind2) 01575 01576 ! Calculate number of wave components to be included in determination of the 01577 ! wave boundary conditions based on the wave record length and width of the 01578 ! wave frequency range 01579 wp%K = ceiling(wp%rtbc*(combspec%f(ind2)-combspec%f(ind1))+1) 01580 ! also include minimum number of components 01581 wp%K = max(wp%K,Kmin) 01582 01583 01584 ! Allocate space in waveparams for all wave train components 01585 allocate(wp%fgen(wp%K)) 01586 allocate(wp%thetagen(wp%K)) 01587 allocate(wp%phigen(wp%K)) 01588 allocate(wp%kgen(wp%K)) 01589 allocate(wp%wgen(wp%K)) 01590 01591 ! Select equidistant wave components between the earlier selected range of 01592 ! frequencies around the peak frequency in the frange subfunction 01593 wp%dfgen = (combspec%f(ind2)-combspec%f(ind1))/(wp%K-1) 01594 do i=1,wp%K 01595 wp%fgen(i)=combspec%f(ind1)+(i-1)*wp%dfgen 01596 enddo 01597 01598 ! This subroutine needs to generate random phase / directions for the individual 01599 ! wave trains. Due to some strange Fortran properties, it is better to select 01600 ! one long vector with random numbers for both the phase and the direction, than 01601 ! one vector for each. 01602 ! Update random seed, if requested 01603 allocate(randnums(2*wp%K)) 01604 if (par%random==1) then 01605 CALL init_seed(seed) 01606 randnums(1) = random(seed) 01607 else 01608 randnums(1) = random(100) 01609 endif 01610 !call random_number(randnums) 01611 do i=2,2*wp%K 01612 randnums(i) = random(0) 01613 enddo 01614 01615 ! Determine a random phase for each wave train component between 0 and 2pi 01616 wp%phigen=randnums(1:wp%K)*2*par%px 01617 01618 ! Determine random directions for each wave train component, based on the CDF of 01619 ! the directional spectrum. For each wave train we will interpolate the directional 01620 ! distribution at that frequency, generate a CDF, and then interpolate the wave 01621 ! train direction from a random number draw and the CDF. 01622 allocate(pdflocal(naint)) 01623 allocate(cdflocal(naint)) 01624 do i=1,wp%K 01625 ! interpolate spectrum at this frequency per angle in the spectrum 01626 do ii=1,naint 01627 call LINEAR_INTERP(combspec%f,combspec%S(:,ii),combspec%nf, wp%fgen(i),pdflocal(ii),dummy) 01628 enddo 01629 ! convert to pdf by ensuring total integral == 1, assuming constant directional bin size 01630 pdflocal = pdflocal/sum(pdflocal) 01631 ! convert to cdf by trapezoidal integration, which is not biased for integration 01632 ! direction. 01633 ! Boundary condition, which may be nonzero: 01634 cdflocal(1) = pdflocal(1) 01635 do ii=2,naint 01636 cdflocal(ii) = cdflocal(ii-1) + (pdflocal(ii)+pdflocal(ii-1))/2 01637 ! Note: this only works if the directional 01638 ! bins are constant in size. Assumed multiplication 01639 ! by one. 01640 enddo 01641 ! interpolate random number draw across the cdf. Note that cdf(1) is not assumed to be zero and 01642 ! there is a posibility of drawing less than cdf(1). Therefore, if the random number is .lt. cdf(1), 01643 ! interpolation should take place across the back of the angle spectrum 01644 ! (i.e., from theta(end)-2pi : theta(1)). 01645 ! Note that we are using the second half of randnums now (K+1:end) 01646 if (randnums(wp%K+i)>=cdflocal(1)) then 01647 call LINEAR_INTERP(cdflocal,combspec%ang,naint,randnums(wp%K+i), wp%thetagen(i),dummy) 01648 else 01649 call LINEAR_INTERP((/0.d0,cdflocal(1)/),(/combspec%ang(naint)-2*par%px,combspec%ang(1)/), & 01650 2,randnums(wp%K+i),wp%thetagen(i),dummy) 01651 if (wp%thetagen(i)<0.d0) then 01652 wp%thetagen(i)=wp%thetagen(i)+2*par%px 01653 endif 01654 endif 01655 enddo 01656 01657 01658 ! determine wave number for each wave train component, using standard dispersion relation 01659 ! solver from wave_functions module. This function returns a negative wave length if the 01660 ! solver did not converge. 01661 do i=1,wp%K 01662 L0 = par%g*(1/wp%fgen(i))**2/2/par%px ! deep water wave length 01663 L = iteratedispersion(L0,L0,par%px,wp%h0) 01664 if (L<0.d0) then 01665 call writelog('lsw','','No dispersion convergence found for wave train ',i, & 01666 ' in boundary condition generation') 01667 L = -L 01668 endif 01669 wp%kgen(i) = 2*par%px/L 01670 enddo 01671 01672 ! Angular frequency 01673 wp%wgen = 2*par%px*wp%fgen 01674 01675 ! In cyclic models, default should adjust wave directions so that an integer number of waves fit in the domain 01676 if(par%cyclicdiradjust==1) then 01677 dx = xe-x0 01678 dy = ye-y0 01679 L_grid = sqrt(dx**2+dy**2) 01680 ! First check: L_grid greater than zero (not a fully circular model) 01681 if (L_grid>1.d-3) then 01682 alfa_grid = atan2(-dx,dy) ! note: not dy,dx 01683 do i=1,wp%K 01684 alfa_wave = wp%thetagen(i)-alfa_grid 01685 L_wave = 2*par%px/wp%kgen(i) 01686 L_proj = L_wave/sin(alfa_wave) 01687 nr = L_grid/L_proj 01688 ni = nint(nr) 01689 ni = min(abs(ni),floor(L_grid/L_wave))*sign(1,ni) 01690 ! cont = .true. 01691 ! do while (cont .and. ni .ne. 0) 01692 ! L_waveadj = L_grid/ni 01693 ! if (abs(L_waveadj)>L_wave) then 01694 ! ni = ni-1*sign(1.d0,nr) 01695 ! else 01696 ! cont = .false. 01697 ! endif 01698 ! enddo 01699 if (ni .ne. 0) then 01700 L_proj_adj = L_grid/ni 01701 alfa_adj = asin(L_wave/L_proj_adj)+alfa_grid 01702 else 01703 alfa_adj = alfa_grid 01704 endif 01705 wp%thetagen(i) = alfa_adj 01706 enddo 01707 else 01708 ! do nothing 01709 endif 01710 endif 01711 01712 01713 ! Free memory 01714 deallocate(pdflocal,cdflocal,randnums) 01715 01716 end subroutine generate_wavetrain_components 01717 01718 ! ----------------------------------------------------------- 01719 ! ----------- Small subroutine to determine ----------------- 01720 ! ------------ representative wave period ------------------- 01721 subroutine tpDcalc(Sf,f,Trep,trepfac,switch) 01722 01723 implicit none 01724 01725 real*8, dimension(:), intent(in) :: Sf, f 01726 real*8, intent(out) :: Trep 01727 real*8, intent(in) :: trepfac 01728 integer, intent(in) :: switch 01729 01730 real*8, dimension(:),allocatable :: temp 01731 01732 allocate(temp(size(Sf))) 01733 temp=0.d0 01734 where (Sf>=trepfac*maxval(Sf)) 01735 temp=1.d0 01736 end where 01737 01738 if (switch == 1) then 01739 Trep=sum(temp*Sf)/sum(temp*Sf*f) ! Tm01 01740 else 01741 Trep = sum(temp*Sf/max(f,0.001d0))/sum(temp*Sf) ! Tm-1,0 01742 endif 01743 01744 deallocate(temp) 01745 01746 end subroutine tpDcalc 01747 01748 01749 ! ----------------------------------------------------------- 01750 ! ---- Small subroutine to determine f-range round peak ----- 01751 ! ----------------------------------------------------------- 01752 subroutine frange(par,f,Sf,fmax,firstp,lastp,findlineout) 01753 01754 use params, only: parameters 01755 01756 implicit none 01757 01758 type(parameters),intent(in) :: par 01759 real*8, dimension(:), intent(in) :: f,Sf 01760 real*8,intent(in) :: fmax 01761 01762 integer, intent(out) :: firstp, lastp 01763 01764 real*8, dimension(:), intent(out),allocatable,optional :: findlineout 01765 real*8, dimension(:),allocatable :: temp, findline 01766 integer :: i = 0 01767 01768 ! find frequency range around peak 01769 allocate(findline(size(Sf))) 01770 findline=0*Sf 01771 where (Sf>par%sprdthr*maxval(Sf)) 01772 findline=1 01773 end where 01774 01775 ! find frequency range below fmax 01776 where(f>fmax) 01777 findline = 0 01778 endwhere 01779 01780 firstp=maxval(maxloc(findline)) ! Picks the first "1" in temp 01781 01782 allocate (temp(size(findline))) 01783 temp=(/(i,i=1,size(findline))/) 01784 lastp=maxval(maxloc(temp*findline)) ! Picks the last "1" in temp 01785 01786 if (present(findlineout)) then 01787 allocate(findlineout(size(Sf))) 01788 findlineout=findline 01789 endif 01790 deallocate(temp, findline) 01791 01792 end subroutine frange 01793 01794 ! ----------------------------------------------------------- 01795 ! --- Small subroutine to reseed random number generator ---- 01796 ! ----------------------------------------------------------- 01797 subroutine init_seed(outseed) 01798 INTEGER :: i, n, clock 01799 INTEGER, DIMENSION(:), ALLOCATABLE :: seed 01800 integer,intent(out) :: outseed 01801 ! RANDOM_SEED does not result in a random seed for each compiler, so we better make sure that we have a pseudo random seed. 01802 ! I am not sure what is a good n 01803 n = 40 01804 i = 1 01805 ! not sure what size means here 01806 ! first call with size 01807 CALL RANDOM_SEED(size = n) 01808 ! define a seed array of size n 01809 ALLOCATE(seed(n)) 01810 ! what time is it 01811 CALL SYSTEM_CLOCK(COUNT=clock) 01812 ! define the seed vector based on a prime, the clock and the set of integers 01813 seed = clock + 37 * (/ (i - 1, i = 1, n) /) 01814 ! if mpi do we need a different seed on each node or the same??? 01815 ! if we do need different seeds on each node 01816 ! seed *= some big prime * rank ? 01817 ! now use the seed 01818 CALL RANDOM_SEED(PUT = seed) 01819 outseed = seed(1) 01820 end subroutine init_seed 01821 01822 01823 ! -------------------------------------------------------------- 01824 ! ---------------- Setup time axis for waves ------------------- 01825 ! -------------------------------------------------------------- 01826 subroutine generate_wave_time_axis(par,wp) 01827 01828 use params, only: parameters 01829 01830 implicit none 01831 ! input/output 01832 type(waveparamsnew),intent(inout) :: wp 01833 type(parameters),intent(in) :: par 01834 ! internal 01835 integer :: ntaper, indend 01836 integer :: i 01837 01838 01839 01840 ! First assume that internal and bc-writing time step is the same 01841 wp%dtin = wp%dtbc 01842 wp%dtchanged = .false. 01843 wp%tslenbc = nint(wp%rtbc/wp%dtbc)+1 01844 01845 ! Check whether the internal frequency is high enough to describe the highest frequency 01846 ! wave train returned from frange (which can be used in the boundary conditions) 01847 if (par%nonhspectrum==0) then 01848 if (wp%dtin>0.5d0/wp%fgen(wp%K)) then 01849 wp%dtin = 0.5d0/wp%fgen(wp%K) 01850 wp%dtchanged = .true. 01851 endif 01852 else 01853 if (wp%dtin>0.1d0/wp%fgen(wp%K)) then 01854 wp%dtin = 0.1d0/wp%fgen(wp%K) 01855 wp%dtchanged = .true. 01856 endif 01857 endif 01858 01859 ! The length of the internal time axis should be even (for Fourier transform) and 01860 ! depends on the internal time step needed and the internal duration (~1/dfgen): 01861 wp%tslen = ceiling(1/wp%dfgen/wp%dtin)+1 01862 if (mod(wp%tslen,2)/=0) then 01863 wp%tslen = wp%tslen +1 01864 end if 01865 01866 ! Now we can make the internal time axis 01867 wp%rtin = wp%tslen * wp%dtin 01868 allocate(wp%tin(wp%tslen)) 01869 do i=1,wp%tslen 01870 wp%tin(i) = (i-1)*wp%dtin 01871 enddo 01872 01873 ! Make a taper function to slowly increase and decrease the boundary condition forcing 01874 ! at the start and the end of the boundary condition file (including any time beyond 01875 ! the external rtbc 01876 allocate(wp%taperf(wp%tslen)) 01877 allocate(wp%taperw(wp%tslen)) 01878 ! fill majority with unity 01879 wp%taperf = 1.d0 01880 wp%taperw = 1.d0 01881 if (par%nonhspectrum==1) then 01882 ! begin taper by building up the wave conditions over 2 wave periods 01883 ntaper = nint((2.0d0*par%Trep)/wp%dtin) 01884 else 01885 ntaper = nint((5.d0*par%Trep)/wp%dtin) 01886 endif 01887 do i=1,min(ntaper,size(wp%taperf)) 01888 wp%taperf(i) = tanh(5.d0*i/ntaper) ! multiplied by five because tanh(5)=~1 01889 wp%taperw(i) = tanh(5.d0*i/ntaper) 01890 enddo 01891 ! We do not want to taperw the end anymore. Instead we pass the wave height at the end of rtbc to 01892 ! the next wave generation iteration. 01893 ! end taper by finding where tin=rtbc, taper before that and set everything to zero after 01894 ! that. 01895 if (wp%tin(wp%tslen)>wp%rtbc) then 01896 indend = minval(minloc(wp%tin,MASK = wp%tin>=wp%rtbc)) 01897 else 01898 indend = wp%tslen 01899 endif 01900 do i=1,min(ntaper,indend) 01901 wp%taperf(indend+1-i) = min(wp%taperf(indend+1-i),tanh(5.d0*i/ntaper)) 01902 enddo 01903 wp%taperf(indend:wp%tslen) = 0.d0 01904 ind_end_taper = indend 01905 end subroutine generate_wave_time_axis 01906 01907 01908 ! -------------------------------------------------------------- 01909 ! -------- Calculate variance at each spectrum location -------- 01910 ! -------------------------------------------------------------- 01911 subroutine generate_wave_train_variance(wp,specinterp) 01912 01913 use interp, only: linear_interp 01914 01915 01916 implicit none 01917 ! input/output 01918 type(waveparamsnew),intent(inout) :: wp 01919 type(spectrum),dimension(nspectra),intent(in):: specinterp 01920 ! internal 01921 integer :: i,ii, dummy 01922 real*8 :: hm0post 01923 01924 ! allocate space for the variance arrays 01925 allocate(wp%vargen(nspectra)) 01926 allocate(wp%vargenq(nspectra)) 01927 01928 ! Determine variance at each spectrum location 01929 do i=1,nspectra 01930 allocate(wp%vargen(i)%Sf(wp%K)) 01931 allocate(wp%vargenq(i)%Sf(wp%K)) 01932 do ii=1,wp%K 01933 ! In order to maintain energy density per frequency, an interpolation 01934 ! is carried out on the input directional variance density spectrum 01935 ! at the frequency and direction locations in fgen. This must be done 01936 ! over the 2D spectrum, not 1D spectrum. 01937 ! call linear_interp_2d(specinterp(i)%f,nfint, & ! input frequency, size 01938 ! specinterp(i)%ang,naint, & ! input angles, size 01939 ! specinterp(i)%S*specinterp(i)%dang, & ! input variance density times grid angle resolution 01940 ! wp%fgen(ii),wp%thetagen(ii), & ! output frquency/angle 01941 ! wp%vargen(i)%Sf(ii), & ! output variance density 01942 ! 'interp',0.d0) ! method and exception return value 01943 call linear_interp(specinterp(i)%f, & 01944 specinterp(i)%Sf, & 01945 nfint, & 01946 wp%fgen(ii), & 01947 wp%vargen(i)%Sf(ii),dummy) 01948 enddo 01949 ! Correct variance to ensure the total variance remains the same as the total variance in the 01950 ! (interpolated) input spectrum 01951 hm0post = 4*sqrt(sum(wp%vargen(i)%Sf)*wp%dfgen) 01952 wp%vargen(i)%Sf = (specinterp(i)%hm0/hm0post)**2*wp%vargen(i)%Sf 01953 ! For the generation of long waves we cannot use wp%vargen%Sf, because it contains an overestimation 01954 ! of energy in the peak frequencies. We can also not use the standard directionally-integrated spectrum 01955 ! because this stores all energy at fgen(K), where it is possible that for this current spectrum, there 01956 ! is no energy at S(f(K),theta(K0)) 01957 ! The current solution is to take the minimum of both methods 01958 do ii=1,wp%K 01959 ! Map Sf input to fgen 01960 call LINEAR_INTERP(specinterp(i)%f,specinterp(i)%Sf,nfint,wp%fgen(ii),wp%vargenq(i)%Sf(ii),dummy) 01961 enddo 01962 wp%vargenq(i)%Sf = min(wp%vargen(i)%Sf,wp%vargenq(i)%Sf) 01963 enddo 01964 01965 end subroutine generate_wave_train_variance 01966 01967 ! -------------------------------------------------------------- 01968 ! ------ Calculate amplitudes at each spectrum location -------- 01969 ! ------------ for every wave train component ------------------ 01970 subroutine generate_wave_train_properties_per_offshore_point(wp,s,nonhspectrum,swkhmin) 01971 01972 use spaceparamsdef, only: spacepars 01973 use interp, only: linear_interp 01974 01975 implicit none 01976 ! input/output 01977 type(waveparamsnew),intent(inout) :: wp 01978 type(spacepars),intent(in) :: s 01979 integer,intent(in) :: nonhspectrum 01980 real*8,intent(in) :: swkhmin 01981 ! internal 01982 integer :: i,ii,dummy 01983 integer,dimension(2) :: interpindex 01984 real*8,dimension(2) :: positions,sf 01985 logical :: interpcase 01986 real*8 :: here,sfnow,sfimp 01987 integer,dimension(size(n_index_loc)) :: temp_index_loc 01988 01989 interpindex = -123 01990 ! allocate space for the amplitude array and representative integration angle 01991 allocate(wp%A(s%ny+1,wp%K)) 01992 ! allocate(wp%danggen(s%ny+1)) 01993 allocate(wp%Sfinterp(s%ny+1,wp%K)) 01994 allocate(wp%Sfinterpq(s%ny+1,wp%K)) 01995 allocate(wp%Hm0interp(s%ny+1)) 01996 01997 ! where necessary, interpolate Sf of each spectrum location 01998 ! to the current grid cell, and use this to calculate A 01999 do i=1,s%ny+1 02000 ! Four possibilities: 02001 ! 1: this exact point has an input spectrum 02002 ! 2: this point lies between two input spectra 02003 ! 3: this point lies before any input spectrum 02004 ! 4: this point lies beyond any input spectrum 02005 if (any(n_index_loc==i)) then ! case 1 02006 interpcase = .false. 02007 do ii=1,nspectra 02008 if(n_index_loc(ii)==i) then 02009 interpindex = ii 02010 endif 02011 enddo 02012 elseif(i<minval(n_index_loc)) then ! case 3 02013 interpcase = .false. 02014 do ii=1,nspectra 02015 if(n_index_loc(ii)==minval(n_index_loc)) then 02016 interpindex = ii 02017 endif 02018 enddo 02019 elseif(i>maxval(n_index_loc)) then ! case 4 02020 interpcase = .false. 02021 do ii=1,nspectra 02022 if(n_index_loc(ii)==maxval(n_index_loc)) then 02023 interpindex = ii 02024 endif 02025 enddo 02026 else ! case 2 02027 interpcase = .true. 02028 temp_index_loc = n_index_loc 02029 where(n_index_loc>i) 02030 temp_index_loc = -huge(0) 02031 endwhere 02032 interpindex(1) = minval(maxloc(temp_index_loc)) 02033 02034 temp_index_loc = n_index_loc 02035 where(n_index_loc<i) 02036 temp_index_loc = huge(0) 02037 endwhere 02038 interpindex(2) = minval(minloc(temp_index_loc)) 02039 endif 02040 02041 ! interpolate Sf 02042 if (interpcase) then 02043 positions(1) = s%ndist(1,n_index_loc(interpindex(1))) 02044 positions(2) = s%ndist(1,n_index_loc(interpindex(2))) 02045 here = s%ndist(1,i) 02046 do ii=1,wp%K 02047 sf(1) = wp%vargen(interpindex(1))%Sf(ii) 02048 sf(2) = wp%vargen(interpindex(2))%Sf(ii) 02049 call LINEAR_INTERP(positions,sf,2,here,wp%Sfinterp(i,ii),dummy) 02050 sf(1) = wp%vargenq(interpindex(1))%Sf(ii) 02051 sf(2) = wp%vargenq(interpindex(2))%Sf(ii) 02052 call LINEAR_INTERP(positions,sf,2,here,wp%Sfinterpq(i,ii),dummy) 02053 enddo 02054 ! Ensure that the total amount of variance has not been reduced/increased 02055 ! during interpolation for short wave energy only. This is not done for 02056 ! Sf for bound long wave generation 02057 sfnow = sum(wp%Sfinterp(i,:)) ! integration over fixed bin df size unnecessary 02058 sf(1) = sum(wp%vargen(interpindex(1))%Sf) 02059 sf(2) = sum(wp%vargen(interpindex(2))%Sf) 02060 call LINEAR_INTERP(positions,sf,2,here,sfimp,dummy) 02061 wp%Sfinterp(i,:)=wp%Sfinterp(i,:)*sfimp/sfnow 02062 else 02063 wp%Sfinterp(i,:) = wp%vargen(interpindex(1))%Sf 02064 wp%Sfinterpq(i,:) = wp%vargenq(interpindex(1))%Sf 02065 endif 02066 02067 wp%A(i,:) = sqrt(2*wp%Sfinterp(i,:)*wp%dfgen) 02068 02069 wp%Hm0interp(i) = 4*sqrt(sum(wp%Sfinterp(i,:))*wp%dfgen) 02070 02071 ! determine if these components will be pahse-resolved, or resolved in the wave energy balance 02072 allocate(wp%PRindex(wp%K)) 02073 if (nonhspectrum==1) then 02074 ! all components phase-resolved 02075 wp%PRindex = 1 02076 else 02077 if(swkhmin>0.d0) then 02078 ! some components resolved, some not 02079 where(wp%kgen*wp%h0 >= swkhmin) 02080 wp%PRindex = 0 02081 elsewhere 02082 wp%PRindex = 1 02083 endwhere 02084 else 02085 ! no components phase-resolved 02086 wp%PRindex = 0 02087 endif 02088 endif 02089 02090 02091 enddo ! i=1,s%ny+1 02092 02093 end subroutine generate_wave_train_properties_per_offshore_point 02094 02095 ! -------------------------------------------------------------- 02096 ! --------- Calculate Fourier componets for each wave ---------- 02097 ! ---------- train component on the offshore boundary ---------- 02098 subroutine generate_wave_train_Fourier(wp,s,par) 02099 02100 use spaceparamsdef, only: spacepars 02101 use math_tools, only: fftkind, flipiv 02102 use params, only: parameters 02103 use logging_module, only: writelog, progress_indicator 02104 use logging_module 02105 02106 implicit none 02107 ! input/output 02108 type(waveparamsnew),intent(inout) :: wp 02109 type(spacepars),intent(in) :: s 02110 type(parameters),intent(in) :: par 02111 ! internal 02112 integer :: i,ii,tempi 02113 complex(fftkind),dimension(:),allocatable :: tempcmplx 02114 02115 ! Determine indices of wave train components in frequency axis and 02116 ! Fourier transform result 02117 allocate(wp%Findex(wp%K)) 02118 tempi = floor(wp%fgen(1)/wp%dfgen) 02119 do i=1,wp%K 02120 wp%Findex(i) = tempi + i 02121 enddo 02122 02123 ! Allocate Fourier coefficients for each y-position at the offshore boundary and 02124 ! each time step 02125 allocate(wp%CompFn(s%ny+1,wp%tslen)) 02126 wp%CompFn=0.d0 02127 allocate(tempcmplx(wp%tslen/2-1)) 02128 call writelog('ls','','Calculating Fourier components') 02129 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02130 do i=1,wp%K 02131 call progress_indicator(.false.,dble(i)/wp%K*100,5.d0,2.d0) 02132 do ii=1,s%ny+1 02133 ! Determine first half of complex Fourier coefficients of wave train 02134 ! components using random phase and amplitudes from sampled spectrum 02135 ! until Nyquist frequency. The amplitudes are represented in a 02136 ! two-sided spectrum, which results in the factor 1/2. 02137 ! Unroll Fourier components along offshore boundary, assuming all wave trains 02138 ! start at x(1,1), y(1,1). 02139 wp%CompFn(ii,wp%Findex(i)) = wp%A(ii,i)/2*exp(compi*wp%phigen(i))* & ! Bas: wp%Findex used in time dimension because t=j*dt in frequency space 02140 exp(-compi*wp%kgen(i)*(dsin(wp%thetagen(i))*(s%yz(1,ii)-s%yz(1,1)) & 02141 +dcos(wp%thetagen(i))*(s%xz(1,ii)-s%xz(1,1))) ) 02142 02143 ! Determine Fourier coefficients beyond Nyquist frequency (equal to 02144 ! coefficients at negative frequency axis) of relevant wave components for 02145 ! first y-coordinate by mirroring 02146 tempcmplx = conjg(wp%CompFn(ii,2:wp%tslen/2)) 02147 call flipiv(tempcmplx,size(tempcmplx)) 02148 wp%CompFn(ii,wp%tslen/2+2:wp%tslen)=tempcmplx 02149 enddo 02150 enddo 02151 02152 ! Free memory 02153 deallocate(tempcmplx) 02154 02155 end subroutine generate_wave_train_Fourier 02156 02157 ! -------------------------------------------------------------- 02158 ! --------- Calculate in which computational wave bin ---------- 02159 ! ------------- each wave train component belongs -------------- 02160 subroutine distribute_wave_train_directions(wp,s,px,nspr,nonhspec) 02161 02162 use spaceparams 02163 use logging_module 02164 02165 implicit none 02166 02167 ! input/output 02168 type(waveparamsnew),intent(inout) :: wp 02169 type(spacepars),intent(in) :: s 02170 real*8,intent(in) :: px 02171 integer,intent(in) :: nspr 02172 logical,intent(in) :: nonhspec 02173 ! internal 02174 integer :: i,ii 02175 real*8,dimension(:),allocatable :: binedges 02176 logical :: toosmall,toolarge 02177 real*8 :: lostvar,keptvar,perclost 02178 integer :: lntheta ! local copies that can be changed without 02179 real*8 :: lthetamin,ldtheta ! damage to the rest of the model 02180 02181 ! Set basic parameters for comparison if using nonh spectrum 02182 if (nonhspec) then 02183 lntheta = 1 02184 lthetamin = -px 02185 ldtheta = 2*px 02186 else 02187 lntheta = s%ntheta 02188 lthetamin = s%thetamin 02189 ldtheta = s%dtheta 02190 endif 02191 02192 ! Calculate the bin edges of all the computational wave bins in the 02193 ! XBeach model (not the input spectrum) 02194 allocate(binedges(lntheta+1)) 02195 do i=1,lntheta+1 02196 binedges(i) = lthetamin+(i-1)*ldtheta 02197 enddo 02198 02199 ! Try to fit as many wave train components as possible between the 02200 ! minimum and maximum bin edges by adding/subtracting 2pi rad. This 02201 ! solves the problem of XBeach not including energy in the computation 02202 ! than should be in there. For instance XBeach does not realize that 02203 ! 90 degrees == -270 degrees == 450 degrees. 02204 ! Note: this does not ensure all energy is included in the wave bins, 02205 ! as wave energy may still fall outside the computational domain. 02206 do i=1,wp%K 02207 if (wp%thetagen(i)>maxval(binedges)) then 02208 toosmall = .false. 02209 toolarge = .true. 02210 elseif (wp%thetagen(i)<minval(binedges)) then 02211 toosmall = .true. 02212 toolarge = .false. 02213 else 02214 toosmall = .false. 02215 toolarge = .false. 02216 endif 02217 ! Continue to reduce/increase wave direction until 02218 ! condition is not met. 02219 do while (toosmall) 02220 wp%thetagen(i) = wp%thetagen(i) + 2*px 02221 toosmall = wp%thetagen(i) < minval(binedges) 02222 enddo 02223 do while (toolarge) 02224 wp%thetagen(i) = wp%thetagen(i) - 2*px 02225 toolarge = wp%thetagen(i) > maxval(binedges) 02226 enddo 02227 enddo 02228 02229 ! distribute the wave train directions among the wave directional bins, 02230 ! store the corresponding bin per wave train in WDindex. 02231 allocate(wp%WDindex(wp%K)) 02232 wp%WDindex = 0 02233 do i=1,wp%K 02234 if (wp%thetagen(i) < minval(binedges)) then 02235 ! remove from computation if wave direction still not within computational 02236 ! bounds 02237 wp%WDindex(i) = 0 02238 elseif (wp%thetagen(i) > maxval(binedges)) then 02239 wp%WDindex(i) = lntheta+1 02240 elseif (wp%thetagen(i) == binedges(lntheta+1) ) then 02241 ! catch upper boundary, and add to upper bin 02242 wp%WDindex(i) = lntheta 02243 else 02244 ! fit this wave train component into one of the lower bins 02245 do ii=1,lntheta 02246 if (wp%thetagen(i)>=binedges(ii) .and. wp%thetagen(i)<binedges(ii+1)) then 02247 wp%WDindex(i) = ii 02248 endif 02249 enddo 02250 endif 02251 enddo 02252 02253 ! If the user has set nspr == 1 then the randomly drawn wave directions 02254 ! should be set to the centres of the wave directional bins. 02255 ! Also move all wave energy falling outside the computational bins, into 02256 ! the computational domain (in the outer wave direction bins) 02257 if (nspr==1) then 02258 if (nonhspec) then 02259 ! We only really want to include all energy between -90 and +90 degrees. 02260 ! Waves outside this direction are sent out of the model 02261 where(wp%WDindex==1) 02262 wp%thetagen = 0.d0 02263 elsewhere 02264 wp%thetagen = -px 02265 endwhere 02266 else 02267 do i=1,wp%K 02268 if (wp%WDindex(i)==0) then 02269 wp%WDindex(i)=1 02270 elseif (wp%WDindex(i)>lntheta) then 02271 wp%WDindex(i)=lntheta 02272 endif 02273 ! reset the direction of this wave train to the centre of the bin 02274 wp%thetagen(i)=s%theta(wp%WDindex(i)) 02275 enddo 02276 endif 02277 endif 02278 02279 ! Check the amount of energy lost to wave trains falling outside the computational 02280 ! domain 02281 lostvar = 0.d0 02282 keptvar = 0.d0 02283 do i=1,wp%K 02284 if (wp%WDindex(i)==0) then 02285 lostvar = lostvar + sum(wp%A(:,i)**2) 02286 else 02287 keptvar = keptvar + sum(wp%A(:,i)**2) 02288 endif 02289 enddo 02290 perclost = 100*(lostvar/(lostvar+keptvar)) 02291 if (perclost>5.0d0) then 02292 call writelog('lsw','(a,f0.1,a)','Large amounts of energy (',perclost, & 02293 '%) fall outside computational domain at the offshore boundary') 02294 call writelog('lsw','','Check specification of input wave angles and wave directional grid') 02295 else 02296 call writelog('ls','(a,f0.1,a)','Wave energy outside computational domain at offshore boundary: ',perclost,'%') 02297 endif 02298 02299 ! Free memory 02300 deallocate(binedges) 02301 02302 end subroutine distribute_wave_train_directions 02303 02304 ! -------------------------------------------------------------- 02305 ! --------- Calculate energy envelope time series from --------- 02306 ! -------- Fourier components, and write to output file -------- 02307 subroutine generate_ebcf(wp,s,par) 02308 02309 use params, only: parameters 02310 use spaceparamsdef, only: spacepars 02311 use math_tools, only: fftkind, fft, flipiv, hilbert 02312 use interp, only: linear_interp 02313 use logging_module, only: writelog, progress_indicator 02314 use filefunctions, only: create_new_fid 02315 02316 implicit none 02317 02318 ! input/output 02319 type(waveparamsnew),intent(inout) :: wp 02320 type(spacepars),intent(in) :: s 02321 type(parameters),intent(in) :: par 02322 ! internal 02323 integer :: itheta,iy,iwc,it,irec 02324 integer :: index,status 02325 integer :: reclen,fid 02326 integer,dimension(:),allocatable :: nwc 02327 real*8,dimension(:,:,:), allocatable :: zeta, Ampzeta, E_tdir, E_interp 02328 real*8,dimension(:,:), allocatable :: eta, Amp 02329 complex(fftkind),dimension(:),allocatable :: Gn, tempcmplx,tempcmplxhalf 02330 integer,dimension(:),allocatable :: tempindex,tempinclude 02331 real*8 :: stdzeta,stdeta,etot,perc,emean 02332 real*8, dimension(:), allocatable :: E_t 02333 02334 if (par%bclwonly==0) then 02335 02336 ! Allocate variables for water level exitation and amplitude with and without 02337 ! directional spreading dependent envelope 02338 allocate(zeta(s%ny+1,wp%tslen,s%ntheta)) 02339 allocate(Ampzeta(s%ny+1,wp%tslen,s%ntheta)) 02340 zeta=0.d0 02341 Ampzeta=0.d0 02342 02343 allocate(eta(s%ny+1,wp%tslen)) 02344 allocate(Amp(s%ny+1,wp%tslen)) 02345 eta=0.d0 02346 Amp=0.d0 02347 02348 ! Calculate wave energy for each y-coordinate along seaside boundary for 02349 ! current computational directional bin 02350 allocate(nwc(s%ntheta)) 02351 allocate(Gn(wp%tslen)) 02352 allocate(tempinclude(wp%K)) 02353 allocate(tempcmplx(wp%tslen)) 02354 allocate(tempcmplxhalf(size(Gn(wp%tslen/2+2:wp%tslen)))) 02355 do itheta=1,s%ntheta 02356 call writelog('ls','(A,I0,A,I0)','Calculating short wave time series for theta bin ',itheta,' of ',s%ntheta) 02357 ! Select wave components that are in the current computational 02358 ! directional bin 02359 tempinclude=0 02360 where (wp%WDindex==itheta .and. wp%PRindex==0) ! only include components in this bin that are not to be phase-resolved 02361 tempinclude=1 02362 end where 02363 02364 ! Determine number of wave components that are in the current 02365 ! computational directional bin 02366 nwc(itheta)=sum(tempinclude) 02367 02368 ! Determine for each wave component in the current computational 02369 ! directional bin its index in the Fourier coefficients array 02370 ! ordered from high to low frequency 02371 allocate(tempindex(nwc(itheta))) 02372 tempindex=0 02373 do iwc=1,nwc(itheta) 02374 ! find highest 02375 index=maxval(maxloc(tempinclude)) 02376 ! reset that one so that the next highest in found in next iteration 02377 tempinclude(index)=0 02378 tempindex(iwc)=wp%Findex(index) 02379 end do 02380 02381 ! Check whether any wave components are in the current computational 02382 ! directional bin 02383 if (nwc(itheta)>0) then 02384 do iy=1,s%ny+1 02385 ! Reset 02386 Gn=0 02387 02388 ! Determine Fourier coefficients of all wave components for current 02389 ! y-coordinate in the current computational directional bin 02390 Gn(tempindex)=wp%CompFn(iy,tempindex) 02391 tempcmplxhalf = conjg(Gn(2:wp%tslen/2)) 02392 call flipiv(tempcmplxhalf,size(tempcmplxhalf)) 02393 Gn(wp%tslen/2+2:wp%tslen)=tempcmplxhalf 02394 02395 ! Inverse Discrete Fourier transformation to transform back to time 02396 ! domain from frequency domain 02397 tempcmplx=Gn 02398 status=0 02399 tempcmplx=fft(tempcmplx,inv=.true.,stat=status) 02400 02401 ! Scale result 02402 tempcmplx=tempcmplx/sqrt(dble(size(tempcmplx))) 02403 02404 ! Superimpose gradual increase and decrease of energy input for 02405 ! current y-coordinate and computational diretional bin on 02406 ! instantaneous water level excitation 02407 ! zeta(iy,:,itheta)=dble(tempcmplx*wp%tslen)*wp%taperw 02408 ! Robert: use final wave elevation from last iteration to startup 02409 ! this boundary condition 02410 zeta(iy,:,itheta)=dble(tempcmplx*wp%tslen) 02411 zeta(iy,:,itheta)=zeta(iy,:,itheta)*wp%taperw+(1-wp%taperw)*lastwaveelevation(iy,itheta) 02412 ! note that taperw is only zero at the start, not the end of the time series, so we 02413 ! can safely copy zeta at the point in time where t=rtbc to lastwaveelevation for use 02414 ! in the generation of the next time series 02415 lastwaveelevation(:,itheta) = zeta(iy,ind_end_taper,itheta) 02416 enddo ! iy=1,ny+1 02417 endif ! nwc>0 02418 deallocate(tempindex) 02419 enddo ! itheta = 1,ntheta 02420 deallocate(tempinclude) 02421 deallocate(Gn) 02422 deallocate(tempcmplxhalf) 02423 ! 02424 ! Calculate energy envelope amplitude 02425 call writelog('ls','(A)','Calculating wave energy envelope at boundary.') 02426 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02427 do iy=1,s%ny+1 02428 ! Integrate instantaneous water level excitation of wave components over directions 02429 eta(iy,:) = sum(zeta(iy,:,:),2) 02430 tempcmplx=eta(iy,:) 02431 ! 02432 ! Hilbert transformation to determine envelope of all total 02433 ! non-directional wave components 02434 call hilbert(tempcmplx,size(tempcmplx)) 02435 ! 02436 ! Determine amplitude of water level envelope by calculating 02437 ! the absolute value of the complex wave envelope descriptions 02438 Amp(iy,:)=abs(tempcmplx) 02439 ! 02440 ! Calculate standard deviation of non-directional 02441 ! instantaneous water level excitation of all 02442 ! wave components to be used as weighing factor 02443 stdeta = sqrt(sum(eta(iy,:)**2)/(size(eta(iy,:))-1)) 02444 do itheta=1,s%ntheta 02445 if (nwc(itheta)>0) then 02446 ! Calculate standard deviations of directional 02447 ! instantaneous water level excitation of all 02448 ! wave components to be used as weighing factor 02449 stdzeta = sqrt(sum(zeta(iy,:,itheta)**2)/(size(zeta(iy,:,itheta))-1)) 02450 ! 02451 ! Calculate amplitude of directional wave envelope 02452 Ampzeta(iy,:,itheta)= Amp(iy,:)*stdzeta/stdeta 02453 else ! nwc==0 02454 ! Current computational directional bin does not contain any wave 02455 ! components, so print message to screen 02456 Ampzeta(iy,:,itheta)=0.d0 02457 endif ! nwc>0 02458 end do ! 1:ntheta 02459 ! Print status message to screen 02460 !call writelog('ls','(A,I0,A,I0,A)','Y-point ',iy,' of ',s%ny+1,' done.') 02461 call progress_indicator(.false.,dble(iy)/(s%ny+1)*100,5.d0,2.d0) 02462 end do ! 1:ny+1 02463 ! 02464 ! free memory 02465 deallocate(tempcmplx) 02466 deallocate(nwc) 02467 ! 02468 else ! only long wave forcing at the boundary 02469 allocate(Ampzeta(s%ny+1,wp%tslen,s%ntheta)) 02470 Ampzeta=0.d0 02471 endif 02472 ! Allocate memory for energy time series 02473 allocate(E_tdir(s%ny+1,wp%tslen,s%ntheta)) 02474 E_tdir=0.0d0 02475 E_tdir=0.5d0*par%rho*par%g*Ampzeta**2 02476 E_tdir=E_tdir/s%dtheta 02477 if (par%wbcEvarreduce<1.d0-1d-10) then 02478 do itheta=1,s%ntheta 02479 do iy=1,s%ny+1 02480 emean = sum(E_tdir(iy,:,itheta))/wp%tslen 02481 E_tdir(iy,:,itheta) = par%wbcEvarreduce*(E_tdir(iy,:,itheta)-emean) + emean 02482 enddo 02483 enddo 02484 endif 02485 02486 ! 02487 ! ! Ensure we scale back to the correct Hm0 02488 ! do iy=1,s%ny+1 02489 ! stdeta = sum(E_tdir(iy,:,:))*s%dtheta ! sum energy 02490 ! stdeta = stdeta/wp%tslen ! mean energy 02491 ! 02492 ! stdzeta = (wp%Hm0interp(iy)/sqrt(2.d0))**2 * par%rhog8 02493 ! 02494 ! E_tdir(iy,:,:) = E_tdir(iy,:,:)*stdzeta/stdeta 02495 ! enddo 02496 02497 ! Print directional energy distribution to screen 02498 if (par%bclwonly==0) then 02499 etot = sum(E_tdir) 02500 do itheta=1,s%ntheta 02501 perc = sum(E_tdir(:,:,itheta))/etot*100 02502 call writelog('ls','(a,i0,a,f0.2,a)','Wave bin ',itheta,' contains ',perc,'% of total energy') 02503 enddo 02504 endif 02505 02506 ! Open file for storage 02507 call writelog('ls','','Writing wave energy to ',trim(wp%Efilename),' ...') 02508 inquire(iolength=reclen) 1.d0 02509 reclen=reclen*(s%ny+1)*(s%ntheta) 02510 fid = create_new_fid() 02511 open(fid,file=trim(wp%Efilename),form='unformatted',access='direct',recl=reclen,status='REPLACE') 02512 02513 ! Write to external file 02514 if (.not. allocated(E_t)) allocate(E_t(wp%tslen)) 02515 if (wp%dtchanged) then 02516 ! Interpolate from internal time axis to output time axis 02517 allocate(E_interp(s%ny+1,wp%tslenbc,s%ntheta)) 02518 do itheta=1,s%ntheta 02519 do iy=1,s%ny+1 02520 E_t=E_tdir(iy,:,itheta) 02521 do it=1,wp%tslenbc 02522 call linear_interp(wp%tin,E_t,wp%tslen,(it-1)*wp%dtbc,E_interp(iy,it,itheta),status) 02523 enddo 02524 enddo 02525 enddo 02526 ! write to file 02527 do irec=1,wp%tslenbc+1 02528 write(fid,rec=irec)E_interp(:,min(irec,wp%tslenbc),:) 02529 end do 02530 deallocate(E_interp) 02531 else 02532 ! no need for interpolation 02533 do irec=1,wp%tslenbc+1 02534 write(fid,rec=irec)E_tdir(:,min(irec,wp%tslenbc),:) 02535 end do 02536 endif 02537 close(fid) 02538 call writelog('sl','','file done') 02539 02540 if (par%bclwonly==0) then 02541 ! Free memory 02542 deallocate(zeta,Ampzeta,E_tdir, Amp, eta) 02543 else 02544 deallocate(Ampzeta,E_tdir) 02545 endif 02546 02547 end subroutine generate_ebcf 02548 02549 02550 ! -------------------------------------------------------------- 02551 ! ------ Calculate time series of short wave at offshore ------- 02552 ! --------------------------- boundary ------------------------- 02553 subroutine generate_swts(wp,s,par) 02554 use params, only: parameters 02555 use logging_module, only: writelog, progress_indicator 02556 use spaceparamsdef, only: spacepars 02557 02558 implicit none 02559 02560 ! input/output 02561 type(waveparamsnew),intent(inout) :: wp 02562 type(spacepars),intent(in) :: s 02563 type(parameters),intent(in) :: par 02564 02565 ! internal 02566 integer :: j,it,ik 02567 real*8 :: u1, u2, z, U, dU 02568 02569 ! allocate memory for time series of data 02570 allocate(wp%zsits(s%ny+1,wp%tslen)) 02571 allocate(wp%uits(s%ny+1,wp%tslen)) 02572 allocate(wp%vits(s%ny+1,wp%tslen)) 02573 allocate(wp%duits(s%ny+1,wp%tslen)) 02574 allocate(wp%dvits(s%ny+1,wp%tslen)) 02575 wp%zsits=0.d0 02576 wp%uits=0.d0 02577 wp%vits=0.d0 02578 wp%duits=0.d0 02579 wp%dvits=0.d0 02580 if (par%bclwonly==0) then 02581 02582 u1 = 0.d0 02583 u2 = 0.d0 02584 U = 0.d0 02585 dU = 0.d0 02586 ! z-level of layer (negative downwards) 02587 z = -1 * (wp%h0 - par%nhlay * wp%h0) 02588 ! total surface elevation 02589 call writelog('ls','','Calculating short wave elevation time series') 02590 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02591 do it=1,wp%tslen 02592 call progress_indicator(.false.,dble(it)/wp%tslen*100,5.d0,2.d0) 02593 do ik=1,wp%K 02594 if (wp%PRindex(ik)==1) then 02595 do j=1,s%ny+1 02596 wp%zsits(j,it)=wp%zsits(j,it)+wp%A(j,ik)*dsin( & 02597 +wp%wgen(ik)*wp%tin(it)& 02598 -wp%kgen(ik)*( dsin(wp%thetagen(ik))*(s%yz(1,j)-s%yz(1,1)) & 02599 +dcos(wp%thetagen(ik))*(s%xz(1,j)-s%xz(1,1))) & 02600 +wp%phigen(ik) & 02601 ) 02602 enddo 02603 endif 02604 enddo 02605 enddo 02606 ! depth-averaged velocity 02607 call writelog('ls','','Calculating short wave velocity time series') 02608 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02609 do it=1,wp%tslen 02610 call progress_indicator(.false.,dble(it)/wp%tslen*100,5.d0,2.d0) 02611 do ik=1,wp%K 02612 if (wp%PRindex(ik)==1) then 02613 do j=1,s%ny+1 02614 02615 if ( par%nonhq3d == 1 .and. par%nhlay>0 ) then 02616 ! Compute layer averaged velocity for layer 1 and 2 based on layer level z. 02617 u1 = wp%wgen(ik) * wp%A(j,ik) / (dsinh(wp%kgen(ik)*wp%h0)*wp%kgen(ik))*(dsinh(wp%kgen(ik)*(z + wp%h0))) * & 02618 dsin(wp%wgen(ik)*wp%tin(it)& 02619 -wp%kgen(ik)*( dsin(wp%thetagen(ik))*(s%yz(1,j)-s%yz(1,1)) & 02620 +dcos(wp%thetagen(ik))*(s%xz(1,j)-s%xz(1,1))) & 02621 +wp%phigen(ik)) 02622 u1 = u1 / (wp%h0+z) 02623 u2 = wp%wgen(ik) * wp%A(j,ik) / (dsinh(wp%kgen(ik)*wp%h0)*wp%kgen(ik))*(dsinh(wp%kgen(ik)*(0 + wp%h0)) - & 02624 dsinh(wp%kgen(ik) * (z + wp%h0))) * & 02625 dsin(wp%wgen(ik)*wp%tin(it) & 02626 -wp%kgen(ik)*( dsin(wp%thetagen(ik))*(s%yz(1,j)-s%yz(1,1)) & 02627 +dcos(wp%thetagen(ik))*(s%xz(1,j)-s%xz(1,1))) & 02628 +wp%phigen(ik)) 02629 u2 = u2/(z*-1) 02630 ! Store U. U = alpha * U1 + (1-alpha) * U2 02631 U = par%nhlay * u1 + (1-par%nhlay) * u2 02632 ! Store dU. Du = U1 - U2 02633 dU = (u1 - u2) 02634 02635 ! Eastward component U: 02636 wp%uits(j,it) = wp%uits(j,it) + dcos(wp%thetagen(ik))*U 02637 ! Northward component V: 02638 wp%vits(j,it) = wp%vits(j,it) + dsin(wp%thetagen(ik))*U 02639 ! Eastward component dU: 02640 wp%duits(j,it) = wp%duits(j,it) + dcos(wp%thetagen(ik))*dU 02641 ! Northward component dV: 02642 wp%dvits(j,it) = wp%dvits(j,it) + dsin(wp%thetagen(ik))*dU 02643 else 02644 ! Depth-average velocity in wave direction: 02645 U = 1.d0/wp%h0*wp%wgen(ik)*wp%A(j,ik) * & 02646 dsin(wp%wgen(ik)*wp%tin(it) & 02647 -wp%kgen(ik)*( dsin(wp%thetagen(ik))*(s%yz(1,j)-s%yz(1,1)) & 02648 +dcos(wp%thetagen(ik))*(s%xz(1,j)-s%xz(1,1))) & 02649 +wp%phigen(ik) & 02650 ) * & 02651 1.d0/wp%kgen(ik) 02652 02653 ! Eastward component: 02654 wp%uits(j,it) = wp%uits(j,it) + dcos(wp%thetagen(ik))*U 02655 ! Northward component: 02656 wp%vits(j,it) = wp%vits(j,it) + dsin(wp%thetagen(ik))*U 02657 endif 02658 enddo 02659 endif 02660 enddo 02661 enddo 02662 ! 02663 ! 02664 ! Apply tapering to time series 02665 do j=1,s%ny+1 02666 wp%uits(j,:)=wp%uits(j,:)*wp%taperf 02667 wp%vits(j,:)=wp%vits(j,:)*wp%taperf 02668 wp%zsits(j,:)=wp%zsits(j,:)*wp%taperf 02669 wp%duits(j,:)=wp%duits(j,:)*wp%taperf 02670 wp%dvits(j,:)=wp%dvits(j,:)*wp%taperf 02671 enddo 02672 else 02673 ! do nothing in this routine 02674 endif 02675 end subroutine generate_swts 02676 02677 02678 ! -------------------------------------------------------------- 02679 ! ----------------------- Bound long wave ---------------------- 02680 ! -------------------------------------------------------------- 02681 subroutine generate_qbcf(wp,s,par) 02682 02683 use params, only: parameters 02684 use spaceparamsdef, only: spacepars 02685 use logging_module, only: writelog, progress_indicator 02686 use math_tools, only: fftkind, fft, flipiv 02687 use interp, only: linear_interp 02688 use filefunctions, only: create_new_fid 02689 02690 implicit none 02691 02692 ! input/output 02693 type(waveparamsnew),intent(inout) :: wp 02694 type(spacepars),intent(in) :: s 02695 type(parameters),intent(in) :: par 02696 ! internal 02697 integer :: j,m,iq,irec,it ! counters 02698 integer :: K ! copy of K 02699 integer :: halflen,reclen,fid,status 02700 logical :: firsttime ! used for output message only 02701 real*8 :: deltaf,qmean ! difference frequency 02702 real*8,dimension(:), allocatable :: term1,term2,term2new,dif,chk1,chk2 02703 real*8,dimension(:,:),allocatable :: Eforc,D,deltheta,KKx,KKy,dphi3,k3,cg3,theta3,Abnd,D_sign 02704 real*8,dimension(:,:,:),allocatable :: q,qinterp 02705 complex(fftkind),dimension(:),allocatable :: Comptemp,Comptemp2,Gn 02706 complex(fftkind),dimension(:,:,:),allocatable:: Ftemp 02707 02708 02709 02710 ! This function has changed with respect to previous versions of XBeach, in that 02711 ! the bound long wave has to be calculated separately at each longshore point, 02712 ! owing to longshore varying incident wave spectra 02713 02714 ! shortcut variable 02715 K = wp%K 02716 02717 ! Print message to screen 02718 call writelog('sl','', 'Calculating primary wave interaction') 02719 02720 ! Allocate two-dimensional variables for all combinations of interacting wave 02721 ! components to be filled triangular 02722 allocate(Eforc(K-1,K)) 02723 allocate(D(K-1,K)) 02724 allocate(deltheta(K-1,K)) 02725 allocate(KKx(K-1,K)) 02726 allocate(KKy(K-1,K)) 02727 allocate(dphi3(K-1,K)) 02728 allocate(k3(K-1,K)) 02729 allocate(cg3(K-1,K)) 02730 allocate(theta3(K-1,K)) 02731 ! Allocate variables for amplitude and Fourier coefficients of bound long wave 02732 allocate(Gn(wp%tslen)) 02733 allocate(Abnd(K-1,K)) 02734 allocate(Ftemp(K-1,K,4)) ! Jaap qx, qy qtot, zeta 02735 ! Storage for output discharge 02736 allocate(q(s%ny+1,wp%tslen,4)) ! qx qy qtot, zeta 02737 ! 02738 ! Initialize variables as zero 02739 Eforc = 0 02740 D = 0 02741 deltheta = 0 02742 KKx = 0 02743 KKy = 0 02744 dphi3 = 0 02745 k3 = 0 02746 cg3 = 0 02747 theta3 = 0 02748 Gn = 0*compi 02749 Abnd = 0 02750 Ftemp = 0*compi 02751 q = 0 02752 02753 ! First time is set true for each time new wave bc are generated 02754 firsttime = .true. 02755 02756 ! upper half of frequency space 02757 halflen = wp%tslen/2 02758 02759 02760 02761 ! Run loop over wave-wave interaction components 02762 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02763 do m=1,K-1 02764 ! call writelog('ls','(a,i0,a,i0)','Wave component ',m,' of ',K-1) 02765 call progress_indicator(.false.,dble(m)/(K-1)*100,5.d0,2.d0) 02766 ! Allocate memory 02767 allocate(term1(K-m),term2(K-m),term2new(K-m),dif(K-m),chk1(K-m),chk2(K-m)) 02768 02769 ! Determine difference frequency 02770 deltaf=m*wp%dfgen 02771 02772 ! Determine difference angles (pi already added) 02773 deltheta(m,1:K-m) = abs(wp%thetagen(m+1:K)-wp%thetagen(1:K-m))+par%px 02774 02775 ! Determine x- and y-components of wave numbers of difference waves 02776 KKy(m,1:K-m)=wp%kgen(m+1:K)*dsin(wp%thetagen(m+1:K))-wp%kgen(1:K-m)*dsin(wp%thetagen(1:K-m)) 02777 KKx(m,1:K-m)=wp%kgen(m+1:K)*dcos(wp%thetagen(m+1:K))-wp%kgen(1:K-m)*dcos(wp%thetagen(1:K-m)) 02778 02779 ! Determine difference wave numbers according to Van Dongeren et al. 2003 02780 ! eq. 19 (+cos due to pi in deltheta) 02781 k3(m,1:K-m) =sqrt(wp%kgen(1:K-m)**2+wp%kgen(m+1:K)**2+ & 02782 2*wp%kgen(1:K-m)*wp%kgen(m+1:K)*dcos(deltheta(m,1:K-m))) 02783 02784 ! Determine group velocity of difference waves 02785 cg3(m,1:K-m)= 2.d0*par%px*deltaf/k3(m,1:K-m) 02786 02787 ! Modification Robert + Jaap: make sure that the bound long wave amplitude does not 02788 ! explode when offshore boundary is too close to shore, 02789 ! by limiting the interaction group velocity 02790 cg3(m,1:K-m) = min(cg3(m,1:K-m),par%nmax*sqrt(par%g/k3(m,1:K-m)*tanh(k3(m,1:K-m)*wp%h0))) 02791 02792 ! Determine difference-interaction coefficient according to Okihiro et al eq. 4a 02793 ! Instead of Herbers 1994 this coefficient is derived for the surface elavtion in stead of the bottom pressure. 02794 term1 = (-wp%wgen(1:K-m))*wp%wgen(m+1:K) 02795 term2 = (-wp%wgen(1:K-m))+wp%wgen(m+1:K) 02796 term2new = cg3(m,1:K-m)*k3(m,1:K-m) 02797 dif = (abs(term2-term2new)) 02798 if (any(dif>0.01*term2) .and. firsttime) then 02799 firsttime = .false. 02800 ! call writelog('lws','','Warning: shallow water so long wave variance is reduced using par%nmax') 02801 endif 02802 chk1 = cosh(wp%kgen(1:K-m)*wp%h0) 02803 chk2 = cosh(wp%kgen(m+1:K)*wp%h0) 02804 02805 D(m,1:K-m) = -par%g*wp%kgen(1:K-m)*wp%kgen(m+1:K)*dcos(deltheta(m,1:K-m))/2.d0/term1+ & 02806 +term2**2/(par%g*2)+par%g*term2/ & 02807 ((par%g*k3(m,1:K-m)*tanh(k3(m,1:K-m)*wp%h0)-(term2new)**2)*term1)* & 02808 (term2*((term1)**2/par%g/par%g - wp%kgen(1:K-m)*wp%kgen(m+1:K)*dcos(deltheta(m,1:K-m))) & 02809 - 0.50d0*((-wp%wgen(1:K-m))*wp%kgen(m+1:K)**2/(chk2**2)+wp%wgen(m+1:K)*wp%kgen(1:K-m)**2/(chk1**2))) 02810 02811 ! Correct for surface elevation input and output instead of bottom pressure 02812 ! so it is consistent with Van Dongeren et al 2003 eq. 18. 02813 !No need to correct for bottom pressure anymore. 02814 !D(m,1:K-m) = D(m,1:K-m)*cosh(k3(m,1:K-m)*wp%h0)/(cosh(wp%kgen(1:K-m)*wp%h0)*cosh(wp%kgen(m+1:K)*wp%h0)) 02815 02816 02817 ! Exclude interactions with components smaller than or equal to current 02818 ! component according to lower limit Herbers 1994 eq. 1 02819 where(wp%fgen<=deltaf) 02820 D(m,:)=0.d0 02821 endwhere 02822 02823 ! Exclude interactions with components that are cut-off by the fcutoff 02824 ! parameter 02825 if (deltaf<=par%fcutoff) D(m,:)=0.d0 02826 02827 ! Determine phase of bound long wave assuming a local equilibrium with 02828 ! forcing of interacting primary waves according to Van Dongeren et al. 02829 ! 2003 eq. 21 (the angle is the imaginary part of the natural log of a 02830 ! complex number as long as the complex number is not zero) 02831 allocate(Comptemp(K-m),Comptemp2(K-m)) 02832 Comptemp=conjg(wp%CompFn(1,wp%Findex(1)+m:wp%Findex(1)+K-1)) 02833 Comptemp2=conjg(wp%CompFn(1,wp%Findex(1):wp%Findex(1)+K-m-1)) 02834 dphi3(m,1:K-m) = par%px+imag(log(Comptemp))-imag(log(Comptemp2)) 02835 ! Menno removed pi, because it is incorporated in the sign of the amplitude 02836 dphi3(m,1:K-m) = imag(log(Comptemp))-imag(log(Comptemp2)) 02837 deallocate (Comptemp,Comptemp2) 02838 ! 02839 ! Determine angle of bound long wave according to Van Dongeren et al. 2003 eq. 22 02840 theta3 = atan2(KKy,KKx) 02841 ! 02842 ! free memory 02843 deallocate(term1,term2,term2new,dif,chk1,chk2) 02844 enddo ! m=1,K-1 02845 ! 02846 ! Output to screen 02847 if (.not. firsttime) then 02848 call writelog('lws','','Warning: shallow water so long wave variance is reduced using par%nmax') 02849 endif 02850 call writelog('sl','', 'Calculating flux at boundary') 02851 ! 02852 ! Allocate temporary arrays for upcoming loop 02853 allocate(Comptemp(halflen-1)) 02854 allocate(Comptemp2(wp%tslen)) 02855 ! 02856 ! Run a loop over the offshore boundary 02857 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02858 do j=1,s%ny+1 02859 ! Determine energy of bound long wave according to Herbers 1994 eq. 1 based 02860 ! on difference-interaction coefficient and energy density spectra of 02861 ! primary waves 02862 ! Robert: E = 2*D**2*S**2*dtheta**2*df can be rewritten as 02863 ! E = 2*D**2*Sf**2*df 02864 Eforc = 0 02865 do m=1,K-1 02866 ! Menno: Use the forced short wave energy to compute the bound long wave energy 02867 ! instead of the corrected energy (Sfinterpq) 02868 if (par%Sfold==1) then 02869 Eforc(m,1:K-m) = 2*D(m,1:K-m)**2*wp%Sfinterpq(j,1:K-m)*wp%Sfinterpq(j,m+1:K)*wp%dfgen 02870 else 02871 Eforc(m,1:K-m) = 2*D(m,1:K-m)**2*wp%Sfinterp(j,1:K-m)*wp%Sfinterp(j,m+1:K)*wp%dfgen 02872 endif 02873 enddo 02874 ! 02875 ! Calculate bound wave amplitude for this offshore grid point 02876 !Abnd = sqrt(2*Eforc*wp%dfgen) 02877 02878 ! Menno: add the sign of the interaction coefficient in the amptlitude. Large dtheta can result in a positive D. 02879 ! The phase of the bound wave is now only determined by phi1+phi2 02880 allocate(D_sign(K-1,K)) 02881 D_sign = 1 02882 ! Menno: put the sign of D in front of D_sign 02883 D_sign = sign(D_sign,D) 02884 ! Multiply amplitude with the sign of D 02885 Abnd = sqrt(2*Eforc*wp%dfgen) * D_sign 02886 deallocate(D_sign) 02887 ! 02888 ! Determine complex description of bound long wave per interaction pair of 02889 ! primary waves for first y-coordinate along seaside boundary 02890 Ftemp(:,:,1) = Abnd/2*exp(-1*compi*dphi3)*cg3*dcos(theta3) ! qx 02891 Ftemp(:,:,2) = Abnd/2*exp(-1*compi*dphi3)*cg3*dsin(theta3) ! qy 02892 Ftemp(:,:,3) = Abnd/2*exp(-1*compi*dphi3)*cg3 ! qtot 02893 Ftemp(:,:,4) = Abnd/2*exp(-1*compi*dphi3) ! eta 02894 ! 02895 ! loop over qx,qy and qtot 02896 do iq=1,4 02897 ! Unroll wave component to correct place along the offshore boundary 02898 Ftemp(:,:,iq) = Ftemp(:,:,iq)* & 02899 exp(-1*compi*(KKy*(s%yz(1,j)-s%yz(1,1))+KKx*(s%xz(1,j)-s%xz(1,1)))) 02900 ! Determine Fourier coefficients 02901 Gn(2:K) = sum(Ftemp(:,:,iq),DIM=2) 02902 Comptemp = conjg(Gn(2:halflen)) 02903 call flipiv(Comptemp,halflen-1) 02904 Gn(halflen+2:wp%tslen) = Comptemp 02905 ! 02906 ! Print status message to screen 02907 if (iq==3) then 02908 !call writelog('ls','(A,I0,A,I0)','Flux ',j,' of ',s%ny+1) 02909 call progress_indicator(.false.,dble(j)/(s%ny+1)*100,5.d0,2.d0) 02910 endif 02911 ! 02912 ! Inverse Discrete Fourier transformation to transform back to time space 02913 ! from frequency space. (Fortran ifft is scaled by sqrt(N) ) 02914 Comptemp2=fft(Gn,inv=.true.) 02915 ! 02916 ! Determine mass flux as function of time and let the flux gradually 02917 ! increase and decrease in and out the wave time record using the earlier 02918 ! specified window. Sqrt(N) is due to scaling of the ifft. 02919 Comptemp2=Comptemp2/sqrt(dble(wp%tslen)) 02920 q(j,:,iq)=dreal(Comptemp2*wp%tslen)*wp%taperf 02921 enddo ! iq=1,4 !Ap Menno 27-06-2018 02922 enddo ! j=1,s%ny+1 02923 ! 02924 ! free memory 02925 deallocate(Comptemp,Comptemp2,Ftemp) 02926 ! 02927 ! 02928 if (par%wbcQvarreduce<1.d0-1d-10) then 02929 do iq=1,4 02930 do j=1,s%ny+1 02931 qmean = sum(q(j,:,iq))/wp%tslen 02932 q(j,:,iq) = par%wbcQvarreduce*(q(j,:,iq)-qmean) + qmean 02933 enddo 02934 enddo 02935 endif 02936 ! 02937 ! 02938 ! 02939 ! Now write data to file 02940 if (par%nonhspectrum==0) then 02941 ! If doing combined wave action balance and swell wave flux with swkhmin>0 then we need to add short wave velocity 02942 ! to time series of q here: 02943 if (par%swkhmin>0.d0) then 02944 do j=1,s%ny+1 02945 q(j,:,1) = q(j,:,1) + wp%uits(j,:)*wp%h0 ! u flux 02946 q(j,:,4) = q(j,:,4) + wp%zsits(j,:) ! eta 02947 enddo 02948 endif 02949 02950 if (par%order==1) then 02951 do j=1,s%ny+1 02952 q(j,:,1) = q(j,:,1) *0 ! u flux 02953 q(j,:,4) = q(j,:,4) *0 ! eta 02954 enddo 02955 endif 02956 ! 02957 ! 02958 ! Open file for storage of bound long wave flux 02959 call writelog('ls','','Writing long wave mass flux to ',trim(wp%qfilename),' ...') 02960 inquire(iolength=reclen) 1.d0 02961 reclen=reclen*((s%ny+1)*4) 02962 fid = create_new_fid() 02963 open(fid,file=trim(wp%qfilename),form='unformatted',access='direct',recl=reclen,status='REPLACE') 02964 ! 02965 ! Write to external file 02966 if (wp%dtchanged) then 02967 ! Interpolate from internal time axis to output time axis 02968 allocate(qinterp(s%ny+1,wp%tslenbc,3)) 02969 do iq=1,3 02970 do it=1,wp%tslenbc 02971 do j=1,s%ny+1 02972 call linear_interp(wp%tin,q(j,:,iq),wp%tslen,(it-1)*wp%dtbc,qinterp(j,it,iq),status) 02973 enddo 02974 enddo 02975 enddo 02976 ! write to file 02977 do irec=1,wp%tslenbc+1 02978 write(fid,rec=irec)qinterp(:,min(irec,wp%tslenbc),:) 02979 end do 02980 deallocate(qinterp) 02981 else 02982 ! no need for interpolation 02983 do irec=1,wp%tslenbc+1 02984 write(fid,rec=irec)q(:,min(irec,wp%tslenbc),:) 02985 end do 02986 endif 02987 close(fid) 02988 call writelog('sl','','file done') 02989 else 02990 do j=1,s%ny+1 02991 ! add to velocity time series 02992 wp%uits(j,:)=wp%uits(j,:)+q(j,:,1)/wp%h0 02993 wp%vits(j,:)=wp%vits(j,:)+q(j,:,2)/wp%h0 02994 ! add to surface elevation time series 02995 wp%zsits(j,:)=wp%zsits(j,:)+q(j,:,4) 02996 enddo 02997 endif ! par%nonhspectrum==1 02998 02999 ! Free memory 03000 deallocate(Eforc,D,deltheta,KKx,KKy,dphi3,k3,cg3,theta3,Gn,Abnd,q) 03001 03002 end subroutine generate_qbcf 03003 03004 03005 ! -------------------------------------------------------------- 03006 ! --------------- Non-hydrostatic wave time -------------------- 03007 ! ---------------- series file generation ---------------------- 03008 subroutine generate_nhtimeseries_file(wp,par) 03009 03010 use params 03011 use filefunctions, only: create_new_fid 03012 use logging_module 03013 03014 implicit none 03015 03016 ! input/output 03017 type(waveparamsnew),intent(inout) :: wp 03018 type(parameters),intent(in) :: par 03019 ! internal 03020 integer :: fid 03021 integer :: it 03022 character(LEN=256) :: rowfmt 03023 03024 03025 call writelog('ls','','Writing short wave time series to ',wp%nhfilename) 03026 03027 fid = create_new_fid() 03028 open(fid,file=trim(wp%nhfilename),status='REPLACE') 03029 write(fid,'(a)')'vector' 03030 write(fid,'(a)')'6' 03031 write(fid,'(a)')'t,U,V,Z,dU,dV' 03032 ! write format descriptor for file to internal string 03033 write(rowfmt,'(A,I4,A)') '(',1+5*(par%ny+1),'(1X,E18.9))' 03034 ! Add one extra point at the start of the time series to deal with interpolation errors 03035 write(fid,rowfmt)par%t-par%dt*par%maxdtfac,wp%uits(:,1),wp%vits(:,1), & 03036 wp%zsits(:,1),wp%duits(:,1),wp%dvits(:,1) 03037 ! write time series to file 03038 do it=1,wp%tslen 03039 write(fid,rowfmt)wp%tin(it)+par%t-par%dt,wp%uits(:,it),wp%vits(:,it),wp%zsits(:,it),wp%duits(:,it),wp%dvits(:,it) 03040 enddo 03041 ! (almost always) add one extra point after the time series to deal with interpolation errors 03042 if (wp%tin(wp%tslen)<wp%rtbc+par%dt*par%maxdtfac) then 03043 write(fid,rowfmt)wp%tin(wp%tslen)+par%t+par%dt*par%maxdtfac,wp%uits(:,wp%tslen),wp%vits(:,wp%tslen), & 03044 wp%zsits(:,wp%tslen),wp%duits(:,wp%tslen),wp%dvits(:,wp%tslen) 03045 endif 03046 close(fid) 03047 03048 end subroutine generate_nhtimeseries_file 03049 03050 03051 subroutine set_stationary_spectrum (s,par,wp,combspec) 03052 use params, only: parameters 03053 use spaceparamsdef, only: spacepars 03054 use interp, only: interp_using_trapez_rule 03055 use filefunctions, only: create_new_fid 03056 use logging_module, only: writelog 03057 type(parameters) :: par 03058 type(spacepars) :: s 03059 type(spectrum) :: combspec 03060 type(waveparamsnew),intent(in) :: wp 03061 real*8 :: xcycle 03062 real*8, dimension(:), allocatable :: angcart,Sdcart 03063 integer :: j 03064 integer :: reclen,fid 03065 03066 allocate(angcart(combspec%nang)) 03067 allocate(Sdcart(combspec%nang)) 03068 03069 xcycle=2*par%px 03070 ! combspec%ang contains nautical directions from 0 to 2 pi; convert to cartesian from -3/2pi to 1/2 pi 03071 ! in reverse order 03072 03073 call interp_using_trapez_rule(combspec%ang,combspec%Sd,combspec%nang,xcycle,s%theta_s,s%ee_s(1,1,:),s%ntheta_s) 03074 do j=1,s%ny+1 03075 s%ee_s(1,j,:)=s%ee_s(1,1,:) 03076 end do 03077 03078 s%ee_s(1,:,:)=s%ee_s(1,:,:)*par%rho*par%g 03079 03080 call writelog('ls','','Writing stationary wave energy directional spread to ',trim(wp%Esfilename),' ...') 03081 inquire(iolength=reclen) 1.d0 03082 reclen=reclen*(s%ny+1)*(s%ntheta_s) 03083 fid = create_new_fid() 03084 open(fid,file=trim(wp%Esfilename),form='unformatted',access='direct',recl=reclen,status='REPLACE') 03085 write(fid,rec=1)s%ee_s(1,:,:) 03086 close(fid) 03087 deallocate(angcart) 03088 deallocate(Sdcart) 03089 03090 end subroutine set_stationary_spectrum 03091 03092 subroutine generate_admin_files(par,wp,rtbc_local,dtbc_local,maindir_local,spectrumendtimeold,spectrumendtime) 03093 03094 use filefunctions, only: create_new_fid 03095 use params, only: parameters 03096 implicit none 03097 03098 type(parameters),intent(in) :: par 03099 type(waveparamsnew),intent(in) :: wp 03100 real*8,intent(in) :: rtbc_local,dtbc_local,maindir_local,spectrumendtimeold,spectrumendtime 03101 integer :: iostat 03102 03103 if (par%nonhspectrum==0) then 03104 if (bccount==1) then 03105 fidelist = create_new_fid() 03106 open(fidelist,file='ebcflist.bcf',form='formatted',status='replace') 03107 else 03108 open(fidelist,file='ebcflist.bcf',form='formatted',position='append') 03109 end if 03110 write(fidelist,'(f12.3,a,f12.3,a,f9.3,a,f9.5,a,f11.5,a)') & 03111 & spectrumendtime,' ',rtbc_local,' ',dtbc_local,' ',par%Trep,' ',maindir_local,' '//trim(wp%Efilename) 03112 close(fidelist) 03113 03114 if (bccount==1) then 03115 fidqlist = create_new_fid() 03116 open(fidqlist,file='qbcflist.bcf',form='formatted',status='replace') 03117 else 03118 open(fidqlist,file='qbcflist.bcf',form='formatted',position='append') 03119 end if 03120 write(fidqlist,'(f12.3,a,f12.3,a,f9.3,a,f9.5,a,f11.5,a)') & 03121 & spectrumendtime,' ',rtbc_local,' ',dtbc_local,' ',par%Trep,' ',maindir_local,' '//trim(wp%qfilename) 03122 close(fidqlist) 03123 03124 if(par%single_dir==1) then 03125 if (bccount==1) then 03126 fideslist = create_new_fid() 03127 open(fideslist,file='esbcflist.bcf',form='formatted',status='replace') 03128 else 03129 open(fideslist,file='esbcflist.bcf',form='formatted',position='append') 03130 end if 03131 write(fideslist,'(f12.3,a,f12.3,a,f9.3,a,f9.5,a,f11.5,a)') & 03132 & spectrumendtime,' ',rtbc_local,' ',dtbc_local,' ',par%Trep,' ',maindir_local,' '//trim(wp%Esfilename) 03133 close(fideslist) 03134 endif 03135 else 03136 if (bccount==1) then 03137 fidnhlist = create_new_fid() 03138 open(fidnhlist,file='nhbcflist.bcf',form='formatted',status='replace',iostat=iostat) 03139 else 03140 open(fidnhlist,file='nhbcflist.bcf',form='formatted',position='append',iostat=iostat) 03141 end if 03142 write(fidnhlist,'(f12.3,a,f12.3,a,f12.3,a)',iostat=iostat) & 03143 & spectrumendtimeold,' ',spectrumendtime,' ',par%Trep,' '//trim(wp%nhfilename) 03144 close(fidnhlist) 03145 end if 03146 end subroutine generate_admin_files 03147 03148 ! Generate time series of bound super harmonics 03149 subroutine generate_secondorder(par,s, wp) 03150 03151 use params, only: parameters 03152 use spaceparamsdef, only: spacepars 03153 use logging_module, only: writelog, progress_indicator 03154 use math_tools, only: fftkind, fft, flipiv 03155 use interp, only: linear_interp 03156 use filefunctions, only: create_new_fid 03157 use spaceparamsdef, only: spacepars 03158 03159 implicit none 03160 03161 ! input/output 03162 type(waveparamsnew),intent(inout) :: wp 03163 type(spacepars),intent(in) :: s 03164 type(parameters),intent(in) :: par 03165 ! internal 03166 integer :: j,m,iq,irec,it,i,count ! counters 03167 integer :: K ! copy of K 03168 integer :: halflen,reclen,fid,status 03169 logical :: firsttime ! used for output message only 03170 real*8 :: deltaf,z 03171 real*8,dimension(:), allocatable :: term1,term2,term2new,dif,chk1,chk2 03172 real*8,dimension(:,:),allocatable :: Eforc,D,deltheta,KKx,KKy,dphi3,k3,w3,c,theta3,Abnd,qx1,qx2,qy1,qy2,D_sign 03173 real*8,dimension(:,:,:),allocatable :: q,qinterp 03174 complex(fftkind),dimension(:),allocatable :: Comptemp,Comptemp2,Gn 03175 complex(fftkind),dimension(:,:,:),allocatable:: Ftemp 03176 integer,dimension(:), allocatable :: ii,jj 03177 03178 03179 03180 03181 ! shortcut variable 03182 K = wp%K 03183 03184 ! Print message to screen 03185 call writelog('sl','', 'Calculating primary wave interaction - super harmonics') 03186 03187 ! Allocate two-dimensional variables for all combinations of interacting wave 03188 ! components to be filled triangular 03189 allocate(Eforc(2*K,2*K)) 03190 allocate(D(2*K,2*K)) 03191 allocate(deltheta(2*K,2*K)) 03192 allocate(KKx(2*K,2*K)) 03193 allocate(KKy(2*K,2*K)) 03194 allocate(dphi3(2*K,2*K)) 03195 allocate(k3(2*K,2*K)) 03196 allocate(w3(2*K,2*K)) 03197 allocate(c(2*K,2*K)) 03198 allocate(theta3(2*K,2*K)) 03199 ! Allocate variables for amplitude and Fourier coefficients of bound wave 03200 allocate(Gn(wp%tslen)) 03201 allocate(Abnd(2*K,2*K)) 03202 if (par%nonhq3d==1 .and. par%nhlay>0) then 03203 allocate(Ftemp(2*K,2*K,6)) ! qx1, qx2, qy1, qy2, qtot, zeta 03204 ! Storage for output discharge 03205 allocate(q(s%ny+1,wp%tslen,6)) ! qx1, qx2, qy1, qy2, qtot, zeta 03206 allocate(qx1(2*K,2*K)) 03207 allocate(qx2(2*K,2*K)) 03208 allocate(qy1(2*K,2*K)) 03209 allocate(qy2(2*K,2*K)) 03210 qx1 = 0 03211 qx2 = 0 03212 qy1 = 0 03213 qy2 = 0 03214 else 03215 allocate(Ftemp(2*K,2*K,4)) ! qx, qy qtot, zeta 03216 ! Storage for output discharge 03217 allocate(q(s%ny+1,wp%tslen,4)) ! qx qy qtot, zeta 03218 endif 03219 03220 ! 03221 ! Initialize variables as zero 03222 Eforc = 0 03223 D = 0 03224 deltheta = 0 03225 KKx = 0 03226 KKy = 0 03227 dphi3 = 0 03228 k3 = 0 03229 w3 = 0 03230 c = 0 03231 theta3 = 0 03232 Gn = 0*compi 03233 Abnd = 0 03234 Ftemp = 0*compi 03235 q = 0 03236 i = 0 03237 count = 0 03238 z = 0 03239 03240 03241 03242 ! upper half of frequency space 03243 halflen = wp%tslen/2 03244 03245 ! Run loop over wave-wave interaction components 03246 call progress_indicator(.true.,0.d0,5.d0,2.d0) 03247 do m=2,2*K 03248 call progress_indicator(.false.,dble(m)/(2*K)*100,5.d0,2.d0) 03249 ! Create indices 03250 if (m .LE. K+1) then 03251 allocate(term1(m-1),term2(m-1),term2new(m-1),dif(m-1),chk1(m-1),chk2(m-1)) 03252 allocate(Comptemp(m-1),Comptemp2(m-1)) 03253 allocate(ii(m-1)) 03254 allocate(jj(m-1)) 03255 ii = 0 03256 jj = 0 03257 ! Indices: ii=1:m-1 03258 ! Indices: jj = flip(ii) 03259 do i=1,m-1 03260 ii(i) = i 03261 jj(m - i) = i 03262 enddo 03263 else 03264 allocate(term1(K+1-(m-K)),term2(K+1-(m-K)),term2new(K+1-(m-K)),dif(K+1-(m-K)),chk1(K+1-(m-K)),chk2(K+1-(m-K))) 03265 allocate(Comptemp(K+1-(m-K)),Comptemp2(K+1-(m-K))) 03266 allocate(ii(K+1-(m-K))) 03267 allocate(jj(K+1-(m-K))) 03268 ii = 0 03269 jj = 0 03270 ! Indices: ii=m-K:K 03271 ! Indices: jj = flip(ii) 03272 do i=1,K+1-(m-K) 03273 ii(i) = m-K + i -1 03274 jj(K+1-(m-K) -i + 1) = m-K + i -1 03275 enddo 03276 endif 03277 ! Determine delta theta: 03278 deltheta(m,ii) = abs(wp%thetagen(ii)-wp%thetagen(jj)) 03279 03280 ! Determine x- and y-components of wave numbers of K3 03281 KKy(m,ii)=wp%kgen(ii)*dsin(wp%thetagen(ii))+wp%kgen(jj)*dsin(wp%thetagen(jj)) 03282 KKx(m,ii)=wp%kgen(ii)*dcos(wp%thetagen(ii))+wp%kgen(jj)*dcos(wp%thetagen(jj)) 03283 03284 ! Determine wave numbers according to Van Dongeren et al. 2003 03285 ! eq. 19 03286 k3(m,ii) =sqrt(wp%kgen(ii)**2+wp%kgen(jj)**2+ & 03287 2*wp%kgen(ii)*wp%kgen(jj)*dcos(deltheta(m,ii))) 03288 ! Determine W3 03289 w3(m,ii) = wp%wgen(ii) + wp%wgen(jj) 03290 ! Determine bound group velocity of waves 03291 c(m,ii) = w3(m,ii)/k3(m,ii) 03292 03293 ! Modification Robert + Jaap: make sure that the bound long wave amplitude does not 03294 ! explode when offshore boundary is too close to shore, 03295 ! by limiting the interaction group velocity 03296 !c(m,ii) = min(c(m,ii),par%nmax*sqrt(par%g/k3(m,ii)*tanh(k3(m,ii)*wp%h0))) 03297 !term2new = cg3(m,1:K-m)*k3(m,1:K-m) 03298 !dif = (abs(term2-term2new)) 03299 !if (any(dif>0.01*term2) .and. firsttime) then 03300 ! firsttime = .false. 03301 ! ! call writelog('lws','','Warning: shallow water so long wave variance is reduced using par%nmax') 03302 !endif 03303 term1 = wp%wgen(ii)*wp%wgen(jj) 03304 term2 = (wp%wgen(ii))+wp%wgen(jj) 03305 chk1 = cosh(wp%kgen(ii)*wp%h0) 03306 chk2 = cosh(wp%kgen(jj)*wp%h0) 03307 03308 ! Determine difference-interaction coefficient according to okihiro 1992 03309 ! eq. 4a 03310 D(m,ii) = -par%g*wp%kgen(ii)*wp%kgen(jj)*dcos(deltheta(m,ii))/2.d0/term1+ & 03311 term2**2/(2*par%g) + & 03312 par%g*term2/ & 03313 ((par%g*k3(m,ii)*tanh(k3(m,ii)*wp%h0)-(term2)**2)*term1)* & 03314 (term2*((term1)**2/par%g/par%g - wp%kgen(ii)*wp%kgen(jj)*dcos(deltheta(m,ii))) & 03315 - 0.50d0*((wp%wgen(ii))*wp%kgen(jj)**2/(chk2**2)+wp%wgen(jj)*wp%kgen(ii)**2/(chk1**2))) 03316 03317 ! Menno: limit the higher components, which can not correcty be resolved. 03318 if (k3(m,1)*wp%h0>5) D(m,:)=0.d0 03319 03320 ! Determine phase of bound long wave assuming a local equilibrium with 03321 ! forcing of interacting primary waves according to Van Dongeren et al. 03322 ! 2003 eq. 21 (the angle is the imaginary part of the natural log of a 03323 ! complex number as long as the complex number is not zero) 03324 Comptemp=conjg(wp%CompFn(1,wp%Findex(1)+ii)) 03325 Comptemp2=conjg(wp%CompFn(1,wp%Findex(1)+jj)) 03326 dphi3(m,ii) = imag(log(Comptemp))+imag(log(Comptemp2)) 03327 deallocate (Comptemp,Comptemp2) 03328 ! 03329 ! Determine angle of bound long wave according to Van Dongeren et al. 2003 eq. 22 03330 theta3 = atan2(KKy,KKx) 03331 ! 03332 ! free memory 03333 deallocate(term1,term2,term2new,dif,chk1,chk2) 03334 deallocate(ii,jj) 03335 enddo ! m=1,K-1 03336 ! 03337 ! Output to screen 03338 !if (.not. firsttime) then 03339 ! call writelog('lws','','Warning: shallow water so long wave variance is reduced using par%nmax') 03340 !endif 03341 call writelog('sl','', 'Calculating flux at boundary') 03342 ! 03343 ! Allocate temporary arrays for upcoming loop 03344 allocate(Comptemp(halflen-1)) 03345 allocate(Comptemp2(wp%tslen)) 03346 ! 03347 ! Run a loop over the offshore boundary 03348 do j=1,s%ny+1 03349 ! Determine energy of bound long wave according to Herbers 1994 eq. 1 based 03350 ! on difference-interaction coefficient and energy density spectra of 03351 ! primary waves 03352 ! Robert: E = 2*D**2*S**2*dtheta**2*df can be rewritten as 03353 ! E = 2*D**2*Sf**2*df 03354 Eforc = 0 03355 do m=2,2*K 03356 if (m .LE. K+1) then 03357 allocate(ii(m-1)) 03358 allocate(jj(m-1)) 03359 ii = 0 03360 jj = 0 03361 do i=1,m-1 03362 ii(i) = i 03363 jj(m -i) = i 03364 enddo 03365 else 03366 allocate(ii(K+1-(m-K))) 03367 allocate(jj(K+1-(m-K))) 03368 ii = 0 03369 jj = 0 03370 do i=1,K+1-(m-K) 03371 ii(i) = m-K + i -1 03372 jj(K+1-(m-K) -i + 1) = m-K + i -1 03373 enddo 03374 endif 03375 ! Use Sfinterp instead of sfinterpq 03376 if (par%Sfold==1) then 03377 Eforc(m,ii) = D(m,ii)**2*wp%Sfinterpq(j,ii)*wp%Sfinterpq(j,jj)*wp%dfgen 03378 else 03379 Eforc(m,ii) = D(m,ii)**2*wp%Sfinterp(j,ii)*wp%Sfinterp(j,jj)*wp%dfgen 03380 endif 03381 03382 ! Compute ratio for nonh. Compute here otherwise divide by zero! 03383 if (par%nonhq3d==1) then 03384 ! z-level of layer 03385 z = -1 * (wp%h0 - par%nhlay * wp%h0) 03386 ! ratio qx1/q 03387 qx1(m,ii) = (dsinh(KKx(m,ii) * (z + wp%h0)))/dsinh(KKx(m,ii)*wp%h0) 03388 ! ratio qx2/q 03389 qx2(m,ii) = (dsinh(KKx(m,ii) * wp%h0) - dsinh(KKx(m,ii)))/dsinh(KKx(m,ii)*wp%h0) 03390 ! ratio qy1/q 03391 qy1(m,ii) = (dsinh(KKy(m,ii) * (z + wp%h0)))/dsinh(KKy(m,ii)*wp%h0) 03392 ! ratio qy2/q 03393 qy2(m,ii) = (dsinh(KKy(m,ii) * wp%h0) - dsinh(KKy(m,ii)))/dsinh(KKy(m,ii)*wp%h0) 03394 endif 03395 deallocate(ii,jj) 03396 enddo 03397 ! 03398 ! Calculate bound wave amplitude for this offshore grid point 03399 ! Menno: add the sign of the interaction coefficient in the amptlitude. Large dtheta can result in a positive D. 03400 ! The phase of the bound wave is now only determined by phi1+phi2 03401 allocate(D_sign(2*K,2*K)) 03402 D_sign = 1 03403 ! Menno: put the sign of D in front of D_sign 03404 D_sign = sign(D_sign,D) 03405 Abnd = sqrt(2*Eforc*wp%dfgen) * D_sign 03406 deallocate(D_sign) 03407 ! 03408 ! Determine complex description of bound long wave per interaction pair of 03409 ! primary waves for first y-coordinate along seaside boundary 03410 if (par%nonhq3d==1) then 03411 count = 6 03412 Ftemp(:,:,1) = Abnd/2*exp(-1*compi*dphi3)*c*dcos(theta3)*qx1 ! qx1 03413 Ftemp(:,:,2) = Abnd/2*exp(-1*compi*dphi3)*c*dcos(theta3)*qx2 ! qx2 03414 Ftemp(:,:,3) = Abnd/2*exp(-1*compi*dphi3)*c*dsin(theta3)*qy1 ! qy1 03415 Ftemp(:,:,4) = Abnd/2*exp(-1*compi*dphi3)*c*dsin(theta3)*qy2 ! qy2 03416 Ftemp(:,:,5) = Abnd/2*exp(-1*compi*dphi3)*c ! qtot 03417 Ftemp(:,:,6) = Abnd/2*exp(-1*compi*dphi3) ! eta 03418 else 03419 count = 4 03420 Ftemp(:,:,1) = Abnd/2*exp(-1*compi*dphi3)*c*dcos(theta3) ! qx 03421 Ftemp(:,:,2) = Abnd/2*exp(-1*compi*dphi3)*c*dsin(theta3) ! qy 03422 Ftemp(:,:,3) = Abnd/2*exp(-1*compi*dphi3)*c ! qtot 03423 Ftemp(:,:,4) = Abnd/2*exp(-1*compi*dphi3) ! eta 03424 endif 03425 ! 03426 ! loop over qx,qy and qtot 03427 do iq=1,count 03428 ! Unroll wave component to correct place along the offshore boundary 03429 Ftemp(:,:,iq) = Ftemp(:,:,iq)* & 03430 exp(-1*compi*(KKy*(s%yz(1,j)-s%yz(1,1))+KKx*(s%xz(1,j)-s%xz(1,1)))) 03431 ! Determine Fourier coefficients 03432 Gn(2*wp%Findex(1):(2*wp%Findex(1)-1)+2*K) = sum(Ftemp(:,:,iq),DIM=2) 03433 Comptemp = conjg(Gn(2:halflen)) 03434 call flipiv(Comptemp,halflen-1) 03435 Gn(halflen+2:wp%tslen) = Comptemp 03436 ! 03437 ! Print status message to screen 03438 if (iq==3) then 03439 call writelog('ls','(A,I0,A,I0)','Flux ',j,' of ',s%ny+1) 03440 endif 03441 ! 03442 ! Inverse Discrete Fourier transformation to transform back to time space 03443 ! from frequency space 03444 Comptemp2=fft(Gn,inv=.true.) 03445 ! 03446 ! Determine mass flux as function of time and let the flux gradually 03447 ! increase and decrease in and out the wave time record using the earlier 03448 ! specified window. Fortran FFT is scaled by squar-roor. 03449 Comptemp2=Comptemp2/sqrt(dble(wp%tslen)) 03450 q(j,:,iq)=dreal(Comptemp2*wp%tslen)*wp%taperf 03451 enddo ! iq=1,3 03452 enddo ! j=1,s%ny+1 03453 03454 ! Now write data to file 03455 ! Menno: I don't know if this is appropriate. Now, nonhspectrum is always 1. 03456 if (par%nonhspectrum==0) then 03457 ! If doing combined wave action balance and swell wave flux with swkhmin>0 then we need to add short wave velocity 03458 ! to time series of q here: 03459 if (par%swkhmin>0.d0) then 03460 do j=1,s%ny+1 03461 q(j,:,1) = q(j,:,1) + wp%uits(j,:)*wp%h0 ! u flux 03462 q(j,:,4) = q(j,:,4) + wp%zsits(j,:) ! eta 03463 enddo 03464 endif 03465 ! 03466 ! 03467 ! Open file for storage of bound long wave flux 03468 call writelog('ls','','Writing long wave mass flux to ',trim(wp%qfilename),' ...') 03469 inquire(iolength=reclen) 1.d0 03470 reclen=reclen*((s%ny+1)*4) 03471 fid = create_new_fid() 03472 open(fid,file=trim(wp%qfilename),form='unformatted',access='direct',recl=reclen,status='REPLACE') 03473 ! 03474 ! Write to external file 03475 if (wp%dtchanged) then 03476 ! Interpolate from internal time axis to output time axis 03477 allocate(qinterp(s%ny+1,wp%tslenbc,3)) 03478 do iq=1,3 03479 do it=1,wp%tslenbc 03480 do j=1,s%ny+1 03481 call linear_interp(wp%tin,q(j,:,iq),wp%tslen,(it-1)*wp%dtbc,qinterp(j,it,iq),status) 03482 enddo 03483 enddo 03484 enddo 03485 ! write to file 03486 do irec=1,wp%tslenbc+1 03487 write(fid,rec=irec)qinterp(:,min(irec,wp%tslenbc),:) 03488 end do 03489 deallocate(qinterp) 03490 else 03491 ! no need for interpolation 03492 do irec=1,wp%tslenbc+1 03493 write(fid,rec=irec)q(:,min(irec,wp%tslenbc),:) 03494 end do 03495 endif 03496 close(fid) 03497 call writelog('sl','','file done') 03498 else 03499 do j=1,s%ny+1 03500 if (par%nonhq3d==1 .and. par%nhlay>0) then 03501 ! add to velocity time series 03502 wp%uits(j,:)=wp%uits(j,:)+ (q(j,:,1)/z + q(j,:,2)/(wp%h0+z))/2.d0 03503 ! add to velocity time series 03504 wp%duits(j,:)=wp%duits(j,:)+ (q(j,:,1)/z - q(j,:,2)/(wp%h0+z)) 03505 ! add to surface elevation time series 03506 wp%zsits(j,:)=wp%zsits(j,:)+q(j,:,6) 03507 if (j==s%ny+1) then 03508 deallocate(Eforc,D,deltheta,KKx,KKy,dphi3,k3,c,theta3,Gn,Abnd,w3,qx1,qx2,qy1,qy2,Ftemp) 03509 endif 03510 03511 else 03512 ! add to velocity time series 03513 wp%uits(j,:)=wp%uits(j,:)+q(j,:,1)/wp%h0 03514 ! add to surface elevation time series 03515 wp%zsits(j,:)=wp%zsits(j,:)+q(j,:,4) 03516 if (j==s%ny+1) then 03517 deallocate(Eforc,D,deltheta,KKx,KKy,dphi3,k3,c,theta3,Gn,Abnd,w3,Ftemp) 03518 endif 03519 endif 03520 enddo 03521 deallocate(q) 03522 endif ! par%nonhspectrum==1 03523 03524 03525 end subroutine generate_secondorder 03526 03527 end module spectral_wave_bc_module
1.7.6.1
