From 5a7b13c068f0eb6133cdd2b8e1c16dbd47ca1f55 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 25 Mar 2026 10:19:17 +0100 Subject: [PATCH 01/25] - As done in PR of jre/snapwave_sfincs, use the SFINCS input routines directly in SnapWave - Made clearer by putting the read functions in new filed 'sfincs_read', and only call those in sfincs_snapwave.f90 --- source/sfincs_lib/sfincs_lib.vfproj | 14 +- source/src/Makefile.am | 5 +- source/src/sfincs_input.f90 | 302 +-------------------------- source/src/sfincs_read.f90 | 303 ++++++++++++++++++++++++++++ source/src/sfincs_snapwave.f90 | 300 +-------------------------- 5 files changed, 320 insertions(+), 604 deletions(-) create mode 100644 source/src/sfincs_read.f90 diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj index 34bf8f520..d73e3675e 100644 --- a/source/sfincs_lib/sfincs_lib.vfproj +++ b/source/sfincs_lib/sfincs_lib.vfproj @@ -34,7 +34,8 @@ - + + @@ -50,14 +51,16 @@ - + + - + + @@ -104,8 +107,11 @@ + + - + + diff --git a/source/src/Makefile.am b/source/src/Makefile.am index 891652f08..f2f427d7d 100644 --- a/source/src/Makefile.am +++ b/source/src/Makefile.am @@ -19,7 +19,8 @@ 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 \ @@ -50,7 +51,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 ba60f5eaa..78b7f5716 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 ! @@ -681,303 +682,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_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..7f400e855 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 ! @@ -729,304 +730,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 From 10dfa1ad19a2fcf637379cba7a8a9f7f8f14b526 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 25 Mar 2026 10:19:17 +0100 Subject: [PATCH 02/25] - As done in PR of jre/snapwave_sfincs, use the SFINCS input routines directly in SnapWave - Made clearer by putting the read functions in new filed 'sfincs_read', and only call those in sfincs_snapwave.f90 --- source/sfincs_lib/sfincs_lib.vfproj | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj index d73e3675e..d9ed21317 100644 --- a/source/sfincs_lib/sfincs_lib.vfproj +++ b/source/sfincs_lib/sfincs_lib.vfproj @@ -61,8 +61,7 @@ - - + From faddb923bf5cada3e676928c7d88f73a8fbaa953 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 25 Mar 2026 10:20:00 +0100 Subject: [PATCH 03/25] - Also in sfincs_spiderweb.f90 replace local read_real_input and read_int_input by use sfincs_read --- source/src/sfincs_spiderweb.f90 | 56 +-------------------------------- 1 file changed, 1 insertion(+), 55 deletions(-) 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) From 16e01ef9ce96def6716983763f8afa8b0a836474 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 25 Mar 2026 10:58:49 +0100 Subject: [PATCH 04/25] - add posibility to inspect snapwave grid as in PR jre/snapwave_sfincs --- source/sfincs_lib/sfincs_lib.vfproj | 5 +- source/src/Makefile.am | 3 +- source/src/sfincs_snapwave.f90 | 1 + source/src/snapwave/snapwave_data.f90 | 1 + source/src/snapwave/snapwave_domain.f90 | 11 ++ source/src/snapwave/snapwave_ncoutput.F90 | 140 ++++++++++++++++++++++ 6 files changed, 159 insertions(+), 2 deletions(-) create mode 100644 source/src/snapwave/snapwave_ncoutput.F90 diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj index d9ed21317..fdf4dac89 100644 --- a/source/sfincs_lib/sfincs_lib.vfproj +++ b/source/sfincs_lib/sfincs_lib.vfproj @@ -36,7 +36,8 @@ - + + @@ -46,6 +47,8 @@ + + diff --git a/source/src/Makefile.am b/source/src/Makefile.am index f2f427d7d..1b1b0fa00 100644 --- a/source/src/Makefile.am +++ b/source/src/Makefile.am @@ -26,7 +26,8 @@ libsfincs_la_SOURCES = \ 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 \ diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index 7f400e855..70d8faf99 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -678,6 +678,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) diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90 index 1ab955526..0c1bfaf44 100644 --- a/source/src/snapwave/snapwave_data.f90 +++ b/source/src/snapwave/snapwave_data.f90 @@ -282,6 +282,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..3e70f33eb 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 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 From 161c8a7891994d3fc69994763df4b979c6b3bade Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 25 Mar 2026 11:00:26 +0100 Subject: [PATCH 05/25] - restart already fixed to true in sfincs_snapwave --- source/src/snapwave/snapwave_solver.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 53cc31d77..b74cbbcd4 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -43,8 +43,6 @@ 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 From ba0148fdc3f7b8155090c05e9fe905dc16b209b8 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 25 Mar 2026 11:12:56 +0100 Subject: [PATCH 06/25] - Remove unused subroutine 'neuboundaries' that is replaced by neuboundaries_light --- source/src/snapwave/snapwave_domain.f90 | 52 ------------------------- 1 file changed, 52 deletions(-) diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90 index 3e70f33eb..38159043f 100644 --- a/source/src/snapwave/snapwave_domain.f90 +++ b/source/src/snapwave/snapwave_domain.f90 @@ -886,58 +886,6 @@ subroutine neuboundaries_light(x,y,msk,no_nodes,tol,neumannconnected) end subroutine neuboundaries_light -subroutine neuboundaries(x,y,no_nodes,xneu,yneu,n_neu,tol,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 - enddo - endif - enddo - ! -end subroutine neuboundaries - subroutine read_snapwave_sfincs_mesh() ! From d5559a9487f3eff6d4a1f74c6f3b0de8ce7bfe72 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 25 Mar 2026 13:15:37 +0100 Subject: [PATCH 07/25] - Preprocess snapwave_ncoutput --- source/sfincs_lib/sfincs_lib.vfproj | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj index fdf4dac89..c4ff9cfdf 100644 --- a/source/sfincs_lib/sfincs_lib.vfproj +++ b/source/sfincs_lib/sfincs_lib.vfproj @@ -38,7 +38,8 @@ - + + @@ -48,10 +49,17 @@ + + + + + + - + + @@ -74,7 +82,8 @@ - + + From 7bedcc17c116dc90ea6845ca5617911a70b2a513 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 26 Mar 2026 15:51:00 +0100 Subject: [PATCH 08/25] - Distinguish this version from the full PR snapwave_domain_updates - Goal is to here add the good stuff to main, without all the reorganisation that killed proper wave growth - All of my relevant changes are already cherry-picked, now onto Maarten's stuff --- source/src/sfincs_lib.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index 5a012f34f..aaaf30ad5 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: 2026-03-20" + build_revision = "$Rev: v2.3.2 mt. Faber+branch:snapwave_domain_updates_core" + build_date = "$Date: 2026-03-26" ! call write_log('', 1) call write_log('------------ Welcome to SFINCS ------------', 1) From 503d6db01fc7b03420f640e15934ebe6849e47a4 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 26 Mar 2026 16:23:32 +0100 Subject: [PATCH 09/25] - Bring over from snapwave_domain_changes of declaring DoverE before snapwave_solver - Bring over update neumann_boundaries_light --- source/src/snapwave/snapwave_data.f90 | 6 +- source/src/snapwave/snapwave_domain.f90 | 139 +++++++++--------------- 2 files changed, 54 insertions(+), 91 deletions(-) diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90 index 0c1bfaf44..4fbb99db5 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 diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90 index 38159043f..551353ad8 100644 --- a/source/src/snapwave/snapwave_domain.f90 +++ b/source/src/snapwave/snapwave_domain.f90 @@ -176,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 ! @@ -217,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. @@ -314,26 +317,9 @@ subroutine initialize_snapwave_domain() if (any(msk == 3)) then ! ! We already have all msk=3 Neumann points, now find each their nearest cell 'neumannconnected' using new 'neuboundaries_light' - call neuboundaries_light(x,y,msk,no_nodes,tol,neumannconnected) - ! - if (ANY(neumannconnected > 0)) then - ! - write(logstr,*)'SnapWave: Neumann connected boundaries found ...' - call write_log(logstr, 0) - ! - do k=1,no_nodes - if (neumannconnected(k)>0) then - if (msk(k)==1) then - ! k is inner and can be neumannconnected - inner(neumannconnected(k))= .false. - msk(neumannconnected(k)) = 3 !TL: should already by 3, but left it like in SnapWave SVN - else - ! we don't allow neumannconnected links if the node is an open boundary - neumannconnected(k) = 0 - endif - endif - enddo - endif + ! + call neuboundaries_light(x, y, msk, no_nodes, neumannconnected) + ! else ! neumannconnected = 0 @@ -807,81 +793,58 @@ subroutine boundaries(x,y,no_nodes,xb,yb,nb,tol,bndpts,nobndpts,bndindx,bndweigh end subroutine boundaries - subroutine neuboundaries_light(x,y,msk,no_nodes,tol,neumannconnected) + subroutine neuboundaries_light(x, y, msk, no_nodes, neumannconnected) ! - ! TL: Based on subroutine find_nearest_depth_for_boundary_points of snapwave_boundaries.f90 - ! implicit none ! integer, intent(in) :: no_nodes - real*8, dimension(no_nodes), intent(in) :: x,y + real*8, dimension(no_nodes), intent(in) :: x, y integer*1, dimension(no_nodes), intent(in) :: msk - real*4, intent(in) :: tol integer, dimension(no_nodes), intent(out) :: neumannconnected ! - real*4 :: h1, h2, fac - ! - real xgb, ygb, dst1, dst2, dst - integer k, ib1, ib2, ic, kmin - ! - ! Loop through all msk=3 cells - ! - do ic = 1, no_nodes - ! Loop through all grid points - ! - if (msk(ic)==3) then ! point ic is on the neumann boundary - ! - dst1 = tol - dst2 = tol - ib1 = 0 - ib2 = 0 - ! - do k = 1, no_nodes - ! - if (msk(k)==1) then - xgb = x(k) - ygb = y(k) - ! - dst = sqrt((x(ic) - xgb)**2 + (y(ic) - ygb)**2) - ! - if (dst 0) then ! - if ( (ib1 > 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 + neumannconnected(ic) = ib1 + ! + endif + ! + endif + ! + enddo ! end subroutine neuboundaries_light From 9a6c4920064b1aec1d74966109e5d548d0594b49 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 26 Mar 2026 16:29:08 +0100 Subject: [PATCH 10/25] - Bring over 'Add triangles around stair case boundaries" --- source/src/snapwave/snapwave_domain.f90 | 130 ++++++++++++++---------- 1 file changed, 78 insertions(+), 52 deletions(-) diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90 index 551353ad8..cd9388875 100644 --- a/source/src/snapwave/snapwave_domain.f90 +++ b/source/src/snapwave/snapwave_domain.f90 @@ -1068,6 +1068,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 @@ -1088,7 +1090,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 ! @@ -1179,6 +1181,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 @@ -1337,47 +1341,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 @@ -1386,7 +1349,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 @@ -1601,14 +1563,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 @@ -1957,6 +1911,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 @@ -2020,7 +2047,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) From bc35896c5c9be05b1d7517753c5b4c28fe555eb7 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 26 Mar 2026 16:33:13 +0100 Subject: [PATCH 11/25] - Last changes besides snapwave_solver.f90 ones I think --- source/src/sfincs_ncoutput.F90 | 4 ++-- source/src/sfincs_snapwave.f90 | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index e3a789070..1c90d4668 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 @@ -4100,7 +4100,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)) ! ! SnapWave IG diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index 70d8faf99..63ecd0ecf 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -624,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) From 132813acb3307ae6994479c7a7ad6d9aede5f515 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 26 Mar 2026 16:34:37 +0100 Subject: [PATCH 12/25] - Add 200m limitation, and dimensionless not > todo --- source/src/snapwave/snapwave_infragravity.f90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/source/src/snapwave/snapwave_infragravity.f90 b/source/src/snapwave/snapwave_infragravity.f90 index 178f32843..502cfb3c1 100644 --- a/source/src/snapwave/snapwave_infragravity.f90 +++ b/source/src/snapwave/snapwave_infragravity.f90 @@ -46,6 +46,8 @@ subroutine determine_ig_bc(x_bwv, y_bwv, hsinc, tpinc, ds, jonswapgam, depth, Ti ! 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 + ! + ! MvO - Is this really the best check? Does not seem very non-dimensional. Maybe something like hsig / depth > 0.5 ? ! write(logstr,*)'ERROR SnapWave - depth at boundary input point ',x_bwv, y_bwv,' dropped below 5 m: ',depth call write_log(logstr, 1) @@ -58,6 +60,12 @@ subroutine determine_ig_bc(x_bwv, y_bwv, hsinc, tpinc, ds, jonswapgam, depth, Ti ! depth = 5.0 ! + 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] From e66aba0d1258471609a259d9dcf5c3e84d0aa6cf Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 26 Mar 2026 17:16:10 +0100 Subject: [PATCH 13/25] - Do many of the visual edits --- source/src/snapwave/snapwave_solver.f90 | 481 +++++++++++++----------- 1 file changed, 261 insertions(+), 220 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index b74cbbcd4..435a42f79 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -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 @@ -48,69 +46,91 @@ subroutine compute_wave_field() ! 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 = sig**2 / g * (1.0 - exp(expon))**(-0.4) + ! else + ! cg_ig = 0.0 kwav_ig = 0.0 + ! endif ! do k = 1, no_nodes + ! sinhkh(k) = sinh(min(kwav(k)*depth(k), 50.0)) + !Hmx(k) = 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 ! @@ -127,18 +147,18 @@ subroutine compute_wave_field() wind, & H,Dw,F,Df,thetam,sinhkh,& Hmx, ee, windspreadfac, u10, niter, crit, & - hmin, baldock_ratio, baldock_ratio_ig, & + hmin, baldock_ratio, baldock_ratio_ig, baldock_exponent, & aa, sig, jadcgdx, sigmin, sigmax,& - c_dispT, WsorE, WsorA, SwE, SwA, Tpini, & - igwaves,kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & + c_dispT, DoverE, WsorE, WsorA, SwE, SwA, Tpini, & + igwaves, kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & beta, srcig, alphaig, Dw_ig, Df_ig, & vegetation, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, Dveg, & zb, nwav, ig_opt, alpha_ig, gamma_ig, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) ! call timer(t3) ! - Fx = F*cos(thetam) - Fy = F*sin(thetam) + Fx = F * cos(thetam) + Fy = F * sin(thetam) ! end subroutine @@ -152,10 +172,10 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & wind, & H,Dw,F,Df,thetam,sinhkh,& Hmx, ee, windspreadfac, u10, niter, crit, & - hmin, baldock_ratio, baldock_ratio_ig, & + hmin, baldock_ratio, baldock_ratio_ig, baldock_exponent, & aa, sig, jadcgdx, sigmin, sigmax,& - c_dispT, WsorE, WsorA, SwE, SwA, Tpini, & - igwaves,kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & + c_dispT, DoverE, WsorE, WsorA, SwE, SwA, Tpini, & + igwaves, kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & 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) @@ -225,6 +245,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 @@ -262,7 +283,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 @@ -272,8 +292,8 @@ 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 @@ -295,21 +315,21 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & 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 @@ -320,7 +340,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 ! @@ -333,7 +352,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 @@ -356,47 +374,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 @@ -418,14 +440,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 @@ -440,23 +455,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 @@ -487,34 +502,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 = 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 + ! + 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) + 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), 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) + 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 @@ -546,7 +566,7 @@ 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 @@ -568,7 +588,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 ! @@ -579,49 +599,49 @@ 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 = 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 ! + 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) , 1, Dwk, Hmx(k)) else @@ -634,7 +654,7 @@ 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) ! if (wind) then ! @@ -644,13 +664,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 - 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)) @@ -659,80 +679,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 @@ -740,63 +764,68 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! IG ! if (igwaves) then - Ek_ig = sum(eeprev_ig)*dtheta + 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 + 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 + 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) = (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 + 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 @@ -811,7 +840,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 ! @@ -820,7 +851,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) @@ -828,15 +861,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 @@ -850,7 +888,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 Date: Thu, 26 Mar 2026 17:22:13 +0100 Subject: [PATCH 15/25] - Updated baldock --- source/src/snapwave/snapwave_solver.f90 | 47 ++++++++++++++++--------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 18df0e516..0d8ddd7ef 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -519,7 +519,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! 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)) + 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 @@ -529,7 +529,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & 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)) + 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 @@ -645,7 +645,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & 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) , 1, Dwk, Hmx(k)) + call baldock(rho, g, alfa, gamma, depth(k), Hk, 2*pi/sig(k), baldock_exponent, Dwk, Hmx(k)) else Dwk = 0. endif @@ -782,7 +782,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! 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)) + call baldock(rho, g, alfa_ig, gamma_ig, depth(k), Hk_ig, T_ig(k), baldock_exponent, Dwk_ig, Hmx_ig(k)) ! else ! @@ -945,10 +945,10 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & Df(k) = 0.28 * rho * fw(k) * uorbi**3 ! if (wind) then - call baldock(rho, g, alfa, gamma, depth(k), H(k), 2.0*pi/sig(k), 1, Dw(k), Hmx(k)) + call baldock(rho, g, alfa, gamma, depth(k), H(k), 2.0*pi/sig(k), baldock_exponent, Dw(k), Hmx(k)) F(k) = Dw(k) * kwav(k) / sig(k) / rho / depth(k) else - call baldock(rho, g, alfa, gamma, depth(k), H(k), Tp(k), 1, Dw(k), Hmx(k)) + call baldock(rho, g, alfa, gamma, depth(k), H(k), Tp(k), baldock_exponent, Dw(k), Hmx(k)) F(k) = Dw(k) * kwav(k) / sig(k) / rho / depth(k) !F(k) = (Dw(k) + Df(k))*kwav(k)/sig(k)/rho/depth(k) !F(k) = (Dw(k) + Df(k))*kwav(k)/sigm ! TODO TL: before was this, now multiplied with rho*depth(k) in sfincs_snapwave.f90 @@ -969,7 +969,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! Df_ig(k) = fw_ig(k) * 0.0361 * (9.81 / depth(k))**1.5 * H(k) * E_ig(k) !org ! - call baldock(rho, g, alfa_ig, gamma_ig, depth(k), H_ig(k), T_ig(k), 1, Dw_ig(k), Hmx_ig(k)) + call baldock(rho, g, alfa_ig, gamma_ig, depth(k), H_ig(k), T_ig(k), baldock_exponent, Dw_ig(k), Hmx_ig(k)) ! ! average beta, alphaig, srcig over directions betamean(k) = sum(beta_local(:,k)) / ntheta ! real mean @@ -1034,7 +1034,7 @@ subroutine solve_tridiag(a, b, c, d, x, n) ! end subroutine solve_tridiag - subroutine baldock (rho,g,alfa,gamma,depth,H,T,opt,Dw,Hmax) + subroutine baldock (rho, g, alfa, gamma, depth, H, T, iexp, Dw, Hmax) ! real*4, intent(in) :: rho real*4, intent(in) :: g @@ -1043,22 +1043,37 @@ subroutine baldock (rho,g,alfa,gamma,depth,H,T,opt,Dw,Hmax) real*4, intent(in) :: depth real*4, intent(in) :: H real*4, intent(in) :: T - integer, intent(in) :: opt + integer, intent(in) :: iexp real*4, intent(out) :: Dw real*4, intent(in) :: Hmax real*4 :: Hloc + real*4 :: f ! ! Compute dissipation according to Baldock ! - Hloc=max(H,1.e-6) - if (opt==1) then - Dw=0.28*alfa*rho*g/T*exp(-(Hmax/Hloc)**2)*(Hmax**2+Hloc**2) - !Dw=0.25*alfa*rho*g/T*exp(-(Hmax/Hloc)**2)*(Hmax**2+Hloc**2) !Old - else - Dw=0.28*alfa*rho*g/T*exp(-(Hmax/Hloc)**2)*(Hmax**3+Hloc**3)/gamma/depth - !Dw=0.25*alfa*rho*g/T*exp(-(Hmax/Hloc)**2)*(Hmax**3+Hloc**3)/gamma/depth !Old + Hloc = max(H, 1.e-6) + ! + if (iexp > 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) From a88f431a9268451c64f366a82ed0ffe2bb2c94dc Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 26 Mar 2026 17:26:28 +0100 Subject: [PATCH 16/25] - Don't limit Hk and Hkig with gamma*depth anymore --- source/src/snapwave/snapwave_solver.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 0d8ddd7ef..9e98a5def 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -512,7 +512,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! ! Make sure DoverE is filled based on previous ee Ek = sum(ee(:, k)) * dtheta - Hk = min(sqrt(Ek / rhog8), gamma * depth(k)) + Hk = sqrt(Ek / rhog8) ! min(sqrt(Ek / rhog8), gamma * depth(k)) Ek = rhog8 * Hk**2 ! if (.not. wind) then @@ -618,7 +618,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & 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)) + Hk = sqrt(Ek / rhog8) !min(sqrt(Ek / rhog8), gamma * depth(k)) Ek = Ek / depthlimfac ! if (wind) then @@ -644,7 +644,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & uorbi = 0.5 * sig(k) * Hk / sinhkh(k) Dfk = 0.28 * rho * fw(k) * uorbi**3 ! - if (Hk>baldock_ratio*Hmx(k)) then + 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. @@ -770,8 +770,8 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! 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? + 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 @@ -938,7 +938,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & thetam(k) = atan2(sum(ee(:, k) * sin(theta)), sum(ee(:, k) * cos(theta))) ! if (wind) then - uorbi = 0.5 * sig(k) * Hk / sinhkh(k) + uorbi = 0.5 * sig(k) * H(k) / sinhkh(k) else uorbi = pi * H(k) / Tp(k) / sinhkh(k) !! to check if we want array endif From 111f78217d1a1aedaee14d56693e8e78be01aa31 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 2 Apr 2026 14:02:37 +0200 Subject: [PATCH 17/25] - Make Herbers IG boundary warning dimensionless, using hsinc / depth - Check further for best ratio later, start at 0.5 --- source/src/snapwave/snapwave_infragravity.f90 | 32 ++++++++----------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/source/src/snapwave/snapwave_infragravity.f90 b/source/src/snapwave/snapwave_infragravity.f90 index 502cfb3c1..bb2c7f590 100644 --- a/source/src/snapwave/snapwave_infragravity.f90 +++ b/source/src/snapwave/snapwave_infragravity.f90 @@ -43,23 +43,19 @@ 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 + ! 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 (depth < 5.0) then - ! - ! MvO - Is this really the best check? Does not seem very non-dimensional. Maybe something like hsig / depth > 0.5 ? - ! - write(logstr,*)'ERROR SnapWave - depth at boundary input point ',x_bwv, y_bwv,' 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 - ! + 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? @@ -77,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) From 6f81cda703febc5df3522ff6f993bc2120fd4700 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 2 Apr 2026 16:35:25 +0200 Subject: [PATCH 18/25] - Already do the bugfix of PR #283 by switching to srcig based on E and redistribute by /E * ee --- source/src/snapwave/snapwave_solver.f90 | 125 +++++++++++++++++------- 1 file changed, 89 insertions(+), 36 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 9e98a5def..05de09b1f 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -487,7 +487,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, 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 @@ -575,7 +575,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! ! 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, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) ! endif ! @@ -1076,7 +1076,15 @@ subroutine baldock (rho, g, alfa, gamma, depth, H, T, iexp, Dw, Hmax) ! 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, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) + ! + ! 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 ! @@ -1095,6 +1103,7 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d 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 ! ! Inout variables real*4, dimension(:,:), intent(inout) :: alphaig_local ! Local infragravity wave shoaling parameter alpha @@ -1107,20 +1116,60 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d integer :: itheta ! directional counter 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 + ! + ! 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 + ! 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 + ! + ! 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 @@ -1152,10 +1201,10 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d 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 ! - Sxxprev(itheta) = w(1, itheta, k)*Sxx(itheta,k1) + w(2, itheta, k)*Sxx(itheta,k2) + Sxxprev(itheta) = w(1, itheta, k) * Sxx(k1) + w(2, itheta, k) * Sxx(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) + 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) ! @@ -1179,26 +1228,30 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d 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) 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 = 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 - ! - dSxx = Sxx_cons - Sxxprev(itheta) - ! - elseif (ig_opt == 2) then ! Option taking actual difference for dSxx/dx - ! - dSxx = Sxx(itheta,k) - Sxxprev(itheta) - ! - 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) - ! + ! + if (ig_opt == 1) 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(k) - Sxxprev(itheta) + ! + endif + ! + dSxx = max(dSxx, 0.0) + ! + ! 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) + ! endif ! else ! TL: option to add future parameterisations here for e.g. coral reef type coasts From b25313bc631e1bacbd9c29ec21bf5c6cb6580f1e Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 15 Apr 2026 13:31:25 +0200 Subject: [PATCH 19/25] - Suggestion of Marlies to keep the original formulation for Hmx_ig --- source/src/snapwave/snapwave_solver.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 05de09b1f..6bbcd17a2 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -98,8 +98,8 @@ subroutine compute_wave_field() if (igwaves) then do k = 1, no_nodes ! - !Hmx_ig(k) = 0.88 / kwav_ig(k) * tanh(gamma_ig * kwav_ig(k) * depth(k) / 0.88) ! Note - uses gamma_ig - Hmx_ig(k) = gamma_ig * depth(k) + 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 endif From c1183ce168b96ebfe223cb042e50dbf4a9d471f0 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 15 Apr 2026 13:30:48 +0200 Subject: [PATCH 20/25] - Bugfix in calculating kwav_ig! --- source/src/snapwave/snapwave_solver.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 6bbcd17a2..15fdc6c39 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -77,8 +77,8 @@ subroutine compute_wave_field() 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 ! From d043e14590a3eed05dc1d21fd9f1576ae05b75ee Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 2 Apr 2026 18:05:17 +0200 Subject: [PATCH 21/25] - Add back Neumann message --- source/src/snapwave/snapwave_domain.f90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90 index cd9388875..746d91102 100644 --- a/source/src/snapwave/snapwave_domain.f90 +++ b/source/src/snapwave/snapwave_domain.f90 @@ -320,6 +320,9 @@ subroutine initialize_snapwave_domain() ! call neuboundaries_light(x, y, msk, no_nodes, neumannconnected) ! + write(logstr,*)'SnapWave: Neumann connected boundaries found ...' + call write_log(logstr, 0) + ! else ! neumannconnected = 0 From 3aa0ed6720121fc1f409a9bdfb531269c6c01d6a Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 29 Apr 2026 10:18:05 +0200 Subject: [PATCH 22/25] - match in and output --- source/src/snapwave/snapwave_solver.f90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 15fdc6c39..6b0cec91d 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -487,7 +487,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, dtheta, 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) ! ! inout: alphaig_local, srcig_local - eeprev, eeprev_ig, cgprev, beta_local ! in: the rest @@ -575,7 +575,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! ! Actual determining of source term - every first sweep of iteration ! - call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, 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) ! endif ! @@ -1108,8 +1108,6 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d ! 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 ! ! Internal variables @@ -1198,10 +1196,9 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d ! 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 + ! Sxx is pre-computed for all nodes before this parallel loop ! - Sxxprev(itheta) = w(1, itheta, k) * Sxx(k1) + w(2, itheta, k) * Sxx(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) From 15a773155ffe1c0cbc7feb1fdff262b1ff8a721d Mon Sep 17 00:00:00 2001 From: Leynse Date: Mon, 9 Mar 2026 09:40:09 +0100 Subject: [PATCH 23/25] - Change default snapwave_gammaig to 0.7 --- source/src/sfincs_snapwave.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index 7f4a7cb51..554e677e2 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -642,7 +642,7 @@ subroutine read_snapwave_input() ! 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.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 From 17390a7b1fa2fc31e0ecedbf16bc7ba0206204e9 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 29 Apr 2026 10:43:36 +0200 Subject: [PATCH 24/25] - Add gamma_fac_br as in PR #281 --- source/src/sfincs_snapwave.f90 | 3 +- source/src/snapwave/snapwave_data.f90 | 1 + source/src/snapwave/snapwave_solver.f90 | 45 ++++++++++++++++++++----- 3 files changed, 40 insertions(+), 9 deletions(-) diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index 554e677e2..ecfd5b81c 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -642,7 +642,8 @@ subroutine read_snapwave_input() ! 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.7) ! 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 diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90 index 4fbb99db5..add49c180 100644 --- a/source/src/snapwave/snapwave_data.f90 +++ b/source/src/snapwave/snapwave_data.f90 @@ -259,6 +259,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 diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 6b0cec91d..8eb704f47 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 @@ -153,7 +153,7 @@ subroutine compute_wave_field() 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) ! @@ -178,7 +178,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & 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 @@ -312,6 +312,7 @@ 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 @@ -487,7 +488,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, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, 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 @@ -575,7 +576,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & ! ! Actual determining of source term - every first sweep of iteration ! - 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) + 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 ! @@ -1077,7 +1078,7 @@ subroutine baldock (rho, g, alfa, gamma, depth, H, T, iexp, Dw, Hmax) end subroutine baldock - subroutine 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) + subroutine 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) ! ! Determining of IG source term as defined in Leijnse et al. 2024 ! @@ -1103,7 +1104,9 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d 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) :: 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 ! ! Inout variables real*4, dimension(:,:), intent(inout) :: alphaig_local ! Local infragravity wave shoaling parameter alpha @@ -1126,6 +1129,12 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d 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 + 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 ! @@ -1136,7 +1145,15 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d ! 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) @@ -1249,6 +1266,18 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d ! 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 From 6fbc1fb3600af50b281029472aa0ccf719f8f9cd Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 2 Apr 2026 14:11:47 +0200 Subject: [PATCH 25/25] - Add user definable snapwave_relax_factor_DoverE and snapwave_relax_factor_DoverA keywords to sfincs_snapwave - By default set to 1.0 for both for now (as in current main) --- source/src/sfincs_snapwave.f90 | 2 ++ source/src/snapwave/snapwave_data.f90 | 2 ++ source/src/snapwave/snapwave_solver.f90 | 39 +++++++++++++------------ 3 files changed, 24 insertions(+), 19 deletions(-) diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index ecfd5b81c..4520995fd 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -637,6 +637,8 @@ 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 ! diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90 index add49c180..e0946ddff 100644 --- a/source/src/snapwave/snapwave_data.f90 +++ b/source/src/snapwave/snapwave_data.f90 @@ -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 diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 8eb704f47..e5f8bf2c8 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -139,17 +139,18 @@ 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, baldock_exponent, & - aa, sig, jadcgdx, sigmin, sigmax,& - c_dispT, DoverE, WsorE, WsorA, SwE, SwA, Tpini, & + 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, & vegetation, no_secveg, veg_ah, veg_bstems, veg_Nstems, veg_Cd, Dveg, & @@ -164,18 +165,19 @@ subroutine compute_wave_field() 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, baldock_exponent, & - aa, sig, jadcgdx, sigmin, sigmax,& - c_dispT, DoverE, WsorE, WsorA, SwE, SwA, Tpini, & - igwaves, kwav_ig, cg_ig,H_ig,ctheta_ig,Hmx_ig, ee_ig,fw_ig, & + 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, gamma_fac_br, eeinc2ig, Tinc2ig, alphaigfac, shinc2ig, iterative_srcig) @@ -237,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 ! @@ -295,9 +299,6 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, & 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 :: relax_factor_DoverA = 0.25 ! underrelaxation factor for DoverA (set to 1.0 to disable) - real*4 :: relax_factor_DoverE = 0.25 ! underrelaxation factor for DoverE (set to 1.0 to disable) real*4 :: oneoverdt real*4 :: oneover2dtheta real*4 :: rhog8