diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj
index 4ba34059..c5982749 100644
--- a/source/sfincs_lib/sfincs_lib.vfproj
+++ b/source/sfincs_lib/sfincs_lib.vfproj
@@ -44,7 +44,11 @@
-
+
+
+
+
+
diff --git a/source/src/Makefile.am b/source/src/Makefile.am
index 2ecb62ea..055454cb 100644
--- a/source/src/Makefile.am
+++ b/source/src/Makefile.am
@@ -37,6 +37,7 @@ libsfincs_la_SOURCES = \
sfincs_boundaries.f90 \
sfincs_continuity.f90 \
sfincs_crosssections.f90 \
+ sfincs_runup_gauges.f90 \
sfincs_discharges.f90 \
sfincs_subgrid.F90 \
sfincs_domain.f90 \
@@ -47,7 +48,7 @@ libsfincs_la_SOURCES = \
sfincs_snapwave.f90 \
deg2utm.f90 \
sfincs_meteo.f90 \
- bicgstab_solver_ilu.f90 \
+ bicgstab_solver_ilu.f90 \
sfincs_nonhydrostatic.f90 \
sfincs_ncoutput.F90 \
sfincs_output.f90 \
diff --git a/source/src/sfincs_boundaries.f90 b/source/src/sfincs_boundaries.f90
index e3684912..b77d2a4d 100644
--- a/source/src/sfincs_boundaries.f90
+++ b/source/src/sfincs_boundaries.f90
@@ -1,6 +1,7 @@
module sfincs_boundaries
use sfincs_log
+ use sfincs_error
contains
@@ -14,6 +15,7 @@ subroutine read_boundary_data()
implicit none
!
integer n, itb, ib, stat, ifreq, iok
+ logical :: ok
!
real*4 dummy,r
!
@@ -38,6 +40,8 @@ subroutine read_boundary_data()
!
call write_log('Info : reading water level boundaries', 0)
!
+ ok = check_file_exists(bndfile, 'Boundary points file', .true.)
+ !
open(500, file=trim(bndfile))
do while(.true.)
read(500,*,iostat = stat)dummy
@@ -54,6 +58,8 @@ subroutine read_boundary_data()
!
! Read water level boundary conditions file
!
+ ok = check_file_exists(bzsfile, 'Boundary time series file', .true.)
+ !
open(500, file=trim(bzsfile))
do while(.true.)
read(500,*,iostat = stat)dummy
@@ -79,6 +85,8 @@ subroutine read_boundary_data()
allocate(zsi_bnd(nbnd,ntbnd))
allocate(zsit_bnd(nbnd))
!
+ ok = check_file_exists(bzifile, 'Boundary infragravity time series file', .true.)
+ !
open(500, file=trim(bzifile))
do itb = 1, ntbnd
read(500,*)dummy,(zsi_bnd(ib, itb), ib = 1, nbnd)
@@ -110,7 +118,7 @@ subroutine read_boundary_data()
!
endif
!
- elseif (include_boundaries) then
+ elseif (water_level_boundaries_in_mask) then
!
write(logstr,'(a)')'Warning! Boundary cells found in mask, without boundary conditions. Using water level of 0.0 m at these points.'
call write_log(logstr, 1)
@@ -161,114 +169,132 @@ subroutine read_boundary_data()
call write_log(logstr, 1)
!
endif
+! !
+! ! Read wave boundaries
+! !
+! nwbnd = 0
+! ntwbnd = 0
+! !
+! if (bwvfile(1:4) /= 'none') then
+! !
+! write(logstr,'(a)')'Info : reading wave boundaries'
+! call write_log(logstr, 0)
+! !
+! ! Locations
+! !
+! open(500, file=trim(bwvfile))
+! do while(.true.)
+! read(500,*,iostat = stat)dummy
+! if (stat<0) exit
+! nwbnd = nwbnd + 1
+! enddo
+! rewind(500)
+! allocate(x_bwv(nwbnd))
+! allocate(y_bwv(nwbnd))
+! do n = 1, nwbnd
+! read(500,*)x_bwv(n),y_bwv(n)
+! enddo
+! close(500)
+! !
+! ! Wave time series
+! !
+! ! First find times in bhs file
+! !
+! open(500, file=trim(bhsfile))
+! do while(.true.)
+! read(500,*,iostat = stat)dummy
+! if (stat<0) exit
+! ntwbnd = ntwbnd + 1
+! enddo
+! close(500)
+! !
+! allocate(t_bwv(ntwbnd))
+! allocate(hst_bwv(nwbnd))
+! allocate(l0t_bwv(nwbnd))
+! allocate(tpt_bwv(nwbnd))
+! allocate(wdt_bwv(nwbnd))
+!! allocate(setupt_bwv(nwbnd))
+! !
+! ! Hs (significant wave height)
+! ! Times in btp and bwd files must be the same as in bhs file!
+! !
+! open(500, file=trim(bhsfile))
+! allocate(hs_bwv(nwbnd,ntwbnd))
+! do itb = 1, ntwbnd
+! read(500,*)t_bwv(itb),(hs_bwv(ib, itb), ib = 1, nwbnd)
+! enddo
+! close(500)
+! !
+! ! Tp (peak period)
+! !
+! open(500, file=trim(btpfile))
+! allocate(tp_bwv(nwbnd,ntwbnd))
+! do itb = 1, ntwbnd
+! read(500,*)t_bwv(itb),(tp_bwv(ib, itb), ib = 1, nwbnd)
+! enddo
+! close(500)
+! !
+! if (bwdfile(1:4) /= 'none') then
+! !
+! ! Wd (wave direction)
+! !
+! open(500, file=trim(bwdfile))
+! allocate(wd_bwv(nwbnd,ntwbnd))
+! do itb = 1, ntwbnd
+! read(500,*)t_bwv(itb),(wd_bwv(ib, itb), ib = 1, nwbnd)
+! enddo
+! close(500)
+! wd_bwv = (270.0 - wd_bwv)*pi/180 ! Convert to cartesian going to
+! !
+! endif
!
- ! Read wave boundaries
+ ! Now the downstream river boundaries. No time series, just points, slope and direction
!
- nwbnd = 0
- ntwbnd = 0
+ nbdr = 0
!
- if (bwvfile(1:4) /= 'none') then
- !
- write(logstr,'(a)')'Info : reading wave boundaries'
- call write_log(logstr, 0)
- !
- ! Locations
- !
- open(500, file=trim(bwvfile))
- do while(.true.)
- read(500,*,iostat = stat)dummy
- if (stat<0) exit
- nwbnd = nwbnd + 1
- enddo
- rewind(500)
- allocate(x_bwv(nwbnd))
- allocate(y_bwv(nwbnd))
- do n = 1, nwbnd
- read(500,*)x_bwv(n),y_bwv(n)
- enddo
- close(500)
- !
- ! Wave time series
- !
- ! First find times in bhs file
- !
- open(500, file=trim(bhsfile))
- do while(.true.)
- read(500,*,iostat = stat)dummy
- if (stat<0) exit
- ntwbnd = ntwbnd + 1
- enddo
- close(500)
- !
- allocate(t_bwv(ntwbnd))
- allocate(hst_bwv(nwbnd))
- allocate(l0t_bwv(nwbnd))
- allocate(tpt_bwv(nwbnd))
- allocate(wdt_bwv(nwbnd))
-! allocate(setupt_bwv(nwbnd))
- !
- ! Hs (significant wave height)
- ! Times in btp and bwd files must be the same as in bhs file!
- !
- open(500, file=trim(bhsfile))
- allocate(hs_bwv(nwbnd,ntwbnd))
- do itb = 1, ntwbnd
- read(500,*)t_bwv(itb),(hs_bwv(ib, itb), ib = 1, nwbnd)
- enddo
- close(500)
- !
- ! Tp (peak period)
- !
- open(500, file=trim(btpfile))
- allocate(tp_bwv(nwbnd,ntwbnd))
- do itb = 1, ntwbnd
- read(500,*)t_bwv(itb),(tp_bwv(ib, itb), ib = 1, nwbnd)
- enddo
- close(500)
- !
- if (bwdfile(1:4) /= 'none') then
+ if (downstream_river_boundaries_in_mask) then
+ !
+ if (bdrfile(1:4) /= 'none') then
!
- ! Wd (wave direction)
+ write(logstr,'(a)')'Info : reading downstream river boundaries'
+ call write_log(logstr, 0)
!
- open(500, file=trim(bwdfile))
- allocate(wd_bwv(nwbnd,ntwbnd))
- do itb = 1, ntwbnd
- read(500,*)t_bwv(itb),(wd_bwv(ib, itb), ib = 1, nwbnd)
+ ok = check_file_exists(bdrfile, 'Downstream boundary points file', .true.)
+ !
+ open(500, file=trim(bdrfile))
+ do while(.true.)
+ read(500,*,iostat = stat)dummy
+ if (stat<0) exit
+ nbdr = nbdr + 1
+ enddo
+ rewind(500)
+ allocate(x_bdr(nbdr))
+ allocate(y_bdr(nbdr))
+ allocate(slope_bdr(nbdr))
+ allocate(azimuth_bdr(nbdr))
+ do n = 1, nbdr
+ read(500,*)x_bdr(n), y_bdr(n), slope_bdr(n), azimuth_bdr(n)
enddo
close(500)
- wd_bwv = (270.0 - wd_bwv)*pi/180 ! Convert to cartesian going to
+ !
+ else
+ !
+ ! This really should not happen
+ !
+ call stop_sfincs('Downstream river points found in mask, without boundary conditions (missing bdrfile in sfincs.inp) !', 1)
!
endif
!
-! if (stufile(1:4) /= 'none') then
-! !
-! ! Set-up
-! !
-! open(500, file=trim(stufile))
-! allocate(setup_bwv(nwbnd, ntwbnd))
-! do itb = 1, ntwbnd
-! read(500,*)t_bwv(itb),(setup_bwv(ib, itb), ib = 1, nwbnd)
-! enddo
-! close(500)
-! !
-! endif
+ else
!
- endif
-! !
-! ! Infragravity frequencies
-! !
-! allocate(freqig(nfreqsig))
-! allocate(costig(nfreqsig))
-! allocate(phiig(nfreqsig))
-! allocate(dphiig(nfreqsig))
-! dfreqig = (freqmaxig - freqminig)/nfreqsig
-! do ifreq = 1, nfreqsig
-! freqig(ifreq) = freqminig + ifreq*dfreqig - 0.5*dfreqig
-! call RANDOM_NUMBER(r)
-! phiig(ifreq) = r*2*3.1416
-!! call RANDOM_NUMBER(r)
-! dphiig(ifreq) = 1.0e-6*2*3.1416/freqig(ifreq)
-! enddo
+ if (bdrfile(1:4) /= 'none') then
+ !
+ write(logstr,'(a)')'Warning : Found bdr file in sfincs.inp, but no downstream river points in were found in the mask!'
+ call write_log(logstr, 0)
+ !
+ endif
+ !
+ endif
!
end subroutine
@@ -279,161 +305,286 @@ subroutine find_boundary_indices()
!
implicit none
!
- ! For each grid boundary point (kcs=2) :
+ ! For each grid boundary point (kcs=2, ) :
!
! Determine indices and weights of boundary points
! For tide and surge, these are the indices and weights of the points in the bnd file
! For waves, these are the indices and weights of the points in the cst file
!
- integer nm, m, n, nb, ib1, ib2, ib, ic
+ integer nm, m, n, nb, ib1, ib2, ib, ic, ibnd, ibdr, iref
+ integer nm_1, nm_2, nm_3, nmi
!
- real x, y, dst1, dst2, dst
+ real x, y, dst1, dst2, dst, phi_riv, phi_r
+ real*4 :: w_1, w_2, w_3, d_1, d_2, d_3
!
- if (nbnd>0 .and. ngbnd>0) then
+ ! Check there are points from bnd file and/or bdr file and that kcs mask contains 2/3/5/6
+ !
+ if (ngbnd == 0) then
!
- ! Count number of boundary points
+ if (nbnd > 0 .or. nbdr > 0) then
+ !
+ write(logstr,'(a)')'Warning : no open boundary points found in mask!'
+ call write_log(logstr, 1)
+ !
+ endif
!
- ! Allocate boundary arrays
+ return
!
- ! First water levels
+ endif
+ !
+ if (nbnd == 0 .and. nbdr == 0) then
!
- ! Water level arrays
+ !write(logstr,'(a)')'Warning : no open boundary points found in mask!'
+ !call write_log(logstr, 1)
!
- allocate(ind1_bnd_gbp(ngbnd))
- allocate(ind2_bnd_gbp(ngbnd))
- allocate(fac_bnd_gbp(ngbnd))
+ return
!
- ! Water levels at boundary
+ endif
+ !
+ ! Allocate boundary arrays
+ !
+ allocate(ind1_bnd_gbp(ngbnd))
+ allocate(ind2_bnd_gbp(ngbnd))
+ allocate(fac_bnd_gbp(ngbnd))
+ !
+ ind1_bnd_gbp = 0
+ ind2_bnd_gbp = 0
+ fac_bnd_gbp = 0.0
+ !
+ if (downstream_river_boundaries_in_mask) then
!
- if (waves) then
- !
- ! Wave arrays
- !
- allocate(ind1_cst_gbp(ngbnd))
- allocate(ind2_cst_gbp(ngbnd))
- allocate(fac_cst_gbp(ngbnd))
- !
- endif
+ ! There are downstream river boundaries, so we need indices, weights and distances of upstream points
!
- ! Find two closest boundary condition points for each boundary point
- ! And the two closest coastline points
+ allocate(slope_gbp(ngbnd))
+ allocate(nm_nbr_gbp(3, ngbnd))
+ allocate(w_nbr_gbp(3, ngbnd))
+ allocate(d_nbr_gbp(3, ngbnd))
!
- nb = 0
+ nm_nbr_gbp = 0
+ w_nbr_gbp = 0.0
+ d_nbr_gbp = 0.0
!
- ! Loop through all grid points
+ endif
+ !
+ if (neumann_boundaries_in_mask) then
!
- do nm = 1, np
- !
- ! Check if this point is a boundary point
+ ! There are Neumann boundaries, so we need nm indices of internal points
+ !
+ allocate(nmi_gbp(ngbnd))
+ nmi_gbp = 0
+ !
+ endif
+ !
+ ! Find two closest boundary condition points for each boundary point
+ !
+ nb = 0
+ !
+ ! Loop through all grid boundary points
+ !
+ do ib = 1, ngbnd
+ !
+ nm = nmindbnd(ib)
+ !
+ x = z_xz(nm)
+ y = z_yz(nm)
+ !
+ if (kcs(nm) == 2) then ! This cell is a water level boundary point
!
- if (kcs(nm) > 1) then
+ if (nbnd > 1) then
!
- nb = nb + 1
+ ! Multiple points in bnd file, so use distance-based weighting
!
- x = z_xz(nm)
- y = z_yz(nm)
+ dst1 = 1.0e10
+ dst2 = 1.0e10
+ ib1 = 0
+ ib2 = 0
!
- ! Indices and weights for water level boundaries
+ ! Loop through all water level boundary points in bnd file
!
- if (nbnd>1) then
+ do ibnd = 1, nbnd
!
- dst1 = 1.0e10
- dst2 = 1.0e10
- ib1 = 0
- ib2 = 0
+ ! Compute distance of this point to grid boundary point
!
- ! Loop through all water level boundary points
+ dst = sqrt((x_bnd(ibnd) - x)**2 + (y_bnd(ibnd) - y)**2)
!
- do ib = 1, nbnd
+ if (dst < dst1) then
!
- ! Compute distance of this point to grid boundary point
+ ! Nearest point found
!
- dst = sqrt((x_bnd(ib) - x)**2 + (y_bnd(ib) - y)**2)
+ dst2 = dst1
+ ib2 = ib1
+ dst1 = dst
+ ib1 = ibnd
!
- if (dst1) then
- !
- dst1 = 1.0e10
- dst2 = 1.0e10
- ib1 = 0
- ib2 = 0
- !
- ! Loop through all water level boundary points
- !
- do ic = 1, ncst
- !
- ! Compute distance of this point to grid boundary point
- !
- dst = sqrt((x_cst(ic) - x)**2 + (y_cst(ic) - y)**2)
- !
- if (dst 0) then
+ if (kcs(nmi) == 1) nmi_gbp(ib) = nmi
+ endif
+ !
+ nmi = find_sfincs_cell(n + 1, m, iref)
+ if (nmi > 0) then
+ if (kcs(nmi) == 1) nmi_gbp(ib) = nmi
+ endif
+ !
+ nmi = find_sfincs_cell(n, m - 1, iref)
+ if (nmi > 0) then
+ if (kcs(nmi) == 1) nmi_gbp(ib) = nmi
+ endif
+ !
+ nmi = find_sfincs_cell(n - 1, m, iref)
+ if (nmi > 0) then
+ if (kcs(nmi) == 1) nmi_gbp(ib) = nmi
+ endif
+ !
+ endif
+ enddo
!
end subroutine
@@ -452,104 +603,108 @@ subroutine update_boundary_points(t)
!
real*4 zstb, tbfac, hs, tp, wd, tb
!
- if (nbnd>0) then
+ if (nbnd == 0) return
+ !
+ ! Interpolate boundary conditions in time
+ !
+ if (t_bnd(1) > (t - 1.0e-3)) then ! use first time in boundary conditions
!
- ! Start with updating values at boundary polylines
+ itb0 = 1
+ itb1 = 1
+ tb = t_bnd(itb0)
!
- ! Water levels
+ elseif (t_bnd(ntbnd) < (t + 1.0e-3)) then ! use last time in boundary conditions
!
- ! Interpolate boundary conditions in time
+ itb0 = ntbnd
+ itb1 = ntbnd
+ tb = t_bnd(itb0)
!
- if (t_bnd(1) > (t - 1.0e-3)) then ! use first time in boundary conditions
- !
- itb0 = 1
- itb1 = 1
- tb = t_bnd(itb0)
- !
- elseif (t_bnd(ntbnd) < (t + 1.0e-3)) then ! use last time in boundary conditions
- !
- itb0 = ntbnd
- itb1 = ntbnd
- tb = t_bnd(itb0)
- !
- else
- !
- do itb = itbndlast, ntbnd ! Loop in time
- if (t_bnd(itb) > (t + 1.0e-6)) then
- itb0 = itb - 1
- itb1 = itb
- tb = t
- itbndlast = itb - 1
- exit
- endif
- enddo
- !
- endif
+ else
!
- tbfac = (tb - t_bnd(itb0))/max(t_bnd(itb1) - t_bnd(itb0), 1.0e-6)
+ do itb = itbndlast, ntbnd ! Loop in time
+ if (t_bnd(itb) > (t + 1.0e-6)) then
+ itb0 = itb - 1
+ itb1 = itb
+ tb = t
+ itbndlast = itb - 1
+ exit
+ endif
+ enddo
!
- do ib = 1, nbnd ! Loop along boundary points
- !
- ! Tide and surge
- !
- zstb = zs_bnd(ib, itb0) + (zs_bnd(ib, itb1) - zs_bnd(ib, itb0))*tbfac
+ endif
+ !
+ tbfac = (tb - t_bnd(itb0))/max(t_bnd(itb1) - t_bnd(itb0), 1.0e-6)
+ !
+ do ib = 1, nbnd ! Loop along boundary points
+ !
+ ! Tide and surge
+ !
+ zstb = zs_bnd(ib, itb0) + (zs_bnd(ib, itb1) - zs_bnd(ib, itb0))*tbfac
+ !
+ zst_bnd(ib) = zstb
+ !
+ if (bzifile(1:4) /= 'none') then
!
- zst_bnd(ib) = zstb
+ ! Incoming infragravity waves
!
- if (bzifile(1:4) /= 'none') then
- !
- ! Incoming infragravity waves
- !
- zsit_bnd(ib) = zsi_bnd(ib, itb0) + (zsi_bnd(ib, itb1) - zsi_bnd(ib, itb0))*tbfac
- !
- endif
+ zsit_bnd(ib) = zsi_bnd(ib, itb0) + (zsi_bnd(ib, itb1) - zsi_bnd(ib, itb0))*tbfac
!
- enddo
- endif
+ endif
+ !
+ enddo
!
end subroutine
-
-
- subroutine update_boundary_conditions(t,dt)
+ subroutine update_boundary_conditions(t, dt)
!
- ! Update values at boundary points
+ ! Update water level at boundary grid points
!
use sfincs_data
!
implicit none
!
- integer ib, ifreq, ic
+ integer ib, ifreq, ic, nm, ibuv, nmi, nmb, indb, iw
!
real*8 :: t
real :: dt
!
- real*4 zst, hst, l0t, wdt, ksi, zsetup, zig, a, angdiff, angfac, hm0ig, tmmin01ig, fp
+ real*8 zst
+ real*4 hst, l0t, wdt, ksi, zsetup, zig, a, angdiff, angfac, hm0ig, tmmin01ig, fp
real*4 zs0act
real*4 smfac
real*4 zs0smooth
+ real*4 sumw
!
! Set water level in all boundary points on grid
!
do ib = 1, ngbnd
!
- if (ibndtype(ib) == 1) then
+ nmb = nmindbnd(ib)
+ !
+ ! kcs = 1 : regular point
+ ! kcs = 2 : water level boundary point
+ ! kcs = 3 : outflow boundary point
+ ! kcs = 4 : wave maker point
+ ! kcs = 5 : river outflow point (dzs/dx = i)
+ ! kcs = 6 : lateral (coastal) boundary point (Neumann dzs/dx = 0.0)
+ !
+ if (kcs(nmb) == 2) then
!
- ! Regular boundary point (otherwise (for kcs==3) zsb and zsb0 have already been initialized at zbmin)
+ ! Regular water level boundary point
!
- ! Water levels (surge+tide)
+ ! Get water levels (surge + tide) from time series boundary conditions
!
- if (nbnd>1) then
+ if (nbnd > 1) then
!
! Interpolation of nearby points
!
- zst = zst_bnd(ind1_bnd_gbp(ib))*fac_bnd_gbp(ib) + zst_bnd(ind2_bnd_gbp(ib))*(1.0 - fac_bnd_gbp(ib))
+ zst = zst_bnd(ind1_bnd_gbp(ib)) * fac_bnd_gbp(ib) + zst_bnd(ind2_bnd_gbp(ib)) * (1.0 - fac_bnd_gbp(ib))
!
else
!
! Just use the value of the one boundary point
!
- zst = zst_bnd(1)
+ zst = zst_bnd(1)
!
endif
!
@@ -557,7 +712,7 @@ subroutine update_boundary_conditions(t,dt)
!
! Barometric pressure correction
!
- zst = zst + ( pavbnd - patmb(ib)) / (rhow*9.81)
+ zst = zst + (pavbnd - patmb(ib)) / (rhow * 9.81)
!
endif
!
@@ -572,7 +727,7 @@ subroutine update_boundary_conditions(t,dt)
!
! Interpolation of nearby points
!
- zig = zsit_bnd(ind1_bnd_gbp(ib))*fac_bnd_gbp(ib) + zsit_bnd(ind2_bnd_gbp(ib))*(1.0 - fac_bnd_gbp(ib))
+ zig = zsit_bnd(ind1_bnd_gbp(ib)) * fac_bnd_gbp(ib) + zsit_bnd(ind2_bnd_gbp(ib)) * (1.0 - fac_bnd_gbp(ib))
!
else
!
@@ -586,7 +741,7 @@ subroutine update_boundary_conditions(t,dt)
!
if (t < (tspinup - 1.0e-3)) then
!
- smfac = 1.0 - (t - t0)/(tspinup - t0)
+ smfac = 1.0 - (t - t0) / (tspinup - t0)
!
zs0act = zst + zsetup
call weighted_average(zini, zs0act, smfac, 1, zs0smooth)
@@ -604,14 +759,79 @@ subroutine update_boundary_conditions(t,dt)
endif
!
if (subgrid) then ! Check on waterlevels minimally equal to z_zmin
- zsb0(ib) = max(zsb0(ib), subgrid_z_zmin(nmindbnd(ib)))
- zsb(ib) = max(zsb(ib), subgrid_z_zmin(nmindbnd(ib)))
+ zsb0(ib) = max(zsb0(ib), subgrid_z_zmin(nmb))
+ zsb(ib) = max(zsb(ib), subgrid_z_zmin(nmb))
else ! Check on waterlevels minimally equal to zb
- zsb0(ib) = max(zsb0(ib), zb(nmindbnd(ib)))
- zsb(ib) = max(zsb(ib), zb(nmindbnd(ib)))
+ zsb0(ib) = max(zsb0(ib), zb(nmb))
+ zsb(ib) = max(zsb(ib), zb(nmb))
+ endif
+ !
+ elseif (kcs(nmb) == 5) then
+ !
+ ! Downstream river point
+ !
+ ! Get water levels from inside model, and adjust for slope.
+ !
+ zst = 0.0
+ sumw = 0.0
+ !
+ do iw = 1, 3
+ !
+! write(*,*)ib,iw,nm_nbr_gbp(iw, ib),w_nbr_gbp(iw, ib)
+ nm = nm_nbr_gbp(iw, ib) ! nm index of internal neighbor
+ !
+ if (nm > 0) then
+ !
+ ! Check that this point is not dry. If so, skip.
+ !
+ if (subgrid) then
+ if (zs(nm) < subgrid_z_zmin(nmb) + 0.01) cycle
+ else
+ if (zs(nm) < zb(nmb) + huthresh) cycle
+ endif
+ !
+ zst = zst + w_nbr_gbp(iw, ib) * (zs(nm) - slope_gbp(ib) * d_nbr_gbp(iw, ib))
+ sumw = sumw + w_nbr_gbp(iw, ib)
+ !
+ endif
+ !
+ enddo
+ !
+ if (sumw > 1.0e-6) then
+ !
+ zst = zst / sumw
+ !
+ else
+ !
+ ! No wet upstream point found, so set water level equal to bed level
+ !
+ zst = -99999.0
+ !
+ endif
+ !
+ ! Make sure water level is not below bed level
+ !
+ if (subgrid) then
+ zst = max(zst, subgrid_z_zmin(nmb))
+ else
+ zst = max(zst, zb(nmb))
endif
+ !
+ zsb(ib) = zst
+ zsb0(ib) = zst
+ !
+ elseif (kcs(nmb) == 6) then
+ !
+ ! Lateral coastal (Neumann) boundary
+ !
+ ! Set water level at boundary point equal to water level inside model.
+ ! No need to set zsb and zsb0, as flux for this type of boundary is solved in sfincs_momentum.f90.
+ ! Lateral boundary u/v points have kcuv=6. They are skipped in update_boundary_fluxes.
+ !
+ zs(nmb) = zs(nmi_gbp(ib)) ! nm index of internal point. Technically there can be more than one internal point. This always uses the last point that was found.
!
endif
+ !
enddo
!
end subroutine
@@ -644,16 +864,23 @@ subroutine update_boundary_fluxes(dt)
!$acc loop independent, private(ib)
do ib = 1, nkcuv2
!
- ip = index_kcuv2(ib) ! Index in uv array of kcuv=2 velocity point
+ indb = ibkcuv2(ib)
+ !
+ nmb = nmbkcuv2(ib) ! nm index of kcs=2/3/5/6 boundary point
+ !
+ if (kcs(nmb) == 6) cycle ! Lateral boundary point. Fluxes computed in sfincs_momentum.f90
!
nmi = nmikcuv2(ib) ! Index of kcs=1 point
- nmb = nmbkcuv2(ib) ! Index of kcs=2/3 boundary point
- indb = ibkcuv2(ib)
+ !
+ ip = index_kcuv2(ib) ! Index in uv array of kcuv=2 velocity point
!
zsnmi = zs(nmi) ! total water level inside model
zsnmb = zsb(indb) ! total water level at boundary
zs0nmb = zsb0(indb) ! average water level inside model (without waves)
!
+ zsnmb = zsb(indb) ! total water level at boundary
+ zs0nmb = zsb0(indb) ! average water level inside model (without waves)
+ !
if (bndtype == 1) then
!
! Weakly reflective boundary (default)
@@ -689,13 +916,13 @@ subroutine update_boundary_fluxes(dt)
!
else
!
- hnmb = max(0.5*(zsnmb + zsnmi) - zbuv(ip), huthresh)
+ hnmb = max(0.5 * (zsnmb + zsnmi) - zbuv(ip), huthresh)
zsnmb = max(zsnmb, zb(nmb))
zs0nmb = max(zs0nmb, zb(nmb))
!
endif
!
- if (hnmb0) then
+ if (nbnd > 0) then
!
- ! Update boundary conditions at boundary points
+ ! Update boundary conditions at boundary points from time series
!
call update_boundary_points(t)
!
- ! Update boundary conditions at grid points (water levels)
- !
- call update_boundary_conditions(t, dt)
- !
endif
!
+ ! Update boundary conditions at grid points (water levels)
+ !
+ call update_boundary_conditions(t, dt)
+ !
! Update boundary fluxes()
!
call update_boundary_fluxes(dt)
@@ -893,4 +1105,27 @@ subroutine weighted_average(val1,val2,fac,iopt,val3)
!
end subroutine
+
+ function find_sfincs_cell(n, m, iref) result (nm)
+ !
+ ! Find nm index for cell n, m, iref
+ !
+ use sfincs_data
+ use quadtree
+ !
+ implicit none
+ !
+ integer, intent(in) :: n
+ integer, intent(in) :: m
+ integer, intent(in) :: iref
+ integer :: nm
+ !
+ integer :: nmq
+ !
+ nmq = find_quadtree_cell_by_index(n, m, iref)
+ nm = index_sfincs_in_quadtree(nmq)
+ !
+ end function
+
+
end module
diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90
index 4d3248da..2a158fb1 100644
--- a/source/src/sfincs_continuity.f90
+++ b/source/src/sfincs_continuity.f90
@@ -328,39 +328,7 @@ subroutine compute_water_levels_subgrid(dt)
!
dvol = 0.0
!
- if (kcs(nm)==1) then
- !
- if (crsgeo) then
- !
- a = cell_area_m2(nm)
- !
- else
- !
- a = cell_area(z_flags_iref(nm))
- !
- endif
- !
- dzsdt = 0.0
- !
- if (precip) then
- !
- ! Add nett rainfall
- !
- dzsdt = dzsdt + netprcp(nm)
- !
- endif
- !
- if (use_qext) then
- !
- ! Add external source (e.g. from XMI coupling)
- !
- dzsdt = dzsdt + qext(nm)
- !
- endif
- !
- ! dvol is still in m/s, so multiply with a * dt to get m^3
- !
- dvol = dzsdt * a * dt
+ if (kcs(nm) == 1) then
!
nmd = z_index_uv_md(nm)
nmu = z_index_uv_mu(nm)
@@ -402,12 +370,12 @@ subroutine compute_water_levels_subgrid(dt)
endif
!
endif
- endif ! kcs==1
+ endif ! kcs==1
!
if (wavemaker .and. kcs(nm) == 4) then
!
! Wave maker point (seaward of wave maker)
- ! Here we use the mean flux at the location of the wave maker
+ ! Here we use the mean flux at the location of the wave maker
!
iwm = z_index_wavemaker(nm)
!
@@ -471,11 +439,51 @@ subroutine compute_water_levels_subgrid(dt)
!
endif
!
- ! We got the volume change dvol in each active cell
- ! Now update the volume and compute new water level
+ ! We got the volume change dvol in each active cell from fluxes
+ ! Now first added precip and qext
+ ! Then adjust for storage volume
+ ! Then update the volume and compute new water level
!
- if (kcs(nm) == 1 .or. kcs(nm) == 4) then
- !
+ if (kcs(nm) == 1 .or. kcs(nm) == 4) then
+ !
+ ! Obtain cell area
+ !
+ if (crsgeo) then
+ !
+ a = cell_area_m2(nm)
+ !
+ else
+ !
+ a = cell_area(z_flags_iref(nm))
+ !
+ endif
+ !
+ if (precip .or. use_qext) then
+ !
+ dzsdt = 0.0
+ !
+ if (precip) then
+ !
+ ! Add nett rainfall
+ !
+ dzsdt = dzsdt + netprcp(nm)
+ !
+ endif
+ !
+ if (use_qext) then
+ !
+ ! Add external source (e.g. from XMI coupling)
+ !
+ dzsdt = dzsdt + qext(nm)
+ !
+ endif
+ !
+ ! dzsdt is still in m/s, so multiply with a * dt to get m^3
+ !
+ dvol = dvol + dzsdt * a * dt
+ !
+ endif
+ !
if (use_storage_volume) then
!
! If water enters the cell through a point discharge, it will NOT end up in storage volume !
@@ -549,6 +557,7 @@ subroutine compute_water_levels_subgrid(dt)
!
endif
!
+ !
if (wiggle_suppression) then
!
zsderv(nm) = zs(nm) - 2 * zs11 + zs00
diff --git a/source/src/sfincs_crosssections.f90 b/source/src/sfincs_crosssections.f90
index cb9e6eaa..0e19057c 100644
--- a/source/src/sfincs_crosssections.f90
+++ b/source/src/sfincs_crosssections.f90
@@ -1,5 +1,8 @@
module sfincs_crosssections
+ use sfincs_log
+ use sfincs_error
+
contains
subroutine read_crs_file()
@@ -27,6 +30,8 @@ subroutine read_crs_file()
integer, dimension(:), allocatable :: uv_indices
integer, dimension(:), allocatable :: vertices
!
+ logical :: ok
+ !
! Read cross sections file
!
ncrs = 0
@@ -36,6 +41,8 @@ subroutine read_crs_file()
!
call write_log('Info : reading cross sections', 0)
!
+ ok = check_file_exists(crsfile, 'Cross sections file', .true.)
+ !
! First count number of polylines
!
open(500, file=trim(crsfile))
diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90
index 22de7978..b5b1ec4d 100644
--- a/source/src/sfincs_data.f90
+++ b/source/src/sfincs_data.f90
@@ -2,8 +2,9 @@ module sfincs_data
!!! Revision data
character*256 :: build_revision, build_date
- !!!!
+ !!!
!!! Time variables
+ !!!
real :: tstart_all, tfinish_all
real*4 :: dtavg
!!!
@@ -102,6 +103,10 @@ module sfincs_data
real*4 nh_tol
real*4 runup_gauge_depth
!
+ real*4 freqmininc
+ real*4 freqmaxinc
+ integer nfreqsinc
+ !
real*4 freqminig
real*4 freqmaxig
integer nfreqsig
@@ -171,6 +176,7 @@ module sfincs_data
character*256 :: wvmfile
character*256 :: qtrfile
character*256 :: volfile
+ character*256 :: bdrfile
!
character*256 :: trefstr_iso8601
character*41 :: treftimefews
@@ -230,12 +236,17 @@ module sfincs_data
logical :: spinup_meteo
logical :: snapwave
logical :: wind_reduction
- logical :: include_boundaries
+ logical :: boundaries_in_mask ! Used when kcs=2, kcs=3, kcs=5, kcs=6
+ logical :: water_level_boundaries_in_mask
+ logical :: outflow_boundaries_in_mask
+ logical :: downstream_river_boundaries_in_mask
+ logical :: neumann_boundaries_in_mask
logical :: use_uv
logical :: wavemaker
logical :: wavemaker_mobile
+ logical :: wavemaker_hinc
logical :: use_quadtree
- LOGICAL :: use_quadtree_output
+ logical :: use_quadtree_output
logical :: interpolate_zst
logical :: advection
logical :: thetasmoothing
@@ -253,6 +264,7 @@ module sfincs_data
logical :: wave_enhanced_roughness
!!!
!!! sfincs_input.f90 switches
+ !!!
integer storevelmax
integer storefluxmax
integer storevel
@@ -403,10 +415,11 @@ module sfincs_data
!
! Boundary velocity points
!
- integer*4, dimension(:), allocatable :: index_kcuv2
- integer*4, dimension(:), allocatable :: nmbkcuv2
- integer*4, dimension(:), allocatable :: nmikcuv2
- integer*1, dimension(:), allocatable :: ibuvdir
+ integer*4, dimension(:), allocatable :: index_kcuv2 ! uv index of kcuv2 point
+ integer*4, dimension(:), allocatable :: nmbkcuv2 ! nm index of boundary cell
+ integer*4, dimension(:), allocatable :: nmikcuv2 ! nm index of internal cell
+ integer*1, dimension(:), allocatable :: ibuvdir ! direction (1 is positive (boundary left or bottom), -1 is negative (boundary right or top))
+ integer*4, dimension(:), allocatable :: ibkcuv2 ! grid boundary index of boundary cell
!
! Mean velocity at boundaries
!
@@ -429,6 +442,12 @@ module sfincs_data
integer*4, dimension(:), allocatable :: wavemaker_ndm
integer*4, dimension(:), allocatable :: z_index_wavemaker
real*4 :: wavemaker_freduv
+ real*4 :: wavemaker_Tinc2ig
+ real*4 :: wavemaker_surfslope
+ real*4 :: wavemaker_hm0_ig_factor
+ real*4 :: wavemaker_hm0_inc_factor
+ real*4 :: wavemaker_gammax
+ real*4 :: wavemaker_tpmin
logical :: wavemaker_spectrum
!
! Wave maker time series
@@ -583,22 +602,23 @@ module sfincs_data
!!!
integer ntb
integer nbnd
+ integer nbdr
integer ngbnd, itbndlast, ntbnd
integer nwbnd, itwbndlast, ntwbnd
!
! Grid boundary points
!
- integer*4, dimension(:), allocatable :: ind1_bnd_gbp
- integer*4, dimension(:), allocatable :: ind2_bnd_gbp
- integer*4, dimension(:), allocatable :: ind1_cst_gbp
- integer*4, dimension(:), allocatable :: ind2_cst_gbp
- real*4, dimension(:), allocatable :: fac_bnd_gbp
- real*4, dimension(:), allocatable :: fac_cst_gbp
- real*4, dimension(:), allocatable :: zsb
- real*4, dimension(:), allocatable :: zsb0
- real*4, dimension(:), allocatable :: patmb
- integer*4, dimension(:), allocatable :: ibkcuv2
- integer*1, dimension(:), allocatable :: ibndtype
+ integer*4, dimension(:), allocatable :: ind1_bnd_gbp ! index of nearest time series point (in bnd file)
+ integer*4, dimension(:), allocatable :: ind2_bnd_gbp ! index of second nearest time series point (in bnd file)
+ real*4, dimension(:), allocatable :: fac_bnd_gbp ! weight factor for nearest neighboring time series points
+ integer, dimension(:,:), allocatable :: nm_nbr_gbp ! nm index of upstream neighbor (for downstream river boundaries)
+ real*4, dimension(:,:), allocatable :: w_nbr_gbp ! weight of upstream neighbor (for downstream river boundaries)
+ real*4, dimension(:,:), allocatable :: d_nbr_gbp ! distance to upstream neighbor (for downstream river boundaries)
+ real*4, dimension(:), allocatable :: slope_gbp ! river slope (for downstream river boundaries)
+ integer, dimension(:), allocatable :: nmi_gbp ! nm index of internal point (for lateral neumann boundary conditions)
+ real*4, dimension(:), allocatable :: zsb ! water level with waves
+ real*4, dimension(:), allocatable :: zsb0 ! water level without waves
+ real*4, dimension(:), allocatable :: patmb ! atmospheric pressure
!
! Water level boundary points
!
@@ -610,33 +630,35 @@ module sfincs_data
real*4, dimension(:), allocatable :: zst_bnd
real*4, dimension(:), allocatable :: zsit_bnd
!
+ ! Water level boundary points
+ !
+ real*4, dimension(:), allocatable :: x_bdr
+ real*4, dimension(:), allocatable :: y_bdr
+ real*4, dimension(:), allocatable :: slope_bdr
+ real*4, dimension(:), allocatable :: azimuth_bdr
+ !
! Wave boundary points
!
- real*4, dimension(:), allocatable :: x_bwv
- real*4, dimension(:), allocatable :: y_bwv
- real*4, dimension(:), allocatable :: t_bwv
- real*4, dimension(:,:), allocatable :: hs_bwv
- real*4, dimension(:,:), allocatable :: tp_bwv
- real*4, dimension(:,:), allocatable :: wd_bwv
- real*4, dimension(:), allocatable :: hst_bwv
- real*4, dimension(:), allocatable :: tpt_bwv
- real*4, dimension(:), allocatable :: l0t_bwv
- real*4, dimension(:), allocatable :: wdt_bwv
- !
- ! cst points (should remove these?)
- !
- integer :: ncst
- real*4, dimension(:), allocatable :: x_cst
- real*4, dimension(:), allocatable :: y_cst
- real*4, dimension(:), allocatable :: slope_cst
- real*4, dimension(:), allocatable :: angle_cst
- real*4, dimension(:), allocatable :: zsetup_cst
- real*4, dimension(:), allocatable :: zig_cst
- real*4, dimension(:), allocatable :: fac_bwv_cst
- integer*4, dimension(:), allocatable :: ind1_bwv_cst
- integer*4, dimension(:), allocatable :: ind2_bwv_cst
- !
- ! IG frequencies
+ !real*4, dimension(:), allocatable :: x_bwv
+ !real*4, dimension(:), allocatable :: y_bwv
+ !real*4, dimension(:), allocatable :: t_bwv
+ !real*4, dimension(:,:), allocatable :: hs_bwv
+ !real*4, dimension(:,:), allocatable :: tp_bwv
+ !real*4, dimension(:,:), allocatable :: wd_bwv
+ !real*4, dimension(:), allocatable :: hst_bwv
+ !real*4, dimension(:), allocatable :: tpt_bwv
+ !real*4, dimension(:), allocatable :: l0t_bwv
+ !real*4, dimension(:), allocatable :: wdt_bwv
+ !
+ ! Incident frequencies (for wave makers)
+ !
+ real*4, dimension(:), allocatable :: freqinc
+ real*4, dimension(:), allocatable :: costinc
+ real*4, dimension(:), allocatable :: phiinc
+ real*4, dimension(:), allocatable :: dphiinc
+ real*4 :: dfreqinc
+ !
+ ! IG frequencies (for wave makers)
!
real*4, dimension(:), allocatable :: freqig
real*4, dimension(:), allocatable :: costig
@@ -999,13 +1021,12 @@ subroutine finalize_parameters()
!
if(allocated(ind1_bnd_gbp)) deallocate(ind1_bnd_gbp)
if(allocated(ind2_bnd_gbp)) deallocate(ind2_bnd_gbp)
- if(allocated(ind1_cst_gbp)) deallocate(ind1_cst_gbp)
- if(allocated(ind2_cst_gbp)) deallocate(ind2_cst_gbp)
+ !if(allocated(ind1_cst_gbp)) deallocate(ind1_cst_gbp)
+ !if(allocated(ind2_cst_gbp)) deallocate(ind2_cst_gbp)
if(allocated(fac_bnd_gbp)) deallocate(fac_bnd_gbp)
- if(allocated(fac_cst_gbp)) deallocate(fac_cst_gbp)
+ !if(allocated(fac_cst_gbp)) deallocate(fac_cst_gbp)
if(allocated(zsb)) deallocate(zsb)
if(allocated(zsb0)) deallocate(zsb0)
- if(allocated(ibndtype)) deallocate(ibndtype)
!
! Water level boundary points
!
@@ -1019,28 +1040,28 @@ subroutine finalize_parameters()
!
! Wave boundary points
!
- if(allocated(x_bwv)) deallocate(x_bwv)
- if(allocated(y_bwv)) deallocate(y_bwv)
- if(allocated(t_bwv)) deallocate(t_bwv)
- if(allocated(hs_bwv)) deallocate(hs_bwv)
- if(allocated(tp_bwv)) deallocate(tp_bwv)
- if(allocated(wd_bwv)) deallocate(wd_bwv)
- if(allocated(hst_bwv)) deallocate(hst_bwv)
- if(allocated(tpt_bwv)) deallocate(tpt_bwv)
- if(allocated(l0t_bwv)) deallocate(l0t_bwv)
- if(allocated(wdt_bwv)) deallocate(wdt_bwv)
+ !if(allocated(x_bwv)) deallocate(x_bwv)
+ !if(allocated(y_bwv)) deallocate(y_bwv)
+ !if(allocated(t_bwv)) deallocate(t_bwv)
+ !if(allocated(hs_bwv)) deallocate(hs_bwv)
+ !if(allocated(tp_bwv)) deallocate(tp_bwv)
+ !if(allocated(wd_bwv)) deallocate(wd_bwv)
+ !if(allocated(hst_bwv)) deallocate(hst_bwv)
+ !if(allocated(tpt_bwv)) deallocate(tpt_bwv)
+ !if(allocated(l0t_bwv)) deallocate(l0t_bwv)
+ !if(allocated(wdt_bwv)) deallocate(wdt_bwv)
!
! cst points
!
- if(allocated(x_cst)) deallocate(x_cst)
- if(allocated(y_cst)) deallocate(y_cst)
- if(allocated(slope_cst)) deallocate(slope_cst)
- if(allocated(angle_cst)) deallocate(angle_cst)
- if(allocated(zsetup_cst)) deallocate(zsetup_cst)
- if(allocated(zig_cst)) deallocate(zig_cst)
- if(allocated(fac_bwv_cst)) deallocate(fac_bwv_cst)
- if(allocated(ind1_bwv_cst)) deallocate(ind1_bwv_cst)
- if(allocated(ind2_bwv_cst)) deallocate(ind2_bwv_cst)
+ !if(allocated(x_cst)) deallocate(x_cst)
+ !if(allocated(y_cst)) deallocate(y_cst)
+ !if(allocated(slope_cst)) deallocate(slope_cst)
+ !if(allocated(angle_cst)) deallocate(angle_cst)
+ !if(allocated(zsetup_cst)) deallocate(zsetup_cst)
+ !if(allocated(zig_cst)) deallocate(zig_cst)
+ !if(allocated(fac_bwv_cst)) deallocate(fac_bwv_cst)
+ !if(allocated(ind1_bwv_cst)) deallocate(ind1_bwv_cst)
+ !if(allocated(ind2_bwv_cst)) deallocate(ind2_bwv_cst)
!
! IG frequencies
!
diff --git a/source/src/sfincs_discharges.f90 b/source/src/sfincs_discharges.f90
index d2710f0b..9e666db3 100644
--- a/source/src/sfincs_discharges.f90
+++ b/source/src/sfincs_discharges.f90
@@ -1,4 +1,7 @@
module sfincs_discharges
+
+ use sfincs_log
+ use sfincs_error
contains
!
@@ -18,6 +21,8 @@ subroutine read_discharges()
!
integer isrc, itsrc, idrn, nm, m, n, stat, j, iref
!
+ logical :: ok
+ !
! Read discharge points
!
nsrc = 0
@@ -26,6 +31,8 @@ subroutine read_discharges()
itsrclast = 1
!
if (srcfile(1:4) /= 'none') then
+ !
+ ok = check_file_exists(srcfile, 'Source points file', .true.)
!
write(logstr,'(a)')'Info : reading discharges'
call write_log(logstr, 0)
@@ -56,6 +63,8 @@ subroutine read_discharges()
write(logstr,'(a)')'Info : reading drainage file'
call write_log(logstr, 0)
!
+ ok = check_file_exists(drnfile, 'Drainage points file', .true.)
+ !
open(501, file=trim(drnfile))
do while(.true.)
read(501,*,iostat = stat)dummy
@@ -88,6 +97,8 @@ subroutine read_discharges()
!
! Read discharge time series
!
+ ok = check_file_exists(disfile, 'Discharge time series file', .true.)
+ !
open(502, file=trim(disfile))
do while(.true.)
read(502,*,iostat = stat)dummy
diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90
index 67f7ed6a..36eff676 100644
--- a/source/src/sfincs_domain.f90
+++ b/source/src/sfincs_domain.f90
@@ -71,7 +71,12 @@ subroutine initialize_processes()
precip = .false.
patmos = .false.
meteo3d = .false.
- include_boundaries = .false.
+ !
+ boundaries_in_mask = .false.
+ water_level_boundaries_in_mask = .false.
+ outflow_boundaries_in_mask = .false.
+ downstream_river_boundaries_in_mask = .false.
+ neumann_boundaries_in_mask = .false.
!
if (spwfile(1:4) /= 'none' .or. wndfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then
!
@@ -767,8 +772,8 @@ subroutine initialize_mesh()
! Flags to describe uv point
!
uv_flags_iref(ip) = z_flags_iref(nm)
- uv_flags_dir(ip) = 0 ! u point
- uv_flags_type(ip) = -1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine
+ uv_flags_dir(ip) = 0 ! u point
+ uv_flags_type(ip) = -1 ! -1 is fine too coarse, 0 is normal, 1 is coarse to fine
!
endif
!
@@ -1362,8 +1367,8 @@ subroutine initialize_mesh()
m = z_index_z_m(nm)
iref = z_flags_iref(nm)
!
- z_xz(nm) = x0 + cosrot*(1.0*(m - 0.5))*dxr(iref) - sinrot*(1.0*(n - 0.5))*dyr(iref)
- z_yz(nm) = y0 + sinrot*(1.0*(m - 0.5))*dxr(iref) + cosrot*(1.0*(n - 0.5))*dyr(iref)
+ z_xz(nm) = x0 + cosrot * (1.0 * (m - 0.5)) * dxr(iref) - sinrot * (1.0 * (n - 0.5)) * dyr(iref)
+ z_yz(nm) = y0 + sinrot * (1.0 * (m - 0.5)) * dxr(iref) + cosrot * (1.0 * (n - 0.5)) * dyr(iref)
!
enddo
!
@@ -1741,13 +1746,24 @@ subroutine initialize_boundaries()
!
! First count number of boundary cells
!
+ ! kcs = 2 : water level
+ ! kcs = 3 : outflow
+ ! kcs = 5 : river
+ ! kcs = 6 : lateral (coastal, dzs/dx = 0.0), in which case we set kcuv = 6 and we do not add this point to the boundary u/v points ! We do add this point to the grid boundary points (later on).
+ !
ngbnd = 0
+ !
do nm = 1, np
- if (kcs(nm)==2 .or. kcs(nm)==3) then
+ if (kcs(nm) == 2 .or. kcs(nm) == 3 .or. kcs(nm) == 5 .or. kcs(nm) == 6) then
!
- include_boundaries = .true.
+ boundaries_in_mask = .true.
ngbnd = ngbnd + 1
!
+ if (kcs(nm) == 2) water_level_boundaries_in_mask = .true.
+ if (kcs(nm) == 3) outflow_boundaries_in_mask = .true.
+ if (kcs(nm) == 5) downstream_river_boundaries_in_mask = .true.
+ if (kcs(nm) == 6) neumann_boundaries_in_mask = .true.
+ !
endif
enddo
!
@@ -1756,7 +1772,12 @@ subroutine initialize_boundaries()
allocate(nmindbnd(ngbnd))
allocate(zsb(ngbnd))
allocate(zsb0(ngbnd))
- allocate(ibndtype(ngbnd)) ! 0 for outflow boundary, 1 for regular boundary
+ !
+ if (neumann_boundaries_in_mask) then
+ !
+ ! Need nm indices of internal points (let's do this in sfincs_boundary.f90)
+ !
+ endif
!
zsb = 0.0 ! Total water level at boundary grid point
zsb0 = 0.0 ! Filtered water level at boundary grid point
@@ -1767,35 +1788,19 @@ subroutine initialize_boundaries()
!
do nm = 1, np
!
- if (kcs(nm)==2 .or. kcs(nm)==3) then
+ if (kcs(nm) == 2 .or. kcs(nm) == 3 .or. kcs(nm) == 5 .or. kcs(nm) == 6) then
!
! This is a boundary point
!
ngbnd = ngbnd + 1
nmindbnd(ngbnd) = nm
!
- ! Determine whether this is a regular or outflow grid point
- !
- if (kcs(nm)==2) then
- !
- ! Regular boundary point
- !
- ibndtype(ngbnd) = 1
- !
- else
- !
- ! Outflow boundary point (set kcs back to 2)
- !
- ibndtype(ngbnd) = 0
-! kcs(nm) = 2
- !
- endif
endif
enddo
!
! UV boundary points
!
- nkcuv2 = 0 ! Number of points with kcuv=2
+ nkcuv2 = 0 ! Number of points with kcuv=2 or kcuv=6
kcuv = 0
!
do ip = 1, npuv
@@ -1803,25 +1808,46 @@ subroutine initialize_boundaries()
nm = uv_index_z_nm(ip)
nmu = uv_index_z_nmu(ip)
!
- if (kcs(nm)>1 .and. kcs(nmu)>1) then
+ if (kcs(nm) > 1 .and. kcs(nmu) > 1) then
!
! Two surrounding boundary points
!
- kcuv(ip) = 0 ! Is this really necessary
+ kcuv(ip) = 0 ! Is this really necessary ?
!
- elseif (kcs(nm)*kcs(nmu) == 0) then
- !
- ! One or two surrounding inactive points, so make this velocity point inactive
- !
- kcuv(ip) = 0
- !
- elseif ((kcs(nm)==1 .and. kcs(nmu)>1) .or. (kcs(nm)>1 .and. kcs(nmu)==1)) then
+ elseif (kcs(nm) * kcs(nmu) == 0) then
+ !
+ ! One or two surrounding inactive points, so make this velocity point inactive
+ !
+ kcuv(ip) = 0
+ !
+ elseif ((kcs(nm) == 1 .and. kcs(nmu) == 2) .or. (kcs(nm) == 2 .and. kcs(nmu) == 1)) then
+ !
+ ! One surrounding water level boundary point, so make this velocity point a boundary point
+ !
+ kcuv(ip) = 2
+ nkcuv2 = nkcuv2 + 1
+ !
+ elseif ((kcs(nm) == 1 .and. kcs(nmu) == 3) .or. (kcs(nm) == 3 .and. kcs(nmu) == 1)) then
+ !
+ ! One surrounding outflow boundary point, so make this velocity point a boundary point
!
- ! One surrounding boundary points, so make this velocity point a boundary point
+ kcuv(ip) = 2
+ nkcuv2 = nkcuv2 + 1
+ !
+ elseif ((kcs(nm) == 1 .and. kcs(nmu) == 5) .or. (kcs(nm) == 5 .and. kcs(nmu) == 1)) then
+ !
+ ! One surrounding downstream river boundary point, so make this velocity point a boundary point
!
kcuv(ip) = 2
nkcuv2 = nkcuv2 + 1
!
+ elseif ((kcs(nm) == 1 .and. kcs(nmu) == 6) .or. (kcs(nm) == 6 .and. kcs(nmu) == 1)) then
+ !
+ ! Lateral coastal boundary
+ !
+ kcuv(ip) = 6
+ nkcuv2 = nkcuv2 + 1
+ !
else
!
kcuv(ip) = 1
@@ -1845,7 +1871,7 @@ subroutine initialize_boundaries()
!
do ip = 1, npuv
!
- if (kcuv(ip)==2) then
+ if (kcuv(ip) == 2 .or. kcuv(ip) == 6) then
!
ikcuv2 = ikcuv2 + 1 ! Counter for kcuv==2 points
index_kcuv2(ikcuv2) = ip
@@ -1875,13 +1901,16 @@ subroutine initialize_boundaries()
!
do ib = 1, ngbnd
!
- if (nmindbnd(ib)==nmbkcuv2(ikcuv2)) then
+ if (nmindbnd(ib) == nmbkcuv2(ikcuv2)) then
!
ibkcuv2(ikcuv2) = ib ! index in grid boundary array of this uv boundary point
!
- ! Set water levels for this point (this only happens here)
+ ! Check if this is an outflow boundary (kcs=3)
!
- if (ibndtype(ib)==0) then
+ if (kcs(nmindbnd(ib)) == 3) then
+ !
+ ! Set water levels for this point (this only happens here, so these do not change throughout the simulation)
+ !
if (subgrid) then
zsb0(ib) = subgrid_z_zmin(nmindbnd(ib))
zsb(ib) = zsb0(ib)
diff --git a/source/src/sfincs_error.f90 b/source/src/sfincs_error.f90
index 36cb5975..4142b02f 100644
--- a/source/src/sfincs_error.f90
+++ b/source/src/sfincs_error.f90
@@ -1,8 +1,23 @@
module sfincs_error
-
-contains
+ !
+ use sfincs_log
+ !
+ contains
+ !
+ subroutine stop_sfincs(message, error_code)
+ !
+ implicit none
+ !
+ integer, intent(in) :: error_code
+ character(len=*), intent(in) :: message
+ !
+ write(logstr,'(a, a)')'Error : ', message
+ call write_log(logstr, 1)
+ stop error_code
+ !
+ end
!
- function check_file_exists(file_name, file_type) result(exists)
+ function check_file_exists(file_name, file_type, stop_on_error) result(exists)
!
! Checks if file exists
! Returns true or false
@@ -13,17 +28,22 @@ function check_file_exists(file_name, file_type) result(exists)
!
character(len=*) :: file_name
character(len=*) :: file_type
+ character(256) :: message
logical :: exists
+ logical :: stop_on_error
!
exists = .true.
!
- inquire( file=trim(file_name), exist=exists )
+ inquire( file=trim(file_name), exist=exists)
!
if (.not. exists) then
!
- error = 2
- !
- write(error_message,'(5a)')' Error! ', trim(file_type), ' "', trim(file_name), '" not found! Simulation has not been started!'
+ if (stop_on_error) then
+ !
+ write(message,'(4a)')trim(file_type), ' "', trim(file_name), '" not found! SFINCS has stopped!'
+ call stop_sfincs(trim(message), 2)
+ !
+ endif
!
endif
!
diff --git a/source/src/sfincs_initial_conditions.F90 b/source/src/sfincs_initial_conditions.F90
index dc826b58..7a9caa8c 100644
--- a/source/src/sfincs_initial_conditions.F90
+++ b/source/src/sfincs_initial_conditions.F90
@@ -3,6 +3,7 @@ module sfincs_initial_conditions
!
use sfincs_data
use sfincs_log
+ use sfincs_error
use netcdf
!
type net_type_ini
@@ -53,6 +54,8 @@ subroutine set_initial_conditions()
write(logstr,'(a,a)')'Info : reading restart file ', trim(rstfile)
call write_log(logstr, 0)
!
+ iok = check_file_exists(rstfile, 'Restart file', .true.)
+ !
call read_binary_restart_file() ! Note - older type real*4 for zs
!
elseif (zsinifile(1:4) /= 'none') then
@@ -63,6 +66,8 @@ subroutine set_initial_conditions()
write(logstr,'(a,a)')'Info : reading initial conditions file ', trim(zsinifile)
call write_log(logstr, 0)
!
+ iok = check_file_exists(zsinifile, 'Initial conditions file', .true.)
+ !
call read_zsini_file()
!
elseif (ncinifile(1:4) /= 'none') then
@@ -75,6 +80,8 @@ subroutine set_initial_conditions()
write(logstr,'(a,a)')'Info : reading NetCDF initial conditions file ', trim(zsinifile)
call write_log(logstr, 0)
!
+ iok = check_file_exists(ncinifile, 'NetCDF initial conditions file', .true.)
+ !
call read_nc_ini_file()
!
else
diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90
index 01c3921a..fec5349c 100644
--- a/source/src/sfincs_input.f90
+++ b/source/src/sfincs_input.f90
@@ -9,6 +9,7 @@ subroutine read_sfincs_input()
use sfincs_data
use sfincs_date
use sfincs_log
+ use sfincs_error
!
implicit none
!
@@ -35,11 +36,14 @@ subroutine read_sfincs_input()
integer iwavemaker
integer iwavemaker_spectrum
integer ispwprecip
- logical iviscosity
+ logical iviscosity
+ logical ok
!
character*256 wmsigstr
character*256 advstr
!
+ ok = check_file_exists('sfincs.inp', 'SFINCS input file', .true.)
+ !
open(500, file='sfincs.inp')
!
call read_int_input(500,'mmax',mmax,0)
@@ -85,6 +89,9 @@ subroutine read_sfincs_input()
call read_int_input(500,'nfreqsig',nfreqsig,100)
call read_real_input(500,'freqminig',freqminig,0.0)
call read_real_input(500,'freqmaxig',freqmaxig,0.1)
+ call read_int_input(500,'nfreqsinc',nfreqsinc,100)
+ call read_real_input(500,'freqmininc',freqmininc,0.04)
+ call read_real_input(500,'freqmaxinc',freqmaxinc,1.0)
call read_real_input(500,'latitude',latitude,0.0)
call read_real_input(500,'pavbnd',pavbnd,0.0)
call read_real_input(500,'gapres',gapres,101200.0)
@@ -112,7 +119,13 @@ subroutine read_sfincs_input()
call read_int_input(500,'dtoutfixed', ioutfixed, 1)
call read_real_input(500,'wmtfilter',wmtfilter,600.0)
call read_real_input(500,'wmfred',wavemaker_freduv,0.99)
- call read_char_input(500,'wmsignal',wmsigstr,'spectrum')
+ call read_char_input(500,'wmsignal',wmsigstr,'spectrum')
+ call read_real_input(500,'wavemaker_tinc2ig', wavemaker_Tinc2ig, -1.0)
+ call read_real_input(500,'wavemaker_surfslope', wavemaker_surfslope, -1.0)
+ call read_real_input(500,'wavemaker_hm0_ig_factor', wavemaker_hm0_ig_factor, 1.0)
+ call read_real_input(500,'wavemaker_hm0_inc_factor', wavemaker_hm0_inc_factor, 1.0)
+ call read_real_input(500,'wavemaker_gammax', wavemaker_gammax, 1.0)
+ call read_real_input(500,'wavemaker_tpmin', wavemaker_tpmin, 1.0)
call read_char_input(500,'advection_scheme',advstr,'upw1')
call read_real_input(500,'btrelax',btrelax,3600.0)
call read_logical_input(500,'wiggle_suppression', wiggle_suppression, .true.)
@@ -133,6 +146,7 @@ subroutine read_sfincs_input()
call read_int_input(500, 'nh_itermax', nh_itermax, 100)
call read_logical_input(500, 'h73table', h73table, .false.)
call read_real_input(500, 'rugdepth', runup_gauge_depth, 0.05)
+ call read_logical_input(500, 'wavemaker_hinc', wavemaker_hinc, .false.)
call read_logical_input(500, 'wave_enhanced_roughness', wave_enhanced_roughness, .false.)
!
! Domain
@@ -154,23 +168,24 @@ subroutine read_sfincs_input()
!
! Forcing
!
- call read_char_input(500,'bndfile',bndfile,'none')
- call read_char_input(500,'bzsfile',bzsfile,'none')
- call read_char_input(500,'bzifile',bzifile,'none')
- call read_char_input(500,'bwvfile',bwvfile,'none')
- call read_char_input(500,'bhsfile',bhsfile,'none')
- call read_char_input(500,'btpfile',btpfile,'none')
- call read_char_input(500,'bwdfile',bwdfile,'none')
- call read_char_input(500,'bdsfile',bdsfile,'none')
- call read_char_input(500,'wfpfile',wfpfile,'none')
- call read_char_input(500,'whifile',whifile,'none')
- call read_char_input(500,'wtifile',wtifile,'none')
- call read_char_input(500,'wstfile',wstfile,'none')
- call read_char_input(500,'srcfile',srcfile,'none')
- call read_char_input(500,'disfile',disfile,'none')
- call read_char_input(500,'spwfile',spwfile,'none')
- call read_char_input(500,'wndfile',wndfile,'none')
- call read_char_input(500,'prcfile',prcpfile,'none')
+ call read_char_input(500, 'bndfile', bndfile, 'none')
+ call read_char_input(500, 'bzsfile', bzsfile, 'none')
+ call read_char_input(500, 'bzifile', bzifile, 'none')
+ call read_char_input(500, 'bdrfile', bdrfile, 'none')
+ !call read_char_input(500, 'bwvfile', bwvfile, 'none')
+ !call read_char_input(500, 'bhsfile', bhsfile, 'none')
+ !call read_char_input(500, 'btpfile', btpfile, 'none')
+ !call read_char_input(500, 'bwdfile', bwdfile, 'none')
+ !call read_char_input(500, 'bdsfile', bdsfile, 'none')
+ call read_char_input(500, 'wfpfile', wfpfile, 'none')
+ call read_char_input(500, 'whifile', whifile, 'none')
+ call read_char_input(500, 'wtifile', wtifile, 'none')
+ call read_char_input(500, 'wstfile', wstfile, 'none')
+ call read_char_input(500, 'srcfile', srcfile, 'none')
+ call read_char_input(500, 'disfile', disfile, 'none')
+ call read_char_input(500, 'spwfile', spwfile, 'none')
+ call read_char_input(500, 'wndfile', wndfile, 'none')
+ call read_char_input(500, 'prcfile', prcpfile, 'none')
if (prcpfile(1:4) == 'none') then
! Try with old keyword
call read_char_input(500,'precipfile',prcpfile,'none')
@@ -204,6 +219,7 @@ subroutine read_sfincs_input()
call read_char_input(500,'netspwfile',netspwfile,'none')
!
! Output
+ !
call read_char_input(500,'obsfile',obsfile,'none')
call read_char_input(500,'crsfile',crsfile,'none')
call read_char_input(500, 'rugfile', rugfile, 'none')
diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90
index 2a5b067d..57d8b7c8 100644
--- a/source/src/sfincs_lib.f90
+++ b/source/src/sfincs_lib.f90
@@ -9,6 +9,7 @@ module sfincs_lib
use sfincs_boundaries
use sfincs_obspoints
use sfincs_crosssections
+ use sfincs_runup_gauges
use sfincs_discharges
use sfincs_meteo
use sfincs_infiltration
@@ -146,16 +147,16 @@ function sfincs_initialize() result(ierr)
!
call read_structures() ! Reads thd files and sets kcuv to zero where necessary
!
- call write_log('Reading boundary data ...', 0)
call read_boundary_data() ! Reads bnd, bzs, etc files
!
call find_boundary_indices()
!
- call write_log('Reading observation points ...', 0)
call read_obs_points() ! Reads obs file
!
call read_crs_file() ! Reads cross sections
!
+ call read_rug_file() ! Reads cross sections
+ !
call read_discharges() ! Reads dis and src file
!
if (nonhydrostatic) then
@@ -249,7 +250,7 @@ function sfincs_initialize() result(ierr)
!
call system_clock(count1, count_rate, count_max)
!
- tinput = 1.0*(count1 - count0)/count_rate
+ tinput = 1.0 * (count1 - count0) / count_rate
!
! Initialize some parameters
!
@@ -600,7 +601,7 @@ function sfincs_update(dtrange) result(ierr)
!
! Stop loop in case of instabilities (make sure time step 'dtmin' does not get too small compared to 'uvmax' flow velocity)
!
- if (dtchk1) then
+ if (dtchk < dtmin .and. nt > 1) then
!
error = 1
!
@@ -624,8 +625,10 @@ function sfincs_update(dtrange) result(ierr)
percdonenext = 1.0 * (int(percdone) + percdoneval)
!
call system_clock(count1, count_rate, count_max)
+ !
trun = 1.0*(count1 - count00)/count_rate
trem = trun / max(0.01*percdone, 1.0e-6) - trun
+ !
if (int(percdone)>0) then
write(logstr,'(i4,a,f7.1,a)')int(percdone),'% complete, ',trem,' s remaining ...'
call write_log(logstr, 1)
@@ -677,7 +680,7 @@ function sfincs_finalize() result(ierr)
write(logstr,'(a,f10.3)') ' Time in input : ',tinput
call write_log(logstr, 1)
!
- if (include_boundaries) then
+ if (boundaries_in_mask) then
write(logstr,'(a,f10.3,a,f5.1,a)') ' Time in boundaries : ',tloopbnd,' (',100*tloopbnd/(tfinish_all - tstart_all),'%)'
call write_log(logstr, 1)
endif
diff --git a/source/src/sfincs_meteo.f90 b/source/src/sfincs_meteo.f90
index 8315329b..178e87b8 100644
--- a/source/src/sfincs_meteo.f90
+++ b/source/src/sfincs_meteo.f90
@@ -2,7 +2,7 @@ module sfincs_meteo
contains
- subroutine read_meteo_data()
+ subroutine read_meteo_data()
!
! Read different meteo input types
!
@@ -10,6 +10,7 @@ subroutine read_meteo_data()
use sfincs_spiderweb
use sfincs_ncinput
use sfincs_log
+ use sfincs_error
!
implicit none
!
@@ -17,12 +18,16 @@ subroutine read_meteo_data()
!
real*4 dummy, wnd, xx, yy
!
+ logical :: ok
+ !
spw_precip = .false.
!
if (spwfile(1:4) /= 'none') then
!
write(logstr,'(a,a)')'Info : reading spiderweb file ', trim(spwfile)
call write_log(logstr, 0)
+ !
+ ok = check_file_exists(spwfile, 'Spiderweb file', .true.)
!
call read_spw_dimensions(spwfile,spw_nt,spw_nrows,spw_ncols,spw_radius,spw_nquant)
!
@@ -84,6 +89,8 @@ subroutine read_meteo_data()
write(logstr,'(a,a)')'Info : reading netcdf spiderweb file ', trim(netspwfile)
call write_log(logstr, 0)
!
+ ok = check_file_exists(netspwfile, 'Spiderweb netCDF file', .true.)
+ !
call read_netcdf_spw_data()
!
endif
@@ -110,6 +117,8 @@ subroutine read_meteo_data()
if (amufile(1:4) /= 'none') then
!
call write_log('Info : reading amu and amv file', 0)
+ !
+ ok = check_file_exists(amufile, 'AMU file', .true.)
!
call read_amuv_dimensions(amufile,amuv_nt,amuv_nrows,amuv_ncols,amuv_x_llcorner,amuv_y_llcorner,amuv_dx,amuv_dy,amuv_nquant)
!
@@ -125,6 +134,8 @@ subroutine read_meteo_data()
call read_amuv_file(amvfile,amuv_nt,amuv_nrows,amuv_ncols,amuv_times,amuv_wv,trefstr)
!
elseif (netamuamvfile(1:4) /= 'none') then ! FEWS compatible Netcdf amu&amv wind spatial input
+ !
+ ok = check_file_exists(amvfile, 'AMV file', .true.)
!
call read_netcdf_amuv_data()
!
@@ -137,6 +148,8 @@ subroutine read_meteo_data()
!
call write_log('Info : reading ampr file', 0)
!
+ ok = check_file_exists(amprfile, 'AMPR file', .true.)
+ !
call read_amuv_dimensions(amprfile,ampr_nt,ampr_nrows,ampr_ncols,ampr_x_llcorner,ampr_y_llcorner,ampr_dx,ampr_dy,ampr_nquant)
!
! Allocate
@@ -148,6 +161,8 @@ subroutine read_meteo_data()
call read_amuv_file(amprfile,ampr_nt,ampr_nrows,ampr_ncols,ampr_times,ampr_pr,trefstr)
!
elseif (netamprfile(1:4) /= 'none') then ! FEWS compatible Netcdf ampr precipitation spatial input
+ !
+ ok = check_file_exists(netamprfile, 'AMPR netCDF file', .true.)
!
call read_netcdf_ampr_data()
!
@@ -159,6 +174,8 @@ subroutine read_meteo_data()
!
call write_log('Info : reading amp file', 0)
!
+ ok = check_file_exists(ampfile, 'AMP file', .true.)
+ !
call read_amuv_dimensions(ampfile,amp_nt,amp_nrows,amp_ncols,amp_x_llcorner,amp_y_llcorner,amp_dx,amp_dy,amp_nquant)
!
! Allocate
@@ -171,6 +188,8 @@ subroutine read_meteo_data()
!
elseif (netampfile(1:4) /= 'none') then ! FEWS compatible Netcdf amp barometric pressure spatial input
!
+ ok = check_file_exists(netampfile, 'AMP netCDF file', .true.)
+ !
call read_netcdf_amp_data()
!
allocate(amp_patm01(amp_nrows, amp_ncols)) ! (has to be allocated somewhere)
@@ -180,9 +199,12 @@ subroutine read_meteo_data()
if (wndfile(1:4) /= 'none') then
!
! Wind in time series file
+ !
write(logstr,'(a,a)')'Info : reading ', trim(wndfile)
call write_log(logstr, 0)
!
+ ok = check_file_exists(wndfile, 'Wind file', .true.)
+ !
ntwnd = 0
itwndlast = 1
!
@@ -213,6 +235,8 @@ subroutine read_meteo_data()
write(logstr,'(a,a)')'Info : reading prcp file ', trim(prcpfile)
call write_log(logstr, 0)
!
+ ok = check_file_exists(prcpfile, 'Precipitation file', .true.)
+ !
ntprcp = 0
itprcplast = 1
!
diff --git a/source/src/sfincs_momentum.f90 b/source/src/sfincs_momentum.f90
index 258c0b6b..36bd831f 100644
--- a/source/src/sfincs_momentum.f90
+++ b/source/src/sfincs_momentum.f90
@@ -134,9 +134,9 @@ subroutine compute_fluxes(dt, min_dt, tloop)
!$acc loop independent, reduction( min : min_dt ), gang, vector
do ip = 1, npuv
!
- if (kcuv(ip)==1) then
+ if (kcuv(ip) == 1 .or. kcuv(ip) == 6) then
!
- ! Regular UV point
+ ! Regular UV point (or a coastal lateral boundary point)
!
! Indices of surrounding water level points
!
@@ -347,9 +347,9 @@ subroutine compute_fluxes(dt, min_dt, tloop)
!
! Interpolation required
!
- dzuv = (zmax - zmin) / (subgrid_nlevels - 1) ! level size (is storing this in memory faster?)
- iuv = min(int((zsu - zmin) / dzuv) + 1, subgrid_nlevels - 1) ! index of level below zsu
- facint = (zsu - (zmin + (iuv - 1) * dzuv) ) / dzuv ! 1d interpolation coefficient
+ dzuv = (zmax - zmin) / (subgrid_nlevels - 1) ! level size (is storing this in memory faster?)
+ iuv = min(int((zsu - zmin) / dzuv) + 1, subgrid_nlevels - 1) ! index of level below zsu
+ facint = (zsu - (zmin + (iuv - 1) * dzuv) ) / dzuv ! 1d interpolation coefficient
!
hu = subgrid_uv_havg(iuv, ip) + (subgrid_uv_havg(iuv + 1, ip) - subgrid_uv_havg(iuv, ip)) * facint ! grid-average depth
gnavg2 = subgrid_uv_nrep(iuv, ip) + (subgrid_uv_nrep(iuv + 1, ip) - subgrid_uv_nrep(iuv, ip)) * facint ! representative g*n^2
@@ -364,6 +364,14 @@ subroutine compute_fluxes(dt, min_dt, tloop)
!
endif
!
+ ! If coupled with a wave model (SnapWave or HurryWave), we may want to set an apparent roughness (which is computed elsewhere for each grid cell)
+ !
+ ! if (use_apparent_roughness) then
+ ! !
+ ! gnavg2 = 9.81 * (0.5 * (napp(nm) + napp(nmu))**2
+ ! !
+ ! endif
+ !
! Compute wet average depth hwet (used in wind and wave forcing)
!
hwet = hu / phi
@@ -575,7 +583,7 @@ subroutine compute_fluxes(dt, min_dt, tloop)
! facmax = 0.25*sqrt(g)*rhow*gammax**2
! fmax = facmax*hu*sqrt(hu)/tp/rhow (we already divided by rhow in sfincs_snapwave)
!
- fwmax = 0.8 * hwet * sqrt(hwet) / 15
+ fwmax = 0.8 * hwet * sqrt(hwet) / 6
!
frc = frc + phi * sign(min(abs(fwuv(ip)), fwmax), fwuv(ip))
!
diff --git a/source/src/sfincs_obspoints.f90 b/source/src/sfincs_obspoints.f90
index dafd6a94..ef6dc84d 100644
--- a/source/src/sfincs_obspoints.f90
+++ b/source/src/sfincs_obspoints.f90
@@ -32,11 +32,7 @@ subroutine read_obs_points()
write(logstr,'(a)')'Info : reading observation points'
call write_log(logstr, 0)
!
- ok = check_file_exists(obsfile, 'obs file')
- !
- if (.not. ok) then
- return
- endif
+ ok = check_file_exists(obsfile, 'Observation points file', .true.)
!
open(500, file=trim(obsfile))
do while(.true.)
@@ -45,6 +41,7 @@ subroutine read_obs_points()
nobs = nobs + 1
enddo
rewind(500)
+ !
allocate(xobs(nobs))
allocate(yobs(nobs))
allocate(zobs(nobs))
@@ -59,7 +56,7 @@ subroutine read_obs_points()
allocate(zbobs(nobs))
!
allocate(nmwindobs(nobs))
- allocate(wobs(4, nobs))
+ ! allocate(wobs(4, nobs))
!
allocate(value(2))
!
@@ -127,6 +124,7 @@ subroutine read_obs_points()
endif
!
iref = z_flags_iref(nm)
+ !
endif
!
write(logstr,'(a,i0,a,a,a,i0,a,i0,a,i0,a,i0,a,f0.3)')'Info : observation point ',iobs,' : "',trim(nameobs(iobs)),'" nm=',nm,' n=',n,' m=',m,' iref=',iref,' z=',zbobs(iobs)
diff --git a/source/src/sfincs_runup_gauges.f90 b/source/src/sfincs_runup_gauges.f90
index 53cdcc0a..44b208bc 100644
--- a/source/src/sfincs_runup_gauges.f90
+++ b/source/src/sfincs_runup_gauges.f90
@@ -22,6 +22,8 @@ subroutine read_rug_file()
!
character(len=256) :: cdummy
!
+ logical :: ok
+ !
! Read cross sections file
!
nr_runup_gauges = 0
@@ -30,6 +32,8 @@ subroutine read_rug_file()
!
call write_log('Info : reading run-up gauges', 0)
!
+ ok = check_file_exists(rugfile, 'Run-up gauge file', .true.)
+ !
! First count number of polylines
!
open(500, file=trim(rugfile))
diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90
index c4384b26..6b8a98b5 100644
--- a/source/src/sfincs_snapwave.f90
+++ b/source/src/sfincs_snapwave.f90
@@ -1,6 +1,7 @@
module sfincs_snapwave
!
use sfincs_log
+ use sfincs_error
!
implicit none
!
@@ -34,6 +35,7 @@ module sfincs_snapwave
integer*4, dimension(:), allocatable :: index_snapwave_in_sfincs
integer*4, dimension(:), allocatable :: index_sfincs_in_snapwave
integer*4, dimension(:), allocatable :: index_sw_in_qt ! used in sfincs_ncoutput (copy of index_snapwave_in_quadtree from snapwave_data)
+ real*4 :: snapwave_hsmean
real*4 :: snapwave_tpmean
real*4 :: snapwave_tpigmean
!
@@ -122,13 +124,14 @@ subroutine couple_snapwave(crsgeo)
snapwave_u10dir = 0.0
endif
!
+ snapwave_hsmean = 0.0
snapwave_tpmean = 0.0
snapwave_tpigmean = 0.0
-
!
call find_matching_cells(index_quadtree_in_snapwave, index_snapwave_in_quadtree)
!
! Copy final snapwave mask from snapwave_domain for output in sfincs_ncoutput
+ !
snapwave_mask = msk
!
call write_log('------------------------------------------', 1)
@@ -194,6 +197,7 @@ subroutine find_matching_cells(index_quadtree_in_snapwave, index_snapwave_in_qua
! Loop through SnapWave points
!
do ipsw = 1, snapwave_no_nodes
+ !
iq = index_quadtree_in_snapwave(ipsw)
ipsf = index_sfincs_in_quadtree(iq)
!
@@ -235,6 +239,7 @@ subroutine find_matching_cells(index_quadtree_in_snapwave, index_snapwave_in_qua
iq = index_quadtree_in_sfincs(ipsf)
ipsw = index_snapwave_in_quadtree(iq)
index_snapwave_in_sfincs(ipsf) = ipsw
+ !if (ipsf==16513) write(*,*)'snapwave',ipsw
enddo
!
! Print warning message
@@ -267,7 +272,7 @@ subroutine update_wave_field(t, tloop)
real*4, dimension(:), allocatable :: dwig0
real*4, dimension(:), allocatable :: dfig0
real*4, dimension(:), allocatable :: cg0
- real*4, dimension(:), allocatable :: qb0
+ !real*4, dimension(:), allocatable :: qb0
real*4, dimension(:), allocatable :: beta0
real*4, dimension(:), allocatable :: srcig0
real*4, dimension(:), allocatable :: alphaig0
@@ -283,7 +288,7 @@ subroutine update_wave_field(t, tloop)
allocate(dwig0(np))
allocate(dfig0(np))
allocate(cg0(np))
- allocate(qb0(np))
+ !allocate(qb0(np))
allocate(beta0(np))
allocate(srcig0(np))
allocate(alphaig0(np))
@@ -295,7 +300,7 @@ subroutine update_wave_field(t, tloop)
dwig0 = 0.0
dfig0 = 0.0
cg0 = 0.0
- qb0 = 0.0
+ !qb0 = 0.0
beta0 = 0.0
srcig0 = 0.0
alphaig0 = 0.0
@@ -306,13 +311,14 @@ subroutine update_wave_field(t, tloop)
!
ip = index_sfincs_in_snapwave(nm) ! matching index in SFINCS mesh
!
- if (ip>0) then
+ if (ip > 0) then
!
! A matching SFINCS point is found
!
+
if (wavemaker) then
!
- snapwave_depth(nm) = max(zsm(ip) - snapwave_z(nm), 0.00001)
+ snapwave_depth(nm) = max(zsm(ip) - snapwave_z(nm), 0.00001)
!
else
!
@@ -386,7 +392,7 @@ subroutine update_wave_field(t, tloop)
dwig0(nm) = snapwave_Dwig(ip)
dfig0(nm) = snapwave_Dfig(ip)
cg0(nm) = snapwave_cg(ip)
- qb0(nm) = snapwave_Qb(ip)
+ !qb0(nm) = snapwave_Qb(ip)
beta0(nm) = snapwave_beta(ip)
srcig0(nm) = snapwave_srcig(ip)
alphaig0(nm) = snapwave_alphaig(ip)
@@ -410,7 +416,7 @@ subroutine update_wave_field(t, tloop)
dwig0(nm) = 0.0
dfig0(nm) = 0.0
cg0(nm) = 0.0
- qb0(nm) = 0.0
+ !qb0(nm) = 0.0
beta0(nm) = 0.0
srcig0(nm) = 0.0
alphaig0(nm) = 0.0
@@ -430,7 +436,7 @@ subroutine update_wave_field(t, tloop)
dwig(nm) = dwig0(nm)
dfig(nm) = dfig0(nm)
cg(nm) = cg0(nm)
- qb(nm) = qb0(nm)
+ !qb(nm) = qb0(nm)
betamean(nm) = beta0(nm)
srcig(nm) = srcig0(nm)
alphaig(nm) = alphaig0(nm)
@@ -439,8 +445,8 @@ subroutine update_wave_field(t, tloop)
!
enddo
!
- hm0 = hm0*sqrt(2.0)
- hm0_ig = hm0_ig*sqrt(2.0)
+ hm0 = hm0 * sqrt(2.0)
+ hm0_ig = hm0_ig * sqrt(2.0)
!
do ip = 1, npuv
!
@@ -516,6 +522,8 @@ subroutine compute_snapwave(t)
snapwave_Fy = Fy * rho * depth
!
! 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
+ !
+ snapwave_hsmean = hsmean_bwv
snapwave_tpmean = tpmean_bwv
!
! Do quick check whether incoming Tpig value seems realistic, before using it:
@@ -550,42 +558,44 @@ subroutine read_snapwave_input()
!
! Input section
!
- call read_real_input(500,'snapwave_gamma',gamma,0.7)
- call read_real_input(500,'snapwave_gammax',gammax,0.6)
- call read_real_input(500,'snapwave_alpha',alpha,1.0)
- call read_real_input(500,'snapwave_hmin',hmin,0.1)
- call read_real_input(500,'snapwave_fw',fw0,0.01)
- call read_real_input(500,'snapwave_fwig',fw0_ig,0.015)
- call read_real_input(500,'snapwave_dt',dt,36000.0)
- call read_real_input(500,'snapwave_tol',tol,1000.0)
- call read_real_input(500,'snapwave_dtheta',dtheta,10.0)
- call read_real_input(500,'snapwave_crit',crit,0.00001) !TL: Old default was 0.01
- call read_int_input(500,'snapwave_nrsweeps',nr_sweeps,4)
- call read_int_input(500,'snapwave_niter',niter, 10) !TL: Old default was 40
- call read_int_input(500,'snapwave_baldock_opt',baldock_opt,1)
- call read_real_input(500,'snapwave_baldock_ratio',baldock_ratio,0.2)
- call read_real_input(500,'rgh_lev_land',rghlevland,0.0)
- call read_real_input(500,'snapwave_fw_ratio',fwratio,1.0)
- call read_real_input(500,'snapwave_fwig_ratio',fwigratio,1.0)
- call read_real_input(500,'snapwave_Tpini',Tpini,1.0)
- call read_int_input (500,'snapwave_mwind',mwind,2)
- call read_real_input(500,'snapwave_sigmin',sigmin,8.0*atan(1.0)/25.0)
- call read_real_input(500,'snapwave_sigmax',sigmax,8.0*atan(1.0)/1.0)
- call read_int_input (500,'snapwave_jadcgdx',jadcgdx,1)
- call read_real_input(500,'snapwave_c_dispT',c_dispT,1.0)
- call read_real_input(500,'snapwave_sector',sector,180.0)
- !
- ! Settings related to IG waves:
- call read_int_input(500,'snapwave_igwaves',igwaves_opt,1)
- call read_real_input(500,'snapwave_alpha_ig',alpha_ig,1.0) !TODO choose whether snapwave_alphaig or snapwave_gamma_ig
- call read_real_input(500,'snapwave_gammaig',gamma_ig,0.2)
- call read_real_input(500,'snapwave_shinc2ig',shinc2ig,1.0) ! Ratio of how much of the calculated IG wave source term, is subtracted from the incident wave energy (0-1, 1=default=all energy as sink)
+ call read_real_input(500, 'snapwave_gamma', gamma, 0.7)
+ call read_real_input(500, 'snapwave_gammax', gammax, 2.0) ! MvO - Changed default gammax from 0.6 to 2.0
+ call read_real_input(500, 'snapwave_alpha', alpha, 1.0)
+ call read_real_input(500, 'snapwave_hmin', hmin, 0.1)
+ call read_real_input(500, 'snapwave_fw', fw0, 0.01)
+ call read_real_input(500, 'snapwave_fwig', fw0_ig, 0.015)
+ call read_real_input(500, 'snapwave_dt', dt, 36000.0)
+ call read_real_input(500, 'snapwave_tol', tol, 1000.0)
+ call read_real_input(500, 'snapwave_dtheta', dtheta, 10.0)
+ call read_real_input(500, 'snapwave_crit', crit, 0.001) !TL: Old default was 0.01
+ call read_int_input(500, 'snapwave_nrsweeps', nr_sweeps, 4)
+ call read_int_input(500, 'snapwave_niter', niter, 10)
+ call read_int_input(500, 'snapwave_baldock_opt', baldock_opt, 1)
+ call read_real_input(500, 'snapwave_baldock_ratio', baldock_ratio, 0.2)
+ call read_real_input(500, 'rgh_lev_land', rghlevland, 0.0)
+ call read_real_input(500, 'snapwave_fw_ratio', fwratio, 1.0)
+ call read_real_input(500, 'snapwave_fwig_ratio', fwigratio, 1.0)
+ call read_real_input(500, 'snapwave_Tpini', Tpini, 1.0)
+ call read_int_input (500, 'snapwave_mwind', mwind, 2)
+ call read_real_input(500, 'snapwave_sigmin', sigmin, 8.0*atan(1.0)/25.0)
+ call read_real_input(500, 'snapwave_sigmax', sigmax, 8.0*atan(1.0)/1.0)
+ call read_int_input (500, 'snapwave_jadcgdx', jadcgdx, 1)
+ call read_real_input(500, 'snapwave_c_dispT', c_dispT, 1.0)
+ call read_real_input(500, 'snapwave_sector', sector, 180.0)
+ !
+ ! Settings related to IG waves
+ !
+ call read_int_input(500,'snapwave_igwaves', igwaves_opt, 1)
+ call read_real_input(500,'snapwave_alpha_ig', alpha_ig, 1.0) !TODO choose whether snapwave_alphaig or snapwave_gamma_ig
+ call read_real_input(500,'snapwave_gammaig', gamma_ig, 0.2)
+ call read_real_input(500,'snapwave_shinc2ig', shinc2ig, 1.0) ! Ratio of how much of the calculated IG wave source term, is subtracted from the incident wave energy (0-1, 1=default=all energy as sink)
call read_real_input(500,'snapwave_alphaigfac',alphaigfac,1.0) ! Multiplication factor for IG shoaling source/sink term
call read_real_input(500,'snapwave_baldock_ratio_ig',baldock_ratio_ig,0.2)
call read_int_input(500,'snapwave_ig_opt',ig_opt,1)
call read_int_input(500,'snapwave_iterative_srcig',iterative_srcig_opt,0) ! 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)
!
- ! IG boundary conditions options:
+ ! IG boundary conditions options
+ !
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)
call read_int_input(500,'snapwave_tpig_opt',tpig_opt,1) ! IG wave period option based on Herbers calculated spectrum, only used if snapwave_use_herbers = 1. Options are: 1=Tm01 (default), 2=Tpsmooth, 3=Tp, 4=Tm-1,0
call read_real_input(500,'snapwave_jonswapgamma',jonswapgam,3.3) ! JONSWAP gamma value for determination offshore spectrum and IG wave conditions using Herbers, default=3.3, only used if snapwave_use_herbers = 1
@@ -601,6 +611,7 @@ subroutine read_snapwave_input()
call read_int_input(500, 'vegetation', vegetation_opt, 0)
!
! Input files
+ !
call read_char_input(500,'snapwave_jonswapfile',snapwave_jonswapfile,'')
call read_char_input(500,'snapwave_bndfile',snapwave_bndfile,'')
call read_char_input(500,'snapwave_encfile',snapwave_encfile,'')
@@ -631,21 +642,21 @@ subroutine read_snapwave_input()
!
if (herbers_opt==0) then
write(logstr,*)'SnapWave: IG bc using use eeinc2ig= ',eeinc2ig,' and snapwave_Tinc2ig= ',Tinc2ig
- call write_log(logstr, 1)
+ call write_log(logstr, 0)
else
igherbers = .true.
endif
!
endif
!
- wind = .true.
+ wind = .true.
if (wind_opt==0) then
- wind = .false.
+ wind = .false.
endif
!
- vegetation = .true.
+ vegetation = .true.
if (vegetation_opt==0) then
- vegetation = .false.
+ vegetation = .false.
endif
!
if (nr_sweeps /= 1 .and. nr_sweeps /= 4) then
diff --git a/source/src/sfincs_structures.f90 b/source/src/sfincs_structures.f90
index 706c8245..01df518b 100644
--- a/source/src/sfincs_structures.f90
+++ b/source/src/sfincs_structures.f90
@@ -1,5 +1,8 @@
module sfincs_structures
+ use sfincs_error
+ use sfincs_log
+
contains
subroutine read_structures()
@@ -96,6 +99,8 @@ subroutine read_structure_file(filename,itype,npars)
write(logstr,'(a)')'Info : reading weir file'
call write_log(logstr, 0)
!
+ okay = check_file_exists(filename, 'Structures file', .true.)
+ !
! Read structures file
!
nthd = 0
@@ -391,6 +396,7 @@ subroutine read_thin_dams()
real :: dst, dstmin, xxx, yyy, dstx, dsty
real :: dummy
character :: cdummy
+ logical :: ok
!
real*4, dimension(:), allocatable :: xthd
real*4, dimension(:), allocatable :: ythd
@@ -404,6 +410,8 @@ subroutine read_thin_dams()
nthd = 0
!
if (thdfile(1:4) /= 'none') then
+ !
+ ok = check_file_exists(thdfile, 'Thin dams file', .true.)
!
! First count number of polylines
!
diff --git a/source/src/sfincs_subgrid.F90 b/source/src/sfincs_subgrid.F90
index b3ac7225..f95b73d3 100644
--- a/source/src/sfincs_subgrid.F90
+++ b/source/src/sfincs_subgrid.F90
@@ -3,6 +3,7 @@ module sfincs_subgrid
!
use netcdf
use sfincs_log
+ use sfincs_error
!
type net_type_subgrid
integer :: ncid
@@ -22,8 +23,12 @@ subroutine read_subgrid_file()
!
implicit none
!
+ logical :: ok
+ !
! Check what sort of file we're dealing with
!
+ ok = check_file_exists(sbgfile, 'Sub-grid file', .true.)
+ !
NF90(nf90_open(trim(sbgfile), NF90_CLOBBER, net_file_sbg%ncid))
!
if (net_file_sbg%ncid > 0) then
diff --git a/source/src/sfincs_wavemaker.f90 b/source/src/sfincs_wavemaker.f90
index 1867d552..e427ad71 100644
--- a/source/src/sfincs_wavemaker.f90
+++ b/source/src/sfincs_wavemaker.f90
@@ -1,5 +1,8 @@
module sfincs_wavemaker
+ use sfincs_log
+ use sfincs_error
+
contains
subroutine read_wavemaker_polylines()
@@ -63,7 +66,7 @@ subroutine read_wavemaker_polylines()
write(logstr,*)'Reading wavemaker polyline file ...'
call write_log(logstr, 0)
!
- ! Loop through all polylines
+ iok = check_file_exists(wvmfile, 'Wave maker file', .true.)
!
open(500, file=trim(wvmfile))
do while(.true.)
@@ -139,7 +142,7 @@ subroutine read_wavemaker_polylines()
!
! Check if these cells have neighbor closer to shore that is also a wavemaker point
!
- if (indwm(ip)==1) then
+ if (indwm(ip) == 1) then
!
iok = .false.
!
@@ -147,7 +150,6 @@ subroutine read_wavemaker_polylines()
!
! Check right and above
!
-! nmu = z_index_uv_mu(ip) ! index of 1st uv neighbor to the right
nmu = z_index_uv_mu1(ip) ! index of 1st uv neighbor to the right
!
if (nmu>0) then
@@ -437,7 +439,7 @@ subroutine read_wavemaker_polylines()
endif
endif
!
- if (iok) then
+ if (iok .and. kcs(ip) == 1) then
!
! This is a valid wave maker point
!
@@ -1124,7 +1126,6 @@ subroutine read_wavemaker_polylines()
!
! Set flags for kcuv points
!
- write(logstr,*)'Setting wave maker flags ...'
call write_log(logstr, 0)
!
do iwm = 1, wavemaker_nr_uv_points
@@ -1309,14 +1310,30 @@ subroutine read_wavemaker_polylines()
allocate(costig(nfreqsig))
allocate(phiig(nfreqsig))
allocate(dphiig(nfreqsig))
- dfreqig = (freqmaxig - freqminig)/nfreqsig
+ dfreqig = (freqmaxig - freqminig) / nfreqsig
do ifreq = 1, nfreqsig
- freqig(ifreq) = freqminig + ifreq*dfreqig - 0.5*dfreqig
+ freqig(ifreq) = freqminig + ifreq * dfreqig - 0.5 * dfreqig
call RANDOM_NUMBER(r)
- phiig(ifreq) = r*2*3.1416
- dphiig(ifreq) = 1.0e-6*2*3.1416/freqig(ifreq)
+ phiig(ifreq) = r * 2 * 3.1416
+ dphiig(ifreq) = 1.0e-6 * 2 * 3.1416 / freqig(ifreq)
enddo
!
+ if (wavemaker_hinc) then
+ !
+ allocate(freqinc(nfreqsinc))
+ allocate(costinc(nfreqsinc))
+ allocate(phiinc(nfreqsinc))
+ allocate(dphiinc(nfreqsinc))
+ dfreqinc = (freqmaxinc - freqmininc) / nfreqsinc
+ do ifreq = 1, nfreqsinc
+ freqinc(ifreq) = freqmininc + ifreq * dfreqinc - 0.5 * dfreqinc
+ call RANDOM_NUMBER(r)
+ phiinc(ifreq) = r * 2 * 3.1416
+ dphiinc(ifreq) = 1.0e-6 * 2 * 3.1416 / freqinc(ifreq)
+ enddo
+ !
+ endif
+ !
end subroutine
@@ -1329,12 +1346,12 @@ subroutine update_wavemaker_fluxes(t, dt, tloop)
!
implicit none
!
- integer ib, nmi, nmb, iuv, ip, ifreq, itb, itb0, itb1
- real*4 hnmb, dt, zsnmi, zsnmb, zs0nmb, zwav
+ integer ib, nmi, nmb, iuv, ip, ifreq, itb, itb0, itb1, kst
+ real*4 hnmb, dt, zsnmi, zsnmb, zs0nmb, zwav, zwav_inc
real*4 alpha, beta
real*8 t, tb
- real*4 a, fp
- real*4 tbfac, hs, tp, tpsum, setup
+ real*4 tbfac, hs, tp_ig, tp_inc, tpsum, setup, fm_ig, a, fm_inc
+ real*4 wave_steepness, betas, zinc, zig, dwvm, ztot, hm0_inc
!
real*4 ui, ub, dzuv, facint, zsuv, depthuv, uvm0
!
@@ -1352,13 +1369,13 @@ subroutine update_wavemaker_fluxes(t, dt, tloop)
!
! Interpolate boundary conditions in time
!
- if (wmf_time(1)>t - 1.0e-3) then ! use first time in boundary conditions
+ if (wmf_time(1) > t - 1.0e-3) then ! use first time in boundary conditions
!
itb0 = 1
itb1 = 1
tb = wmf_time(itb0)
!
- elseif (wmf_time(ntwmfp)t + 1.0e-6) then
+ if (wmf_time(itb) > t + 1.0e-6) then
itb0 = itb - 1
itb1 = itb
tb = t
@@ -1378,69 +1395,136 @@ subroutine update_wavemaker_fluxes(t, dt, tloop)
!
endif
!
- tbfac = (tb - wmf_time(itb0))/max(wmf_time(itb1) - wmf_time(itb0), 1.0e-6)
+ tbfac = (tb - wmf_time(itb0)) / max(wmf_time(itb1) - wmf_time(itb0), 1.0e-6)
!
tpsum = 0.0
!
do ib = 1, nwmfp ! Loop along forcing points
!
- hs = wmf_hm0_ig(ib, itb0) + (wmf_hm0_ig(ib, itb1) - wmf_hm0_ig(ib, itb0))*tbfac
- tp = wmf_tp_ig(ib, itb0) + (wmf_tp_ig(ib, itb1) - wmf_tp_ig(ib, itb0))*tbfac
- setup = wmf_setup(ib, itb0) + (wmf_setup(ib, itb1) - wmf_setup(ib, itb0))*tbfac
+ hs = wmf_hm0_ig(ib, itb0) + (wmf_hm0_ig(ib, itb1) - wmf_hm0_ig(ib, itb0)) * tbfac
+ tp_ig = wmf_tp_ig(ib, itb0) + (wmf_tp_ig(ib, itb1) - wmf_tp_ig(ib, itb0)) * tbfac
+ setup = wmf_setup(ib, itb0) + (wmf_setup(ib, itb1) - wmf_setup(ib, itb0)) * tbfac
!
wmf_hm0_ig_t(ib) = hs
wmf_setup_t(ib) = setup
!
- tpsum = tpsum + tp
+ tpsum = tpsum + tp_ig
!
enddo
!
- tp = tpsum / nwmfp ! Take average Tp from boundary points
+ tp_ig = tpsum / nwmfp ! Take average Tp from boundary points
+ tp_inc = 10.0 ! update this!
!
else
!
! Use mean peak period from SnapWave boundary conditions
!
- tp = snapwave_tpigmean ! TL: Now calculated in SnapWave, different options for using a period based on Herbers spectrum (snapwave_tpig_opt, if snapwave_use_herbers=1, or user defined snapwave_Tinc2ig ratio (if snapwave_use_herbers = 0)
+ tp_ig = snapwave_tpigmean ! TL: Now calculated in SnapWave, different options for using a period based on Herbers spectrum (snapwave_tpig_opt, if snapwave_use_herbers=1, or user defined snapwave_Tinc2ig ratio (if snapwave_use_herbers = 0)
!
+ ! We may want to use Herbers for computation of IG waves in SnapWave, but we want to have control over peak IG period at wave makers.
+ !
+ if (wavemaker_Tinc2ig > 0.0) then
+ !
+ ! Use factor on mean Tp_inc at boundaries
+ !
+ tp_ig = snapwave_tpmean * wavemaker_Tinc2ig
+ !
+ elseif (wavemaker_surfslope > 0.0) then ! Dean a
+ !
+ ! Estimate surfzone slope from Dean's a, using gambr = 1.0
+ !
+ betas = snapwave_hsmean / (snapwave_hsmean / (1.0 * wavemaker_surfslope))**(3.0 / 2.0)
+ !
+ wave_steepness = snapwave_hsmean / (1.56 * snapwave_tpmean**2)
+ !
+ ! From empirical run-up equation (van Ormondt et al., 2021), but slightly adjusted
+ !
+ tp_ig = snapwave_tpmean * max(1.86 * betas**-0.43 * wave_steepness**0.07, 5.0)
+ !
+ endif
+ !
+ tp_inc = max(snapwave_tpmean, wavemaker_tpmin)
+ !
endif
!
- alpha = min(dt/wmtfilter, 1.0)
- beta = min(dt/(0.2*wmtfilter), 1.0)
+ ! Factors for double-exponential filtering
+ !
+ alpha = min(dt / wmtfilter, 1.0)
+ beta = min(dt / (0.2 * wmtfilter), 1.0)
+ !
+ zwav = 0.0
+ zwav_inc = 0.0
!
if (wavemaker_spectrum) then
!
- ! First update phases
- !
- do ifreq = 1, nfreqsig
- phiig(ifreq) = phiig(ifreq) + dphiig(ifreq)*dt
- costig(ifreq) = cos(2*pi*t*freqig(ifreq) + phiig(ifreq))
- enddo
- !
- fp = 1.0/tp ! Wave period
+ fm_ig = 1.0 / tp_ig ! Wave period
!
! Now spectrum and wave excitation
!
- zwav = 0.0
- !
do ifreq = 1, nfreqsig
- a = 0.125*(1.0**2)*(fp**(-2))*freqig(ifreq)*exp(-freqig(ifreq)/fp)
- zwav = zwav + costig(ifreq)*sqrt(a*dfreqig)
+ !
+ ! Update phase
+ !
+ phiig(ifreq) = phiig(ifreq) + dphiig(ifreq) * dt
+ costig(ifreq) = cos(2 * pi * t * freqig(ifreq) + phiig(ifreq))
+ !
+ ! Use this spectral shape instead
+ !
+ a = 0.125 * (fm_ig**-2) * freqig(ifreq) * (exp(-freqig(ifreq) / fm_ig))
+ !
+ zwav = zwav + costig(ifreq) * sqrt(a * dfreqig)
+ !
enddo
!
+ if (wavemaker_hinc) then
+ !
+ fm_inc = 1.0 / tp_inc ! Wave period
+ !
+ do ifreq = 1, nfreqsig
+ !
+ phiinc(ifreq) = phiinc(ifreq) + dphiinc(ifreq) * dt
+ costinc(ifreq) = cos(2 * pi * t * freqinc(ifreq) + phiinc(ifreq))
+ !
+ ! The ISSC spectrum (also known as Bretschneider or modified Pierson-Moskowitz)
+ !
+ a = 0.625 * (fm_inc**4) * (freqinc(ifreq)**-5) * (exp(-1.25 * (freqinc(ifreq) / fm_inc)**-4))
+ !
+ zwav_inc = zwav_inc + costinc(ifreq) * sqrt(a * dfreqinc)
+ !
+ enddo
+ !
+ !zwav_inc = 0.5 * sin(2 * pi * t / tp_inc)
+ !
+ ! Saw tooth
+ !
+ !zwav_inc = - mod(t, tp_inc) / tp_inc + 0.5
+ !
+ ! Let zwav_inc be modulated by zwav_ig (i.e. higher incident waves at the peaks of the IG wave)
+ !
+ !zwav_inc = zwav_inc * sqrt(max(zwav + 1.0, 0.0)) ! this assumes zwav is somewhere between -0.5 and +0.5
+ !
+ endif
+ !
else
!
! Monochromatic signal
!
- zwav = 0.5*sin(2*pi*t/tp)
+ zwav = 0.5 * sin(2 * pi * t / tp_ig)
+ !
+ if (wavemaker_hinc) then
+ !
+ zwav_inc = 0.5 * sin(2 * pi * t / tp_inc)
+ !
+ endif
!
endif
!
if (t=subgrid_uv_zmax(ip) - 1.0e-3) then
+ if (zsuv >= subgrid_uv_zmax(ip) - 1.0e-3) then
!
! Entire cell is wet, no interpolation from table needed
!
depthuv = subgrid_uv_havg_zmax(ip) + zsuv
!
- elseif (zsuv>subgrid_uv_zmin(ip)) then
+ elseif (zsuv > subgrid_uv_zmin(ip)) then
!
! Interpolation required
!
dzuv = (subgrid_uv_zmax(ip) - subgrid_uv_zmin(ip)) / (subgrid_nlevels - 1)
- iuv = int((zsuv - subgrid_uv_zmin(ip))/dzuv) + 1
- facint = (zsuv - (subgrid_uv_zmin(ip) + (iuv - 1)*dzuv) ) / dzuv
- depthuv = subgrid_uv_havg(iuv, ip) + (subgrid_uv_havg(iuv + 1, ip) - subgrid_uv_havg(iuv, ip))*facint
+ iuv = int((zsuv - subgrid_uv_zmin(ip)) / dzuv) + 1
+ facint = (zsuv - (subgrid_uv_zmin(ip) + (iuv - 1) * dzuv) ) / dzuv
+ depthuv = subgrid_uv_havg(iuv, ip) + (subgrid_uv_havg(iuv + 1, ip) - subgrid_uv_havg(iuv, ip)) * facint
!
else
!
@@ -1503,7 +1608,7 @@ subroutine update_wavemaker_fluxes(t, dt, tloop)
!
endif
!
- hnmb = max(depthuv, huthresh)
+ hnmb = depthuv
zsnmb = max(zsnmb, subgrid_z_zmin(nmb))
zs0nmb = max(zs0nmb, subgrid_z_zmin(nmb))
!
@@ -1517,7 +1622,7 @@ subroutine update_wavemaker_fluxes(t, dt, tloop)
!
! Weakly reflective boundary
!
- if (hnmb 25.0)) then
- write(logstr,*)'DEBUG SnapWave - input wave height is outside acceptable range of 0-25 m: ',hs, ' and is therefore limited back to this range, please check whether your input is realistic!'
- call write_log(logstr, 1)
- !
- hs = max(min(hs, 25.0), 0.0)
+ !
+ write(logstr,*)'DEBUG SnapWave - input wave height is outside acceptable range of 0-25 m: ',hs, ' and is therefore limited back to this range, please check whether your input is realistic!'
+ call write_log(logstr, 1)
+ !
+ hs = max(min(hs, 25.0), 0.0)
+ !
endif
!
! Limit directional spreading (1 < ds < 60)
!
if ((dsp < 3*pi/180) .or. (dsp > 90*pi/180)) then
- write(logstr,*)'DEBUG SnapWave - input wave spreading is outside acceptable range of 3-90 degrees: ',dsp/pi*180, ' and is therefore limited back to this range, please check whether your input is realistic!'
- call write_log(logstr, 1)
- !
- dsp = max(min(dsp, 90*pi/180), 3*pi/180)
+ !
+ write(logstr,*)'DEBUG SnapWave - input wave spreading is outside acceptable range of 3-90 degrees: ',dsp/pi*180, ' and is therefore limited back to this range, please check whether your input is realistic!'
+ call write_log(logstr, 1)
+ !
+ dsp = max(min(dsp, 90*pi/180), 3*pi/180)
+ !
endif
!
! Limit period (0.1 < tps < 25)
!
if ((tps < 0.1) .or. (tps > 25.0)) then
- write(logstr,*)'DEBUG SnapWave - input wave period is outside acceptable range of 0.1-25 s: ',tps, ' and is therefore limited back to this range, please check whether your input is realistic!'
- call write_log(logstr, 1)
- !
- tps = max(min(tps, 25.0), 0.1)
+ !
+ write(logstr,*)'DEBUG SnapWave - input wave period is outside acceptable range of 0.1-25 s: ',tps, ' and is therefore limited back to this range, please check whether your input is realistic!'
+ call write_log(logstr, 1)
+ !
+ tps = max(min(tps, 25.0), 0.1)
+ !
endif
!
call weighted_average(wd_bwv(ib, itb0), wd_bwv(ib, itb1), 1.0 - tbfac, 2, wd) !wavdir
@@ -594,48 +601,19 @@ subroutine update_boundary_points(t)
tpt_bwv(ib) = tps
wdt_bwv(ib) = wd
dst_bwv(ib) = dsp
- !zst_bwv(ib) = zst
!
enddo
!
-! do itb = itwbndlast, ntwbnd ! Loop in time
-! !
-! if (t_bwv(itb)>t) then
-! !
-! tbfac = (t - t_bwv(itb - 1))/(t_bwv(itb) - t_bwv(itb - 1))
-! !
-! do ib = 1, nwbnd ! Loop along boundary points
-! !
-! hs = hs_bwv(ib, itb - 1) + (hs_bwv(ib, itb) - hs_bwv(ib, itb - 1))*tbfac
-! tps = tp_bwv(ib, itb - 1) + (tp_bwv(ib, itb) - tp_bwv(ib, itb - 1))*tbfac
-! dsp = ds_bwv(ib, itb - 1) + (ds_bwv(ib, itb) - ds_bwv(ib, itb - 1))*tbfac !dirspr
-! zst = zs_bwv(ib, itb - 1) + (zs_bwv(ib, itb) - zs_bwv(ib, itb - 1))*tbfac
-! !
-! call weighted_average(wd_bwv(ib, itb - 1), wd_bwv(ib, itb), 1.0 - tbfac, 2, wd) !wavdir
-! !
-! hst_bwv(ib) = hs
-! tpt_bwv(ib) = tps
-! wdt_bwv(ib) = wd
-! dst_bwv(ib) = dsp
-! zst_bwv(ib) = zst
-! !
-! enddo
-! !
-! itwbndlast = itb
-! exit
-! !
-! endif
-! enddo
- !
! Now generate wave spectra at the boundary points
!
! Average wave period and direction to determine theta grid
!
- tpmean_bwv = sum(tpt_bwv)/size(tpt_bwv)
+ tpmean_bwv = sum(tpt_bwv) / size(tpt_bwv)
+ hsmean_bwv = sum(hst_bwv) / size(hst_bwv)
!zsmean_bwv = sum(zst_bwv)/size(zst_bwv)
!depth = max(zsmean_bwv - zb,hmin) ! TL: For SFINCS we don't want this, because it overrides the real updated water depth we're inserting
depth = max(depth,hmin)
- wdmean_bwv = atan2(sum(sin(wdt_bwv)*hst_bwv)/sum(hst_bwv), sum(cos(wdt_bwv)*hst_bwv)/sum(hst_bwv))
+ wdmean_bwv = atan2(sum(sin(wdt_bwv) * hst_bwv) / sum(hst_bwv), sum(cos(wdt_bwv) * hst_bwv) / sum(hst_bwv))
!
! Determine IG boundary conditions
!
@@ -643,7 +621,8 @@ subroutine update_boundary_points(t)
!
if (igherbers) then
!
- ! Get local water depth at boundary points (can change in time)
+ ! Get local water depth at boundary points (can change in time)
+ !
call find_nearest_depth_for_boundary_points() ! Output is: deptht_bwv
!
do ib = 1, nwbnd ! Loop along boundary points
@@ -674,12 +653,26 @@ subroutine update_boundary_points(t)
!
thetamean = wdmean_bwv
!
- if (ntwbnd > 0) then
- call make_theta_grid(wdmean_bwv)
- else
- thetamean=u10dmean
- call make_theta_grid(u10dmean)
- endif
+ ind = nint(thetamean / dtheta) + 1
+ !
+ do itheta = 1, ntheta
+ ! i360(itheta) = mod2(itheta + ind - 10, 36)
+ i360(itheta) = mod2(itheta + ind - (1 + ntheta / 2), ntheta * 2)
+ enddo
+ !
+ do itheta = 1, ntheta
+ !
+ theta(itheta) = theta360(i360(itheta))
+ !
+ do k = 1, no_nodes
+ w(1, itheta, k) = w360(1, i360(itheta), k)
+ w(2, itheta, k) = w360(2, i360(itheta), k)
+ prev(1, itheta, k) = prev360(1, i360(itheta), k)
+ prev(2, itheta, k) = prev360(2, i360(itheta), k)
+ ds(itheta, k) = ds360(i360(itheta), k)
+ enddo
+ !
+ enddo
!
! Build spectra on wave boundary support points
!
@@ -695,6 +688,7 @@ subroutine update_boundary_points(t)
enddo
!
! Build IG spectra on wave boundary support points
+ !
if (igwaves) then
if (igherbers) then
do ib = 1, nwbnd ! Loop along boundary points
diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90
index 670fd0b5..f6b8561a 100644
--- a/source/src/snapwave/snapwave_data.f90
+++ b/source/src/snapwave/snapwave_data.f90
@@ -62,6 +62,7 @@ module snapwave_data
real*4, dimension(:), allocatable :: Hmx_ig
real*4, dimension(:,:), allocatable :: ee ! directional energy density
real*4, dimension(:,:), allocatable :: ee_ig ! directional infragravity energy density
+ real*4, dimension(:), allocatable :: DoverE
!
real*4, dimension(:,:), allocatable :: aa ! directional action density
real*4, dimension(:), allocatable :: sig ! mean frequency
@@ -100,6 +101,7 @@ module snapwave_data
!
integer :: nwbnd ! number of support points wave boundary
integer :: ntwbnd ! number of time points wave boundary
+ real*4 :: hsmean_bwv ! mean wave height over boundary points for given time
real*4 :: tpmean_bwv ! mean tp over boundary points for given time
real*4 :: wdmean_bwv ! mean wave direction for given time, used to make theta grid
real*4 :: zsmean_bwv ! mean water level for given time, used to make theta grid
diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90
index 54dc1040..bcf595e0 100644
--- a/source/src/snapwave/snapwave_domain.f90
+++ b/source/src/snapwave/snapwave_domain.f90
@@ -46,9 +46,13 @@ subroutine initialize_snapwave_domain()
! Quadtree file
!
if (coupled_to_sfincs) then
+ !
call read_snapwave_quadtree_mesh(.false.) ! no need to read qtr file again
+ !
else
+ !
call read_snapwave_quadtree_mesh(.true.)
+ !
endif
!
elseif (ext=='tx') then
@@ -110,6 +114,7 @@ subroutine initialize_snapwave_domain()
allocate(Cg(no_nodes))
allocate(sinhkh(no_nodes))
allocate(Hmx(no_nodes))
+ allocate(DoverE(no_nodes))
allocate(fw(no_nodes))
allocate(H(no_nodes))
allocate(Tp(no_nodes))
@@ -207,6 +212,7 @@ subroutine initialize_snapwave_domain()
WsorA = 0.0
SwE = 0.0
SwA = 0.0
+ DoverE = 0.0
windspreadfac = 0.0
!
generate_upw = .true.
@@ -304,26 +310,9 @@ subroutine initialize_snapwave_domain()
if (any(msk == 3)) then
!
! We already have all msk=3 Neumann points, now find each their nearest cell 'neumannconnected' using new 'neuboundaries_light'
- call neuboundaries_light(x,y,msk,no_nodes,tol,neumannconnected)
- !
- if (ANY(neumannconnected > 0)) then
- !
- write(logstr,*)'SnapWave: Neumann connected boundaries found ...'
- call write_log(logstr, 0)
- !
- do k=1,no_nodes
- if (neumannconnected(k)>0) then
- if (msk(k)==1) then
- ! k is inner and can be neumannconnected
- inner(neumannconnected(k))= .false.
- msk(neumannconnected(k)) = 3 !TL: should already by 3, but left it like in SnapWave SVN
- else
- ! we don't allow neumannconnected links if the node is an open boundary
- neumannconnected(k) = 0
- endif
- endif
- enddo
- endif
+ !
+ call neuboundaries_light(x, y, msk, no_nodes, neumannconnected)
+ !
else
!
neumannconnected = 0
@@ -387,7 +376,7 @@ subroutine find_upwind_neighbours(x,y,no_nodes,sferic,theta,ntheta,kp,np,w,prev,
pi = 4*atan(1.0)
!
do k = 1, no_nodes
- call findloc(kp(:,k), np, 0, nploc)
+ call findloc(kp(:,k), np, 0, nploc)
nploc = nploc - 1
if (kp(1,k)/=0) then
do itheta = 1, ntheta
@@ -402,6 +391,7 @@ subroutine find_upwind_neighbours(x,y,no_nodes,sferic,theta,ntheta,kp,np,w,prev,
xsect=[x(ind1)-x(k), x(ind2)-x(k)]*circumf_eq/360.0*cos(y(k)*180.0/pi)
ysect=[y(ind1)-y(k), y(ind2)-y(k)]*circumf_pole/360.0
call intersect_angle(0.d0, 0.d0, theta(itheta) + pi, xsect, ysect, ww, dss, xi, yi)
+
endif
if (dss/=0) then
w(1, itheta,k) = ww(1)
@@ -602,7 +592,7 @@ subroutine fm_surrounding_points(xn,yn,zn,no_nodes,sferic,face_nodes,no_faces,kp
dhdx(kn)=0.
dhdy(kn)=0.
endif
- !write(*,*)'fmgsp 2',kn, no_nodes,isp,isp2
+ !
enddo
!
end subroutine
@@ -797,81 +787,58 @@ subroutine boundaries(x,y,no_nodes,xb,yb,nb,tol,bndpts,nobndpts,bndindx,bndweigh
end subroutine boundaries
- subroutine neuboundaries_light(x,y,msk,no_nodes,tol,neumannconnected)
+ subroutine neuboundaries_light(x, y, msk, no_nodes, neumannconnected)
!
! TL: Based on subroutine find_nearest_depth_for_boundary_points of snapwave_boundaries.f90
!
implicit none
!
integer, intent(in) :: no_nodes
- real*8, dimension(no_nodes), intent(in) :: x,y
+ real*8, dimension(no_nodes), intent(in) :: x, y
integer*1, dimension(no_nodes), intent(in) :: msk
- real*4, intent(in) :: tol
integer, dimension(no_nodes), intent(out) :: neumannconnected
!
- real*4 :: h1, h2, fac
- !
- real xgb, ygb, dst1, dst2, dst
- integer k, ib1, ib2, ic, kmin
- !
- ! Loop through all msk=3 cells
- !
- do ic = 1, no_nodes
- ! Loop through all grid points
- !
- if (msk(ic)==3) then ! point ic is on the neumann boundary
- !
- dst1 = tol
- dst2 = tol
- ib1 = 0
- ib2 = 0
- !
- do k = 1, no_nodes
- !
- if (msk(k)==1) then
- xgb = x(k)
- ygb = y(k)
- !
- dst = sqrt((x(ic) - xgb)**2 + (y(ic) - ygb)**2)
- !
- if (dst 0) .and. (ib2 > 0) ) then
- !
- ! Determine the index of the minimum value, if points found within 'tol' distance
- !
- if (dst1 < dst2) then
- kmin = ib1
- else
- kmin = ib2
- endif
- !
- neumannconnected(kmin)=ic
- !
- !write(*,*)kmin,ic
- !
- endif
- !
- endif
- enddo
+ real xgb, ygb, dst1, dst
+ integer k, ib1, ic
+ !
+ ! Loop through all msk=3 cells
+ !
+ do ic = 1, no_nodes
+ !
+ ! Loop through all grid points
+ !
+ if (msk(ic) == 3) then ! point ic is on the neumann boundary
+ !
+ dst1 = 1.0e9
+ ib1 = 0
+ !
+ do k = 1, no_nodes
+ !
+ if (msk(k) == 1) then
+ !
+ xgb = x(k)
+ ygb = y(k)
+ !
+ dst = sqrt((x(ic) - xgb)**2 + (y(ic) - ygb)**2)
+ !
+ if (dst < dst1) then
+ !
+ ! Nearest point found
+ !
+ dst1 = dst
+ ib1 = k
+ !
+ endif
+ endif
+ enddo
+ !
+ if (ib1 > 0) then
+ neumannconnected(ic) = ib1
+ endif
+ !
+ endif
+ !
+ enddo
!
end subroutine neuboundaries_light
@@ -1144,13 +1111,19 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
integer :: n
integer :: nu1
integer :: nu2
+ integer :: nd1
+ integer :: nd2
integer :: m
integer :: mu1
integer :: mu2
+ integer :: md1
+ integer :: md2
integer :: mnu1
integer :: nra
integer*1 :: mu
integer*1 :: nu
+ integer*1 :: md
+ integer*1 :: nd
integer*1 :: mnu
!
logical :: load_quadtree
@@ -1191,6 +1164,7 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
allocate(msk_tmp(quadtree_nr_points)) ! Make temporary mask with all quadtree points
allocate(msk_tmp2(quadtree_nr_points)) ! Make second temporary mask with all quadtree points
allocate(zb_tmp(quadtree_nr_points)) ! Make temporary array with bed level on all quadtree points
+ !
zb_tmp = -10.0
!
msk_tmp = 1 ! Without mask file, all points will be active
@@ -1256,7 +1230,9 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
!
! STEP 4 - Make faces
!
- allocate(faces(4, 4*quadtree_nr_points)) ! max 4 nodes per faces, and max 4 faces per node
+ allocate(faces(4, 4 * quadtree_nr_points)) ! max 4 nodes per faces, and max 4 faces per node
+ !
+ faces = 0
!
nfaces = 0
!
@@ -1465,7 +1441,6 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
!
! Type 1 - Most normal cell possible
!
-! write(*,'(a,20i6)')'ip,mu,nu,mnu,mu1,nu1,mnu1',ip,mu,nu,mnu,mu1,nu1,mnu1
if (mu1>0 .and. nu1>0 .and. mnu1>0) then
nfaces = nfaces + 1
faces(1, nfaces) = ip
@@ -2034,7 +2009,80 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
msk_tmp2(nu1) = 1
endif
endif
- endif
+ endif
+ !
+ ! Add triangles around stair case boundaries
+ !
+ ! Inactive cell top left
+ !
+ mu = quadtree_mu(ip)
+ mu1 = quadtree_mu1(ip)
+ mu2 = quadtree_mu2(ip)
+ md = quadtree_md(ip)
+ md1 = quadtree_md1(ip)
+ md2 = quadtree_md2(ip)
+ nu = quadtree_nu(ip)
+ nu1 = quadtree_nu1(ip)
+ nu2 = quadtree_nu2(ip)
+ nd = quadtree_nd(ip)
+ nd1 = quadtree_nd1(ip)
+ nd2 = quadtree_nd2(ip)
+ !
+ ! Check for inactive cell top left
+ !
+ if (md == 0 .and. nu == 0 .and. md1 > 0 .and. nu1 > 0) then
+ if (msk_tmp(md1) == 2 .and. msk_tmp(nu1) == 2 .and. msk_tmp(ip) == 1) then
+ !
+ nfaces = nfaces + 1
+ faces(1, nfaces) = md1
+ faces(2, nfaces) = ip
+ faces(3, nfaces) = mu1
+ !
+ endif
+ !
+ endif
+ !
+ ! Check for inactive cell bottom left
+ !
+ if (md == 0 .and. nd == 0 .and. md1 > 0 .and. nd1 > 0) then
+ if (msk_tmp(md1) == 2 .and. msk_tmp(nd1) == 2 .and. msk_tmp(ip) == 1) then
+ !
+ nfaces = nfaces + 1
+ faces(1, nfaces) = md1
+ faces(2, nfaces) = nd1
+ faces(3, nfaces) = ip
+ !
+ endif
+ !
+ endif
+ !
+ ! Check for inactive cell top right
+ !
+ if (mu == 0 .and. nu == 0 .and. mu1 > 0 .and. nu1 > 0) then
+ if (msk_tmp(mu1) == 2 .and. msk_tmp(nu1) == 2 .and. msk_tmp(ip) == 1) then
+ !
+ nfaces = nfaces + 1
+ faces(1, nfaces) = ip
+ faces(2, nfaces) = mu1
+ faces(3, nfaces) = nu1
+ !
+ endif
+ !
+ endif
+ !
+ ! Check for inactive cell bottom right
+ !
+ if (mu == 0 .and. nd == 0 .and. mu1 > 0 .and. nd1 > 0) then
+ if (msk_tmp(mu1) == 2 .and. msk_tmp(nd1) == 2 .and. msk_tmp(ip) == 1) then
+ !
+ nfaces = nfaces + 1
+ faces(1, nfaces) = ip
+ faces(2, nfaces) = nd1
+ faces(3, nfaces) = mu1
+ !
+ endif
+ !
+ endif
!
enddo
!
@@ -2099,7 +2147,6 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
!
! Set node values
!
-! zb(nac) = zb_tmp(ip)
zb(nac) = quadtree_zz(ip)
x(nac) = quadtree_xz(ip)
y(nac) = quadtree_yz(ip)
@@ -2115,7 +2162,7 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
!
do iface = 1, no_faces
do j = 1, 4
- if (faces(j, iface)>0) then
+ if (faces(j, iface) > 0) then
ip0 = faces(j, iface) ! index in full quadtree
ip1 = index_snapwave_in_quadtree(ip0) ! index in reduced quadtree
face_nodes(j, iface) = ip1 ! set index to that of reduced mesh
diff --git a/source/src/snapwave/snapwave_infragravity.f90 b/source/src/snapwave/snapwave_infragravity.f90
index f38ce782..cf6eea97 100644
--- a/source/src/snapwave/snapwave_infragravity.f90
+++ b/source/src/snapwave/snapwave_infragravity.f90
@@ -37,37 +37,49 @@ subroutine determine_ig_bc(x_bwv, y_bwv, hsinc, tpinc, ds, jonswapgam, depth, Ti
correctHm0 = 0 ! Choice between correcting Hm0 in build_jonswap if build 2D Vardens spectrum too low (1) or not (default, 0)
! TL: Question - should we make this user definable?
!
- pi = 4.*atan(1.)
+ pi = 4.0 * atan(1.0)
!
! Convert wave spreading in degrees (input) to S
- scoeff = (2/ds**2) - 1
+ !
+ scoeff = (2 / ds**2) - 1.0
!
! Call function that calculates Hig0 following Herbers, as also implemented in XBeach and secordspec2 in Matlab
! Loosely based on 3 step calculation in waveparams.F90 of XBeach (build_jonswap, build_etdir, build_boundw), here all in 1 subroutine calculate_herbers
!
+ !
if (depth < 5.0) then
- write(logstr,*)'ERROR SnapWave - depth at boundary input point ',x_bwv, y_bwv,' dropped below 5 m: ',depth, ' which might lead to large values of Hm0ig as bc, especially when directional spreading is low! Please specify input in deeper water. '
- call write_log(logstr, 1)
- !
- write(logstr,*)'Depth set back to 5 meters for stability, simulation will continue.'
- call write_log(logstr, 1)
- !
- depth = 5.0
+ !
+ ! MvO - Is this really the best check? Does not seem very non-dimensional. Maybe something like hsig / depth > 0.5 ?
+ !
+ write(logstr, *)'ERROR SnapWave - depth at boundary input point ', x_bwv, ',', y_bwv,' below 5 m: ', depth
+ call write_log(logstr, 0)
+ write(logstr, *)'This may lead to large values of Hm0ig as bc, especially when directional spreading is low! Please specify input in deeper water.'
+ call write_log(logstr, 0)
+ write(logstr,*)'Depth set back to 5 meters for stability, simulation will continue.'
+ call write_log(logstr, 0)
+ !
+ depth = 5.0
+ !
endif
!
call compute_herbers(hsig, Tm01, Tm10, Tp, Tpsmooth, hsinc, tpinc, scoeff, jonswapgam, depth, correctHm0) ![out,out,out,out,out, in,in,in,in,in,in]
!
! Catch NaN values (if depth=0 probably) or unrealistically large values above 2 meters
if (hsig < 0.0) then
+ !
write(logstr,*)'DEBUG SnapWave - computed hm0ig at boundary dropped below 0 m: ',hsig, ' and is therefore limited back to 0 m!'
- call write_log(logstr, 1)
+ call write_log(logstr, 1)
+ !
hsig = max(hsig, 0.0)
+ !
endif
+ !
if (hsig > 3.0) then
+ !
write(logstr,*)'DEBUG SnapWave - computed hm0ig at boundary exceeds 3 meter: ',hsig, ' - please check whether this might be realistic!'
- call write_log(logstr, 1)
-
- !write(*,*)'DEBUG - computed hm0ig at boundary exceeds 3 meter: ',hsig, ' and is therefore limited back to 3 m!'
+ call write_log(logstr, 1)
+ !
+ !write(*,*)'DEBUG - computed hm0ig at boundary exceeds 3 meter: ',hsig, ' and is therefore limited back to 3 m!'
!hsig = min(hsig, 3.0)
endif
!
@@ -93,12 +105,12 @@ subroutine determine_ig_bc(x_bwv, y_bwv, hsinc, tpinc, ds, jonswapgam, depth, Ti
endif
!
! Check on ratio tpig/tpinc whether it is deemed realistic
- if (tpig/tpinc < 2.0) then
+ if (tpig / tpinc < 2.0) then
write(logstr,*)'DEBUG SnapWave - computed tpig/tpinc ratio at offshore boundary dropped below 2 and might be unrealistic! value: ',tpig/tpinc
- call write_log(logstr, 1)
+ call write_log(logstr, 1)
elseif (tpig/tpinc > 20.0) then
write(logstr,*)'DEBUG SnapWave - computed tpig/tpinc ratio at offshore boundary increased above 20 and might be unrealistic! value: ',tpig/tpinc
- call write_log(logstr, 1)
+ call write_log(logstr, 1)
endif
!
end subroutine
diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90
index d8c05226..7750123a 100644
--- a/source/src/snapwave/snapwave_solver.f90
+++ b/source/src/snapwave/snapwave_solver.f90
@@ -3,6 +3,7 @@ module snapwave_solver
use sfincs_log
implicit none
+
contains
subroutine compute_wave_field()
@@ -13,10 +14,7 @@ subroutine compute_wave_field()
!
!real*8, intent(in) :: time > TL: not used in this implementation
!
- real*4 :: tpb
- !
real*4, parameter :: waveps = 1e-5
- !real*4, dimension(:), allocatable :: sig
real*4, dimension(:), allocatable :: sigm_ig
real*4, dimension(:), allocatable :: expon
!
@@ -27,7 +25,7 @@ subroutine compute_wave_field()
allocate(sigm_ig(no_nodes))
!
g = 9.81
- pi = 4.*atan(1.)
+ pi = 4 * atan(1.0)
!
call timer(t0)
!
@@ -37,83 +35,102 @@ subroutine compute_wave_field()
!
do k = 1, no_nodes
if (inner(k)) then
- ee(:,k) = waveps
+ ee(:, k) = waveps
endif
enddo
!
ee_ig = waveps
!
- restart=1 !TODO TL: CHECK > we need this turned on right now for IG...
+ restart = 1 !TODO TL: CHECK > we need this turned on right now for IG...
!
endif
!
! Initialize wave period
!
do k = 1, no_nodes
+ !
if (inner(k)) then
+ !
Tp(k) = Tpini
+ !
endif
- if (neumannconnected(k)>0) then
- Tp(neumannconnected(k))=Tpini
+ !
+ if (neumannconnected(k) > 0) then
+ !
+ Tp(neumannconnected(k)) = Tpini
+ !
endif
+ !
enddo
!
! Compute celerities and refraction speed
!
- Tp = max(tpmean_bwv,Tpini) ! to check voor windgroei
- sig = 2.0*pi/Tp
+ Tp = max(tpmean_bwv, Tpini) ! to check voor windgroei
+ sig = 2.0 * pi / Tp
Tp_ig = tpmean_bwv_ig! TL: now determined in snapwave_boundaries.f90 instead of Tinc2ig*Tp
- sigm_ig = 2.0*pi/Tp_ig !TODO - TL: Question do we want Tp_ig now as contant, or also spatially varying like Tp ?
- !
- expon = -(sig*sqrt(depth/g))**(2.5)
- kwav = sig**2/g*(1.0-exp(expon))**(-0.4)
- C = sig/kwav
- nwav = 0.5+kwav*depth/sinh(min(2*kwav*depth,50.0))
- Cg = nwav*C
+ sigm_ig = 2.0 * pi / Tp_ig !TODO - TL: Question do we want Tp_ig now as contant, or also spatially varying like Tp ?
+ expon = - (sig * sqrt(depth / g))**2.5
+ kwav = sig**2 / g * (1.0 - exp(expon))**-0.4
+ C = sig / kwav
+ nwav = 0.5 + kwav * depth / sinh(min(2 * kwav * depth, 50.0))
+ Cg = nwav * C
!
if (igwaves) then
cg_ig = Cg
- expon = -(sigm_ig*sqrt(depth/g))**(2.5)
- kwav_ig = sig**2/g*(1.0-exp(expon))**(-0.4)
+ expon = -(sigm_ig * sqrt(depth / g))**2.5
+ kwav_ig = sig**2 / g * (1.0 - exp(expon))**-0.4
else
cg_ig = 0.0
kwav_ig = 0.0
endif
!
do k = 1, no_nodes
- sinhkh(k) = sinh(min(kwav(k)*depth(k), 50.0))
- Hmx(k) = gamma*depth(k)
+ !
+ sinhkh(k) = sinh(min(kwav(k) * depth(k), 50.0))
+ Hmx(k) = gamma * depth(k)
+ !
enddo
+ !
if (igwaves) then
+ !
do k = 1, no_nodes
- Hmx_ig(k) = 0.88/kwav_ig(k)*tanh(gamma_ig*kwav_ig(k)*depth(k)/0.88) ! Note - uses gamma_ig
+ Hmx_ig(k) = 0.88 / kwav_ig(k) * tanh(gamma_ig * kwav_ig(k) * depth(k) / 0.88) ! Note - uses gamma_ig
enddo
+ !
else
- Hmx_ig = 0.0
+ !
+ Hmx_ig = 0.0
+ !
endif
!
do itheta = 1, ntheta
- ctheta(itheta,:) = sig/sinh(min(2.0*kwav*depth, 50.0))*(dhdx*sin(theta(itheta)) - dhdy*cos(theta(itheta)))
+ !
+ ctheta(itheta,:) = sig / sinh(min(2 * kwav * depth, 50.0)) * (dhdx * sin(theta(itheta)) - dhdy * cos(theta(itheta)))
+ !
enddo
!
if (igwaves) then
+ !
do itheta = 1, ntheta
- ctheta_ig(itheta,:) = sigm_ig/sinh(min(2.0*kwav_ig*depth, 50.0))*(dhdx*sin(theta(itheta)) - dhdy*cos(theta(itheta)))
+ ctheta_ig(itheta,:) = sigm_ig / sinh(min(2 * kwav_ig * depth, 50.0)) * (dhdx * sin(theta(itheta)) - dhdy * cos(theta(itheta)))
enddo
+ !
else
+ !
ctheta_ig = 0.0
+ !
endif
!
! Limit unrealistic refraction speed to 1/2 pi per wave period
!
do k = 1, no_nodes
- ctheta(:,k) = sign(1.0, ctheta(:,k))*min(abs(ctheta(:, k)), sig(k)/4)
+ ctheta(:,k) = sign(1.0, ctheta(:,k)) * min(abs(ctheta(:, k)), sig(k) / 4)
enddo
!
if (igwaves) then
- do k=1, no_nodes
- ctheta_ig(:,k) = sign(1.0, ctheta_ig(:,k))*min(abs(ctheta_ig(:, k)), sigm_ig(k)/4.0)
- enddo
+ do k=1, no_nodes
+ ctheta_ig(:,k) = sign(1.0, ctheta_ig(:,k)) * min(abs(ctheta_ig(:, k)), sigm_ig(k) / 4)
+ enddo
endif
!
! Solve the directional wave energy balance on an unstructured grid
@@ -131,8 +148,8 @@ subroutine compute_wave_field()
Hmx, ee, windspreadfac, u10, niter, crit, &
hmin, baldock_ratio, baldock_ratio_ig, &
aa, sig, jadcgdx, sigmin, sigmax,&
- c_dispT, WsorE, WsorA, SwE, SwA, Tpini, &
- igwaves,kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, &
+ c_dispT, DoverE, WsorE, WsorA, SwE, SwA, Tpini, &
+ igwaves, kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, &
beta, srcig, alphaig, Dw_ig, Df_ig, &
vegetation, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, Dveg, &
zb, nwav, ig_opt, alpha_ig, gamma_ig, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig)
@@ -145,6 +162,7 @@ subroutine compute_wave_field()
end subroutine
+
subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
w, ds, prev, &
neumannconnected, &
@@ -156,14 +174,13 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
Hmx, ee, windspreadfac, u10, niter, crit, &
hmin, baldock_ratio, baldock_ratio_ig, &
aa, sig, jadcgdx, sigmin, sigmax,&
- c_dispT, WsorE, WsorA, SwE, SwA, Tpini, &
+ c_dispT, DoverE, WsorE, WsorA, SwE, SwA, Tpini, &
igwaves,kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, &
betamean, srcig, alphaig, Dw_ig, Df_ig, &
vegetation, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, Dveg, &
zb, nwav, ig_opt, alfa_ig, gamma_ig, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig)
!
use snapwave_windsource
- !use snapwave_ncoutput ! TL: removed, we don't use this in SF+SW
!
implicit none
!
@@ -203,7 +220,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
real*4, intent(in) :: alfa,gamma, gammax ! coefficients in Baldock wave breaking dissipation
real*4, intent(in) :: baldock_ratio ! option controlling from what depth wave breaking should take place: (Hk>baldock_ratio*Hmx(k)), default baldock_ratio=0.2
real*4, intent(in) :: baldock_ratio_ig ! option controlling from what depth wave breaking should take place for IG waves: (Hk_ig>baldock_ratio_ig*Hmx_ig(k)), default baldock_ratio_ig=0.2
- real*4, dimension(no_nodes), intent(inout) :: H ! wave height - TODO - TL - CHECK > inout needed to have updated 'H' for determining srcig
+ real*4, dimension(no_nodes), intent(out) :: H ! wave height
real*4, dimension(no_nodes), intent(out) :: H_ig ! wave height
real*4, dimension(no_nodes), intent(out) :: Dw ! wave breaking dissipation
real*4, dimension(no_nodes), intent(out) :: Dw_ig ! wave breaking dissipation IG
@@ -227,6 +244,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
real*4, intent(in) :: sigmin, sigmax, c_dispT
real*4, dimension(ntheta, no_nodes), intent(in) :: windspreadfac !< [-] distribution array for wind input
real*4, dimension(ntheta,no_nodes), intent(inout) :: aa
+ real*4, dimension(no_nodes), intent(inout) :: DoverE
real*4, dimension(ntheta,no_nodes), intent(out) :: WsorE, WsorA
real*4, dimension(no_nodes), intent(out) :: SwE, SwA
real*4, dimension(no_nodes), intent(inout) :: sig
@@ -264,7 +282,6 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
real*4, dimension(:), allocatable :: A,B,C,R ! coefficients in the tridiagonal matrix solved per point
real*4, dimension(:), allocatable :: B_aa,R_aa,aaprev ! coefficients in the tridiagonal matrix solved per point
real*4, dimension(:), allocatable :: A_ig,B_ig,C_ig,R_ig ! coefficients in the tridiagonal matrix solved per point
- real*4, dimension(:), allocatable :: DoverE ! ratio of mean wave dissipation over mean wave energy
real*4, dimension(:), allocatable :: DoverA ! ratio of mean wave dissipation over mean wave energy
real*4, dimension(:), allocatable :: DoverE_ig ! ratio of mean wave dissipation over mean wave energy
real*4, dimension(:), allocatable :: E ! mean wave energy
@@ -277,7 +294,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
real*4 :: pi = 4.*atan(1.0)
real*4 :: g=9.81
real*4 :: hmin ! minimum water depth! TL: make user changeable also here according to 'snapwave_hmin' in sfincs.inp
- real*4 :: fac=1.0 ! underrelaxation factor for DoverA
+ real*4 :: fac=0.25 ! underrelaxation factor for DoverA
real*4 :: oneoverdt
real*4 :: oneover2dtheta
real*4 :: rhog8
@@ -296,22 +313,32 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
real*4 :: Tinc2ig ! ratio compared to period Tinc to estimate Tig
real*4 :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0
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)
+ real*4 :: fdrspr ! Inc-IG reduction factor to account for directional spreading
integer, save :: callno=1
!
- real*4, dimension(ntheta) :: sinth,costh ! distribution of wave angles and offshore wave energy density
+ integer :: baldock_option ! 1 or 2
+ real*4 :: baldock_hrms2hs
+ !
+ real*4, dimension(ntheta) :: sinth, costh ! distribution of wave angles and offshore wave energy density
!
- !local wind source vars
+ ! Local wind source vars
!
real*4 :: Ak
real*4 :: DwT
real*4 :: DwAk
- real*4 :: ndissip !
- real*4 :: depthlimfac=1.0
+ real*4 :: ndissip
+ real*4 :: depthlimfac
real*4 :: waveps=0.0001
!
! Allocate local arrays
!
- waveps = 0.0001
+ waveps = 0.0001
+ baldock_option = 1 ! 1 or 2, 2 avoids insufficient dissipation at steep coasts
+ !baldock_hrms2hs = sqrt(2.0) ! or use 1.0 for original implementation
+ baldock_hrms2hs = 1.0
+ !
+ fdrspr = 0.65
+ !
allocate(ok(no_nodes)); ok=0
allocate(indx(no_nodes,4)); indx=0
allocate(eeold(ntheta,no_nodes)); eeold=0.0
@@ -322,7 +349,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
allocate(B(ntheta)); B=0.0
allocate(C(ntheta)); C=0.0
allocate(R(ntheta)); R=0.0
- allocate(DoverE(no_nodes)); DoverE=0.0
+ !allocate(DoverE(no_nodes)); DoverE=0.0
allocate(E(no_nodes)); E=waveps
allocate(Eold(no_nodes)); Eold=0.0
!
@@ -335,7 +362,6 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
allocate(cgprev_ig(ntheta)); cgprev_ig=0.0
allocate(DoverE_ig(no_nodes)); DoverE_ig=0.0
allocate(E_ig(no_nodes)); E_ig=waveps
- !allocate(T_ig(no_nodes)); T_ig=0.0
allocate(sigm_ig(no_nodes)); sigm_ig=0.0
allocate(depthprev(ntheta,no_nodes)); depthprev=0.0
allocate(beta_local(ntheta,no_nodes)); beta_local=0.0
@@ -358,55 +384,61 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
costh(itheta) = cos(theta(itheta))
enddo
!
- df = 0.0
- dw = 0.0
+ Df = 0.0
+ Dw = 0.0
+ Dveg = 0.0
F = 0.0
!
- ok = 0
- indx = 0
- eemax = maxval(ee)
- dtheta = theta(2) - theta(1)
- if (dtheta<0.) dtheta = dtheta + 2.*pi
+ ok = 0
+ indx = 0
+ eemax = maxval(ee)
+ dtheta = theta(2) - theta(1)
+ !
+ if (dtheta < 0.0) dtheta = dtheta + 2*pi
+ !
if (wind) then
- sig = 2*pi/Tpini
+ sig = 2 * pi / Tpini
else
- sig = 2*pi/Tp
+ sig = 2 * pi / Tp
endif
- oneoverdt = 1.0/dt
- oneover2dtheta = 1.0/2.0/dtheta
- rhog8 = 0.125*rho*g
+ !
+ oneoverdt = 1.0 / dt
+ oneover2dtheta = 1.0 / 2.0 / dtheta
+ rhog8 = 0.125 * rho * g
thetam = 0.0
- !H = 0.0 ! TODO - TL: CHeck > needed for restart for IG > set to 0 now in snapwave_domain.f90
- Dveg = 0.0
!
if (igwaves) then
- !T_ig = Tinc2ig*Tp
- sigm_ig = 2*pi/T_ig
+ !
+ sigm_ig = 2 * pi / T_ig
DoverE_ig = 0.0
+ !
endif
!
if (wind) then
+ !
DoverA = 0.0
ndissip = 3.0
WsorE = 0.0
WsorA = 0.0
- Ak = waveps/sigmax
+ Ak = waveps / sigmax
+ !
endif
!
- ! Sort coordinates in sweep directions
+ ! Sort coordinates in sweep directions (can we not do this already in snapwave_domain?)
+ !
+ shift = [0, 1, -1, 2]
!
- shift = [0,1,-1,2]
do sweep = 1, 4
!
- ra = x*cos(thetamean + 0.5*pi*shift(sweep)) + y*sin(thetamean + 0.5*pi*shift(sweep))
+ ra = x * cos(thetamean + 0.5 * pi * shift(sweep)) + y * sin(thetamean + 0.5 * pi * shift(sweep))
call hpsort_eps_epw(no_nodes, ra , indx(:, sweep), 1.0e-6)
!
enddo
- !
- ! Set inner to false for all points at grid edge or adjacent to dry point
!
- do k=1,no_nodes
- !
+ ! Set inner to false for all points at grid edge
+ !
+ do k = 1, no_nodes
+ !
do itheta = 1, ntheta
!
k1 = prev(1, itheta, k)
@@ -427,136 +459,117 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!exit
!
endif
+ !
enddo
enddo
!
+ ! Start iteration
!
- ! 0-a) Set boundary and initial conditions
- !
- do k = 1, no_nodes
+ do iter = 1, niter * 4
!
- ! Boundary condition at sea side (uniform)
+ ! write(*,*)'iter=',iter
!
- if (.not.inner(k)) then
+ sweep = mod(iter, 4) !TODO - TL: problem that we don't have option for sweep = 1 anymore?
+ !
+ if (sweep == 0) then
+ sweep = 4
+ endif
+ !
+ ! At start of each sweep, compute E, H, E_ig and H_ig, and aa
+ !
+ ! This loop can be parallelized !
+ !
+ do k = 1, no_nodes
!
- ee(:,k)=max(ee(:,k),waveps)
- E(k) = sum(ee(:, k))*dtheta
- H(k) = sqrt(8*E(k)/rho/g)
- thetam(k) = atan2(sum(ee(:, k)*sin(theta)), sum(ee(:, k)*cos(theta)))
+ if (ok(k) == 1) cycle
!
- !ee_ig(:, k) = eeinc2ig*ee(:,k) !TODO TL: determined in snapwave_boundaries.f90
+ ee(:,k) = max(ee(:, k), waveps)
!
- if (igwaves) then
- E_ig(k) = sum(ee_ig(:, k))*dtheta
- H_ig(k) = sqrt(8*E_ig(k)/rho/g)
- endif
- !
- if (wind) then
- sig(k) = 2*pi/Tp(k)
- !aa(:,k) = max(aa(:,k),waveps/sig(k))
- aa(:,k) = max(ee(:,k),waveps)/sig(k)
- Ak = E(k)/sig(k)
- endif
+ ! Limit energy with gammax
!
- endif
- enddo
- !
- ! 0-b) Determine IG source/sink term
- !
- if (igwaves) then
- !
- ! As defined in Leijnse, van Ormondt, van Dongeren, Aerts & Muis et al. 2024
- !
- ! Actual determining of source term:
- !
- call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
- !
- ! inout: alphaig_local, srcig_local - eeprev, eeprev_ig, cgprev, beta_local
- ! in: the rest
- !
- ! NOTE - This is based on the energy in the precious SnapWave timestep 'ee' and 'ee_ig', and waveheight 'H', which should therefore be made available.
- !
- endif
- !
- ! 0-c) Set initial condition at inner cells
- !
- do k = 1, no_nodes
- !
- if (inner(k)) then
+ depthlimfac = max(1.0, (sqrt(sum(ee(:, k)) * dtheta / rhog8) / (gammax * depth(k)) )**2.0)
+ ee(:,k) = ee(:,k) / depthlimfac
!
+ E(k) = sum(ee(:, k)) * dtheta
+ H(k) = sqrt(8 * E(k) / rho / g)
+ !
+ if (igwaves) then
+ !
+ depthlimfac = max(1.0, (sqrt(sum(ee_ig(:, k)) * dtheta / rhog8) / (gammax * depth(k)) )**2.0)
+ ee_ig(:,k) = ee_ig(:,k) / depthlimfac
+ !
+ E_ig(k) = sum(ee_ig(:, k)) * dtheta
+ H_ig(k) = sqrt(8 * E_ig(k) / rho / g)
+ !
+ endif
+ !
if (wind) then
- ee(:,k) = waveps
- sig(k) = 2*pi/Tpini
- aa(:,k) = ee(:,k)/sig(k)
- else
- ee(:,k) = waveps
+ !
+ sig(k) = 2 * pi / Tp(k)
+ sig(k) = max(sig(k), sigmin)
+ sig(k) = min(sig(k), sigmax)
+ aa(:, k) = max(ee(:, k), waveps) / sig(k)
+ !
endif
!
- ! Make sure DoverE is filled based on previous ee
- Ek = sum(ee(:, k))*dtheta
- Hk = min(sqrt(Ek/rhog8), gamma*depth(k))
- Ek = rhog8*Hk**2
- if (.not. wind) then
- uorbi = 0.5*sig(k)*Hk/sinhkh(k)
- Dfk = 0.28*rho*fw(k)*uorbi**3
- call baldock(rho, g, alfa, gamma, depth(k), Hk, Tp(k), 1, Dwk, Hmx(k))
- DoverE(k) = (Dwk + Dfk)/max(Ek, 1.0e-6)
- endif
+ ! Set Neumann boundaries
!
- if (wind) then
- call compute_celerities(depth(k), sig(k), sinth, costh, ntheta, gamma, dhdx(k), dhdy(k), sinhkh(k), Hmx(k), kwav(k), cg(k), ctheta(:,k))
- uorbi = 0.5*sig(k)*Hk/sinhkh(k)
- Dfk = 0.28*rho*fw(k)*uorbi**3
- call baldock(rho, g, alfa, gamma, depth(k), Hk, 2.0*pi/sig(k), 1, Dwk, Hmx(k))
- DoverE(k) = (Dwk + Dfk)/max(Ek, 1.0e-6)
+ if (neumannconnected(k) /= 0) then
+ !
+ ! Flip indices around compared to orginal implementation
+ !
+ ! Do we really need all of these?
+ !
+ kn = neumannconnected(k) ! Index of internal point
+ !
+ sinhkh(k) = sinhkh(kn)
+ kwav(k) = kwav(kn)
+ Hmx(k) = Hmx(kn)
+ ee(:, k) = ee(:, kn)
+ ee_ig(:, k) = ee_ig(:, kn)
+ ctheta(:, k) = ctheta(:, kn)
+ cg(k) = cg(kn)
+ !
+ if (wind) then
+ !
+ sig(k) = sig(kn)
+ Tp(k) = 2 * pi / sig(kn)
+ WsorE(:,k) = WsorE(:,kn)
+ WsorA(:,k) = WsorA(:,kn)
+ aa(:,k) = aa(:,kn)
+ !
+ endif
+ !
+ Df(kn) = Df(kn)
+ Dw(kn) = Dw(kn)
!
- ! initial conditions are not equal to bc conditions
- DwT = - c_dispT/(1.0 -ndissip)*(2.*pi)/sig(k)**2*cg(k)*kwav(k) * DoverE(k)
- DwAk = 0.5/pi * (E(k)*DwT+2.0*pi*Ak*DoverE(k) )
- DoverA(k) = DwAk/max(Ak,1e-6)
endif
!
- endif
- !
- enddo
- !
- ! Start iteration
- !
- do iter=1,niter
+ enddo
!
- sweep = mod(iter, 4) !TODO - TL: problem that we don't have option for sweep = 1 anymore?
- !
- if (sweep==0) then
- sweep = 4
- endif
- !
- if (sweep==1) then
+ if (sweep == 1) then
+ !
eeold = ee
+ !
do k = 1, no_nodes
Eold(k) = sum(eeold(:, k))
enddo
!
- if (igwaves) then
- !
- if (iterative_srcig) then
- ! Update H(k) based on updated ee(:,k), as used in IG source term to determine alphaig
- !
- do k = 1, no_nodes
- !
- if (inner(k)) then
- !
- H(k) = sqrt(8*sum(ee(:, k))*dtheta/rho/g)
- !
- endif
- enddo
- !
- ! Actual determining of source term - every first sweep of iteration
- !
- call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
- !
- endif
+ endif
+ !
+ ! Update IG source and sink terms
+ !
+ if (igwaves) then
+ !
+ ! Do this in each first sweep, or (in case of iterative_srcig) in every sweep
+ !
+ if (sweep == 1 .or. iterative_srcig) then
!
- endif
+ call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, &
+ ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, &
+ alphaig_local, beta_local, srcig_local)
+ !
+ endif
!
endif
!
@@ -564,233 +577,269 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
do count = 1, no_nodes
!
- k=indx(count, sweep)
+ k = indx(count, sweep)
!
- if (inner(k)) then
- if (depth(k)>1.1*hmin) then
+ if (inner(k)) then ! Regular point
+ !
+ ! Set Ek, Hk, Ek_ig, Hk_ig and Ak
+ !
+ Ek = E(k)
+ Hk = H(k)
+ !
+ if (igwaves) then
+ !
+ Ek_ig = E_ig(k)
+ Hk_ig = H_ig(k)
+ !
+ endif
+ !
+ if (wind) then
+ !
+ Ak = Ek / sig(k)
!
- if (ok(k) == 0) then
+ endif
+ !
+ if (depth(k) > hmin) then
+ !
+ if (ok(k) == 0) then ! Only perform computations on wet inner points that are not yet converged (ok=1)
!
- ! Only perform computations on wet inner points that are not yet converged (ok)
+ ! Get upwind data
!
do itheta = 1, ntheta
!
k1 = prev(1, itheta, k)
k2 = prev(2, itheta, k)
!
- eeprev(itheta) = w(1, itheta, k)*ee(itheta, k1) + w(2, itheta, k)*ee(itheta, k2)
- cgprev(itheta) = w(1, itheta, k)*cg(k1) + w(2, itheta, k)*cg(k2)
+ eeprev(itheta) = w(1, itheta, k) * ee(itheta, k1) + w(2, itheta, k) * ee(itheta, k2)
+ cgprev(itheta) = w(1, itheta, k) * cg(k1) + w(2, itheta, k) * cg(k2)
!
if (igwaves) then
- eeprev_ig(itheta) = w(1, itheta, k)*ee_ig(itheta, k1) + w(2, itheta, k)*ee_ig(itheta, k2)
- cgprev_ig(itheta) = w(1, itheta, k)*cg_ig(k1) + w(2, itheta, k)*cg_ig(k2)
+ !
+ eeprev_ig(itheta) = w(1, itheta, k) * ee_ig(itheta, k1) + w(2, itheta, k) * ee_ig(itheta, k2)
+ cgprev_ig(itheta) = w(1, itheta, k) * cg_ig(k1) + w(2, itheta, k) * cg_ig(k2)
+ !
endif
!
if (wind) then
- aaprev(itheta) = w(1, itheta, k)*aa(itheta, k1) + w(2, itheta, k)*aa(itheta, k2)
+ !
+ aaprev(itheta) = w(1, itheta, k) * aa(itheta, k1) + w(2, itheta, k) * aa(itheta, k2)
+ aaprev(itheta) = min(aaprev(itheta), eeprev(itheta) / sigmin)
+ aaprev(itheta) = max(aaprev(itheta), eeprev(itheta) / sigmax)
+ !
endif
!
enddo
!
- Ek = sum(eeprev)*dtheta ! to check
+ ! Fill DoverE
!
- depthlimfac = max(1.0, (sqrt(Ek/rhog8)/(gammax*depth(k)))**2.0)
- Hk = min(sqrt(Ek/rhog8), gamma*depth(k))
- Ek = Ek/depthlimfac
+ ! Bottom friction
!
- if (wind) then
- !
- Ak = sum(aaprev)*dtheta
- !
- Ak = Ak/depthlimfac
- ee(:,k) = ee(:,k) / depthlimfac
- aa(:,k) = aa(:,k) / depthlimfac
- sig(k) = Ek/Ak
- sig(k) = max(sig(k),sigmin)
- sig(k) = min(sig(k),sigmax)
- Ak = Ek/sig(k) ! to avoid small T in windinput
- if (wind) then
- aaprev=min(aaprev,eeprev/sigmin)
- aaprev=max(aaprev,eeprev/sigmax)
- endif
- !
- call compute_celerities(depth(k), sig(k), sinth, costh, ntheta, gamma, dhdx(k), dhdy(k), sinhkh(k), Hmx(k), kwav(k), cg(k), ctheta(:,k))
- endif
- !
- ! Fill DoverE
- uorbi = 0.5*sig(k)*Hk/sinhkh(k)
- Dfk = 0.28*rho*fw(k)*uorbi**3
- !if (Hk>0.) then !
- if (Hk>baldock_ratio*Hmx(k)) then
- call baldock(rho, g, alfa, gamma, depth(k), Hk, 2*pi/sig(k) , 1, Dwk, Hmx(k))
- else
- Dwk = 0.
- endif
- !
- if (vegetation) then
- call vegatt(sig(k), no_nodes, kwav(k), no_secveg, veg_ah(k,:), veg_bstems(k,:), veg_Nstems(k,:), veg_Cd(k,:), depth(k), rho, g, Hk, Dvegk)
+ uorbi = 0.5 * sig(k) * Hk / sinhkh(k)
+ Dfk = 0.28 * rho * fw(k) * uorbi**3
+ !
+ ! Wave breaking
+ !
+ if (Hk > baldock_ratio * Hmx(k)) then
+ !
+ ! Baldock may expect Hs so multiply H and Hmax with baldock_hrms2hs (sqrt(2))
+ !
+ call baldock(rho, g, alfa, gamma, depth(k), Hk * baldock_hrms2hs, 2 * pi / sig(k) , baldock_option, Dwk, Hmx(k) * baldock_hrms2hs)
+ !
else
- Dvegk = 0.
+ !
+ Dwk = 0.0
+ !
endif
!
- DoverE(k) = (Dwk + Dfk + Dvegk)/max(Ek, 1.0e-6)
+ ! Vegetation
!
- if (wind) then
+ if (vegetation) then
!
- if (iter==1) then
- call windinput(u10(k), rho, g, depth(k), ntheta, windspreadfac(:,k), Ek, Ak, cg(k), eeprev, aaprev, ds(:,k), WsorE(:,k), WsorA(:,k), jadcgdx)
- else
- call windinput(u10(k), rho, g, depth(k), ntheta, windspreadfac(:,k), Ek, Ak, cg(k), ee(:,k), aa(:,k), ds(:,k), WsorE(:,k), WsorA(:,k), jadcgdx)
- endif
+ call vegatt(sig(k), no_nodes, kwav(k), no_secveg, veg_ah(k,:), veg_bstems(k,:), veg_Nstems(k,:), veg_Cd(k,:), depth(k), rho, g, Hk, Dvegk)
!
- DwT = - c_dispT/(1.0 -ndissip)*(2.0*pi)/sig(k)**2*cg(k)*kwav(k) * DoverE(k)
- DwAk = 1/2.0/pi * (E(k)*DwT+2.0*pi*Ak*DoverE(k) )
+ else
!
- if (iter==1) then
- DoverA(k) = DwAk/max(Ak,1e-6)
- else
- DoverA(k) = (1.0-fac)*DoverA(k)+fac*DwAk/max(Ak,1.e-6)
- endif
- !
- call compute_celerities(depth(k), sig(k), sinth, costh, ntheta, gamma, dhdx(k), dhdy(k), sinhkh(k), Hmx(k), kwav(k), cg(k), ctheta(:,k))
+ Dvegk = 0.0
!
endif
!
- do itheta = 1, ntheta
+ !DoverE(k) = (Dwk + Dfk + Dvegk) / max(Ek, 1.0e-6)
+ !
+ ! Re-introduce relaxation
+ !
+ DoverE(k) = (1.0 - fac) * DoverE(k) + fac * (Dwk + Dfk + Dvegk) / max(Ek, 1.0e-6)
+ !
+ ! Store dissipation terms for output
+ !
+ Df(k) = Dfk
+ Dw(k) = Dwk
+ Dveg(k) = Dvegk
+ !
+ do itheta = 1, ntheta
!
- R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta)*eeprev(itheta)/ds(itheta, k) - srcig_local(itheta, k) * shinc2ig
+ R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta) * eeprev(itheta) / ds(itheta, k) - srcig_local(itheta, k) * shinc2ig * fdrspr
!
enddo
!
do itheta = 2, ntheta - 1
!
- A(itheta) = -ctheta(itheta - 1, k)*oneover2dtheta
- B(itheta) = oneoverdt + cg(k)/ds(itheta,k) + DoverE(k)
- C(itheta) = ctheta(itheta + 1, k)*oneover2dtheta
+ A(itheta) = -ctheta(itheta - 1, k) * oneover2dtheta
+ B(itheta) = oneoverdt + cg(k) / ds(itheta,k) + DoverE(k)
+ C(itheta) = ctheta(itheta + 1, k) * oneover2dtheta
!
enddo
!
- A(1) = -ctheta(ntheta, k)*oneover2dtheta
- B(1) = oneoverdt + cg(k)/ds(1,k) + DoverE(k)
- C(1) = ctheta(2, k)*oneover2dtheta
+ A(1) = -ctheta(ntheta, k) * oneover2dtheta
+ B(1) = oneoverdt + cg(k) / ds(1,k) + DoverE(k)
+ C(1) = ctheta(2, k) * oneover2dtheta
!
- A(ntheta) = -ctheta(ntheta - 1, k)*oneover2dtheta
- B(ntheta) = oneoverdt + cg(k)/ds(ntheta,k) + DoverE(k)
- C(ntheta) = ctheta(1, k)*oneover2dtheta
+ A(ntheta) = -ctheta(ntheta - 1, k) * oneover2dtheta
+ B(ntheta) = oneoverdt + cg(k) / ds(ntheta,k) + DoverE(k)
+ C(ntheta) = ctheta(1, k) * oneover2dtheta
!
! Solve tridiagonal system per point
!
if (wind) then
+ !
+ call compute_celerities(depth(k), sig(k), sinth, costh, ntheta, gamma, dhdx(k), dhdy(k), sinhkh(k), Hmx(k), kwav(k), cg(k), ctheta(:,k))
+ !
+ call windinput(u10(k), rho, g, depth(k), ntheta, windspreadfac(:,k), Ek, Ak, cg(k), ee(:,k), aa(:,k), ds(:,k), WsorE(:,k), WsorA(:,k), jadcgdx)
+ !
+ DwT = - c_dispT / (1.0 - ndissip) * (2 * pi) / sig(k)**2 * cg(k) * kwav(k) * DoverE(k)
+ DwAk = 1.0 / 2.0 / pi * (E(k) * DwT + 2.0 * pi * Ak * DoverE(k) )
+ !
+ if (iter == 1) then
+ !
+ DoverA(k) = DwAk / max(Ak, 1e-6)
+ !
+ else
+ !
+ DoverA(k) = (1.0 - fac) * DoverA(k) + fac * DwAk / max(Ak, 1.0e-6)
+ !
+ endif
+ !
do itheta = 2, ntheta - 1
- B_aa(itheta) = oneoverdt + cg(k)/ds(itheta,k) + DoverA(k)
- R_aa(itheta) = (oneoverdt)*aa(itheta, k) + cgprev(itheta)*aaprev(itheta)/ds(itheta, k)
+ !
+ B_aa(itheta) = oneoverdt + cg(k) / ds(itheta,k) + DoverA(k)
+ R_aa(itheta) = (oneoverdt) * aa(itheta, k) + cgprev(itheta) * aaprev(itheta) / ds(itheta, k)
+ !
enddo
!
- if (ctheta(1,k)<0) then
- B_aa(1) = oneoverdt - ctheta(1, k)/dtheta + cg(k)/ds(1, k) + DoverA(k)
- R_aa(1) = (oneoverdt)*aa(1, k) + cgprev(1)*aaprev(1)/ds(1, k)
+ if (ctheta(1, k) < 0) then
+ !
+ B_aa(1) = oneoverdt - ctheta(1, k) / dtheta + cg(k) / ds(1, k) + DoverA(k)
+ R_aa(1) = (oneoverdt) * aa(1, k) + cgprev(1) * aaprev(1) / ds(1, k)
+ !
else
- B_aa(1) = oneoverdt + cg(k)/ds(1, k) + DoverA(k)
- R_aa(1) = (oneoverdt)*aa(1, k) + cgprev(1)*aaprev(1)/ds(1, k)
+ !
+ B_aa(1) = oneoverdt + cg(k) / ds(1, k) + DoverA(k)
+ R_aa(1) = (oneoverdt) * aa(1, k) + cgprev(1) * aaprev(1) / ds(1, k)
+ !
endif
!
- if (ctheta(ntheta, k)>0) then
- B_aa(ntheta) = oneoverdt + ctheta(ntheta, k)/dtheta + cg(k)/ds(ntheta, k) + DoverA(k)
- R_aa(ntheta) = (oneoverdt )*aa(ntheta,k) + cgprev(ntheta)*aaprev(ntheta)/ds(ntheta, k)
+ if (ctheta(ntheta, k) > 0) then
+ !
+ B_aa(ntheta) = oneoverdt + ctheta(ntheta, k) / dtheta + cg(k) / ds(ntheta, k) + DoverA(k)
+ R_aa(ntheta) = oneoverdt * aa(ntheta,k) + cgprev(ntheta) * aaprev(ntheta) / ds(ntheta, k)
+ !
else
- B_aa(ntheta) = oneoverdt + cg(k)/ds(ntheta, k) + DoverA(k)
- R_aa(ntheta) = (oneoverdt)*aa(ntheta,k) + cgprev(ntheta)*aaprev(ntheta)/ds(ntheta, k)
+ !
+ B_aa(ntheta) = oneoverdt + cg(k) / ds(ntheta, k) + DoverA(k)
+ R_aa(ntheta) = (oneoverdt) * aa(ntheta,k) + cgprev(ntheta) * aaprev(ntheta) / ds(ntheta, k)
+ !
endif
+ !
R(:) = R(:) + WsorE(:,k)
R_aa(:) = R_aa(:) + WsorA(:,k)
!
call solve_tridiag(A, B, C, R, ee(:,k), ntheta)
- call solve_tridiag(A,B_aa,C,R_aa,aa(:,k),ntheta)
- ee(:, k) = max(ee(:, k), waveps)
- aa(:,k) = max(aa(:,k),waveps/sigmax)
- aa(:,k) = max(aa(:,k),waveps/sig(k))
- !
- Ek = sum(ee(:, k))*dtheta
- Ak = sum(aa(:,k))*dtheta
+ call solve_tridiag(A, B_aa, C, R_aa, aa(:,k), ntheta)
!
- depthlimfac = max(1.0, (sqrt(Ek/rhog8)/(gammax*depth(k)))**2.0)
- Hk = sqrt(Ek/rhog8/depthlimfac)
- Ek = Ek/depthlimfac
- Ak = Ak/depthlimfac
- ee(:,k) = ee(:,k)/depthlimfac
- aa(:,k) = aa(:,k)/depthlimfac
+ ee(:, k) = max(ee(:, k), waveps)
+ aa(:,k) = max(aa(:,k), waveps / sigmax)
+ aa(:,k) = max(aa(:,k), waveps / sig(k))
!
- sig(k) = Ek/Ak
- sig(k) = max(sig(k),sigmin)
- sig(k) = min(sig(k),sigmax)
- call compute_celerities(depth(k), sig(k), sinth, costh, ntheta, gamma, dhdx(k), dhdy(k), sinhkh(k), Hmx(k), kwav(k), cg(k), ctheta(:,k))
- if (sig(k)<0.1) then
- a=1
- endif
else
!
! Solve tridiagonal system per point
!
- call solve_tridiag(A, B, C, R, ee(:,k), ntheta)
- ee(:, k) = max(ee(:, k),waveps)
+ call solve_tridiag(A, B, C, R, ee(:, k), ntheta)
+ !
+ ee(:, k) = max(ee(:, k), waveps)
!
endif !wind
!
! IG
!
if (igwaves) then
- Ek_ig = sum(eeprev_ig)*dtheta
- !Hk_ig = sqrt(Ek_ig/rhog8) !org trunk
- Hk_ig = min(sqrt(Ek_ig/rhog8), gamma_ig*depth(k)) !TL: Question - why not this one?
- Ek_ig = rhog8*Hk_ig**2
+ !
+ ! Update Hk
+ !
+ H(k) = sqrt(8 * sum(ee(:, k)) * dtheta / rho / g)
!
! Bottom friction Henderson and Bowen (2002) - D = 0.015*rhow*(9.81/depth(k))**1.5*(Hk/sqrt(8.0))*Hk_ig**2/8
!
- Dfk_ig = fw_ig(k)*0.0361*(9.81/depth(k))**1.5*Hk*Ek_ig
+ Dfk_ig = fw_ig(k) * 0.0361 * (9.81 / depth(k))**1.5 * Hk * Ek_ig
!
! Dissipation of infragravity waves
!
- if (Hk_ig>baldock_ratio_ig*Hmx_ig(k)) then
- call baldock(rho, g, alfa_ig, gamma_ig, depth(k), Hk_ig, T_ig(k), 1, Dwk_ig, Hmx_ig(k))
+ if (Hk_ig > baldock_ratio_ig * Hmx_ig(k)) then
+ !
+ call baldock(rho, g, alfa_ig, gamma_ig, depth(k), Hk_ig * baldock_hrms2hs, T_ig(k), baldock_option, Dwk_ig, Hmx_ig(k) * baldock_hrms2hs)
+ !
else
- Dwk_ig = 0.
+ !
+ Dwk_ig = 0.0
+ !
endif
!
- DoverE_ig(k) = (Dwk_ig + Dfk_ig)/max(Ek_ig, 1.0e-6) ! org trunk
- !DoverE_ig(k) = (1.0 - fac)*DoverE_ig(k) + fac*(Dwk_ig + Dfk_ig)/max(Ek_ig, 1.0e-6) ! TODO - TL CHECK - why not with relaxation anymore?
+ ! Store dissipation terms for output
+ !
+ Df_ig(k) = Dfk_ig
+ Dw_ig(k) = Dwk_ig
+ !
+ DoverE_ig(k) = (Dwk_ig + Dfk_ig) / max(Ek_ig, 1.0e-6)
!
do itheta = 1, ntheta
!
- R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta)*eeprev_ig(itheta)/ds(itheta, k) + srcig_local(itheta, k) !TL: new version
+ R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta) * eeprev_ig(itheta) / ds(itheta, k) + srcig_local(itheta, k) * fdrspr
!
enddo
!
do itheta = 2, ntheta - 1
!
- A_ig(itheta) = -ctheta_ig(itheta - 1, k)*oneover2dtheta
- B_ig(itheta) = oneoverdt + cg_ig(k)/ds(itheta,k) + DoverE_ig(k)
- C_ig(itheta) = ctheta_ig(itheta + 1, k)*oneover2dtheta
+ A_ig(itheta) = - ctheta_ig(itheta - 1, k) * oneover2dtheta
+ B_ig(itheta) = oneoverdt + cg_ig(k) / ds(itheta,k) + DoverE_ig(k)
+ C_ig(itheta) = ctheta_ig(itheta + 1, k) * oneover2dtheta
!
enddo
!
- if (ctheta_ig(1,k)<0) then
+ if (ctheta_ig(1,k) < 0) then
+ !
A_ig(1) = 0.0
- B_ig(1) = oneoverdt - ctheta_ig(1, k)/dtheta + cg_ig(k)/ds(1, k) + DoverE_ig(k)
- C_ig(1) = ctheta_ig(2, k)/dtheta
+ B_ig(1) = oneoverdt - ctheta_ig(1, k) / dtheta + cg_ig(k) / ds(1, k) + DoverE_ig(k)
+ C_ig(1) = ctheta_ig(2, k) / dtheta
+ !
else
- A_ig(1)=0.0
- B_ig(1)=1.0/dt + cg_ig(k)/ds(1, k) + DoverE_ig(k)
- C_ig(1)=0.0
+ !
+ A_ig(1) = 0.0
+ B_ig(1) = 1.0 / dt + cg_ig(k) / ds(1, k) + DoverE_ig(k)
+ C_ig(1) = 0.0
+ !
endif
!
- if (ctheta_ig(ntheta, k)>0) then
- A_ig(ntheta) = -ctheta_ig(ntheta - 1, k)/dtheta
- B_ig(ntheta) = oneoverdt + ctheta_ig(ntheta, k)/dtheta + cg_ig(k)/ds(ntheta, k) + DoverE_ig(k)
+ if (ctheta_ig(ntheta, k) > 0) then
+ !
+ A_ig(ntheta) = -ctheta_ig(ntheta - 1, k) / dtheta
+ B_ig(ntheta) = oneoverdt + ctheta_ig(ntheta, k) / dtheta + cg_ig(k) / ds(ntheta, k) + DoverE_ig(k)
C_ig(ntheta) = 0.0
+ !
else
+ !
A_ig(ntheta) = 0.0
- B_ig(ntheta) = oneoverdt + cg_ig(k)/ds(ntheta, k) + DoverE_ig(k)
+ B_ig(ntheta) = oneoverdt + cg_ig(k) / ds(ntheta, k) + DoverE_ig(k)
C_ig(ntheta) = 0.0
+ !
endif
!
! Solve tridiagonal system per point
@@ -809,35 +858,17 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
else
!
ee(:, k) = 0.0
+ !
if (wind) then
aa(:,k) = 0.0
endif
+ !
ee_ig(:, k) = 0.0
!
endif
!
endif
!
- if (neumannconnected(k)/=0) then
- kn = neumannconnected(k)
- sinhkh(kn) = sinhkh(k)
- kwav(kn) = kwav(k)
- Hmx(kn) = Hmx(k)
- ee(:, kn) = ee(:, k)
- ee_ig(:, kn) = ee_ig(:, k) ! TL: Added Neumann option for IG
- ctheta(:, kn) = ctheta(:, k)
- cg(kn) = cg(k)
- if (wind) then
- sig(kn) = sig(k)
- Tp(kn) = 2.0*pi/sig(k)
- WsorE(:,kn) = WsorE(:,k)
- WsorA(:,kn) = WsorA(:,k)
- aa(:,kn) = aa(:,k)
- endif
- Df(kn) = Df(k)
- Dw(kn) = Dw(k)
- endif
- !
enddo
!
if (sweep==4) then
@@ -849,103 +880,87 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
dee = ee(:, k) - eeold(:, k)
diff(k) = maxval(abs(dee))
!
- if (diff(k)/eemax Ozeren et al. (2013); 2 => Mendez and Losada (2004)
!
myflag = 2
- pi = 4.d0*atan(1.d0)
+ pi = 4.d0*atan(1.d0)
!
! Representative wave period
+ !
Tp = 2*pi/sigm
!
! Coefficient alfa
+ !
if (ahh>=depth) then
alfav = 1.d0
else
@@ -1518,6 +1565,7 @@ subroutine bulkdragcoeff(ahh, m, Cdterm, no_nodes, no_secveg, depth, H, kwav, ve
KC = um*Tp/veg_bstems
!
! Bulk drag coefficient
+ !
if (myflag == 1) then
!
! Approach from Ozeren et al. (2013), eq?
@@ -1527,17 +1575,20 @@ subroutine bulkdragcoeff(ahh, m, Cdterm, no_nodes, no_secveg, depth, H, kwav, ve
else
Cdterm = 0.036d0+50.d0/(10.d0**0.926d0)
endif
+ !
elseif (myflag == 2) then
!
! Approach from Mendez and Losada (2004), eq. 40
! Only applicable for Laminaria Hyperborea (kelp)???
!
Q = KC/(alfav**0.76d0)
+ !
if (Q>=7) then
Cdterm = exp(-0.0138*Q)/(Q**0.3d0)
else
Cdterm = exp(-0.0138*7)/(7**0.3d0)
endif
+ !
endif
!
end subroutine bulkdragcoeff