|
XBeach
|
00001 module wave_boundary_update_module 00002 use wave_boundary_datastore 00003 ! 00004 implicit none 00005 save 00006 private 00007 public generate_wave_boundary_surfbeat 00008 ! 00009 type spectrum ! These are related to input spectra 00010 real*8,dimension(:,:),pointer :: S ! 2D variance density spectrum 00011 real*8,dimension(:),pointer :: f,ang ! 1D frequency and direction vectors for S 00012 real*8,dimension(:),pointer :: Sf ! S integrated over all directions 00013 integer :: nf,nang ! number of frequencies and angles 00014 real*8 :: df,dang ! frequency and angle step size 00015 real*8 :: hm0,fp,dir0,scoeff ! imposed significant wave height, peak frequency, 00016 ! main wave angle and spreading coefficient 00017 real*8 :: trep,dirm ! representative period and mean wave direction 00018 endtype spectrum 00019 type shortspectrum 00020 real*8,dimension(:),pointer :: Sf ! S integrated over all directions 00021 endtype shortspectrum 00022 type waveparamsnew ! These are related to the generated bc series 00023 real*8 :: h0 ! average water depth on offshore boundary 00024 integer :: K ! number of wave train components 00025 real*8 :: rtbc,dtbc ! duration and time step to be written to boundary condition file 00026 real*8 :: rtin,dtin ! duration and time step for the internal time axis, based on Fourier 00027 ! limitations and the number of wave train components (wp%K) 00028 real*8,dimension(:),pointer :: tin ! internal time axis 00029 real*8,dimension(:),pointer :: taperf,taperw ! internal taper function for flow and waves 00030 integer :: tslen ! internal time axis length (is always even, not odd) 00031 integer :: tslenbc ! time axis length for boundary condition file 00032 logical :: dtchanged ! quick check to see if dtin == dtbc (useful for interpolation purpose) 00033 real*8,dimension(:),pointer :: fgen,thetagen,phigen,kgen,wgen ! frequency, angle, phase, wave number, radian frequency 00034 ! of wavetrain components for boundary signal 00035 real*8 :: dfgen ! frequency grid size in the generated components 00036 type(shortspectrum),dimension(:),pointer :: vargen ! This is where the variance for each wave train at each spectrum location 00037 ! is stored 00038 type(shortspectrum),dimension(:),pointer :: vargenq ! This is where the variance for each wave train at each spectrum location 00039 ! is stored, which is not scaled and is used for the generation of bound waves 00040 real*8,dimension(:,:),pointer :: A ! Amplitude, per wave train component, per offshore grid point A(ny+1,K) 00041 real*8,dimension(:,:),pointer :: Sfinterp ! S integrated over all directions at frequency locations of fgen, 00042 ! per offshore grid point Sfinterp(ny+1,K) 00043 real*8,dimension(:,:),pointer :: Sfinterpq ! S integrated over all directions at frequency locations of fgen, 00044 ! per offshore grid point Sfinterpq(ny+1,K), uncorrected for generation of bound waves 00045 real*8,dimension(:),pointer :: Hm0interp ! Hm0 per offshore point, based on intergration of Sfinterp and used to scale 00046 ! final time series 00047 integer, dimension(:),pointer :: Findex ! Index of wave train component locations on frequency/Fourier axis 00048 integer, dimension(:),pointer :: WDindex ! Index of wave train component locations on wave directional bin axis 00049 double complex,dimension(:,:),pointer :: CompFn ! Fourier components of the wave trains 00050 character(1024) :: Efilename,qfilename,nhfilename 00051 real*8,dimension(:,:),pointer :: zsits ! time series of total surface elevation for nonhspectrum==1 00052 real*8,dimension(:,:),pointer :: uits ! time series of depth-averaged horizontal velocity nonhspectrum==1 00053 real*8,dimension(:,:),pointer :: wits ! time series of depth-averaged vertical velocity for nonhspectrum==1 ?? 00054 endtype waveparamsnew 00055 ! 00056 ! These parameters control a lot how the spectra are handled. They could be put in params.txt, 00057 ! but most users will want to keep these at their default values anyway 00058 integer,parameter :: nfint = 801 ! size of standard 2D spectrum in frequency dimension 00059 integer,parameter :: naint = 401 ! size of standard 2D spectrum in angular dimension 00060 integer,parameter :: Kmin = 200 ! minimum number of wave train components 00061 real*8,parameter :: wdmax = 5.d0 ! maximum depth*reliable angular wave frequency that can be resolved by 00062 ! nonhydrostatic wave model. All frequencies above this are removed 00063 ! from nonhspectrum generation 00064 ! Shortcut pointers to commonly used parameters 00065 integer :: nspectra,bccount,npb 00066 real*8 :: hb0 00067 ! Physical constants 00068 real*8,parameter :: par_pi = 4.d0*atan(1.d0) 00069 real*8,parameter :: par_g = 9.81d0 00070 complex(kind(0.0d0)),parameter :: par_compi = (0.0d0,1.0d0) 00071 contains 00072 00073 subroutine generate_wave_boundary_surfbeat(durationlength) 00074 #ifdef BUILDXBEACH 00075 use logging_module 00076 #endif 00077 use wave_boundary_datastore 00078 implicit none 00079 ! Input and output variables 00080 real*8,intent(out) :: durationlength 00081 ! 00082 ! 00083 ! Internal variables 00084 type(spectrum),dimension(:),allocatable :: specin,specinterp 00085 type(spectrum) :: combspec 00086 type(waveparamsnew) :: wp ! Most will be deallocated, but some variables useful to keep? 00087 integer :: iloc 00088 integer :: iostat 00089 real*8 :: spectrumendtimeold,fmax 00090 00091 ! Shortcuts, could be pointers? 00092 nspectra = waveSpectrumAdministration%nspectra 00093 bccount = waveSpectrumAdministration%bccount 00094 ! Offshore water depth, which is used in various computations in this module 00095 hb0 = waveBoundaryParameters%hboundary 00096 npb = waveBoundaryParameters%np 00097 ! 00098 ! 00099 ! Start Wave boundary condition time series generation 00100 if (waveSpectrumAdministration%repeatwbc) then 00101 ! Return wave boundary conditions that have already been computed in a previous 00102 ! call to this subroutine. Modify time axis to reflect shift in time since the 00103 ! previous call 00104 waveBoundaryTimeSeries%tbc = min(waveBoundaryTimeSeries%tbc, & 00105 huge(0.d0)-waveBoundaryAdministration%startComputeNewSeries) + & 00106 waveBoundaryAdministration%startComputeNewSeries 00107 00108 else 00109 #ifdef BUILDXBEACH 00110 call writelog('l','','--------------------------------') 00111 call writelog('ls','','Calculating spectral wave boundary conditions ') 00112 call writelog('l','','--------------------------------') 00113 #endif 00114 ! allocate temporary input storage 00115 if (.not. allocated(specin)) then 00116 allocate(specin(nspectra)) 00117 allocate(specinterp(nspectra)) 00118 endif 00119 00120 ! Read through input spectra files 00121 fmax = 1.d0 ! assume 1Hz as maximum frequency. Increase in loop below if needed. 00122 do iloc = 1,nspectra 00123 #ifdef BUILDXBEACH 00124 call writelog('sl','(a,i0)','Reading spectrum at location ',iloc) 00125 #endif 00126 ! Read input file 00127 call read_spectrum_input(wp,bcfiles(iloc),specin(iloc)) 00128 fmax = max(fmax,maxval(specin(iloc)%f)) 00129 enddo 00130 do iloc = 1,nspectra 00131 #ifdef BUILDXBEACH 00132 call writelog('sl','(a,i0)','Interpreting spectrum at location ',iloc) 00133 #endif 00134 ! Interpolate input 2D spectrum to standard 2D spectrum 00135 call interpolate_spectrum(specin(iloc),specinterp(iloc),par,fmax) 00136 #ifdef BUILDXBEACH 00137 call writelog('sl','','Values calculated from interpolated spectrum:') 00138 call writelog('sl','(a,f0.2,a)','Hm0 = ',specinterp(iloc)%hm0,' m') 00139 call writelog('sl','(a,f0.2,a)','Trep = ',specinterp(iloc)%trep,' s') 00140 call writelog('sl','(a,f0.2,a)','Mean dir = ',specinterp(iloc)%dirm,' degN') 00141 #endif 00142 enddo 00143 00144 ! Determine whether all the spectra are to be reused, which implies that the global repeatwbc should be 00145 ! set to true (no further computations required in future calls) 00146 call set_repeatwbc 00147 ! 00148 ! calculate the mean combined spectra (used for combined Trep, determination of wave components, etc.) 00149 ! now still uses simple averaging, but could be improved to use weighting for distance etc. 00150 call generate_combined_spectrum(specinterp,combspec) 00151 ! Store these data in wave administration for exchange with outer 00152 ! XBeach or D3D/FM models 00153 waveSpectrumAdministration%Hbc = combspec%hm0 00154 waveSpectrumAdministration%Tbc = combspec%trep 00155 waveSpectrumAdministration%Dbc = combspec%dirm 00156 #ifdef BUILDXBEACH 00157 call writelog('sl','(a,f0.2,a)','Overall Trep from all spectra calculated: ',& 00158 waveSpectrumAdministration%Tbc,' s') 00159 #endif 00160 ! Wave trains that are used by XBeach. The number of wave trains, their frequencies and directions 00161 ! are based on the combined spectra of all the locations to ensure all wave conditions are 00162 ! represented in the XBeach model 00163 call generate_wavetrain_components(combspec,wp) 00164 00165 ! We can now apply a correction to the wave train components if necessary. This section can be 00166 ! improved later 00167 if (nspectra==1 .and. specin(1)%scoeff>1000.d0) then 00168 ! this can be used both for Jonswap and vardens input 00169 wp%thetagen = mod(specin(1)%dir0,2*par_pi) 00170 endif 00171 00172 ! Set up time axis, including the time axis for output to boundary condition files and an 00173 ! internal time axis, which may differ in length to the output time axis 00174 call generate_wave_time_axis(wp) 00175 00176 ! Determine the variance for each wave train component, at every spectrum location point 00177 call generate_wave_train_variance(wp,specinterp) 00178 00179 ! Determine the amplitude of each wave train component, at every point along the 00180 ! offshore boundary 00181 call generate_wave_train_properties_per_offshore_point(wp,s) 00182 00183 ! Generate Fourier components for all wave train component, at every point along 00184 ! the offshore boundary 00185 call generate_wave_train_Fourier(wp) 00186 00187 ! Time series of short wave energy or surface elevation 00188 ! Distribute all wave train components among the wave direction bins. Also rearrage 00189 ! the randomly drawn wave directions to match the centres of the wave bins if the 00190 ! user-defined nspr is set on. 00191 call distribute_wave_train_directions(wp,waveBoundaryParameters%nspr) 00192 if (.not.waveBoundaryParameters%nonhspectrum) then 00193 ! Calculate the wave energy envelope per offshore grid point and write to output file 00194 call generate_ebcf(wp) 00195 else 00196 ! Generate time series of surface elevation, horizontal velocity and vertical velocity 00197 call generate_swts(wp,s) 00198 endif 00199 00200 ! Calculate the bound long wave from the wave train components and write to output file 00201 call generate_qbcf(wp) 00202 00203 ! Write non-hydrostatic time series of combined short and long waves if necessary 00204 if (waveBoundaryParameters%nonhspectrum) then 00205 call generate_nhtimeseries_file(wp) 00206 endif 00207 00208 ! Deallocate a lot of memory 00209 deallocate(specin,specinterp) 00210 deallocate(wp%tin,wp%taperf,wp%taperw) 00211 deallocate(wp%fgen,wp%thetagen,wp%phigen,wp%kgen,wp%wgen) 00212 deallocate(wp%vargen) 00213 deallocate(wp%vargenq) 00214 deallocate(wp%Sfinterp) 00215 deallocate(wp%Sfinterpq) 00216 deallocate(wp%Hm0interp) 00217 deallocate(wp%A) 00218 deallocate(wp%Findex) 00219 deallocate(wp%CompFn) 00220 if (waveBoundaryParameters%nonhspectrum) then 00221 deallocate(wp%zsits) 00222 deallocate(wp%uits) 00223 else 00224 deallocate(wp%WDindex) 00225 endif 00226 ! Send message to screen and log 00227 #if BUILDXBEACH 00228 call writelog('l','','--------------------------------') 00229 call writelog('ls','','Spectral wave boundary conditions complete ') 00230 call writelog('l','','--------------------------------') 00231 #endif 00232 endif ! repeatwbc 00233 ! 00234 ! Return the duration length of the current boundary conditions 00235 durantionlength = wp%rtbc 00236 ! 00237 ! Update the number of boundary conditions generated by this module 00238 bccount = bccount+1 00239 00240 endif 00241 00242 endsubroutine generate_wave_boundary_surfbeat 00243 00244 ! -------------------------------------------------------------- 00245 ! ---------------- Read input spectra files -------------------- 00246 ! -------------------------------------------------------------- 00247 subroutine read_spectrum_input(wp,fn,specin) 00248 00249 use filefunctions 00250 use logging_module 00251 00252 implicit none 00253 ! Interface 00254 type(waveparamsnew),intent(inout) :: wp 00255 type(filenames),intent(inout) :: fn 00256 type(spectrum),intent(inout) :: specin 00257 ! internal 00258 integer :: fid 00259 character(8) :: testline 00260 character(1024) :: readfile 00261 logical :: filelist 00262 integer :: i,ier 00263 00264 ! check the first line of the boundary condition file for FILELIST keyword 00265 fid = create_new_fid() 00266 open(fid,file=fn%fname,status='old',form='formatted') 00267 read(fid,*,iostat=ier)testline 00268 if (ier .ne. 0) then 00269 call report_file_read_error(fn%fname) 00270 endif 00271 if (trim(testline)=='FILELIST') then 00272 filelist = .true. 00273 ! move listline off its default position of zero to the first row 00274 if (fn%listline==0) then 00275 fn%listline = 1 00276 endif 00277 fn%repeat = .false. 00278 else 00279 filelist = .false. 00280 if (par%instat /= INSTAT_JONS_TABLE) then 00281 fn%repeat = .true. 00282 else 00283 fn%repeat = .false. 00284 endif 00285 endif 00286 close(fid) 00287 ! If file has list of spectra, read through the lines to find the correct 00288 ! spectrum file name and rtbc and dtbc. Else the filename and rtbc, dtbc are in 00289 ! params.txt 00290 if (filelist) then 00291 fid = create_new_fid() 00292 open(fid,file=fn%fname,status='old',form='formatted') 00293 do i=1,fn%listline 00294 read(fid,*)testline ! old stuff, not needed anymore 00295 enddo 00296 read(fid,*,iostat=ier)wp%rtbc,wp%dtbc,readfile ! new boundary condition 00297 if (ier .ne. 0) then 00298 call report_file_read_error(fn%fname) 00299 endif 00300 ! we have to adjust this to morphological time, as done in params.txt 00301 if (par%morfacopt==1) then 00302 wp%rtbc = wp%rtbc/max(par%morfac,1.d0) 00303 endif 00304 fn%listline = fn%listline + 1 ! move one on from the last time we opened this file 00305 close(fid) 00306 else 00307 wp%rtbc = par%rt ! already set to morphological time 00308 wp%dtbc = par%dtbc 00309 readfile = fn%fname 00310 endif 00311 00312 ! based on the value of instat, we need to read either Jonswap, Swan or vardens files 00313 ! note: jons_table is also handeled by read_jonswap_file subroutine 00314 !select case (par%instat(1:4)) 00315 select case(par%instat) 00316 !case ('jons') 00317 case (INSTAT_JONS, INSTAT_JONS_TABLE) 00318 ! wp type sent in to receive rtbc and dtbc from jons_table file 00319 ! fn%listline sent in to find correct row in jons_table file 00320 ! pfff.. 00321 call read_jonswap_file(par,wp,readfile,fn%listline,specin) 00322 !case ('swan') 00323 case (INSTAT_SWAN) 00324 call read_swan_file(par,readfile,specin) 00325 !case ('vard') 00326 case (INSTAT_VARDENS) 00327 call read_vardens_file(par,readfile,specin) 00328 endselect 00329 00330 endsubroutine read_spectrum_input 00331 00332 ! -------------------------------------------------------------- 00333 ! ------------------- Read JONSWAP files ----------------------- 00334 ! -------------------------------------------------------------- 00335 subroutine read_jonswap_file(par,wp,readfile,listline,specin) 00336 00337 use readkey_module 00338 use params 00339 use logging_module 00340 use filefunctions 00341 00342 IMPLICIT NONE 00343 00344 ! Input / output variables 00345 type(parameters), INTENT(IN) :: par 00346 type(waveparamsnew),intent(inout) :: wp 00347 character(len=*), intent(IN) :: readfile 00348 integer, intent(INOUT) :: listline 00349 type(spectrum),intent(inout) :: specin 00350 00351 ! Internal variables 00352 integer :: i,ii,ier,ip,ind 00353 integer :: nmodal 00354 integer :: fid 00355 integer :: forcepartition 00356 integer,dimension(2) :: indvec 00357 real*8,dimension(:),allocatable :: x, y, Dd, tempdir 00358 real*8,dimension(:),allocatable :: Hm0,fp,gam,mainang,scoeff 00359 real*8 :: dfj, fnyq 00360 real*8 :: Tp 00361 character(len=80) :: dummystring 00362 type(spectrum),dimension(:),allocatable :: multinomalspec,scaledspec 00363 logical :: cont 00364 real*8,dimension(:),allocatable :: scalefac1,scalefac2,tempmax,avgscale 00365 real*8,dimension(:),allocatable :: oldvariance,newvariance 00366 real*8 :: newconv,oldconv 00367 00368 ! First part: read JONSWAP parameter data 00369 00370 ! Check whether spectrum characteristics or table should be used 00371 if (par%instat /= INSTAT_JONS_TABLE) then 00372 ! Use spectrum characteristics 00373 call writelog('sl','','waveparams: Reading from ',trim(readfile),' ...') 00374 ! 00375 ! First read if the spectrum is multinodal, and how many peaks there should be 00376 ! 00377 nmodal = readkey_int (readfile, 'nmodal', 1, 1, 4, bcast=.false.) 00378 if (nmodal<1) then 00379 call writelog('lswe','','Error: number of spectral partions may not be less than 1 in ',trim(readfile)) 00380 call halt_program 00381 endif 00382 ! 00383 ! Allocate space for all spectral parameters 00384 ! 00385 allocate(Hm0(nmodal)) 00386 allocate(fp(nmodal)) 00387 allocate(gam(nmodal)) 00388 allocate(mainang(nmodal)) 00389 allocate(scoeff(nmodal)) 00390 ! 00391 ! Read the spectral parameters for all spectrum components 00392 ! 00393 ! Wave height (required) 00394 Hm0 = readkey_dblvec(readfile, 'Hm0',nmodal,nmodal, 0.0d0, 0.0d0, 5.0d0, bcast=.false.,required=.true. ) 00395 ! 00396 ! Wave period (required) 00397 ! allow both Tp and fp specification to bring in line with params.txt 00398 if (isSetParameter(readfile,'Tp',bcast=.false.) .and. .not. isSetParameter(readfile,'fp',bcast=.false.)) then 00399 fp = 1.d0/readkey_dblvec(readfile, 'Tp',nmodal,nmodal, 12.5d0, 2.5d0, 20.0d0, bcast=.false.) 00400 elseif (isSetParameter(readfile,'fp',bcast=.false.) .and. .not. isSetParameter(readfile,'Tp',bcast=.false.)) then 00401 fp = readkey_dblvec(readfile, 'fp',nmodal,nmodal, 0.08d0,0.0625d0, 0.4d0, bcast=.false.) 00402 elseif (.not. isSetParameter(readfile,'fp',bcast=.false.) .and. .not. isSetParameter(readfile,'Tp',bcast=.false.)) then 00403 call writelog('lswe','','Error: missing required value for parameter ''Tp'' or ''fp'' in ',trim(readfile)) 00404 call halt_program 00405 else 00406 fp = 1.d0/readkey_dblvec(readfile, 'Tp',nmodal,nmodal, 12.5d0, 2.5d0, 20.0d0, bcast=.false.) 00407 call writelog('lsw','','Warning: selecting to read peak period (Tp) instead of frequency (fp) in ',trim(readfile)) 00408 endif 00409 ! 00410 ! Wave spreading in frequency domain (peakedness) 00411 ! 00412 gam = readkey_dblvec(readfile, 'gammajsp',nmodal,nmodal, 3.3d0, 1.0d0, 5.0d0, bcast=.false.) 00413 ! 00414 ! Wave spreading in directional domain 00415 ! 00416 scoeff = readkey_dblvec(readfile, 's',nmodal,nmodal, 10.0d0, 1.0d0, 1000.0d0, bcast=.false.) 00417 ! 00418 ! Main wave direction 00419 ! allow both mainang and dir0 specification to bring in line with params.txt 00420 if (isSetParameter(readfile,'mainang',bcast=.false.) .and. & 00421 .not. isSetParameter(readfile,'dir0',bcast=.false.)) then 00422 mainang = readkey_dblvec(readfile, 'mainang',nmodal,nmodal, & 00423 270.0d0, 0.0d0, 360.0d0, bcast=.false.) 00424 elseif (isSetParameter(readfile,'dir0',bcast=.false.) .and. & 00425 .not. isSetParameter(readfile,'mainang',bcast=.false.)) then 00426 mainang = readkey_dblvec(readfile, 'dir0',nmodal,nmodal, & 00427 270.0d0, 0.0d0, 360.0d0, bcast=.false.) 00428 elseif (.not. isSetParameter(readfile,'dir0',bcast=.false.) .and. & 00429 .not. isSetParameter(readfile,'mainang',bcast=.false.)) then 00430 mainang = 270.d0 00431 else 00432 mainang = readkey_dblvec(readfile, 'mainang',nmodal,nmodal, 270.0d0, 0.0d0, 360.0d0, bcast=.false.) 00433 call writelog('lsw','','Warning: selecting to read ''mainang'' instead of ''dir0'' in ',trim(readfile)) 00434 endif 00435 ! 00436 ! Nyquist parameters used only in this subroutine 00437 ! are not read individually for each spectrum partition 00438 if (par%oldnyq==1) then 00439 fnyq = readkey_dbl (readfile, 'fnyq', 0.3d0, 0.2d0, 1.0d0, bcast=.false. ) 00440 else 00441 fnyq = readkey_dbl (readfile, 'fnyq',max(0.3d0,3.d0*maxval(fp)), 0.2d0, 1.0d0, bcast=.false. ) 00442 endif 00443 dfj = readkey_dbl (readfile, 'dfj', fnyq/200, fnyq/1000, fnyq/20, bcast=.false. ) 00444 ! 00445 ! Finally check if XBeach should accept even the most stupid partioning (sets error level to warning 00446 ! level when computing partition overlap 00447 if (nmodal>1) then 00448 forcepartition = readkey_int (readfile, 'forcepartition', 0, 0, 1, bcast=.false.) 00449 endif 00450 ! check for other strange values in this file 00451 call readkey(readfile,'checkparams',dummystring) 00452 else 00453 nmodal = 1 00454 allocate(Hm0(nmodal)) 00455 allocate(fp(nmodal)) 00456 allocate(gam(nmodal)) 00457 allocate(mainang(nmodal)) 00458 allocate(scoeff(nmodal)) 00459 ! Use spectrum table 00460 fid = create_new_fid() 00461 call writelog('sl','','waveparams: Reading from table ',trim(readfile),' ...') 00462 open(fid,file=readfile,status='old',form='formatted') 00463 ! read junk up to the correct line in the file 00464 do i=1,listline 00465 read(fid,*,iostat=ier)dummystring 00466 if (ier .ne. 0) then 00467 call report_file_read_error(readfile) 00468 endif 00469 enddo 00470 read(fid,*,iostat=ier)Hm0(1),Tp,mainang(1),gam(1),scoeff(1),wp%rtbc,wp%dtbc 00471 if (ier .ne. 0) then 00472 call writelog('lswe','','Error reading file ',trim(readfile)) 00473 close(fid) 00474 call halt_program 00475 endif 00476 ! move the line pointer in the file 00477 listline = listline+1 00478 ! convert to morphological time 00479 if (par%morfacopt==1) then 00480 wp%rtbc = wp%rtbc/max(par%morfac,1.d0) 00481 endif 00482 fp(1)=1.d0/Tp 00483 fnyq=3.d0*fp(1) 00484 dfj=fp(1)/50 00485 close(fid) 00486 endif 00487 ! 00488 ! Second part: generate 2D spectrum from input parameters 00489 ! 00490 ! Define number of frequency bins by defining an array of the necessary length 00491 ! using the Nyquist frequency and frequency step size 00492 specin%nf = ceiling((fnyq-dfj)/dfj) 00493 specin%df = dfj 00494 ! 00495 ! Define array with actual eqidistant frequency bins 00496 allocate(specin%f(specin%nf)) 00497 do i=1,specin%nf 00498 specin%f(i)=i*dfj 00499 enddo 00500 ! 00501 ! we need a normalised frequency and variance vector for JONSWAP generation 00502 allocate(x(size(specin%f))) 00503 allocate(y(size(specin%f))) 00504 ! 00505 ! Define 200 directions relative to main angle running from 0 to 2*pi 00506 ! we also need a temporary vector for direction distribution 00507 specin%nang = naint 00508 allocate(tempdir(specin%nang)) 00509 allocate(Dd(specin%nang)) 00510 allocate(specin%ang(specin%nang)) 00511 specin%dang = 2*par_pi/(naint-1) 00512 do i=1,specin%nang 00513 specin%ang(i)=(i-1)*specin%dang 00514 enddo 00515 ! 00516 ! Generate density spectrum for each spectrum partition 00517 ! 00518 allocate(multinomalspec(nmodal)) 00519 ! 00520 do ip=1,nmodal 00521 ! relative frequenct vector 00522 x=specin%f/fp(ip) 00523 ! Calculate unscaled and non-directional JONSWAP spectrum using 00524 ! peak-enhancement factor and pre-determined frequency bins 00525 call jonswapgk(x,gam(ip),y) 00526 ! Determine scaled and non-directional JONSWAP spectrum using the JONSWAP 00527 ! characteristics 00528 y=(Hm0(ip)/(4.d0*sqrt(sum(y)*dfj)))**2*y 00529 ! Convert main angle from degrees to radians and from nautical convention to 00530 ! internal grid 00531 mainang(ip)=(1.5d0*par_pi)-mainang(ip)*par_pi/180 00532 ! Make sure the main angle is defined between 0 and 2*pi 00533 do while (mainang(ip)>2*par_pi .or. mainang(ip)<0.d0) !Robert en Ap 00534 if (mainang(ip)>2*par_pi) then 00535 mainang(ip)=mainang(ip)-2*par_pi 00536 elseif (mainang(ip)<0.d0) then 00537 mainang(ip)=mainang(ip)+2*par_pi 00538 endif 00539 enddo 00540 ! Convert 200 directions relative to main angle to directions relative to 00541 ! internal grid 00542 ! Bas: apparently division by 2 for cosine law happens already here 00543 tempdir = (specin%ang-mainang(ip))/2 00544 ! Make sure all directions around the main angle are defined between 0 and 2*pi 00545 do while (any(tempdir>2*par_pi) .or. any(tempdir<0.d0)) 00546 where (tempdir>2*par_pi) 00547 tempdir=tempdir-2*par_pi 00548 elsewhere (tempdir<0.d0) 00549 tempdir=tempdir+2*par_pi 00550 endwhere 00551 enddo 00552 ! Calculate directional spreading based on cosine law 00553 Dd = dcos(tempdir)**(2*nint(scoeff(ip))) ! Robert: apparently nint is needed here, else MATH error 00554 ! Scale directional spreading to have a surface of unity by dividing by its 00555 ! own surface 00556 Dd = Dd / (sum(Dd)*specin%dang) 00557 ! Define two-dimensional variance density spectrum array and distribute 00558 ! variance density for each frequency over directional bins 00559 allocate(multinomalspec(ip)%S(specin%nf,specin%nang)) 00560 do i=1,specin%nang 00561 do ii=1,specin%nf 00562 multinomalspec(ip)%S(ii,i)=y(ii)*Dd(i) 00563 end do 00564 end do 00565 enddo ! 1,nmodal 00566 do ip=1,nmodal 00567 write(*,*)'Hm0 = ',4*sqrt(sum(multinomalspec(ip)%S)*specin%dang*specin%df) 00568 enddo 00569 ! 00570 ! Combine spectrum partitions so that the total spectrum is correct 00571 ! 00572 if (nmodal==1) then 00573 ! Set all the useful parameters and arrays 00574 specin%Hm0 = Hm0(1) 00575 specin%fp = fp(1) 00576 specin%dir0 = mainang(1) 00577 specin%scoeff = scoeff(1) 00578 allocate(specin%S(specin%nf,specin%nang)) 00579 specin%S = multinomalspec(1)%S 00580 else 00581 specin%Hm0 = sqrt(sum(Hm0**2)) ! total wave height 00582 ind = maxval(maxloc(Hm0)) ! spectrum that is largest 00583 specin%fp = fp(ind) ! not really used in further calculation 00584 specin%dir0 = mainang(ind) ! again not really used in further calculation 00585 ! if all scoeff>1000 then all waves should be in the same direction exactly 00586 if (all(scoeff>=1024.d0) .and. all(mainang==mainang(ind))) then 00587 specin%scoeff = 1024.d0 00588 else 00589 specin%scoeff = min(scoeff(ind),999.d0) 00590 endif 00591 ! Now we have to loop over all partitioned spectra. Where two or more spectra 00592 ! overlap, only the largest is counted, and all others are set to zero. Afterwards 00593 ! all spectra are scaled so that the total energy is maintained. Since the scaling 00594 ! affects where spectra overlap, this loop is repeated until only minor changes 00595 ! in the scaling occur. Warnings or errors are given if the spectrum is scaled too 00596 ! much, i.e. spectra overlap too much. 00597 ! 00598 ! Allocate space 00599 allocate(scaledspec(nmodal)) ! used to store scaled density spectra 00600 do ip=1,nmodal 00601 allocate(scaledspec(ip)%S(specin%nf,specin%nang)) 00602 scaledspec(ip)%S = multinomalspec(ip)%S 00603 enddo 00604 allocate(scalefac1(nmodal)) ! this is the scaling factor required to maintain Hm0 00605 ! including that parts of the spectrum are set to zero 00606 scalefac1 = 1.d0 00607 allocate(scalefac2(nmodal)) ! this is the scaling factor required to maintain Hm0 00608 ! in the previous iteration 00609 scalefac2 = 1.d0 00610 allocate(avgscale(nmodal)) 00611 avgscale = 1.d0 00612 ! these are convergence criteria 00613 newconv = 0.d0 00614 oldconv = huge(0.d0) 00615 cont = .true. 00616 allocate(tempmax(nmodal)) ! used to store maximum value in f,theta space 00617 allocate(oldvariance(nmodal)) ! used to store the sum of variance in the original partitions 00618 do ip=1,nmodal 00619 oldvariance(ip) = sum(multinomalspec(ip)%S) 00620 enddo 00621 allocate(newvariance(nmodal)) ! used to store the sum of variance in the new partitions 00622 ! 00623 ! Start convergence loop 00624 do while (cont) 00625 avgscale = (avgscale+scalefac1)/2 00626 ! First scale the spectra 00627 do ip=1,nmodal 00628 ! scale the spectrum using less than half the additional scale factor, this to 00629 ! ensure the method does not become unstable 00630 ! scaledspec(ip)%S = multinomalspec(ip)%S*(0.51d0+scalefac1(ip)*0.49d0) 00631 ! scaledspec(ip)%S = multinomalspec(ip)%S*(scalefac1(ip)+scalefac2(ip))/2 00632 scaledspec(ip)%S = multinomalspec(ip)%S*avgscale(ip) 00633 00634 enddo 00635 do i=1,specin%nang 00636 do ii=1,specin%nf 00637 ! vector of variance densities at this point in f,theta space 00638 do ip=1,nmodal 00639 tempmax(ip) = scaledspec(ip)%S(ii,i) 00640 end do 00641 ! All spectra other than the one that is greatest in this point 00642 ! is set to zero. In case multiple spectra are equal largest, 00643 ! the first spectrum in the list is chosen (minval(maxloc)) 00644 ind = minval(maxloc(tempmax)) 00645 do ip=1,nmodal 00646 if (ip/=ind) then 00647 scaledspec(ip)%S(ii,i) = 0.d0 00648 endif 00649 enddo 00650 enddo 00651 end do 00652 ! 00653 ! Now rescale adjusted partition spectra so that they match the incident Hm0 00654 scalefac2 = scalefac1 ! keep previous results 00655 do ip=1,nmodal 00656 newvariance(ip) = sum(scaledspec(ip)%S) 00657 if (newvariance(ip)>0.01d0*oldvariance(ip)) then 00658 scalefac1(ip) = oldvariance(ip)/newvariance(ip) ! want to maximise to a factor 2 00659 ! else can generate rediculous results 00660 else 00661 scalefac1(ip) = 0.d0 ! completely remove this spectrum 00662 endif 00663 enddo 00664 ! 00665 ! check convergence criteria (if error is increasing, we have passed best fit) 00666 newconv = maxval(abs(scalefac2-scalefac1)/scalefac1) 00667 if (newconv<0.0001d0 .or. abs(newconv-oldconv)<0.0001d0) then 00668 cont = .false. 00669 endif 00670 oldconv = newconv 00671 end do 00672 ! Ensure full energy conservation by scaling now according to full scalefac, not just the 00673 ! half used in the iteration loop 00674 do ip=1,nmodal 00675 scaledspec(ip)%S = multinomalspec(ip)%S*scalefac1(ip) 00676 enddo 00677 ! 00678 ! Now check the total scaling that has taken place across the spectrum, also accounting 00679 ! for the fact that some parts of the spectra are set to zero. This differs from the 00680 ! scalefac1 and scalefac2 vectors in the loop that only adjust the part of the spectrum 00681 ! that is not zero. 00682 do ip=1,nmodal 00683 write(*,*)'Hm0m = ',4*sqrt(sum(scaledspec(ip)%S)*specin%dang*specin%df) 00684 enddo 00685 do ip=1,nmodal 00686 indvec = maxloc(scaledspec(ip)%S) 00687 scalefac1(ip) = scaledspec(ip)%S(indvec(1),indvec(2)) / & 00688 multinomalspec(ip)%S(indvec(1),indvec(2))-1.d0 00689 enddo 00690 ! 00691 ! Warning and/or error criteria here if spectra overlap each other too much 00692 do ip=1,nmodal 00693 if(scalefac1(ip)>0.5d0) then 00694 if (forcepartition==1) then 00695 call writelog('lsw','(a,f0.0,a,i0,a,a,a)', & 00696 'Warning: ',scalefac1(ip)*100,'% of energy in spectrum partition ''',ip, & 00697 ''' in ', trim(readfile),' is overlapped by other partitions') 00698 call writelog('lsw','',' Check spectral partitioning in ',trim(readfile)) 00699 else 00700 call writelog('lswe','(a,f0.0,a,i0,a,a,a)', & 00701 'Error: ',scalefac1(ip)*100,'% of energy in spectrum partition ''',ip, & 00702 ''' in ', trim(readfile),' is overlapped by other partitions') 00703 call writelog('lswe','','This spectrum overlaps too much with another spectrum partition.', & 00704 ' Check spectral partitioning in ',trim(readfile)) 00705 call writelog('lswe','','If the partitioning should be carried out regardles of energy loss, ', & 00706 'set ''forcepartition = 1'' in ',trim(readfile)) 00707 call halt_program 00708 endif 00709 elseif (scalefac1(ip)>0.2d0 .and. scalefac1(ip)<=0.5d0) then 00710 call writelog('lsw','(a,f0.0,a,i0,a,a,a)', & 00711 'Warning: ',scalefac1(ip)*100,'% of energy in spectrum partition ''',ip, & 00712 ''' in ', trim(readfile),' is overlapped by other partitions') 00713 call writelog('lsw','',' Check spectral partitioning in ',trim(readfile)) 00714 elseif (scalefac1(ip)<0.d0) then 00715 call writelog('lsw','(a,i0,a,a,a)','Warning: spectrum partition ''',ip,''' in ',trim(readfile), & 00716 ' has been removed') 00717 call writelog('lsw','','This spectrum is entirely overlapped by another spectrum partition.',& 00718 ' Check spectral partitioning in ',trim(readfile)) 00719 endif 00720 enddo 00721 ! 00722 ! Now set total spectrum 00723 allocate(specin%S(specin%nf,specin%nang)) 00724 specin%S = 0.d0 00725 do ip=1,nmodal 00726 specin%S = specin%S+scaledspec(ip)%S 00727 enddo 00728 00729 deallocate(scaledspec) 00730 deallocate(tempmax) 00731 deallocate(scalefac1,scalefac2) 00732 deallocate(newvariance,oldvariance) 00733 endif 00734 00735 ! We need frequency spectrum to ensure Sf remains correct between interpoloation 00736 ! routines 00737 allocate(specin%Sf(specin%nf)) 00738 specin%Sf=0.d0 00739 do i=1,specin%nf 00740 do ii=1,specin%nang 00741 if (ii==1) then 00742 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(2)-specin%ang(1)) 00743 elseif (ii==specin%nang) then 00744 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii)-specin%ang(ii-1)) 00745 else 00746 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii+1)-specin%ang(ii-1))/2 00747 endif 00748 enddo 00749 enddo 00750 00751 deallocate (Dd) 00752 deallocate(Hm0,fp,gam,mainang,scoeff) 00753 deallocate (x,y) 00754 deallocate(tempdir) 00755 ! TODO dereference pointers first.... 00756 do ip=1,nmodal 00757 deallocate(multinomalspec(ip)%S) 00758 end do 00759 deallocate(multinomalspec) 00760 00761 00762 end subroutine read_jonswap_file 00763 00764 00765 00766 ! ----------------------------------------------------------- 00767 ! --------- JONSWAP unscaled JONSWAP spectrum -------------- 00768 ! -------------(used by read_jonswap_files)------------------ 00769 subroutine jonswapgk(x,gam,y) 00770 00771 IMPLICIT NONE 00772 ! Required input: - x : nondimensional frequency, divided by the peak frequency 00773 ! - gam : peak enhancement factor, optional parameter (DEFAULT 3.3) 00774 ! - y is output : nondimensional relative spectral density, equal to one at the peak 00775 00776 real*8, INTENT(IN) :: gam 00777 real*8,dimension(:), INTENT(IN) :: x 00778 real*8,dimension(:), INTENT(INOUT) :: y 00779 00780 ! Internal variables 00781 real*8,dimension(size(x)) :: xa, sigma, fac1, fac2, fac3, temp 00782 00783 xa=abs(x) 00784 00785 where (xa==0) 00786 xa=1e-20 00787 end where 00788 00789 sigma=xa 00790 00791 where (sigma<1.) 00792 sigma=0.07 00793 end where 00794 00795 where (sigma>=1.) 00796 sigma=0.09 00797 end where 00798 00799 temp=0*xa+1 00800 00801 fac1=xa**(-5) 00802 fac2=exp(-1.25*(xa**(-4))) 00803 fac3=(gam*temp)**(exp(-((xa-1)**2)/(2*(sigma**2)))) 00804 00805 y=fac1*fac2*fac3 00806 y=y/maxval(y) 00807 00808 return 00809 00810 end subroutine jonswapgk 00811 00812 00813 ! -------------------------------------------------------------- 00814 ! ----------------------Read SWAN files ------------------------ 00815 ! -------------------------------------------------------------- 00816 subroutine read_swan_file(par,readfile,specin) 00817 00818 use params 00819 use logging_module 00820 use filefunctions 00821 use xmpi_module, only: Halt_Program 00822 use math_tools, only: flipv,flipa 00823 00824 IMPLICIT NONE 00825 00826 ! Input / output variables 00827 type(parameters), INTENT(IN) :: par 00828 character(len=*), intent(IN) :: readfile 00829 type(spectrum),intent(inout) :: specin 00830 00831 ! Internal variables 00832 character(6) :: rtext 00833 real*8 :: factor,exc 00834 integer :: fid,switch 00835 integer :: i,ii,ier,ier2,ier3 00836 logical :: flipped 00837 integer :: nt,Ashift 00838 real*8, dimension(:),allocatable :: temp 00839 real*8, dimension(:,:),allocatable :: tempA 00840 00841 flipped =.false. 00842 switch = 0 00843 00844 call writelog('sl','','Reading from SWAN file ',trim(readfile),' ...') 00845 fid = create_new_fid() 00846 open(fid,file=readfile,form='formatted',status='old') 00847 00848 ! Read file until RFREQ or AFREQ is found 00849 do while (switch==0) 00850 read(fid,'(a)',iostat=ier)rtext 00851 if (ier .ne. 0) then 00852 call report_file_read_error(readfile) 00853 endif 00854 if (rtext == 'RFREQ ') then 00855 switch = 1 00856 elseif (rtext == 'AFREQ ') then 00857 switch = 2 00858 end if 00859 end do 00860 00861 ! Read nfreq and f 00862 ! Note f is not monotonically increasing in most simulations 00863 read(fid,*,iostat=ier)specin%nf 00864 if (ier .ne. 0) then 00865 call report_file_read_error(readfile) 00866 endif 00867 allocate(specin%f(specin%nf)) 00868 do i=1,specin%nf 00869 read(fid,*,iostat=ier)specin%f(i) 00870 if (ier .ne. 0) then 00871 call report_file_read_error(readfile) 00872 endif 00873 end do 00874 00875 ! Convert to absolute frequencies: 00876 ! STILL TO BE DONE 00877 if (switch == 1) then 00878 specin%f = specin%f 00879 else 00880 specin%f = specin%f 00881 end if 00882 00883 ! Read CDIR or NDIR 00884 read(fid,'(a)',iostat=ier)rtext 00885 if (ier .ne. 0) then 00886 call report_file_read_error(readfile) 00887 endif 00888 if (rtext == 'NDIR ') then 00889 switch = 1 00890 elseif (rtext == 'CDIR ') then 00891 switch = 2 00892 else 00893 call writelog('ewls','', 'SWAN directional bins keyword not found') 00894 call halt_program 00895 endif 00896 00897 ! Read ndir, theta 00898 read(fid,*,iostat=ier)specin%nang 00899 if (ier .ne. 0) then 00900 call report_file_read_error(readfile) 00901 endif 00902 allocate(specin%ang(specin%nang)) 00903 do i=1,specin%nang 00904 read(fid,*,iostat=ier)specin%ang(i) 00905 if (ier .ne. 0) then 00906 call report_file_read_error(readfile) 00907 endif 00908 end do 00909 00910 ! Convert angles to cartesian degrees relative to East 00911 if (switch == 1) then 00912 ! nautical to cartesian East 00913 specin%ang = 270.d0-specin%ang 00914 else 00915 ! cartesian to cartesian East 00916 specin%ang = specin%ang-par%dthetaS_XB 00917 ! dthetaS_XB is the (counter-clockwise) angle in the degrees to rotate from the x-axis in SWAN to the 00918 ! x-axis pointing East 00919 end if 00920 00921 ! Ensure angles are increasing instead of decreasing 00922 if (specin%ang(2)<specin%ang(1)) then 00923 call flipv(specin%ang,size(specin%ang)) 00924 flipped=.true. 00925 end if 00926 00927 nt = 0 00928 Ashift = 0 00929 ! Make sure that all angles are in range of 0 to 360 degrees 00930 if(minval(specin%ang)<0.d0)then 00931 allocate (temp(specin%nang)) 00932 Ashift=-1 00933 temp=0.d0 00934 do i=1,specin%nang 00935 if (specin%ang(i)<0.d0) then 00936 specin%ang(i)=specin%ang(i)+360.0d0 00937 nt = nt+1 00938 endif 00939 enddo 00940 temp(1:specin%nang-nt)=specin%ang(nt+1:specin%nang) 00941 temp(specin%nang-nt+1:specin%nang)=specin%ang(1:nt) 00942 specin%ang=temp 00943 deallocate(temp) 00944 elseif(maxval(specin%ang)>360.0d0)then 00945 allocate (temp(specin%nang)) 00946 Ashift=1 00947 temp=0.d0 00948 do i=1,specin%nang 00949 if (specin%ang(i)>360.d0) then 00950 specin%ang(i)=specin%ang(i)-360.0d0 00951 nt = nt+1 00952 endif 00953 enddo 00954 temp(nt+1:specin%nang)=specin%ang(1:specin%nang-nt) 00955 temp(1:nt)=specin%ang(specin%nang-nt+1:specin%nang) 00956 specin%ang=temp 00957 deallocate(temp) 00958 endif 00959 00960 ! convert to radians 00961 specin%ang=specin%ang*par_pi/180 00962 specin%dang=specin%ang(2)-specin%ang(1) 00963 00964 ! Skip Quant, next line, read VaDens or EnDens 00965 read(fid,'(a)',iostat=ier)rtext 00966 read(fid,'(a)',iostat=ier2)rtext 00967 read(fid,'(a)',iostat=ier3)rtext 00968 if (ier+ier2+ier3 .ne. 0) then 00969 call report_file_read_error(readfile) 00970 endif 00971 if (rtext == 'VaDens') then 00972 switch = 1 00973 elseif (rtext == 'EnDens') then 00974 switch = 2 00975 else 00976 call writelog('slwe','', 'SWAN VaDens/EnDens keyword not found') 00977 call halt_program 00978 end if 00979 read(fid,'(a)',iostat=ier)rtext 00980 read(fid,*,iostat=ier2)exc 00981 if (ier+ier2 .ne. 0) then 00982 call report_file_read_error(readfile) 00983 endif 00984 00985 i=0 00986 ! Find FACTOR keyword 00987 do while (i==0) 00988 read(fid,'(a)',iostat=ier)rtext 00989 if (ier .ne. 0) then 00990 call report_file_read_error(readfile) 00991 endif 00992 if (rtext == 'FACTOR') then 00993 i=1 00994 elseif (rtext == 'ZERO ') then 00995 call writelog('lswe','','Zero energy density input for this point') 00996 call halt_program 00997 elseif (rtext == 'NODATA') then 00998 call writelog('lwse','','SWAN file has no data for this point') 00999 call halt_program 01000 end if 01001 end do 01002 read(fid,*,iostat=ier)factor 01003 if (ier .ne. 0) then 01004 call report_file_read_error(readfile) 01005 endif 01006 01007 ! Read 2D S array 01008 allocate(specin%S(specin%nf,specin%nang)) 01009 do i=1,specin%nf 01010 read(fid,*,iostat=ier)(specin%S(i,ii),ii=1,specin%nang) 01011 if (ier .ne. 0) then 01012 call report_file_read_error(readfile) 01013 endif 01014 end do 01015 01016 ! Finished reading file 01017 close(fid) 01018 01019 ! Replace exception value 01020 where (specin%S == exc) 01021 specin%S = 0.d0 01022 endwhere 01023 01024 ! If angles were decreasing, flip S_array as also dir is flipped 01025 if (flipped) then 01026 call flipa(specin%S,specin%nf,specin%nang,2) 01027 end if 01028 01029 ! If the order of the angles in specin%ang was reordered, so the same in 01030 ! specin%S array 01031 if(Ashift==-1)then 01032 allocate(tempA(specin%nf,specin%nang)) 01033 tempA=0 01034 tempA(:,1:specin%nang-nt)=specin%S(:,nt+1:specin%nang) 01035 tempA(:,specin%nang-nt+1:specin%nang)=specin%S(:,1:nt) 01036 specin%S=tempA 01037 deallocate(tempA) 01038 elseif (Ashift==1) then 01039 allocate(tempA(specin%nf,specin%nang)) 01040 tempA=0 01041 tempA(:,nt+1:specin%nang)=specin%S(:,1:specin%nang-nt) 01042 tempA(:,1:nt)=specin%S(:,specin%nang-nt+1:specin%nang) 01043 specin%S=tempA 01044 deallocate(tempA) 01045 endif 01046 01047 ! multiply by SWAN output factor 01048 specin%S=specin%S*factor 01049 01050 ! Convert from energy density to variance density 01051 if (switch == 2) then 01052 specin%S=specin%S/(waveBoundaryParameters%rho*par_g) 01053 end if 01054 01055 ! Convert to m2/Hz/rad 01056 specin%S=specin%S*180/par_pi 01057 01058 ! We need a value for spreading. The assumption is that it is less than 1000 01059 ! This way, wp%fgen will not be set to just one angle. 01060 specin%scoeff = -1.d0 01061 ! We need to know if hm0 was set explicitly, not the case for Swan files 01062 specin%hm0 = -1.d0 01063 01064 ! We need frequency spectrum to ensure Sf remains correct between interpoloation 01065 ! routines 01066 allocate(specin%Sf(specin%nf)) 01067 specin%Sf=0.d0 01068 do i=1,specin%nf 01069 do ii=1,specin%nang 01070 if (ii==1) then 01071 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(2)-specin%ang(1)) 01072 elseif (ii==specin%nang) then 01073 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii)-specin%ang(ii-1)) 01074 else 01075 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii+1)-specin%ang(ii-1))/2 01076 endif 01077 enddo 01078 enddo 01079 01080 end subroutine read_swan_file 01081 01082 01083 ! -------------------------------------------------------------- 01084 ! -----------------Read variance density files ----------------- 01085 ! -------------------------------------------------------------- 01086 subroutine read_vardens_file(par,readfile,specin) 01087 01088 use params 01089 use logging_module 01090 use filefunctions 01091 use xmpi_module, only: Halt_Program 01092 use math_tools, only: flipv,flipa 01093 01094 IMPLICIT NONE 01095 01096 ! Input / output variables 01097 type(parameters), INTENT(IN) :: par 01098 character(len=*), intent(IN) :: readfile 01099 type(spectrum),intent(inout) :: specin 01100 01101 ! Internal variables 01102 integer :: fid,i,ii,nnz,ier,nt,Ashift 01103 real*8,dimension(:),allocatable :: Sd,temp 01104 real*8, dimension(:,:),allocatable :: tempA 01105 01106 ! Open file to start read 01107 call writelog('sl','','Reading from vardens file ',trim(readfile),' ...') 01108 fid = create_new_fid() 01109 open(fid,file=readfile,form='formatted',status='old') 01110 01111 ! Read number of frequencies and frequency vector 01112 read(fid,*,iostat=ier)specin%nf 01113 if (ier .ne. 0) then 01114 call report_file_read_error(readfile) 01115 endif 01116 allocate(specin%f(specin%nf)) 01117 do i=1,specin%nf 01118 read(fid,*,iostat=ier)specin%f(i) 01119 if (ier .ne. 0) then 01120 call report_file_read_error(readfile) 01121 endif 01122 end do 01123 01124 ! Read number of angles and angles vector 01125 read(fid,*,iostat=ier)specin%nang 01126 if (ier .ne. 0) then 01127 call report_file_read_error(readfile) 01128 endif 01129 allocate(specin%ang(specin%nang)) 01130 do i=1,specin%nang 01131 read(fid,*,iostat=ier)specin%ang(i) 01132 if (ier .ne. 0) then 01133 call report_file_read_error(readfile) 01134 endif 01135 end do 01136 ! Directions should span 0-360 degrees, not negative 01137 nt = 0 01138 Ashift = 0 01139 ! Make sure that all angles are in range of 0 to 360 degrees 01140 if(minval(specin%ang)<0.d0)then 01141 allocate (temp(specin%nang)) 01142 Ashift=-1 01143 temp=0.d0 01144 do i=1,specin%nang 01145 if (specin%ang(i)<0.d0) then 01146 specin%ang(i)=specin%ang(i)+360.0d0 01147 nt = nt+1 01148 endif 01149 enddo 01150 temp(1:specin%nang-nt)=specin%ang(nt+1:specin%nang) 01151 temp(specin%nang-nt+1:specin%nang)=specin%ang(1:nt) 01152 specin%ang=temp 01153 deallocate(temp) 01154 elseif(maxval(specin%ang)>360.0d0)then 01155 allocate (temp(specin%nang)) 01156 Ashift=1 01157 temp=0.d0 01158 do i=1,specin%nang 01159 if (specin%ang(i)>360.d0) then 01160 specin%ang(i)=specin%ang(i)-360.0d0 01161 nt = nt+1 01162 endif 01163 enddo 01164 temp(nt+1:specin%nang)=specin%ang(1:specin%nang-nt) 01165 temp(1:nt)=specin%ang(specin%nang-nt+1:specin%nang) 01166 specin%ang=temp 01167 deallocate(temp) 01168 endif 01169 01170 ! Convert from degrees to rad 01171 specin%ang=specin%ang*par_pi/180 01172 specin%dang=specin%ang(2)-specin%ang(1) 01173 01174 ! Read 2D S array 01175 allocate(specin%S(specin%nf,specin%nang)) 01176 do i=1,specin%nf 01177 read(fid,*,iostat=ier)(specin%S(i,ii),ii=1,specin%nang) 01178 if (ier .ne. 0) then 01179 call report_file_read_error(readfile) 01180 endif 01181 end do 01182 01183 ! If the order of the angles in specin%ang was reordered, so the same in 01184 ! specin%S array 01185 if(Ashift==-1)then 01186 allocate(tempA(specin%nf,specin%nang)) 01187 tempA=0 01188 tempA(:,1:specin%nang-nt)=specin%S(:,nt+1:specin%nang) 01189 tempA(:,specin%nang-nt+1:specin%nang)=specin%S(:,1:nt) 01190 specin%S=tempA 01191 deallocate(tempA) 01192 elseif (Ashift==1) then 01193 allocate(tempA(specin%nf,specin%nang)) 01194 tempA=0 01195 tempA(:,nt+1:specin%nang)=specin%S(:,1:specin%nang-nt) 01196 tempA(:,1:nt)=specin%S(:,specin%nang-nt+1:specin%nang) 01197 specin%S=tempA 01198 deallocate(tempA) 01199 endif 01200 01201 ! Finished reading file 01202 close(fid) 01203 01204 ! Convert to m2/Hz/rad 01205 specin%S=specin%S*180/par_pi 01206 01207 ! We need a value for spreading. The assumption is that it is less than 1000 01208 ! if more than one direction column has variance. If only one column has variance 01209 ! we assume the model requires long crested waves, and wp%thetagen should be exactly 01210 ! equal to the input direction. We then also need to find the dominant angle. 01211 ! First flatten spectrum to direction only, then determine if all but one are zero, 01212 ! use the direction of the non-zero as the main wave direction, set spreading to 01213 ! a high value (greater than 1000). Else set spreading to a negative value. 01214 allocate(Sd(specin%nang)) 01215 Sd = sum(specin%S,DIM = 1) 01216 nnz = count(Sd>0.d0) 01217 if (nnz == 1) then 01218 specin%scoeff = 1024.d0 01219 specin%dir0 = specin%ang(minval(maxloc(Sd))) 01220 else 01221 specin%scoeff = -1.d0 01222 endif 01223 01224 ! We need frequency spectrum to ensure Sf remains correct between interpoloation 01225 ! routines 01226 allocate(specin%Sf(specin%nf)) 01227 specin%Sf=0.d0 01228 do i=1,specin%nf 01229 do ii=1,specin%nang 01230 if (ii==1) then 01231 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(2)-specin%ang(1)) 01232 elseif (ii==specin%nang) then 01233 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii)-specin%ang(ii-1)) 01234 else 01235 specin%Sf(i) = specin%Sf(i)+specin%S(i,ii)*abs(specin%ang(ii+1)-specin%ang(ii-1))/2 01236 endif 01237 enddo 01238 enddo 01239 01240 01241 ! We need to know if hm0 was set explicitly, not the case for vardens files 01242 specin%hm0 = -1.d0 01243 01244 ! Free memory 01245 deallocate(Sd) 01246 end subroutine read_vardens_file 01247 01248 01249 ! -------------------------------------------------------------- 01250 ! ------------- Interpolate to standard spectrum --------------- 01251 ! -------------------------------------------------------------- 01252 subroutine interpolate_spectrum(specin,specinterp,par,fmax) 01253 01254 use interp 01255 use params 01256 01257 IMPLICIT NONE 01258 01259 ! input/output 01260 type(spectrum),intent(in) :: specin 01261 type(spectrum),intent(inout) :: specinterp 01262 type(parameters),intent(in) :: par 01263 real*8, intent(in) :: fmax 01264 ! internal 01265 integer :: i,j,dummy 01266 real*8 :: m0,df,dang 01267 real*8,dimension(naint) :: Sd 01268 real*8 :: hm0pre,hm0post,Sfnow,factor,tempt0 01269 01270 ! allocate size of f,ang,Sf and S arrays in specinterp 01271 allocate(specinterp%f(nfint)) 01272 allocate(specinterp%ang(naint)) 01273 allocate(specinterp%Sf(nfint)) 01274 allocate(specinterp%S(nfint,naint)) 01275 01276 ! fill f and ang arrays 01277 specinterp%nf = nfint 01278 specinterp%df = fmax/(nfint-1) ! this usually gives a range of 0-1Hz, 01279 ! with 1/800Hz resolution 01280 do i=1,nfint 01281 specinterp%f(i)=(i-1)*specinterp%df 01282 enddo 01283 specinterp%nang = naint 01284 specinterp%dang = 2*par_pi/(naint-1) ! this is exactly the same as in the JONSWAP construction 01285 do i=1,specinterp%nang 01286 specinterp%ang(i)=(i-1)*specinterp%dang 01287 enddo 01288 01289 ! If hm0 was set explicitly, then use that, else calculate hm0 01290 if (specin%hm0>0.d0) then 01291 hm0pre = specin%hm0 01292 else 01293 ! pre-interpolation hm0 value (can be on a non-monotonic f,ang grid) 01294 m0 = 0 01295 do j=1,specin%nang 01296 do i=1,specin%nf 01297 df = specin%f(max(2,i))-specin%f(max(1,i-1)) 01298 dang = specin%ang(max(2,j))-specin%ang(max(1,j-1)) 01299 m0 = m0 + specin%S(i,j)*df*dang 01300 enddo 01301 enddo 01302 hm0pre = 4*sqrt(m0) 01303 endif 01304 01305 ! interpolation (no extrapolation) of input 2D spectrum to standard 2D spectrum 01306 do j=1,naint 01307 do i=1,nfint 01308 call linear_interp_2d(specin%f,specin%nf,specin%ang,specin%nang,specin%S, & 01309 specinterp%f(i),specinterp%ang(j),specinterp%S(i,j),'interp',0.d0) 01310 enddo 01311 enddo 01312 ! hm0 post (is always on a monotonic f,ang grid 01313 m0 = sum(specinterp%S)*specinterp%df*specinterp%dang 01314 hm0post = 4*sqrt(m0) 01315 01316 ! correct the wave energy to ensure hm0post == hm0pre 01317 ! specinterp%S = (hm0pre/hm0post)**2*specinterp%S 01318 ! This is done later using Sf 01319 01320 ! calculate 1D spectrum, summed over directions 01321 specinterp%Sf = sum(specinterp%S, DIM = 2)*specinterp%dang 01322 01323 ! correct the wave energy by setting Sfpost == Sfpre 01324 do i=1,nfint 01325 if (specinterp%f(i)>=minval(specin%f) .and. specinterp%f(i)<=maxval(specin%f)) then 01326 call LINEAR_INTERP(specin%f,specin%Sf,specin%nf,specinterp%f(i),Sfnow,dummy) 01327 if (specinterp%Sf(i)>0.d0 .and. Sfnow>0.d0) then 01328 factor = Sfnow/specinterp%Sf(i) 01329 specinterp%Sf(i) = specinterp%Sf(i)*factor 01330 specinterp%S(i,:) = specinterp%S(i,:)*factor 01331 elseif (Sfnow==0.d0) then 01332 specinterp%Sf(i) = 0.d0 01333 specinterp%S(i,:) = 0.d0 01334 else 01335 specinterp%Sf(i) = Sfnow 01336 dummy = maxval(maxloc(specin%S(i,:))) 01337 tempt0 = specin%ang(dummy) 01338 dummy = maxval(minloc(abs(specinterp%ang-tempt0))) 01339 specinterp%S(i,dummy) = Sfnow/specinterp%dang 01340 endif 01341 endif 01342 enddo 01343 01344 ! calculate other wave statistics from interpolated spectrum 01345 m0 = sum(specinterp%Sf)*specinterp%df 01346 specinterp%hm0 = 4*sqrt(m0) 01347 Sd = sum(specinterp%S, DIM = 1)*specinterp%df 01348 i = maxval(maxloc(Sd)) 01349 specinterp%dir0 = 270.d0 - specinterp%ang(i)*180/par_pi ! converted back into nautical degrees 01350 i = maxval(maxloc(specinterp%Sf)) 01351 specinterp%fp = specinterp%f(i) 01352 call tpDcalc(specinterp%Sf,specinterp%f,specinterp%trep,par%trepfac,par%Tm01switch) 01353 specinterp%dirm = 270.d0-180.d0/par_pi*atan2( sum(sin(specinterp%ang)*Sd)/sum(Sd),& 01354 sum(cos(specinterp%ang)*Sd)/sum(Sd) ) 01355 01356 end subroutine interpolate_spectrum 01357 01358 ! -------------------------------------------------------------- 01359 ! ----------- Small subroutine to determine if the ------------ 01360 ! ---------- global repeatwbc should be true or false ---------- 01361 subroutine set_repeatwbc 01362 01363 implicit none 01364 ! internal 01365 integer :: i 01366 01367 ! set default to repeat all boundary conditions 01368 waveSpectrumAdministration%repeatwbc = .true. 01369 ! find any spectrum that cannot be repeated, then change global 01370 ! repeatwbc to false 01371 do i=1,nspectra 01372 if (.not. bcfiles(i)%repeat) then 01373 waveSpectrumAdministration%repeatwbc = .false. 01374 endif 01375 enddo 01376 01377 end subroutine set_repeatwbc 01378 01379 ! -------------------------------------------------------------- 01380 ! ----------- Small subroutine to set the filenames ------------ 01381 ! ----------- of the boundary condition output files ----------- 01382 subroutine set_bcfilenames(wp) 01383 01384 implicit none 01385 ! input/output 01386 type(waveparamsnew),intent(inout) :: wp 01387 ! internal 01388 integer :: i1,i2,i3,i4,i5 01389 01390 if (waveSpectrumAdministration%repeatwbc) then 01391 wp%Efilename = 'E_reuse.bcf' 01392 wp%qfilename = 'q_reuse.bcf' 01393 wp%nhfilename = 'nh_reuse.bcf' 01394 else 01395 i1=floor(real(bccount)/10000) 01396 i2=floor(real(bccount-i1*10000)/1000) 01397 i3=floor(real(bccount-i1*10000-i2*1000)/100) 01398 i4=floor(real(bccount-i1*10000-i2*1000-i3*100)/10) 01399 i5=bccount-i1*10000-i2*1000-i3*100-i4*10 01400 wp%Efilename='E_series'//char(48+i1)//char(48+i2)//char(48+i3)//char(48+i4)//char(48+i5)//'.bcf' 01401 wp%qfilename='q_series'//char(48+i1)//char(48+i2)//char(48+i3)//char(48+i4)//char(48+i5)//'.bcf' 01402 wp%nhfilename='nh_series'//char(48+i1)//char(48+i2)//char(48+i3)//char(48+i4)//char(48+i5)//'.bcf' 01403 endif 01404 01405 end subroutine set_bcfilenames 01406 01407 ! -------------------------------------------------------------- 01408 ! ----------- Merge all separate spectra into one -------------- 01409 ! -------------- average spectrum for other use ---------------- 01410 subroutine generate_combined_spectrum(specinterp,combspec) 01411 01412 implicit none 01413 type(spectrum),dimension(nspectra),intent(in) :: specinterp 01414 type(spectrum),intent(inout) :: combspec 01415 integer :: iloc 01416 real*8,dimension(naint) :: Sd 01417 real*8,dimension(3) :: peakSd,peakang 01418 real*8,dimension(:),allocatable :: tempSf 01419 01420 allocate(combspec%f(nfint)) 01421 allocate(combspec%Sf(nfint)) 01422 allocate(combspec%ang(naint)) 01423 allocate(combspec%S(nfint,naint)) 01424 combspec%f=specinterp(1)%f 01425 combspec%nf=nfint 01426 combspec%df=specinterp(1)%df 01427 combspec%ang=specinterp(1)%ang 01428 combspec%nang=naint 01429 combspec%dang=specinterp(1)%dang 01430 combspec%trep = 0.d0 01431 combspec%S=0.d0 01432 combspec%Sf=0.d0 01433 combspec%hm0=0.d0 01434 01435 do iloc = 1,nspectra 01436 combspec%S = combspec%S +specinterp(iloc)%S /nspectra 01437 combspec%Sf = combspec%Sf +specinterp(iloc)%Sf /nspectra 01438 enddo 01439 ! 01440 ! Calculate peak wave angle 01441 ! 01442 ! frequency integrated variance array 01443 Sd = sum(combspec%S,DIM=1) 01444 ! peak location of f-int array 01445 iloc = maxval(maxloc(Sd)) 01446 ! pick two neighbouring directional bins, including effect of closing circle at 01447 ! 0 and 2pi rad 01448 if (iloc>1 .and. iloc<naint) then 01449 peakSd = (/Sd(iloc-1),Sd(iloc),Sd(iloc+1)/) 01450 peakang = (/combspec%ang(iloc-1),combspec%ang(iloc),combspec%ang(iloc+1)/) 01451 elseif (iloc==1) then 01452 peakSd = (/Sd(naint),Sd(1),Sd(2)/) 01453 peakang = (/combspec%ang(naint)-2*par_pi,combspec%ang(1),combspec%ang(2)/) 01454 elseif (iloc==naint) then 01455 peakSd = (/Sd(naint-1),Sd(naint),Sd(1)/) 01456 peakang = (/combspec%ang(naint-1),combspec%ang(naint),combspec%ang(1)+2*par_pi/) 01457 endif 01458 ! dir0 calculated as mean over peak and two neighbouring cells 01459 combspec%dir0 = sum(peakSd*peakang)/sum(peakSd) 01460 ! return to 0<=dir0<=2*par_pi 01461 if (combspec%dir0>2*par_pi) then 01462 combspec%dir0 = combspec%dir0-2*par_pi 01463 elseif (combspec%dir0<0) then 01464 combspec%dir0 = combspec%dir0+2*par_pi 01465 endif 01466 ! 01467 ! Compute peak wave frequency 01468 ! 01469 allocate(tempSf(size(combspec%Sf))) 01470 tempSf=0*combspec%Sf 01471 ! find frequency range of 95% wave energy around peak 01472 where (combspec%Sf>0.95d0*maxval(combspec%Sf)) 01473 tempSf=combspec%Sf 01474 end where 01475 ! smoothed peak is weighted mean frequency of 95% range 01476 combspec%fp = sum(tempSf*combspec%f)/sum(tempSf) 01477 ! 01478 ! Compute representative wave period 01479 ! 01480 call tpDcalc(combspec%Sf,combspec%f,combspec%trep,par%trepfac,par%Tm01switch) 01481 ! 01482 ! Compute combined wave height 01483 do iloc = 1,nspectra 01484 combspec%hm0 = combspec%hm0+specinterp(iloc)%hm0**2/nspectra 01485 enddo 01486 combspec%hm0 = sqrt(combspec%hm0) 01487 01488 end subroutine generate_combined_spectrum 01489 01490 01491 ! -------------------------------------------------------------- 01492 ! ----------- Choose wave train components based on ------------ 01493 ! --------------------- combined spectrum ---------------------- 01494 subroutine generate_wavetrain_components(combspec,wp) 01495 01496 #ifdef BUILDXBEACH 01497 use logging_module 01498 use interp 01499 use wave_functions_module, only: iteratedispersion 01500 #endif 01501 use wave_boundary_datastore 01502 use math_tools 01503 01504 implicit none 01505 ! input/output 01506 type(spectrum),intent(in) :: combspec 01507 type(waveparamsnew),intent(inout) :: wp 01508 ! internal 01509 integer :: i,ii 01510 integer :: ind1,ind2,dummy 01511 real*8,dimension(:),allocatable :: randnums,pdflocal,cdflocal 01512 real*8 :: L0,L,kmax,fmax,dummy_real 01513 01514 01515 ! If we are running non-hydrostatic boundary conditions, we want to remove 01516 ! very high frequency components, as these will not be computed well anyway 01517 if (waveBoundaryParameters%nonhspectrum) then 01518 kmax = wdmax/hb0 01519 fmax = par_g*kmax*tanh(kmax*hb0)/2/par_pi 01520 else 01521 ! this is really already taken into account 01522 ! by specifying waveBoundaryParameters%sprdthr 01523 fmax = 2*maxval(combspec%f) 01524 endif 01525 01526 ! Determine frequencies around peak frequency of one-dimensional 01527 ! non-directional variance density spectrum, based on factor sprdthr, which 01528 ! should be included in the determination of the wave boundary conditions 01529 call frange(combspec%f,combspec%Sf,fmax,ind1,ind2) 01530 01531 ! Calculate number of wave components to be included in determination of the 01532 ! wave boundary conditions based on the wave record length and width of the 01533 ! wave frequency range 01534 wp%K = ceiling(wp%rtbc*(combspec%f(ind2)-combspec%f(ind1))+1) 01535 ! also include minimum number of components 01536 wp%K = max(wp%K,Kmin) 01537 01538 ! Allocate space in waveparams for all wave train components 01539 allocate(wp%fgen(wp%K)) 01540 allocate(wp%thetagen(wp%K)) 01541 allocate(wp%phigen(wp%K)) 01542 allocate(wp%kgen(wp%K)) 01543 allocate(wp%wgen(wp%K)) 01544 01545 ! Select equidistant wave components between the earlier selected range of 01546 ! frequencies around the peak frequency in the frange subfunction 01547 wp%dfgen = (combspec%f(ind2)-combspec%f(ind1))/(wp%K-1) 01548 do i=1,wp%K 01549 wp%fgen(i)=combspec%f(ind1)+(i-1)*wp%dfgen 01550 enddo 01551 01552 ! This subroutine needs to generate random phase / directions for the individual 01553 ! wave trains. Due to some strange Fortran properties, it is better to select 01554 ! one long vector with random numbers for both the phase and the direction, than 01555 ! one vector for each. 01556 ! Update random seed, if requested 01557 allocate(randnums(2*wp%K)) 01558 ! Spin up random number generator, else all low seed numbers mean random sequence 01559 ! start at ~0. Don't know how much spin up is required (strictly no spin-up needed 01560 ! if random integer given as seed) 01561 dummy_real = random(waveBoundaryParameters%randomseed) 01562 do i=1,10 01563 dummy_real = random(0) 01564 enddo 01565 do i=1,2*wp%K 01566 randnums(i) = random(0) 01567 enddo 01568 01569 ! Determine a random phase for each wave train component between 0 and 2pi 01570 wp%phigen=randnums(1:wp%K)*2*par_pi 01571 01572 ! Determine random directions for each wave train component, based on the CDF of 01573 ! the directional spectrum. For each wave train we will interpolate the directional 01574 ! distribution at that frequency, generate a CDF, and then interpolate the wave 01575 ! train direction from a random number draw and the CDF. 01576 allocate(pdflocal(naint)) 01577 allocate(cdflocal(naint)) 01578 do i=1,wp%K 01579 ! interpolate spectrum at this frequency per angle in the spectrum 01580 do ii=1,naint 01581 #ifdef BUILDXBEACH 01582 call LINEAR_INTERP(combspec%f,combspec%S(:,ii),combspec%nf, wp%fgen(i),pdflocal(ii),dummy) 01583 #else 01584 ! interpolate in other systems 01585 #endif 01586 enddo 01587 ! convert to pdf by ensuring total integral == 1, assuming constant directional bin size 01588 pdflocal = pdflocal/sum(pdflocal) 01589 ! convert to cdf by trapezoidal integration, which is not biased for integration 01590 ! direction. 01591 ! Boundary condition, which may be nonzero: 01592 cdflocal(1) = pdflocal(1) 01593 do ii=2,naint 01594 cdflocal(ii) = cdflocal(ii-1) + (pdflocal(ii)+pdflocal(ii-1))/2 01595 ! Note: this only works if the directional 01596 ! bins are constant in size. Assumed multiplication 01597 ! by one. 01598 enddo 01599 ! interpolate random number draw across the cdf. Note that cdf(1) is not assumed to be zero and 01600 ! there is a posibility of drawing less than cdf(1). Therefore, if the random number is .lt. cdf(1), 01601 ! interpolation should take place across the back of the angle spectrum 01602 ! (i.e., from theta(end)-2pi : theta(1)). 01603 ! Note that we are using the second half of randnums now (K+1:end) 01604 if (randnums(wp%K+i)>=cdflocal(1)) then 01605 #ifdef BUILDXBEACH 01606 call LINEAR_INTERP(cdflocal,combspec%ang,naint,randnums(wp%K+i), wp%thetagen(i),dummy) 01607 #else 01608 ! interpolate in other systems 01609 #endif 01610 else 01611 #ifdef BUILDXBEACH 01612 call LINEAR_INTERP((/0.d0,cdflocal(1)/),(/combspec%ang(naint)-2*par_pi,combspec%ang(1)/), & 01613 2,randnums(wp%K+i),wp%thetagen(i),dummy) 01614 #else 01615 ! interpolate in other systems 01616 #endif 01617 endif 01618 ! ensure wave direction 0<=theta<2pi 01619 wp%thetagen(i) = mod(wp%thetagen(i),2*par_pi) 01620 enddo 01621 01622 ! determine wave number for each wave train component, using standard dispersion relation 01623 ! solver from wave_functions module. This function returns a negative wave length if the 01624 ! solver did not converge. 01625 #ifdef BUILDXBEACH 01626 do i=1,wp%K 01627 L0 = par_g*(1/wp%fgen(i))**2/2/par_pi ! deep water wave length 01628 L = iteratedispersion(L0,L0,par_pi,hb0) 01629 if (L<0.d0) then 01630 call writelog('lsw','','No dispersion convergence found for wave train ',i, & 01631 ' in boundary condition generation') 01632 L = -L 01633 endif 01634 wp%kgen(i) = 2*par_pi/L 01635 enddo 01636 #else 01637 ! do dispersion relation in other model 01638 #endif 01639 01640 ! Angular frequency 01641 wp%wgen = 2*par_pi*wp%fgen 01642 01643 ! Free memory 01644 deallocate(pdflocal,cdflocal,randnums) 01645 01646 end subroutine generate_wavetrain_components 01647 01648 ! ----------------------------------------------------------- 01649 ! ----------- Small subroutine to determine ----------------- 01650 ! ------------ representative wave period ------------------- 01651 subroutine tpDcalc(Sf,f,Trep,trepfac,switch) 01652 01653 implicit none 01654 01655 real*8, dimension(:), intent(in) :: Sf, f 01656 real*8, intent(out) :: Trep 01657 real*8, intent(in) :: trepfac 01658 integer, intent(in) :: switch 01659 01660 real*8, dimension(:),allocatable :: temp 01661 01662 allocate(temp(size(Sf))) 01663 temp=0.d0 01664 where (Sf>=trepfac*maxval(Sf)) 01665 temp=1.d0 01666 end where 01667 01668 if (switch == 1) then 01669 Trep=sum(temp*Sf)/sum(temp*Sf*f) ! Tm01 01670 else 01671 Trep = sum(temp*Sf/max(f,0.001d0))/sum(temp*Sf) ! Tm-1,0 01672 endif 01673 01674 deallocate(temp) 01675 01676 end subroutine tpDcalc 01677 01678 01679 ! ----------------------------------------------------------- 01680 ! ---- Small subroutine to determine f-range round peak ----- 01681 ! ----------------------------------------------------------- 01682 subroutine frange(f,Sf,fmax,firstp,lastp,findlineout) 01683 01684 use wave_boundary_datastore 01685 01686 implicit none 01687 01688 real*8, dimension(:), intent(in) :: f,Sf 01689 real*8,intent(in) :: fmax 01690 01691 integer, intent(out) :: firstp, lastp 01692 01693 real*8, dimension(:), intent(out),allocatable,optional :: findlineout 01694 real*8, dimension(:),allocatable :: temp, findline 01695 integer :: i = 0 01696 01697 ! find frequency range around peak 01698 allocate(findline(size(Sf))) 01699 findline=0*Sf 01700 where (Sf>waveBoundaryParameters%sprdthr*maxval(Sf)) 01701 findline=1 01702 end where 01703 01704 ! find frequency range below fmax 01705 where(f>fmax) 01706 findline = 0 01707 endwhere 01708 01709 firstp=maxval(maxloc(findline)) ! Picks the first "1" in temp 01710 01711 allocate (temp(size(findline))) 01712 temp=(/(i,i=1,size(findline))/) 01713 lastp=maxval(maxloc(temp*findline)) ! Picks the last "1" in temp 01714 01715 if (present(findlineout)) then 01716 allocate(findlineout(size(Sf))) 01717 findlineout=findline 01718 endif 01719 deallocate(temp, findline) 01720 01721 end subroutine frange 01722 01723 01724 ! -------------------------------------------------------------- 01725 ! ---------------- Setup time axis for waves ------------------- 01726 ! -------------------------------------------------------------- 01727 subroutine generate_wave_time_axis(wp) 01728 01729 use wave_boundary_datastore 01730 01731 implicit none 01732 ! input/output 01733 type(waveparamsnew),intent(inout) :: wp 01734 ! internal 01735 integer :: ntaper, indend 01736 integer :: i 01737 01738 01739 01740 ! First assume that internal and bc-writing time step is the same 01741 wp%dtin = wp%dtbc 01742 wp%dtchanged = .false. 01743 wp%tslenbc = nint(wp%rtbc/wp%dtbc)+1 01744 01745 ! Check whether the internal frequency is high enough to describe the highest frequency 01746 ! wave train returned from frange (which can be used in the boundary conditions) 01747 if (.not.waveBoundaryParameters%nonhspectrum) then 01748 if (wp%dtin>0.5d0/wp%fgen(wp%K)) then 01749 wp%dtin = 0.5d0/wp%fgen(wp%K) 01750 wp%dtchanged = .true. 01751 endif 01752 else 01753 if (wp%dtin>0.1d0/wp%fgen(wp%K)) then 01754 wp%dtin = 0.1d0/wp%fgen(wp%K) 01755 wp%dtchanged = .true. 01756 endif 01757 endif 01758 01759 ! The length of the internal time axis should be even (for Fourier transform) and 01760 ! depends on the internal time step needed and the internal duration (~1/dfgen): 01761 wp%tslen = ceiling(1/wp%dfgen/wp%dtin)+1 01762 if (mod(wp%tslen,2)/=0) then 01763 wp%tslen = wp%tslen +1 01764 end if 01765 01766 ! Now we can make the internal time axis 01767 wp%rtin = wp%tslen * wp%dtin 01768 allocate(wp%tin(wp%tslen)) 01769 do i=1,wp%tslen 01770 wp%tin(i) = (i-1)*wp%dtin 01771 enddo 01772 01773 ! Make a taper function to slowly increase and decrease the boundary condition forcing 01774 ! at the start and the end of the boundary condition file (including any time beyond 01775 ! the external rtbc 01776 allocate(wp%taperf(wp%tslen)) 01777 allocate(wp%taperw(wp%tslen)) 01778 ! fill majority with unity 01779 wp%taperf = 1.d0 01780 wp%taperw = 1.d0 01781 if (waveBoundaryParameters%nonhspectrum) then 01782 ! begin taper by building up the wave conditions over 2 wave periods 01783 ntaper = nint((2.0d0*waveSpectrumAdministration%Tbc)/wp%dtin) 01784 else 01785 ntaper = nint((5.d0*waveSpectrumAdministration%Tbc)/wp%dtin) 01786 endif 01787 do i=1,min(ntaper,size(wp%taperf)) 01788 wp%taperf(i) = tanh(5.d0*i/ntaper) ! multiplied by five because tanh(5)=~1 01789 wp%taperw(i) = tanh(5.d0*i/ntaper) 01790 enddo 01791 ! We do not want to taperw the end anymore. Instead we pass the wave height at the end of rtbc to 01792 ! the next wave generation iteration. 01793 ! end taper by finding where tin=rtbc, taper before that and set everything to zero after 01794 ! that. 01795 if (wp%tin(wp%tslen)>wp%rtbc) then 01796 indend = minval(minloc(wp%tin,MASK = wp%tin>=wp%rtbc)) 01797 else 01798 indend = wp%tslen 01799 endif 01800 do i=1,min(ntaper,indend) 01801 wp%taperf(indend+1-i) = min(wp%taperf(indend+1-i),tanh(5.d0*i/ntaper)) 01802 enddo 01803 wp%taperf(indend:wp%tslen) = 0.d0 01804 ind_end_taper = indend 01805 end subroutine generate_wave_time_axis 01806 01807 01808 ! -------------------------------------------------------------- 01809 ! -------- Calculate variance at each spectrum location -------- 01810 ! -------------------------------------------------------------- 01811 subroutine generate_wave_train_variance(wp,specinterp) 01812 01813 #ifdef BUILDXBEACH 01814 use interp 01815 #endif 01816 01817 implicit none 01818 ! input/output 01819 type(waveparamsnew),intent(inout) :: wp 01820 type(spectrum),dimension(nspectra),intent(in):: specinterp 01821 ! internal 01822 integer :: i,ii, dummy 01823 real*8 :: hm0post 01824 01825 ! allocate space for the variance arrays 01826 allocate(wp%vargen(nspectra)) 01827 allocate(wp%vargenq(nspectra)) 01828 01829 ! Determine variance at each spectrum location 01830 do i=1,nspectra 01831 allocate(wp%vargen(i)%Sf(wp%K)) 01832 allocate(wp%vargenq(i)%Sf(wp%K)) 01833 do ii=1,wp%K 01834 ! In order to maintain energy density per frequency, an interpolation 01835 ! is carried out on the input directional variance density spectrum 01836 ! at the frequency and direction locations in fgen. 01837 call linear_interp(specinterp(i)%f, & 01838 specinterp(i)%Sf, & 01839 nfint, & 01840 wp%fgen(ii), & 01841 wp%vargen(i)%Sf(ii),dummy) 01842 enddo 01843 ! Correct variance to ensure the total variance remains the same as the total variance in the 01844 ! (interpolated) input spectrum 01845 hm0post = 4*sqrt(sum(wp%vargen(i)%Sf)*wp%dfgen) 01846 wp%vargen(i)%Sf = (specinterp(i)%hm0/hm0post)**2*wp%vargen(i)%Sf 01847 ! For the generation of long waves we cannot use wp%vargen%Sf, because it contains an overestimation 01848 ! of energy in the peak frequencies. We can also not use the standard directionally-integrated spectrum 01849 ! because this stores all energy at fgen(K), where it is possible that for this current spectrum, there 01850 ! is no energy at S(f(K),theta(K0)) 01851 ! The current solution is to take the minimum of both methods 01852 do ii=1,wp%K 01853 ! Map Sf input to fgen 01854 call LINEAR_INTERP(specinterp(i)%f,specinterp(i)%Sf,specinterp(i)%nf,wp%fgen(ii),wp%vargenq(i)%Sf(ii),dummy) 01855 enddo 01856 wp%vargenq(i)%Sf = min(wp%vargen(i)%Sf,wp%vargenq(i)%Sf) 01857 enddo 01858 01859 end subroutine generate_wave_train_variance 01860 01861 ! -------------------------------------------------------------- 01862 ! ------ Calculate amplitudes at each spectrum location -------- 01863 ! ------------ for every wave train component ------------------ 01864 subroutine generate_wave_train_properties_per_offshore_point(wp,s) 01865 01866 use spaceparams 01867 use interp 01868 01869 implicit none 01870 ! input/output 01871 type(waveparamsnew),intent(inout) :: wp 01872 type(spacepars),intent(in) :: s 01873 ! internal 01874 integer :: i,ii,dummy 01875 integer,dimension(2) :: interpindex 01876 real*8,dimension(2) :: positions,sf 01877 logical :: interpcase 01878 real*8 :: here,sfnow,sfimp 01879 integer,dimension(size(n_index_loc)) :: temp_index_loc 01880 01881 ! allocate space for the amplitude array and representative integration angle 01882 allocate(wp%A(s%ny+1,wp%K)) 01883 ! allocate(wp%danggen(s%ny+1)) 01884 allocate(wp%Sfinterp(s%ny+1,wp%K)) 01885 allocate(wp%Sfinterpq(s%ny+1,wp%K)) 01886 allocate(wp%Hm0interp(s%ny+1)) 01887 01888 ! where necessary, interpolate Sf of each spectrum location 01889 ! to the current grid cell, and use this to calculate A 01890 do i=1,s%ny+1 01891 ! Four possibilities: 01892 ! 1: this exact point has an input spectrum 01893 ! 2: this point lies between two input spectra 01894 ! 3: this point lies before any input spectrum 01895 ! 4: this point lies beyond any input spectrum 01896 if (any(n_index_loc==i)) then ! case 1 01897 interpcase = .false. 01898 do ii=1,nspectra 01899 if(n_index_loc(ii)==i) then 01900 interpindex = ii 01901 endif 01902 enddo 01903 elseif(i<minval(n_index_loc)) then ! case 3 01904 interpcase = .false. 01905 do ii=1,nspectra 01906 if(n_index_loc(ii)==minval(n_index_loc)) then 01907 interpindex = ii 01908 endif 01909 enddo 01910 elseif(i>maxval(n_index_loc)) then ! case 4 01911 interpcase = .false. 01912 do ii=1,nspectra 01913 if(n_index_loc(ii)==maxval(n_index_loc)) then 01914 interpindex = ii 01915 endif 01916 enddo 01917 else ! case 2 01918 interpcase = .true. 01919 temp_index_loc = n_index_loc 01920 where(n_index_loc>i) 01921 temp_index_loc = -huge(0) 01922 endwhere 01923 interpindex(1) = minval(maxloc(temp_index_loc)) 01924 01925 temp_index_loc = n_index_loc 01926 where(n_index_loc<i) 01927 temp_index_loc = huge(0) 01928 endwhere 01929 interpindex(2) = minval(minloc(temp_index_loc)) 01930 endif 01931 01932 ! interpolate Sf 01933 if (interpcase) then 01934 positions(1) = s%ndist(1,n_index_loc(interpindex(1))) 01935 positions(2) = s%ndist(1,n_index_loc(interpindex(2))) 01936 here = s%ndist(1,i) 01937 do ii=1,wp%K 01938 sf(1) = wp%vargen(interpindex(1))%Sf(ii) 01939 sf(2) = wp%vargen(interpindex(2))%Sf(ii) 01940 call LINEAR_INTERP(positions,sf,2,here,wp%Sfinterp(i,ii),dummy) 01941 sf(1) = wp%vargenq(interpindex(1))%Sf(ii) 01942 sf(2) = wp%vargenq(interpindex(2))%Sf(ii) 01943 call LINEAR_INTERP(positions,sf,2,here,wp%Sfinterpq(i,ii),dummy) 01944 enddo 01945 ! Ensure that the total amount of variance has not been reduced/increased 01946 ! during interpolation for short wave energy only. This is not done for 01947 ! Sf for bound long wave generation 01948 sfnow = sum(wp%Sfinterp(i,:)) ! integration over fixed bin df size unnecessary 01949 sf(1) = sum(wp%vargen(interpindex(1))%Sf) 01950 sf(2) = sum(wp%vargen(interpindex(2))%Sf) 01951 call LINEAR_INTERP(positions,sf,2,here,sfimp,dummy) 01952 wp%Sfinterp(i,:)=wp%Sfinterp(i,:)*sfimp/sfnow 01953 else 01954 wp%Sfinterp(i,:) = wp%vargen(interpindex(1))%Sf 01955 wp%Sfinterpq(i,:) = wp%vargenq(interpindex(1))%Sf 01956 endif 01957 01958 wp%A(i,:) = sqrt(2*wp%Sfinterp(i,:)*wp%dfgen) 01959 01960 wp%Hm0interp(i) = 4*sqrt(sum(wp%Sfinterp(i,:))*wp%dfgen) 01961 01962 enddo ! i=1,s%ny+1 01963 01964 end subroutine generate_wave_train_properties_per_offshore_point 01965 01966 ! -------------------------------------------------------------- 01967 ! --------- Calculate Fourier componets for each wave ---------- 01968 ! ---------- train component on the offshore boundary ---------- 01969 subroutine generate_wave_train_Fourier(wp) 01970 01971 #ifdef BUILDXBEACH 01972 use interp 01973 use logging_module 01974 #endif 01975 use math_tools 01976 use wave_boundary_datastore 01977 01978 implicit none 01979 ! input/output 01980 type(waveparamsnew),intent(inout) :: wp 01981 ! internal 01982 integer :: i,ii,tempi 01983 complex(fftkind),dimension(:),allocatable :: tempcmplx 01984 01985 ! Determine indices of wave train components in frequency axis and 01986 ! Fourier transform result 01987 allocate(wp%Findex(wp%K)) 01988 tempi = floor(wp%fgen(1)/wp%dfgen) 01989 do i=1,wp%K 01990 wp%Findex(i) = tempi + i 01991 enddo 01992 01993 ! Allocate Fourier coefficients for each y-position at the offshore boundary and 01994 ! each time step 01995 allocate(wp%CompFn(npb,wp%tslen)) 01996 wp%CompFn=0.d0 01997 allocate(tempcmplx(wp%tslen/2-1)) 01998 #ifdef BUILDXBEACH 01999 call writelog('ls','','Calculating Fourier components') 02000 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02001 #endif 02002 do i=1,wp%K 02003 #ifdef BUILDXBEACH 02004 call progress_indicator(.false.,dble(i)/wp%K*100,5.d0,2.d0) 02005 #endif 02006 do ii=1,npb 02007 ! Determine first half of complex Fourier coefficients of wave train 02008 ! components using random phase and amplitudes from sampled spectrum 02009 ! until Nyquist frequency. The amplitudes are represented in a 02010 ! two-sided spectrum, which results in the factor 1/2. 02011 ! Unroll Fourier components along offshore boundary, assuming all wave trains 02012 ! start at x(1,1), y(1,1). 02013 wp%CompFn(ii,wp%Findex(i)) = wp%A(ii,i)/2*exp(par_compi*wp%phigen(i))* & ! Bas: wp%Findex used in time dimension because t=j*dt in frequency space 02014 exp(-par_compi*wp%kgen(i)* & 02015 (dsin(wp%thetagen(i))*(waveBoundaryParameters%yb(ii)-waveBoundaryParameters%y0) & 02016 +dcos(wp%thetagen(i))*(waveBoundaryParameters%xb(ii)-waveBoundaryParameters%x0)) ) 02017 02018 ! Determine Fourier coefficients beyond Nyquist frequency (equal to 02019 ! coefficients at negative frequency axis) of relevant wave components for 02020 ! first y-coordinate by mirroring 02021 tempcmplx = conjg(wp%CompFn(ii,2:wp%tslen/2)) 02022 tempcmplx = flipiv(tempcmplx,size(tempcmplx)) 02023 wp%CompFn(ii,wp%tslen/2+2:wp%tslen)=tempcmplx 02024 enddo 02025 enddo 02026 02027 ! Free memory 02028 deallocate(tempcmplx) 02029 02030 end subroutine generate_wave_train_Fourier 02031 02032 ! -------------------------------------------------------------- 02033 ! --------- Calculate in which computational wave bin ---------- 02034 ! ------------- each wave train component belongs -------------- 02035 subroutine distribute_wave_train_directions(wp,nspr) 02036 02037 #ifdef BUILDXBEACH 02038 use logging_module 02039 #endif 02040 use wave_boundary_datastore 02041 02042 implicit none 02043 02044 ! input/output 02045 type(waveparamsnew),intent(inout) :: wp 02046 integer,intent(in) :: nspr 02047 ! internal 02048 integer :: i,ii,itheta 02049 real*8,dimension(:,:),allocatable :: binedges 02050 logical :: toosmall,toolarge 02051 real*8 :: lostvar,keptvar,perclost 02052 integer :: lntheta ! local copies that can be changed without 02053 real*8 :: ldtheta ! damage to the rest of the model 02054 real*8 :: dang 02055 02056 ! Calculate the bin edges of all the computational wave bins in the 02057 ! XBeach model (not the input spectrum) 02058 allocate(binedges(lntheta,2)) 02059 do itheta=1,lntheta 02060 binedges(itheta,1) = mod(waveBoundaryParameters%theta(itheta)-ldtheta/2,2*par_pi) 02061 binedges(itheta,2) = mod(waveBoundaryParameters%theta(itheta)+ldtheta/2,2*par_pi) 02062 enddo 02063 02064 ! All generated wave components are in the rang 0<=theta<2pi. 02065 ! We link wave components to a wave direction bin if the direction falls 02066 ! within the bin boundaries. Note the >= and <=, ranther than >= and <. This 02067 ! is not necessarily a problem, but solves having to make an exception for the 02068 ! highest wave direction bin, in which >= and <= should be applicable. 02069 ! In the case of a central bin and a wave direction exactly (!) on the bin 02070 ! interface, the wave will be linked to the first wave bin in the ordering, 02071 ! rather than the higher of the two bins. 02072 ! 02073 ! Initially set WDindex to zero. This marks a wave direction outside the 02074 ! computational wave bins. In case it does not fit in a directional wave 02075 ! bin, it remains zero at the end of the loops. 02076 ! Note: this does not ensure all energy is included in the wave bins, 02077 ! as wave energy may still fall outside the computational domain. 02078 wp%WDindex = 0 02079 do i=1,wp%K 02080 do itheta=1,lntheta 02081 ! special case if this bin spans 0 degrees 02082 if (binedges(i,1)>binedges(i,2)) then 02083 ! Need to check both above and below zero degrees 02084 if((wp%thetagen(i)>=binedges(itheta,1) .and. wp%thetagen(i)<=2*par_px) .or. & 02085 (wp%thetagen(i)>=0.d0 .and. wp%thetagen(i)<=binedges(itheta,2) ) ) then 02086 wp%WDindex(i) = itheta 02087 ! We now have the correct wave bin, move to next wave component K 02088 exit 02089 endif 02090 else 02091 if(wp%thetagen(i)>=binedges(itheta,1) .and. wp%thetagen(i)<=binedges(itheta,2)) then 02092 wp%WDindex(i) = itheta 02093 ! We now have the correct wave bin, move to next wave component K 02094 exit 02095 endif 02096 endif 02097 enddo 02098 enddo 02099 02100 ! If the user has set nspr == 1 then the randomly drawn wave directions 02101 ! should be set to the centres of the wave directional bins. 02102 ! Also move all wave energy falling outside the computational bins, into 02103 ! the computational domain (in the outer wave direction bins) 02104 if (nspr==1) then 02105 do i=1,wp%K 02106 if (wp%WDindex(i)>0) then 02107 ! reset the direction of this wave train to the centre of the bin 02108 wp%thetagen(i)=waveBoundaryParameters%theta(wp%WDindex(i)) 02109 endif 02110 enddo 02111 endif 02112 02113 ! Check the amount of energy lost to wave trains falling outside the computational 02114 ! domain 02115 lostvar = 0.d0 02116 keptvar = 0.d0 02117 do i=1,wp%K 02118 if (wp%WDindex(i)==0) then 02119 lostvar = lostvar + sum(wp%A(:,i)**2) 02120 else 02121 keptvar = keptvar + sum(wp%A(:,i)**2) 02122 endif 02123 enddo 02124 perclost = 100*(lostvar/(lostvar+keptvar)) 02125 #ifdef BUILDXBEACH 02126 if (perclost>5.0d0) then 02127 call writelog('lsw','(a,f0.1,a)','Large amounts of energy (',perclost, & 02128 '%) fall outside computational domain at the offshore boundary') 02129 call writelog('lsw','','Check specification of input wave angles and wave directional grid') 02130 else 02131 call writelog('ls','(a,f0.1,a)','Wave energy outside computational domain at offshore boundary: ',perclost,'%') 02132 endif 02133 #endif 02134 ! Free memory 02135 deallocate(binedges) 02136 02137 end subroutine distribute_wave_train_directions 02138 02139 ! -------------------------------------------------------------- 02140 ! --------- Calculate energy envelope time series from --------- 02141 ! -------- Fourier components, and write to output file -------- 02142 subroutine generate_ebcf(wp) 02143 02144 02145 #ifdef BUILDXBEACH 02146 use interp 02147 use logging_module 02148 #endif 02149 use math_tools 02150 02151 implicit none 02152 02153 ! input/output 02154 type(waveparamsnew),intent(inout) :: wp 02155 ! internal 02156 integer :: itheta,iy,iwc,it,irec 02157 integer :: index,status 02158 integer :: reclen,fid 02159 integer,dimension(:),allocatable :: nwc 02160 real*8,dimension(:,:,:), allocatable :: zeta, Ampzeta, E_tdir, E_interp 02161 real*8,dimension(:,:), allocatable :: eta, Amp 02162 complex(fftkind),dimension(:),allocatable :: Gn, tempcmplx,tempcmplxhalf 02163 integer,dimension(:),allocatable :: tempindex,tempinclude 02164 real*8 :: stdzeta,stdeta,etot,perc 02165 02166 02167 ! Allocate variables for water level exitation and amplitude with and without 02168 ! directional spreading dependent envelope 02169 allocate(zeta(npb,wp%tslen,ntheta)) 02170 allocate(Ampzeta(npb,wp%tslen,ntheta)) 02171 zeta=0.d0 02172 Ampzeta=0.d0 02173 02174 allocate(eta(npb,wp%tslen)) 02175 allocate(Amp(npb,wp%tslen)) 02176 eta=0.d0 02177 Amp=0.d0 02178 02179 ! Calculate wave energy for each y-coordinate along seaside boundary for 02180 ! current computational directional bin 02181 allocate(nwc(ntheta)) 02182 allocate(Gn(wp%tslen)) 02183 allocate(tempinclude(wp%K)) 02184 allocate(tempcmplx(wp%tslen)) 02185 allocate(tempcmplxhalf(size(Gn(wp%tslen/2+2:wp%tslen)))) 02186 do itheta=1,ntheta 02187 #ifdef BUILDXBEACH 02188 call writelog('ls','(A,I0,A,I0)','Calculating short wave time series for theta bin ',itheta,' of ',s%ntheta) 02189 #endif 02190 ! Select wave components that are in the current computational 02191 ! directional bin 02192 tempinclude=0 02193 where (wp%WDindex==itheta) 02194 tempinclude=1 02195 end where 02196 02197 ! Determine number of wave components that are in the current 02198 ! computational directional bin 02199 nwc(itheta)=sum(tempinclude) 02200 02201 ! Determine for each wave component in the current computational 02202 ! directional bin its index in the Fourier coefficients array 02203 ! ordered from high to low frequency 02204 allocate(tempindex(nwc(itheta))) 02205 tempindex=0 02206 do iwc=1,nwc(itheta) 02207 ! find highest 02208 index=maxval(maxloc(tempinclude)) 02209 ! reset that one so that the next highest in found in next iteration 02210 tempinclude(index)=0 02211 tempindex(iwc)=wp%Findex(index) 02212 end do 02213 02214 ! Check whether any wave components are in the current computational 02215 ! directional bin 02216 if (nwc(itheta)>0) then 02217 do iy=1,npb 02218 ! Reset 02219 Gn=0 02220 02221 ! Determine Fourier coefficients of all wave components for current 02222 ! y-coordinate in the current computational directional bin 02223 Gn(tempindex)=wp%CompFn(iy,tempindex) 02224 tempcmplxhalf = conjg(Gn(2:wp%tslen/2)) 02225 call flipiv(tempcmplxhalf,size(tempcmplxhalf)) 02226 Gn(wp%tslen/2+2:wp%tslen)=tempcmplxhalf 02227 02228 ! Inverse Discrete Fourier transformation to transform back to time 02229 ! domain from frequency domain 02230 tempcmplx=Gn 02231 status=0 02232 tempcmplx=fft(tempcmplx,inv=.true.,stat=status) 02233 02234 ! Scale result 02235 tempcmplx=tempcmplx/sqrt(dble(size(tempcmplx))) 02236 02237 ! Superimpose gradual increase and decrease of energy input for 02238 ! current y-coordinate and computational diretional bin on 02239 ! instantaneous water level excitation 02240 ! 02241 ! Robert: use final wave elevation from last iteration to startup 02242 ! this boundary condition 02243 zeta(iy,:,itheta)=dble(tempcmplx*wp%tslen) 02244 ! The first time this function is called in a simulation, lastwaveelevation is unknown, 02245 ! so would be set to zero. However, this artificially increases the taper time, and is 02246 ! not useful if repeatwbc = .true., as this zero level is repeated every boundary condition 02247 ! Instead, we select zeta at the end of the first boundary condition time series to 02248 ! initialise lastwaveelevation 02249 if (bccount==0) then 02250 lastwaveelevation(:,itheta) = zeta(iy,ind_end_taper,itheta) 02251 endif 02252 zeta(iy,:,itheta)=zeta(iy,:,itheta)*wp%taperw+(1-wp%taperw)*lastwaveelevation(iy,itheta) 02253 ! note that taperw is only zero at the start, not the end of the time series, so we 02254 ! can safely copy zeta at the point in time where t=rtbc to lastwaveelevation for use 02255 ! in the generation of the next time series 02256 lastwaveelevation(:,itheta) = zeta(iy,ind_end_taper,itheta) 02257 enddo ! iy=1,ny+1 02258 endif ! nwc>0 02259 deallocate(tempindex) 02260 enddo ! itheta = 1,ntheta 02261 deallocate(tempinclude) 02262 deallocate(Gn) 02263 deallocate(tempcmplxhalf) 02264 ! 02265 ! Calculate energy envelope amplitude 02266 02267 ! Robert: note oldwbc has been removed 02268 do iy=1,npb 02269 ! Integrate instantaneous water level excitation of wave 02270 ! components over directions 02271 eta(iy,:) = sum(zeta(iy,:,:),2) 02272 tempcmplx=eta(iy,:) 02273 02274 ! Hilbert transformation to determine envelope of all total 02275 ! non-directional wave components 02276 call hilbert(tempcmplx,size(tempcmplx)) 02277 02278 ! Determine amplitude of water level envelope by calculating 02279 ! the absolute value of the complex wave envelope descriptions 02280 Amp(iy,:)=abs(tempcmplx) 02281 02282 ! Calculate standard deviation of non-directional 02283 ! instantaneous water level excitation of all 02284 ! wave components to be used as weighing factor 02285 stdeta = sqrt(sum(eta(iy,:)**2)/(size(eta(iy,:))-1)) 02286 do itheta=1,ntheta 02287 if (nwc(itheta)>0) then 02288 ! Calculate standard deviations of directional 02289 ! instantaneous water level excitation of all 02290 ! wave components to be used as weighing factor 02291 stdzeta = sqrt(sum(zeta(iy,:,itheta)**2)/(size(zeta(iy,:,itheta))-1)) 02292 02293 ! Calculate amplitude of directional wave envelope 02294 Ampzeta(iy,:,itheta)= Amp(iy,:)*stdzeta/stdeta 02295 else ! nwc==0 02296 ! Current computational directional bin does not contain any wave 02297 ! components 02298 Ampzeta(iy,:,itheta)=0.d0 02299 endif ! nwc>0 02300 end do ! 1:ntheta 02301 ! Print status message to screen 02302 #ifdef BUILDXBEACH 02303 call writelog('ls','(A,I0,A,I0,A)','Y-point ',iy,' of ',s%ny+1,' done.') 02304 #endif 02305 end do ! 1:ny+1 02306 ! 02307 02308 ! free memory 02309 deallocate(tempcmplx) 02310 deallocate(nwc) 02311 02312 ! Allocate memory for energy time series 02313 allocate(E_tdir(npb,wp%tslen,ntheta)) 02314 E_tdir=0.0d0 02315 E_tdir=0.5d0*waveBoundaryParameters%rho*par_g*Ampzeta**2 02316 E_tdir=E_tdir/waveBoundaryParameters%dtheta 02317 02318 ! Print directional energy distribution to screen 02319 etot = sum(E_tdir) 02320 do itheta=1,ntheta 02321 perc = sum(E_tdir(:,:,itheta))/etot*100 02322 #ifdef BUILDXBEACH 02323 call writelog('ls','(a,i0,a,f0.2,a)','Wave bin ',itheta,' contains ',perc,'% of total energy') 02324 #endif 02325 enddo 02326 ! 02327 ! 02328 ! Store time series in internal memory 02329 ! 02330 ! First deallocate arrays if necessary 02331 if (allocated(waveBoundaryTimeSeries%eebct)) deallocate(waveBoundaryTimeSeries%eebct) 02332 if (allocated(waveBoundaryTimeSeries%tbc)) deallocate(waveBoundaryTimeSeries%tbc) 02333 ! 02334 ! Allocate in the correct size 02335 allocate(waveBoundaryTimeSeries%eebct(npb,wp%tslenbc+2,ntheta)) 02336 allocate(waveBoundaryTimeSeries%tbc(wp%tslenbc+2)) 02337 ! 02338 ! Store time vector 02339 do it=2,wp%tslenbc+1 02340 waveBoundaryTimeSeries%tbc(2:wp%tslenbc+1) = (it-2)*wp%dtbc + waveBoundaryAdministration%startCurrentSeries 02341 enddo 02342 ! add 'inifinite ends to the time series in case of mismatch at the end or start of the time 02343 ! series generation and interpolation at each time step 02344 waveBoundaryTimeSeries%tbc(1) = -1.d0*huge(0.d0) 02345 waveBoundaryTimeSeries%tbc(wp%tslenbc+2) = 1.d0*huge(0.d0) 02346 if (wp%dtchanged) then 02347 ! Interpolate from internal time axis to output time axis 02348 do itheta=1,ntheta 02349 do it=2,wp%tslenbc+1 02350 do iy=1,npb 02351 call linear_interp(wp%tin,E_tdir(iy,:,itheta),wp%tslen, & 02352 (it-1)*wp%dtbc,waveBoundaryTimeSeries%eebct(iy,it,itheta),status) 02353 enddo 02354 enddo 02355 enddo 02356 else 02357 ! no need for interpolation 02358 waveBoundaryTimeSeries%eebct(:,2:wp%tslenbc+1,:) = E_tdir 02359 endif 02360 ! 02361 ! add 'inifinite ends to the time series in case of mismatch at the end or start of the time 02362 ! series generation and interpolation at each time step 02363 waveBoundaryTimeSeries%eebct(:,1,:) = waveBoundaryTimeSeries%eebct(:,2,:) 02364 waveBoundaryTimeSeries%eebct(:,wp%tslenbc+2,:) = waveBoundaryTimeSeries%eebct(:,wp%tslenbc+1,:) 02365 ! 02366 ! 02367 ! Free memory 02368 deallocate(zeta,Ampzeta,E_tdir, Amp, eta) 02369 02370 end subroutine generate_ebcf 02371 02372 02373 ! -------------------------------------------------------------- 02374 ! ------ Calculate time series of short wave at offshore ------- 02375 ! --------------------------- boundary ------------------------- 02376 subroutine generate_swts(wp) 02377 #ifdef BUILDXBEACH 02378 use logging_module 02379 #endif 02380 use wave_boundary_datastore 02381 02382 implicit none 02383 02384 ! input/output 02385 type(waveparamsnew),intent(inout) :: wp 02386 ! internal 02387 integer :: j,it,ik 02388 real*8,dimension(npb) :: distx,disty 02389 02390 ! allocate memory for time series of data 02391 allocate(wp%zsits(npb,wp%tslen)) 02392 allocate(wp%uits(npb,wp%tslen)) 02393 wp%zsits=0.d0 02394 wp%uits=0.d0 02395 ! distance of each grid point to reference point 02396 do j=1,npb 02397 distx(j) = waveBoundaryParameters%xb(j)-waveBoundaryParameters%x0 02398 disty(j) = waveBoundaryParameters%yb(j)-waveBoundaryParameters%y0 02399 enddo 02400 ! 02401 ! total surface elevation 02402 #ifdef BUILDXBEACH 02403 call writelog('ls','','Calculating short wave elevation time series') 02404 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02405 #endif 02406 do it=1,wp%tslen 02407 #ifdef BUILDXBEACH 02408 call progress_indicator(.false.,dble(it)/wp%tslen*100,5.d0,2.d0) 02409 #endif 02410 do ik=1,wp%K 02411 do j=1,npb 02412 wp%zsits(j,it)=wp%zsits(j,it)+wp%A(j,ik)*dsin( & 02413 +wp%wgen(ik)*wp%tin(it)& 02414 -wp%kgen(ik)*( dsin(wp%thetagen(ik))*disty(j) & 02415 +dcos(wp%thetagen(ik))*distx(j) & 02416 ) & 02417 +wp%phigen(ik) & 02418 ) 02419 enddo 02420 enddo 02421 enddo 02422 ! depth-averaged velocity 02423 #ifdef BUILDXBEACH 02424 call writelog('ls','','Calculating short wave velocity time series') 02425 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02426 #endif 02427 do it=1,wp%tslen 02428 #ifdef BUILDXBEACH 02429 call progress_indicator(.false.,dble(it)/wp%tslen*100,5.d0,2.d0) 02430 #endif 02431 do ik=1,wp%K 02432 do j=1,npb 02433 wp%uits(j,it) = wp%uits(j,it) + & 02434 1.d0/hb0*wp%wgen(ik)*wp%A(j,ik)/sinh(wp%kgen(ik)*hb0) * & 02435 dsin( & 02436 +wp%wgen(ik)*wp%tin(it)& 02437 -wp%kgen(ik)*( dsin(wp%thetagen(ik))*disty(j) & 02438 +dcos(wp%thetagen(ik))*distx(j) & 02439 ) & 02440 +wp%phigen(ik) & 02441 ) * & 02442 1.d0/wp%kgen(ik)*sinh(wp%kgen(ik)*hb0) 02443 02444 enddo 02445 enddo 02446 enddo 02447 ! 02448 ! 02449 ! Apply tapering to time series 02450 do j=1,npb 02451 wp%uits(j,:)=wp%uits(j,:)*wp%taperf 02452 wp%zsits(j,:)=wp%zsits(j,:)*wp%taperf 02453 enddo 02454 end subroutine generate_swts 02455 02456 02457 ! -------------------------------------------------------------- 02458 ! ----------------------- Bound long wave ---------------------- 02459 ! -------------------------------------------------------------- 02460 subroutine generate_qbcf(wp) 02461 02462 #ifdef BUILDXBEACH 02463 use logging_module 02464 use interp 02465 #endif 02466 use math_tools 02467 use wave_boundary_datastore 02468 02469 implicit none 02470 02471 ! input/output 02472 type(waveparamsnew),intent(inout) :: wp 02473 ! internal 02474 integer :: j,m,iq,irec,it ! counters 02475 integer :: K ! copy of K 02476 integer :: halflen,reclen,fid,status 02477 logical :: firsttime ! used for output message only 02478 real*8 :: deltaf ! difference frequency 02479 real*8,dimension(:), allocatable :: term1,term2,term2new,dif,chk1,chk2 02480 real*8,dimension(npb) :: distx,disty 02481 real*8,dimension(:,:),allocatable :: Eforc,D,deltheta,KKx,KKy,dphi3,k3,cg3,theta3,Abnd 02482 real*8,dimension(:,:,:),allocatable :: q,qinterp 02483 complex(fftkind),dimension(:),allocatable :: Comptemp,Comptemp2,Gn 02484 complex(fftkind),dimension(:,:,:),allocatable:: Ftemp 02485 02486 02487 02488 ! This function has changed with respect to previous versions of XBeach, in that 02489 ! the bound long wave has to be calculated separately at each longshore point, 02490 ! owing to longshore varying incident wave spectra 02491 02492 ! shortcut variable 02493 K = wp%K 02494 02495 ! Print message to screen 02496 #ifdef BUILDXBEACH 02497 call writelog('sl','', 'Calculating primary wave interaction') 02498 #endif 02499 02500 ! Allocate two-dimensional variables for all combinations of interacting wave 02501 ! components to be filled triangular 02502 allocate(Eforc(K-1,K)) 02503 allocate(D(K-1,K)) 02504 allocate(deltheta(K-1,K)) 02505 allocate(KKx(K-1,K)) 02506 allocate(KKy(K-1,K)) 02507 allocate(dphi3(K-1,K)) 02508 allocate(k3(K-1,K)) 02509 allocate(cg3(K-1,K)) 02510 allocate(theta3(K-1,K)) 02511 ! Allocate variables for amplitude and Fourier coefficients of bound long wave 02512 allocate(Gn(wp%tslen)) 02513 allocate(Abnd(K-1,K)) 02514 allocate(Ftemp(K-1,K,3)) ! Jaap qx, qy, zeta 02515 ! Storage for output discharge 02516 allocate(q(s%ny+1,wp%tslen,4)) ! qx qy qtot, zeta 02517 ! 02518 ! Initialize variables as zero 02519 Eforc = 0 02520 D = 0 02521 deltheta = 0 02522 KKx = 0 02523 KKy = 0 02524 dphi3 = 0 02525 k3 = 0 02526 cg3 = 0 02527 theta3 = 0 02528 Gn = 0*par_compi 02529 Abnd = 0 02530 Ftemp = 0*par_compi 02531 q = 0 02532 02533 ! First time is set true for each time new wave bc are generated 02534 firsttime = .true. 02535 02536 ! upper half of frequency space 02537 halflen = wp%tslen/2 02538 02539 ! Run loop over wave-wave interaction components 02540 #ifdef BUILDXBEACH 02541 call progress_indicator(.true.,0.d0,5.d0,2.d0) 02542 #endif 02543 do m=1,K-1 02544 #ifdef BUILDXBEACH 02545 call progress_indicator(.false.,dble(m)/(K-1)*100,5.d0,2.d0) 02546 #endif 02547 ! Allocate memory 02548 allocate(term1(K-m),term2(K-m),term2new(K-m),dif(K-m),chk1(K-m),chk2(K-m)) 02549 02550 ! Determine difference frequency 02551 deltaf=m*wp%dfgen 02552 02553 ! Determine difference angles (pi already added) 02554 deltheta(m,1:K-m) = abs(wp%thetagen(m+1:K)-wp%thetagen(1:K-m))+par_pi 02555 02556 ! Determine x- and y-components of wave numbers of difference waves 02557 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)) 02558 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)) 02559 02560 ! Determine difference wave numbers according to Van Dongeren et al. 2003 02561 ! eq. 19 02562 k3(m,1:K-m) =sqrt(wp%kgen(1:K-m)**2+wp%kgen(m+1:K)**2+ & 02563 2*wp%kgen(1:K-m)*wp%kgen(m+1:K)*dcos(deltheta(m,1:K-m))) 02564 02565 ! Determine group velocity of difference waves 02566 cg3(m,1:K-m)= 2.d0*par_pi*deltaf/k3(m,1:K-m) 02567 02568 ! Modification Robert + Jaap: make sure that the bound long wave amplitude does not 02569 ! explode when offshore boundary is too close to shore, 02570 ! by limiting the interaction group velocity 02571 cg3(m,1:K-m) = min(cg3(m,1:K-m),waveBoundaryParameters%nmax*sqrt(par_g/k3(m,1:K-m)*tanh(k3(m,1:K-m)*hb0))) 02572 02573 ! Determine difference-interaction coefficient according to Herbers 1994 02574 ! eq. A5 02575 term1 = (-wp%wgen(1:K-m))*wp%wgen(m+1:K) 02576 term2 = (-wp%wgen(1:K-m))+wp%wgen(m+1:K) 02577 term2new = cg3(m,1:K-m)*k3(m,1:K-m) 02578 dif = (abs(term2-term2new)) 02579 if (any(dif>0.01*term2) .and. firsttime) then 02580 firsttime = .false. 02581 endif 02582 chk1 = cosh(wp%kgen(1:K-m)*hb0) 02583 chk2 = cosh(wp%kgen(m+1:K)*hb0) 02584 02585 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+par_g*term2*(chk1*chk2)/ & 02586 ((par_g*k3(m,1:K-m)*tanh(k3(m,1:K-m)*hb0)-(term2new)**2)*term1*cosh(k3(m,1:K-m)*hb0))* & 02587 (term2*((term1)**2/par_g/par_g - wp%kgen(1:K-m)*wp%kgen(m+1:K)*dcos(deltheta(m,1:K-m))) & 02588 - 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))) 02589 02590 ! Correct for surface elevation input and output instead of bottom pressure 02591 ! so it is consistent with Van Dongeren et al 2003 eq. 18 02592 D(m,1:K-m) = D(m,1:K-m)*cosh(k3(m,1:K-m)*hb0)/(cosh(wp%kgen(1:K-m)*hb0)*cosh(wp%kgen(m+1:K)*hb0)) 02593 02594 ! Exclude interactions with components smaller than or equal to current 02595 ! component according to lower limit Herbers 1994 eq. 1 02596 where(wp%fgen<=deltaf) D(m,:)=0.d0 02597 02598 ! Exclude interactions with components that are cut-off by the fcutoff 02599 ! parameter 02600 if (deltaf<=waveBoundaryParameters%fcutoff) D(m,:)=0.d0 02601 02602 ! Determine phase of bound long wave assuming a local equilibrium with 02603 ! forcing of interacting primary waves according to Van Dongeren et al. 02604 ! 2003 eq. 21 (the angle is the imaginary part of the natural log of a 02605 ! complex number as long as the complex number is not zero) 02606 allocate(Comptemp(K-m),Comptemp2(K-m)) 02607 Comptemp=conjg(wp%CompFn(1,wp%Findex(1)+m:wp%Findex(1)+K-1)) 02608 Comptemp2=conjg(wp%CompFn(1,wp%Findex(1):wp%Findex(1)+K-m-1)) 02609 dphi3(m,1:K-m) = par_pi+imag(log(Comptemp))-imag(log(Comptemp2)) 02610 deallocate (Comptemp,Comptemp2) 02611 ! 02612 ! Determine angle of bound long wave according to Van Dongeren et al. 2003 eq. 22 02613 theta3 = atan2(KKy,KKx) 02614 ! 02615 ! free memory 02616 deallocate(term1,term2,term2new,dif,chk1,chk2) 02617 enddo ! m=1,K-1 02618 ! 02619 ! Output to screen 02620 #ifdef BUILDXBEACH 02621 if (.not. firsttime) then 02622 call writelog('lws','','Warning: shallow water so long wave variance is reduced using par%nmax') 02623 endif 02624 call writelog('sl','', 'Calculating flux at boundary') 02625 #endif 02626 ! 02627 ! Allocate temporary arrays for upcoming loop 02628 allocate(Comptemp(halflen-1)) 02629 allocate(Comptemp2(wp%tslen)) 02630 ! 02631 ! 02632 ! distance of each grid point to reference point 02633 do j=1,npb 02634 distx(j) = waveBoundaryParameters%xb(j)-waveBoundaryParameters%x0 02635 disty(j) = waveBoundaryParameters%yb(j)-waveBoundaryParameters%y0 02636 enddo 02637 ! 02638 ! Run a loop over the offshore boundary 02639 do j=1,npb 02640 ! Determine energy of bound long wave according to Herbers 1994 eq. 1 based 02641 ! on difference-interaction coefficient and energy density spectra of 02642 ! primary waves 02643 ! Robert: E = 2*D**2*S**2*dtheta**2*df can be rewritten as 02644 ! E = 2*D**2*Sf**2*df 02645 Eforc = 0 02646 do m=1,K-1 02647 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 02648 enddo 02649 ! 02650 ! Calculate bound wave amplitude for this offshore grid point 02651 Abnd = sqrt(2*Eforc*wp%dfgen) 02652 ! 02653 ! Determine complex description of bound long wave per interaction pair of 02654 ! primary waves for first y-coordinate along seaside boundary 02655 Ftemp(:,:,1) = Abnd/2*exp(-1*par_compi*dphi3)*cg3*dcos(theta3) ! qx 02656 Ftemp(:,:,2) = Abnd/2*exp(-1*par_compi*dphi3)*cg3*dsin(theta3) ! qy 02657 !Ftemp(:,:,3) = Abnd/2*exp(-1*par_compi*dphi3)*cg3 ! qtot 02658 Ftemp(:,:,3) = Abnd/2*exp(-1*par_compi*dphi3) ! eta 02659 ! 02660 ! loop over qx,qy and qtot 02661 do iq=1,3 02662 ! Unroll wave component to correct place along the offshore boundary 02663 Ftemp(:,:,iq) = Ftemp(:,:,iq)* & 02664 exp(-1*par_compi*(KKy*disty+KKx*distx)) 02665 ! Determine Fourier coefficients 02666 Gn(2:K) = sum(Ftemp(:,:,iq),DIM=2) 02667 Comptemp = conjg(Gn(2:halflen)) 02668 call flipiv(Comptemp,halflen-1) 02669 Gn(halflen+2:wp%tslen) = Comptemp 02670 ! 02671 ! Print status message to screen 02672 #ifdef BUILDXBEACH 02673 if (iq==3) then 02674 call writelog('ls','(A,I0,A,I0)','Flux ',j,' of ',s%ny+1) 02675 endif 02676 #endif 02677 ! 02678 ! Inverse Discrete Fourier transformation to transform back to time space 02679 ! from frequency space 02680 Comptemp2=fft(Gn,inv=.true.) 02681 ! 02682 ! Determine mass flux as function of time and let the flux gradually 02683 ! increase and decrease in and out the wave time record using the earlier 02684 ! specified window 02685 Comptemp2=Comptemp2/sqrt(dble(wp%tslen)) 02686 q(j,:,iq)=dreal(Comptemp2*wp%tslen)*wp%taperf 02687 enddo ! iq=1,3 02688 enddo ! j=1,s%ny+1 02689 ! 02690 ! free memory 02691 deallocate(Comptemp,Comptemp2,Ftemp) 02692 ! 02693 if (.not. waveBoundaryParameters%nonhspectrum) then 02694 if(allocated(waveBoundaryTimeSeries%qxbct)) deallocate(waveBoundaryTimeSeries%qxbct) 02695 if(allocated(waveBoundaryTimeSeries%qybct)) deallocate(waveBoundaryTimeSeries%qybct) 02696 allocate(waveBoundaryTimeSeries%qxbct(npb,tslenbc+2)) 02697 allocate(waveBoundaryTimeSeries%qxbct(npb,tslenbc+2)) 02698 if (wp%dtchanged) then 02699 ! Interpolate from internal time axis to output time axis 02700 do it=2,wp%tslenbc+1 02701 do j=1,npb 02702 call linear_interp(wp%tin,q(j,:,1),wp%tslen, & 02703 (it-1)*wp%dtbc,waveBoundaryTimeSeries%qxbct(j,it),status) 02704 call linear_interp(wp%tin,q(j,:,2),wp%tslen, & 02705 (it-1)*wp%dtbc,waveBoundaryTimeSeries%qybct(j,it),status) 02706 enddo 02707 enddo 02708 else 02709 ! no need for interpolation 02710 waveBoundaryTimeSeries%qxbct(:,2:wp%tslenbc+1) = q(:,:,1) 02711 waveBoundaryTimeSeries%qybct(:,2:wp%tslenbc+1) = q(:,:,2) 02712 endif 02713 waveBoundaryTimeSeries%qxbct(:,1) = waveBoundaryTimeSeries%qxbct(:,2) 02714 waveBoundaryTimeSeries%qxbct(:,wp%tslenbc+2) = waveBoundaryTimeSeries%qxbct(:,wp%tslenbc+1) 02715 waveBoundaryTimeSeries%qybct(:,1) = waveBoundaryTimeSeries%qybct(:,2) 02716 waveBoundaryTimeSeries%qybct(:,wp%tslenbc+2) = waveBoundaryTimeSeries%qybct(:,wp%tslenbc+1) 02717 else 02718 do j=1,npb 02719 ! add to velocity time series 02720 wp%uits(j,:)=wp%uits(j,:)+q(j,:,1)/hb0 02721 ! add to surface elevation time series 02722 wp%zsits(j,:)=wp%zsits(j,:)+q(j,:,3) 02723 enddo 02724 endif ! par%nonhspectrum==1 02725 02726 ! Free memory 02727 deallocate(Eforc,D,deltheta,KKx,KKy,dphi3,k3,cg3,theta3,Gn,Abnd,q) 02728 02729 end subroutine generate_qbcf 02730 02731 02732 ! -------------------------------------------------------------- 02733 ! --------------- Non-hydrostatic wave time -------------------- 02734 ! ---------------- series file generation ---------------------- 02735 subroutine generate_nhtimeseries_file(wp) 02736 02737 #ifdef BUILDXBEACH 02738 use logging_module 02739 #endif 02740 use wave_boundary_datastore 02741 02742 implicit none 02743 02744 ! input/output 02745 type(waveparamsnew),intent(inout) :: wp 02746 ! internal 02747 02748 #ifdef BUILDXBEACH 02749 call writelog('ls','','Writing short wave time series to ',wp%nhfilename) 02750 #endif 02751 if(allocated(waveBoundaryTimeSeries%zsbct)) deallocate(waveBoundaryTimeSeries%zsbct) 02752 if(allocated(waveBoundaryTimeSeries%ubct)) deallocate(waveBoundaryTimeSeries%ubct) 02753 allocate(waveBoundaryTimeSeries%zsbct(npb,wp%tslen+2)) 02754 allocate(waveBoundaryTimeSeries%ubct(npb,wp%tslen+2)) 02755 02756 waveBoundaryTimeSeries%zsbct(:,2:wp%tslen+1) = wp%zsits 02757 waveBoundaryTimeSeries%ubct(:,2:wp%tslen+1) = wp%uits 02758 02759 waveBoundaryTimeSeries%zsbct(:,1) = waveBoundaryTimeSeries%zsbct(:,2) 02760 waveBoundaryTimeSeries%zsbct(:,wp%tslen+2) = waveBoundaryTimeSeries%zsbct(:,wp%tslen+1) 02761 waveBoundaryTimeSeries%ubct(:,1) = waveBoundaryTimeSeries%ubct(:,2) 02762 waveBoundaryTimeSeries%ubct(:,wp%tslen+2) = waveBoundaryTimeSeries%ubct(:,wp%tslen+1) 02763 02764 end subroutine generate_nhtimeseries_file 02765 02766 02767 02768 end module
1.7.6.1
