diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj
index 34bf8f520..c4ff9cfdf 100644
--- a/source/sfincs_lib/sfincs_lib.vfproj
+++ b/source/sfincs_lib/sfincs_lib.vfproj
@@ -34,9 +34,12 @@
-
-
-
+
+
+
+
+
+
@@ -45,21 +48,31 @@
+
+
+
+
+
+
+
+
-
+
+
-
+
+
-
-
+
+
@@ -69,7 +82,8 @@
-
+
+
@@ -104,8 +118,11 @@
+
+
-
+
+
diff --git a/source/src/Makefile.am b/source/src/Makefile.am
index 891652f08..1b1b0fa00 100644
--- a/source/src/Makefile.am
+++ b/source/src/Makefile.am
@@ -19,13 +19,15 @@ 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 \
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 +52,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_input.f90 b/source/src/sfincs_input.f90
index efbad2e9a..25d228bd3 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
!
@@ -734,303 +735,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 ca3441bf4..42d4f602b 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+branch:318"
- build_date = "$Date: 2025-04-13"
+ build_revision = "$Rev: v2.3.2 mt. Faber+branch:snapwave_domain_updates_core"
+ build_date = "$Date: 2026-04-29"
!
call write_log('', 1)
call write_log('------------ Welcome to SFINCS ------------', 1)
diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90
index e3a1895e1..34f380549 100644
--- a/source/src/sfincs_ncoutput.F90
+++ b/source/src/sfincs_ncoutput.F90
@@ -3913,7 +3913,7 @@ 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, &
+ 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, 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,7 +4115,7 @@ 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_factor',waveforces_factor))
!
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 81d5015e0..4520995fd 100644
--- a/source/src/sfincs_snapwave.f90
+++ b/source/src/sfincs_snapwave.f90
@@ -604,6 +604,7 @@ subroutine read_snapwave_input()
! Reads snapwave data from sfincs.inp
!
use snapwave_data
+ use sfincs_read
!
implicit none
!
@@ -623,8 +624,9 @@ subroutine read_snapwave_input()
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_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)
@@ -635,12 +637,15 @@ subroutine read_snapwave_input()
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,1.0) ! underrelaxation factor for DoverA (set to 1.0 to disable)
+ call read_real_input(500,'snapwave_relax_factor_DoverE',relax_factor_DoverE,1.0) ! 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) ! 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_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) ! ! option controlling from what depth wave breaking should take place for IG waves, default baldock_ratio_ig=0.2
@@ -677,6 +682,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 +735,7 @@ subroutine read_snapwave_input()
restart = .true.
coupled_to_sfincs = .true.
!
- end subroutine
+ end subroutine read_snapwave_input
-
- 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_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/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90
index 1ab955526..e0946ddff 100644
--- a/source/src/snapwave/snapwave_data.f90
+++ b/source/src/snapwave/snapwave_data.f90
@@ -62,6 +62,7 @@ module snapwave_data
real*4, dimension(:), allocatable :: Hmx_ig
real*4, dimension(:,:), allocatable :: ee ! directional energy density
real*4, dimension(:,:), allocatable :: ee_ig ! directional infragravity energy density
+ real*4, dimension(:), allocatable :: DoverE
!
real*4, dimension(:,:), allocatable :: aa ! directional action density
real*4, dimension(:), allocatable :: sig ! mean frequency
@@ -172,10 +173,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 +216,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 +261,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
@@ -282,6 +285,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..746d91102 100644
--- a/source/src/snapwave/snapwave_domain.f90
+++ b/source/src/snapwave/snapwave_domain.f90
@@ -8,6 +8,7 @@ subroutine initialize_snapwave_domain()
!
use snapwave_data
use snapwave_boundaries
+ use snapwave_ncoutput
use interp
use sfincs_error
!
@@ -21,6 +22,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,6 +85,15 @@ subroutine initialize_snapwave_domain()
if (face_nodes(4,k)==0) face_nodes(4,k) = -999
enddo
!
+ ! write mesh to file
+ if (storesnapwavegrid) then
+ !
+ snapwave_ncfname = 'snapwavegrid.nc'
+ !
+ call write_snapwave_mesh(snapwave_ncfname, sferic == 1)
+ !
+ endif
+ !
! Done with the mesh
!
! keep on also if ja_vegetation==0, so array Dveg is initialized with zeroes
@@ -165,6 +176,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
!
@@ -206,7 +218,9 @@ subroutine initialize_snapwave_domain()
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 +317,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 +796,61 @@ subroutine boundaries(x,y,no_nodes,xb,yb,nb,tol,bndpts,nobndpts,bndindx,bndweigh
end subroutine boundaries
- subroutine neuboundaries_light(x,y,msk,no_nodes,tol,neumannconnected)
+ subroutine neuboundaries_light(x, y, msk, no_nodes, neumannconnected)
!
- ! TL: Based on subroutine find_nearest_depth_for_boundary_points of snapwave_boundaries.f90
- !
implicit none
!
integer, intent(in) :: no_nodes
- real*8, dimension(no_nodes), intent(in) :: x,y
+ real*8, dimension(no_nodes), intent(in) :: x, y
integer*1, dimension(no_nodes), intent(in) :: msk
- real*4, intent(in) :: tol
integer, dimension(no_nodes), intent(out) :: neumannconnected
!
- real*4 :: h1, h2, fac
- !
- real xgb, ygb, dst1, dst2, dst
- integer k, ib1, ib2, ic, kmin
- !
- ! Loop through all msk=3 cells
- !
- do ic = 1, no_nodes
- ! Loop through all grid points
- !
- if (msk(ic)==3) then ! point ic is on the neumann boundary
- !
- dst1 = tol
- dst2 = tol
- ib1 = 0
- ib2 = 0
- !
- do k = 1, no_nodes
- !
- if (msk(k)==1) then
- xgb = x(k)
- ygb = y(k)
- !
- dst = sqrt((x(ic) - xgb)**2 + (y(ic) - ygb)**2)
- !
- if (dst 0) .and. (ib2 > 0) ) then
- !
- ! Determine the index of the minimum value, if points found within 'tol' distance
- !
- if (dst1 < dst2) then
- kmin = ib1
- else
- kmin = ib2
- endif
- !
- neumannconnected(kmin)=ic
- !
- !write(*,*)kmin,ic
- !
- endif
- !
- endif
- enddo
+ real :: xgb, ygb, dst1, dst
+ integer :: k, ib1, ic
!
- end subroutine neuboundaries_light
-
-
-subroutine neuboundaries(x,y,no_nodes,xneu,yneu,n_neu,tol,neumannconnected)
+ ! Loop through all msk=3 cells and find their nearest msk=1 cell, save in 'neumannconnected'
!
- implicit none
- !
- 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()
@@ -1146,6 +1071,8 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
integer :: m
integer :: mu1
integer :: mu2
+ integer :: md1
+ integer :: md2
integer :: mnu1
integer :: nra
integer*1 :: mu
@@ -1166,7 +1093,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
!
@@ -1257,6 +1184,8 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
!
allocate(faces(4, 4*quadtree_nr_points)) ! max 4 nodes per faces, and max 4 faces per node
!
+ faces = 0
+ !
nfaces = 0
!
do ip = 1, quadtree_nr_points
@@ -1415,47 +1344,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 +1352,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 +1566,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
@@ -2035,6 +1914,79 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
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
!
! STEP 5 - set back snapwave_mask = 2&3 values
@@ -2098,7 +2050,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)
diff --git a/source/src/snapwave/snapwave_infragravity.f90 b/source/src/snapwave/snapwave_infragravity.f90
index 178f32843..bb2c7f590 100644
--- a/source/src/snapwave/snapwave_infragravity.f90
+++ b/source/src/snapwave/snapwave_infragravity.f90
@@ -43,21 +43,25 @@ 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
- !
+ ! 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..e5f8bf2c8 100644
--- a/source/src/snapwave/snapwave_solver.f90
+++ b/source/src/snapwave/snapwave_solver.f90
@@ -1,4 +1,4 @@
-module snapwave_solver
+module snapwave_solver
use sfincs_log
@@ -16,14 +16,12 @@ subroutine compute_wave_field()
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
!
integer :: itheta
integer :: k
!
- !allocate(sig(no_nodes))
allocate(sigm_ig(no_nodes))
!
g = 9.81
@@ -43,76 +41,96 @@ subroutine compute_wave_field()
!
ee_ig = waveps
!
- restart=1 !TODO TL: CHECK > we need this turned on right now for IG...
- !
endif
!
! Initialize wave period
!
do k = 1, no_nodes
+ !
if (inner(k)) then
+ !
Tp(k) = Tpini
+ !
endif
- if (neumannconnected(k)>0) then
- Tp(neumannconnected(k))=Tpini
+ !
+ if (neumannconnected(k) > 0) then
+ !
+ Tp(neumannconnected(k)) = Tpini
+ !
endif
+ !
enddo
!
! Compute celerities and refraction speed
!
Tp = max(tpmean_bwv,Tpini) ! to check voor windgroei
- sig = 2.0*pi/Tp
+ 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 ?
+ 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
+ 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
!
do k = 1, no_nodes
+ !
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)
+ !
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
+ do k = 1, no_nodes
+ !
+ Hmx_ig(k) = 0.88 / kwav_ig(k) * tanh(gamma_ig * kwav_ig(k) * depth(k) / 0.88) ! Note - uses gamma_ig
+ !Hmx_ig(k) = gamma_ig * depth(k)
+ !
enddo
- else
- Hmx_ig = 0.0
endif
!
do itheta = 1, ntheta
- ctheta(itheta,:) = sig/sinh(min(2.0*kwav*depth, 50.0))*(dhdx*sin(theta(itheta)) - dhdy*cos(theta(itheta)))
+ !
+ ctheta(itheta,:) = sig / sinh(min(2 * kwav * depth, 50.0)) * (dhdx * sin(theta(itheta)) - dhdy * cos(theta(itheta)))
+ !
enddo
!
if (igwaves) then
do itheta = 1, ntheta
- ctheta_ig(itheta,:) = sigm_ig/sinh(min(2.0*kwav_ig*depth, 50.0))*(dhdx*sin(theta(itheta)) - dhdy*cos(theta(itheta)))
+ !
+ ctheta_ig(itheta,:) = sigm_ig / sinh(min(2 * kwav_ig * depth, 50.0)) * (dhdx * sin(theta(itheta)) - dhdy * cos(theta(itheta)))
+ !
enddo
else
+ !
ctheta_ig = 0.0
+ !
endif
!
! Limit unrealistic refraction speed to 1/2 pi per wave period
!
do k = 1, no_nodes
- ctheta(:,k) = sign(1.0, ctheta(:,k))*min(abs(ctheta(:, k)), sig(k)/4)
+ ctheta(:,k) = sign(1.0, ctheta(:,k)) * min(abs(ctheta(:, k)), sig(k)/4)
enddo
!
if (igwaves) then
do k=1, no_nodes
- ctheta_ig(:,k) = sign(1.0, ctheta_ig(:,k))*min(abs(ctheta_ig(:, k)), sigm_ig(k)/4.0)
+ ctheta_ig(:,k) = sign(1.0, ctheta_ig(:,k)) * min(abs(ctheta_ig(:, k)), sigm_ig(k)/4)
enddo
endif
!
@@ -121,46 +139,48 @@ 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, &
+ igwaves, kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, &
beta, srcig, alphaig, Dw_ig, Df_ig, &
vegetation, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, Dveg, &
- zb, nwav, ig_opt, alpha_ig, gamma_ig, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig)
+ 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, &
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
@@ -219,7 +239,9 @@ 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)
!
! wind source vars
!
@@ -227,6 +249,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 ! ratio of mean wave dissipation over mean wave energy
real*4, dimension(ntheta,no_nodes), intent(out) :: WsorE, WsorA
real*4, dimension(no_nodes), intent(out) :: SwE, SwA
real*4, dimension(no_nodes), intent(inout) :: sig
@@ -264,7 +287,6 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
real*4, dimension(:), allocatable :: A,B,C,R ! coefficients in the tridiagonal matrix solved per point
real*4, dimension(:), allocatable :: B_aa,R_aa,aaprev ! coefficients in the tridiagonal matrix solved per point
real*4, dimension(:), allocatable :: A_ig,B_ig,C_ig,R_ig ! coefficients in the tridiagonal matrix solved per point
- real*4, dimension(:), allocatable :: DoverE ! ratio of mean wave dissipation over mean wave energy
real*4, dimension(:), allocatable :: DoverA ! ratio of mean wave dissipation over mean wave energy
real*4, dimension(:), allocatable :: DoverE_ig ! ratio of mean wave dissipation over mean wave energy
real*4, dimension(:), allocatable :: E ! mean wave energy
@@ -274,10 +296,9 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!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
@@ -292,26 +313,27 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
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
!
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
allocate(ok(no_nodes)); ok=0
allocate(indx(no_nodes,4)); indx=0
allocate(eeold(ntheta,no_nodes)); eeold=0.0
@@ -322,7 +344,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,7 +356,6 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
allocate(cgprev_ig(ntheta)); cgprev_ig=0.0
allocate(DoverE_ig(no_nodes)); DoverE_ig=0.0
allocate(E_ig(no_nodes)); E_ig=waveps
- !allocate(T_ig(no_nodes)); T_ig=0.0
allocate(sigm_ig(no_nodes)); sigm_ig=0.0
allocate(depthprev(ntheta,no_nodes)); depthprev=0.0
allocate(beta_local(ntheta,no_nodes)); beta_local=0.0
@@ -358,47 +378,51 @@ 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
+ !
+ if (dtheta < 0.0) dtheta = dtheta + 2 * pi
+ !
if (wind) then
- sig = 2*pi/Tpini
+ sig = 2 * pi / Tpini
else
- sig = 2*pi/Tp
+ sig = 2 * pi / Tp
endif
- oneoverdt = 1.0/dt
- oneover2dtheta = 1.0/2.0/dtheta
- rhog8 = 0.125*rho*g
+ oneoverdt = 1.0 / dt
+ oneover2dtheta = 1.0 / 2.0 / dtheta
+ rhog8 = 0.125 * rho * g
thetam = 0.0
- !H = 0.0 ! TODO - TL: CHeck > needed for restart for IG > set to 0 now in snapwave_domain.f90
- Dveg = 0.0
!
if (igwaves) then
- !T_ig = Tinc2ig*Tp
- sigm_ig = 2*pi/T_ig
+ !
+ sigm_ig = 2 * pi / T_ig
DoverE_ig = 0.0
+ !
endif
!
if (wind) then
+ !
DoverA = 0.0
ndissip = 3.0
WsorE = 0.0
WsorA = 0.0
- Ak = waveps/sigmax
+ Ak = waveps / sigmax
+ !
endif
!
! Sort coordinates in sweep directions
!
- 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
@@ -420,14 +444,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
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.
- !
- !exit
+ exit
!
endif
enddo
@@ -442,23 +459,23 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
if (.not.inner(k)) then
!
- 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)))
+ 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)))
!
!ee_ig(:, k) = eeinc2ig*ee(:,k) !TODO TL: determined in snapwave_boundaries.f90
!
if (igwaves) then
- E_ig(k) = sum(ee_ig(:, k))*dtheta
- H_ig(k) = sqrt(8*E_ig(k)/rho/g)
+ 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)
+ aa(:,k) = max(ee(:,k),waveps) / sig(k)
+ Ak = E(k) / sig(k)
endif
!
endif
@@ -472,7 +489,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
! 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)
+ call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local, gamma, gamma_fac_br)
!
! inout: alphaig_local, srcig_local - eeprev, eeprev_ig, cgprev, beta_local
! in: the rest
@@ -489,34 +506,39 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
if (wind) then
ee(:,k) = waveps
- sig(k) = 2*pi/Tpini
- aa(:,k) = ee(:,k)/sig(k)
+ sig(k) = 2 * pi / Tpini
+ aa(:,k) = ee(:,k) / sig(k)
else
ee(:,k) = waveps
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
+ Ek = sum(ee(:, k)) * dtheta
+ Hk = sqrt(Ek / rhog8) ! 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)
+ !
+ 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), baldock_exponent, Dwk, Hmx(k))
+ DoverE(k) = (Dwk + Dfk) / max(Ek, 1.0e-6)
+ !
endif
!
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)
+ 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), baldock_exponent, 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)
+ DwT = - c_dispT / (1.0 - ndissip) * (2. * pi) / sig(k)**2 * cg(k) * kwav(k) * DoverE(k)
+ DwAk = 0.5 / pi * (E(k) * DwT + 2.0 * pi * Ak * DoverE(k) )
+ DoverA(k) = DwAk / max(Ak,1e-6)
+ !
endif
!
endif
@@ -548,14 +570,14 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
if (inner(k)) then
!
- H(k) = sqrt(8*sum(ee(:, k))*dtheta/rho/g)
+ 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)
+ call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local, gamma, gamma_fac_br)
!
endif
!
@@ -570,7 +592,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
k=indx(count, sweep)
!
if (inner(k)) then
- if (depth(k)>1.1*hmin) then
+ if (depth(k) > 1.1 * hmin) then
!
if (ok(k) == 0) then
!
@@ -581,51 +603,51 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
k1 = prev(1, itheta, k)
k2 = prev(2, itheta, k)
!
- eeprev(itheta) = w(1, itheta, k)*ee(itheta, k1) + w(2, itheta, k)*ee(itheta, k2)
- cgprev(itheta) = w(1, itheta, k)*cg(k1) + w(2, itheta, k)*cg(k2)
+ eeprev(itheta) = w(1, itheta, k) * ee(itheta, k1) + w(2, itheta, k) * ee(itheta, k2)
+ cgprev(itheta) = w(1, itheta, k) * cg(k1) + w(2, itheta, k) * cg(k2)
!
if (igwaves) then
- eeprev_ig(itheta) = w(1, itheta, k)*ee_ig(itheta, k1) + w(2, itheta, k)*ee_ig(itheta, k2)
- cgprev_ig(itheta) = w(1, itheta, k)*cg_ig(k1) + w(2, itheta, k)*cg_ig(k2)
+ eeprev_ig(itheta) = w(1, itheta, k) * ee_ig(itheta, k1) + w(2, itheta, k) * ee_ig(itheta, k2)
+ cgprev_ig(itheta) = w(1, itheta, k) * cg_ig(k1) + w(2, itheta, k) * cg_ig(k2)
endif
!
if (wind) then
- aaprev(itheta) = w(1, itheta, k)*aa(itheta, k1) + w(2, itheta, k)*aa(itheta, k2)
+ aaprev(itheta) = w(1, itheta, k) * aa(itheta, k1) + w(2, itheta, k) * aa(itheta, k2)
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
+ depthlimfac = max(1.0, (sqrt(Ek / rhog8) / (gammax * depth(k)))**2.0)
+ Hk = sqrt(Ek / rhog8) !min(sqrt(Ek / rhog8), gamma * depth(k))
+ Ek = Ek / depthlimfac
!
if (wind) then
!
- Ak = sum(aaprev)*dtheta
+ Ak = sum(aaprev) * dtheta
!
Ak = Ak/depthlimfac
ee(:,k) = ee(:,k) / depthlimfac
aa(:,k) = aa(:,k) / depthlimfac
- sig(k) = Ek/Ak
+ 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
+ Ak = Ek / sig(k) ! to avoid small T in windinput
if (wind) then
- aaprev=min(aaprev,eeprev/sigmin)
- aaprev=max(aaprev,eeprev/sigmax)
+ 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))
+ uorbi = 0.5 * sig(k) * Hk / sinhkh(k)
+ Dfk = 0.28 * rho * fw(k) * uorbi**3
+ !
+ 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
Dwk = 0.
endif
@@ -636,7 +658,10 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
Dvegk = 0.
endif
!
- DoverE(k) = (Dwk + Dfk + Dvegk)/max(Ek, 1.0e-6)
+ !DoverE(k) = (Dwk + Dfk + Dvegk) / max(Ek, 1.0e-6)
+ ! Use some under relaxation (unless relax_factor_DoverE is set to 1.0)
+ !
+ DoverE(k) = (1.0 - relax_factor_DoverE) * DoverE(k) + relax_factor_DoverE * (Dwk + Dfk + Dvegk) / max(Ek, 1.0e-6)
!
if (wind) then
!
@@ -646,13 +671,13 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
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) )
+ 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)
+ DoverA(k) = DwAk / max(Ak,1e-6)
else
- DoverA(k) = (1.0-fac)*DoverA(k)+fac*DwAk/max(Ak,1.e-6)
+ DoverA(k) = (1.0 - relax_factor_DoverA) * DoverA(k) + relax_factor_DoverA * 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))
@@ -661,80 +686,84 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
do itheta = 1, ntheta
!
- R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta)*eeprev(itheta)/ds(itheta, k) - srcig_local(itheta, k) * shinc2ig
+ R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta) * eeprev(itheta) / ds(itheta, k) - srcig_local(itheta, k) * shinc2ig
!
enddo
!
do itheta = 2, ntheta - 1
!
- A(itheta) = -ctheta(itheta - 1, k)*oneover2dtheta
- B(itheta) = oneoverdt + cg(k)/ds(itheta,k) + DoverE(k)
- C(itheta) = ctheta(itheta + 1, k)*oneover2dtheta
+ A(itheta) = -ctheta(itheta - 1, k) * oneover2dtheta
+ B(itheta) = oneoverdt + cg(k) / ds(itheta,k) + DoverE(k)
+ C(itheta) = ctheta(itheta + 1, k) * oneover2dtheta
!
enddo
!
- A(1) = -ctheta(ntheta, k)*oneover2dtheta
- B(1) = oneoverdt + cg(k)/ds(1,k) + DoverE(k)
- C(1) = ctheta(2, k)*oneover2dtheta
+ A(1) = -ctheta(ntheta, k) * oneover2dtheta
+ B(1) = oneoverdt + cg(k) / ds(1,k) + DoverE(k)
+ C(1) = ctheta(2, k) * oneover2dtheta
!
- A(ntheta) = -ctheta(ntheta - 1, k)*oneover2dtheta
- B(ntheta) = oneoverdt + cg(k)/ds(ntheta,k) + DoverE(k)
- C(ntheta) = ctheta(1, k)*oneover2dtheta
+ A(ntheta) = -ctheta(ntheta - 1, k) * oneover2dtheta
+ B(ntheta) = oneoverdt + cg(k) / ds(ntheta,k) + DoverE(k)
+ C(ntheta) = ctheta(1, k) * oneover2dtheta
!
! Solve tridiagonal system per point
!
if (wind) then
do itheta = 2, ntheta - 1
- B_aa(itheta) = oneoverdt + cg(k)/ds(itheta,k) + DoverA(k)
- R_aa(itheta) = (oneoverdt)*aa(itheta, k) + cgprev(itheta)*aaprev(itheta)/ds(itheta, k)
+ B_aa(itheta) = oneoverdt + cg(k) / ds(itheta,k) + DoverA(k)
+ R_aa(itheta) = (oneoverdt) * aa(itheta, k) + cgprev(itheta) * aaprev(itheta) / ds(itheta, k)
enddo
!
if (ctheta(1,k)<0) then
- B_aa(1) = oneoverdt - ctheta(1, k)/dtheta + cg(k)/ds(1, k) + DoverA(k)
- R_aa(1) = (oneoverdt)*aa(1, k) + cgprev(1)*aaprev(1)/ds(1, k)
+ B_aa(1) = oneoverdt - ctheta(1, k) / dtheta + cg(k) / ds(1, k) + DoverA(k)
+ R_aa(1) = (oneoverdt) * aa(1, k) + cgprev(1) * aaprev(1) / ds(1, k)
else
- B_aa(1) = oneoverdt + cg(k)/ds(1, k) + DoverA(k)
- R_aa(1) = (oneoverdt)*aa(1, k) + cgprev(1)*aaprev(1)/ds(1, k)
+ B_aa(1) = oneoverdt + cg(k) / ds(1, k) + DoverA(k)
+ R_aa(1) = (oneoverdt) * aa(1, k) + cgprev(1) * aaprev(1) / ds(1, k)
endif
!
if (ctheta(ntheta, k)>0) then
- B_aa(ntheta) = oneoverdt + ctheta(ntheta, k)/dtheta + cg(k)/ds(ntheta, k) + DoverA(k)
- R_aa(ntheta) = (oneoverdt )*aa(ntheta,k) + cgprev(ntheta)*aaprev(ntheta)/ds(ntheta, k)
+ B_aa(ntheta) = oneoverdt + ctheta(ntheta, k) / dtheta + cg(k) / ds(ntheta, k) + DoverA(k)
+ R_aa(ntheta) = (oneoverdt) * aa(ntheta,k) + cgprev(ntheta) * aaprev(ntheta) / ds(ntheta, k)
else
- B_aa(ntheta) = oneoverdt + cg(k)/ds(ntheta, k) + DoverA(k)
- R_aa(ntheta) = (oneoverdt)*aa(ntheta,k) + cgprev(ntheta)*aaprev(ntheta)/ds(ntheta, k)
+ B_aa(ntheta) = oneoverdt + cg(k) / ds(ntheta, k) + DoverA(k)
+ R_aa(ntheta) = (oneoverdt) * aa(ntheta,k) + cgprev(ntheta) * aaprev(ntheta) / ds(ntheta, k)
endif
R(:) = R(:) + WsorE(:,k)
R_aa(:) = R_aa(:) + WsorA(:,k)
!
call solve_tridiag(A, B, C, R, ee(:,k), ntheta)
call solve_tridiag(A,B_aa,C,R_aa,aa(:,k),ntheta)
+ !
ee(:, k) = max(ee(:, k), waveps)
- aa(:,k) = max(aa(:,k),waveps/sigmax)
- aa(:,k) = max(aa(:,k),waveps/sig(k))
+ 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
+ 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
+ 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) = 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
+ 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
@@ -742,63 +771,67 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
! 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
+ 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
+ Dfk_ig = fw_ig(k) * 0.0361 * (9.81 / depth(k))**1.5 * Hk * Ek_ig
!
! Dissipation of infragravity waves
!
- if (Hk_ig>baldock_ratio_ig*Hmx_ig(k)) then
- call baldock(rho, g, alfa_ig, gamma_ig, depth(k), Hk_ig, T_ig(k), 1, Dwk_ig, Hmx_ig(k))
+ if (Hk_ig > baldock_ratio_ig * Hmx_ig(k)) then
+ !
+ call baldock(rho, g, alfa_ig, gamma_ig, depth(k), Hk_ig, T_ig(k), baldock_exponent, 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?
+ DoverE_ig(k) = (Dwk_ig + Dfk_ig) / max(Ek_ig, 1.0e-6) ! Not using underrelaxation for IG dissipation for now, but we could add this if needed (relax_factor_DoverE_ig)
!
do itheta = 1, ntheta
!
- R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta)*eeprev_ig(itheta)/ds(itheta, k) + srcig_local(itheta, k) !TL: new version
+ R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta) * eeprev_ig(itheta) / ds(itheta, k) + srcig_local(itheta, k) !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
+ 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
+ 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)
+ 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)
+ A_ig(ntheta) = -ctheta_ig(ntheta - 1, k) / dtheta
+ B_ig(ntheta) = oneoverdt + ctheta_ig(ntheta, k) / dtheta + cg_ig(k) / ds(ntheta, k) + DoverE_ig(k)
C_ig(ntheta) = 0.0
else
A_ig(ntheta) = 0.0
- B_ig(ntheta) = oneoverdt + cg_ig(k)/ds(ntheta, k) + DoverE_ig(k)
+ B_ig(ntheta) = oneoverdt + cg_ig(k) / ds(ntheta, k) + DoverE_ig(k)
C_ig(ntheta) = 0.0
endif
!
! Solve tridiagonal system per point
!
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
@@ -813,7 +846,9 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
!
ee(:, k) = 0.0
if (wind) then
+ !
aa(:,k) = 0.0
+ !
endif
ee_ig(:, k) = 0.0
!
@@ -822,7 +857,9 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
endif
!
if (neumannconnected(k)/=0) then
- kn = neumannconnected(k)
+ !
+ kn = neumannconnected(k) ! Index of internal point
+ !
sinhkh(kn) = sinhkh(k)
kwav(kn) = kwav(k)
Hmx(kn) = Hmx(k)
@@ -830,15 +867,20 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
ee_ig(:, kn) = ee_ig(:, k) ! TL: Added Neumann option for IG
ctheta(:, kn) = ctheta(:, k)
cg(kn) = cg(k)
+ !
if (wind) then
+ !
sig(kn) = sig(k)
- Tp(kn) = 2.0*pi/sig(k)
+ Tp(kn) = 2.0 * pi / sig(k)
WsorE(:,kn) = WsorE(:,k)
WsorA(:,kn) = WsorA(:,k)
aa(:,kn) = aa(:,k)
+ !
endif
+ !
Df(kn) = Df(k)
Dw(kn) = Dw(k)
+ !
endif
!
enddo
@@ -852,7 +894,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
dee = ee(:, k) - eeold(:, k)
diff(k) = maxval(abs(dee))
!
- if (diff(k)/eemax 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)
+ !
+ ! Limit srcig to 0 after waves start (significantly) breaking, as defined here as gam=Hrms,inc / h > (gamma_fac_br * gamma)
+ ! 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
+ !
+ ! 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
!
else ! TL: option to add future parameterisations here for e.g. coral reef type coasts