|
XBeach
|
00001 module wave_boundary_main_module 00002 implicit none 00003 save 00004 private 00005 public create_incident_waves 00006 ! This module is the main entry to generating wave boundary conditions in XBeach 00007 ! and other models. This module is accessed by wave_boundary_update_module. 00008 ! The interface 'create_incident_wave', which returns wave boundary conditions at 00009 ! the correct time step can be called from outside the module. 00010 ! 00011 ! NOTE! 00012 ! 'wave_boundary' should only be called by processes with a boundary 00013 ! for which wave boundary conditions are required (in XBeach only if xmpi_istop). 00014 ! Else processes will waste I/O resources and computational time generating 00015 ! useless information. 00016 ! 00017 ! TO FIX: 00018 ! 00019 ! - generate_wave_train_properties_per_offshore_point 00020 ! - read all spectrum files 00021 ! - continue at line 2117 wave_boundary_update 00022 ! 00023 ! 00024 ! 00025 ! generate an interface so we don't have to pass unnecessary vectors 00026 ! when using different types of boundary conditions 00027 ! 00028 ! 00029 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00030 ! ! 00031 ! INTERFACE TO OTHER MODELS ! 00032 ! ! 00033 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00034 ! 00035 interface create_incident_waves 00036 module procedure create_incident_waves_surfbeat 00037 !module procedure generate_wave_boundary_nonhydrostatic 00038 end interface create_incident_waves 00039 ! 00040 ! 00041 ! 00042 contains 00043 ! 00044 ! 00045 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00046 ! ! 00047 ! SUBROUTINES CALLED BY INTERFACE ! 00048 ! ! 00049 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00050 ! 00051 ! 00052 subroutine create_incident_waves_surfbeat(np,xb,yb,ntheta,dtheta,theta,t, & 00053 bctype,bcfile, & 00054 x0,y0,hboundary, & 00055 randomseed, & 00056 eebc,qsbc,qnbc, & 00057 Hbc,Tbc,Dbc,isRecomputed, & 00058 nonhspectrum, & 00059 sprdthr,trepfac,nmax,fcutoff,rho, & 00060 Tm01switch,nspr) 00061 ! This subroutine handles all calls for surf-beat wave boundary conditions. 00062 ! The subroutine automatically initialises all variables if needed, and returns 00063 ! boundary condition information at all offshore points at the required point 00064 ! in time. 00065 ! 00066 ! Input variables 00067 ! np : number of offshore grid points (-) 00068 ! xb,yb : vectors of x and y coordinates of offshore grid points (m) 00069 ! ntheta : number of computational wave bins in directional space (-) 00070 ! dtheta : (constant) grid size of wave direction bins (rad) 00071 ! theta : vector of centre of wave direction bins (rad). Angles are in 00072 ! cartesian system relative to the coordinate system x (UTM East) axis 00073 ! t : current time (s) 00074 ! bctype : integer specifying the type of boundary conditions to produce, 00075 ! equal to par%instat in XBeach (-), see paramsconst.F90 for 00076 ! options 00077 ! bcfile : name of main wave boundary condition file to read, equal to 00078 ! : par%bcfile in XBeach (-) 00079 ! x0,y0 : reference point coordinates for wave conditions (for instance 00080 ! origin of model grid). Used to determine phase of wave 00081 ! components in boundary conditions. Must be identical on all 00082 ! processes (m) 00083 ! hboundary : average water depth along the entire offshore boundary (m). This 00084 ! value should be constant, or not vary strongly, across all 00085 ! processes. 00086 ! randomseed : integer containing initial random seed value between 1 and 2^31-2. 00087 ! This should have an identical value on all processes, and is used 00088 ! to generate identical random numbers sequences on all processors (-). 00089 ! If set to same integer throughout simulation, identical random numbers 00090 ! will be generated. Use int(clocktime) to have new random numbers for 00091 ! each time series generation. 00092 ! 00093 ! Output variables 00094 ! eebc : array of size (np,ntheta) containing wave energy density per 00095 ! offshore point and wave direction at time=t (J/m2/rad) 00096 ! qsbc,qnbc : vector of size (np) containing depth-avaraged discharge per along-boundary 00097 ! meter in the s (landward) and n (longshore) direction per 00098 ! offshore point at time=t (m2/s) 00099 ! Hbc : Hm0 computed for the entire boundary based on the input spectra, valid for 00100 ! the duration of the entire spectrum (m) 00101 ! Tbc : Trep computed for the entire boundary based on the input spectra, valid for 00102 ! the duration of the entire spectrum (s) 00103 ! Dbc : mean wave direction computed for the entire boundary based on the input spectra, 00104 ! valid for the duration of the entire spectrum (rad) 00105 ! isRecomputed : logical, indicating whether a new spectrum has been read and computed and therfore 00106 ! showing new values for Hbc, Tbc and Dbc 00107 ! 00108 ! Optional input variables 00109 ! nonhspectrum : generate a non-hydrostatic time series instead of a surf-beat time 00110 ! series. Default = .false. 00111 ! sprdthr : Threshold ratio to maxval of S above which spec dens are read in, see XBeach. 00112 ! Default = 0.08 00113 ! trepfac : Compute mean wave period over energy band: par%trepfac*maxval(Sf); converges to Tm01 00114 ! for trepfac = 0.0. Default = 0.01 00115 ! nmax : maximum ratio of cg/c fro computing long wave boundary conditions. Default = 0.8d0 00116 ! fcutoff : Low-frequency cutoff frequency for long-wave interaction components. Default = 0.d0 00117 ! rho : Density of water (kg/m3). Default = 1025.0 00118 ! Tm01switch : Turn on Tm01 (1) or Tm-10 (0) to compute Trep. Default = 0 00119 ! nspr : nspr = 1 long wave direction forced into centres of short wave bins, 00120 ! nspr = 0 regular long wave spreading. Default = 0 00121 use wave_boundary_datastore 00122 use wave_boundary_update_module, only: generate_wave_boundary_surfbeat 00123 #ifdef BUILDXBEACH 00124 use interp 00125 #endif 00126 ! 00127 implicit none 00128 ! 00129 ! Input variables 00130 integer,intent(in) :: np,ntheta,bctype 00131 real*8,intent(in) :: t,x0,y0,hboundary,dtheta 00132 character(len=*),intent(in) :: bcfile 00133 real*8,dimension(np),intent(in) :: xb,yb 00134 real*8,dimension(ntheta),intent(in) :: theta 00135 integer,intent(in) :: randomseed 00136 ! output variables 00137 real*8,intent(out) :: Hbc,Tbc,Dbc 00138 logical,intent(out) :: isRecomputed 00139 real*8,dimension(np),intent(out) :: qsbc,qnbc 00140 real*8,dimension(np,ntheta),intent(out):: eebc 00141 ! Optional variables 00142 logical,optional,intent(in) :: nonhspectrum 00143 real*8,optional,intent(in) :: sprdthr,trepfac,nmax,rho,fcutoff 00144 integer,optional,intent(in) :: Tm01switch 00145 ! internal variables 00146 integer :: i,itheta,l,dummy 00147 real*8 :: durationlength 00148 ! 00149 ! 00150 ! Check function input arguments and set defaults 00151 if(.not.present(nonhspectrum)) then 00152 waveBoundaryParameters%nonhspectrum = .false. 00153 else 00154 waveBoundaryParameters%nonhspectrum = nonhspectrum 00155 endif 00156 if(.not.present(sprdthr)) then 00157 waveBoundaryParameters%sprdthr = 0.08d0 00158 else 00159 waveBoundaryParameters%sprdthr = sprdthr 00160 endif 00161 if(.not.present(trepfac)) then 00162 waveBoundaryParameters%trepfac = 0.01d0 00163 else 00164 waveBoundaryParameters%trepfac = trepfac 00165 endif 00166 if(.not.present(Tm01switch)) then 00167 waveBoundaryParameters%Tm01switch = 0 00168 else 00169 waveBoundaryParameters%Tm01switch = Tm01switch 00170 endif 00171 if(.not.present(nspr)) then 00172 waveBoundaryParameters%nspr = 0 00173 else 00174 waveBoundaryParameters%nspr = nspr 00175 endif 00176 if(.not.present(nmax)) then 00177 waveBoundaryParameters%nmax = 0.8d0 00178 else 00179 waveBoundaryParameters%nmax = nmax 00180 endif 00181 if(.not.present(fcutoff)) then 00182 waveBoundaryParameters%fcutoff = 0.0d0 00183 else 00184 waveBoundaryParameters%fcutoff = fcutoff 00185 endif 00186 if(.not.present(rho)) then 00187 waveBoundaryParameters%rho = 1025.d0 00188 else 00189 waveBoundaryParameters%rho = rho 00190 endif 00191 ! 00192 ! 00193 ! Check if the wave boundary conditions have been initialised 00194 if (.not.waveBoundaryAdministration%initialized) then 00195 ! 00196 ! 00197 ! Allocate copies of main grid properties in this module. Note, 00198 ! this may be possible with pointers instead of allocated arrays 00199 waveBoundaryParameters%masterFileName = bcfile 00200 waveBoundaryParameters%np = np 00201 waveBoundaryParameters%ntheta = ntheta 00202 waveBoundaryParameters%x0 = x0 00203 waveBoundaryParameters%y0 = y0 00204 waveBoundaryParameters%hboundary = hboundary 00205 ! Potentially, initialized can be set outside this module, so ensure all these 00206 ! arrays are deallocated if necessary 00207 if(allocated(waveBoundaryParameters%xb)) deallocate(waveBoundaryParameters%xb) 00208 if(allocated(waveBoundaryParameters%yb)) deallocate(waveBoundaryParameters%yb) 00209 if(allocated(waveBoundaryParameters%theta)) deallocate(waveBoundaryParameters%theta) 00210 ! Now allocate arrays to the correct size and set values 00211 allocate(waveBoundaryParameters%xb(np)) 00212 allocate(waveBoundaryParameters%yb(np)) 00213 allocate(waveBoundaryParameters%theta(ntheta)) 00214 waveBoundaryParameters%xb = xb 00215 waveBoundaryParameters%yb = yb 00216 waveBoundaryParameters%theta = theta 00217 ! Ensure all theta directions are between 0 and 2pi, required for some trig. on some compilers 00218 do itheta=1,ntheta 00219 waveBoundaryParameters%theta(itheta) = mod(waveBoundaryParameters%theta(itheta),8.d0*atan(1.d0)) 00220 enddo 00221 ! Allocate space for the random seed. This seed is set to 40 integers and 00222 ! should be identical on all processes 00223 if(allocated(waveBoundaryParameters%randomseed)) deallocate(waveBoundaryParameters%randomseed) 00224 waveBoundaryParameters%randomseed = randomseed 00225 ! 00226 ! 00227 ! Start initialization subroutines 00228 call initialise_wave_spectrum_parameters 00229 waveBoundaryAdministration%initialized = .true. 00230 ! 00231 ! 00232 ! Set time to recompute new boundary condition time series to 00233 ! now so boundary conditions are generated in first time step 00234 waveBoundaryAdministration%startComputeNewSeries = t 00235 endif ! initialized 00236 ! 00237 ! 00238 ! Generate or interpolate boundary condition time series 00239 if (t>=waveBoundaryAdministration%startComputeNewSeries) then 00240 ! The start of the current boundary condition should be the end of the previous 00241 ! boundary condition 00242 waveBoundaryAdministration%startCurrentSeries = waveBoundaryAdministration%startComputeNewSeries 00243 ! Call subroutine to generate wave boundary condition time series from spectral 00244 ! input 00245 call generate_wave_boundary_surfbeat(durationlength) 00246 ! 00247 ! 00248 ! Update time administration 00249 waveBoundaryAdministration%startComputeNewSeries = waveBoundaryAdministration%startComputeNewSeries + & 00250 durationlength 00251 00252 00253 isRecomputed = .true. 00254 endif 00255 00256 ! Interpolate energy and discharge at all locations in time 00257 #ifdef BUILDXBEACH 00258 l = size(waveBoundaryTimeSeries%tbc) 00259 do itheta=1,ntheta 00260 do i=1,np 00261 call linear_interp(waveBoundaryTimeSeries%tbc,waveBoundaryTimeSeries%eebct(i,itheta,:),l,& 00262 t,eebc(i,itheta),dummy) 00263 enddo 00264 enddo 00265 do i=1,np 00266 call linear_interp(waveBoundaryTimeSeries%tbc,waveBoundaryTimeSeries%qxbct(i,:),l,& 00267 t,qsbc(i),dummy) 00268 call linear_interp(waveBoundaryTimeSeries%tbc,waveBoundaryTimeSeries%qybct(i,:),l,& 00269 t,qnbc(i),dummy) 00270 enddo 00271 #else 00272 ! Interpolate time series in other models. 00273 #endif 00274 Hbc = waveSpectrumAdministration%Hbc 00275 Tbc = waveSpectrumAdministration%Tbc 00276 Dbc = waveSpectrumAdministration%Dbc 00277 00278 end subroutine create_incident_waves_surfbeat 00279 00280 subroutine initialise_wave_spectrum_parameters 00281 ! This subroutine reads 00282 use wave_boundary_datastore 00283 #ifdef BUILDXBEACH 00284 use filefunctions 00285 use logging_module 00286 use xmpi_module, only: Halt_Program 00287 #endif 00288 00289 !internal variables 00290 integer :: fid,err 00291 integer :: i,nspectra 00292 character(1024) :: testline 00293 character(1) :: testchar 00294 00295 #ifdef BUILDXBEACH 00296 call writelog('l','','--------------------------------') 00297 call writelog('l','','Initializing spectral wave boundary conditions ') 00298 #endif 00299 ! Initialize that wave boundary conditions need to be calculated (first time at least) 00300 ! Stored and defined in wave_boundary_main_module 00301 waveSpectrumAdministration%repeatwbc = .false. 00302 ! Initialize the number of times wave boundary conditions have been generated. 00303 ! Stored and defined in wave_boundary_main_module 00304 waveSpectrumAdministration%bccount = 0 00305 ! Initialize bcendtime to zero. 00306 ! Stored and defined in wave_boundary_main_module 00307 waveSpectrumAdministration%spectrumendtime = 0.d0 00308 ! Initialise lastwaveheight to zero 00309 ! Stored and defined in wave_boundary_main_module 00310 allocate(waveSpectrumAdministration%lastwaveelevation(waveBoundaryParameters%np,& 00311 waveBoundaryParameters%ntheta)) 00312 ! 00313 ! open location list file 00314 #ifdef BUILDXBEACH 00315 fid = create_new_fid() 00316 #endif 00317 open(fid,file=waveBoundaryParameters%masterFileName,status='old',form='formatted') 00318 ! check for LOCLIST 00319 read(fid,*)testline 00320 if (trim(testline)=='LOCLIST') then 00321 ! check the length of this file 00322 i = 0 00323 do while (err==0) 00324 i=i+1 00325 read(fid,*,IOSTAT=err)testchar 00326 enddo 00327 rewind(fid) 00328 nspectra = i-1 00329 ! store this information in the main module 00330 waveSpectrumAdministration%nspectra = nspectra 00331 ! We need at least 1 spectrum location 00332 if (nspectra<1) then 00333 #ifdef BUILDXBEACH 00334 call writelog('ewls','(a,a)','Error reading file ', & 00335 trim(waveBoundaryParameters%masterFileName)) 00336 call writelog('ewls','(a,a)','Ensure that if LOCLIST option is used, ', & 00337 'at least one spectrum location is specified') 00338 call halt_program 00339 #endif 00340 endif 00341 ! read first line again to set cursor at correct point in the file 00342 read(fid,*,IOSTAT=err)testline 00343 ! 00344 ! 00345 ! Now allocate space for variables 00346 allocate(waveSpectrumAdministration%bcfiles(nspectra)) 00347 allocate(waveSpectrumAdministration%xspec(nspectra)) 00348 allocate(waveSpectrumAdministration%yspec(nspectra)) 00349 do i=1,nspectra 00350 ! read x,y and file name per location 00351 read(fid,*,IOSTAT=err)waveSpectrumAdministration%xspec(i), & 00352 waveSpectrumAdministration%yspec(i), & 00353 waveSpectrumAdministration%bcfiles(i)%fname 00354 waveSpectrumAdministration%bcfiles(i)%listline = 0 00355 if (err /= 0) then 00356 #ifdef BUILDXBEACH 00357 ! something has gone wrong during the read of this file 00358 call writelog('ewls','(a,i0,a,a)','error reading line ',i+1,' of file ', & 00359 trim(waveBoundaryParameters%masterFileName)) 00360 call writelog('ewls','','check file for format errors') 00361 call halt_program 00362 #endif 00363 endif 00364 enddo 00365 else 00366 nspectra = 1 00367 allocate(waveSpectrumAdministration%bcfiles(nspectra)) 00368 allocate(waveSpectrumAdministration%xspec(nspectra)) 00369 allocate(waveSpectrumAdministration%yspec(nspectra)) 00370 waveSpectrumAdministration%bcfiles(1)%fname = waveBoundaryParameters%masterFileName 00371 waveSpectrumAdministration%bcfiles(1)%listline = 0 00372 waveSpectrumAdministration%xspec = waveBoundaryParameters%x0 00373 waveSpectrumAdministration%yspec = waveBoundaryParameters%y0 00374 endif 00375 close(fid) 00376 #ifdef BUILDXBEACH 00377 call writelog('l','','--------------------------------') 00378 #endif 00379 end subroutine initialise_wave_spectrum_parameters 00380 00381 00382 end module wave_boundary_main_module 00383 00384 00385 00386
1.7.6.1
