diff --git a/.gitignore b/.gitignore index 22752d026..13f9c8823 100644 --- a/.gitignore +++ b/.gitignore @@ -65,3 +65,12 @@ source/third_party_open/netcdf/x64 source/sfincs/sfincs.opt.yaml /source/sfincs_lib/*.yaml /source/third_party_open/netcdf/netcdf-fortran-4.6.1/Debug +/docs/_build +/source/build_nvfortran_gpu_h7.sh +/source/build_nvfortran_gpu.sh +/source/Singularityfile-gpu.def +/source/Dockerfile.xpu +/source/Dockerfile.gpu.test +/source/Dockerfile.gpu.update01 +/source/Dockerfile.gpu +/source/Dockerfile.gpu.25.5.ccall diff --git a/source/code_docs/vegetation_lookup_figure.png b/source/code_docs/vegetation_lookup_figure.png new file mode 100644 index 000000000..d54d87bca Binary files /dev/null and b/source/code_docs/vegetation_lookup_figure.png differ diff --git a/source/code_docs/vegetation_lookup_figure.py b/source/code_docs/vegetation_lookup_figure.py new file mode 100644 index 000000000..55d32f888 --- /dev/null +++ b/source/code_docs/vegetation_lookup_figure.py @@ -0,0 +1,151 @@ +""" +Visualisation of SFINCS vegetation drag lookup table concept. +Run with: python vegetation_lookup_figure.py +""" + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.patches as mpatches +from matplotlib.gridspec import GridSpec + +# -------------------------------------------------------------------------- +# Example vegetation: 3 vertical sections stacked from bed +# -------------------------------------------------------------------------- +sections = [ + {"ah": 0.5, "cd_wd": 0.8, "label": "section 1\n(dense canopy base)", "color": "#4CAF50"}, + {"ah": 0.8, "cd_wd": 0.4, "label": "section 2\n(mid canopy)", "color": "#8BC34A"}, + {"ah": 0.4, "cd_wd": 0.15, "label": "section 3\n(sparse canopy top)", "color": "#CDDC39"}, +] + +# Cumulative section boundaries +bottoms = [0.0] +for s in sections: + bottoms.append(bottoms[-1] + s["ah"]) +tops = bottoms[1:] +h_veg_total = bottoms[-1] # = 1.7 m + +# -------------------------------------------------------------------------- +# Build lookup table (vegetation_nlookup = 20) +# -------------------------------------------------------------------------- +nlookup = 20 +dh = h_veg_total / nlookup +h_nodes = np.arange(0, nlookup + 1) * dh # shape (21,) + +table = np.zeros(nlookup + 1) +for k in range(1, nlookup + 1): + h_k = k * dh + sec_bot = 0.0 + for s in sections: + sec_top = sec_bot + s["ah"] + table[k] += s["cd_wd"] * max(0.0, min(sec_top, h_k) - sec_bot) + sec_bot = sec_top + +slope_table = np.diff(table) # length nlookup + +# -------------------------------------------------------------------------- +# Runtime lookup example: water depth hu +# -------------------------------------------------------------------------- +hu = 1.1 # example water depth at a timestep +hu_eff = min(hu, h_veg_total) +frac = hu_eff / dh +ik = min(int(frac), nlookup - 1) +frac_r = frac - ik +veg_cd_eff = table[ik] + frac_r * slope_table[ik] + +# -------------------------------------------------------------------------- +# Plot +# -------------------------------------------------------------------------- +fig = plt.figure(figsize=(13, 6)) +fig.suptitle("SFINCS vegetation drag — lookup table concept", fontsize=13, fontweight="bold") +gs = GridSpec(1, 3, figure=fig, wspace=0.40, left=0.06, right=0.97, top=0.88, bottom=0.10) + +# ── Panel A: vegetation sections ───────────────────────────────────────── +ax1 = fig.add_subplot(gs[0]) +ax1.set_title("A Vegetation sections\n(single uv-point)", fontsize=10) + +for i, s in enumerate(sections): + ax1.barh( + bottoms[i] + s["ah"] / 2, 1.0, + height=s["ah"], left=0, + color=s["color"], edgecolor="k", linewidth=0.8, alpha=0.85 + ) + ax1.text( + 0.5, bottoms[i] + s["ah"] / 2, + f'cd·b·N = {s["cd_wd"]:.2f}\nah = {s["ah"]} m', + ha="center", va="center", fontsize=7.5 + ) + # section boundary lines + ax1.axhline(bottoms[i], color="gray", lw=0.7, ls="--") + +ax1.axhline(h_veg_total, color="gray", lw=0.7, ls="--") +ax1.set_xlim(0, 1); ax1.set_xticks([]) +ax1.set_ylim(-0.15, 2.1) +ax1.set_ylabel("Height above bed (m)") +ax1.set_xlabel("← stem density schematic →") + +# water surface at hu +ax1.axhline(hu, color="#1565C0", lw=2, ls="-", label=f"water depth hu = {hu} m") +ax1.fill_betweenx([0, hu], 0, 1, color="#90CAF9", alpha=0.25) +ax1.legend(fontsize=7.5, loc="upper right") + +# ── Panel B: lookup table ───────────────────────────────────────────────── +ax2 = fig.add_subplot(gs[1]) +ax2.set_title("B Pre-computed lookup table\n(built once in initialize_vegetation)", fontsize=10) + +ax2.step(table, h_nodes, where="post", color="k", lw=1.2, label="table values") +ax2.plot(table, h_nodes, "o", color="k", ms=4, zorder=5) + +# shade each section contribution band +sec_bot = 0.0 +for s in sections: + sec_top = sec_bot + s["ah"] + ax2.axhspan(sec_bot, sec_top, color=s["color"], alpha=0.20) + sec_bot = sec_top + +# interpolation at hu +ax2.axhline(hu_eff, color="#1565C0", lw=1.5, ls="--", label=f"hu = {hu} m") +ax2.plot(veg_cd_eff, hu_eff, "*", color="red", ms=12, zorder=10, + label=f"veg_cd_eff = {veg_cd_eff:.3f} m²/s² (interpolated)") + +# show interpolation bracket +ax2.plot([table[ik], table[ik+1]], [h_nodes[ik], h_nodes[ik+1]], + "r-", lw=1.5, label="linear interpolation") +ax2.plot([table[ik], table[ik+1]], [h_nodes[ik], h_nodes[ik+1]], + "rs", ms=6) + +ax2.set_xlabel("Cumulative drag integral\n∑ cd·b·N · submerged thickness (m²/s²)") +ax2.set_ylabel("Depth level h_k (m)") +ax2.set_ylim(-0.15, 2.1) +ax2.legend(fontsize=7.5, loc="lower right") + +# ── Panel C: flux update (implicit) ────────────────────────────────────── +ax3 = fig.add_subplot(gs[2]) +ax3.set_title("C Explicit flux update\n(every timestep in compute_fluxes)", fontsize=10) +ax3.axis("off") + +textblock = ( + r"$\bf{Runtime\ lookup\ (O(1),\ no\ inner\ loop):}$" + "\n\n" + r"$\mathrm{frac} = \min(h_u,\; h_{veg}) \;/\; \Delta h$" + "\n" + r"$ik = \lfloor \mathrm{frac} \rfloor$" + "\n" + r"$cd_{eff} = \mathrm{table}[ik] + (\mathrm{frac}-ik)\times\mathrm{slope}[ik]$" + "\n\n" + r"$\bf{Explicit\ momentum\ update:}$" + "\n\n" + r"$F_{veg} = -\phi\;cd_{eff}\;u_0\;|u_0|$" + "\n\n" + r"$q^{n+1} = \dfrac{q^n + (F_{ext} + F_{veg})\,\Delta t}" + r"{1 + \dfrac{g\,n^2\,|q|}{h_u^{7/3}}\,\Delta t}$" + + "\n\n" + r"$\bf{Key\ properties:}$" + "\n" + "• No inner loop over sections at runtime\n" + "• Linear interpolation between table bins\n" + "• Sections stacked from bed upward\n" + " (consistent with SnapWave swvegatt)\n" + r"• $\phi$ = wet fraction (subgrid correction)" +) + +ax3.text(0.03, 0.97, textblock, transform=ax3.transAxes, + fontsize=9, va="top", ha="left", + bbox=dict(boxstyle="round,pad=0.5", fc="#F5F5F5", ec="#BDBDBD"), + linespacing=1.6) + +plt.savefig("vegetation_lookup_figure.png", dpi=150) +print("Saved: vegetation_lookup_figure.png") +plt.show() diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj index 34bf8f520..7ffd6464c 100644 --- a/source/sfincs_lib/sfincs_lib.vfproj +++ b/source/sfincs_lib/sfincs_lib.vfproj @@ -31,11 +31,15 @@ + - - - + + + + + + @@ -45,34 +49,48 @@ + + + + + + + + + - + + - + + - + + + - + + - - - + + + @@ -93,9 +111,9 @@ - + + - @@ -104,8 +122,11 @@ + + - + + @@ -116,6 +137,8 @@ + + diff --git a/source/src/Makefile.am b/source/src/Makefile.am index 891652f08..7c0a0698c 100644 --- a/source/src/Makefile.am +++ b/source/src/Makefile.am @@ -19,13 +19,16 @@ libsfincs_la_SOURCES = \ sfincs_date.f90 \ sfincs_spiderweb.f90 \ sfincs_data.f90 \ - ../third_party_open/Delft3D/astro.f90 \ + sfincs_read.f90 \ + ../third_party_open/Delft3D/astro.f90 \ ../third_party_open/utils/geometry.f90 \ sfincs_error.f90 \ sfincs_quadtree.f90 \ + sfincs_vegetation.f90 \ snapwave/interp.F90 \ snapwave/snapwave_data.f90 \ - snapwave/snapwave_ncinput.F90 \ + snapwave/snapwave_ncinput.F90 \ + snapwave/snapwave_ncoutput.F90 \ snapwave/snapwave_infragravity.f90 \ snapwave/snapwave_boundaries.f90 \ snapwave/snapwave_date.f90 \ @@ -50,7 +53,7 @@ libsfincs_la_SOURCES = \ sfincs_snapwave.f90 \ ../third_party_open/utils/deg2utm.f90 \ sfincs_meteo.f90 \ - ../third_party_open/bicgstab/bicgstab_solver_ilu.f90 \ + ../third_party_open/bicgstab/bicgstab_solver_ilu.f90 \ sfincs_nonhydrostatic.f90 \ sfincs_ncoutput.F90 \ sfincs_output.f90 \ diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 14066306a..90e766588 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -106,6 +106,7 @@ module sfincs_data real*4 factor_pres real*4 factor_prcp real*4 factor_spw_size + real*4 waveforces_ratio ! integer mmax integer nmax @@ -164,6 +165,7 @@ module sfincs_data character*256 :: z0lfile character*256 :: qtrfile character*256 :: volfile + character*256 :: veggiefile ! character*256 :: trefstr_iso8601 character*41 :: treftimefews @@ -192,6 +194,8 @@ module sfincs_data logical :: subgrid logical :: manning2d ! spatially-varying roughness logical :: coriolis + logical :: vegetation + logical :: snapwave_vegetation logical :: store_cumulative_precipitation logical :: store_maximum_waterlevel logical :: store_maximum_waterdepth @@ -206,6 +210,7 @@ module sfincs_data logical :: store_zvolume logical :: store_storagevolume logical :: store_meteo + logical :: store_vegetation logical :: store_wind logical :: store_wind_max logical :: store_wave_forces @@ -419,6 +424,26 @@ module sfincs_data ! real*4, dimension(:), allocatable :: uvmean ! + ! Vegetation + ! + integer :: vegetation_vertical_segments ! nr of vegetation sections in vertical + real*4, dimension(:,:), allocatable :: vegetation_cd + real*4, dimension(:,:), allocatable :: vegetation_stems_height + real*4, dimension(:,:), allocatable :: vegetation_stems_height_uv + real*4, dimension(:,:), allocatable :: vegetation_stems_width + real*4, dimension(:,:), allocatable :: vegetation_stems_density + real*4, dimension(:,:), allocatable :: vegetation_stems_cd_width_density_uv + ! Lookup table arrays (pre-computed in initialize_vegetation, used in compute_fluxes without inner loop) + integer :: vegetation_nlookup ! number of equidistant vertical sections in lookup table (default 20, set via sfincs.inp) + real*4, dimension(:), allocatable :: vegetation_lookup_hmin_uv ! minimum vegetation height on uv points + real*4, dimension(:), allocatable :: vegetation_lookup_hmax_uv ! maximum vegetation height on uv points + real*4, dimension(:), allocatable :: vegetation_lookup_dh_uv ! bin width of lookup table: hmax / vegetation_nlookup, per uv point + real*4, dimension(:,:), allocatable :: vegetation_cd_sum_table ! cumulative sum of cd*width*density at vegetation_nlookup equidistant depth levels, (npuv, 0:vegetation_nlookup) + real*4, dimension(:,:), allocatable :: vegetation_cd_slope_table ! slope between consecutive table entries: table(k+1)-table(k), (npuv, 0:vegetation_nlookup-1) + ! + ! Wave forces limiter determined in sfincs_snapwave + real*4 :: fwmaxfac + ! !!! Wave makers ! character*256 :: wavemaker_wvmfile ! polylines @@ -624,7 +649,10 @@ module sfincs_data real*4, dimension(:), allocatable :: df real*4, dimension(:), allocatable :: dwig real*4, dimension(:), allocatable :: dfig - real*4, dimension(:), allocatable :: cg + real*4, dimension(:), allocatable :: cg + real*4, dimension(:), allocatable :: cgig + real*4, dimension(:), allocatable :: qb + real*4, dimension(:), allocatable :: gam real*4, dimension(:), allocatable :: betamean real*4, dimension(:), allocatable :: srcig real*4, dimension(:), allocatable :: alphaig diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 65bd1015e..fb9da5a6e 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -10,6 +10,7 @@ subroutine initialize_domain() use sfincs_data use quadtree use sfincs_infiltration + use sfincs_vegetation use sfincs_timestep_analysis ! implicit none @@ -28,6 +29,8 @@ subroutine initialize_domain() ! call initialize_storage_volume() ! + call initialize_vegetation() + ! call initialize_hydro() ! if (timestep_analysis) then @@ -2246,6 +2249,12 @@ subroutine initialize_hydro() dfig = 0.0 allocate(cg(np)) cg = 0.0 + allocate(cgig(np)) + cgig = 0.0 + allocate(qb(np)) + qb = 0.0 + allocate(gam(np)) + gam = 0.0 allocate(betamean(np)) betamean = 0.0 allocate(srcig(np)) diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index 5aab1f182..8b447c2a5 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -10,6 +10,7 @@ subroutine read_sfincs_input() use sfincs_date use sfincs_log use sfincs_error + use sfincs_read ! implicit none ! @@ -203,6 +204,7 @@ subroutine read_sfincs_input() call read_logical_input(500, 'bathtub', bathtub, .false.) call read_real_input(500, 'bathtub_fachs', bathtub_fac_hs, 0.2) call read_real_input(500, 'bathtub_dt', bathtub_dt, -999.0) + call read_logical_input(500,'vegetation',vegetation,.false.) ! ! Domain ! @@ -219,6 +221,8 @@ subroutine read_sfincs_input() call read_char_input(500,'manningfile',manningfile,'none') call read_char_input(500,'drnfile',drnfile,'none') call read_char_input(500,'volfile',volfile,'none') + call read_char_input(500,'vegetationfile',veggiefile,'none') + call read_int_input(500,'vegetation_nlookup',vegetation_nlookup,20) ! ! Forcing ! @@ -299,7 +303,9 @@ subroutine read_sfincs_input() percdoneval = max(min(percdoneval,100), 0) ! ! Coupled SnapWave solver related - call read_int_input(500,'snapwave_wind',iwind,0) + call read_int_input(500,'snapwave_wind',iwind,0) + call read_logical_input(500,'snapwave_vegetation',snapwave_vegetation,.false.) + call read_real_input(500,'snapwave_waveforces_ratio',waveforces_ratio,1.0) ! ! Wind drag ! @@ -484,6 +490,14 @@ subroutine read_sfincs_input() endif endif ! + store_vegetation = .false. + if (vegetation==.true. .or. snapwave_vegetation==.true.) then + ! + store_vegetation = .true. + ! vegetation can be used in SnapWave and/or SFINCS calculations + ! + endif + ! store_twet = .false. if (storetwet==1) then store_twet = .true. @@ -733,303 +747,4 @@ subroutine read_sfincs_input() end subroutine - - subroutine read_real_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - real*4, intent(out) :: value - real*4, intent(in) :: default - integer j,stat,ilen - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - read(valstr,*)value - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - subroutine read_real_array_input(fileid,keyword,value,default,nr) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - integer, intent(in) :: nr - real*4, dimension(:), intent(out), allocatable :: value - real*4, intent(in) :: default - integer j,stat, m,ilen - ! - allocate(value(nr)) - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - read(valstr,*)(value(m), m = 1, nr) - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - - subroutine read_int_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - integer, intent(out) :: value - integer, intent(in) :: default - integer j,stat,ilen - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - read(valstr,*)value - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - - subroutine read_char_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr0 - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - character(*), intent(in) :: default - character(*), intent(out) :: value - integer j,stat,ilen,jn - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - value = valstr - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - - subroutine read_logical_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr0 - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - logical, intent(in) :: default - logical, intent(out) :: value - integer j,stat,ilen - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - if (valstr(1:1) == '1' .or. valstr(1:1) == 'y' .or. valstr(1:1) == 'Y' .or. valstr(1:1) == 't' .or. valstr(1:1) == 'T') then - value = .true. - else - value = .false. - endif - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - subroutine read_line(line0, keystr, valstr) - ! - ! Reads line from input file, returns keyword and value strings - ! - character(*), intent(in) :: line0 - character(len=256) :: line - character(*), intent(out) :: keystr - character(*), intent(out) :: valstr - integer j, ilen, jn - ! - keystr = '' - valstr = '' - ! - ! Change tabs into spaces. - ! - call notabs(line0, line, ilen) - ! - ! Look for line ending character. Remove it if it exists. - ! - jn = index(line, '\r') - ! - if (jn > 0) then - ! - ! New line character detected (probably sfincs.inp with windows line endings, running in linux) - ! - line = line(1 : jn - 1) - ! - endif - ! - ! Remove leading and trailing spaces. - ! - line = trim(line) - ! - if (line(1:1) == '#' .or. line(1:1) == '!' .or. line(1:1) == '@') return - ! - ! Find "=" - ! - j = index(line, '=') - ! - if (j == 0) return - ! - keystr = trim(line(1:j-1)) - ! - valstr = trim(line(j+1:)) - ! - ! Remove comments - ! - jn = index(valstr, '#') - ! - if (jn > 0) then - ! - valstr = trim(valstr(1 : jn - 1)) - ! - endif - ! - valstr = adjustl(trim(valstr)) - ! - end subroutine - - - subroutine notabs(INSTR,OUTSTR,ILEN) - ! @(#) convert tabs in input to spaces in output while maintaining columns, assuming a tab is set every 8 characters - ! - ! USES: - ! It is often useful to expand tabs in input files to simplify further processing such as tokenizing an input line. - ! Some FORTRAN compilers hate tabs in input files; some printers; some editors will have problems with tabs - ! AUTHOR: - ! John S. Urban - ! - ! SEE ALSO: - ! GNU/Unix commands expand(1) and unexpand(1) - ! - use ISO_FORTRAN_ENV, only : ERROR_UNIT ! get unit for standard error. if not supported yet, define ERROR_UNIT for your system (typically 0) - character(len=*),intent(in) :: INSTR ! input line to scan for tab characters - character(len=*),intent(out) :: OUTSTR ! tab-expanded version of INSTR produced - integer,intent(out) :: ILEN ! column position of last character put into output string - - integer,parameter :: TABSIZE=8 ! assume a tab stop is set every 8th column - character(len=1) :: c ! character read from stdin - integer :: ipos ! position in OUTSTR to put next character of INSTR - integer :: lenin ! length of input string trimmed of trailing spaces - integer :: lenout ! number of characters output string can hold - integer :: i10 ! counter that advances thru input string INSTR one character at a time - ! - IPOS=1 ! where to put next character in output string OUTSTR - lenin=len(INSTR) ! length of character variable INSTR - lenin=len_trim(INSTR(1:lenin)) ! length of INSTR trimmed of trailing spaces - lenout=len(OUTSTR) ! number of characters output string OUTSTR can hold - OUTSTR=" " ! this SHOULD blank-fill string, a buggy machine required a loop to set all characters - ! - do i10=1,lenin ! look through input string one character at a time - c=INSTR(i10:i10) - if(ichar(c) == 9)then ! test if character is a tab (ADE (ASCII Decimal Equivalent) of tab character is 9) - IPOS = IPOS + (TABSIZE - (mod(IPOS-1,TABSIZE))) - else ! c is anything else other than a tab insert it in output string - if(IPOS > lenout)then - write(ERROR_UNIT,*)"*notabs* output string overflow" - exit - else - OUTSTR(IPOS:IPOS)=c - IPOS=IPOS+1 - endif - endif - enddo - ! - ILEN=len_trim(OUTSTR(:IPOS)) ! trim trailing spaces - return - ! - end subroutine notabs - - -end module +end module \ No newline at end of file diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index 8e3707475..607cdb22f 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -94,8 +94,8 @@ function sfincs_initialize() result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = "$Rev: v2.3.2 mt. Faber+" - build_date = "$Date: 2025-04-10" + build_revision = "$Rev: v2.3.2 mt. Faber+branch:336-added_onto_301" + build_date = "$Date: 2026-05-21" ! call write_log('', 1) call write_log('------------ Welcome to SFINCS ------------', 1) @@ -255,6 +255,11 @@ function sfincs_initialize() result(ierr) else ! call write_log('Non-hydrostatic : no', 1) endif + if (vegetation) then + call write_log('Vegetation : yes', 1) + else + ! call write_log('Vegetation : no', 1) + endif if (bathtub) then call write_log('Bathtub : yes', 1) else diff --git a/source/src/sfincs_momentum.f90 b/source/src/sfincs_momentum.f90 index 4e19f72fa..6d33454c9 100644 --- a/source/src/sfincs_momentum.f90 +++ b/source/src/sfincs_momentum.f90 @@ -21,9 +21,10 @@ subroutine compute_fluxes(dt, tloop) integer :: ip integer :: nm integer :: nmu + integer :: ik_veg integer :: n integer :: m - + ! integer :: idir integer :: iref integer :: itype @@ -89,6 +90,8 @@ subroutine compute_fluxes(dt, tloop) ! real*4 :: min_dt_ip ! + real*4 :: frac_veg + ! real*4, parameter :: expo = 1.0 / 3.0 !integer, parameter :: expo = 1 ! @@ -137,6 +140,7 @@ subroutine compute_fluxes(dt, tloop) !$omp end do !$omp end parallel ! + ! ! Copy flux and velocity from previous time step ! !$acc parallel, present( kcuv, kfuv, zs, q, q0, uv, uv0, zsderv, & @@ -146,11 +150,11 @@ subroutine compute_fluxes(dt, tloop) !$acc uv_index_z_nm, uv_index_z_nmu, uv_index_u_nmd, uv_index_u_nmu, uv_index_u_ndm, uv_index_u_num, & !$acc uv_index_v_ndm, uv_index_v_ndmu, uv_index_v_nm, uv_index_v_nmu, cuv_index_uv, cuv_index_uv1, cuv_index_uv2, & !$acc zb, zbuv, zbuvmx, tauwu, tauwv, patm, fwuv, gn2uv, dxminv, dxrinv, dyrinv, dxm2inv, dxr2inv, dyr2inv, & - !$acc dxrinvc, dyrinvc, fcorio2d, nuvisc, z_volume, gnapp2, x73, timestep_analysis_required_timestep ) num_gangs( 1024 ) vector_length( 128 ) + !$acc dxrinvc, dyrinvc, fcorio2d, nuvisc, z_volume, gnapp2, x73, timestep_analysis_required_timestep, vegetation_cd_sum_table, vegetation_cd_slope_table, vegetation_lookup_hmax_uv, vegetation_lookup_dh_uv ) num_gangs( 1024 ) vector_length( 128 ) !$omp parallel & !$omp private ( ip,hu,qfr,qsm,qx_nm,nm,nmu,dzdx,frc,idir,itype,iref,dxuvinv,dxuv2inv,dyuvinv,dyuv2inv, & !$omp qx_nmd,qx_nmu,qy_nm,qy_ndm,qy_nmu,qy_ndmu,uu_nm,uu_nmd,uu_nmu,uu_num,uu_ndm,vu, & - !$omp fcoriouv,gnavg2,iok,zsu,dzuv,iuv,facint,fwmax,zmax,zmin,one_minus_facint,dqxudx,dqyudy,uu,ud,qu,qd,qy,hwet,phi,adv,mdrv,hu73,min_dt_ip ) & + !$omp fcoriouv,gnavg2,iok,zsu,dzuv,iuv,facint,fwmax,zmax,zmin,one_minus_facint,dqxudx,dqyudy,uu,ud,qu,qd,qy,hwet,phi,adv,mdrv,hu73,min_dt_ip,ik_veg,frac_veg ) & !$omp reduction ( min : min_dt ) !$omp do schedule ( dynamic, 256 ) !$acc loop, reduction( min : min_dt ), gang, vector @@ -597,12 +601,35 @@ subroutine compute_fluxes(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 + ! old: fwmax = 0.8 * hwet * sqrt(hwet) / 15 + ! fix for lab cases: fwmax = 999 + ! + fwmax = fwmaxfac * hwet * sqrt(hwet) + ! Note, fwmaxfac is determined in sfincs_snapwave every 'update_wave_field' call ! - frc = frc + phi * sign(min(abs(fwuv(ip)), fwmax), fwuv(ip)) + frc = frc + phi * sign(min(abs(fwuv(ip)), fwmax), fwuv(ip)) ! endif ! + if (vegetation) then + ! + ! Vegetation drag due to mean flow + ! Direct lookup in pre-computed table - no inner loop over vegetation layers + ! + if (vegetation_lookup_hmax_uv(ip) > 0.0 .and. hwet > 0.0) then + ! + frac_veg = min(hwet, vegetation_lookup_hmax_uv(ip)) / vegetation_lookup_dh_uv(ip) + ! + ik_veg = min(int(frac_veg), vegetation_nlookup - 1) + ! + frac_veg = frac_veg - real(ik_veg) + ! + frc = frc - phi * (vegetation_cd_sum_table(ip, ik_veg) + frac_veg * vegetation_cd_slope_table(ip, ik_veg)) * uv0(ip) * abs(uv0(ip)) + ! + endif + ! + endif + ! ! Compute flux qfr used for friction term ! if (kfuv(ip) == 0) then diff --git a/source/src/sfincs_ncinput.F90 b/source/src/sfincs_ncinput.F90 index 40238e1b6..eadecddd9 100644 --- a/source/src/sfincs_ncinput.F90 +++ b/source/src/sfincs_ncinput.F90 @@ -244,6 +244,46 @@ subroutine read_netcdf_storage_volume() ! end subroutine + subroutine read_netcdf_quadtree_get_dimension(ncfile, varname, var) + ! For instance: vegetationfile, nsec, vegetation_vertical_segments + ! + use netcdf + use sfincs_data + use quadtree + ! + implicit none + ! + integer :: nm, ip, nrcells, status + ! + character*256 :: ncfile + character*256 :: varname + ! + integer, intent(inout) :: var ! variable that we are mapping to + ! + real*4, dimension(:), allocatable :: vartmp + ! + ! Open netcdf file + ! + NF90(nf90_open(trim(ncfile), NF90_CLOBBER, net_file_generic%ncid)) + ! + ! Get dimensions id's: nr points + ! + NF90(nf90_inq_dimid(net_file_generic%ncid, varname, net_file_generic%np_dimid)) + ! + ! Get dimensions sizes + ! + status = nf90_inquire_dimension(net_file_generic%ncid, net_file_generic%np_dimid, len = var) + ! + ! Stop SFINCS if wanted variable was not found + if (status /= nf90_noerr) then + write(logstr,'(a,a,a,a,a)')'Error : netcdf input file ',trim(ncfile),' does not contain needed variable: ',trim(varname),' !' + call stop_sfincs(trim(logstr), 1) + endif + ! + NF90(nf90_close(net_file_generic%ncid)) + ! + end subroutine + subroutine read_netcdf_quadtree_to_sfincs(ncfile, varname, var) ! For instance: storage_volume.nc, vol, storage_volume ! diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index 33f0ba8e3..388bf5405 100644 --- a/source/src/sfincs_ncoutput.F90 +++ b/source/src/sfincs_ncoutput.F90 @@ -26,6 +26,9 @@ module sfincs_ncoutput integer :: manning_varid integer :: pnonh_varid integer :: subgridslope_varid + ! Vegetation + integer :: nsec_dimid + integer :: veg_cd_varid, veg_ah_varid, veg_bstems_varid, veg_Nstems_varid ! integer :: mesh2d_varid integer :: mesh2d_node_x_varid, mesh2d_node_y_varid @@ -55,7 +58,7 @@ module sfincs_ncoutput integer :: patm_varid, wind_speed_varid, wind_dir_varid integer :: inp_varid, total_runtime_varid, average_dt_varid, status_varid integer :: hm0_varid, hm0ig_varid, zsm_varid, tp_varid, tpig_varid, wavdir_varid, dirspr_varid - integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, beta_varid, srcig_varid, alphaig_varid + integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, cgig_varid, beta_varid, srcig_varid, alphaig_varid, qb_varid, gam_varid integer :: runup_gauge_name_varid, runup_gauge_zs_varid ! end type @@ -869,7 +872,7 @@ subroutine ncoutput_quadtree_map_init() ! implicit none ! - integer :: nm, nmq, nmu1, num1, n, m, nn, ntmx, n_nodes, n_faces, iref + integer :: nm, nmq, nmu1, num1, n, m, nn, ntmx, n_nodes, n_faces, iref, isec real*4 :: dxx, dyy ! real, dimension(:), allocatable :: nodes_x @@ -877,6 +880,7 @@ subroutine ncoutput_quadtree_map_init() integer*4, dimension(:,:), allocatable :: face_nodes real*4, dimension(:), allocatable :: vtmp integer*4, dimension(:), allocatable :: vtmpi + real*4, dimension(:,:), allocatable :: vtmp2d ! ! Very lazy for now ! @@ -943,11 +947,15 @@ subroutine ncoutput_quadtree_map_init() ! Time ! NF90(nf90_def_dim(map_file%ncid, 'time', NF90_UNLIMITED, map_file%time_dimid)) ! time - ntmx = max(ceiling((t1out - t0out)/dtmaxout), 1) + ntmx = max(ceiling((t1out - t0out)/dtmaxout), 1) NF90(nf90_def_dim(map_file%ncid, 'timemax', ntmx, map_file%timemax_dimid)) ! time - NF90(nf90_def_dim(map_file%ncid, 'runtime', 1, map_file%runtime_dimid)) ! total_runtime, average_dt + NF90(nf90_def_dim(map_file%ncid, 'runtime', 1, map_file%runtime_dimid)) ! total_runtime, average_dt ! - ! Some metadata attributes + if (store_vegetation) then + NF90(nf90_def_dim(map_file%ncid, 'vegetation_vertical_segments', vegetation_vertical_segments, map_file%nsec_dimid)) ! number of vegetation vertical sections + endif + ! + ! Some metadata attributes ! NF90(nf90_put_att(map_file%ncid,nf90_global, "Conventions", "Conventions = 'CF-1.8 UGRID-1.0 Deltares-0.10'")) NF90(nf90_put_att(map_file%ncid,nf90_global, "Build-Revision-Date-Netcdf-library", trim(nf90_inq_libvers()))) ! version of netcdf library @@ -1064,9 +1072,41 @@ subroutine ncoutput_quadtree_map_init() NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'units', '-')) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'standard_name', 'mask')) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'long_name', 'msk_active_cells')) - NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'description', 'inactive=0, active=1, normal_boundary=2, outflow_boundary=3, wavemaker=4')) + NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'description', 'inactive=0, active=1, normal_boundary=2, outflow_boundary=3, wavemaker=4')) + ! + if (store_vegetation) then + ! + NF90(nf90_def_var(map_file%ncid, 'vegetation_cd', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%nsec_dimid/), map_file%veg_cd_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%veg_cd_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_cd_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_cd_varid, 'units', '-')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_cd_varid, 'standard_name', 'vegetation_cd')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_cd_varid, 'long_name', 'bulk_drag_coefficient_per_vegetation_section')) + ! + NF90(nf90_def_var(map_file%ncid, 'vegetation_stems_height', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%nsec_dimid/), map_file%veg_ah_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%veg_ah_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_ah_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_ah_varid, 'units', 'm')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_ah_varid, 'standard_name', 'vegetation_stems_height')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_ah_varid, 'long_name', 'vegetation_section_thickness')) + ! + NF90(nf90_def_var(map_file%ncid, 'vegetation_stems_width', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%nsec_dimid/), map_file%veg_bstems_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%veg_bstems_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_bstems_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_bstems_varid, 'units', 'm')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_bstems_varid, 'standard_name', 'vegetation_stems_width')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_bstems_varid, 'long_name', 'width_of_individual_vegetation_stems_per_section')) + ! + NF90(nf90_def_var(map_file%ncid, 'vegetation_stems_density', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%nsec_dimid/), map_file%veg_Nstems_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%veg_Nstems_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_Nstems_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%veg_Nstems_varid, 'units', 'm-2')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_Nstems_varid, 'standard_name', 'vegetation_stems_density')) + NF90(nf90_put_att(map_file%ncid, map_file%veg_Nstems_varid, 'long_name', 'number_of_stems_per_unit_horizontal_area_per_section')) + ! + endif ! - ! Time variables + ! Time variables ! trefstr_iso8601 = date_to_iso8601(trefstr) ! @@ -1599,13 +1639,67 @@ subroutine ncoutput_quadtree_map_init() ! endif ! + ! Write vegetation fields (static, written once at init) + ! + if (store_vegetation) then + ! + allocate(vtmp2d(n_faces, vegetation_vertical_segments)) + ! + vtmp2d = FILL_VALUE + do nmq = 1, quadtree_nr_points + nm = index_sfincs_in_quadtree(nmq) + if (nm > 0) then + do isec = 1, vegetation_vertical_segments + vtmp2d(nmq, isec) = vegetation_cd(nm, isec) + enddo + endif + enddo + NF90(nf90_put_var(map_file%ncid, map_file%veg_cd_varid, vtmp2d)) + ! + vtmp2d = FILL_VALUE + do nmq = 1, quadtree_nr_points + nm = index_sfincs_in_quadtree(nmq) + if (nm > 0) then + do isec = 1, vegetation_vertical_segments + vtmp2d(nmq, isec) = vegetation_stems_height(nm, isec) + enddo + endif + enddo + NF90(nf90_put_var(map_file%ncid, map_file%veg_ah_varid, vtmp2d)) + ! + vtmp2d = FILL_VALUE + do nmq = 1, quadtree_nr_points + nm = index_sfincs_in_quadtree(nmq) + if (nm > 0) then + do isec = 1, vegetation_vertical_segments + vtmp2d(nmq, isec) = vegetation_stems_width(nm, isec) + enddo + endif + enddo + NF90(nf90_put_var(map_file%ncid, map_file%veg_bstems_varid, vtmp2d)) + ! + vtmp2d = FILL_VALUE + do nmq = 1, quadtree_nr_points + nm = index_sfincs_in_quadtree(nmq) + if (nm > 0) then + do isec = 1, vegetation_vertical_segments + vtmp2d(nmq, isec) = vegetation_stems_density(nm, isec) + enddo + endif + enddo + NF90(nf90_put_var(map_file%ncid, map_file%veg_Nstems_varid, vtmp2d)) + ! + deallocate(vtmp2d) + ! + endif + ! ! write away intermediate data ! NF90(nf90_sync(map_file%ncid)) !write away intermediate data ! end subroutine - - + + subroutine ncoutput_his_init() ! ! 1. Initialise dimensions/variables/attributes @@ -1959,6 +2053,27 @@ subroutine ncoutput_his_init() NF90(nf90_put_att(his_file%ncid, his_file%cg_varid, 'long_name', 'wave group velocity')) NF90(nf90_put_att(his_file%ncid, his_file%cg_varid, 'coordinates', 'station_id station_name point_x point_y')) ! + NF90(nf90_def_var(his_file%ncid, 'point_cgig', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%cgig_varid)) ! time-varying water level point + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'units', 'm/s')) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'standard_name', 'infragravity_wave_velocity')) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'long_name', 'infragravity wave velocity')) + NF90(nf90_put_att(his_file%ncid, his_file%cgig_varid, 'coordinates', 'station_id station_name point_x point_y')) + ! + NF90(nf90_def_var(his_file%ncid, 'point_qb', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%qb_varid)) ! time-varying water level point + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'units', '-')) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'standard_name', 'fraction_breaking_waves')) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'long_name', 'fraction breaking incident waves')) + NF90(nf90_put_att(his_file%ncid, his_file%qb_varid, 'coordinates', 'station_id station_name point_x point_y')) + ! + NF90(nf90_def_var(his_file%ncid, 'point_gam', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%gam_varid)) ! time-varying water level point + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'units', '-')) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'standard_name', 'local_wave_height_water_depth_ratio')) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'long_name', 'local wave height water depth ratio')) + NF90(nf90_put_att(his_file%ncid, his_file%gam_varid, 'coordinates', 'station_id station_name point_x point_y')) + ! NF90(nf90_def_var(his_file%ncid, 'point_beta', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%beta_varid)) ! time-varying water level point NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, 'units', '-')) @@ -3038,6 +3153,9 @@ subroutine ncoutput_update_his(t,nthisout) real*4, dimension(nobs) :: dwigobs real*4, dimension(nobs) :: dfigobs real*4, dimension(nobs) :: cgobs + real*4, dimension(nobs) :: cgigobs + real*4, dimension(nobs) :: qbobs + real*4, dimension(nobs) :: gamobs real*4, dimension(nobs) :: betaobs real*4, dimension(nobs) :: srcigobs real*4, dimension(nobs) :: alphaigobs @@ -3062,6 +3180,9 @@ subroutine ncoutput_update_his(t,nthisout) dwobs = FILL_VALUE dfobs = FILL_VALUE cgobs = FILL_VALUE + cgigobs = FILL_VALUE + qbobs = FILL_VALUE + gamobs = FILL_VALUE betaobs = FILL_VALUE srcigobs = FILL_VALUE alphaigobs = FILL_VALUE @@ -3173,7 +3294,10 @@ subroutine ncoutput_update_his(t,nthisout) dfobs(iobs) = df(nm) dwigobs(iobs) = dwig(nm) dfigobs(iobs) = dfig(nm) - cgobs(iobs) = cg(nm) + cgobs(iobs) = cg(nm) + cgigobs(iobs) = cgig(nm) + qbobs(iobs) = qb(nm) + gamobs(iobs) = gam(nm) betaobs(iobs) = betamean(nm) srcigobs(iobs) = srcig(nm) alphaigobs(iobs) = alphaig(nm) @@ -3236,6 +3360,9 @@ subroutine ncoutput_update_his(t,nthisout) NF90(nf90_put_var(his_file%ncid, his_file%dfig_varid, dfigobs, (/1, nthisout/))) ! NF90(nf90_put_var(his_file%ncid, his_file%cg_varid, cgobs, (/1, nthisout/))) + NF90(nf90_put_var(his_file%ncid, his_file%cgig_varid, cgigobs, (/1, nthisout/))) + NF90(nf90_put_var(his_file%ncid, his_file%qb_varid, qbobs, (/1, nthisout/))) + NF90(nf90_put_var(his_file%ncid, his_file%gam_varid, gamobs, (/1, nthisout/))) ! NF90(nf90_put_var(his_file%ncid, his_file%beta_varid, betaobs, (/1, nthisout/))) NF90(nf90_put_var(his_file%ncid, his_file%srcig_varid, srcigobs, (/1, nthisout/))) @@ -3913,8 +4040,8 @@ subroutine ncoutput_add_params(ncid, varid) use sfincs_data ! ! Because of overlapping names, only important specific values from snapwave_data - use snapwave_data, only: gamma, gammax, alpha, hmin, fw0, fw0_ig, dt, tol, dtheta, crit, nr_sweeps, baldock_opt, baldock_ratio, & - igwaves_opt, alpha_ig, gamma_ig, shinc2ig, alphaigfac, baldock_ratio_ig, ig_opt, herbers_opt, tpig_opt, eeinc2ig, tinc2ig, & + use snapwave_data, only: gamma, gammax, alpha, hmin, fw0, fw0_ig, dt, tol, dtheta, crit, nr_sweeps, baldock_exponent, baldock_ratio, & + igwaves_opt, alpha_ig, gamma_ig, gamma_fac_br, shinc2ig, alphaigfac, baldock_ratio_ig, ig_opt, herbers_opt, tpig_opt, eeinc2ig, tinc2ig, & snapwave_jonswapfile, snapwave_encfile, snapwave_bndfile, snapwave_bhsfile, snapwave_btpfile, snapwave_bwdfile, snapwave_bdsfile, upwfile, gridfile ! @@ -4115,14 +4242,16 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'snapwave_dtheta',dtheta)) NF90(nf90_put_att(ncid, varid, 'snapwave_crit',crit)) NF90(nf90_put_att(ncid, varid, 'snapwave_nrsweeps',nr_sweeps)) - NF90(nf90_put_att(ncid, varid, 'snapwave_baldock_opt',baldock_opt)) + NF90(nf90_put_att(ncid, varid, 'snapwave_baldock_exponent',baldock_exponent)) NF90(nf90_put_att(ncid, varid, 'snapwave_baldock_ratio',baldock_ratio)) + NF90(nf90_put_att(ncid, varid, 'snapwave_waveforces_ratio',waveforces_ratio)) ! ! SnapWave IG ! NF90(nf90_put_att(ncid, varid, 'snapwave_igwaves',igwaves_opt)) NF90(nf90_put_att(ncid, varid, 'snapwave_alpha_ig',alpha_ig)) - NF90(nf90_put_att(ncid, varid, 'snapwave_gammaig',gamma_ig)) + NF90(nf90_put_att(ncid, varid, 'snapwave_gammaig',gamma_ig)) + NF90(nf90_put_att(ncid, varid, 'snapwave_gamma_fac_br',gamma_fac_br)) NF90(nf90_put_att(ncid, varid, 'snapwave_shinc2ig',shinc2ig)) NF90(nf90_put_att(ncid, varid, 'snapwave_alphaigfac',alphaigfac)) NF90(nf90_put_att(ncid, varid, 'snapwave_baldock_ratio_ig',baldock_ratio_ig)) diff --git a/source/src/sfincs_openacc.f90 b/source/src/sfincs_openacc.f90 index 857260203..fbd3e3875 100644 --- a/source/src/sfincs_openacc.f90 +++ b/source/src/sfincs_openacc.f90 @@ -36,7 +36,8 @@ subroutine initialize_openacc() !$acc gnapp2, & !$acc timestep_analysis_required_timestep, timestep_analysis_average_required_timestep, timestep_analysis_times_wet, timestep_analysis_times_limiting, & !$acc qinffield, qinfmap, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, & - !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0 ) + !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0, & + !$acc vegetation_cd_sum_table, vegetation_cd_slope_table, vegetation_lookup_hmin_uv, vegetation_lookup_hmax_uv, vegetation_lookup_dh_uv ) ! end subroutine ! @@ -66,7 +67,8 @@ subroutine finalize_openacc() !$acc gnapp2, & !$acc timestep_analysis_required_timestep, timestep_analysis_average_required_timestep, timestep_analysis_times_wet, timestep_analysis_times_limiting, & !$acc qinffield, qinfmap, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, & - !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0 ) + !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0, & + !$acc vegetation_cd_sum_table, vegetation_cd_slope_table, vegetation_lookup_hmin_uv, vegetation_lookup_hmax_uv, vegetation_lookup_dh_uv ) ! end subroutine finalize_openacc ! diff --git a/source/src/sfincs_quadtree.F90 b/source/src/sfincs_quadtree.F90 index 72a5397a9..28e197a8d 100644 --- a/source/src/sfincs_quadtree.F90 +++ b/source/src/sfincs_quadtree.F90 @@ -47,9 +47,10 @@ module quadtree integer*1, dimension(:), allocatable :: quadtree_snapwave_mask integer*1, dimension(:), allocatable :: quadtree_nonh_mask ! + ! type net_type_qtr integer :: ncid - integer :: np_dimid + integer :: np_dimid, nsec_dimid integer :: n_varid, m_varid integer :: level_varid integer :: nu_varid, mu_varid, nd_varid, md_varid @@ -301,7 +302,7 @@ subroutine quadtree_read_file_netcdf(qtrfile, snapwave, nonhydrostatic) logical, intent(in) :: snapwave, nonhydrostatic ! integer*1 :: iversion - integer :: np, ip, iepsg, status + integer :: np, nm, ip, iepsg, status ! write(logstr,'(a,a)')'Info : reading QuadTree netCDF file ', trim(qtrfile) call write_log(logstr, 0) @@ -343,7 +344,7 @@ subroutine quadtree_read_file_netcdf(qtrfile, snapwave, nonhydrostatic) NF90(nf90_inq_varid(net_file_qtr%ncid, 'snapwave_mask', net_file_qtr%snapwave_mask_varid)) ! allocate(quadtree_snapwave_mask(np)) - ! + ! endif ! ! Allocate variables @@ -391,7 +392,9 @@ subroutine quadtree_read_file_netcdf(qtrfile, snapwave, nonhydrostatic) NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%mask_varid, quadtree_mask(:))) ! if (snapwave) then + ! NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%snapwave_mask_varid, quadtree_snapwave_mask(:))) + ! endif ! ! Nonhydrostatic mask diff --git a/source/src/sfincs_read.f90 b/source/src/sfincs_read.f90 new file mode 100644 index 000000000..553c82fac --- /dev/null +++ b/source/src/sfincs_read.f90 @@ -0,0 +1,303 @@ +module sfincs_read + +contains + + subroutine read_real_input(fileid,keyword,value,default) + ! + character(*), intent(in) :: keyword + character(len=256) :: keystr + character(len=256) :: valstr + character(len=256) :: line + integer, intent(in) :: fileid + real*4, intent(out) :: value + real*4, intent(in) :: default + integer j,stat,ilen + ! + value = default + ! + rewind(fileid) + ! + do while(.true.) + ! + read(fileid,'(a)',iostat = stat)line + ! + if (stat==-1) exit + ! + call read_line(line, keystr, valstr) + ! + if (trim(keystr)==trim(keyword)) then + ! + read(valstr,*)value + ! + exit + ! + endif + ! + enddo + ! + end subroutine + + subroutine read_real_array_input(fileid,keyword,value,default,nr) + ! + character(*), intent(in) :: keyword + character(len=256) :: keystr + character(len=256) :: valstr + character(len=256) :: line + integer, intent(in) :: fileid + integer, intent(in) :: nr + real*4, dimension(:), intent(out), allocatable :: value + real*4, intent(in) :: default + integer j,stat, m,ilen + ! + allocate(value(nr)) + ! + value = default + ! + rewind(fileid) + ! + do while(.true.) + ! + read(fileid,'(a)',iostat = stat)line + ! + if (stat==-1) exit + ! + call read_line(line, keystr, valstr) + ! + if (trim(keystr)==trim(keyword)) then + ! + read(valstr,*)(value(m), m = 1, nr) + ! + exit + ! + endif + ! + enddo + ! + end subroutine + + + subroutine read_int_input(fileid,keyword,value,default) + ! + character(*), intent(in) :: keyword + character(len=256) :: keystr + character(len=256) :: valstr + character(len=256) :: line + integer, intent(in) :: fileid + integer, intent(out) :: value + integer, intent(in) :: default + integer j,stat,ilen + ! + value = default + ! + rewind(fileid) + ! + do while(.true.) + ! + read(fileid,'(a)',iostat = stat)line + ! + if (stat==-1) exit + ! + call read_line(line, keystr, valstr) + ! + if (trim(keystr)==trim(keyword)) then + ! + read(valstr,*)value + ! + exit + ! + endif + ! + enddo + ! + end subroutine + + + subroutine read_char_input(fileid,keyword,value,default) + ! + character(*), intent(in) :: keyword + character(len=256) :: keystr0 + character(len=256) :: keystr + character(len=256) :: valstr + character(len=256) :: line + integer, intent(in) :: fileid + character(*), intent(in) :: default + character(*), intent(out) :: value + integer j,stat,ilen,jn + ! + value = default + ! + rewind(fileid) + ! + do while(.true.) + ! + read(fileid,'(a)',iostat = stat)line + ! + if (stat==-1) exit + ! + call read_line(line, keystr, valstr) + ! + if (trim(keystr)==trim(keyword)) then + ! + value = valstr + ! + exit + ! + endif + ! + enddo + ! + end subroutine + + + subroutine read_logical_input(fileid,keyword,value,default) + ! + character(*), intent(in) :: keyword + character(len=256) :: keystr0 + character(len=256) :: keystr + character(len=256) :: valstr + character(len=256) :: line + integer, intent(in) :: fileid + logical, intent(in) :: default + logical, intent(out) :: value + integer j,stat,ilen + ! + value = default + ! + rewind(fileid) + ! + do while(.true.) + ! + read(fileid,'(a)',iostat = stat)line + ! + if (stat==-1) exit + ! + call read_line(line, keystr, valstr) + ! + if (trim(keystr)==trim(keyword)) then + ! + if (valstr(1:1) == '1' .or. valstr(1:1) == 'y' .or. valstr(1:1) == 'Y' .or. valstr(1:1) == 't' .or. valstr(1:1) == 'T') then + value = .true. + else + value = .false. + endif + ! + exit + ! + endif + ! + enddo + ! + end subroutine + + subroutine read_line(line0, keystr, valstr) + ! + ! Reads line from input file, returns keyword and value strings + ! + character(*), intent(in) :: line0 + character(len=256) :: line + character(*), intent(out) :: keystr + character(*), intent(out) :: valstr + integer j, ilen, jn + ! + keystr = '' + valstr = '' + ! + ! Change tabs into spaces. + ! + call notabs(line0, line, ilen) + ! + ! Look for line ending character. Remove it if it exists. + ! + jn = index(line, '\r') + ! + if (jn > 0) then + ! + ! New line character detected (probably sfincs.inp with windows line endings, running in linux) + ! + line = line(1 : jn - 1) + ! + endif + ! + ! Remove leading and trailing spaces. + ! + line = trim(line) + ! + if (line(1:1) == '#' .or. line(1:1) == '!' .or. line(1:1) == '@') return + ! + ! Find "=" + ! + j = index(line, '=') + ! + if (j == 0) return + ! + keystr = trim(line(1:j-1)) + ! + valstr = trim(line(j+1:)) + ! + ! Remove comments + ! + jn = index(valstr, '#') + ! + if (jn > 0) then + ! + valstr = trim(valstr(1 : jn - 1)) + ! + endif + ! + valstr = adjustl(trim(valstr)) + ! + end subroutine + + + subroutine notabs(INSTR,OUTSTR,ILEN) + ! @(#) convert tabs in input to spaces in output while maintaining columns, assuming a tab is set every 8 characters + ! + ! USES: + ! It is often useful to expand tabs in input files to simplify further processing such as tokenizing an input line. + ! Some FORTRAN compilers hate tabs in input files; some printers; some editors will have problems with tabs + ! AUTHOR: + ! John S. Urban + ! + ! SEE ALSO: + ! GNU/Unix commands expand(1) and unexpand(1) + ! + use ISO_FORTRAN_ENV, only : ERROR_UNIT ! get unit for standard error. if not supported yet, define ERROR_UNIT for your system (typically 0) + character(len=*),intent(in) :: INSTR ! input line to scan for tab characters + character(len=*),intent(out) :: OUTSTR ! tab-expanded version of INSTR produced + integer,intent(out) :: ILEN ! column position of last character put into output string + + integer,parameter :: TABSIZE=8 ! assume a tab stop is set every 8th column + character(len=1) :: c ! character read from stdin + integer :: ipos ! position in OUTSTR to put next character of INSTR + integer :: lenin ! length of input string trimmed of trailing spaces + integer :: lenout ! number of characters output string can hold + integer :: i10 ! counter that advances thru input string INSTR one character at a time + ! + IPOS=1 ! where to put next character in output string OUTSTR + lenin=len(INSTR) ! length of character variable INSTR + lenin=len_trim(INSTR(1:lenin)) ! length of INSTR trimmed of trailing spaces + lenout=len(OUTSTR) ! number of characters output string OUTSTR can hold + OUTSTR=" " ! this SHOULD blank-fill string, a buggy machine required a loop to set all characters + ! + do i10=1,lenin ! look through input string one character at a time + c=INSTR(i10:i10) + if(ichar(c) == 9)then ! test if character is a tab (ADE (ASCII Decimal Equivalent) of tab character is 9) + IPOS = IPOS + (TABSIZE - (mod(IPOS-1,TABSIZE))) + else ! c is anything else other than a tab insert it in output string + if(IPOS > lenout)then + write(ERROR_UNIT,*)"*notabs* output string overflow" + exit + else + OUTSTR(IPOS:IPOS)=c + IPOS=IPOS+1 + endif + endif + enddo + ! + ILEN=len_trim(OUTSTR(:IPOS)) ! trim trailing spaces + return + ! + end subroutine notabs + + +end module diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index 66825a0f5..8d1ccf284 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -27,6 +27,9 @@ module sfincs_snapwave real*4, dimension(:), allocatable :: snapwave_Dwig real*4, dimension(:), allocatable :: snapwave_Dfig real*4, dimension(:), allocatable :: snapwave_cg + real*4, dimension(:), allocatable :: snapwave_cgig + real*4, dimension(:), allocatable :: snapwave_qb + real*4, dimension(:), allocatable :: snapwave_gam real*4, dimension(:), allocatable :: snapwave_beta real*4, dimension(:), allocatable :: snapwave_srcig real*4, dimension(:), allocatable :: snapwave_alphaig @@ -37,6 +40,7 @@ module sfincs_snapwave real*4 :: snapwave_hsmean real*4 :: snapwave_tpmean real*4 :: snapwave_tpigmean + real*4 :: snapwave_fwmaxfac ! contains ! @@ -305,8 +309,10 @@ subroutine update_wave_field(t, tloop) real*4, dimension(:), allocatable :: df0 real*4, dimension(:), allocatable :: dwig0 real*4, dimension(:), allocatable :: dfig0 - real*4, dimension(:), allocatable :: cg0 - !real*4, dimension(:), allocatable :: qb0 + real*4, dimension(:), allocatable :: cg0 + real*4, dimension(:), allocatable :: cgig0 + real*4, dimension(:), allocatable :: qb0 + real*4, dimension(:), allocatable :: gam0 real*4, dimension(:), allocatable :: beta0 real*4, dimension(:), allocatable :: srcig0 real*4, dimension(:), allocatable :: alphaig0 @@ -321,8 +327,10 @@ subroutine update_wave_field(t, tloop) allocate(df0(np)) allocate(dwig0(np)) allocate(dfig0(np)) - allocate(cg0(np)) - !allocate(qb0(np)) + allocate(cg0(np)) + allocate(cgig0(np)) + allocate(qb0(np)) + allocate(gam0(np)) allocate(beta0(np)) allocate(srcig0(np)) allocate(alphaig0(np)) @@ -334,7 +342,9 @@ subroutine update_wave_field(t, tloop) dwig0 = 0.0 dfig0 = 0.0 cg0 = 0.0 - !qb0 = 0.0 + cgig0 = 0.0 + qb0 = 0.0 + gam0 = 0.0 beta0 = 0.0 srcig0 = 0.0 alphaig0 = 0.0 @@ -430,7 +440,9 @@ 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) + cgig0(nm) = snapwave_cgig(ip) + qb0(nm) = snapwave_qb(ip) + gam0(nm) = snapwave_gam(ip) beta0(nm) = snapwave_beta(ip) srcig0(nm) = snapwave_srcig(ip) alphaig0(nm) = snapwave_alphaig(ip) @@ -454,7 +466,9 @@ subroutine update_wave_field(t, tloop) dwig0(nm) = 0.0 dfig0(nm) = 0.0 cg0(nm) = 0.0 - !qb0(nm) = 0.0 + cgig0(nm) = 0.0 + qb0(nm) = 0.0 + gam0(nm) = 0.0 beta0(nm) = 0.0 srcig0(nm) = 0.0 alphaig0(nm) = 0.0 @@ -473,8 +487,10 @@ subroutine update_wave_field(t, tloop) df(nm) = df0(nm) dwig(nm) = dwig0(nm) dfig(nm) = dfig0(nm) - cg(nm) = cg0(nm) - !qb(nm) = qb0(nm) + cg(nm) = cg0(nm) + cgig(nm) = cgig0(nm) + qb(nm) = qb0(nm) + gam(nm) = gam0(nm) betamean(nm) = beta0(nm) srcig(nm) = srcig0(nm) alphaig(nm) = alphaig0(nm) @@ -498,13 +514,13 @@ subroutine update_wave_field(t, tloop) ! ! U point ! - fwuv(ip) = (0.5*(cosrot*fwx0(nm) + sinrot*fwy0(nm)) + 0.5*( cosrot*fwx0(nmu) + sinrot*fwy0(nmu)))/rhow - ! + fwuv(ip) = waveforces_ratio * (0.5 * (cosrot * fwx0(nm) + sinrot * fwy0(nm)) + 0.5 * ( cosrot * fwx0(nmu) + sinrot * fwy0(nmu))) / rhow + ! waveforces_ratio = 1.0 by default, but can be set to 0 to avoid double counting incident setup if wavemaker_hinc true else ! ! V point ! - fwuv(ip) = (0.5*(-sinrot*fwx0(nm) + cosrot*fwy0(nm)) + 0.5*(-sinrot*fwx0(nmu) + cosrot*fwy0(nmu)))/rhow + fwuv(ip) = waveforces_ratio * (0.5 * (-sinrot * fwx0(nm) + cosrot * fwy0(nm)) + 0.5 * (-sinrot * fwx0(nmu) + cosrot * fwy0(nmu))) / rhow ! endif ! @@ -512,6 +528,9 @@ subroutine update_wave_field(t, tloop) ! !$acc update device(fwuv) ! + ! Set wave forces fwmaxfac factor + fwmaxfac = snapwave_fwmaxfac + ! call system_clock(count1, count_rate, count_max) tloop = tloop + 1.0*(count1 - count0)/count_rate ! @@ -551,14 +570,22 @@ subroutine compute_snapwave(t) snapwave_Dwig = Dw_ig snapwave_Dfig = Df_ig snapwave_cg = cg + snapwave_cgig = cg_ig + snapwave_qb = qb + snapwave_gam = gam snapwave_beta = beta snapwave_srcig = srcig snapwave_alphaig = alphaig - ! + ! ! Convert wave force to correct unit [Dw/C] as expected by SFINCS, assumed to be piecewise (seems to work) snapwave_Fx = Fx * rho * depth snapwave_Fy = Fy * rho * depth ! + ! Pre-alculate wave forces limiter factor + snapwave_fwmaxfac = 0.25 * sqrt(g) * rho * gammax**2 / tpmean_bwv + ! + ! FIXME - should we limit snapwave_fwmaxfac to a certain range? + ! ! Loop over points and set Tp, cg, direction, spreading to 0 where H and/or H_ig are zero ! TL: needed because e.g. Tp is set to Tpini initially, so shows values even if cell remains dry with H=0 do k = 1, no_nodes @@ -567,10 +594,13 @@ subroutine compute_snapwave(t) snapwave_mean_direction(k) = 0.0 snapwave_directional_spreading(k) = 0.0 snapwave_cg(k) = 0.0 + snapwave_qb(k) = 0.0 + snapwave_gam(k) = 0.0 endif ! if (snapwave_H_ig(k) <= 0.0) then - snapwave_Tp_ig(k) = 0.0 + snapwave_Tp_ig(k) = 0.0 + snapwave_cgig(k) = 0.0 endif enddo ! @@ -583,17 +613,8 @@ subroutine compute_snapwave(t) if (igwaves) then ! snapwave_tpigmean = tpmean_bwv_ig - ! - if (snapwave_tpigmean < 10.0) then - ! These warnings should not occur here - write(logstr,*)'DEBUG SFINCS_SnapWave - incoming tp for IG wave at wavemaker might be unrealistically small! value: ',snapwave_tpigmean - call write_log(logstr, 0) - elseif (snapwave_tpigmean > 250.0) then - write(logstr,*)'DEBUG SFINCS_SnapWave - incoming tp for IG wave at wavemaker might be unrealistically large! value: ',snapwave_tpigmean - call write_log(logstr, 0) - endif + ! endif - ! TL: NOTE - in first timestep run of SnapWave tp = 0, therefore excluded that case from the check ! end subroutine @@ -604,6 +625,7 @@ subroutine read_snapwave_input() ! Reads snapwave data from sfincs.inp ! use snapwave_data + use sfincs_read ! implicit none ! @@ -611,49 +633,58 @@ subroutine read_snapwave_input() ! ! Input section ! - 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) ! Compute IG waves (1=default), or not (0) - 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) ! Wave breaking parameter for IG waves, default=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) ! ! option controlling from what depth wave breaking should take place for IG waves, default baldock_ratio_ig=0.2 - call read_int_input(500, 'snapwave_ig_opt', ig_opt, 1) ! option of IG wave settings (1 = default = conservative shoaling based dSxx as in Leijnse et al. 2024) - 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), or effectively based on previous timestep (faster, potential mismatch, =0=default) - ! - ! 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 - call read_real_input(500, 'snapwave_eeinc2ig', eeinc2ig, 0.01) ! Only used if snapwave_use_herbers = 0 - call read_real_input(500, 'snapwave_Tinc2ig', Tinc2ig, 7.0) ! Only used if snapwave_use_herbers = 0 + 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) + 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_int_input(500,'snapwave_baldock_exponent',baldock_exponent,0) ! Exponent for multiplying the Baldock dissipation with a factor 'f = (Hloc / Hmax)**iexp' to enhance breaking when H > Hmax, with iexp = 0 (default, means unused), 1 or 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) + call read_real_input(500,'snapwave_relax_factor_DoverA',relax_factor_DoverA,0.25) ! underrelaxation factor for DoverA (set to 1.0 to disable) + call read_real_input(500,'snapwave_relax_factor_DoverE',relax_factor_DoverE,0.25) ! underrelaxation factor for DoverE (set to 1.0 to disable) + ! + ! 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.7) ! Wave breaking parameter for IG waves, default=0.7 + call read_real_input(500,'snapwave_gamma_fac_br',gamma_fac_br,0.45) ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point + 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 steep slope related: + call read_real_input(500,'snapwave_steep_fac1',steep_fac1,0.0) ! Cut-off gamma, below this alphaig_steep = 0 + call read_real_input(500,'snapwave_steep_fac2',steep_fac2,0.1) ! Multiplication factor + call read_real_input(500,'snapwave_steep_fac3',steep_fac3,0.07) ! Cut-off beta, below this alphaig_steep = 0, and above this it increases with beta + call read_real_input(500,'snapwave_steep_fac4',steep_fac4,0.6) ! Exponent + call read_real_input(500,'snapwave_steep_fac5',steep_fac5,1.0) ! Cut-off gamma, above this alphaig_steep = 0 + ! + ! 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 + call read_real_input(500,'snapwave_eeinc2ig',eeinc2ig,0.01) ! Only used if snapwave_use_herbers = 0 + call read_real_input(500,'snapwave_Tinc2ig',Tinc2ig,7.0) ! Only used if snapwave_use_herbers = 0 ! ! Wind ! @@ -661,7 +692,7 @@ subroutine read_snapwave_input() ! ! Vegetation input ! - call read_int_input(500, 'vegetation', vegetation_opt, 0) + call read_int_input(500, 'snapwave_vegetation', vegetation_opt, 0) ! ! Input files ! @@ -677,6 +708,7 @@ subroutine read_snapwave_input() call read_char_input(500, 'snapwave_depfile', depfile, 'none') call read_char_input(500, 'snapwave_ncfile', gridfile, 'snapwave_net.nc') call read_char_input(500, 'netsnapwavefile', netsnapwavefile, 'none') + call read_logical_input(500,'storesnapwavegrid',storesnapwavegrid,.false.) call read_char_input(500, 'tref', trefstr, '20000101 000000') ! Read again > needed in sfincs_ncinput.F90 ! close(500) @@ -729,304 +761,7 @@ subroutine read_snapwave_input() restart = .true. coupled_to_sfincs = .true. ! - end subroutine - - - - subroutine read_real_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - real*4, intent(out) :: value - real*4, intent(in) :: default - integer j,stat,ilen - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - read(valstr,*)value - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - subroutine read_real_array_input(fileid,keyword,value,default,nr) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - integer, intent(in) :: nr - real*4, dimension(:), intent(out), allocatable :: value - real*4, intent(in) :: default - integer j,stat, m,ilen - ! - allocate(value(nr)) - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - read(valstr,*)(value(m), m = 1, nr) - ! - exit - ! - endif - ! - enddo - ! - end subroutine + end subroutine read_snapwave_input - subroutine read_int_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - integer, intent(out) :: value - integer, intent(in) :: default - integer j,stat,ilen - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - read(valstr,*)value - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - - subroutine read_char_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr0 - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - character(*), intent(in) :: default - character(*), intent(out) :: value - integer j,stat,ilen,jn - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - value = valstr - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - subroutine read_logical_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr0 - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - logical, intent(in) :: default - logical, intent(out) :: value - integer j,stat,ilen - ! - value = default - ! - rewind(fileid) - ! - do while(.true.) - ! - read(fileid,'(a)',iostat = stat)line - ! - if (stat==-1) exit - ! - call read_line(line, keystr, valstr) - ! - if (trim(keystr)==trim(keyword)) then - ! - if (valstr(1:1) == '1' .or. valstr(1:1) == 'y' .or. valstr(1:1) == 'Y' .or. valstr(1:1) == 't' .or. valstr(1:1) == 'T') then - value = .true. - else - value = .false. - endif - ! - exit - ! - endif - ! - enddo - ! - end subroutine - - subroutine read_line(line0, keystr, valstr) - ! - ! Reads line from input file, returns keyword and value strings - ! - character(*), intent(in) :: line0 - character(len=256) :: line - character(*), intent(out) :: keystr - character(*), intent(out) :: valstr - integer j, ilen, jn - ! - keystr = '' - valstr = '' - ! - ! Change tabs into spaces. - ! - call notabs(line0, line, ilen) - ! - ! Look for line ending character. Remove it if it exists. - ! - jn = index(line, '\r') - ! - if (jn > 0) then - ! - ! New line character detected (probably sfincs.inp with windows line endings, running in linux) - ! - line = line(1 : jn - 1) - ! - endif - ! - ! Remove leading and trailing spaces. - ! - line = trim(line) - ! - if (line(1:1) == '#' .or. line(1:1) == '!' .or. line(1:1) == '@') return - ! - ! Find "=" - ! - j = index(line, '=') - ! - if (j == 0) return - ! - keystr = trim(line(1:j-1)) - ! - valstr = trim(line(j+1:)) - ! - ! Remove comments - ! - jn = index(valstr, '#') - ! - if (jn > 0) then - ! - valstr = trim(valstr(1 : jn - 1)) - ! - endif - ! - valstr = adjustl(trim(valstr)) - ! - end subroutine - - - subroutine notabs(INSTR,OUTSTR,ILEN) - ! @(#) convert tabs in input to spaces in output while maintaining columns, assuming a tab is set every 8 characters - ! - ! USES: - ! It is often useful to expand tabs in input files to simplify further processing such as tokenizing an input line. - ! Some FORTRAN compilers hate tabs in input files; some printers; some editors will have problems with tabs - ! AUTHOR: - ! John S. Urban - ! - ! SEE ALSO: - ! GNU/Unix commands expand(1) and unexpand(1) - ! - use ISO_FORTRAN_ENV, only : ERROR_UNIT ! get unit for standard error. if not supported yet, define ERROR_UNIT for your system (typically 0) - character(len=*),intent(in) :: INSTR ! input line to scan for tab characters - character(len=*),intent(out) :: OUTSTR ! tab-expanded version of INSTR produced - integer,intent(out) :: ILEN ! column position of last character put into output string - - integer,parameter :: TABSIZE=8 ! assume a tab stop is set every 8th column - character(len=1) :: c ! character read from stdin - integer :: ipos ! position in OUTSTR to put next character of INSTR - integer :: lenin ! length of input string trimmed of trailing spaces - integer :: lenout ! number of characters output string can hold - integer :: i10 ! counter that advances thru input string INSTR one character at a time - ! - IPOS=1 ! where to put next character in output string OUTSTR - lenin=len(INSTR) ! length of character variable INSTR - lenin=len_trim(INSTR(1:lenin)) ! length of INSTR trimmed of trailing spaces - lenout=len(OUTSTR) ! number of characters output string OUTSTR can hold - OUTSTR=" " ! this SHOULD blank-fill string, a buggy machine required a loop to set all characters - ! - do i10=1,lenin ! look through input string one character at a time - c=INSTR(i10:i10) - if(ichar(c) == 9)then ! test if character is a tab (ADE (ASCII Decimal Equivalent) of tab character is 9) - IPOS = IPOS + (TABSIZE - (mod(IPOS-1,TABSIZE))) - else ! c is anything else other than a tab insert it in output string - if(IPOS > lenout)then - write(ERROR_UNIT,*)"*notabs* output string overflow" - exit - else - OUTSTR(IPOS:IPOS)=c - IPOS=IPOS+1 - endif - endif - enddo - ! - ILEN=len_trim(OUTSTR(:IPOS)) ! trim trailing spaces - return - ! - end subroutine notabs - end module diff --git a/source/src/sfincs_spiderweb.f90 b/source/src/sfincs_spiderweb.f90 index 2c0abb618..7e9e6818b 100644 --- a/source/src/sfincs_spiderweb.f90 +++ b/source/src/sfincs_spiderweb.f90 @@ -1,6 +1,7 @@ module sfincs_spiderweb use sfincs_log + use sfincs_read contains @@ -404,61 +405,6 @@ subroutine read_amuv_dimensions(filename,nt,nrows,ncols,x_llcorner,y_llcorner,dx ! end subroutine - - subroutine read_real_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - real*4, intent(out) :: value - real*4, intent(in) :: default - integer j,stat - ! - value = default - rewind(fileid) - do while(.true.) - read(fileid,'(a)',iostat = stat)line - if (stat<0) exit - j=index(line,'=') - keystr = trim(line(1:j-1)) - if (trim(keystr)==trim(keyword)) then - valstr = trim(line(j+1:256)) - read(valstr,*)value - exit - endif - enddo - ! - end subroutine - - - subroutine read_int_input(fileid,keyword,value,default) - ! - character(*), intent(in) :: keyword - character(len=256) :: keystr - character(len=256) :: valstr - character(len=256) :: line - integer, intent(in) :: fileid - integer, intent(out) :: value - integer, intent(in) :: default - integer j,stat - ! - value = default - rewind(fileid) - do while(.true.) - read(fileid,'(a)',iostat = stat)line - if (stat<0) exit - j=index(line,'=') - keystr = trim(line(1:j-1)) - if (trim(keystr)==trim(keyword)) then - valstr = trim(line(j+1:256)) - read(valstr,*)value - exit - endif - enddo - ! - end subroutine subroutine compute_time_in_seconds(line,trefstr,dtsec) diff --git a/source/src/sfincs_vegetation.f90 b/source/src/sfincs_vegetation.f90 new file mode 100644 index 000000000..f6d2d8681 --- /dev/null +++ b/source/src/sfincs_vegetation.f90 @@ -0,0 +1,215 @@ +module sfincs_vegetation + + use sfincs_log + use sfincs_error + +contains + + subroutine initialize_vegetation() + ! + use sfincs_data + use sfincs_ncinput + ! + implicit none + ! + integer :: nm, nmu, ip, iveg, k + real*4 :: dh_veg, h_k, section_bottom, section_top + ! + logical :: ok + ! + character*256 :: varname + ! + if (use_quadtree .eqv. .false.) then + ! + call stop_sfincs('Error ! Netcdf vegetation input format can only be specified for quadtree mesh model !', 1) + ! + endif + ! + if (store_vegetation) then !either SFINCS and/or SnapWave needs veggie input + ! + write(logstr,'(a,a)')'Info : reading vegetation file ',trim(veggiefile) + call write_log(logstr, 0) + ! + ok = check_file_exists(veggiefile, 'Vegetation file', .true.) + ! + ! Get dimension of vertical sections + ! + ! Call the generic quadtree nc file reader function + varname = 'nsec' + !varname = 'vegetation_vertical_segments' ! TODO: change naming netcdf file into this + call read_netcdf_quadtree_get_dimension(veggiefile, varname, vegetation_vertical_segments) !ncfile, varname, varout) + ! + if (vegetation_vertical_segments > 4) then + ! + call stop_sfincs('Error ! Maximum allowed vertical sections in vegetationfile is 4 !', 1) + elseif(vegetation_vertical_segments == 0) then + ! + call stop_sfincs('Error ! Prescribed vertical sections in vegetationfile is 0 !', 1) + ! + endif + ! + ! allocate variables + allocate(vegetation_cd(np, vegetation_vertical_segments)) + allocate(vegetation_stems_height(np, vegetation_vertical_segments)) !=vegetation_ah + allocate(vegetation_stems_width(np, vegetation_vertical_segments)) !=vegetation_bstems + allocate(vegetation_stems_density(np, vegetation_vertical_segments)) !=vegetation_Nstems + ! + vegetation_cd = 0.0 + vegetation_stems_height = 0.0 + vegetation_stems_width = 0.0 + vegetation_stems_density = 0.0 + ! + ! Call the generic quadtree nc file reader function + varname = 'snapwave_veg_Cd' + !varname = 'vegegation_cd' ! TODO: change naming netcdf file into this + call read_netcdf_quadtree_to_sfincs(veggiefile, varname, vegetation_cd) !ncfile, varname, varout) + ! + ! Call the generic quadtree nc file reader function + varname = 'snapwave_veg_ah' + !varname = 'vegetation_stems_height' ! TODO: change naming netcdf file into this + call read_netcdf_quadtree_to_sfincs(veggiefile, varname, vegetation_stems_height) !ncfile, varname, varout) + ! + ! Call the generic quadtree nc file reader function + varname = 'snapwave_veg_bstems' + !varname = 'vegetation_stems_width' ! TODO: change naming netcdf file into this + call read_netcdf_quadtree_to_sfincs(veggiefile, varname, vegetation_stems_width) !ncfile, varname, varout) + ! + ! Call the generic quadtree nc file reader function + varname = 'snapwave_veg_Nstems' + !varname = 'vegetation_stems_density' ! TODO: change naming netcdf file into this + call read_netcdf_quadtree_to_sfincs(veggiefile, varname, vegetation_stems_density) !ncfile, varname, varout) + ! + endif + ! + ! For SFINCS precalculate cd * bstems * Nstems for each vertical section, as well as the vegetation height on uv points + ! + ! TODO - do this now as pre-processing table + ! + if (vegetation) then + ! + allocate(vegetation_stems_cd_width_density_uv(npuv, vegetation_vertical_segments)) ! product of cd, width and density on uv points + ! + allocate(vegetation_stems_height_uv(npuv, vegetation_vertical_segments)) ! vegetation height on uv points + ! + vegetation_stems_cd_width_density_uv = 0.0 + vegetation_stems_height_uv = 0.0 + ! + ! Interpolate vegetation properties from z-points to uv-points + ! + do ip = 1, npuv + ! + nm = uv_index_z_nm(ip) + nmu = uv_index_z_nmu(ip) + ! + do iveg = 1, vegetation_vertical_segments + ! + vegetation_stems_height_uv(ip,iveg) = 0.5*(vegetation_stems_height(nm,iveg)+vegetation_stems_height(nmu,iveg)) + ! + vegetation_stems_cd_width_density_uv(ip,iveg) = 0.5 * (0.5 * (vegetation_cd(nm,iveg) + vegetation_cd(nmu,iveg))) * (0.5 * (vegetation_stems_width(nm,iveg) + vegetation_stems_width(nmu,iveg))) * (0.5 * (vegetation_stems_density(nm,iveg) + vegetation_stems_density(nmu,iveg))) / rhow + ! + ! vegetation_stems_cd_width_density = 0.5 * cd * stems_width * stems_density / rhow, so everything that is precalculatable + ! + enddo + ! + enddo + ! + ! Pre-compute lookup table: cumulative sum of cd*width*density at vegetation_nlookup equidistant depth levels + ! Sections are stacked from the bed (consistent with swvegatt in snapwave_solver.f90): + ! vegetation_cd_sum_table(ip, k) = sum_iveg( cd_wd(ip,iveg) * max(0, min(section_top_iveg, h_k) - section_bottom_iveg) ) + ! In compute_fluxes: fvm = table_lookup(ip, hu) * uv0 * |uv0| (no inner do-loop needed) + ! + allocate(vegetation_lookup_hmin_uv(npuv)) + allocate(vegetation_lookup_hmax_uv(npuv)) + allocate(vegetation_lookup_dh_uv(npuv)) + allocate(vegetation_cd_sum_table(npuv, 0:vegetation_nlookup)) + allocate(vegetation_cd_slope_table(npuv, 0:vegetation_nlookup-1)) + ! + vegetation_lookup_hmin_uv = 0.0 + vegetation_lookup_hmax_uv = 0.0 + vegetation_lookup_dh_uv = 0.0 + vegetation_cd_sum_table = 0.0 + vegetation_cd_slope_table = 0.0 + ! + !$omp parallel do private( ip, k, iveg, dh_veg, h_k, section_bottom, section_top ) schedule( static ) + do ip = 1, npuv + ! + ! Sections are stacked from the bed upward (consistent with swvegatt in snapwave_solver): + ! section_bottom(iveg) = sum of all previous section heights + ! section_top(iveg) = section_bottom + vegetation_stems_height_uv(ip, iveg) + ! hmax = total vegetation height = sum of all section thicknesses + ! hmin = height of the bottom of the lowest section = 0 (all sections start at the bed) + ! + vegetation_lookup_hmin_uv(ip) = 0.0 + vegetation_lookup_hmax_uv(ip) = sum(vegetation_stems_height_uv(ip,:)) + ! + if (vegetation_lookup_hmax_uv(ip) > 0.0) then + dh_veg = vegetation_lookup_hmax_uv(ip) / real(vegetation_nlookup) + vegetation_lookup_dh_uv(ip) = dh_veg + do k = 1, vegetation_nlookup + h_k = k * dh_veg + section_bottom = 0.0 + do iveg = 1, vegetation_vertical_segments + section_top = section_bottom + vegetation_stems_height_uv(ip, iveg) + vegetation_cd_sum_table(ip, k) = vegetation_cd_sum_table(ip, k) + vegetation_stems_cd_width_density_uv(ip, iveg) * max(0.0, min(section_top, h_k) - section_bottom) + section_bottom = section_top + enddo + enddo + endif + ! + enddo + !$omp end parallel do + ! + ! Pre-compute slope between consecutive table entries to avoid the subtraction in compute_fluxes + ! + !$omp parallel do private( ip, k ) schedule( static ) + do ip = 1, npuv + do k = 0, vegetation_nlookup - 1 + vegetation_cd_slope_table(ip, k) = vegetation_cd_sum_table(ip, k+1) - vegetation_cd_sum_table(ip, k) + enddo + enddo + !$omp end parallel do + ! + endif + ! + ! ----------------------------------------------------------------------- + ! Summary: vegetation drag pre-computation + ! ----------------------------------------------------------------------- + ! + ! INPUT (read from vegetation NetCDF file, on z-points): + ! vegetation_cd : bulk drag coefficient [-] (np, nsec) + ! vegetation_stems_height : section thickness, stacked bed upward [m] (np, nsec) + ! vegetation_stems_width : stem diameter [m] (np, nsec) + ! vegetation_stems_density : stem density [m-2] (np, nsec) + ! + ! STEP 1 - interpolate to uv-points and pre-multiply constants: + ! vegetation_stems_height_uv(ip,iveg) = average of nm and nmu cell heights + ! vegetation_stems_cd_width_density_uv(ip,iveg) = 0.5 * cd * b * N / rho_w + ! (factor 0.5 and rho_w division folded in once; never recomputed at runtime) + ! + ! STEP 2 - build lookup table (vegetation_nlookup equidistant depth bins): + ! hmax(ip) = sum of all section heights (total vegetation height) + ! dh(ip) = hmax / vegetation_nlookup (bin width, stored for runtime use) + ! For each bin k (depth level h_k = k * dh): + ! table(ip,k) = sum_iveg [ cd_wd(ip,iveg) * max(0, min(sec_top, h_k) - sec_bot) ] + ! This is the cumulative drag integral up to depth h_k. + ! Sections are stacked from the bed (sec_bot = sum of previous section heights), + ! consistent with the swvegatt convention in snapwave_solver.f90. + ! + ! STEP 3 - pre-compute slope table (avoids subtraction at runtime): + ! slope_table(ip,k) = table(ip,k+1) - table(ip,k) + ! + ! RUNTIME USE (compute_fluxes in sfincs_momentum.f90): + ! Given water depth hu at a uv-point: + ! frac = min(hu, hmax(ip)) / dh(ip) ! fractional bin index + ! ik = floor(frac) ! integer bin + ! frac = frac - ik ! remainder for interpolation + ! cd_eff = table(ip,ik) + frac * slope(ip,ik) ! O(1) lookup, no section loop + ! Vegetation drag force (explicit, added to frc): + ! F_veg = -phi * cd_eff * u0 * |u0| + ! Flux update (Manning friction handled implicitly in denominator): + ! q_new = (q_old + (F_ext + F_veg) * dt) / (1 + g*n^2*|q|/hu^(7/3) * dt) + ! ----------------------------------------------------------------------- + ! + end subroutine + +end module \ No newline at end of file diff --git a/source/src/sfincs_wavemaker.f90 b/source/src/sfincs_wavemaker.f90 index 1bf5a6e04..28951f3ee 100644 --- a/source/src/sfincs_wavemaker.f90 +++ b/source/src/sfincs_wavemaker.f90 @@ -51,7 +51,7 @@ subroutine initialize_wavemakers() real*4, dimension(:), allocatable :: wavemaker_xfp real*4, dimension(:), allocatable :: wavemaker_yfp ! - logical :: iok, ok + logical :: iok, ok, refinement_warning ! integer ib1, ib2, ib, ic, nmb, nrwvm ! @@ -304,7 +304,7 @@ subroutine initialize_wavemakers() ! if (nmu>0) then ! - iz = uv_index_z_nm(ip) + iz = uv_index_z_nm(nmu) ! if (indwm(iz)==0) then ! @@ -759,6 +759,8 @@ subroutine initialize_wavemakers() write(logstr,*)'Setting wave makers ...' call write_log(logstr, 0) ! + refinement_warning = .false. ! set to true if we find a wavemaker point that has refinemed neighbor + ! do ip = 1, np ! if (kcs(ip)==4) then @@ -798,6 +800,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nmu(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -807,6 +811,8 @@ subroutine initialize_wavemakers() wavemaker_idir(iwm) = 1 wavemaker_angfac(iwm) = max(cos(phi(ip) - 0.0), 0.0) ! + wavemaker_nmu(nok) = iwm + ! endif ! endif @@ -842,6 +848,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nmu(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -851,6 +859,8 @@ subroutine initialize_wavemakers() wavemaker_idir(iwm) = 1 wavemaker_angfac(iwm) = max(sin(phi(ip) - 0.0), 0.0) ! + wavemaker_num(nok) = iwm + ! endif ! endif @@ -888,6 +898,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nm(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -897,6 +909,8 @@ subroutine initialize_wavemakers() wavemaker_idir(iwm) = -1 wavemaker_angfac(iwm) = max(cos(pi - phi(ip)), 0.0) ! + wavemaker_nmd(nok) = iwm + ! endif ! endif @@ -932,6 +946,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nmu(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -941,6 +957,8 @@ subroutine initialize_wavemakers() wavemaker_idir(iwm) = 1 wavemaker_angfac(iwm) = max(sin(phi(ip) - 0.0), 0.0) ! + wavemaker_num(nok) = iwm + ! endif ! endif @@ -978,6 +996,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nm(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -987,6 +1007,8 @@ subroutine initialize_wavemakers() wavemaker_idir(iwm) = -1 wavemaker_angfac(iwm) = max(cos(pi - phi(ip)), 0.0) ! + wavemaker_nmd(nok) = iwm + ! endif ! endif @@ -1007,7 +1029,7 @@ subroutine initialize_wavemakers() wavemaker_index_nmi(iwm) = iz wavemaker_index_nmb(iwm) = ip wavemaker_idir(iwm) = -1 - wavemaker_angfac(iwm) = max(sin(pi - phi(ip)), 0.0) + wavemaker_angfac(iwm) = max(-sin(phi(ip)), 0.0) ! wavemaker_ndm(nok) = iwm ! @@ -1022,6 +1044,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nm(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -1029,7 +1053,9 @@ subroutine initialize_wavemakers() wavemaker_index_nmi(iwm) = iz wavemaker_index_nmb(iwm) = ip wavemaker_idir(iwm) = -1 - wavemaker_angfac(iwm) = max(sin(pi - phi(ip)), 0.0) + wavemaker_angfac(iwm) = max(-sin(phi(ip)), 0.0) + ! + wavemaker_ndm(nok) = iwm ! endif ! @@ -1067,6 +1093,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nmu(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -1076,6 +1104,8 @@ subroutine initialize_wavemakers() wavemaker_idir(iwm) = 1 wavemaker_angfac(iwm) = max(cos(phi(ip) - 0.0), 0.0) ! + wavemaker_nmu(nok) = iwm + ! endif ! endif @@ -1096,7 +1126,7 @@ subroutine initialize_wavemakers() wavemaker_index_nmi(iwm) = iz wavemaker_index_nmb(iwm) = ip wavemaker_idir(iwm) = -1 - wavemaker_angfac(iwm) = max(sin(pi - phi(ip)), 0.0) + wavemaker_angfac(iwm) = max(-sin(phi(ip)), 0.0) ! wavemaker_ndm(nok) = iwm ! @@ -1111,6 +1141,8 @@ subroutine initialize_wavemakers() iz = uv_index_z_nm(nmu) ! if (kcs(iz) == 1) then + ! + refinement_warning = .true. ! iwm = iwm + 1 ! @@ -1118,7 +1150,9 @@ subroutine initialize_wavemakers() wavemaker_index_nmi(iwm) = iz wavemaker_index_nmb(iwm) = ip wavemaker_idir(iwm) = -1 - wavemaker_angfac(iwm) = max(sin(pi - phi(ip)), 0.0) + wavemaker_angfac(iwm) = max(-sin(phi(ip)), 0.0) + ! + wavemaker_ndm(nok) = iwm ! endif ! @@ -1128,6 +1162,15 @@ subroutine initialize_wavemakers() endif enddo ! + ! Give warning if we found a wavemaker point that has refined neighbor + ! + if (refinement_warning) then + ! + write(logstr,'(a)')' WARNING! Found wavemaker point along quadtree refinement boundary, this is not recommended! The simulation will continue.' + call write_log(logstr, 1) + ! + endif + ! ! Set flags for kcuv points ! do iwm = 1, wavemaker_nr_uv_points @@ -1437,10 +1480,6 @@ subroutine update_wavemaker_fluxes(t, dt, tloop) tp_inc = 10.0 ! Later make it possible to also specify Tp_inc in time series forcing, but for now just add a fixed value (that is not used) ! else - ! - ! Use mean peak period from SnapWave boundary conditions - ! - 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. ! @@ -1465,13 +1504,28 @@ subroutine update_wavemaker_fluxes(t, dt, tloop) ! ! ! tp_ig = snapwave_tpmean * max(1.86 * betas**-0.43 * wave_steepness**0.07, 5.0) ! ! + else + ! + ! Use mean peak period from SnapWave boundary conditions + ! + 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) + ! + if (tp_ig < 10.0) then + ! These warnings should not occur here + write(logstr,*)'DEBUG SFINCS_SnapWave - incoming tp for IG wave at wavemaker might be unrealistically small! value: ',tp_ig + call write_log(logstr, 0) + elseif (tp_ig > 250.0) then + write(logstr,*)'DEBUG SFINCS_SnapWave - incoming tp for IG wave at wavemaker might be unrealistically large! value: ',tp_ig + call write_log(logstr, 0) + endif + ! endif ! tp_inc = max(snapwave_tpmean, wavemaker_tpmin) ! tp_ig = max(tp_ig, wavemaker_tpmin) ! - endif + endif ! ! Now determine zwav_ig and zwav_inc based on spectrum or monochromatic signal. ! Time series of zwav_ig and zwav_inc will be used to modulate water level at wave maker points. @@ -1511,7 +1565,7 @@ subroutine update_wavemaker_fluxes(t, dt, tloop) ! fm_inc = 1.0 / tp_inc ! Wave period ! - do ifreq = 1, wavemaker_nfreqs_ig + do ifreq = 1, wavemaker_nfreqs_inc ! wavemaker_phi_inc(ifreq) = modulo(wavemaker_phi_inc(ifreq) + wavemaker_dphi_inc(ifreq) * dt, 2 * pi) wavemaker_cost_inc(ifreq) = cos(2 * pi * t * wavemaker_freq_inc(ifreq) + wavemaker_phi_inc(ifreq)) @@ -1613,6 +1667,11 @@ subroutine update_wavemaker_fluxes(t, dt, tloop) ! zsnmb = zs0nmb + min(zinc + zig, wavemaker_gammax * dwvm) ! total water level in wave maker (i.e. mean water level plus wave) ! + if (( zinc + zig) > wavemaker_gammax * dwvm) then + write(*,*)'WARNING! Incident wave height at wave maker exceeds maximum allowed value based on local water depth! Value: ', zinc + zig, ' Max allowed: ', wavemaker_gammax * dwvm + endif + + ! endif ! if (subgrid) then diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90 index 1ab955526..2925314b7 100644 --- a/source/src/snapwave/snapwave_data.f90 +++ b/source/src/snapwave/snapwave_data.f90 @@ -49,8 +49,8 @@ module snapwave_data real*4, dimension(:,:), allocatable :: ctheta360 ! refraction speed, per grid point and direction ! real*4, dimension(:), allocatable :: xn,yn,zn ! coordinates of nodes of unstructured grid real*4, dimension(:), allocatable :: dzdx,dzdy ! bed slopes at nodes of unstructured grid - integer, dimension(:,:), allocatable :: face_nodes ! node numbers connected to each cell - integer, dimension(:,:), allocatable :: edge_nodes ! node numbers connected to each edge + integer, dimension(:,:), allocatable :: face_nodes ! node numbers connected to each cell + integer, dimension(:,:), allocatable :: edge_nodes ! node numbers connected to each edge real*4, dimension(:), allocatable :: bndindx real*4, dimension(:), allocatable :: tau ! real*4, dimension(:,:), allocatable :: Fluxtab @@ -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 @@ -73,6 +74,8 @@ module snapwave_data real*4, dimension(:), allocatable :: beta real*4, dimension(:), allocatable :: srcig real*4, dimension(:), allocatable :: alphaig + real*4, dimension(:), allocatable :: qb + real*4, dimension(:), allocatable :: gam ! integer*4, dimension(:), allocatable :: index_snapwave_in_quadtree integer*4, dimension(:), allocatable :: index_quadtree_in_snapwave @@ -172,10 +175,9 @@ module snapwave_data real*4 :: fwcutoff ! depth below which to apply space-varying fw real*4 :: alpha,gamma ! coefficients in Baldock wave breaking dissipation model real*4 :: gammax ! max wave height/water depth ratio - integer :: baldock_opt ! option of Baldock wave breaking dissipation model (opt=1 is without gamma&depth, else is including) + !integer :: baldock_opt ! option of Baldock wave breaking dissipation model (opt=1 is without gamma&depth, else is including) real*4 :: baldock_ratio ! option controlling from what depth wave breaking should take place: (Hk>baldock_ratio*Hmx(k)), default baldock_ratio=0.2 - ! TODO - TL: bring back baldock_ratio? - + integer :: baldock_exponent! Exponent for multiplying the Baldock dissipation with a factor 'f = (Hloc / Hmax)**iexp' to enhance breaking when H > Hmax, with iexp = 0 (default, means unused), 1 or 2 real*4 :: hmin ! minimum water depth character*256 :: gridfile ! name of gridfile (Delft3D .grd format) integer :: sferic ! sferical (1) or cartesian (0) grid @@ -216,6 +218,8 @@ module snapwave_data real*4 :: rghlevland ! Elevation separation as in SFINCS for simple elevation varying roughness real*4 :: fwratio ! Above 'rghlevland' elevation of zb, the friction for incident waves is multiplied with value 'fwratio' real*4 :: fwigratio ! Above 'rghlevland' elevation of zb, the friction for IG waves is multiplied with value 'fwratio' + real*4 :: relax_factor_DoverA ! underrelaxation factor for DoverA (set to 1.0 to disable) + real*4 :: relax_factor_DoverE ! underrelaxation factor for DoverE (set to 1.0 to disable) ! character*3 :: outputformat integer :: ja_save_each_iter ! logical to save output after each iteration or not @@ -259,6 +263,7 @@ module snapwave_data ! integer :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) real*4 :: alpha_ig,gamma_ig ! coefficients in Baldock wave breaking dissipation model for IG waves + real*4 :: gamma_fac_br ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point 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 :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0 real*4 :: eeinc2ig ! ratio of incident wave energy as first estimate of IG wave energy at boundary @@ -270,6 +275,7 @@ module snapwave_data ! ... or just a priori based on effectively incident wave energy from previous timestep only integer :: herbers_opt ! 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) integer :: tpig_opt ! IG wave period option based on Herbers calculated spectrum, only used if herbers_opt = 1. Options are: 1=Tm01 (default), 2=Tpsmooth, 3=Tp, 4=Tm-1,0 + real*4 :: steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5 ! ! Switches logical :: igwaves ! switch whether include IG or not @@ -282,6 +288,7 @@ module snapwave_data ! logical :: restart logical :: coupled_to_sfincs + logical :: storesnapwavegrid ! integer :: nr_quadtree_points ! diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90 index 7c1992f2f..994270184 100644 --- a/source/src/snapwave/snapwave_domain.f90 +++ b/source/src/snapwave/snapwave_domain.f90 @@ -8,8 +8,10 @@ subroutine initialize_snapwave_domain() ! use snapwave_data use snapwave_boundaries + use snapwave_ncoutput use interp use sfincs_error + use quadtree ! ! Local input variables ! @@ -21,6 +23,7 @@ subroutine initialize_snapwave_domain() integer*4 :: idummy character*2 :: ext logical :: generate_upw, exists + character(len=256) :: snapwave_ncfname ! real*8 :: xmn, ymn ! ! First set some constants @@ -83,17 +86,16 @@ subroutine initialize_snapwave_domain() if (face_nodes(4,k)==0) face_nodes(4,k) = -999 enddo ! - ! Done with the mesh + ! write mesh to file + if (storesnapwavegrid) then + ! + snapwave_ncfname = 'snapwavegrid.nc' + ! + call write_snapwave_mesh(snapwave_ncfname, sferic == 1) + ! + endif ! - ! keep on also if ja_vegetation==0, so array Dveg is initialized with zeroes - !if (ja_vegetation==1) then - ! call veggie_init() - !else - allocate(veg_Cd(no_nodes, no_secveg)) - allocate(veg_ah(no_nodes, no_secveg)) - allocate(veg_bstems(no_nodes, no_secveg)) - allocate(veg_Nstems(no_nodes, no_secveg)) - !endif + ! Done with the mesh ! ntheta360 = nint(360./dtheta) ntheta = nint(sector/dtheta) @@ -136,6 +138,8 @@ subroutine initialize_snapwave_domain() allocate(beta(no_nodes)) allocate(srcig(no_nodes)) allocate(alphaig(no_nodes)) + allocate(qb(no_nodes)) + allocate(gam(no_nodes)) ! allocate(uorb(no_nodes)) allocate(ctheta(ntheta,no_nodes)) allocate(ctheta_ig(ntheta,no_nodes)) @@ -165,6 +169,7 @@ subroutine initialize_snapwave_domain() allocate(WsorA (ntheta,no_nodes)) allocate(SwE (no_nodes)) allocate(SwA (no_nodes)) + allocate(DoverE(no_nodes)) ! ! Spatially-uniform bottom friction coefficients ! @@ -172,10 +177,10 @@ subroutine initialize_snapwave_domain() fw_ig = fw0_ig ! do k=1,no_nodes - if (zb(k) > rghlevland) then - fw(k) = fw0 * fwratio - fw_ig(k) = fw0_ig * fwigratio - endif + if (zb(k) > rghlevland) then + fw(k) = fw0 * fwratio + fw_ig(k) = fw0_ig * fwigratio + endif enddo ! ! Initialization of reference tables @@ -200,13 +205,16 @@ subroutine initialize_snapwave_domain() prev360 = 0 H = 0.0 H_ig = 0.0 + Dw = 0.0 aa = 0.0 sig = 0.0 WsorE = 0.0 WsorA = 0.0 SwE = 0.0 SwA = 0.0 - windspreadfac = 0.0 + DoverE = 0.0 + windspreadfac = 0.0 + Hmx_ig = 0.0 ! generate_upw = .true. exists = .true. @@ -303,26 +311,12 @@ 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) + ! + write(logstr,*)'SnapWave: Neumann connected boundaries found ...' + call write_log(logstr, 0) + ! else ! neumannconnected = 0 @@ -796,136 +790,62 @@ 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 - ! - end subroutine neuboundaries_light - - -subroutine neuboundaries(x,y,no_nodes,xneu,yneu,n_neu,tol,neumannconnected) + real :: xgb, ygb, dst1, dst + integer :: k, ib1, ic ! - implicit none + ! Loop through all msk=3 cells and find their nearest msk=1 cell, save in 'neumannconnected' ! - integer, intent(in) :: no_nodes - integer, intent(in) :: n_neu - real*8, dimension(no_nodes), intent(in) :: x,y - real*8, dimension(n_neu), intent(in) :: xneu,yneu - real*4, intent(in) :: tol - integer, dimension(no_nodes), intent(out) :: neumannconnected - ! - integer :: ib,k,kmin, k2 - real*8 :: alpha, cosa,sina, distmin, x1,y1,x2,y2, xend - ! - neumannconnected=0 - do ib=1,n_neu-1 - if (xneu(ib).ne.-999.and.xneu(ib+1).ne.-999) then - alpha=atan2(yneu(ib+1)-yneu(ib),xneu(ib+1)-xneu(ib)) - cosa=cos(alpha) - sina=sin(alpha) - xend=(xneu(ib+1)-xneu(ib))*cosa+(yneu(ib+1)-yneu(ib))*sina - do k=1,no_nodes - x1= (x(k)-xneu(ib))*cosa+(y(k)-yneu(ib))*sina - y1=-(x(k)-xneu(ib))*sina+(y(k)-yneu(ib))*cosa - if (x1>=0.d0 .and. x1<=xend) then - if (abs(y1)0) then - neumannconnected(kmin)=k - write(logstr,*)kmin,k - call write_log(logstr, 0) - endif - endif - endif + 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 + ! + enddo ! -end subroutine neuboundaries + end subroutine neuboundaries_light + subroutine read_snapwave_sfincs_mesh() @@ -1125,6 +1045,9 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) ! use snapwave_data use quadtree + use sfincs_data, only: vegetation_cd, vegetation_stems_height, & + vegetation_stems_width, vegetation_stems_density, & + vegetation_vertical_segments ! ! Local input variables ! @@ -1143,13 +1066,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 @@ -1166,7 +1095,7 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) ! 4) Loop through all points and make cells for points where msk==1. ! The node indices in the cells will point to the indices of the entire quadtree. ! In a second temporary mask array msk_tmp2, determine which nodes are actually active (being part a cell) - ! 5) Set back snapwave_mask = 2&3 values of wave boudnary and neumann cells + ! 5) Set back snapwave_mask = 2&3 values of wave boundary and neumann cells ! 6) Count actual number of active nodes and cells, and allocate arrays ! 7) Set node data and re-map indices ! @@ -1174,12 +1103,12 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) ! ! Check if qtr file has already been loaded by other model (sfincs) ! -! if (load_quadtree) then -! ! -! write(*,*)'Reading SnapWave quadtree file ', trim(gridfile), ' ...' -! call quadtree_read_file(gridfile) -! ! -! endif + ! if (load_quadtree) then + ! ! + ! write(*,*)'Reading SnapWave quadtree file ', trim(gridfile), ' ...' + ! call quadtree_read_file(gridfile) + ! ! + ! endif ! allocate(index_snapwave_in_quadtree(quadtree_nr_points)) ! Needed for mapping to sfincs ! @@ -1190,6 +1119,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 @@ -1221,16 +1151,6 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) ! Count number of active points ! This is also the number of points in the dep file ! -! nra = 0 -! do ip = 1, quadtree_nr_points -! if (msk_tmp(ip)>0) then -! nra = nra + 1 -! endif -! enddo - ! -! allocate(zb_tmp2(nra)) ! Make (very) temporary array with bed level on all active quadtree points -! zb_tmp2 = -10.0 - ! if (depfile /= 'none') then ! write(logstr,*)'Reading SnapWave depth file ',trim(depfile),' ...' @@ -1243,19 +1163,11 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) ! ! Now loop through all quadtree points and set depth ! -! nra = 0 -! do ip = 1, quadtree_nr_points -! if (msk_tmp(ip)>0) then -! nra = nra + 1 -! zb_tmp(ip) = zb_tmp2(nra) -! endif -! enddo -! ! -! deallocate(zb_tmp2) - ! ! 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 ! @@ -1415,47 +1327,6 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) endif endif ! - if (mnu1==0) then - ! Didn't find it going to the right, try via above -! if nu==0 -! ! same level above -! if nu1>0 -! if buq.mu(nu1)==0 -! ! same level above right -! if buq.mu1(nu1)>0 -! ! and it exists -! mnu=0; -! mnu1=buq.mu1(nu1); -! end -! end -! end -! elseif mu==-1 -! ! coarser above -! if nu1>0 -! if buq.mu(nu1)==0 -! ! same level above right -! if buq.mu1(nu1)>0 -! ! and it exists -! mnu=-1; -! mnu1=buq.mu1(nu1); -! end -! end -! end -! else -! ! finer above -! if nu2>0 -! if buq.mu(nu2)==0 -! ! same level above right -! if buq.mu1(nu2)>0 -! ! and it exists -! mnu=1; -! mnu1=buq.mu1(nu2); -! end -! end -! end -! end - endif - ! ! Okay, found all the neighbors! ! ! Now let's see what sort of cells we need @@ -1464,7 +1335,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 @@ -1679,14 +1549,6 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) msk_tmp2(nu1) = 1 endif ! -!% elseif (mu==-1 .and. nu==0 .and. mnu==0 .and. odd(buq.n(ip))) -!% % Type 9 -!% if mu1>0 .and. nu1>0 -!% nfaces=nfaces+1; -!% faces(1, nfaces) = ip; -!% faces(2, nfaces) = mu1; -!% faces(3, nfaces) = nu1; -!% end elseif (mu==-1 .and. nu==-1 .and. mnu==-1) then ! ! Type 10 @@ -2033,7 +1895,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 .and. quadtree_nu1(md1) == 0) then + ! + nfaces = nfaces + 1 + faces(1, nfaces) = md1 + faces(2, nfaces) = ip + faces(3, nfaces) = nu1 + ! + 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 .and. quadtree_nd1(md1) == 0) 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 .and. quadtree_nu1(mu1) == 0) 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 .and. quadtree_nd1(mu1) == 0) then + ! + nfaces = nfaces + 1 + faces(1, nfaces) = ip + faces(2, nfaces) = nd1 + faces(3, nfaces) = mu1 + ! + endif + ! + endif ! enddo ! @@ -2098,7 +2033,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) @@ -2110,11 +2044,11 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) ! enddo ! - ! Loop through cells to re-maps the face nodes + ! STEP 8 - Loop through cells to re-maps the face nodes ! 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 @@ -2122,6 +2056,44 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree) enddo enddo ! + ! STEP 9 - if vegetation, re-map veggie input from quadtree netcdf vegetationfile + ! Set 'no_secveg' from sfincs_vegetation.f90 for use in snapwave_data + ! + no_secveg = vegetation_vertical_segments + ! + allocate(veg_Cd(no_nodes, no_secveg)) + allocate(veg_ah(no_nodes, no_secveg)) + allocate(veg_bstems(no_nodes, no_secveg)) + allocate(veg_Nstems(no_nodes, no_secveg)) + ! + veg_Cd = 0.0 + veg_ah = 0.0 + veg_bstems = 0.0 + veg_Nstems = 0.0 + ! + if (vegetation) then + ! copy from the quadtree snapwave_veg + nac = 0 + ! + do ip = 1, quadtree_nr_points + ! + if (msk_tmp2(ip)>0) then + ! + nac = nac + 1 + ! + ! Set node values for all points in the vertical + do iq = 1, no_secveg + veg_Cd(nac,iq) = vegetation_cd(ip,iq) + veg_ah(nac,iq) = vegetation_stems_height(ip,iq) + veg_bstems(nac,iq) = vegetation_stems_width(ip,iq) + veg_Nstems(nac,iq) = vegetation_stems_density(ip,iq) + enddo + ! + endif + enddo + ! + endif + ! end subroutine end module diff --git a/source/src/snapwave/snapwave_infragravity.f90 b/source/src/snapwave/snapwave_infragravity.f90 index 178f32843..741ce4515 100644 --- a/source/src/snapwave/snapwave_infragravity.f90 +++ b/source/src/snapwave/snapwave_infragravity.f90 @@ -43,22 +43,26 @@ subroutine determine_ig_bc(x_bwv, y_bwv, hsinc, tpinc, ds, jonswapgam, depth, Ti scoeff = (2/ds**2) - 1 ! ! 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 - call write_log(logstr, 1) - ! - write(logstr,*)'This 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 - ! - endif + ! Loosely based on 3 step calculation in waveparams.F90 of XBeach (build_jonswap, build_etdir, build_boundw), here all in 1 subroutine compute_herbers + ! + if (hsinc / depth > 0.5) then + ! + write(logstr, *)'ERROR SnapWave - Hs over depth at boundary input point ', x_bwv, ',', y_bwv,' is larger then 0.5: ', hsinc / 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 2.0 * hsinc meters for stability, simulation will continue.' + call write_log(logstr, 0) + ! + depth = 2.0 * hsinc + ! + elseif (depth > 200.0) then + ! + ! Limit depth to 200 m. Larger depth can result in NaNs. @Tim, why? + ! + depth = 200.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] ! @@ -69,9 +73,7 @@ subroutine determine_ig_bc(x_bwv, y_bwv, hsinc, tpinc, ds, jonswapgam, depth, Ti call write_log(logstr, 1) hsig = max(hsig, 0.0) ! - endif - ! - if (hsig > 3.0) then + elseif (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) diff --git a/source/src/snapwave/snapwave_ncoutput.F90 b/source/src/snapwave/snapwave_ncoutput.F90 new file mode 100644 index 000000000..bd2710fb0 --- /dev/null +++ b/source/src/snapwave/snapwave_ncoutput.F90 @@ -0,0 +1,140 @@ +#define NF90(nf90call) call handle_err(nf90call,__FILE__,__LINE__) +module snapwave_ncoutput + ! + use sfincs_log + use netcdf + ! + implicit none + ! + contains + ! + subroutine write_snapwave_mesh(fname, crsgeo) + ! + use snapwave_data + ! + implicit none + ! + character(len=256), intent(in) :: fname + logical, intent(in) :: crsgeo + ! + integer :: ncid + integer :: nmesh2d_node_dimid, nmesh2d_face_dimid, max_nmesh2d_face_nodes_dimid + integer :: mesh2d_varid + integer :: mesh2d_node_x_varid, mesh2d_node_y_varid, crs_varid + integer :: mesh2d_face_nodes_varid + integer :: zb_varid + ! + integer, parameter :: nc_deflate_level = 2 + real*4, parameter :: FILL_VALUE = -99999.0 + ! + ! dimensions + NF90(nf90_create(trim(fname), ior(NF90_CLOBBER, NF90_NETCDF4), ncid)) + NF90(nf90_def_dim(ncid, 'nmesh2d_node', no_nodes, nmesh2d_node_dimid)) + NF90(nf90_def_dim(ncid, 'nmesh2d_face', no_faces, nmesh2d_face_dimid)) + NF90(nf90_def_dim(ncid, 'max_nmesh2d_face_nodes', 4, max_nmesh2d_face_nodes_dimid)) + ! + ! global attributes + NF90(nf90_put_att(ncid,nf90_global, "Conventions", "Conventions = 'CF-1.8 UGRID-1.0 Deltares-0.10'")) + NF90(nf90_put_att(ncid,nf90_global, "Build-Revision-Date-Netcdf-library", trim(nf90_inq_libvers()))) + NF90(nf90_put_att(ncid,nf90_global, "Producer", "SFINCS model: Super-Fast INundation of CoastS")) + NF90(nf90_put_att(ncid,nf90_global, "Build-Revision", trim(build_revision))) + NF90(nf90_put_att(ncid,nf90_global, "Build-Date", trim(build_date))) + NF90(nf90_put_att(ncid,nf90_global, "title", "Snapwave grid")) + ! + ! mesh topology + NF90(nf90_def_var(ncid, 'mesh2d', NF90_INT, mesh2d_varid)) + NF90(nf90_put_att(ncid, mesh2d_varid, 'cf_role', 'mesh_topology')) + NF90(nf90_put_att(ncid, mesh2d_varid, 'long_name', 'Topology data of 2D network')) + NF90(nf90_put_att(ncid, mesh2d_varid, 'topology_dimension', 2)) + NF90(nf90_put_att(ncid, mesh2d_varid, 'node_coordinates', 'mesh2d_node_x mesh2d_node_y')) + NF90(nf90_put_att(ncid, mesh2d_varid, 'node_dimension', 'nmesh2d_node')) + NF90(nf90_put_att(ncid, mesh2d_varid, 'max_face_nodes_dimension', 'max_nmesh2d_face_nodes')) + NF90(nf90_put_att(ncid, mesh2d_varid, 'face_node_connectivity', 'mesh2d_face_nodes')) + NF90(nf90_put_att(ncid, mesh2d_varid, 'face_dimension', 'nmesh2d_face')) + ! + if (crsgeo) then + ! + NF90(nf90_def_var(ncid, 'mesh2d_node_x', NF90_FLOAT, (/nmesh2d_node_dimid/), mesh2d_node_x_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(ncid, mesh2d_node_x_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'units', 'degrees')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'standard_name', 'longitude')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'long_name', 'longitude')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'mesh', 'mesh2d')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'location', 'node')) + ! + NF90(nf90_def_var(ncid, 'mesh2d_node_y', NF90_FLOAT, (/nmesh2d_node_dimid/), mesh2d_node_y_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(ncid, mesh2d_node_y_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'units', 'degrees')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'standard_name', 'latitude')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'long_name', 'latitude')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'mesh', 'mesh2d')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'location', 'node')) + ! + else + ! + NF90(nf90_def_var(ncid, 'mesh2d_node_x', NF90_DOUBLE, (/nmesh2d_node_dimid/), mesh2d_node_x_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(ncid, mesh2d_node_x_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'units', 'm')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'standard_name', 'projection_x_coordinate')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'long_name', 'x-coordinate of mesh nodes')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'mesh', 'mesh2d')) + NF90(nf90_put_att(ncid, mesh2d_node_x_varid, 'location', 'node')) + ! + NF90(nf90_def_var(ncid, 'mesh2d_node_y', NF90_DOUBLE, (/nmesh2d_node_dimid/), mesh2d_node_y_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(ncid, mesh2d_node_y_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'units', 'm')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'standard_name', 'projection_y_coordinate')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'long_name', 'y-coordinate of mesh nodes')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'mesh', 'mesh2d')) + NF90(nf90_put_att(ncid, mesh2d_node_y_varid, 'location', 'node')) + ! + endif + ! + NF90(nf90_def_var(ncid, 'mesh2d_face_nodes', NF90_INT, (/max_nmesh2d_face_nodes_dimid, nmesh2d_face_dimid/), mesh2d_face_nodes_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(ncid, mesh2d_face_nodes_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(ncid, mesh2d_face_nodes_varid, 'cf_role', 'face_node_connectivity')) + NF90(nf90_put_att(ncid, mesh2d_face_nodes_varid, 'mesh', 'mesh2d')) + NF90(nf90_put_att(ncid, mesh2d_face_nodes_varid, 'location', 'face')) + NF90(nf90_put_att(ncid, mesh2d_face_nodes_varid, 'long_name', 'Mapping from every face to its corner nodes (counterclockwise)')) + NF90(nf90_put_att(ncid, mesh2d_face_nodes_varid, 'start_index', 1)) + NF90(nf90_put_att(ncid, mesh2d_face_nodes_varid, '_FillValue', -999)) + ! + NF90(nf90_def_var(ncid, 'crs', NF90_INT, crs_varid)) ! For EPSG code + NF90(nf90_put_att(ncid, crs_varid, 'EPSG', '-')) + ! + NF90(nf90_def_var(ncid, 'mesh2d_node_z', NF90_FLOAT, (/nmesh2d_node_dimid/), zb_varid)) ! bed level in cell centre + NF90(nf90_def_var_deflate(ncid, zb_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(ncid, zb_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(ncid, zb_varid, 'units', 'm')) + NF90(nf90_put_att(ncid, zb_varid, 'standard_name', 'altitude')) + NF90(nf90_put_att(ncid, zb_varid, 'long_name', 'bed_level_above_reference_level')) + ! + NF90(nf90_enddef(ncid)) + ! + ! put variables + NF90(nf90_put_var(ncid, mesh2d_node_x_varid, x)) ! write node x + NF90(nf90_put_var(ncid, mesh2d_node_y_varid, y)) ! write node y + NF90(nf90_put_var(ncid, mesh2d_face_nodes_varid, face_nodes)) + NF90(nf90_put_var(ncid, zb_varid, zb)) + + ! close file + NF90(nf90_close(ncid)) + + end subroutine write_snapwave_mesh + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine handle_err(status,file,line) + ! + integer, intent ( in) :: status + character(*), intent(in) :: file + integer, intent ( in) :: line + integer :: status2 + + if(status /= nf90_noerr) then + ! !UNIT=6 for stdout and UNIT=0 for stderr. + write(0,'("NETCDF ERROR: ",a,i6,":",a)') file,line,trim(nf90_strerror(status)) + end if + end subroutine handle_err + ! + end module \ No newline at end of file diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 53cc31d77..16af218f0 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -1,9 +1,10 @@ -module snapwave_solver +module snapwave_solver - use sfincs_log + use sfincs_log - implicit none - contains + 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) ! @@ -35,85 +33,120 @@ subroutine compute_wave_field() ! ! Set energies to 0.0; note that boundary values have been set in update_boundaries ! + !$omp parallel do schedule(static) do k = 1, no_nodes if (inner(k)) then - ee(:,k) = waveps + ee(:, k) = waveps endif enddo + !$omp end parallel do ! ee_ig = waveps ! - restart=1 !TODO TL: CHECK > we need this turned on right now for IG... - ! endif ! ! Initialize wave period ! + !$omp parallel do schedule(static) 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 + !$omp end parallel do ! ! 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 = sigm_ig**2 / g * (1.0 - exp(expon))**-0.4 + ! else + ! cg_ig = 0.0 kwav_ig = 0.0 + ! endif ! + ! Set Hmx and sinh(kh) for regular waves, and Hmx_ig for IG waves. Note that we use the same gamma for regular and IG waves, but this can be easily changed if needed. + ! + !$omp parallel do schedule(static) 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) = 0.88 / kwav(k) * tanh(gamma * kwav(k) * depth(k) / 0.88) + Hmx(k) = gamma * depth(k) + ! + if (igwaves) then + ! + ! Why is this different from Hmx for regular waves where we use gamma * h? + ! + 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) = gamma_ig * depth(k) + ! + endif + ! 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 - enddo - else - Hmx_ig = 0.0 - endif + !$omp end parallel do ! 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 ! + !$omp parallel do schedule(static) 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 + !$omp end parallel do ! 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 + !$omp parallel do schedule(static) + do k=1, no_nodes + ! + ctheta_ig(:,k) = sign(1.0, ctheta_ig(:,k)) * min(abs(ctheta_ig(:, k)), sigm_ig(k) / 4) + ! + enddo + !$omp end parallel do endif ! ! Solve the directional wave energy balance on an unstructured grid @@ -121,49 +154,53 @@ subroutine compute_wave_field() call timer(t2) ! call solve_energy_balance2Dstat (x,y,dhdx, dhdy, no_nodes,inner, & - w, ds, prev, & - neumannconnected, & + w, ds, prev, & + neumannconnected, & theta,ntheta,thetamean, & - depth,kwav,cg,ctheta,fw, & + depth,kwav,cg,ctheta,fw, & Tp,Tp_ig,dt,rho,alpha,gamma, gammax, & - wind, & - H,Dw,F,Df,thetam,sinhkh,& + wind, & + H,Dw,F,Df,thetam,sinhkh, & Hmx, ee, windspreadfac, u10, niter, crit, & - hmin, baldock_ratio, baldock_ratio_ig, & - aa, sig, jadcgdx, sigmin, sigmax,& + hmin, baldock_ratio, baldock_ratio_ig, baldock_exponent, & + aa, sig, jadcgdx, sigmin, sigmax, & + DoverE, relax_factor_DoverE, relax_factor_DoverA, & c_dispT, 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, & + igwaves, kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & + beta, srcig, alphaig, Dw_ig, Df_ig, qb, gam, & + steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5, & 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) + zb, nwav, ig_opt, alpha_ig, gamma_ig, gamma_fac_br, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) ! call timer(t3) ! - Fx = F*cos(thetam) - Fy = F*sin(thetam) + Fx = F * cos(thetam) + Fy = F * sin(thetam) ! end subroutine + subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & - w, ds, prev, & - neumannconnected, & + w, ds, prev, & + neumannconnected, & theta,ntheta,thetamean, & - depth,kwav,cg,ctheta,fw, & + depth,kwav,cg,ctheta,fw, & Tp,T_ig,dt,rho,alfa,gamma, gammax, & - wind, & - H,Dw,F,Df,thetam,sinhkh,& + wind, & + H,Dw,F,Df,thetam,sinhkh, & Hmx, ee, windspreadfac, u10, niter, crit, & - hmin, baldock_ratio, baldock_ratio_ig, & - aa, sig, jadcgdx, sigmin, sigmax,& + hmin, baldock_ratio, baldock_ratio_ig, baldock_exponent, & + aa, sig, jadcgdx, sigmin, sigmax, & + DoverE, relax_factor_DoverE, relax_factor_DoverA, & c_dispT, 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, & + betamean, srcig, alphaig, Dw_ig, Df_ig, qb, gam, & + steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5, & 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) + zb, nwav, ig_opt, alfa_ig, gamma_ig, gamma_fac_br, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) ! use snapwave_windsource - !use snapwave_ncoutput ! TL: removed, we don't use this in SF+SW ! implicit none ! @@ -188,13 +225,15 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4, dimension(no_nodes), intent(in) :: kwav_ig ! wave number real*4, dimension(no_nodes), intent(inout) :: cg ! group velocity real*4, dimension(ntheta,no_nodes), intent(inout):: ctheta ! refractioon speed - real*4, dimension(no_nodes), intent(in) :: cg_ig ! group velocity + real*4, dimension(no_nodes), intent(inout) :: cg_ig ! group velocity real*4, dimension(no_nodes), intent(in) :: nwav ! wave number n real*4, dimension(ntheta,no_nodes), intent(inout):: ee ! real*4, dimension(ntheta,no_nodes), intent(inout):: ee_ig ! real*4, dimension(ntheta,no_nodes), intent(in) :: ctheta_ig ! refractioon speed real*4, dimension(no_nodes), intent(in) :: fw ! wave friction factor real*4, dimension(no_nodes), intent(in) :: fw_ig ! wave friction factor + real*4, dimension(no_nodes), intent(out) :: qb ! Fraction of breaking waves according to Baldock's formulation + real*4, dimension(no_nodes), intent(out) :: gam ! Local incident wave height water depth ratio real*4, dimension(no_nodes), intent(out) :: betamean ! Mean local bed slope parameter real*4, dimension(no_nodes), intent(out) :: srcig ! Directionally averaged incident wave sink/infragravity source term real*4, dimension(no_nodes), intent(out) :: alphaig ! Mean IG shoaling parameter alpha @@ -203,9 +242,9 @@ 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(inout) :: Dw ! wave breaking dissipation real*4, dimension(no_nodes), intent(out) :: Dw_ig ! wave breaking dissipation IG real*4, dimension(no_nodes), intent(out) :: F ! wave force Dw/C/rho/h real*4, dimension(no_nodes), intent(out) :: Df ! wave friction dissipation @@ -219,7 +258,10 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4, dimension(no_nodes), intent(in) :: u10 ! wind speed and direction integer, intent(in) :: niter ! max number of iterations real*4, intent(in) :: crit ! relative accuracy for stopping criterion - integer :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx as in Leijnse et al. 2024) + integer, intent(in) :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx as in Leijnse et al. 2024) + real*4, intent(in) :: relax_factor_DoverA ! underrelaxation factor for DoverA (set to 1.0 to disable) + real*4, intent(in) :: relax_factor_DoverE ! underrelaxation factor for DoverE (set to 1.0 to disable) + real*4, intent(in) :: steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5 ! ! wind source vars ! @@ -227,6 +269,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 @@ -236,48 +279,55 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! logical, intent(in) :: vegetation ! logical yes/no real*4, dimension(no_nodes), intent(out) :: Dveg ! dissipation by vegetation: N.B. spatial field! - integer, intent(in) :: no_secveg + integer, intent(in) :: no_secveg ! number of sections in the vertical real*4, dimension(no_nodes,no_secveg), intent(in) :: veg_ah ! Height of vertical sections used in vegetation schematization [m wrt zb_ini (zb0)] real*4, dimension(no_nodes,no_secveg), intent(in) :: veg_bstems ! Width/diameter of individual vegetation stems [m] real*4, dimension(no_nodes,no_secveg), intent(in) :: veg_Nstems ! Number of vegetation stems per unit horizontal area [m-2] real*4, dimension(no_nodes,no_secveg), intent(in) :: veg_Cd ! Bulk drag coefficient [-] real*4 :: Dvegk ! dissipation by vegetation: N.B. scalar value! + real*4, dimension(no_nodes) :: Fvw ! vegetation wave drag force + real*4, dimension(no_nodes,50) :: unl ! non-linear wave orbital velocity time series, in 50 points per wave length + real*4, dimension(no_nodes,50) :: etaw0 ! non-linear sea surface time series, in 50 points per wave length ! ! ! Local variables and arrays ! integer, dimension(:), allocatable :: ok ! mask for fully iterated points - real*4 :: eemax,dtheta ! maximum wave energy density, directional resolution + integer, dimension(:), allocatable :: ok_ig ! mask for fully iterated IG points + real*4 :: dtheta ! directional resolution + real*4 :: eemax ! maximum wave energy density + real*4 :: eemax_ig ! maximum IG wave energy density real*4 :: uorbi integer :: sweep,iter ! sweep number, number of iterations integer :: k,k1,k2,count,kn,itheta ! counters (k is grid index) integer, dimension(:,:), allocatable :: indx ! index for grid sorted per sweep direction real*4, dimension(:,:), allocatable :: eeold ! wave energy density, energy density previous iteration + real*4, dimension(:,:), allocatable :: eeold_ig ! IG wave energy density, energy density previous iteration real*4, dimension(:), allocatable :: Eold ! mean wave energy, previous iteration real*4, dimension(:,:), allocatable :: srcig_local ! Energy source/sink term because of IG wave energy transfer from incident waves real*4, dimension(:,:), allocatable :: beta_local ! Local bed slope based on bed level per direction real*4, dimension(:,:), allocatable :: alphaig_local ! Local infragravity wave shoaling parameter alpha - real*4, dimension(:,:), allocatable :: depthprev ! water depth at upwind intersection point per direction + real*4, dimension(:,:), allocatable :: depthprev ! water depth at upwind intersection point per direction + real*4, dimension(:,:), allocatable :: qb_local ! local percentage of breaking waves Qb + real*4, dimension(:,:), allocatable :: gam_local ! local incident wave height over water depth ratio real*4, dimension(:), allocatable :: dee ! difference with energy previous iteration real*4, dimension(:), allocatable :: eeprev, cgprev ! energy density and group velocity at upwind intersection point real*4, dimension(:), allocatable :: eeprev_ig, cgprev_ig ! energy density and group velocity at upwind intersection point 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 real*4, dimension(:), allocatable :: E_ig ! mean wave energy real*4, dimension(:), allocatable :: diff ! maximum difference of wave energy relative to previous iteration + real*4, dimension(:), allocatable :: diff_ig ! maximum difference of IG wave energy relative to previous iteration real*4, dimension(:), allocatable :: ra ! coordinate in sweep direction - !real*4, dimension(:), allocatable :: sig real*4, dimension(:), allocatable :: sigm_ig integer, dimension(4) :: shift - real*4 :: pi = 4.*atan(1.0) - real*4 :: g=9.81 + real*4 :: pi = 4.0 * 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 :: oneoverdt real*4 :: oneover2dtheta real*4 :: rhog8 @@ -286,32 +336,37 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & real*4 :: Ek real*4 :: Hk real*4 :: percok + real*4 :: percok_ig ! percentage of converged IG points real*4 :: error + real*4 :: error_ig ! relative maximum IG wave error real*4 :: Dfk_ig real*4 :: Dwk_ig real*4 :: Ek_ig real*4 :: Hk_ig real*4 :: alfa_ig,gamma_ig ! coefficients in Baldock wave breaking dissipation model for IG waves + real*4 :: gamma_fac_br ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point real*4 :: eeinc2ig ! ratio of incident wave energy as first estimate of IG wave energy at boundary 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) integer, save :: callno=1 + integer, intent(in) :: baldock_exponent ! Exponent for multiplying the Baldock dissipation with a factor 'f = (Hloc / Hmax)**iexp' to enhance breaking when H > Hmax, with iexp = 0 (default, means unused), 1 or 2 ! - real*4, dimension(ntheta) :: sinth,costh ! distribution of wave angles and offshore wave energy density + 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 + ! allocate(ok(no_nodes)); ok=0 allocate(indx(no_nodes,4)); indx=0 allocate(eeold(ntheta,no_nodes)); eeold=0.0 @@ -322,7 +377,6 @@ 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(E(no_nodes)); E=waveps allocate(Eold(no_nodes)); Eold=0.0 ! @@ -335,11 +389,15 @@ 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 + allocate(depthprev(ntheta,no_nodes)); depthprev=0.0 + allocate(beta_local(ntheta,no_nodes)); beta_local=0.0 allocate(alphaig_local(ntheta,no_nodes)); alphaig_local=0.0 + allocate(qb_local(ntheta,no_nodes)); qb_local=0.0 + allocate(gam_local(ntheta,no_nodes)); gam_local=0.0 + allocate(eeold_ig(ntheta,no_nodes)); eeold_ig=0.0 + allocate(diff_ig(no_nodes)); diff_ig=0.0 + allocate(ok_ig(no_nodes)); ok_ig=0 endif ! if (wind) then @@ -358,55 +416,99 @@ 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) + if (igwaves) then + ok_ig = 0 + eemax_ig = maxval(ee_ig) + endif + 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 + Fvw = 0.0 + unl = 0.0 + etaw0 = 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 + ! + DoverA = 0.0 ndissip = 3.0 WsorE = 0.0 WsorA = 0.0 - Ak = waveps/sigmax + Ak = waveps / sigmax + ! + ! Re-initialise Tp at inner / neumann-connected cells only; + ! boundary cells must keep their prescribed Tp (set by + ! update_boundaries) so the wind iteration starts from the + ! correct boundary forcing. + ! + do k = 1, no_nodes + ! + if (inner(k)) Tp(k) = Tpini + if (neumannconnected(k) > 0) Tp(neumannconnected(k)) = Tpini + ! + enddo + ! + do k = 1, no_nodes + ! + sig(k) = 2 * pi / Tp(k) + sig(k) = min(max(sig(k), sigmin), sigmax) + ! + if (.not. inner(k)) then + aa(:, k) = max(ee(:, k), waveps) / sig(k) + 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)) + ! + enddo + ! 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) @@ -418,537 +520,616 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! inner(k) = .false. ! - elseif ((k1==1 .and. k2==1)) then ! TL: for now still needed for a working IG solver - inner(k)=.false. - exit - !elseif (depth(k1) < hmin .or. depth(k2) < hmin .or. (k1 == 1 .and. k2 == 1)) then - ! - ! Do not change inner here! It should be static! In a next update of the wave fields, these points may be wet. - ! - !inner(k) = .false. + elseif (k1==1 .and. k2==1) then ! TL: for now still needed for a working IG solver ! - !exit + inner(k) = .false. + exit ! endif + ! enddo enddo ! + ! Start iteration ! - ! 0-a) Set boundary and initial conditions - ! - do k = 1, no_nodes + do iter = 1, niter * 4 + ! + sweep = mod(iter, 4) !TODO - TL: problem that we don't have option for sweep = 1 anymore? + ! + if (sweep == 0) then + sweep = 4 + endif + ! + !write(*,*)'iter:', iter, 'sweep:', sweep ! - ! Boundary condition at sea side (uniform) + ! At start of each sweep, compute E, H, E_ig and H_ig, and aa ! - if (.not.inner(k)) then + !$omp parallel do private(kn) schedule(static) + 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) + E(k) = sum(ee(:, k)) * dtheta + H(k) = sqrt(8 * E(k) / rho / g) + ! + 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) + ! + ! Set Neumann boundaries + ! + if (neumannconnected(k) /= 0) then + ! + ! Do we really need all of these? Hmx? + ! + 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(k) = Df(kn) + Dw(k) = Dw(kn) + ! + endif + ! + enddo + !$omp end parallel do + ! + if (sweep == 1) then + ! + eeold = ee + ! + if (igwaves) then + ! + eeold_ig = ee_ig + ! endif ! + !$omp parallel do schedule(static) + do k = 1, no_nodes + ! + Eold(k) = sum(eeold(:, k)) + ! + enddo + !$omp end parallel do + ! 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. + ! Update IG source and sink terms ! - endif - ! - ! 0-c) Set initial condition at inner cells - ! - do k = 1, no_nodes - ! - if (inner(k)) then + 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 + ! + ! Determining of IG source term as defined in Leijnse et al. 2024 + ! + call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local, Dw, Hmx, qb_local, gam_local, gamma, gamma_fac_br, steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5) + ! + endif + ! + endif + ! + ! Loop over all points depending on sweep direction + ! + do count = 1, no_nodes + ! + k = indx(count, sweep) + ! + ! Skip non-inner (boundary) nodes + ! + if (.not. inner(k)) cycle + ! + ! For nodes below minimum depth: zero out energy and skip + ! + if (depth(k) <= hmin) then + ! + ee(:, k) = 0.0 + if (wind) aa(:, k) = 0.0 + ee_ig(:, k) = 0.0 + cycle + ! + endif + ! + ! Skip nodes that have already converged + ! + if (ok(k) == 1) cycle + ! + ! Retrieve integrated quantities computed in the start-of-sweep pre-loop + ! + Ek = E(k) + Hk = H(k) + ! + if (igwaves) then + ! + Ek_ig = E_ig(k) + Hk_ig = H_ig(k) + ! + endif + ! + ! --- Step 1: Upwind energy, group velocity (and action for wind) ------------ + ! + 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) + ! + 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) + endif + ! + if (wind) then + 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 + ! + ! --- Step 2: Pre-solve sig and celerities from upwind Ek/Ak (wind only) ---- + ! + ! The upwind Ek/Ak provides the best pre-solve estimate of local sig, so that + ! both the ee and aa matrices are assembled with a consistent cg(k). + ! Post-solve, sig is updated from the solved Ek/Ak in Step 6. ! if (wind) then - ee(:,k) = waveps - sig(k) = 2*pi/Tpini - aa(:,k) = ee(:,k)/sig(k) + ! + Ek = sum(eeprev) * dtheta + Ak = sum(aaprev) * dtheta + sig(k) = max(min(Ek / Ak, sigmax), sigmin) + Ak = Ek / sig(k) + aaprev = min(aaprev, eeprev / sigmin) + aaprev = max(aaprev, eeprev / 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)) + ! + endif + ! + ! --- Step 3: Source and sink terms ------------------------------------------ + ! + ! Bottom friction + ! + uorbi = 0.5 * sig(k) * Hk / sinhkh(k) + Dfk = 0.28 * rho * fw(k) * uorbi**3 + ! + ! Wave breaking (Baldock) + ! First check if wave breaking could occur based on Baldock criterion (Hk > baldock_ratio * Hmx(k)) + ! + if (Hk > baldock_ratio * Hmx(k)) then + ! + call baldock(rho, g, alfa, gamma, depth(k), Hk, 2 * pi / sig(k), baldock_exponent, Dwk, Hmx(k)) + ! else - ee(:,k) = waveps + ! + ! No wave breaking according to Baldock criterion + ! + Dwk = 0.0 + ! 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) + ! Vegetation + ! + 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) + ! + else + ! + Dvegk = 0.0 + ! endif ! + ! Energy dissipation ratio (with under-relaxation) + ! + DoverE(k) = (1.0 - relax_factor_DoverE) * DoverE(k) & + + relax_factor_DoverE * (Dwk + Dfk + Dvegk) / max(Ek, 1.0e-6) + ! + Df(k) = Dfk + Dw(k) = Dwk + Dveg(k) = Dvegk + ! + ! Wind source terms and action dissipation ratio DoverA (wind only) + ! 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) - ! - ! 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) + ! + 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 + ! + DwT = -c_dispT / (1.0 - ndissip) * (2.0 * 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, 1.0e-6) + else + DoverA(k) = (1.0 - relax_factor_DoverA) * DoverA(k) & + + relax_factor_DoverA * DwAk / max(Ak, 1.0e-6) + endif + ! endif ! - endif - ! - enddo - ! - ! Start iteration - ! - do iter=1,niter - ! - 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 - eeold = ee - do k = 1, no_nodes - Eold(k) = sum(eeold(:, k)) + ! --- Step 4: Assemble and solve energy balance (ee) ------------------------- + ! + do itheta = 1, ntheta + ! + R(itheta) = oneoverdt * ee(itheta, k) + cgprev(itheta) * eeprev(itheta) / ds(itheta, k) & + - srcig_local(itheta, k) * shinc2ig + ! + 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 + ! enddo ! + 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 + ! + if (wind) R(:) = R(:) + WsorE(:, k) + ! + call solve_tridiag(A, B, C, R, ee(:, k), ntheta) + ee(:, k) = max(ee(:, k), waveps) + ! + ! --- Step 5: Assemble and solve action balance (aa, wind only) --------------- + ! + ! A and C are the same as for ee (refraction terms don't change). + ! Only B_aa differs: DoverA instead of DoverE, plus upwind BC for ctheta endpoints. + ! + if (wind) then + ! + 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) + enddo + ! + if (ctheta(1, k) < 0.0) then + B_aa(1) = oneoverdt - ctheta(1, k) / dtheta + cg(k) / ds(1, k) + DoverA(k) + else + B_aa(1) = oneoverdt + cg(k) / ds(1, k) + DoverA(k) + endif + ! + R_aa(1) = oneoverdt * aa(1, k) + cgprev(1) * aaprev(1) / ds(1, k) + ! + if (ctheta(ntheta, k) > 0.0) then + B_aa(ntheta) = oneoverdt + ctheta(ntheta, k) / dtheta + cg(k) / ds(ntheta, k) + DoverA(k) + else + B_aa(ntheta) = oneoverdt + cg(k) / ds(ntheta, k) + DoverA(k) + endif + ! + R_aa(ntheta) = oneoverdt * aa(ntheta, k) + cgprev(ntheta) * aaprev(ntheta) / ds(ntheta, k) + R_aa(:) = R_aa(:) + WsorA(:, k) + ! + call solve_tridiag(A, B_aa, C, R_aa, aa(:, k), ntheta) + ! + aa(:, k) = max(aa(:, k), waveps / sigmax) + aa(:, k) = max(aa(:, k), waveps / sig(k)) + ! + endif + ! + ! --- Step 6: Depth-limit energy (and action), update sig and celerities ----- + ! + Ek = sum(ee(:, k)) * dtheta + depthlimfac = max(1.0, (sqrt(Ek / rhog8) / (gammax * depth(k)))**2) + ee(:, k) = ee(:, k) / depthlimfac + ! +! if (wind) then +! ! +! Ek = Ek / depthlimfac +! Ak = sum(aa(:, k)) * dtheta +! Ak = Ak / depthlimfac +! aa(:, k) = aa(:, k) / depthlimfac +! sig(k) = max(min(Ek / Ak, sigmax), sigmin) +! ! +! 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 + ! + ! --- Step 7: IG wave balance (optional) ------------------------------------- + ! 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) - ! + ! + ! Update incident Hk from post-solve ee (needed for IG bottom friction) + ! + Hk = sqrt(8.0 * sum(ee(:, k)) * dtheta / rho / g) + Dfk_ig = fw_ig(k) * 0.0361 * (9.81 / depth(k))**1.5 * Hk * Ek_ig + ! + ! IG wave breaking (Baldock) + ! + 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), baldock_exponent, Dwk_ig, Hmx_ig(k)) + else + Dwk_ig = 0.0 + endif + ! + Df_ig(k) = Dfk_ig + Dw_ig(k) = Dwk_ig + ! + ! Not using underrelaxation for IG dissipation for now, but we could add this if needed (relax_factor_DoverE_ig) + DoverE_ig(k) = (Dwk_ig + Dfk_ig) / max(Ek_ig, 1.0e-6) + ! + ! + ! IG RHS + ! + 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) + enddo + ! + ! IG matrix with directional boundary conditions + ! + 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 + enddo + ! + if (ctheta_ig(1, k) < 0.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 + else + A_ig(1) = 0.0 + B_ig(1) = oneoverdt + cg_ig(k) / ds(1, k) + DoverE_ig(k) + C_ig(1) = 0.0 endif ! - endif + if (ctheta_ig(ntheta, k) > 0.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) + C_ig(ntheta) = 0.0 + endif + ! + call solve_tridiag(A_ig, B_ig, C_ig, R_ig, ee_ig(:, k), ntheta) + ee_ig(:, k) = max(ee_ig(:, k), 0.0) + ! + ! Depth-limit IG energy + ! + depthlimfac = max(1.0, (sqrt(sum(ee_ig(:, k)) * dtheta / rhog8) / (gammax * depth(k)))**2) + ee_ig(:, k) = ee_ig(:, k) / depthlimfac + ! + else + ! + ee_ig(:, k) = 0.0 + ! + endif ! - endif - ! - ! Loop over all points depending on sweep direction + enddo ! - do count = 1, no_nodes + if (sweep==4) then ! - k=indx(count, sweep) + ! Check convergence after all 4 sweeps ! - if (inner(k)) then - if (depth(k)>1.1*hmin) then + !$omp parallel do private(dee) schedule(static) + do k = 1, no_nodes + ! + dee = ee(:, k) - eeold(:, k) + diff(k) = maxval(abs(dee)) + ! + if (diff(k) / eemax < crit) then + ok(k) = 1 + endif + ! + enddo + !$omp end parallel do + ! + ! Percentage of converged points + ! + percok = sum(ok) / dble(no_nodes) * 100.0 + eemax = maxval(ee) + ! + ! Relative maximum error + ! + error = maxval(diff) / eemax + ! + ! Check convergence of IG waves + ! + if (igwaves) then + ! + !$omp parallel do private(dee) schedule(static) + do k = 1, no_nodes ! - if (ok(k) == 0) then - ! - ! Only perform computations on wet inner points that are not yet converged (ok) - ! - 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) - ! - 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) - endif - ! - if (wind) then - aaprev(itheta) = w(1, itheta, k)*aa(itheta, k1) + w(2, itheta, k)*aa(itheta, k2) - endif - ! - enddo - ! - Ek = sum(eeprev)*dtheta ! to check - ! - depthlimfac = max(1.0, (sqrt(Ek/rhog8)/(gammax*depth(k)))**2.0) - Hk = min(sqrt(Ek/rhog8), gamma*depth(k)) - Ek = Ek/depthlimfac - ! - 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) - else - Dvegk = 0. - endif - ! - DoverE(k) = (Dwk + Dfk + Dvegk)/max(Ek, 1.0e-6) - ! - if (wind) 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 - ! - 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) ) - ! - 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)) - ! - endif - ! - do itheta = 1, ntheta - ! - R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta)*eeprev(itheta)/ds(itheta, k) - srcig_local(itheta, k) * shinc2ig - ! - 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 - ! - enddo - ! - 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 - ! - ! Solve tridiagonal system per point - ! - if (wind) then - 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) - 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) - 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) - 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) - 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) - 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 - ! - 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 - ! - 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) - ! - 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 - ! - ! 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 - ! - ! 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)) - else - Dwk_ig = 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? - ! - 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 - ! - 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 - ! - enddo - ! - 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 - 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 - 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) - C_ig(ntheta) = 0.0 - else - A_ig(ntheta) = 0.0 - B_ig(ntheta) = oneoverdt + cg_ig(k)/ds(ntheta, k) + DoverE_ig(k) - C_ig(ntheta) = 0.0 - endif - ! - ! Solve tridiagonal system per point - ! - call solve_tridiag(A_ig, B_ig, C_ig, R_ig, ee_ig(:,k), ntheta) - ee_ig(:, k) = max(ee_ig(:, k), 0.0) - ! - else - ! - ee_ig(:, k) = 0.0 - ! - endif - ! + dee = ee_ig(:, k) - eeold_ig(:, k) + diff_ig(k) = maxval(abs(dee)) + ! + if (diff_ig(k) / eemax_ig < crit) then + ok_ig(k) = 1 endif ! - else + enddo + !$omp end parallel do + ! + percok_ig = sum(ok_ig) / dble(no_nodes) * 100.0 + eemax_ig = maxval(ee_ig) + error_ig = maxval(diff_ig) / eemax_ig + ! + write(logstr,'(a,i6,a,f10.5,a,f7.2,a,f10.5,a,f7.2)')' iteration ', iter / 4 , & + ' error = ', error,' %ok = ', percok,' error_ig = ', error_ig,' %ok_ig = ', percok_ig + call write_log(logstr, 0) + ! + if ((error < crit .or. percok > 99.0)) then!.and. (error_ig < crit .or. percok_ig > 99.0)) then + ! + write(logstr,'(a,i6,a,f10.5,a,f7.2,a,f10.5,a,f7.2)')' converged at iteration ', iter / 4 , & + ' error = ', error,' %ok = ', percok,' error_ig = ', error_ig,' %ok_ig = ', percok_ig + call write_log(logstr, 0) + exit + ! + elseif (iter == niter * 4) then ! Made it to the end without reaching 'error 0 .and. Hloc > Hmax) then + ! + ! Add extra dissipation when Hloc exceeds Hmax. + ! This is needed at very steep coast lines, where Baldock dissipation cannot always keep up with + ! the wave height increase due to shoaling. The extra dissipation is added by multiplying + ! the Baldock dissipation with a factor f, which is larger than 1 when Hloc > Hmax. + ! + f = (Hloc / Hmax)**iexp + ! + else + ! + f = 1.0 + ! endif ! + Dw = 0.28 * alfa * rho * g / T * exp( - (Hmax / Hloc)**2) * (Hmax**2 + Hloc**2) * f + ! + ! Other options for wave breaking dissipation (not used, but left here for reference) + ! + ! Dw = 0.28 * alfa * rho * g / T * exp( - (Hmax / Hloc)**2) * (Hmax**3 + Hloc**3) / gamma / depth + ! end subroutine baldock - - subroutine 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) + + subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local, Dw, Hmx, qb_local, gam_local, gamma, gamma_fac_br, steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5) ! + ! + ! Determining of IG source term as defined in Leijnse et al. 2024 + ! + ! inout: alphaig_local, srcig_local, beta_local + ! in: the rest + ! + ! NOTE - This is based on the energy in the previous SnapWave timestep 'ee' and 'ee_ig', and waveheight 'H', which should therefore be made available. + ! implicit none ! ! Incoming variables + ! logical, dimension(no_nodes), intent(in) :: inner ! mask of inner grid points (not on boundary) integer, intent(in) :: no_nodes,ntheta ! number of grid points, number of directions real*4, dimension(2,ntheta,no_nodes),intent(in) :: w ! weights of upwind grid points, 2 per grid point and per wave direction real*4, dimension(ntheta,no_nodes), intent(in) :: ds ! distance to interpolated upwind point, per grid point and direction integer, dimension(2,ntheta,no_nodes),intent(in) :: prev ! two upwind grid points per grid point and wave direction - real*4, dimension(no_nodes), intent(in) :: cg_ig ! group velocity + real*4, dimension(no_nodes), intent(inout) :: cg_ig ! group velocity real*4, dimension(no_nodes), intent(in) :: nwav ! wave number n real*4, dimension(no_nodes), intent(in) :: depth ! water depth real*4, dimension(no_nodes), intent(in) :: zb ! actual bed level - real*4, dimension(no_nodes), intent(in) :: H ! wave height + real*4, dimension(no_nodes), intent(in) :: H ! wave height real*4, dimension(ntheta,no_nodes), intent(in) :: ee ! energy density - real*4, dimension(ntheta,no_nodes), intent(in) :: ee_ig ! energy density infragravity waves + real*4, dimension(ntheta,no_nodes), intent(in) :: ee_ig ! energy density infragravity waves integer, intent(in) :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) real*4, intent(in) :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0 + real*4, intent(in) :: dtheta ! directional resolution + real*4, intent(in) :: gamma ! coefficients in Baldock wave breaking dissipation + real*4, intent(in) :: gamma_fac_br ! factor times gamma that is used to determine the maximum incident wave breaking point in the surf zone using local incident wave height over water depth ratio, among others used to set the IG source term to 0 shallower than this point + real*4, dimension(no_nodes), intent(in) :: Dw ! wave breaking dissipation + real*4, dimension(no_nodes), intent(in) :: Hmx ! Hmax + real*4, intent(in) :: steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5 ! ! Inout variables + ! real*4, dimension(:,:), intent(inout) :: alphaig_local ! Local infragravity wave shoaling parameter alpha real*4, dimension(:,:), intent(inout) :: srcig_local ! Energy source/sink term because of IG wave shoaling - real*4, dimension(:), intent(inout) :: eeprev, cgprev ! energy density and group velocity at upwind intersection point - real*4, dimension(:), intent(inout) :: eeprev_ig ! energy density at upwind intersection point - real*4, dimension(ntheta,no_nodes), intent(inout):: beta_local ! Local bed slope based on bed level per direction + real*4, dimension(ntheta,no_nodes), intent(inout):: beta_local ! Local bed slope based on bed level per direction + real*4, dimension(ntheta,no_nodes), intent(inout):: qb_local ! local percentage of breaking waves Qb + real*4, dimension(ntheta,no_nodes), intent(inout):: gam_local ! local incident wave height over water depth ratio ! ! Internal variables + ! integer :: itheta ! directional counter - integer :: k ! counters (k is grid index) + integer :: k ! counters (k is grid index) integer :: k1,k2 ! upwind counters (k is grid index) - real*4 :: gam ! local gamma (Hinc / depth ratio) - real*4, dimension(ntheta,no_nodes) :: depthprev ! water depth at upwind intersection point - real*4, dimension(ntheta,no_nodes) :: Sxx ! Radiation Stress - real*4, dimension(:), allocatable :: Sxxprev ! radiation stress at upwind intersection point - real*4, dimension(:), allocatable :: Hprev ! Incident wave height at upwind intersection point + real*4 :: gam ! local gamma (Hinc / depth ratio) + real*4, dimension(ntheta,no_nodes) :: depthprev ! water depth at upwind intersection point + real*4, dimension(no_nodes) :: Sxx ! radiation Stress + real*4, dimension(ntheta) :: Sxxprev ! radiation stress at upwind point + real*4, dimension(ntheta) :: Hprev ! wave height at upwind point + real*4, dimension(ntheta) :: cgprev ! group velocity at upwind point + real*4, dimension(ntheta) :: Eprev ! Mean incident wave energy at upwind intersection point + real*4, dimension(ntheta) :: Eprev_ig ! Mean infragravity wave energy at upwind intersection point + real*4, dimension(no_nodes) :: E_local ! mean wave energy waves - just local + real*4, dimension(no_nodes) :: E_ig_local ! mean wave energy infragravity waves - just local real*4 :: dSxx ! difference in Radiation stress - real*4 :: Sxx_cons ! conservative estimate of radiation stress using conservative shoaling - ! - ! Allocate internal variables - allocate(Sxxprev(ntheta)) - allocate(Hprev(ntheta)) + real*4 :: Sxx_cons ! conservative estimate of radiation stress using conservative shoaling + real*4 :: delta_Dw ! difference of Dw compared to upwind point, to get sign for max breaking point + real*4 :: Qb ! Percentage of breaking incident waves + real*4 :: transition_factor ! Transition factor for letting srcig go to zero smoothly, around gamma*gamma_fac_br + real*4 :: transition_factor_width_1 ! Width factor of generalized (Fermi�Dirac style) transfer function with adjustable midpoint and width + real*4 :: transition_factor_width_2 ! Width factor of generalized (Fermi�Dirac style) transfer function with adjustable midpoint and width + real*4 :: gamma_fac_br_transition ! Transitioned version of gamma_fac_br, so that for steep slopes it remains 1.0 + real*4 :: beta_limit_1 ! Cut-off beta_local for end of validity alphaig formulation of Leijnse et al. 2024 + real*4 :: beta_limit_2 ! Beta_local limit for transition function + ! + ! Set internal variables ! Sxx = 0.0 + Hprev = 0.0 + Eprev = 0.0 + Eprev_ig = 0.0 + ! + E_local = 0.0 + E_ig_local = 0.0 + Sxx = 0.0 + ! + ! Used is generalized (Fermi�Dirac style) transfer function with adjustable midpoint and width + ! + transition_factor_width_1 = 0.005 + transition_factor_width_2 = 0.002 + beta_limit_1 = 0.07 + !beta_limit_2 = beta_limit_1 - 0.01 + beta_limit_2 = beta_limit_1 - 0.02 + ! + ! Pre-compute Sxx for all nodes + ! + !$omp parallel do schedule(static) + do k = 1, no_nodes + ! + if (inner(k)) then !TODO: check whether should be on only 'inner' or not + ! + ! Update E (not saved from previous timestep) + ! + E_local(k) = sum(ee(:,k)) * dtheta + ! + ! Update E_ig (not saved from previous timestep) + ! + E_ig_local(k) = sum(ee_ig(:, k)) * dtheta + ! + endif + ! + Sxx(k) = ((2.0 * max(0.0, min(1.0, nwav(k)))) - 0.5) * E_local(k) + ! + enddo + !$omp end parallel do + ! + ! Actual computation of srcig + ! + ! Main loop: compute IG source/sink term per node. + ! All writes target the column (itheta, k), so the loop is data-independent across k. + ! Per-k scratch arrays (cgprev, Eprev, Eprev_ig, Sxxprev, Hprev) are + ! listed as private so each thread gets its own copy on the stack. ! + !$omp parallel do & + !$omp& private(itheta, k1, k2, gam, dSxx, Sxx_cons, & + !$omp& cgprev, Eprev, Eprev_ig, Sxxprev, Hprev) & + !$omp& schedule(static) do k = 1, no_nodes ! - if (inner(k)) then + if (inner(k)) then ! - ! Compute exchange source term inc to ig waves - per direction + ! Compute exchange source term inc to ig waves - per direction ! do itheta = 1, ntheta ! k1 = prev(1, itheta, k) k2 = prev(2, itheta, k) ! - if (k1>0 .and. k2>0) then ! IMPORTANT - for some reason (k1*k2)>0 is not reliable always, resulting in directions being uncorrectly skipped!!! + if (k1 > 0 .and. k2 > 0) then ! IMPORTANT - for some reason (k1*k2)>0 is not reliable always, resulting in directions being uncorrectly skipped!!! ! ! First calculate upwind direction dependent variables - depthprev(itheta,k) = w(1, itheta, k)*depth(k1) + w(2, itheta, k)*depth(k2) - ! - beta_local(itheta,k) = max((w(1, itheta, k)*(zb(k) - zb(k1)) + w(2, itheta, k)*(zb(k) - zb(k2)))/ds(itheta, k), 0.0) + ! + depthprev(itheta,k) = w(1, itheta, k) * depth(k1) + w(2, itheta, k) * depth(k2) + ! + beta_local(itheta,k) = max((w(1, itheta, k) * (zb(k) - zb(k1)) + w(2, itheta, k) * (zb(k) - zb(k2))) / ds(itheta, k), 0.0) + ! + ! FIXME - shorter, but also same result? + !beta_local(itheta,k) = max(zb(k) - (w(1, itheta, k) * zb(k1) + w(2, itheta, k) * zb(k2))/ds(itheta, k), 0.0) ! ! Notes: ! - use actual bed level now for slope, because depth changes because of wave setup/tide/surge ! - in zb, depth is negative > therefore zb(k) minus zb(k1) ! - beta=0 means a horizontal or decreasing slope > need alphaig=0 then in IG src/sink term ! - !betan_local(itheta,k) = (beta/sigm_ig)*sqrt(9.81/max(depth(k), hmin)) ! TL: in case in the future we would need the normalised bed slope again + !betan_local(itheta,k) = (beta/sigm_ig)*sqrt(9.81/max(depth(k), hmin)) ! TL: in case in the future we would need the normalised bed slope again ! - ! TL - Note: cg_ig = cg - cgprev(itheta) = w(1, itheta, k)*cg_ig(k1) + w(2, itheta, k)*cg_ig(k2) - ! - Sxx(itheta,k1) = ((2.0 * max(0.0,min(1.0,nwav(k1)))) - 0.5) * ee(itheta, k1) ! limit so value of nwav is between 0 and 1 - Sxx(itheta,k2) = ((2.0 * max(0.0,min(1.0,nwav(k2)))) - 0.5) * ee(itheta, k2) ! limit so value of nwav is between 0 and 1 + ! Fraction of breaking waves, based on H(k) + !Qb = min(max(exp(-(Hmx(k)/H(k))**2), 0.0), 1.0) ! Qb percentage of breaking waves according to Baldock's formulation, between 0 and 1 + ! Base on upwind point: + Qb = min(max(exp(-((w(1, itheta, k)*Hmx(k1) + w(2, itheta, k)*Hmx(k2)) / Hprev(itheta))**2), 0.0), 1.0) ! Qb percentage of breaking waves according to Baldock's formulation, between 0 and 1 ! - Sxxprev(itheta) = w(1, itheta, k)*Sxx(itheta,k1) + w(2, itheta, k)*Sxx(itheta,k2) + qb_local(itheta, k) = Qb + ! + cgprev(itheta) = w(1, itheta, k) * cg_ig(k1) + w(2, itheta, k) * cg_ig(k2) ! - eeprev(itheta) = w(1, itheta, k)*ee(itheta, k1) + w(2, itheta, k)*ee(itheta, k2) - eeprev_ig(itheta) = w(1, itheta, k)*ee_ig(itheta, k1) + w(2, itheta, k)*ee_ig(itheta, k2) + Sxxprev(itheta) = w(1, itheta, k) * Sxx(k1) + w(2, itheta, k) * Sxx(k2) + ! + Eprev(itheta) = w(1, itheta, k) * E_local(k1) + w(2, itheta, k) * E_local(k2) + Eprev_ig(itheta) = w(1, itheta, k) * E_ig_local(k1) + w(2, itheta, k) * E_ig_local(k2) + ! + Hprev(itheta) = w(1, itheta, k) * H(k1) + w(2, itheta, k) * H(k2) ! - Hprev(itheta) = w(1, itheta, k)*H(k1) + w(2, itheta, k)*H(k2) - ! ! Determine relative waterdepth 'gam' ! - gam = max(0.5*(Hprev(itheta)/depthprev(itheta,k) + H(k)/depth(k)), 0.0) ! mean gamma over current and upwind point + gam = max(0.5 * (Hprev(itheta) / depthprev(itheta,k) + H(k) / depth(k)), 0.0) ! mean gamma over current and upwind point + ! + gam_local(itheta, k) = gam + ! + ! Free waves and no IG source/sink term if incident waves start breaking + ! + ! Adjust cg_ig for free infragravity waves release in surfzone + ! TL - Note: cg_ig = cg + !if (ig_opt == X) then + ! ! + ! if (gam > (gamma_fac_br * gamma)) then + ! ! + ! cg_ig(k) = sqrt(9.81 * depth(k)) + ! ! + ! endif + ! ! + !endif ! ! Determine dSxx and IG source/sink term 'srcig' ! - if (ig_opt == 1 .or. ig_opt == 2) then + if (ig_opt == 1 .or. ig_opt == 2 .or. ig_opt == 11 .or. ig_opt == 12 .or. ig_opt == 13 .or. ig_opt == 14 .or. ig_opt == 15) then ! ! Calculate shoaling parameter alpha_ig following Leijnse et al. (2024) ! + if (ig_opt == 11 .or. ig_opt == 12 .or. ig_opt == 13) then + ! + ! Limit beta to max 0.07 (=beta_limit_1) before going into alphaig parametrisation + ! + beta_local(itheta,k) = min(beta_local(itheta,k), beta_limit_1) + ! + endif + ! call estimate_shoaling_parameter_alphaig(beta_local(itheta,k), gam, alphaig_local(itheta,k)) ! [input, input, output] + ! + ! Steep slope addition + ! + if (ig_opt == 14 .or. ig_opt == 15) then + ! + call estimate_shoaling_parameter_alphaig_steep_slopes(beta_local(itheta, k), gam, alphaig_local(itheta, k), steep_fac1, steep_fac2, steep_fac3, steep_fac4, steep_fac5) + ! [input, input, inout, input, input, input, input, input] + ! + endif ! ! Now calculate source term component ! ! Newest dSxx/dx based method, using estimate of Sxx(k) using conservative shoaling - if (Sxxprev(itheta)<=0.0) then + ! + if (Sxxprev(itheta) <= 0.0) then + ! + srcig_local(itheta, k) = 0.0 !Avoid big jumps in dSxx that can happen if a upwind point is a boundary point with Hinc=0 + ! + else + ! + if (ig_opt == 1 .or. ig_opt == 11 .or. ig_opt == 12 .or. ig_opt == 13 .or. ig_opt == 14 .or. ig_opt == 15) then ! Option using conservative shoaling for dSxx/dx + ! + ! Calculate Sxx based on conservative shoaling of upwind point's energy: + ! Sxx_cons = E(i-1) * Cg(i-1) / Cg * (2 * n(i) - 0.5) + Sxx_cons = Eprev(itheta) * cgprev(itheta) / cg_ig(k) * ((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5) + ! Note - limit so value of nwav is between 0 and 1, and Sxx therefore doesn't become NaN for nwav=Infinite + ! + dSxx = Sxx_cons - Sxxprev(itheta) + ! + elseif (ig_opt == 2) then ! Option taking actual difference for dSxx/dx + ! + dSxx = Sxx(itheta) - Sxxprev(itheta) + endif ! - srcig_local(itheta, k) = 0.0 !Avoid big jumps in dSxx that can happen if a upwind point is a boundary point with Hinc=0 + dSxx = max(dSxx, 0.0) ! - else - ! - if (ig_opt == 1) then ! Option using conservative shoaling for dSxx/dx + !if (ig_opt == 1 .or. ig_opt == 2.or. ig_opt == 12 .or. ig_opt == 13) then ! - ! Calculate Sxx based on conservative shoaling of upwind point's energy: - ! Sxx_cons = E(i-1) * Cg(i-1) / Cg * (2 * n(i) - 0.5) - Sxx_cons = eeprev(itheta) * cgprev(itheta) / cg_ig(k) * ((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5) - ! Note - limit so value of nwav is between 0 and 1, and Sxx therefore doesn't become NaN for nwav=Infinite + ! Base on E_prev_ig instead of eeprev_ig(itheta) > no bins but total energy + ! NOTE - already here multiplied with ee(itheta,k), for direct inclusion in 'R'-term + srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(Eprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k) /max(E_local(k), 1.0e-6) * ee(itheta,k) ! - dSxx = Sxx_cons - Sxxprev(itheta) + !elseif (ig_opt == 20) then + ! + ! NOTE - in main script this is multiplied with ee(itheta,k) to get directional energy, for direct inclusion in 'B'-term + ! + !srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(Eprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k) /max(E_local(k), 1.0e-6) !* ee(itheta,k) + !endif ! - elseif (ig_opt == 2) then ! Option taking actual difference for dSxx/dx + ! Limit srcig to 0 after waves start (significantly) breaking, as defined here as gam=Hrms,inc / h > (gamma_fac_br * gamma) ! - dSxx = Sxx(itheta,k) - Sxxprev(itheta) + ! Ergo, it is assumed that after this point IG waves are free, and no bound wave forcing is happening anymore, so srcig should be 0 from here on + ! + if (ig_opt == 12) then + ! + ! Let srcig transition to 0 more smoothly using fac_transition that reduced from 1 to 0 around gamma_fac_br * snapwave_gamma + ! Similar as before, but then smooth: + ! ! Note - gam is in Hrms + ! if (gam > (gamma_fac_br * gamma)) then + ! ! + ! srcig_local(itheta, k) = 0.0 + ! ! + ! endif + ! + transition_factor = 1.0 - (1.0 / (1.0 + exp(- (gam - (gamma_fac_br * gamma)) / transition_factor_width_1))) + ! + srcig_local(itheta, k) = transition_factor * srcig_local(itheta, k) + ! + elseif (ig_opt == 13 .or. ig_opt == 15) then + ! + ! Let srcig transition to 0 more smoothly using fac_transition that reduced from 1 to 0 around gamma_fac_br * snapwave_gamma + ! + ! But, only for beta_local < 0.07, so adjust based on beta_local so that transition_factor = 1.0 for Beta_local = 0.07 + ! + gamma_fac_br_transition = gamma_fac_br + ((1-gamma_fac_br) / (1 + exp(- (beta_local(itheta,k) - beta_limit_2) / transition_factor_width_2))) + ! + transition_factor = 1.0 - (1.0 / (1.0 + exp(- (gam - (gamma_fac_br_transition * gamma)) / transition_factor_width_1))) + ! + srcig_local(itheta, k) = transition_factor * srcig_local(itheta, k) + ! + endif ! - endif - ! - dSxx = max(dSxx, 0.0) - ! - srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(eeprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k) - ! endif ! else ! TL: option to add future parameterisations here for e.g. coral reef type coasts @@ -1153,15 +1514,17 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d ! endif ! - enddo + enddo ! endif ! - enddo - ! + enddo + !$omp end parallel do + ! end subroutine determine_infragravity_source_sink_term subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) + ! real*4, intent(in) :: beta real*4, intent(in) :: gam real*4, intent(out) :: alphaig @@ -1171,11 +1534,13 @@ subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) ! Estimate shoaling parameter alphaig - as in Leijnse et al. (2024) ! ! Determine constants + ! beta1 = 0.016993 beta2 = 0.5 beta3 = 17.7104 beta4 = 1 beta5 = 0.7 + !beta5 = 0.5 beta6 = 0.11841 beta7 = 0.34037 ! @@ -1190,8 +1555,10 @@ subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) alphaig = exp(-beta3 * beta ** beta4) * ((beta5 - gam) * beta6 + (beta7 - gam) * (beta1 / beta ** beta2)) ! elseif (gam >= beta7) then ! shallow water - for gam>0.7 the fit automatically goes to 0 + !elseif (gam >= beta7 .and. gam < 0.5) then ! shallow water - for gam>0.7 the fit automatically goes to 0 ! - alphaig = exp(-beta3 * beta ** beta4) * (max(beta5 - gam, 0.0)) * beta6 + alphaig = exp(-beta3 * beta ** beta4) * (max(beta5 - gam, 0.0)) * beta6 + !alphaig = exp(-beta3 * beta ** beta4) * (max(0.5 - gam, 0.0)) * beta6 ! else ! for safety, but negative gamma should not occur ! @@ -1208,6 +1575,72 @@ subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) ! end subroutine estimate_shoaling_parameter_alphaig + + subroutine estimate_shoaling_parameter_alphaig_steep_slopes(beta, gam, alphaig, fac1, fac2, fac3, fac4, fac5) + ! [input, input, inout, input, input, input, input, input] + real*4, intent(in) :: beta + real*4, intent(in) :: gam + real*4, intent(inout) :: alphaig + ! + real*4, intent(in) :: fac1, fac2, fac3, fac4, fac5 + real*4 :: alphaig_steep + ! + ! Estimate shoaling parameter alphaig - for steep slopes with beta > 0.07 + ! These were not covered in training dataset of Leijnse et al. 2024 + ! + !alphaig_total = alphaig Eq11 + alphaig_steep + !alphaig_steep = 0.1 * max(beta-0.07, 0)^0.6 * max(1.0-gamma, 0) + !alphaig_steep = fac2 * max(beta-fac3, 0)**fac4 * max(fac5-gam, 0) + ! + ! Determine constants + ! + !fac1 = 0.3 !Cut-off gamma, below this alphaig_steep = 0 + !fac1 = 0.0 !Cut-off gamma, below this alphaig_steep = 0 + !fac2 = 0.1 ! Multiplication factor + !fac3 = 0.07 ! Cut-off beta, below this alphaig_steep = 0, and above this it increases with beta + !fac4 = 0.6 ! Exponent + !fac5 = 1.0 ! Cut-off gamma, above this alphaig_steep = 0 + ! + ! For deep water or negative slope, alphaig_steep = 0 + ! + alphaig_steep = 0.0 + ! + !if (beta > 0.0) then + ! write(*,*)'Starting alphaig steep slope addition, beta, gam, alphaig before', beta, gam, alphaig + !endif + ! + ! If positively increasing local bed slope beta + ! + if (beta > 0.0) then + ! + if (gam >= fac1) then ! shallow(er) water - for gam>1.0 (=fac5) the fit automatically goes to 0 + ! + alphaig_steep = fac2 * (max(beta - fac3, 0.0) ** fac4) * (max(fac5 - gam, 0.0)) + ! + else ! for gam < fac1 + ! + alphaig_steep = 0.0 + ! + endif + ! + endif + ! + !if (beta > 0.0) then + ! write(*,*)'Calculated alphaig_steep addition', alphaig_steep + !endif + ! + ! Combine + ! + alphaig = alphaig + alphaig_steep + ! + ! Limit total alphaig between [0, 1] to prevent large overshoots in case of low gamma and very small beta + ! + alphaig = max(alphaig, 0.0) + alphaig = min(alphaig, 1.0) + ! + end subroutine estimate_shoaling_parameter_alphaig_steep_slopes + + subroutine hpsort_eps_epw (n, ra, ind, eps) !--------------------------------------------------------------------- ! sort an array ra(1:n) into ascending order using heapsort algorithm, @@ -1352,8 +1785,8 @@ end subroutine hpsort_eps_epw subroutine timer(t) real*4,intent(out) :: t integer*4 :: count,count_rate,count_max - call system_clock (count,count_rate,count_max) - t = real(count)/count_rate + call system_clock (count,count_rate, count_max) + t = real(count) / count_rate end subroutine timer subroutine vegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, depth, rho, g, H, Dveg) @@ -1364,7 +1797,7 @@ subroutine vegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nstem ! declare variables real*4, intent(in) :: sigm ! wave frequency (per cell) integer, intent(in) :: no_nodes ! number of unstructured grid nodes - integer, intent(in) :: no_secveg + integer, intent(in) :: no_secveg ! number of sections in the vertical real*4, dimension(no_secveg), intent(in) :: veg_ah ! Height of vertical sections used in vegetation schematization [m wrt zb_ini (zb0)] (per cell) real*4, dimension(no_secveg), intent(in) :: veg_bstems ! Width/diameter of individual vegetation stems [m] (per cell) real*4, dimension(no_secveg), intent(in) :: veg_Nstems ! Number of vegetation stems per unit horizontal area [m-2] (per cell) @@ -1377,7 +1810,7 @@ subroutine vegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nstem integer :: m real*4, intent(in) :: kwav ! wave number (per cell) real*4, intent(out) :: Dveg ! dissipation by vegetation (per cell) - + ! Set dissipation in vegetation to zero everywhere for a start Dveg = 0.d0 @@ -1389,16 +1822,19 @@ subroutine vegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nstem if (no_secveg > 0) then ! only in case vegetation is present do m=1,no_secveg ! for each vertical vegetation section if (veg_Cd(m) < 0.d0) then ! If Cd is not user specified: call subroutine of M. Bendoni (see below) - write(logstr,*)'Cd is not user specified: using subroutine bulkdragcoeff to compute Cd' - call write_log(logstr, 0) ! - call bulkdragcoeff(veg_ah(m),m,Cdterm,no_nodes,no_secveg,depth,H,kwav,veg_bstems(m),sigm) ! bulkdragcoeff(ahveg(k,m)+zb0(k)-zb(k),m,k,Cdterm) <- no bed level change implemented in Snapwave - !write(*,*)'Cd is not user specified: putting default value of 0.7' - !veg_Cd(k,m) = 0.7 + !call bulkdragcoeff(veg_ah(m),m,Cdterm,no_nodes,no_secveg,depth,H,kwav,veg_bstems(m),sigm) ! bulkdragcoeff(ahveg(k,m)+zb0(k)-zb(k),m,k,Cdterm) <- no bed level change implemented in Snapwave + !write(logstr,*)'Cd is not user specified: using m. bendoni bulkdragcoefficient to compute cd: ',cdterm + !veg_Cd(m) = Cdterm + ! + write(logstr,*)'SnapWave ERROR - Cd is not specified for layer: ',m + call write_log(logstr, 0) + ! + ! endif enddo endif - + ! ! Attenuation by vegetation is computed in wave action balance (swvegatt) and the momentum balance (momeqveg); ! 1) Short wave dissipation by vegetation @@ -1412,9 +1848,9 @@ end subroutine vegatt subroutine swvegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, depth, rho, g, H, Dveg)! Short wave dissipation by vegetation !use snapwave_data !use snapwave_domain - + ! implicit none - + ! ! declare variables integer, intent(in) :: no_nodes ! number of unstructured grid nodes integer, intent(in) :: no_secveg @@ -1427,20 +1863,20 @@ subroutine swvegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nst real*4, intent(in) :: rho real*4, intent(in) :: g real*4, intent(in) :: H ! wave height - + ! ! local variables real*4 :: pi ! 3.14159 integer :: k,m ! indices of actual x,y point - + ! real*4 :: aht,hterm,htermold,Dvgt,ahtold real*4 :: Dvg,kmr!,kwav real*4, intent(in) :: kwav!,k - + ! real*4, intent(out) :: Dveg - + ! pi = 4.d0*atan(1.d0) kmr = min(max(kwav, 0.01d0), 100.d0) - + ! ! Set dissipation in vegetation to zero everywhere for a start Dvg = 0.d0 Dvgt = 0.d0 @@ -1448,24 +1884,24 @@ subroutine swvegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nst ahtold = 0.d0 if (no_secveg>0) then ! only if vegetation is present do m=1,no_secveg - + ! ! Determine height of vegetation section (restricted to current bed level) !aht = veg(ind)%ah(m)+ahtold !+s%zb0(k,j)-s%zb(k,j)!(max(veg(ind)%zv(m)+s%zb0(k,j),s%zb(k,j))) aht = veg_ah(m)+ahtold - + ! ! restrict vegetation height to local water depth aht = min(aht, depth) - + ! ! compute hterm based on ah hterm = (sinh(kmr*aht)**3+3*sinh(kmr*aht))/(3.d0*kmr*cosh(kmr* depth)**3) ! - + ! ! compute dissipation based on aht and correct for lower elevated dissipation layers (following Suzuki et al. 2012) Dvgt = 0.5d0/sqrt(pi)*rho*veg_Cd(m)*veg_bstems(m)*veg_Nstems(m)*(0.5d0*kmr*g/sigm)**3*(hterm-htermold)*H**3 - + ! ! save hterm to htermold to correct possibly in next vegetation section htermold = hterm ahtold = aht - + ! ! add dissipation current vegetation section Dvg = Dvg + Dvgt enddo @@ -1473,27 +1909,29 @@ subroutine swvegatt(sigm, no_nodes, kwav, no_secveg, veg_ah, veg_bstems, veg_Nst Dveg = Dvg end subroutine swvegatt - subroutine bulkdragcoeff(ahh, m, Cdterm, no_nodes, no_secveg, depth, H, kwav, veg_bstems, sigm)!(ahh,m,i,Cdterm) + subroutine bulkdragcoeff(ahh, m, Cdterm, no_nodes, no_secveg, depth, H, kwav, veg_bstems, sigm) !(ahh,m,i,Cdterm) + ! ! Michele Bendoni: subroutine to calculate bulk drag coefficient for short wave ! energy dissipation based on the Keulegan-Carpenter number (adapted from XBeach) ! Ozeren et al. (2013) or Mendez and Losada (2004) - ! + ! implicit none - ! + ! real*4, intent(out) :: Cdterm real*4, intent(in) :: ahh ! [m] plant (total) height integer, intent(in) :: m - integer, intent(in) :: no_nodes ! number of unstructured grid nodes - integer, intent(in) :: no_secveg - real*4, intent(in) :: depth ! bed level, water depth - real*4, intent(in) :: H ! wave height + integer, intent(in) :: no_nodes ! number of unstructured grid nodes + integer, intent(in) :: no_secveg + real*4, intent(in) :: depth ! bed level, water depth + real*4, intent(in) :: H ! wave height real*4, intent(in) :: kwav ! wave number - real*4, intent(in) :: veg_bstems ! Width/diameter of individual vegetation stems [m] - real*4, intent(in) :: sigm ! [rad/s] mean frequency - ! + real*4, intent(in) :: veg_bstems ! Width/diameter of individual vegetation stems [m] + real*4, intent(in) :: sigm ! [rad/s] mean frequency + ! ! Local variables + ! real*4 :: pi ! 3.14159 - real*4 :: alfav ! [-] ratio between plant height and water depth + real*4 :: alfav ! [-] ratio between plant height and water depth real*4 :: um ! [m/s] typical velocity acting on the plant real*4 :: Tp ! [s] reference wave period real*4 :: KC ! [-] Keulegan-Carpenter number @@ -1501,12 +1939,14 @@ subroutine bulkdragcoeff(ahh, m, Cdterm, no_nodes, no_secveg, depth, H, kwav, ve integer :: myflag ! 1 => 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 @@ -1521,6 +1961,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? @@ -1530,19 +1971,196 @@ 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 + +subroutine momeqveg(no_nodes, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, depth, rho, H, Trep, unl, Fvw) + ! INput: no_nodes, no_secveg, veg_ah(k,:), veg_bstems(k,:), veg_Nstems(k,:), veg_Cd(k,:), depth(k), rho, H(k), Tp(k), unl(k,:), Fvw(k) + ! + implicit none + ! + ! Inputs + integer, intent(in) :: no_nodes, no_secveg + real*4, intent(in) :: depth ,rho, H, Trep + real*4, dimension(no_secveg), intent(in) :: veg_ah, veg_bstems, veg_Nstems, veg_Cd + real*4, dimension(50), intent(in) :: unl + ! + ! Output + real*4, intent(out) :: Fvw + ! + ! Local variables + integer :: m, t + real*4 :: dt, hvegeff, Fvgnlt, integral + real*4 :: Cd, b, N + ! + ! Initialize output force + ! + Fvw = 0.0 + ! + ! Time step within wave period + ! + dt = Trep / 50.0 + ! + ! Loop over vertical vegetation sections + do m = 1 , no_secveg + ! Effective submerged height of vegetation section + hvegeff = min(veg_ah(m), depth) + ! Read vegetation parameters + Cd = veg_Cd(m) + b = veg_bstems(m) + N = veg_Nstems(m) + ! Integrate vegetation drag over wave period using unl + integral = 0.0 + do t = 1, 50 !50=PPWL + integral = integral + (0.5 * Cd * b * N * hvegeff * unl(t) * abs(unl(t) ) ) * dt + enddo + ! Convert to force per unit mass and sum + Fvgnlt = -integral / depth / Trep !> units match with F(k) m/s2 + + Fvw = Fvw + Fvgnlt + enddo + ! +end subroutine momeqveg +subroutine swvegnonlin(no_nodes, kwav, depth, H, g, Trep, unl, etaw0) + use snapwave_RFtable + ! + implicit none + ! + integer, intent(in) :: no_nodes + real*4, dimension(no_nodes), intent(in) :: kwav + real*4, dimension(no_nodes), intent(in) :: depth + real*4, dimension(no_nodes), intent(in) :: H + real*4, intent(in) :: g + real*4, dimension(no_nodes), intent(in) :: Trep + real*4, dimension(no_nodes, 50),intent(out) :: unl + real*4, dimension(no_nodes, 50),intent(out) :: etaw0 + + real*4, dimension(:), save , allocatable :: h0, t0 + integer, save :: nh , nt ! save as it only needs to be done at first call + real*4, save :: dh , dt + real*4, dimension(50 ,8), save :: cs , sn ! MvdL: what is this fixed dimension 8 and 50? + + real*4, dimension(8) :: urf0 + real*4, dimension(50) :: urf2 , urf + real*4, dimension(50, 8) :: urf1 + + real*4, dimension(no_nodes) :: kmr , Urs , phi , w1 , w2 + real*4 :: p ,q , f0 , f1 , f2 , f3 + + integer :: k, irf, ih0, it0, jrf, ih1, it1 + ! + real*4 :: pi = 4.*atan(1.0) + + real*4, dimension(:,:,:), allocatable :: RFveg + ! + allocate(RFveg(11,18,20)) + ! + ! Based on Deltares' XBeach SurfBeat' subroutine: swvegnonlin + ! Subroutine to compute a net drag force due to wave skewness. Based on (matlab based) roller model with veggies by Ad. + ! + ! Background: + ! The drag force (Fveg) is a function of u*abs(u), which is zero for linear waves. For non-linear, skewed waves the + ! depth-averaged velocity integrated over the wave period is zero. However, due to the sharp peaks and flat troughs + ! the integral of u*abs(u) is non-zero, and can significantly reduce wave setup, or even lead to set-down (e.g. Dean & Bender,2006). + ! + ! Here we use a method based on Rienecker & Fenton (1981), similar to the method used for onshore sediment transport due to wave asymmetry/ + ! skewness (see also morphevolution.F90 + Van Thiel de Vries Phd thesis par 6.2.3). + ! + ! load Ad's RF-table (update for depth averaged velocities?) + call load_RFtable(RFveg) + ! + ! Initialize/Prepare for interpolation of RF-value from RFveg-table + if (.not. allocated(h0)) then + allocate(h0(no_nodes)) + allocate(t0(no_nodes)) + dh = 0.03 + dt = 1.25 + nh = floor(0.54/ dh) + nt = floor(25 / dt ) + !construct velocity profile based on cosine/sine functions / Fourier components + do irf =1 ,8 + do jrf =1 ,50 + cs ( jrf , irf ) = cos (( jrf * 2 * pi / 50) * irf ) + sn ( jrf , irf ) = sin (( jrf * 2 * pi / 50) * irf ) + enddo + enddo + endif + ! + h0 = min(nh * dh, max(dh, min(H, depth) / depth) ) + t0 = min(nt * dt, max(dt, Trep * sqrt (g / depth) ) ) + ! + ! Initialize + urf0 = 0 + urf1 = 0 + urf2 = 0 + urf = 0 + w1 = 0 + w2 = 0 + phi = 0 + Urs = 0 + kmr = 0 + ! + ! Now compute weight factors (w1,w2) for relative contribution of cosine and sine functions (for w1 = 1: only cosines -> + ! fully skewed Stokes wave, for w2 = 1: only sines -> fully asymmetric wave) based on Ruessink. + kmr = min(max(kwav, 0.01), 100.0) + Urs = H / (kmr * kmr * (depth **3) ) + + ! Compute phase and weight factors + phi = pi /2 * (1 - tanh (0.815/(Urs **0.672) ) ) + w1 = 1 - phi /( pi /2) + w2 = 1 - w1 + ! + ! Interpolate RieneckerFenton velocity from RFveg table from Ad + ! in ftab-dimension, only read 4:11 and sum later + do k =1, no_nodes + ! + ih0 = floor( h0(k) / dh) + it0 = floor( t0(k) / dt) + ih1 = min(ih0 + 1, nh) + it1 = min(it0 + 1, nt) + p = ( h0(k) - ih0 * dh) / dh + q = ( t0(k) - it0 * dt) / dt + f0 = (1 - p) * (1 - q) + f1 = p * (1 - q) + f2 = q * (1 - p) + f3 = p * q + ! + ! Compute velocity amplitude per component + do irf = 1, 8 + urf0(irf) = f0 * RFveg(irf + 3, ih0, it0) + f1 * RFveg(irf + 3, ih1, it0) + f2 * RFveg(irf+3, ih0, it1) + f3 * RFveg(irf + 3, ih1, it1) + enddo + ! fill velocity amplitude matrix urf1([50 time points, 8 components]) + do irf = 1, 8 + urf1(:, irf) = urf0(irf) + enddo + ! + ! Compute velocity profile matrix per component + urf1 = urf1 * (w1(k) * cs + w2(k) * sn ) + ! + ! Add velocity components + urf2 = sum(urf1, 2) + ! + ! Scale the results to get velocity profile over wave period + unl(k,:) = urf2 * sqrt(g * depth(k) ) + etaw0(k,:) = unl(k,:)*sqrt(max(depth(k),0.d0)/g) + enddo + ! +end subroutine swvegnonlin + end module snapwave_solver