Skip to content

Commit 701b541

Browse files
authored
82 overcome dependency snapwave ig source term on previous timestep (#83)
* - First version without previous timestep dependency - Note ; had to change definition of gamma, so not anymore mean gamma over current and upwind point > changes gamma and thus source term and IG growth, have to check still whether that is a problem or not - todo: clean up * - Solved so we can use the mean gamma again, issue was het H(k) is set to 0 everytime snapwave_solver is called, so now determined again - Code can be cleaned up by doing a H inout definition probably * - Bit cleanup, and move back 'Make sure DoverE is filled based on previous ee' part - and result now only little bit different anymore, because of removal of previous loop now? * - Hprev should be based on full E, not ee(itheta), now result is same as in paper - Thereby not using H_inc_old and H_ig_old anymore from wave heights of last time step (but implicitly we still do) * - Clean up by changing H to inout, so that we don't have to determine it again at the start of solve_energy_balance2Dstat - Remove Sxx outside of subroutine 'solve_energy_balance2Dstat', since that is not needed anymore * - Change water depth for warning limited depth in Herbers calculation to determine Hig,0 from 10 m to 5 m * - Put the whole determination of the IG source term in a new subroutine 'determine_infragravity_source_sink_term' for clarity - Also makes it easier to experiment puting the source term determination in the ' do iter=1,niter '-loop * - Some style cleanup * - New implementation where IG source/sink term is first determined prior to iteration loop, to fill 'srcsh', based on effectively incident energy from previous SnapWave timestep - Then depending on new option, 'snapwave_iterative_srcsh=1, is new default', the subroutine is called again to determine srcsh now based on just updated incident wave energy. - If nr_sweeps = 4, this is slower because subroutine 'determine_infragravity_source_sink_term' is called multiple times - But, Hm0 IG converges now immediately, instead of after a few timesteps! - However, to speed up you can with snapwave_iterative_srcsh = 0 still keep old behaviour of only using the initial filling of 'srcsh'. Depends on application whether slower convergence IG wave height is problem - TODO: check whether slower performance can still be counteracted some way * - First significant speedup by putting ' do itheta = 1, ntheta' loop inside the subroutine * save intermediate attempts: - ! Update H(k), used in subroutine to calculate relative waterdepth 'gam' > TL: this creates a variability over sweeps/iterations/timesteps since 'gam' changes every time > seems more stable not to do this - OMP loop in solver did not seem to provide speedup - OMP loop in subroutine 'determine_infragravity_source_sink_term' made things slower - Even without these 3 attempts, performance loss of more times calling subroutine 'determine_infragravity_source_sink_term' in iteration loop seems negligible for small model * - In iteration loop, save a little bit of computation time by updating srcsh only for the first sweep (but still every iteration) * - Rename all references to 'srcsh' to more descriptive 'srcig' * - Do actually limit the depth to 5m, not only warning * Added text * Combine loops * - depth needed to be inout to be able to change it * Revert "Combine loops" This reverts commit 04a8e61. * - Determinen above, so that we can correctly allocatie nmindbnd * - Update H(k) in first sweep of first iteration, so that in first timestep we get already a good estimate of srcig and thereby Hig, and not just only conservative shoaling from the boundary - Only done once to prevent changing H(k) > and thereby gam and alphaig > every time, to prevent oscillations * - Fix quadtree hmax output * - set/check random seed for random_number - Todo: check whether this is best implementation and results in actual random seed if wmrandom = 1 * - SnapWave doesn't output directional spreading, and in sfincs_snapwave the mean wave direction is actually added, which gives wrong values in hisfile. Turn off for now * - After Maarten's suggestion take out determinatioon of depthprev and beta_local from every time determine_infragravity_source_sink_term is called - Did not take out beta_local completely, because it depends on the grid direction, that is updated in boundary_points for every timestep. - So now depthprev and beta_local only determined once, and not every iteration - Note - could be done still using w360 etc, but more work
1 parent f0f3adc commit 701b541

File tree

10 files changed

+302
-205
lines changed

10 files changed

+302
-205
lines changed

source/src/sfincs_data.f90

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -221,6 +221,7 @@ module sfincs_data
221221
logical :: friction2d
222222
logical :: advection_limiter
223223
logical :: advection_mask
224+
logical :: wmrandom
224225
!!!
225226
!!! sfincs_input.f90 switches
226227
integer storevelmax
@@ -532,7 +533,7 @@ module sfincs_data
532533
real*4, dimension(:), allocatable :: cg
533534
real*4, dimension(:), allocatable :: qb
534535
real*4, dimension(:), allocatable :: betamean
535-
real*4, dimension(:), allocatable :: srcsh
536+
real*4, dimension(:), allocatable :: srcig
536537
real*4, dimension(:), allocatable :: alphaig
537538

538539
! real*4, dimension(:), allocatable :: tauwavv

source/src/sfincs_domain.f90

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2291,8 +2291,8 @@ subroutine initialize_hydro()
22912291
qb = 0.0
22922292
allocate(betamean(np))
22932293
betamean = 0.0
2294-
allocate(srcsh(np))
2295-
srcsh = 0.0
2294+
allocate(srcig(np))
2295+
srcig = 0.0
22962296
allocate(alphaig(np))
22972297
alphaig = 0.0
22982298
endif

source/src/sfincs_input.f90

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -114,6 +114,7 @@ subroutine read_sfincs_input()
114114
call read_real_input(500,'wmtfilter',wmtfilter,600.0)
115115
call read_real_input(500,'wmfred',wavemaker_freduv,0.99)
116116
call read_char_input(500,'wmsignal',wmsigstr,'spectrum')
117+
call read_logical_input(500,'wmrandom',wmrandom,.true.)
117118
call read_char_input(500,'advection_scheme',advstr,'upw1')
118119
call read_real_input(500,'btrelax',btrelax,3600.0)
119120
!

source/src/sfincs_ncoutput.F90

Lines changed: 30 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ module sfincs_ncoutput
5050
integer :: patm_varid, wind_speed_varid, wind_dir_varid
5151
integer :: inp_varid, total_runtime_varid, average_dt_varid, status_varid
5252
integer :: hm0_varid, hm0ig_varid, zsm_varid, tp_varid, tpig_varid, wavdir_varid, dirspr_varid
53-
integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, qb_varid, beta_varid, srcsh_varid, alphaig_varid
53+
integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, qb_varid, beta_varid, srcig_varid, alphaig_varid
5454
!
5555
end type
5656
!
@@ -1500,12 +1500,12 @@ subroutine ncoutput_his_init()
15001500
NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, 'long_name', 'directionally averaged normalised bed slope'))
15011501
NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, 'coordinates', 'station_id station_name point_x point_y'))
15021502
!
1503-
NF90(nf90_def_var(his_file%ncid, 'srcsh', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%srcsh_varid)) ! time-varying water level point
1504-
NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, '_FillValue', FILL_VALUE))
1505-
NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'units', '-'))
1506-
NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'standard_name', 'directionally_averaged_ig_energy_source'))
1507-
NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'long_name', 'directionally averaged ig energy source'))
1508-
NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'coordinates', 'station_id station_name point_x point_y'))
1503+
NF90(nf90_def_var(his_file%ncid, 'srcig', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%srcig_varid)) ! time-varying water level point
1504+
NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, '_FillValue', FILL_VALUE))
1505+
NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'units', '-'))
1506+
NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'standard_name', 'directionally_averaged_ig_energy_source'))
1507+
NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'long_name', 'directionally averaged ig energy source'))
1508+
NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'coordinates', 'station_id station_name point_x point_y'))
15091509
!
15101510
NF90(nf90_def_var(his_file%ncid, 'alphaig', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%alphaig_varid)) ! time-varying water level point
15111511
NF90(nf90_put_att(his_file%ncid, his_file%alphaig_varid, '_FillValue', FILL_VALUE))
@@ -2224,7 +2224,7 @@ subroutine ncoutput_update_his(t,nthisout)
22242224
real*4, dimension(nobs) :: cgobs
22252225
real*4, dimension(nobs) :: qbobs
22262226
real*4, dimension(nobs) :: betaobs
2227-
real*4, dimension(nobs) :: srcshobs
2227+
real*4, dimension(nobs) :: srcigobs
22282228
real*4, dimension(nobs) :: alphaigobs
22292229
real*4, dimension(:), allocatable :: qq
22302230
!
@@ -2248,7 +2248,7 @@ subroutine ncoutput_update_his(t,nthisout)
22482248
cgobs = FILL_VALUE
22492249
qbobs = FILL_VALUE
22502250
betaobs = FILL_VALUE
2251-
srcshobs = FILL_VALUE
2251+
srcigobs = FILL_VALUE
22522252
alphaigobs = FILL_VALUE
22532253
!
22542254
do iobs = 1, nobs ! determine zs and prcp of obervation points at required timestep
@@ -2350,7 +2350,7 @@ subroutine ncoutput_update_his(t,nthisout)
23502350
cgobs(iobs) = cg(nm)
23512351
qbobs(iobs) = qb(nm)
23522352
betaobs(iobs) = betamean(nm)
2353-
srcshobs(iobs) = srcsh(nm)
2353+
srcigobs(iobs) = srcig(nm)
23542354
alphaigobs(iobs) = alphaig(nm)
23552355
!
23562356
endif
@@ -2392,7 +2392,7 @@ subroutine ncoutput_update_his(t,nthisout)
23922392
if (store_wave_direction) then
23932393
!
23942394
NF90(nf90_put_var(his_file%ncid, his_file%wavdir_varid, wavdirobs, (/1, nthisout/)))
2395-
NF90(nf90_put_var(his_file%ncid, his_file%dirspr_varid, dirsprobs, (/1, nthisout/)))
2395+
!NF90(nf90_put_var(his_file%ncid, his_file%dirspr_varid, dirsprobs, (/1, nthisout/)))
23962396
!
23972397
endif
23982398
!
@@ -2414,7 +2414,7 @@ subroutine ncoutput_update_his(t,nthisout)
24142414
!
24152415
NF90(nf90_put_var(his_file%ncid, his_file%qb_varid, qbobs, (/1, nthisout/)))
24162416
NF90(nf90_put_var(his_file%ncid, his_file%beta_varid, betaobs, (/1, nthisout/)))
2417-
NF90(nf90_put_var(his_file%ncid, his_file%srcsh_varid, srcshobs, (/1, nthisout/)))
2417+
NF90(nf90_put_var(his_file%ncid, his_file%srcig_varid, srcigobs, (/1, nthisout/)))
24182418
NF90(nf90_put_var(his_file%ncid, his_file%alphaig_varid, alphaigobs, (/1, nthisout/)))
24192419
!
24202420
endif
@@ -2687,27 +2687,26 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout)
26872687
!
26882688
! Write maximum water depth
26892689
if (subgrid .eqv. .false. .or. store_hsubgrid .eqv. .true.) then
2690+
!
26902691
zstmp = FILL_VALUE
26912692
!
2692-
if (subgrid) then
2693-
do nm = 1, np
2694-
!
2695-
if ( (zsmax(nm) - subgrid_z_zmin(nm)) > huthresh) then
2696-
zstmp(nm) = zsmax(nm) - subgrid_z_zmin(nm)
2697-
endif
2698-
!
2699-
enddo
2700-
else
2701-
do nm = 1, np
2702-
!
2703-
if ( (zsmax(nm) - zb(nm)) > huthresh) then
2704-
zstmp(nm) = zsmax(nm) - zb(nm)
2705-
endif
2706-
enddo
2707-
endif
2708-
endif
2709-
!
2710-
if (subgrid .eqv. .false. .or. store_hsubgrid .eqv. .true.) then
2693+
do nmq = 1, quadtree_nr_points
2694+
!
2695+
nm = index_sfincs_in_quadtree(nmq)
2696+
!
2697+
if (kcs(nm)>0) then
2698+
if (subgrid) then
2699+
if ( (zsmax(nm) - subgrid_z_zmin(nm)) > huthresh) then
2700+
zstmp(nmq) = zsmax(nm) - subgrid_z_zmin(nm)
2701+
endif
2702+
else
2703+
if ( (zsmax(nm) - zb(nm)) > huthresh) then
2704+
zstmp(nmq) = zsmax(nm) - zb(nm)
2705+
endif
2706+
endif
2707+
endif
2708+
enddo
2709+
!
27112710
NF90(nf90_put_var(map_file%ncid, map_file%hmax_varid, zstmp, (/1, ntmaxout/))) ! write hmax
27122711
endif
27132712
!

source/src/sfincs_snapwave.f90

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ module sfincs_snapwave
2222
real*4, dimension(:), allocatable :: snapwave_cg
2323
real*4, dimension(:), allocatable :: snapwave_Qb
2424
real*4, dimension(:), allocatable :: snapwave_beta
25-
real*4, dimension(:), allocatable :: snapwave_srcsh
25+
real*4, dimension(:), allocatable :: snapwave_srcig
2626
real*4, dimension(:), allocatable :: snapwave_alphaig
2727
integer, dimension(:,:), allocatable :: snapwave_connected_nodes
2828
integer*4, dimension(:), allocatable :: index_snapwave_in_sfincs
@@ -140,7 +140,7 @@ subroutine update_wave_field(t, tloop)
140140
real*4, dimension(:), allocatable :: cg0
141141
real*4, dimension(:), allocatable :: qb0
142142
real*4, dimension(:), allocatable :: beta0
143-
real*4, dimension(:), allocatable :: srcsh0
143+
real*4, dimension(:), allocatable :: srcig0
144144
real*4, dimension(:), allocatable :: alphaig0
145145
integer :: ip, ii, m, n, nm, nmu, idir
146146
real*4 :: f
@@ -157,7 +157,7 @@ subroutine update_wave_field(t, tloop)
157157
allocate(cg0(np))
158158
allocate(qb0(np))
159159
allocate(beta0(np))
160-
allocate(srcsh0(np))
160+
allocate(srcig0(np))
161161
allocate(alphaig0(np))
162162
!
163163
fwx0 = 0.0
@@ -169,7 +169,7 @@ subroutine update_wave_field(t, tloop)
169169
cg0 = 0.0
170170
qb0 = 0.0
171171
beta0 = 0.0
172-
srcsh0 = 0.0
172+
srcig0 = 0.0
173173
alphaig0 = 0.0
174174
!
175175
! Determine SnapWave water depth
@@ -221,7 +221,7 @@ subroutine update_wave_field(t, tloop)
221221
cg0(nm) = snapwave_cg(ip)
222222
qb0(nm) = snapwave_Qb(ip)
223223
beta0(nm) = snapwave_beta(ip)
224-
srcsh0(nm) = snapwave_srcsh(ip)
224+
srcig0(nm) = snapwave_srcig(ip)
225225
alphaig0(nm) = snapwave_alphaig(ip)
226226
if (store_wave_direction) then
227227
mean_wave_direction(nm) = 270.0 - snapwave_mean_direction(ip)*180/pi
@@ -243,7 +243,7 @@ subroutine update_wave_field(t, tloop)
243243
cg0(nm) = 0.0
244244
qb0(nm) = 0.0
245245
beta0(nm) = 0.0
246-
srcsh0(nm) = 0.0
246+
srcig0(nm) = 0.0
247247
alphaig0(nm) = 0.0
248248
if (store_wave_direction) then
249249
mean_wave_direction(nm) = 0.0
@@ -263,7 +263,7 @@ subroutine update_wave_field(t, tloop)
263263
cg(nm) = cg0(nm)
264264
qb(nm) = qb0(nm)
265265
betamean(nm) = beta0(nm)
266-
srcsh(nm) = srcsh0(nm)
266+
srcig(nm) = srcig0(nm)
267267
alphaig(nm) = alphaig0(nm)
268268
!
269269
endif
@@ -328,7 +328,7 @@ subroutine compute_snapwave(t)
328328
snapwave_H = H
329329
snapwave_H_ig = H_ig
330330
snapwave_mean_direction = thetam
331-
snapwave_directional_spreading = thetam ! L: CORRECT? > is not spreading but mean direction?
331+
snapwave_directional_spreading = thetam ! TL: CORRECT? > is not spreading but mean direction?
332332
snapwave_Fx = Fx
333333
snapwave_Fy = Fy
334334
snapwave_Dw = Dw
@@ -338,7 +338,7 @@ subroutine compute_snapwave(t)
338338
snapwave_cg = cg
339339
snapwave_Qb = Qb
340340
snapwave_beta = beta
341-
snapwave_srcsh = srcsh
341+
snapwave_srcig = srcig
342342
snapwave_alphaig = alphaig
343343
!
344344
! Wave periods from SnapWave, used in e.g. wavemakers - TL: moved behind call update_boundary_conditions & compute_wave_field so values at first timestep are not 0
@@ -399,6 +399,7 @@ subroutine read_snapwave_input()
399399
call read_real_input(500,'snapwave_alphaigfac',alphaigfac,1.0) ! Multiplication factor for IG shoaling source/sink term
400400
call read_real_input(500,'snapwave_baldock_ratio_ig',baldock_ratio_ig,0.2)
401401
call read_int_input(500,'snapwave_ig_opt',ig_opt,1)
402+
call read_int_input(500,'snapwave_iterative_srcig',iterative_srcig,1) ! Option whether to calculate IG source/sink term in iterative lower (better, but potentially slower, 1=default), or effectively based on previous timestep (faster, potential mismatch, =0)
402403
!
403404
! IG boundary conditions options:
404405
call read_int_input(500,'snapwave_use_herbers',herbers_opt,1) ! Choice whether you want IG Hm0&Tp be calculated by herbers (=1, default), or want to specify user defined values (0> then snapwave_eeinc2ig & snapwave_Tinc2ig are used)
@@ -416,7 +417,7 @@ subroutine read_snapwave_input()
416417
call read_char_input(500,'snapwave_btpfile',btpfile,'')
417418
call read_char_input(500,'snapwave_bwdfile',bwdfile,'')
418419
call read_char_input(500,'snapwave_bdsfile',bdsfile,'')
419-
call read_char_input(500,'snapwave_upwfile',upwfile,'')
420+
call read_char_input(500,'snapwave_upwfile',upwfile,'snapwave.upw')
420421
call read_char_input(500,'snapwave_mskfile',mskfile,'')
421422
call read_char_input(500,'snapwave_depfile',depfile,'none')
422423
call read_char_input(500,'snapwave_ncfile', gridfile,'snapwave_net.nc')

source/src/sfincs_wavemaker.f90

Lines changed: 14 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,8 @@ subroutine read_wavemaker_polylines()
4848
real*4 :: phip
4949
!
5050
real*4 :: r
51+
integer, allocatable :: seed(:)
52+
integer :: nseed
5153
! real*4 :: dxcross
5254
! real*4 :: dst
5355
! real*4 :: xxx
@@ -1287,8 +1289,19 @@ subroutine read_wavemaker_polylines()
12871289
!
12881290
endif
12891291
!
1290-
! Infragravity frequencies
1292+
! Random seed wavemaker
1293+
!
1294+
call random_seed(size = nseed)
1295+
if (.NOT. wmrandom) then
1296+
call random_seed(put = [ 2147483562, 2147483398])
1297+
endif
12911298
!
1299+
allocate(seed(nseed))
1300+
call random_seed(get=seed)
1301+
write (*, *)'Wavemaker random seed: ',seed
1302+
!
1303+
! Infragravity frequencies
1304+
!
12921305
allocate(freqig(nfreqsig))
12931306
allocate(costig(nfreqsig))
12941307
allocate(phiig(nfreqsig))

source/src/snapwave/snapwave_data.f90

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -11,12 +11,9 @@ module snapwave_data
1111
real*4, dimension(:), allocatable :: kwav_ig, nwav_ig ! wave number, ratio Cg/C
1212
real*4, dimension(:), allocatable :: C, Cg ! wave celerity, group velocity0
1313
real*4, dimension(:), allocatable :: C_ig, Cg_ig ! wave celerity, group velocity0
14-
real*4, dimension(:,:), allocatable :: Sxx ! directional radiation stress Sxx
1514
real*4, dimension(:), allocatable :: fw ! friction coefficient
1615
real*4, dimension(:), allocatable :: fw_ig ! friction coefficient
1716
real*4, dimension(:), allocatable :: H, H_ig ! rms wave height
18-
real*4, dimension(:), allocatable :: H_ig_old ! IG wave height at previous timestep
19-
real*4, dimension(:), allocatable :: H_inc_old ! Incident wave height at previous timestep
2017
real*4, dimension(:), allocatable :: Dw,Df ! dissipation due to breaking, bed friction
2118
real*4, dimension(:), allocatable :: Dw_ig,Df_ig ! dissipation due to breaking, bed friction for IG
2219
real*4, dimension(:), allocatable :: F ! wave force Dw/C/rho/depth
@@ -59,7 +56,7 @@ module snapwave_data
5956
!
6057
real*4, dimension(:), allocatable :: Qb
6158
real*4, dimension(:), allocatable :: beta
62-
real*4, dimension(:), allocatable :: srcsh
59+
real*4, dimension(:), allocatable :: srcig
6360
real*4, dimension(:), allocatable :: alphaig
6461
!
6562
integer*4, dimension(:), allocatable :: index_snapwave_in_quadtree
@@ -162,6 +159,8 @@ module snapwave_data
162159
real*4 :: baldock_ratio_ig ! option controlling from what depth wave breaking should take place for IG waves: (Hk>baldock_ratio*Hmx(k)), default baldock_ratio_ig=0.2
163160
integer :: igwaves_opt ! option of IG waves on (1) or off (0)
164161
integer :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking)
162+
integer :: iterative_srcig ! option whether IG source/sink term should be calculated in the iterative loop again (iterative_srcig = 1, is a bit slower),
163+
! ... or just a priori based on effectively incident wave energy from previous timestep only
165164
real*4 :: snapwave_alpha_ig,gamma_ig ! coefficients in Baldock wave breaking dissipation model for IG waves
166165
real*4 :: shinc2ig ! Ratio of how much of the calculated IG wave source term, is subtracted from the incident wave energy (0-1, 0=default)
167166
real*4 :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0

source/src/snapwave/snapwave_domain.f90

Lines changed: 20 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -114,13 +114,10 @@ subroutine initialize_snapwave_domain()
114114
allocate(nwav_ig(no_nodes))
115115
allocate(C_ig(no_nodes))
116116
allocate(Cg_ig(no_nodes))
117-
allocate(Sxx(ntheta,no_nodes))
118117
allocate(sinhkh_ig(no_nodes))
119118
allocate(Hmx_ig(no_nodes))
120119
allocate(fw_ig(no_nodes))
121-
allocate(H_ig(no_nodes))
122-
allocate(H_ig_old(no_nodes))
123-
allocate(H_inc_old(no_nodes))
120+
allocate(H_ig(no_nodes))
124121
allocate(Dw(no_nodes))
125122
allocate(Dw_ig(no_nodes))
126123
allocate(F(no_nodes))
@@ -131,7 +128,7 @@ subroutine initialize_snapwave_domain()
131128
allocate(thetam(no_nodes))
132129
allocate(Qb(no_nodes))
133130
allocate(beta(no_nodes))
134-
allocate(srcsh(no_nodes))
131+
allocate(srcig(no_nodes))
135132
allocate(alphaig(no_nodes))
136133
! allocate(uorb(no_nodes))
137134
allocate(ctheta(ntheta,no_nodes))
@@ -196,7 +193,8 @@ subroutine initialize_snapwave_domain()
196193
ds360d0 = 0.d0
197194
w360d0 = 0.d0
198195
prev360 = 0
199-
H_ig_old = 0.0
196+
H = 0.0
197+
H_ig = 0.0
200198
!
201199
generate_upw = .false.
202200
exists = .true.
@@ -274,6 +272,8 @@ subroutine initialize_snapwave_domain()
274272
!
275273
endif
276274
!
275+
nb = 0
276+
!
277277
do k=1,no_nodes
278278
!if (msk(k)==3) msk(k) = 1 ! Set outflow points to regular points > now should become neumann, so don't do this
279279
do itheta=1,ntheta360
@@ -282,25 +282,27 @@ subroutine initialize_snapwave_domain()
282282
if (inout>0) msk(k) = 2
283283
endif
284284
enddo
285-
enddo
286-
!
287-
nb = 0
288-
do k = 1, no_nodes
289-
if (msk(k)==1) then
290-
inner(k) = .true.
291-
else
292-
inner(k) = .false.
285+
!
286+
if (msk(k)>1) then
287+
nb = nb + 1
293288
endif
294-
if (msk(k)==2) nb = nb + 1
289+
!
295290
enddo
296291
!
297292
allocate(nmindbnd(nb))
298293
!
299-
nb = 0
300294
do k = 1, no_nodes
295+
!
296+
if (msk(k)==1) then
297+
inner(k) = .true.
298+
else
299+
inner(k) = .false.
300+
endif
301+
!
301302
if (msk(k)>1) then
302-
nb = nb + 1
303-
nmindbnd(nb) = k
303+
!
304+
nmindbnd(nb) = k
305+
!
304306
endif
305307
enddo
306308
!

0 commit comments

Comments
 (0)