diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj index 26372260..4e5db456 100644 --- a/source/sfincs_lib/sfincs_lib.vfproj +++ b/source/sfincs_lib/sfincs_lib.vfproj @@ -78,7 +78,15 @@ - + + + + + + + + + diff --git a/source/src/Makefile.am b/source/src/Makefile.am index 718eccf2..34b2103a 100644 --- a/source/src/Makefile.am +++ b/source/src/Makefile.am @@ -48,7 +48,7 @@ libsfincs_la_SOURCES = \ sfincs_output.f90 \ sfincs_momentum.f90 \ sfincs_wavemaker.f90 \ - sfincs_infiltration.f90 \ + sfincs_infiltration.F90 \ sfincs_lib.f90 \ sfincs_bmi.f90 diff --git a/source/src/quadtree.F90 b/source/src/quadtree.F90 index 998ca03a..31d34000 100644 --- a/source/src/quadtree.F90 +++ b/source/src/quadtree.F90 @@ -852,13 +852,19 @@ subroutine find_uv_points_intersected_by_polyline(uv_indices,vertices,nint,xpol, do m = m0, m1 do n = n0, n1 ! + ! First check whether point is on the quadtree grid at all nm = find_quadtree_cell_by_index(n, m, iref) ! if (nm==0) then cycle endif ! + ! Second check is whether it is on the active SFINCS mask nm = index_sfincs_in_quadtree(nm) + ! + if (nm==0) then + cycle + endif ! ! Right (same level or coarser) ! diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 7d115e18..06714d37 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -136,6 +136,7 @@ module sfincs_data character*256 :: netampfile character*256 :: netamprfile character*256 :: netspwfile + character*256 :: netinfiltrationfile character*256 :: scsfile character*256 :: smaxfile character*256 :: sefffile @@ -532,7 +533,7 @@ module sfincs_data real*4, dimension(:), allocatable :: cg real*4, dimension(:), allocatable :: qb real*4, dimension(:), allocatable :: betamean - real*4, dimension(:), allocatable :: srcsh + real*4, dimension(:), allocatable :: srcig real*4, dimension(:), allocatable :: alphaig ! real*4, dimension(:), allocatable :: tauwavv diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 421994a5..3b39ee85 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -88,7 +88,7 @@ subroutine initialize_processes() ! ! Check for time-spatially varying meteo ! - if (spwfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. amprfile(1:4) /= 'none' .or. ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. netamprfile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then + if (prcpfile(1:4) /= 'none' .or. spwfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. amprfile(1:4) /= 'none' .or. ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. netamprfile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then meteo3d = .true. write(*,*)'Turning on flag: meteo3d' endif @@ -1804,18 +1804,10 @@ subroutine initialize_roughness() subroutine initialize_infiltration() ! use sfincs_data + use sfincs_infiltration ! implicit none ! - real*4, dimension(:), allocatable :: rghfield - ! - integer :: idummy - integer :: ip - integer :: nm - integer :: nmu - integer :: n - integer :: m - ! ! INFILTRATION ! ! Infiltration only works when rainfall is activated ! If you want infiltration without rainfall, use a precip file with 0.0s @@ -1824,296 +1816,46 @@ subroutine initialize_infiltration() ! infiltration = .false. ! - ! Four options for infiltration: - ! - ! 1) Spatially-uniform constant infiltration - ! Requires: - - ! 2) Spatially-varying constant infiltration - ! Requires: qinfmap (does not require qinffield !) - ! 3) Spatially-varying infiltration with CN numbers (old) - ! Requires: cumprcp, cuminf, qinfmap, qinffield - ! 4) Spatially-varying infiltration with CN numbers (new) - ! Requires: qinfmap, qinffield, qinffield, ksfield, scs_P1, scs_F1, scs_Se and scs_rain (but not necessarily cuminf and cumprcp) - ! - ! cumprcp and cuminf are stored in the netcdf output if store_cumulative_precipitation == .true. which is the default - ! - ! We need to keep cumprcp and cuminf in memory when: - ! a) store_cumulative_precipitation == .true. - ! or: - ! b) inftype == 'cna' or inftype == 'cnb' + ! Two input types for infiltration: + ! a) Netcdf input file for values on quadtree grid (CNA, CNB, GA, Horton). Required for quadtree input! Check for existence first, if not check: + ! b) Original input type with simple input keyword (qinf) or binary input files (qinffile, CNA, CNB, GA, Horton). Only works for original regular grid! ! - ! First we determine precipitation type - ! - if (precip) then + if (precip) then !TL: CHECK - for Horton this should not be a requirement anymore !? ! - if (qinf > 0.0) then - ! - ! Spatially-uniform constant infiltration (specified as +mm/hr) - ! - inftype = 'con' - infiltration = .true. - ! - elseif (qinffile /= 'none') then - ! - ! Spatially-varying constant infiltration - ! - inftype = 'c2d' - infiltration = .true. - ! - elseif (scsfile /= 'none') then - ! - ! Spatially-varying infiltration with CN numbers (old) - ! - inftype = 'cna' - infiltration = .true. - ! - elseif (sefffile /= 'none') then - ! - ! Spatially-varying infiltration with CN numbers (new) - ! - inftype = 'cnb' - infiltration = .true. - ! - elseif (psifile /= 'none') then - ! - ! The Green-Ampt (GA) model for infiltration - ! - inftype = 'gai' - infiltration = .true. - ! - elseif (f0file /= 'none') then - ! - ! The Horton Equation model for infiltration - ! - inftype = 'hor' - infiltration = .true. - store_meteo = .true. - ! - endif + if (netinfiltrationfile /= 'none') then + ! + if (use_quadtree) then + call read_infiltration_file_netcdf() + else + write(*,*)'Warning : new Netcdf infiltration input format can only be specified for quadtree mesh model!' + endif + ! + else! Original + ! + if (use_quadtree .eqv. .false.) then + call read_infiltration_file_original() + else + write(*,*)'Warning : infiltration input for quadtree mesh can only be specified using the new Netcdf infiltration input format!' + endif + ! + endif ! - if (precip) then - ! - ! We need cumprcp and cuminf - ! - allocate(cumprcp(np)) - cumprcp = 0.0 - ! - allocate(cuminf(np)) - cuminf = 0.0 - ! - endif + ! We need cumprcp and cuminf ! - ! Now allocate and read spatially-varying inputs + allocate(cumprcp(np)) + cumprcp = 0.0 ! - if (infiltration) then - ! - allocate(qinfmap(np)) - qinfmap = 0.0 - ! - endif + allocate(cuminf(np)) + cuminf = 0.0 ! - if (inftype == 'con') then - ! - ! Spatially-uniform constant infiltration (specified as +mm/hr) - ! - write(*,*)'Turning on process: Spatially-uniform constant infiltration' - ! - allocate(qinffield(np)) - do nm = 1, np - if (subgrid) then - if (subgrid_z_zmin(nm) > qinf_zmin) then - qinffield(nm) = qinf - else - qinffield(nm) = 0.0 - endif - else - if (zb(nm) > qinf_zmin) then - qinffield(nm) = qinf - else - qinffield(nm) = 0.0 - endif - endif - enddo - ! - elseif (inftype == 'c2d') then - ! - ! Spatially-varying constant infiltration - ! - write(*,*)'Turning on process: Spatially-varying constant infiltration' - ! - ! Read spatially-varying infiltration (only binary, specified in +mm/hr) - ! - write(*,*)'Reading ', trim(qinffile), ' ...' - allocate(qinffield(np)) - open(unit = 500, file = trim(qinffile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) - ! - do nm = 1, np - qinffield(nm) = qinffield(nm)/3.6e3/1.0e3 ! convert to +m/s - enddo - ! - elseif (inftype == 'cna') then - ! - ! Spatially-varying infiltration with CN numbers (old) - ! - write(*,*)'Turning on process: Infiltration (via Curve Number method - A)' - ! - allocate(qinffield(np)) - qinffield = 0.0 - ! - write(*,*)'Reading ',trim(scsfile) - open(unit = 500, file = trim(scsfile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) - ! - ! already convert qinffield from inches to m here - do nm = 1, np - qinffield(nm) = qinffield(nm)*0.0254 !to m - enddo - ! - elseif (inftype == 'cnb') then - ! - ! Spatially-varying infiltration with CN numbers (new) - ! - write(*,*)'Turning on process: Infiltration (via Curve Number method - B)' - ! - ! Allocate Smax - allocate(qinffield(np)) - qinffield = 0.0 - write(*,*)'Reading ',trim(smaxfile) - open(unit = 500, file = trim(smaxfile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) - ! - ! Allocate Se - allocate(scs_Se(np)) - scs_Se = 0.0 - write(*,*)'Reading ',trim(sefffile) - open(unit = 501, file = trim(sefffile), form = 'unformatted', access = 'stream') - read(501)scs_Se - close(501) - ! - ! Compute recovery ! Equation 4-36 - ! Allocate Ks - allocate(ksfield(np)) - ksfield = 0.0 - write(*,*)'Reading ',trim(ksfile) - open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') - read(502)ksfield - close(502) - ! - ! Compute recovery ! Equation 4-36 - allocate(inf_kr(np)) - inf_kr = sqrt(ksfield/25.4) / 75 ! Note that we assume ksfield to be in mm/hr, convert it here to inch/hr (/25.4) - ! /75 is conversion to recovery rate (in days) - ! - ! Allocate support variables - allocate(scs_P1(np)) - scs_P1 = 0.0 - allocate(scs_F1(np)) - scs_F1 = 0.0 - allocate(rain_T1(np)) - rain_T1 = 0.0 - allocate(scs_S1(np)) - scs_S1 = 0.0 - allocate(scs_rain(np)) - scs_rain = 0 - ! - elseif (inftype == 'gai') then - ! - ! Spatially-varying infiltration with the Green-Ampt (GA) model - ! - write(*,*)'Turning on process: Infiltration (via Green-Ampt)' - ! - ! Allocate suction head at the wetting front - allocate(GA_head(np)) - GA_head = 0.0 - write(*,*)'Reading ',trim(psifile) - open(unit = 500, file = trim(psifile), form = 'unformatted', access = 'stream') - read(500)GA_head - close(500) - ! - ! Allocate maximum soil moisture deficit - allocate(GA_sigma_max(np)) - GA_sigma_max = 0.0 - write(*,*)'Reading ',trim(sigmafile) - open(unit = 501, file = trim(sigmafile), form = 'unformatted', access = 'stream') - read(501)GA_sigma_max - close(501) - ! - ! Allocate saturated hydraulic conductivity - allocate(ksfield(np)) - ksfield = 0.0 - write(*,*)'Reading ',trim(ksfile) - open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') - read(502)ksfield - close(502) - ! - ! Compute recovery ! Equation 4-36 - allocate(inf_kr(np)) - inf_kr = sqrt(ksfield/25.4) / 75 ! Note that we assume ksfield to be in mm/hr, convert it here to inch/hr (/25.4) - ! /75 is conversion to recovery rate (in days) - - allocate(rain_T1(np)) ! minimum amount of time that a soil must remain in recovery - rain_T1 = 0.0 - ! Allocate support variables - allocate(GA_sigma(np)) ! variable for sigma_max_du - GA_sigma = GA_sigma_max - allocate(GA_F(np)) ! total infiltration - GA_F = 0.0 - allocate(GA_Lu(np)) ! depth of upper soil recovery zone - GA_Lu = 4 *sqrt(25.4) * sqrt(ksfield) ! Equation 4-33 - ! - ! Input values for green-ampt are in mm and mm/hr, but computation is in m a m/s - GA_head = GA_head/1000 ! from mm to m - GA_Lu = GA_Lu/1000 ! from mm to m - ksfield = ksfield/1000/3600 ! from mm/hr to m/s - ! - ! First time step doesnt have an estimate yet - allocate(qinffield(np)) - qinffield(nm) = 0.00 - ! - elseif (inftype == 'hor') then - ! - ! Spatially-varying infiltration with the modified Horton Equation - ! - write(*,*)'Turning on process: Infiltration (via modified Horton)' - ! - ! Horton: final infiltration capacity (fc) - ! Note that qinffield = horton_fc - allocate(horton_fc(np)) - horton_fc = 0.0 - write(*,*)'Reading ',trim(fcfile) - open(unit = 500, file = trim(fcfile), form = 'unformatted', access = 'stream') - read(500)horton_fc - close(500) - ! - ! Horton: initial infiltration capacity (f0) - allocate(horton_f0(np)) - horton_f0 = 0.0 - write(*,*)'Reading ',trim(f0file) - open(unit = 501, file = trim(f0file), form = 'unformatted', access = 'stream') - read(501)horton_f0 - close(501) - ! - ! Prescribe the current estimate (for output only; initial capacity) - qinffield = horton_f0/3600/1000 - ! - ! Empirical constant (1/hr) k => note that this is different than ks used in Curve Number and Green-Ampt - allocate(horton_kd(np)) - horton_kd = 0.0 - write(*,*)'Reading ',trim(kdfile) - open(unit = 502, file = trim(kdfile), form = 'unformatted', access = 'stream') - read(502)horton_kd - close(502) - write(*,*)'Using constant recovery rate that is based on constant factor relative to ',trim(kdfile) - ! - ! Estimate of time - allocate(rain_T1(np)) - rain_T1 = 0.0 - ! - endif + ! Now allocatespatially-varying inputs + ! + if (infiltration) then + ! + allocate(qinfmap(np)) + qinfmap = 0.0 + ! + endif ! else ! @@ -2291,8 +2033,8 @@ subroutine initialize_hydro() qb = 0.0 allocate(betamean(np)) betamean = 0.0 - allocate(srcsh(np)) - srcsh = 0.0 + allocate(srcig(np)) + srcig = 0.0 allocate(alphaig(np)) alphaig = 0.0 endif diff --git a/source/src/sfincs_infiltration.F90 b/source/src/sfincs_infiltration.F90 new file mode 100644 index 00000000..4defd658 --- /dev/null +++ b/source/src/sfincs_infiltration.F90 @@ -0,0 +1,802 @@ +#define NF90(nf90call) call handle_err(nf90call,__FILE__,__LINE__) +module sfincs_infiltration + ! + use netcdf + ! + type net_type_inf + integer :: ncid + integer :: np_dimid + integer :: mask_varid, type_varid + integer :: qinffield_varid, scs_varid, scs_se_varid, ksfield_varid + end type + type(net_type_inf) :: net_file_inf + ! +contains + + subroutine read_infiltration_file_original() + ! + ! Six options for infiltration: + ! + ! 1) Spatially-uniform constant infiltration + ! Requires: - + ! 2) Spatially-varying constant infiltration + ! Requires: qinfmap (does not require qinffield !) + ! 3) Spatially-varying infiltration with CN numbers (old) + ! Requires: cumprcp, cuminf, qinfmap, qinffield + ! 4) Spatially-varying infiltration with CN numbers (new) + ! Requires: qinfmap, qinffield, ksfield, scs_P1, scs_F1, scs_Se and scs_rain (but not necessarily cuminf and cumprcp) + ! 5) Spatially-varying infiltration with the Green-Ampt (GA) model + ! Requires: qinfmap, qinffield, ksfield, GA_head, GA_sigma_max, GA_Lu + ! 6) Spatially-varying infiltration with the modified Horton Equation + ! Requires: qinfmap, qinffield, horton_fc, horton_f0 + ! + ! cumprcp and cuminf are stored in the netcdf output if store_cumulative_precipitation == .true. which is the default + ! + ! We need to keep cumprcp and cuminf in memory when: + ! a) store_cumulative_precipitation == .true. + ! or: + ! b) inftype == 'cna' or inftype == 'cnb' + ! + ! First we determine infiltration type + ! + use sfincs_data + ! + implicit none + ! + integer :: nm + ! + if (qinf > 0.0) then + ! + ! Spatially-uniform constant infiltration (specified as +mm/hr) + ! + inftype = 'con' + infiltration = .true. + ! + elseif (qinffile /= 'none') then + ! + ! Spatially-varying constant infiltration + ! + inftype = 'c2d' + infiltration = .true. + ! + elseif (scsfile /= 'none') then + ! + ! Spatially-varying infiltration with CN numbers (old) + ! + inftype = 'cna' + infiltration = .true. + ! + elseif (sefffile /= 'none') then + ! + ! Spatially-varying infiltration with CN numbers (new) + ! + inftype = 'cnb' + infiltration = .true. + ! + elseif (psifile /= 'none') then + ! + ! The Green-Ampt (GA) model for infiltration + ! + inftype = 'gai' + infiltration = .true. + ! + elseif (f0file /= 'none') then + ! + ! The Horton Equation model for infiltration + ! + inftype = 'hor' + infiltration = .true. + store_meteo = .true. + ! + endif + ! + if (inftype == 'con') then + ! + ! Spatially-uniform constant infiltration (specified as +mm/hr) + ! + write(*,*)'Turning on process: Spatially-uniform constant infiltration' + ! + allocate(qinffield(np)) + do nm = 1, np + if (subgrid) then + if (subgrid_z_zmin(nm) > qinf_zmin) then + qinffield(nm) = qinf + else + qinffield(nm) = 0.0 + endif + else + if (zb(nm) > qinf_zmin) then + qinffield(nm) = qinf + else + qinffield(nm) = 0.0 + endif + endif + enddo + ! + elseif (inftype == 'c2d') then + ! + ! Spatially-varying constant infiltration + ! + write(*,*)'Turning on process: Spatially-varying constant infiltration' + ! + ! Read spatially-varying infiltration (only binary, specified in +mm/hr) + ! + write(*,*)'Reading ', trim(qinffile), ' ...' + allocate(qinffield(np)) + open(unit = 500, file = trim(qinffile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + do nm = 1, np + qinffield(nm) = qinffield(nm)/3.6e3/1.0e3 ! convert to +m/s + enddo + ! + elseif (inftype == 'cna') then + ! + ! Spatially-varying infiltration with CN numbers (old) + ! + write(*,*)'Turning on process: Infiltration (via Curve Number method - A)' + ! + allocate(qinffield(np)) + qinffield = 0.0 + ! + write(*,*)'Reading ',trim(scsfile) + open(unit = 500, file = trim(scsfile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + ! already convert qinffield from inches to m here + do nm = 1, np + qinffield(nm) = qinffield(nm)*0.0254 !to m + enddo + ! + elseif (inftype == 'cnb') then + ! + ! Spatially-varying infiltration with CN numbers (new) + ! + write(*,*)'Turning on process: Infiltration (via Curve Number method - B)' + ! + ! Allocate Smax + allocate(qinffield(np)) + qinffield = 0.0 + write(*,*)'Reading ',trim(smaxfile) + open(unit = 500, file = trim(smaxfile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + ! Allocate Se + allocate(scs_Se(np)) + scs_Se = 0.0 + write(*,*)'Reading ',trim(sefffile) + open(unit = 501, file = trim(sefffile), form = 'unformatted', access = 'stream') + read(501)scs_Se + close(501) + ! + ! Compute recovery ! Equation 4-36 + ! Allocate Ks + allocate(ksfield(np)) + ksfield = 0.0 + write(*,*)'Reading ',trim(ksfile) + open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') + read(502)ksfield + close(502) + ! + ! Compute recovery ! Equation 4-36 + allocate(inf_kr(np)) + inf_kr = sqrt(ksfield/25.4) / 75 ! Note that we assume ksfield to be in mm/hr, convert it here to inch/hr (/25.4) + ! /75 is conversion to recovery rate (in days) + ! + ! Allocate support variables + allocate(scs_P1(np)) + scs_P1 = 0.0 + allocate(scs_F1(np)) + scs_F1 = 0.0 + allocate(rain_T1(np)) + rain_T1 = 0.0 + allocate(scs_S1(np)) + scs_S1 = 0.0 + allocate(scs_rain(np)) + scs_rain = 0 + ! + elseif (inftype == 'gai') then + ! + ! Spatially-varying infiltration with the Green-Ampt (GA) model + ! + write(*,*)'Turning on process: Infiltration (via Green-Ampt)' + ! + ! Allocate suction head at the wetting front + allocate(GA_head(np)) + GA_head = 0.0 + write(*,*)'Reading ',trim(psifile) + open(unit = 500, file = trim(psifile), form = 'unformatted', access = 'stream') + read(500)GA_head + close(500) + ! + ! Allocate maximum soil moisture deficit + allocate(GA_sigma_max(np)) + GA_sigma_max = 0.0 + write(*,*)'Reading ',trim(sigmafile) + open(unit = 501, file = trim(sigmafile), form = 'unformatted', access = 'stream') + read(501)GA_sigma_max + close(501) + ! + ! Allocate saturated hydraulic conductivity + allocate(ksfield(np)) + ksfield = 0.0 + write(*,*)'Reading ',trim(ksfile) + open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') + read(502)ksfield + close(502) + ! + ! Compute recovery ! Equation 4-36 + allocate(inf_kr(np)) + inf_kr = sqrt(ksfield/25.4) / 75 ! Note that we assume ksfield to be in mm/hr, convert it here to inch/hr (/25.4) + ! /75 is conversion to recovery rate (in days) + + allocate(rain_T1(np)) ! minimum amount of time that a soil must remain in recovery + rain_T1 = 0.0 + ! Allocate support variables + allocate(GA_sigma(np)) ! variable for sigma_max_du + GA_sigma = GA_sigma_max + allocate(GA_F(np)) ! total infiltration + GA_F = 0.0 + allocate(GA_Lu(np)) ! depth of upper soil recovery zone + GA_Lu = 4 *sqrt(25.4) * sqrt(ksfield) ! Equation 4-33 + ! + ! Input values for green-ampt are in mm and mm/hr, but computation is in m a m/s + GA_head = GA_head/1000 ! from mm to m + GA_Lu = GA_Lu/1000 ! from mm to m + ksfield = ksfield/1000/3600 ! from mm/hr to m/s + ! + ! First time step doesnt have an estimate yet + allocate(qinffield(np)) + qinffield(nm) = 0.00 + ! + elseif (inftype == 'hor') then + ! + ! Spatially-varying infiltration with the modified Horton Equation + ! + write(*,*)'Turning on process: Infiltration (via modified Horton)' + ! + ! Horton: final infiltration capacity (fc) + ! Note that qinffield = horton_fc + allocate(horton_fc(np)) + horton_fc = 0.0 + write(*,*)'Reading ',trim(fcfile) + open(unit = 500, file = trim(fcfile), form = 'unformatted', access = 'stream') + read(500)horton_fc + close(500) + ! + ! Horton: initial infiltration capacity (f0) + allocate(horton_f0(np)) + horton_f0 = 0.0 + write(*,*)'Reading ',trim(f0file) + open(unit = 501, file = trim(f0file), form = 'unformatted', access = 'stream') + read(501)horton_f0 + close(501) + ! + ! Prescribe the current estimate (for output only; initial capacity) + qinffield = horton_f0/3600/1000 + ! + ! Empirical constant (1/hr) k => note that this is different than ks used in Curve Number and Green-Ampt + allocate(horton_kd(np)) + horton_kd = 0.0 + write(*,*)'Reading ',trim(kdfile) + open(unit = 502, file = trim(kdfile), form = 'unformatted', access = 'stream') + read(502)horton_kd + close(502) + write(*,*)'Using constant recovery rate that is based on constant factor relative to ',trim(kdfile) + ! + ! Estimate of time + allocate(rain_T1(np)) + rain_T1 = 0.0 + ! + endif + ! + end subroutine read_infiltration_file_original + + + subroutine read_infiltration_file_netcdf() + ! + use sfincs_data + ! + implicit none + ! + !character*256, intent(in) :: netinfiltrationfile + ! + integer :: np_nc,ip, nm + ! + integer*4 :: quadtree_nr_points + ! + integer*1, dimension(:), allocatable :: quadtree_mask + integer*4, dimension(:), allocatable :: z_index + real*4, dimension(:), allocatable :: rtmpz + ! + write(*,*)'Reading Infiltration netCDF file on quadtree mesh...' + ! + NF90(nf90_open(trim(netinfiltrationfile), NF90_CLOBBER, net_file_inf%ncid)) + ! + ! Get dimensions id's: nr points + ! + NF90(nf90_inq_dimid(net_file_inf%ncid, "mesh2d_nFaces", net_file_inf%np_dimid)) + ! + ! Get dimensions sizes + ! + NF90(nf90_inquire_dimension(net_file_inf%ncid, net_file_inf%np_dimid, len = np_nc)) ! nr of cells + ! + ! Get variable id's + NF90(nf90_inq_varid(net_file_inf%ncid, 'mask', net_file_inf%mask_varid)) + ! + !NF90(nf90_inq_varid(net_file_inf%ncid, 'type', net_file_inf%type_varid)) + inftype = 'cnb' !TL: later determine this based on the netcdf file + ! + ! Allocate variables + allocate(z_index(np_nc)) + allocate(rtmpz(np_nc)) + !allocate(kcs(np)) + allocate(quadtree_mask(np)) + ! + ! Netcdf quadtree infiltrationfile contains data for the entire quadtree. So also for points with kcs==0 ! + ! This means that the data needs to be re-mapped to the active cell indices: + if (use_quadtree) then + ! + do ip = 1, np_nc !np_nc~quadtree_nr_points + ! + nm = index_sfincs_in_quadtree(ip) + ! + if (nm>0) then + z_index(nm) = ip + endif + ! + enddo + ! + !else + ! ! Regular grid > Tl: only keep regular if we want to support this in the end + ! do nm = 1, np_nc + ! z_index(nm) = nm + ! enddo + endif + ! + ! Read Z points + ! + NF90(nf90_get_var(net_file_inf%ncid, net_file_inf%mask_varid, rtmpz(:))) + ! + do ip = 1, np + quadtree_mask(ip) = rtmpz(z_index(ip)) + enddo + ! + ! Check whether number of cells matches + !if (np_nc /= np) then + ! write(*,'(a,i8,a,i8,a)')'Error! Number of cells in netinfiltrationfile ',np_nc,' does not match number of active cells in mesh ',np,' ! Abort reading infiltration input...',np_nc + ! continue + !endif + ! + ! Check whether incoming mask is matches the one as retrieved from netcdf quadtree grid file (our general 'msk') + if (all(quadtree_mask == kcs)) then + write(*,*)'Check succesfull: The mask arrays are the same.' + else + write(*,*)'Error! Mask of netinfiltrationfile does not match with actiave mask in mesh! Abort reading infiltration input...' + continue + endif + ! + ! Now read in infiltration type specific variables + ! + if (inftype == 'cna') then + ! + ! Spatially-varying infiltration with CN numbers (old) + ! + write(*,*)'Turning on process: Infiltration (via Curve Number method - A)' + ! + ! Allocate qinffield + allocate(qinffield(np)) + ! + ! Read variables for cna method + ! + ! Get variable id's + NF90(nf90_inq_varid(net_file_inf%ncid, 'scs', net_file_inf%scs_varid)) + ! Read Z points + ! + NF90(nf90_get_var(net_file_inf%ncid, net_file_inf%scs_varid, rtmpz(:))) + do ip = 1, np + qinffield(ip) = rtmpz(z_index(ip)) + enddo + ! + ! already convert qinffield from inches to m here + do nm = 1, np + qinffield(nm) = qinffield(nm)*0.0254 !to m + enddo + ! + elseif (inftype == 'cnb') then + ! + ! Spatially-varying infiltration with CN numbers (with recovery) + ! + write(*,*)'Turning on process: Infiltration (via Curve Number method - B)' + ! + ! Allocate Smax, Se, Ks + allocate(qinffield(np)) + allocate(scs_Se(np)) + allocate(ksfield(np)) + allocate(inf_kr(np)) + ! + ! Read variables for cnb method + ! + ! Get variable id's + NF90(nf90_inq_varid(net_file_inf%ncid, 'smax', net_file_inf%qinffield_varid)) + NF90(nf90_inq_varid(net_file_inf%ncid, 'seff', net_file_inf%scs_se_varid)) + NF90(nf90_inq_varid(net_file_inf%ncid, 'ks', net_file_inf%ksfield_varid)) + ! From Kees below: note that scs_Se is S and qinffield is Smax + + ! Read Z points + ! + NF90(nf90_get_var(net_file_inf%ncid, net_file_inf%qinffield_varid, rtmpz(:))) + do ip = 1, np + qinffield(ip) = rtmpz(z_index(ip)) + enddo + ! + NF90(nf90_get_var(net_file_inf%ncid, net_file_inf%scs_se_varid, rtmpz(:))) + do ip = 1, np + scs_Se(ip) = rtmpz(z_index(ip)) + enddo + ! + NF90(nf90_get_var(net_file_inf%ncid, net_file_inf%ksfield_varid, rtmpz(:))) + do ip = 1, np + ksfield(ip) = rtmpz(z_index(ip)) + enddo + ! + ! Compute recovery ! Equation 4-36 + inf_kr = sqrt(ksfield/25.4) / 75 ! Note that we assume ksfield to be in mm/hr, convert it here to inch/hr (/25.4) + ! ! /75 is conversion to recovery rate (in days) + ! + ! Allocate support variables + allocate(scs_P1(np)) + scs_P1 = 0.0 + allocate(scs_F1(np)) + scs_F1 = 0.0 + allocate(rain_T1(np)) + rain_T1 = 0.0 + allocate(scs_S1(np)) + scs_S1 = 0.0 + allocate(scs_rain(np)) + scs_rain = 0 + ! + endif + ! + infiltration = .true. + ! + NF90(nf90_close(net_file_inf%ncid)) + ! + ! Deallocate vars + deallocate(rtmpz) + deallocate(z_index) + deallocate(quadtree_mask) + ! + end subroutine read_infiltration_file_netcdf + + + subroutine update_infiltration_map(dt) + ! + ! Update infiltration rates in each grid cell + ! + use sfincs_data + ! + implicit none + ! + integer nm + ! + real*4 :: Qq + real*4 :: I + real*4 :: hh_local, a + real*4 :: dt + ! + if (inftype == 'con' .or. inftype == 'c2d') then + ! + ! Infiltration rate map stays constant + ! + !$omp parallel & + !$omp private ( nm ) + !$omp do + do nm = 1, np + ! + qinfmap(nm) = qinffield(nm) + ! + ! No infiltration if there is no water + ! + if (subgrid) then + if (z_volume(nm)<=0.0) then + qinfmap(nm) = 0.0 + endif + else + if (zs(nm)<=zb(nm)) then + qinfmap(nm) = 0.0 + endif + endif + ! + cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt + netprcp(nm) = netprcp(nm) - qinfmap(nm) + ! + enddo + !$omp end do + !$omp end parallel + ! + elseif (inftype == 'cna') then + ! + ! Determine infiltration rate with Curve Number (old method; no recovery) + ! +! !$acc update host(cumprcp), async(1) + ! + !$omp parallel & + !$omp private ( Qq,I,nm ) + !$omp do + do nm = 1, np + ! + ! Check if Ia (0.2 x S) is larger than cumulative rainfall + ! + if (cumprcp(nm) > sfacinf*qinffield(nm)) then ! qinffield is S + ! + ! Compute runoff as function of rain + ! + Qq = (cumprcp(nm) - sfacinf*qinffield(nm))**2 / (cumprcp(nm) + (1.0 - sfacinf)*qinffield(nm)) ! cumulative runoff in m + I = cumprcp(nm) - Qq ! cumulative infiltration in m + qinfmap(nm) = (I - cuminf(nm))/dt ! infiltration in m/s + ! + else + ! + ! Everything still infiltrating + ! + qinfmap(nm) = prcp(nm) + ! + endif + ! + cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt + netprcp(nm) = netprcp(nm) - qinfmap(nm) + ! + enddo + !$omp end do + !$omp end parallel + ! + !$acc update device(qinfmap), async(1) + ! + ! For now, curve number infiltration is done on the CPU +! !$acc update device(cuminf), async(1) + ! + ! Provide update to user + !write(*,'(a,f5.3,a,f10.3,a)') ' update from SCS-CN method: Average cumulative rainfall of ',sum(cumprcp)/size(cumprcp),' meter and infiltration rate of ',sum(qinfmap*3.6e3*1.0e3)/size(qinfmap) ,' mm/hour' + ! + elseif (inftype == 'cnb') then + ! + ! Determine infiltration rate with Curve Number with recovery + ! + !$omp parallel & + !$omp private ( Qq,I,nm ) + !$omp do + do nm = 1, np + ! + ! If there is precip in this grid cell for this time step + ! + if (prcp(nm) > 0.0) then + ! + ! Is raining now + ! + if (scs_rain(nm) == 1) then + ! if it was raining before; do nothing + else + ! initalise these variables + scs_P1(nm) = 0.0 ! cumulative rainfall for this 'event' + scs_F1(nm) = 0.0 ! cumulative infiltration for this 'event' + scs_S1(nm) = scs_Se(nm) ! S for this 'event' + scs_rain(nm) = 1 ! logic used to determine if there is an event ongoing + endif + ! + ! Compute cum rainfall + ! + scs_P1(nm) = scs_P1(nm) + prcp(nm)*dt + ! + ! Compute runoff + ! + if (scs_P1(nm) > (sfacinf* scs_S1(nm)) ) then ! scs_S1 is S + Qq = (scs_P1(nm) - (sfacinf*scs_S1(nm)))**2 / (scs_P1(nm) + (1.0 - sfacinf)*scs_S1(nm)) ! cumulative runoff in m + I = scs_P1(nm) - Qq ! cum infiltration this event + qinfmap(nm) = (I - scs_F1(nm))/dt ! infiltration in m/s + scs_F1(nm) = I ! cum infiltration this event + else + Qq = 0.0 ! no runoff + scs_F1(nm) = scs_P1(nm) ! all rainfall is infiltrated + qinfmap(nm) = prcp(nm) ! infiltration rate = rainfall rate + endif + ! + ! Compute "remaining S", but note that scs_Se is not used in computation + ! + scs_Se(nm) = scs_Se(nm) - qinfmap(nm)*dt + scs_Se(nm) = max(scs_Se(nm), 0.0) + qinfmap(nm) = max(qinfmap(nm), 0.0) + ! + else + ! + ! It is not raining here + if (scs_rain(nm) == 1) then + ! + ! if it was raining before; cange logic and set rate to 0 + scs_rain(nm) = 0 + qinfmap(nm) = 0.0 + rain_T1(nm) = 0.0 + ! + endif + ! + ! Add to recovery time + rain_T1(nm) = rain_T1(nm) + dt / 3600 + ! + ! compute recovery of S if time is larger than this + if (rain_T1(nm) > (0.06 / inf_kr(nm)) ) then ! Equation 4-37 from SWMM + ! note that scs_Se is S and qinffield is Smax + scs_Se(nm) = scs_Se(nm) + (inf_kr(nm) * qinffield(nm) * dt / 3600) ! scs_kr is recovery in hours + scs_Se(nm) = min(scs_Se(nm), qinffield(nm)) + ! + endif + ! + endif + ! + ! Compute cumulative values + cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt + netprcp(nm) = netprcp(nm) - qinfmap(nm) + ! + enddo + !$omp end do + !$omp end parallel + ! + !$acc update device(qinfmap), async(1) + ! + elseif (inftype == 'gai') then + ! + ! Determine infiltration rate with with the Green-Ampt (GA) model + !$omp parallel & + !$omp private ( nm ) + !$omp do + do nm = 1, np + ! + ! If there is precip in this grid cell for this time step? + if (prcp(nm) > 0.0) then + ! + ! Is raining now + if ( (prcp(nm)) < ksfield(nm)) then + ! + ! Small amounts of rainfall - infiltration is same as soil + qinfmap(nm) = prcp(nm) ! infiltration is same as rainfall + ! + else + ! + ! Larger amounts of rainfall - Equation 4-27 from SWMM manual + qinfmap(nm) = (ksfield(nm) * (1 + (GA_head(np)*GA_sigma(np)) / GA_F(nm))) + qinfmap(nm) = min(qinfmap(nm), prcp(nm)) ! never more than rainfall + qinfmap(nm) = max(qinfmap(nm), 0.0) ! and never negative + ! + endif + ! + ! Update sigma + GA_sigma(nm) = GA_sigma(nm) - (qinfmap(nm)*dt/GA_Lu(nm)) + GA_sigma(nm) = max(GA_sigma(nm),0.00) + ! + ! Update others + GA_F(nm) = GA_F(nm) + qinfmap(nm)*dt ! internal cumulative rainfall from Green-Ampt + rain_T1(nm) = 0.0 ! recovery time not started + ! + else + ! + ! Not raining here + ! + ! Add to recovery time + rain_T1(nm) = rain_T1(nm) + dt / 3600 + ! + ! compute recovery of S if time is larger than this + if (rain_T1(nm) > (0.06 / inf_kr(nm)) ) then ! Equation 4-37 from SWMM + ! + ! Update sigma + GA_sigma(nm) = GA_sigma(nm) + (inf_kr(nm) * GA_sigma_max(nm) * dt/3600) ! Equation 4-35 + GA_sigma(nm) = min(GA_sigma(nm), GA_sigma_max(nm)) ! never more than max + ! + ! Update internal cumulative rainfall + GA_F(nm) = GA_F(nm) - (inf_kr(nm) * GA_sigma_max(nm) * dt/3600*GA_Lu(nm)) ! Page 112 SWMM + GA_F(nm) = max(GA_F(nm), 0.0) ! never negative + endif + endif + ! + ! Compute cumulative values + qinffield(nm) = qinfmap(nm) + cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt + netprcp(nm) = netprcp(nm) - qinfmap(nm) + ! + enddo + !$omp end do + !$omp end parallel + ! + !$acc update device(qinfmap), async(1) + ! + elseif (inftype == 'hor') then + ! + ! Determine infiltration rate with with the Horton model + !$omp parallel & + !$omp private ( Qq,I,nm,a,hh_local) + !$omp do + do nm = 1, np + ! + ! Get local water depth estimate + if (subgrid) then + if (crsgeo) then + a = cell_area_m2(nm) + else + a = cell_area(z_flags_iref(nm)) + endif + hh_local = z_volume(nm)/a + else + hh_local = zs(nm) - zb(nm) + endif + ! + ! Check if there is water + if (hh_local> 0.0 .OR. prcp(nm) > 0.0) then + ! + ! Infiltrating here + ! + ! Count how long this is already going + if (rain_T1(nm) > 0) then + rain_T1 = 0.0 + endif + rain_T1(nm) = rain_T1(nm) - dt ! negative amount of how long it is infiltrating + ! + ! Compute estimate of infiltration ! Note that qinffield = horton_fc + I = exp(horton_kd(nm)*rain_T1(nm)/3600) ! note that horton_kd is factor in hours while dt is seconds + ! + ! Stop keeping track of this when less than 1% left (same for time) + if (I < 0.01) then + I = 0.0 ! which reduces qinfmap to horton_fc + rain_T1(nm) = rain_T1(nm) + dt ! also make sure time doesnt further decrease + qinfmap(nm) = (horton_fc(nm))/3600/1000 ! from mm/hr to m/s + else + qinfmap(nm) = (horton_fc(nm) + (horton_f0(nm) - horton_fc(nm))*I)/3600/1000 ! from mm/hr to m/s + endif + ! + ! Check how much there can infiltrate + if (hh_local > 0.00) then + Qq = prcp(nm)*dt + (zs(nm) - zb(nm)) ! Qq is estimate in meter of how much water there is + else + Qq = prcp(nm)*dt ! if no water; only compare with rainfall + endif + ! + ! Compare how much Horton want to infiltrate + I = qinfmap(nm) * dt ! I is estimate in meter of how much Horton allows + if (I > Qq) then + qinfmap(nm) = qinfmap(nm) * Qq/I ! scale Horton if capacity > available + endif + ! + else + ! + ! Not raining here NOR ponding + rain_T1(nm) = rain_T1(nm) + dt/horton_kr_kd ! positive amount of how long it is infiltrating + qinfmap(nm) = 0.0 + ! + endif + ! + ! Compute cumulative values + cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt + netprcp(nm) = netprcp(nm) - qinfmap(nm) + ! + enddo + !$omp end do + !$omp end parallel + ! + endif + ! + end subroutine + + + 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 + write(0,'("NETCDF ERROR: ",a,i6,":",a)') file,line,trim(nf90_strerror(status)) + end if + end subroutine handle_err + ! +end module diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 deleted file mode 100644 index 76897f2d..00000000 --- a/source/src/sfincs_infiltration.f90 +++ /dev/null @@ -1,318 +0,0 @@ -module sfincs_infiltration - -contains - - subroutine update_infiltration_map(dt) - ! - ! Update infiltration rates in each grid cell - ! - use sfincs_data - ! - implicit none - ! - integer nm - ! - real*4 :: Qq - real*4 :: I - real*4 :: hh_local, a - real*4 :: dt - ! - if (inftype == 'con' .or. inftype == 'c2d') then - ! - ! Infiltration rate map stays constant - ! - !$omp parallel & - !$omp private ( nm ) - !$omp do - do nm = 1, np - ! - qinfmap(nm) = qinffield(nm) - ! - ! No infiltration if there is no water - ! - if (subgrid) then - if (z_volume(nm)<=0.0) then - qinfmap(nm) = 0.0 - endif - else - if (zs(nm)<=zb(nm)) then - qinfmap(nm) = 0.0 - endif - endif - ! - cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! - enddo - !$omp end do - !$omp end parallel - ! - elseif (inftype == 'cna') then - ! - ! Determine infiltration rate with Curve Number (old method; no recovery) - ! -! !$acc update host(cumprcp), async(1) - ! - !$omp parallel & - !$omp private ( Qq,I,nm ) - !$omp do - do nm = 1, np - ! - ! Check if Ia (0.2 x S) is larger than cumulative rainfall - ! - if (cumprcp(nm) > sfacinf*qinffield(nm)) then ! qinffield is S - ! - ! Compute runoff as function of rain - ! - Qq = (cumprcp(nm) - sfacinf*qinffield(nm))**2 / (cumprcp(nm) + (1.0 - sfacinf)*qinffield(nm)) ! cumulative runoff in m - I = cumprcp(nm) - Qq ! cumulative infiltration in m - qinfmap(nm) = (I - cuminf(nm))/dt ! infiltration in m/s - ! - else - ! - ! Everything still infiltrating - ! - qinfmap(nm) = prcp(nm) - ! - endif - ! - cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! - enddo - !$omp end do - !$omp end parallel - ! - !$acc update device(qinfmap), async(1) - ! - ! For now, curve number infiltration is done on the CPU -! !$acc update device(cuminf), async(1) - ! - ! Provide update to user - !write(*,'(a,f5.3,a,f10.3,a)') ' update from SCS-CN method: Average cumulative rainfall of ',sum(cumprcp)/size(cumprcp),' meter and infiltration rate of ',sum(qinfmap*3.6e3*1.0e3)/size(qinfmap) ,' mm/hour' - ! - elseif (inftype == 'cnb') then - ! - ! Determine infiltration rate with Curve Number with recovery - ! - !$omp parallel & - !$omp private ( Qq,I,nm ) - !$omp do - do nm = 1, np - ! - ! If there is precip in this grid cell for this time step - ! - if (prcp(nm) > 0.0) then - ! - ! Is raining now - ! - if (scs_rain(nm) == 1) then - ! if it was raining before; do nothing - else - ! initalise these variables - scs_P1(nm) = 0.0 ! cumulative rainfall for this 'event' - scs_F1(nm) = 0.0 ! cumulative infiltration for this 'event' - scs_S1(nm) = scs_Se(nm) ! S for this 'event' - scs_rain(nm) = 1 ! logic used to determine if there is an event ongoing - endif - ! - ! Compute cum rainfall - ! - scs_P1(nm) = scs_P1(nm) + prcp(nm)*dt - ! - ! Compute runoff - ! - if (scs_P1(nm) > (sfacinf* scs_S1(nm)) ) then ! scs_S1 is S - Qq = (scs_P1(nm) - (sfacinf*scs_S1(nm)))**2 / (scs_P1(nm) + (1.0 - sfacinf)*scs_S1(nm)) ! cumulative runoff in m - I = scs_P1(nm) - Qq ! cum infiltration this event - qinfmap(nm) = (I - scs_F1(nm))/dt ! infiltration in m/s - scs_F1(nm) = I ! cum infiltration this event - else - Qq = 0.0 ! no runoff - scs_F1(nm) = scs_P1(nm) ! all rainfall is infiltrated - qinfmap(nm) = prcp(nm) ! infiltration rate = rainfall rate - endif - ! - ! Compute "remaining S", but note that scs_Se is not used in computation - ! - scs_Se(nm) = scs_Se(nm) - qinfmap(nm)*dt - scs_Se(nm) = max(scs_Se(nm), 0.0) - qinfmap(nm) = max(qinfmap(nm), 0.0) - ! - else - ! - ! It is not raining here - if (scs_rain(nm) == 1) then - ! - ! if it was raining before; cange logic and set rate to 0 - scs_rain(nm) = 0 - qinfmap(nm) = 0.0 - rain_T1(nm) = 0.0 - ! - endif - ! - ! Add to recovery time - rain_T1(nm) = rain_T1(nm) + dt / 3600 - ! - ! compute recovery of S if time is larger than this - if (rain_T1(nm) > (0.06 / inf_kr(nm)) ) then ! Equation 4-37 from SWMM - ! note that scs_Se is S and qinffield is Smax - scs_Se(nm) = scs_Se(nm) + (inf_kr(nm) * qinffield(nm) * dt / 3600) ! scs_kr is recovery in hours - scs_Se(nm) = min(scs_Se(nm), qinffield(nm)) - ! - endif - ! - endif - ! - ! Compute cumulative values - cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! - enddo - !$omp end do - !$omp end parallel - ! - !$acc update device(qinfmap), async(1) - ! - elseif (inftype == 'gai') then - ! - ! Determine infiltration rate with with the Green-Ampt (GA) model - !$omp parallel & - !$omp private ( nm ) - !$omp do - do nm = 1, np - ! - ! If there is precip in this grid cell for this time step? - if (prcp(nm) > 0.0) then - ! - ! Is raining now - if ( (prcp(nm)) < ksfield(nm)) then - ! - ! Small amounts of rainfall - infiltration is same as soil - qinfmap(nm) = prcp(nm) ! infiltration is same as rainfall - ! - else - ! - ! Larger amounts of rainfall - Equation 4-27 from SWMM manual - qinfmap(nm) = (ksfield(nm) * (1 + (GA_head(np)*GA_sigma(np)) / GA_F(nm))) - qinfmap(nm) = min(qinfmap(nm), prcp(nm)) ! never more than rainfall - qinfmap(nm) = max(qinfmap(nm), 0.0) ! and never negative - ! - endif - ! - ! Update sigma - GA_sigma(nm) = GA_sigma(nm) - (qinfmap(nm)*dt/GA_Lu(nm)) - GA_sigma(nm) = max(GA_sigma(nm),0.00) - ! - ! Update others - GA_F(nm) = GA_F(nm) + qinfmap(nm)*dt ! internal cumulative rainfall from Green-Ampt - rain_T1(nm) = 0.0 ! recovery time not started - ! - else - ! - ! Not raining here - ! - ! Add to recovery time - rain_T1(nm) = rain_T1(nm) + dt / 3600 - ! - ! compute recovery of S if time is larger than this - if (rain_T1(nm) > (0.06 / inf_kr(nm)) ) then ! Equation 4-37 from SWMM - ! - ! Update sigma - GA_sigma(nm) = GA_sigma(nm) + (inf_kr(nm) * GA_sigma_max(nm) * dt/3600) ! Equation 4-35 - GA_sigma(nm) = min(GA_sigma(nm), GA_sigma_max(nm)) ! never more than max - ! - ! Update internal cumulative rainfall - GA_F(nm) = GA_F(nm) - (inf_kr(nm) * GA_sigma_max(nm) * dt/3600*GA_Lu(nm)) ! Page 112 SWMM - GA_F(nm) = max(GA_F(nm), 0.0) ! never negative - endif - endif - ! - ! Compute cumulative values - qinffield(nm) = qinfmap(nm) - cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! - enddo - !$omp end do - !$omp end parallel - ! - !$acc update device(qinfmap), async(1) - ! - elseif (inftype == 'hor') then - ! - ! Determine infiltration rate with with the Horton model - !$omp parallel & - !$omp private ( Qq,I,nm,a,hh_local) - !$omp do - do nm = 1, np - ! - ! Get local water depth estimate - if (subgrid) then - if (crsgeo) then - a = cell_area_m2(nm) - else - a = cell_area(z_flags_iref(nm)) - endif - hh_local = z_volume(nm)/a - else - hh_local = zs(nm) - zb(nm) - endif - ! - ! Check if there is water - if (hh_local> 0.0 .OR. prcp(nm) > 0.0) then - ! - ! Infiltrating here - ! - ! Count how long this is already going - if (rain_T1(nm) > 0) then - rain_T1 = 0.0 - endif - rain_T1(nm) = rain_T1(nm) - dt ! negative amount of how long it is infiltrating - ! - ! Compute estimate of infiltration ! Note that qinffield = horton_fc - I = exp(horton_kd(nm)*rain_T1(nm)/3600) ! note that horton_kd is factor in hours while dt is seconds - ! - ! Stop keeping track of this when less than 1% left (same for time) - if (I < 0.01) then - I = 0.0 ! which reduces qinfmap to horton_fc - rain_T1(nm) = rain_T1(nm) + dt ! also make sure time doesnt further decrease - qinfmap(nm) = (horton_fc(nm))/3600/1000 ! from mm/hr to m/s - else - qinfmap(nm) = (horton_fc(nm) + (horton_f0(nm) - horton_fc(nm))*I)/3600/1000 ! from mm/hr to m/s - endif - ! - ! Check how much there can infiltrate - if (hh_local > 0.00) then - Qq = prcp(nm)*dt + (zs(nm) - zb(nm)) ! Qq is estimate in meter of how much water there is - else - Qq = prcp(nm)*dt ! if no water; only compare with rainfall - endif - ! - ! Compare how much Horton want to infiltrate - I = qinfmap(nm) * dt ! I is estimate in meter of how much Horton allows - if (I > Qq) then - qinfmap(nm) = qinfmap(nm) * Qq/I ! scale Horton if capacity > available - endif - ! - else - ! - ! Not raining here NOR ponding - rain_T1(nm) = rain_T1(nm) + dt/horton_kr_kd ! positive amount of how long it is infiltrating - qinfmap(nm) = 0.0 - ! - endif - ! - ! Compute cumulative values - cuminf(nm) = cuminf(nm) + qinfmap(nm)*dt - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! - enddo - !$omp end do - !$omp end parallel - ! - endif - ! - end subroutine - -end module diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index 138dad69..ef381749 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -179,6 +179,7 @@ subroutine read_sfincs_input() call read_char_input(500,'netamprfile',netamprfile,'none') call read_char_input(500,'netampfile',netampfile,'none') call read_char_input(500,'netspwfile',netspwfile,'none') + call read_char_input(500,'netinfiltrationfile',netinfiltrationfile,'none') ! ! Output call read_char_input(500,'obsfile',obsfile,'none') diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index e2c8028a..1b937b99 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -85,8 +85,8 @@ function sfincs_initialize(config_file) result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = '$Rev: v2.0.6-alpha' - build_date = '$Date: 2024-03-28' + build_revision = '$Rev: v2.0.6-alpha:branch-80-add-option-for-netcdf_infiltration_quadtree-file-input' + build_date = '$Date: 2024-04-25' ! write(*,'(a)')'' write(*,*)'----------- Welcome to SFINCS -----------' diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index 2ebf7df4..e1e58c33 100644 --- a/source/src/sfincs_ncoutput.F90 +++ b/source/src/sfincs_ncoutput.F90 @@ -50,7 +50,7 @@ module sfincs_ncoutput integer :: patm_varid, wind_speed_varid, wind_dir_varid integer :: inp_varid, total_runtime_varid, average_dt_varid, status_varid integer :: hm0_varid, hm0ig_varid, zsm_varid, tp_varid, tpig_varid, wavdir_varid, dirspr_varid - integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, qb_varid, beta_varid, srcsh_varid, alphaig_varid + integer :: dw_varid, df_varid, dwig_varid, dfig_varid, cg_varid, qb_varid, beta_varid, srcig_varid, alphaig_varid ! end type ! @@ -1500,12 +1500,12 @@ subroutine ncoutput_his_init() NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, 'long_name', 'directionally averaged normalised bed slope')) NF90(nf90_put_att(his_file%ncid, his_file%beta_varid, 'coordinates', 'station_id station_name point_x point_y')) ! - NF90(nf90_def_var(his_file%ncid, 'srcsh', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%srcsh_varid)) ! time-varying water level point - NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, '_FillValue', FILL_VALUE)) - NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'units', '-')) - NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'standard_name', 'directionally_averaged_ig_energy_source')) - NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'long_name', 'directionally averaged ig energy source')) - NF90(nf90_put_att(his_file%ncid, his_file%srcsh_varid, 'coordinates', 'station_id station_name point_x point_y')) + NF90(nf90_def_var(his_file%ncid, 'srcig', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%srcig_varid)) ! time-varying water level point + NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'units', '-')) + NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'standard_name', 'directionally_averaged_ig_energy_source')) + NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'long_name', 'directionally averaged ig energy source')) + NF90(nf90_put_att(his_file%ncid, his_file%srcig_varid, 'coordinates', 'station_id station_name point_x point_y')) ! NF90(nf90_def_var(his_file%ncid, 'alphaig', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%alphaig_varid)) ! time-varying water level point NF90(nf90_put_att(his_file%ncid, his_file%alphaig_varid, '_FillValue', FILL_VALUE)) @@ -2224,7 +2224,7 @@ subroutine ncoutput_update_his(t,nthisout) real*4, dimension(nobs) :: cgobs real*4, dimension(nobs) :: qbobs real*4, dimension(nobs) :: betaobs - real*4, dimension(nobs) :: srcshobs + real*4, dimension(nobs) :: srcigobs real*4, dimension(nobs) :: alphaigobs real*4, dimension(:), allocatable :: qq ! @@ -2248,7 +2248,7 @@ subroutine ncoutput_update_his(t,nthisout) cgobs = FILL_VALUE qbobs = FILL_VALUE betaobs = FILL_VALUE - srcshobs = FILL_VALUE + srcigobs = FILL_VALUE alphaigobs = FILL_VALUE ! do iobs = 1, nobs ! determine zs and prcp of obervation points at required timestep @@ -2350,7 +2350,7 @@ subroutine ncoutput_update_his(t,nthisout) cgobs(iobs) = cg(nm) qbobs(iobs) = qb(nm) betaobs(iobs) = betamean(nm) - srcshobs(iobs) = srcsh(nm) + srcigobs(iobs) = srcig(nm) alphaigobs(iobs) = alphaig(nm) ! endif @@ -2414,7 +2414,7 @@ subroutine ncoutput_update_his(t,nthisout) ! NF90(nf90_put_var(his_file%ncid, his_file%qb_varid, qbobs, (/1, nthisout/))) NF90(nf90_put_var(his_file%ncid, his_file%beta_varid, betaobs, (/1, nthisout/))) - NF90(nf90_put_var(his_file%ncid, his_file%srcsh_varid, srcshobs, (/1, nthisout/))) + NF90(nf90_put_var(his_file%ncid, his_file%srcig_varid, srcigobs, (/1, nthisout/))) NF90(nf90_put_var(his_file%ncid, his_file%alphaig_varid, alphaigobs, (/1, nthisout/))) ! endif @@ -2575,7 +2575,7 @@ subroutine ncoutput_update_max(t,ntmaxout) n = z_index_z_n(nm) m = z_index_z_m(nm) ! - zstmp(m, n) = cumprcp(nm)*1000 + zstmp(m, n) = cumprcp(nm) ! enddo ! @@ -2669,17 +2669,19 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) ! nm = index_sfincs_in_quadtree(nmq) ! - if (kcs(nm)>0) then - if (subgrid) then - if ( (zsmax(nm) - subgrid_z_zmin(nm)) > huthresh) then - zstmp(nmq) = zsmax(nm) + if (nm>0) then + if (kcs(nm)>0) then + if (subgrid) then + if ( (zsmax(nm) - subgrid_z_zmin(nm)) > huthresh) then + zstmp(nmq) = zsmax(nm) + endif + else + if ( (zsmax(nm) - zb(nm)) > huthresh) then + zstmp(nmq) = zsmax(nm) + endif endif - else - if ( (zsmax(nm) - zb(nm)) > huthresh) then - zstmp(nmq) = zsmax(nm) - endif endif - endif + endif enddo ! NF90(nf90_put_var(map_file%ncid, map_file%timemax_varid, t, (/ntmaxout/))) ! write time_max @@ -2687,27 +2689,28 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) ! ! Write maximum water depth if (subgrid .eqv. .false. .or. store_hsubgrid .eqv. .true.) then + ! zstmp = FILL_VALUE ! - if (subgrid) then - do nm = 1, np - ! - if ( (zsmax(nm) - subgrid_z_zmin(nm)) > huthresh) then - zstmp(nm) = zsmax(nm) - subgrid_z_zmin(nm) - endif - ! - enddo - else - do nm = 1, np - ! - if ( (zsmax(nm) - zb(nm)) > huthresh) then - zstmp(nm) = zsmax(nm) - zb(nm) - endif - enddo - endif - endif - ! - if (subgrid .eqv. .false. .or. store_hsubgrid .eqv. .true.) then + do nmq = 1, quadtree_nr_points + ! + nm = index_sfincs_in_quadtree(nmq) + ! + if (nm>0) then + if (kcs(nm)>0) then + if (subgrid) then + if ( (zsmax(nm) - subgrid_z_zmin(nm)) > huthresh) then + zstmp(nmq) = zsmax(nm) - subgrid_z_zmin(nm) + endif + else + if ( (zsmax(nm) - zb(nm)) > huthresh) then + zstmp(nmq) = zsmax(nm) - zb(nm) + endif + endif + endif + endif + enddo + ! NF90(nf90_put_var(map_file%ncid, map_file%hmax_varid, zstmp, (/1, ntmaxout/))) ! write hmax endif ! @@ -2719,9 +2722,11 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) zstmp = FILL_VALUE do nmq = 1, quadtree_nr_points nm = index_sfincs_in_quadtree(nmq) - if (kcs(nm)>0) then - zstmp(nmq) = cumprcp(nm)*1000 - endif + if (nm>0) then + if (kcs(nm)>0) then + zstmp(nmq) = cumprcp(nm) + endif + endif enddo NF90(nf90_put_var(map_file%ncid, map_file%cumprcp_varid, zstmp, (/1, ntmaxout/))) ! write cumprcp ! @@ -2730,9 +2735,11 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) zstmp = FILL_VALUE do nmq = 1, quadtree_nr_points nm = index_sfincs_in_quadtree(nmq) - if (kcs(nm)>0) then - zstmp(nmq) = cuminf(nm) - endif + if (nm>0) then + if (kcs(nm)>0) then + zstmp(nmq) = cuminf(nm) + endif + endif enddo NF90(nf90_put_var(map_file%ncid, map_file%cuminf_varid, zstmp, (/1, ntmaxout/))) ! write cuminf ! @@ -2743,9 +2750,11 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) zstmp = FILL_VALUE do nmq = 1, quadtree_nr_points nm = index_sfincs_in_quadtree(nmq) - if (kcs(nm)>0) then - zstmp(nmq) = vmax(nm) - endif + if (nm>0) then + if (kcs(nm)>0) then + zstmp(nmq) = vmax(nm) + endif + endif enddo NF90(nf90_put_var(map_file%ncid, map_file%vmax_varid, zstmp, (/1, ntmaxout/))) ! write vmax endif @@ -2755,9 +2764,11 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) zstmp = FILL_VALUE do nmq = 1, quadtree_nr_points nm = index_sfincs_in_quadtree(nmq) - if (kcs(nm)>0) then - zstmp(nmq) = qmax(nm) - endif + if (nm>0) then + if (kcs(nm)>0) then + zstmp(nmq) = qmax(nm) + endif + endif enddo NF90(nf90_put_var(map_file%ncid, map_file%qmax_varid, zstmp, (/1, ntmaxout/))) ! write qmax endif @@ -2767,8 +2778,10 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) zstmp = FILL_VALUE do nmq = 1, quadtree_nr_points nm = index_sfincs_in_quadtree(nmq) - if (kcs(nm)>0) then - zstmp(nmq) = twet(nm) + if (nm>0) then + if (kcs(nm)>0) then + zstmp(nmq) = twet(nm) + endif endif enddo NF90(nf90_put_var(map_file%ncid, map_file%tmax_varid, zstmp, (/1, ntmaxout/))) ! write tmax @@ -2779,24 +2792,14 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) zstmp = FILL_VALUE do nmq = 1, quadtree_nr_points nm = index_sfincs_in_quadtree(nmq) - if (kcs(nm)>0) then - zstmp(nmq) = windmax(nm) + if (nm>0) then + if (kcs(nm)>0) then + zstmp(nmq) = windmax(nm) + endif endif enddo NF90(nf90_put_var(map_file%ncid, map_file%windmax_varid, zstmp, (/1, ntmaxout/))) ! write windmax endif - ! - ! Cumulative infiltration - if (infiltration) then - zstmp = FILL_VALUE - do nmq = 1, quadtree_nr_points - nm = index_sfincs_in_quadtree(nmq) - if (kcs(nm)>0) then - zstmp(nmq) = cuminf(nm) - endif - enddo - NF90(nf90_put_var(map_file%ncid, map_file%cuminf_varid, zstmp, (/1, ntmaxout/))) ! write cuminf - endif ! end subroutine ! diff --git a/source/src/sfincs_obspoints.f90 b/source/src/sfincs_obspoints.f90 index b547bac3..4882b78f 100644 --- a/source/src/sfincs_obspoints.f90 +++ b/source/src/sfincs_obspoints.f90 @@ -107,20 +107,24 @@ subroutine read_obs_points() if (nmq>0) then ! nm = index_sfincs_in_quadtree(nmq) - nmindobs(iobs) = nm - n = z_index_z_n(nm) - m = z_index_z_m(nm) ! - xgobs(iobs) = z_xz(nm) - ygobs(iobs) = z_yz(nm) - ! - if (subgrid) then - zbobs(iobs) = subgrid_z_zmin(nm) - else - zbobs(iobs) = zb(nm) - endif - ! - iref = z_flags_iref(nm) + if (nm > 0) then + ! + nmindobs(iobs) = nm + n = z_index_z_n(nm) + m = z_index_z_m(nm) + ! + xgobs(iobs) = z_xz(nm) + ygobs(iobs) = z_yz(nm) + ! + if (subgrid) then + zbobs(iobs) = subgrid_z_zmin(nm) + else + zbobs(iobs) = zb(nm) + endif + ! + iref = z_flags_iref(nm) + endif ! ! write(*,'(a,i0,a,a,a,i0,a,i0,a,i0,a,i0,a,f0.3)')' Observation point ',iobs,' : "',trim(nameobs(iobs)),'" nm=',nm,' n=',n,' m=',m,' iref=',iref,' z=',zbobs(iobs) ! diff --git a/source/src/sfincs_snapwave.f90 b/source/src/sfincs_snapwave.f90 index 7a0e997e..6a624738 100644 --- a/source/src/sfincs_snapwave.f90 +++ b/source/src/sfincs_snapwave.f90 @@ -22,7 +22,7 @@ module sfincs_snapwave real*4, dimension(:), allocatable :: snapwave_cg real*4, dimension(:), allocatable :: snapwave_Qb real*4, dimension(:), allocatable :: snapwave_beta - real*4, dimension(:), allocatable :: snapwave_srcsh + real*4, dimension(:), allocatable :: snapwave_srcig real*4, dimension(:), allocatable :: snapwave_alphaig integer, dimension(:,:), allocatable :: snapwave_connected_nodes integer*4, dimension(:), allocatable :: index_snapwave_in_sfincs @@ -140,7 +140,7 @@ subroutine update_wave_field(t, tloop) real*4, dimension(:), allocatable :: cg0 real*4, dimension(:), allocatable :: qb0 real*4, dimension(:), allocatable :: beta0 - real*4, dimension(:), allocatable :: srcsh0 + real*4, dimension(:), allocatable :: srcig0 real*4, dimension(:), allocatable :: alphaig0 integer :: ip, ii, m, n, nm, nmu, idir real*4 :: f @@ -157,7 +157,7 @@ subroutine update_wave_field(t, tloop) allocate(cg0(np)) allocate(qb0(np)) allocate(beta0(np)) - allocate(srcsh0(np)) + allocate(srcig0(np)) allocate(alphaig0(np)) ! fwx0 = 0.0 @@ -169,7 +169,7 @@ subroutine update_wave_field(t, tloop) cg0 = 0.0 qb0 = 0.0 beta0 = 0.0 - srcsh0 = 0.0 + srcig0 = 0.0 alphaig0 = 0.0 ! ! Determine SnapWave water depth @@ -221,7 +221,7 @@ subroutine update_wave_field(t, tloop) cg0(nm) = snapwave_cg(ip) qb0(nm) = snapwave_Qb(ip) beta0(nm) = snapwave_beta(ip) - srcsh0(nm) = snapwave_srcsh(ip) + srcig0(nm) = snapwave_srcig(ip) alphaig0(nm) = snapwave_alphaig(ip) if (store_wave_direction) then mean_wave_direction(nm) = 270.0 - snapwave_mean_direction(ip)*180/pi @@ -243,7 +243,7 @@ subroutine update_wave_field(t, tloop) cg0(nm) = 0.0 qb0(nm) = 0.0 beta0(nm) = 0.0 - srcsh0(nm) = 0.0 + srcig0(nm) = 0.0 alphaig0(nm) = 0.0 if (store_wave_direction) then mean_wave_direction(nm) = 0.0 @@ -263,7 +263,7 @@ subroutine update_wave_field(t, tloop) cg(nm) = cg0(nm) qb(nm) = qb0(nm) betamean(nm) = beta0(nm) - srcsh(nm) = srcsh0(nm) + srcig(nm) = srcig0(nm) alphaig(nm) = alphaig0(nm) ! endif @@ -338,7 +338,7 @@ subroutine compute_snapwave(t) snapwave_cg = cg snapwave_Qb = Qb snapwave_beta = beta - snapwave_srcsh = srcsh + snapwave_srcig = srcig snapwave_alphaig = alphaig ! ! Wave periods from SnapWave, used in e.g. wavemakers - TL: moved behind call update_boundary_conditions & compute_wave_field so values at first timestep are not 0 @@ -399,6 +399,7 @@ subroutine read_snapwave_input() call read_real_input(500,'snapwave_alphaigfac',alphaigfac,1.0) ! Multiplication factor for IG shoaling source/sink term call read_real_input(500,'snapwave_baldock_ratio_ig',baldock_ratio_ig,0.2) call read_int_input(500,'snapwave_ig_opt',ig_opt,1) + call read_int_input(500,'snapwave_iterative_srcig',iterative_srcig,1) ! Option whether to calculate IG source/sink term in iterative lower (better, but potentially slower, 1=default), or effectively based on previous timestep (faster, potential mismatch, =0) ! ! IG boundary conditions options: call read_int_input(500,'snapwave_use_herbers',herbers_opt,1) ! Choice whether you want IG Hm0&Tp be calculated by herbers (=1, default), or want to specify user defined values (0> then snapwave_eeinc2ig & snapwave_Tinc2ig are used) @@ -416,7 +417,7 @@ subroutine read_snapwave_input() call read_char_input(500,'snapwave_btpfile',btpfile,'') call read_char_input(500,'snapwave_bwdfile',bwdfile,'') call read_char_input(500,'snapwave_bdsfile',bdsfile,'') - call read_char_input(500,'snapwave_upwfile',upwfile,'') + call read_char_input(500,'snapwave_upwfile',upwfile,'snapwave.upw') call read_char_input(500,'snapwave_mskfile',mskfile,'') call read_char_input(500,'snapwave_depfile',depfile,'none') call read_char_input(500,'snapwave_ncfile', gridfile,'snapwave_net.nc') diff --git a/source/src/sfincs_wavemaker.f90 b/source/src/sfincs_wavemaker.f90 index aeb7bff2..d08ca7e5 100644 --- a/source/src/sfincs_wavemaker.f90 +++ b/source/src/sfincs_wavemaker.f90 @@ -1467,13 +1467,13 @@ subroutine update_wavemaker_fluxes(t, dt, tloop) ! zsuv = max(zsnmb, zsnmi) ! - if (zsuv>=subgrid_uv_zmax(nm) - 1.0e-3) then + if (zsuv>=subgrid_uv_zmax(ip) - 1.0e-3) then ! ! Entire cell is wet, no interpolation from table needed ! depthuv = subgrid_uv_havg_zmax(ip) + zsuv ! - elseif (zsuv>subgrid_uv_zmin(nm)) then + elseif (zsuv>subgrid_uv_zmin(ip)) then ! ! Interpolation required ! diff --git a/source/src/snapwave/snapwave_data.f90 b/source/src/snapwave/snapwave_data.f90 index 5f8e5deb..7659715d 100644 --- a/source/src/snapwave/snapwave_data.f90 +++ b/source/src/snapwave/snapwave_data.f90 @@ -11,12 +11,9 @@ module snapwave_data real*4, dimension(:), allocatable :: kwav_ig, nwav_ig ! wave number, ratio Cg/C real*4, dimension(:), allocatable :: C, Cg ! wave celerity, group velocity0 real*4, dimension(:), allocatable :: C_ig, Cg_ig ! wave celerity, group velocity0 - real*4, dimension(:,:), allocatable :: Sxx ! directional radiation stress Sxx real*4, dimension(:), allocatable :: fw ! friction coefficient real*4, dimension(:), allocatable :: fw_ig ! friction coefficient real*4, dimension(:), allocatable :: H, H_ig ! rms wave height - real*4, dimension(:), allocatable :: H_ig_old ! IG wave height at previous timestep - real*4, dimension(:), allocatable :: H_inc_old ! Incident wave height at previous timestep real*4, dimension(:), allocatable :: Dw,Df ! dissipation due to breaking, bed friction real*4, dimension(:), allocatable :: Dw_ig,Df_ig ! dissipation due to breaking, bed friction for IG real*4, dimension(:), allocatable :: F ! wave force Dw/C/rho/depth @@ -59,7 +56,7 @@ module snapwave_data ! real*4, dimension(:), allocatable :: Qb real*4, dimension(:), allocatable :: beta - real*4, dimension(:), allocatable :: srcsh + real*4, dimension(:), allocatable :: srcig real*4, dimension(:), allocatable :: alphaig ! integer*4, dimension(:), allocatable :: index_snapwave_in_quadtree @@ -162,6 +159,8 @@ module snapwave_data real*4 :: baldock_ratio_ig ! option controlling from what depth wave breaking should take place for IG waves: (Hk>baldock_ratio*Hmx(k)), default baldock_ratio_ig=0.2 integer :: igwaves_opt ! option of IG waves on (1) or off (0) integer :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) + integer :: iterative_srcig ! option whether IG source/sink term should be calculated in the iterative loop again (iterative_srcig = 1, is a bit slower), + ! ... or just a priori based on effectively incident wave energy from previous timestep only real*4 :: snapwave_alpha_ig,gamma_ig ! coefficients in Baldock wave breaking dissipation model for IG waves 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 diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90 index 484eac1a..8d864cec 100644 --- a/source/src/snapwave/snapwave_domain.f90 +++ b/source/src/snapwave/snapwave_domain.f90 @@ -114,13 +114,10 @@ subroutine initialize_snapwave_domain() allocate(nwav_ig(no_nodes)) allocate(C_ig(no_nodes)) allocate(Cg_ig(no_nodes)) - allocate(Sxx(ntheta,no_nodes)) allocate(sinhkh_ig(no_nodes)) allocate(Hmx_ig(no_nodes)) allocate(fw_ig(no_nodes)) - allocate(H_ig(no_nodes)) - allocate(H_ig_old(no_nodes)) - allocate(H_inc_old(no_nodes)) + allocate(H_ig(no_nodes)) allocate(Dw(no_nodes)) allocate(Dw_ig(no_nodes)) allocate(F(no_nodes)) @@ -131,7 +128,7 @@ subroutine initialize_snapwave_domain() allocate(thetam(no_nodes)) allocate(Qb(no_nodes)) allocate(beta(no_nodes)) - allocate(srcsh(no_nodes)) + allocate(srcig(no_nodes)) allocate(alphaig(no_nodes)) ! allocate(uorb(no_nodes)) allocate(ctheta(ntheta,no_nodes)) @@ -196,7 +193,8 @@ subroutine initialize_snapwave_domain() ds360d0 = 0.d0 w360d0 = 0.d0 prev360 = 0 - H_ig_old = 0.0 + H = 0.0 + H_ig = 0.0 ! generate_upw = .false. exists = .true. @@ -274,6 +272,8 @@ subroutine initialize_snapwave_domain() ! endif ! + nb = 0 + ! do k=1,no_nodes !if (msk(k)==3) msk(k) = 1 ! Set outflow points to regular points > now should become neumann, so don't do this do itheta=1,ntheta360 @@ -282,25 +282,27 @@ subroutine initialize_snapwave_domain() if (inout>0) msk(k) = 2 endif enddo - enddo - ! - nb = 0 - do k = 1, no_nodes - if (msk(k)==1) then - inner(k) = .true. - else - inner(k) = .false. + ! + if (msk(k)>1) then + nb = nb + 1 endif - if (msk(k)==2) nb = nb + 1 + ! enddo ! allocate(nmindbnd(nb)) ! - nb = 0 do k = 1, no_nodes + ! + if (msk(k)==1) then + inner(k) = .true. + else + inner(k) = .false. + endif + ! if (msk(k)>1) then - nb = nb + 1 - nmindbnd(nb) = k + ! + nmindbnd(nb) = k + ! endif enddo ! diff --git a/source/src/snapwave/snapwave_infragravity.f90 b/source/src/snapwave/snapwave_infragravity.f90 index 6fbc48da..1df9a55d 100644 --- a/source/src/snapwave/snapwave_infragravity.f90 +++ b/source/src/snapwave/snapwave_infragravity.f90 @@ -21,8 +21,9 @@ subroutine determine_ig_bc(hsinc, tpinc, ds, jonswapgam, depth, Tinc2ig, tpig_op ! implicit none ! - real*4, intent(in) :: hsinc, tpinc, ds, jonswapgam, depth, Tinc2ig + real*4, intent(in) :: hsinc, tpinc, ds, jonswapgam, Tinc2ig integer, intent(in) :: tpig_opt + real*4, intent(inout) :: depth real*4, intent(out) :: hsig, tpig ! real*4 :: pi, scoeff @@ -40,8 +41,10 @@ subroutine determine_ig_bc(hsinc, tpinc, ds, jonswapgam, depth, Tinc2ig, tpig_op ! Call function that calculates Hig0 following Herbers, as also implemented in XBeach and secordspec2 in Matlab ! Loosely based on 3 step calculation in waveparams.F90 of XBeach (build_jonswap, build_etdir, build_boundw), here all in 1 subroutine calculate_herbers ! - if (depth < 10.0) then - write(*,*)'DEBUG SnapWave - depth at boundary input point dropped below 10 m: ',depth, ' which might lead to large values of Hm0ig as bc, especially when directional spreading is low!' + if (depth < 5.0) then + write(*,*)'ERROR SnapWave - depth at boundary input point dropped below 5 m: ',depth, ' which might lead to large values of Hm0ig as bc, especially when directional spreading is low! Please specify input in deeper water. ' + write(*,*)'Depth set back to 5 meters for stability, simulation will continue.' + depth = 5.0 endif ! call compute_herbers(hsig, Tm01, Tm10, Tp, Tpsmooth, hsinc, tpinc, scoeff, jonswapgam, depth, correctHm0) ![out,out,out,out,out, in,in,in,in,in,in] diff --git a/source/src/snapwave/snapwave_solver.f90 b/source/src/snapwave/snapwave_solver.f90 index 606a1335..00745eda 100644 --- a/source/src/snapwave/snapwave_solver.f90 +++ b/source/src/snapwave/snapwave_solver.f90 @@ -39,7 +39,6 @@ subroutine compute_wave_field() call disper_approx(depth, Tpb, kwav, nwav, C, Cg, no_nodes) call disper_approx(depth, Tpb_ig, kwav_ig, nwav_ig, C_ig, Cg_ig, no_nodes) cg_ig = cg - Sxx = 0 ! do k = 1, no_nodes sinhkh(k) = sinh(min(kwav(k)*depth(k), 100.0)) @@ -47,11 +46,6 @@ subroutine compute_wave_field() sinhkh_ig(k) = sinh(min(kwav_ig(k)*depth(k), 50.0)) Hmx_ig(k) = 0.88/kwav_ig(k)*tanh(gamma_ig*kwav_ig(k)*depth(k)/0.88) ! - ! Calculate radiation stress Sxx = ((2 .* n) - 0.5) .* Einc > is from energy of previous timestep, and directional! - do itheta = 1, ntheta - ! old, non-directional: Sxx(itheta,k) = ((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5) * sum(ee(:, k)) * dtheta - Sxx(itheta,k) = ((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5) * ee(itheta, k) ! limit so value of nwav is between 0 and 1, and Sxx therefore doesn't become NaN for nwav=Infinite - enddo enddo ! do itheta = 1, ntheta @@ -90,17 +84,13 @@ subroutine compute_wave_field() call solve_energy_balance2Dstat (x,y,no_nodes,w,ds,inner,prev,neumannconnected, & theta,ntheta,thetamean, & depth,zb,kwav,kwav_ig,cg,cg_ig,ctheta,ctheta_ig,fw,fw_ig,Tpb,Tpb_ig,50000.,rho,snapwave_alpha,snapwave_alpha_ig,gamma,& - H,H_ig,Dw,Dw_ig,F,Df,Df_ig,thetam,sinhkh,sinhkh_ig,Hmx,Hmx_ig, ee, ee_ig, igwaves, nr_sweeps, crit, hmin, gamma_ig, shinc2ig, ig_opt, baldock_opt, baldock_ratio, baldock_ratio_ig, alphaigfac, Qb, beta, srcsh, alphaig, Sxx, H_ig_old, H_inc_old, nwav) + H,H_ig,Dw,Dw_ig,F,Df,Df_ig,thetam,sinhkh,sinhkh_ig,Hmx,Hmx_ig, ee, ee_ig, igwaves, nr_sweeps, crit, hmin, gamma_ig, shinc2ig, ig_opt, iterative_srcig, baldock_opt, baldock_ratio, baldock_ratio_ig, alphaigfac, Qb, beta, srcig, alphaig, nwav) ! call timer(t3) ! Fx = F*cos(thetam) Fy = F*sin(thetam) ! - ! Incident and IG wave height after solving energy balance - H_inc_old = H - H_ig_old = H_ig - ! write(*,*)'Computation SnapWave timestep took: ', t3 - t2, ' seconds' ! end subroutine @@ -111,8 +101,8 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec theta,ntheta,thetamean, & depth,zb,kwav,kwav_ig,cg,cg_ig,ctheta,ctheta_ig,fw,fw_ig,T,T_ig,dt,rho,alfa,alfa_ig,gamma, & H,H_ig,Dw,Dw_ig,F,Df,Df_ig,thetam,sinhkh,sinhkh_ig,Hmx,Hmx_ig, ee, ee_ig, igwaves, nr_sweeps, crit, hmin, & - gamma_ig, shinc2ig, ig_opt, baldock_opt, baldock_ratio, baldock_ratio_ig, alphaigfac, Qb, betamean, srcsh, & - alphaig, Sxx, H_ig_old, H_inc_old, nwav) + gamma_ig, shinc2ig, ig_opt, iterative_srcig, baldock_opt, baldock_ratio, baldock_ratio_ig, alphaigfac, Qb, betamean, srcig, & + alphaig, nwav) ! implicit none ! @@ -136,7 +126,6 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec real*4, dimension(ntheta,no_nodes), intent(in) :: ctheta ! refractioon speed real*4, dimension(no_nodes), intent(in) :: cg_ig ! group velocity real*4, dimension(no_nodes), intent(in) :: nwav ! wave number n - real*4, dimension(ntheta,no_nodes), intent(in) :: Sxx ! Radiation Stress real*4, dimension(ntheta,no_nodes), intent(inout) :: ee ! real*4, dimension(ntheta,no_nodes), intent(inout) :: ee_ig ! real*4, dimension(ntheta,no_nodes), intent(in) :: ctheta_ig ! refractioon speed @@ -144,7 +133,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec real*4, dimension(no_nodes), intent(in) :: fw_ig ! wave friction factor real*4, dimension(no_nodes), intent(out) :: Qb ! Fraction of breaking waves according to Baldock's formulation real*4, dimension(no_nodes), intent(out) :: betamean ! Mean local bed slope parameter - real*4, dimension(no_nodes), intent(out) :: srcsh ! Directionally averaged incident wave sink/infragravity source term + real*4, dimension(no_nodes), intent(out) :: srcig ! Directionally averaged incident wave sink/infragravity source term real*4, dimension(no_nodes), intent(out) :: alphaig ! Mean IG shoaling parameter alpha real*4, intent(in) :: T ! wave period real*4, intent(in) :: T_ig ! IG wave period @@ -155,8 +144,10 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec real*4, intent(in) :: baldock_ratio ! option controlling from what depth wave breaking should take place: (Hk>baldock_ratio*Hmx(k)), default baldock_ratio=0.2 real*4, intent(in) :: baldock_ratio_ig ! option controlling from what depth wave breaking should take place for IG waves: (Hk>baldock_ratio*Hmx(k)), default baldock_ratio_ig=0.2 integer :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking) - real*4, dimension(no_nodes), intent(out) :: H ! wave height - real*4, dimension(no_nodes), intent(out) :: H_ig ! wave height + integer :: iterative_srcig ! option whether IG source/sink term should be calculated in the iterative loop again (iterative_srcig = 1, is a bit slower), + ! ... or just a priori based on effectively incident wave energy from previous timestep only + real*4, dimension(no_nodes), intent(inout) :: H ! wave height + real*4, dimension(no_nodes), intent(inout) :: H_ig ! wave height real*4, dimension(no_nodes), intent(out) :: Dw ! wave breaking dissipation real*4, dimension(no_nodes), intent(out) :: Dw_ig ! wave breaking dissipation IG real*4, dimension(no_nodes), intent(out) :: F ! wave force Dw/C/rho/h @@ -167,8 +158,6 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec real*4, dimension(no_nodes), intent(in) :: Hmx ! Hmax real*4, dimension(no_nodes), intent(in) :: sinhkh_ig ! sinh(k*depth) real*4, dimension(no_nodes), intent(in) :: Hmx_ig ! Hmax - real*4, dimension(no_nodes), intent(in) :: H_ig_old ! wave height of previous timestep - real*4, dimension(no_nodes), intent(in) :: H_inc_old ! wave height of previous timestep ! ! Local variables and arrays ! @@ -179,20 +168,13 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec integer :: k,k1,k2,k12,i,ind(1),count,kn,itheta ! counters (k is grid index) integer :: indint integer, dimension(:,:), allocatable :: indx ! index for grid sorted per sweep direction -! real*4, dimension(:,:), allocatable :: ee,eeold ! wave energy density, energy density previous iteration -! real*4, dimension(:,:), allocatable :: ee_ig ! wave energy density - real*4, dimension(:,:), allocatable :: eeold ! wave energy density, energy density previous iteration - real*4, dimension(:,:), allocatable :: srcsh_local ! Energy source/sink term because of IG wave shoaling + real*4, dimension(:,:), allocatable :: eeold ! wave energy density, energy density previous iteration + real*4, dimension(:,:), allocatable :: srcig_local ! Energy source/sink term because of IG wave shoaling real*4, dimension(:,:), allocatable :: beta_local ! Local bed slope based on bed level per direction - real*4, dimension(:,:), allocatable :: alphaig_local ! + real*4, dimension(:,:), allocatable :: alphaig_local ! Local infragravity wave shoaling parameter alpha real*4, dimension(:), allocatable :: dee ! difference with energy previous iteration real*4, dimension(:), allocatable :: eeprev, cgprev ! energy density and group velocity at upwind intersection point real*4, dimension(:), allocatable :: eeprev_ig, cgprev_ig ! energy density and group velocity at upwind intersection point - real*4, dimension(:), allocatable :: Sxxprev ! radiation stress at upwind intersection point - real*4, dimension(:), allocatable :: Hprev ! Incident wave height at upwind intersection point - real*4, dimension(:), allocatable :: H_igprev ! IG wave height at upwind intersection point - real*4, dimension(:), allocatable :: H_incprev ! Incident wave height at upwind intersection point - real*4, dimension(:), allocatable :: depthprev ! water depth at upwind intersection point real*4, dimension(:), allocatable :: A,B,C,R ! 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 @@ -230,8 +212,6 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec real*4 :: Sxx_cons real*4 :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0 real*4 :: sigm_ig - real*4 :: beta - real*4 :: gam ! local gamma (Hinc / depth ratio) character*20 :: fname integer, save :: callno=1 ! @@ -261,7 +241,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec allocate(indx(no_nodes,nr_sweeps)) ! allocate(ee(ntheta,no_nodes)) ! allocate(ee_ig(ntheta,no_nodes)) - allocate(eeold(ntheta,no_nodes)) + allocate(eeold(ntheta,no_nodes)) allocate(dee(ntheta)) allocate(eeprev(ntheta)) allocate(cgprev(ntheta)) @@ -273,7 +253,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec allocate(E(no_nodes)) allocate(diff(no_nodes)) allocate(ra(no_nodes)) - allocate(srcsh_local(ntheta,no_nodes)) + allocate(srcig_local(ntheta,no_nodes)) ! if (igwaves) then ! @@ -287,11 +267,6 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec allocate(E_ig(no_nodes)) allocate(beta_local(ntheta,no_nodes)) allocate(alphaig_local(ntheta,no_nodes)) - allocate(Sxxprev(ntheta)) - allocate(Hprev(ntheta)) - allocate(H_igprev(ntheta)) - allocate(H_incprev(ntheta)) - allocate(depthprev(ntheta)) ! endif ! @@ -305,7 +280,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec df = 0.0 dw = 0.0 F = 0.0 - srcsh_local = 0.0 + srcig_local = 0.0 alphaig_local = 0.0 ! ok = 0 @@ -318,8 +293,6 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec oneover2dtheta = 1.0/2.0/dtheta rhog8 = 0.125*rho*g thetam = 0.0 - H = 0.0 - H_ig = 0.0 ! !gamma_ig = gamma --> now user defineable ! @@ -367,83 +340,15 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec ! if (igwaves) then ! - ! Compute exchange source term inc to ig waves - per direction - do itheta = 1, ntheta - ! - k1 = prev(1, itheta, k) - k2 = prev(2, itheta, k) - ! - if (k1>0 .and. k2>0) then ! IMPORTANT - for some reason (k1*k2)>0 is not reliable always, resulting in directions being uncorrectly skipped!!! - ! - ! Calculate upwind direction dependent variables - ! TL - Note: cg_ig = cg - cgprev(itheta) = w(1, itheta, k)*cg_ig(k1) + w(2, itheta, k)*cg_ig(k2) - Sxxprev(itheta) = w(1, itheta, k)*Sxx(itheta,k1) + w(2, itheta, k)*Sxx(itheta,k2) - Hprev(itheta) = w(1, itheta, k)*H(k1) + w(2, itheta, k)*H(k2) - H_igprev(itheta) = w(1, itheta, k)*H_ig_old(k1) + w(2, itheta, k)*H_ig_old(k2) - H_incprev(itheta) = w(1, itheta, k)*H_inc_old(k1) + w(2, itheta, k)*H_inc_old(k2) - depthprev(itheta) = w(1, itheta, k)*depth(k1) + w(2, itheta, k)*depth(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) - ! - beta = max((w(1, itheta, k)*(zb(k) - zb(k1)) + w(2, itheta, k)*(zb(k) - zb(k2)))/ds(itheta, k), 0.0) - ! Notes: - ! - use actual bed level now for slope, because depth changes because of wave setup/tide/surge - ! - in zb, depth is negative > therefore zb(k) minus zb(k1) - ! - beta=0 means a horizontal or decreasing slope > need alphaig=0 then - ! - beta_local(itheta,k) = beta - !betan_local(itheta,k) = (beta/sigm_ig)*sqrt(9.81/max(depth(k), hmin)) ! TL: in case in the future we would need the normalised bed slope again - ! - gam = max(0.5*(H_incprev(itheta)/depthprev(itheta) + H_inc_old(k)/depth(k)), 0.0) ! mean gamma over current and upwind point - ! - if (ig_opt == 1 .or. ig_opt == 2) then - ! - ! Calculate shoaling parameter alpha_ig following Leijnse et al. (2024) - ! - call estimate_shoaling_parameter_alphaig(beta, gam, alphaig_local(itheta,k)) ! [input, input, output] - ! - ! Now calculate source term component - ! - ! Newest dSxx/dx based method, using estimate of Sxx(k) using conservative shoaling - if (Sxxprev(itheta)<=0.0) then - ! - srcsh_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) - ! - srcsh_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(eeprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta) * dSxx / ds(itheta, k) - ! - endif - ! - else ! TL: option to add future parameterisations here for e.g. coral reef type coasts - ! - srcsh_local(itheta, k) = 0.0 - ! - endif - ! - srcsh_local(itheta, k) = max(srcsh_local(itheta, k), 0.0) - ! - endif - ! - enddo + ! Determine IG source/sink term as in Leijnse, van Ormondt, van Dongeren, Aerts & Muis et al. 2024 + ! + call determine_infragravity_source_sink_term(no_nodes, ntheta, k, w, ds, prev, cg_ig, nwav, depth, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, zb, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) + ! inout: alphaig_local, beta_local, srcig_local - eeprev, eeprev_ig, cgprev + ! in: the rest + ! + ! NOTE - this is now only used to fill 'srcig_local', so that we have a good starting point for the source term in incident wave energy balance in the 'do iter=1,niter' loop + ! THis is therefore based on the energy in the precious SnapWave timestep. + ! The real source term is calculated now in the iteration loop ! endif ! @@ -471,7 +376,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec ! DoverE_ig(k) = (Dwk_ig + Dfk_ig)/max(Ek_ig, 1.0e-6) ! - endif + endif ! endif ! @@ -524,7 +429,6 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec if (igwaves) then !TL: now calculated double? 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) !TL: needed to calculate? Or is same as cgprev(itheta)? - !depthprev(itheta) = w(1, itheta, k)*depth(k1) + w(2, itheta, k)*depth(k2) endif ! enddo @@ -534,7 +438,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec A(itheta) = -ctheta(itheta - 1, k)*oneover2dtheta B(itheta) = oneoverdt + cg(k)/ds(itheta,k) + DoverE(k) C(itheta) = ctheta(itheta + 1, k)*oneover2dtheta - R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta)*eeprev(itheta)/ds(itheta, k) - srcsh_local(itheta, k) * shinc2ig + R(itheta) = oneoverdt*ee(itheta, k) + cgprev(itheta)*eeprev(itheta)/ds(itheta, k) - srcig_local(itheta, k) * shinc2ig ! enddo ! @@ -542,7 +446,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec A(1) = 0.0 B(1) = oneoverdt - ctheta(1, k)/dtheta + cg(k)/ds(1, k) + DoverE(k) C(1) = ctheta(2, k)/dtheta - R(1) = oneoverdt*ee(1, k) + cgprev(1)*eeprev(1)/ds(1, k) - srcsh_local(1, k) * shinc2ig + R(1) = oneoverdt*ee(1, k) + cgprev(1)*eeprev(1)/ds(1, k) - srcig_local(1, k) * shinc2ig else A(1) = 0.0 B(1) = 1.0/dt @@ -554,7 +458,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec A(ntheta) = -ctheta(ntheta - 1, k)/dtheta B(ntheta) = oneoverdt + ctheta(ntheta, k)/dtheta + cg(k)/ds(ntheta, k) + DoverE(k) C(ntheta) = 0.0 - R(ntheta) = oneoverdt*ee(ntheta,k) + cgprev(ntheta)*eeprev(ntheta)/ds(ntheta, k) - srcsh_local(ntheta, k) * shinc2ig + R(ntheta) = oneoverdt*ee(ntheta,k) + cgprev(ntheta)*eeprev(ntheta)/ds(ntheta, k) - srcig_local(ntheta, k) * shinc2ig else A(ntheta) = 0.0 B(ntheta) = oneoverdt @@ -618,6 +522,39 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec ! endif ! + ! Update IG source/sink term based on just updated incident wave energy balance + ! + if (igwaves) then + ! + if (iterative_srcig == 1) then + ! + ! Update H in first iteration - first sweep - ONLY: + if ((iter == 1) .and. (sweep==1)) then + ! + H(k) = sqrt(8*sum(ee(:, k))*dtheta/rho/g) ! is combined: E(k) = sum(ee(:, k))*dtheta And H(k) = sqrt(8*E(k)/rho/g) + ! + endif + ! + ! Compute exchange source term inc to ig waves - per direction + ! + ! Do only on first sweep to save computation time: + if (sweep==1) then + ! Determine IG source/sink term as in Leijnse, van Ormondt, van Dongeren, Aerts & Muis et al. 2024 + ! + call determine_infragravity_source_sink_term(no_nodes, ntheta, k, w, ds, prev, cg_ig, nwav, depth, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, zb, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) + ! inout: alphaig_local, beta_local, srcig_local - eeprev, eeprev_ig, cgprev + ! in: the rest + ! + ! NOTE - this is now only used to fill 'srcig_local', so that we have a good starting point for the source term in incident wave energy balance in the 'do iter=1,niter' loop + ! THis is therefore based on the energy in the precious SnapWave timestep. + ! The real source term is calculated now in the iteration loop + ! + endif + ! + endif + ! + endif + ! if (igwaves) then ! ! IG @@ -627,7 +564,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec 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 - R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta)*eeprev_ig(itheta)/ds(itheta, k) + srcsh_local(itheta, k) + R_ig(itheta) = oneoverdt*ee_ig(itheta, k) + cgprev_ig(itheta)*eeprev_ig(itheta)/ds(itheta, k) + srcig_local(itheta, k) ! enddo ! @@ -635,7 +572,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec 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 - R_ig(1) = oneoverdt*ee_ig(1, k) + cgprev_ig(1)*eeprev_ig(1)/ds(1, k) + srcsh_local(1, k) + R_ig(1) = oneoverdt*ee_ig(1, k) + cgprev_ig(1)*eeprev_ig(1)/ds(1, k) + srcig_local(1, k) else A_ig(1)=0.0 B_ig(1)=1.0/dt @@ -647,7 +584,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec 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 - R_ig(ntheta) = oneoverdt*ee_ig(ntheta,k) + cgprev_ig(ntheta)*eeprev_ig(ntheta)/ds(ntheta,k) + srcsh_local(ntheta, k) + R_ig(ntheta) = oneoverdt*ee_ig(ntheta,k) + cgprev_ig(ntheta)*eeprev_ig(ntheta)/ds(ntheta,k) + srcig_local(ntheta, k) else A_ig(ntheta) = 0.0 B_ig(ntheta) = oneoverdt @@ -764,7 +701,7 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec endif ! enddo - ! + ! do k=1,no_nodes ! ! Compute directionally integrated parameters for output @@ -799,22 +736,22 @@ subroutine solve_energy_balance2Dstat(x,y,no_nodes,w,ds,inner,prev,neumannconnec call baldock(g, rho, alfa_ig, gamma_ig, kwav_ig(k), depth(k), H_ig(k), T_ig, baldock_opt, Dw_ig(k), Hmx_ig(k)) endif ! - ! average beta, alphaig, srcsh over directions + ! average beta, alphaig, srcig over directions betamean(k) = sum(beta_local(:,k))/ntheta ! real mean !betamean(k) = maxval(beta_local(:,k)) ! alphaig(k) = sum(alphaig_local(:,k))/ntheta ! real mean ! - srcsh(k) = sum(srcsh_local(:,k)) /ntheta ! real mean - !srcsh(k) = maxval(srcsh_local(:,k)) - !srcsh(k) = sum(srcsh_local(:,k))*dtheta + srcig(k) = sum(srcig_local(:,k)) /ntheta ! real mean + !srcig(k) = maxval(srcig_local(:,k)) + !srcig(k) = sum(srcig_local(:,k))*dtheta ! endif ! endif endif ! - enddo + enddo ! callno=callno+1 ! @@ -911,6 +848,136 @@ subroutine baldock (rho,g,alfa,gamma,k,depth,H,T,opt,Dw,Hmax) ! end subroutine baldock + subroutine determine_infragravity_source_sink_term(no_nodes, ntheta, k, w, ds, prev, cg_ig, nwav, depth, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, zb, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local) + ! + ! Incoming variables + integer, intent(in) :: no_nodes,ntheta ! number of grid points, number of directions + integer, intent(in) :: k ! counters (k is grid index) + real*4, dimension(2,ntheta,no_nodes),intent(in) :: w ! weights of upwind grid points, 2 per grid point and per wave direction + real*4, dimension(ntheta,no_nodes), intent(in) :: ds ! distance to interpolated upwind point, per grid point and direction + integer, dimension(2,ntheta,no_nodes),intent(in) :: prev ! two upwind grid points per grid point and wave direction + real*4, dimension(no_nodes), intent(in) :: cg_ig ! group velocity + real*4, dimension(no_nodes), intent(in) :: nwav ! wave number n + real*4, dimension(no_nodes), intent(in) :: depth ! water depth + real*4, dimension(no_nodes), intent(in) :: H ! wave height + real*4, dimension(ntheta,no_nodes), intent(in) :: ee ! energy density + real*4, dimension(ntheta,no_nodes), intent(in) :: ee_ig ! energy density infragravity waves + real*4, dimension(no_nodes), intent(in) :: zb ! actual bed level + 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 + ! + ! 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(:,:), intent(inout) :: beta_local ! Local bed slope based on bed level per direction + ! + ! Internal variables + integer :: itheta ! directional counter + integer :: k1,k2 ! upwind counters (k is grid index) + + real*4, dimension(ntheta,no_nodes) :: Sxx ! Radiation Stress + real*4, dimension(:), allocatable :: Sxxprev ! radiation stress at upwind intersection point + real*4, dimension(:), allocatable :: depthprev ! water depth at upwind intersection point + real*4, dimension(:), allocatable :: Hprev ! Incident wave height at upwind intersection point + real*4 :: beta ! local bedslope + real*4 :: gam ! local gamma (Hinc / depth ratio) + ! + ! Allocate internal variables + allocate(Sxxprev(ntheta)) + allocate(depthprev(ntheta)) + allocate(Hprev(ntheta)) + ! + ! + ! Compute exchange source term inc to ig waves - per direction + do itheta = 1, ntheta + ! + k1 = prev(1, itheta, k) + k2 = prev(2, itheta, k) + ! + if (k1>0 .and. k2>0) then ! IMPORTANT - for some reason (k1*k2)>0 is not reliable always, resulting in directions being uncorrectly skipped!!! + ! + ! First calculate upwind direction dependent variables + ! 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 + ! + Sxxprev(itheta) = w(1, itheta, k)*Sxx(itheta,k1) + w(2, itheta, k)*Sxx(itheta,k2) + ! + depthprev(itheta) = w(1, itheta, k)*depth(k1) + w(2, itheta, k)*depth(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) + ! + Hprev(itheta) = w(1, itheta, k)*H(k1) + w(2, itheta, k)*H(k2) + ! + ! Determine local bed slope and relative waterdepth 'gam' + beta = max((w(1, itheta, k)*(zb(k) - zb(k1)) + w(2, itheta, k)*(zb(k) - zb(k2)))/ds(itheta, k), 0.0) + ! Notes: + ! - use actual bed level now for slope, because depth changes because of wave setup/tide/surge + ! - in zb, depth is negative > therefore zb(k) minus zb(k1) + ! - beta=0 means a horizontal or decreasing slope > need alphaig=0 then + ! + beta_local(itheta,k) = beta + !betan_local(itheta,k) = (beta/sigm_ig)*sqrt(9.81/max(depth(k), hmin)) ! TL: in case in the future we would need the normalised bed slope again + ! + gam = max(0.5*(Hprev(itheta)/depthprev(itheta) + H(k)/depth(k)), 0.0) ! mean gamma over current and upwind point + ! + ! Determine dSxx and IG source/sink term 'srcig' + ! + if (ig_opt == 1 .or. ig_opt == 2) then + ! + ! Calculate shoaling parameter alpha_ig following Leijnse et al. (2024) + ! + call estimate_shoaling_parameter_alphaig(beta, gam, alphaig_local(itheta,k)) ! [input, input, output] + ! + ! Now calculate source term component + ! + ! Newest dSxx/dx based method, using estimate of Sxx(k) using conservative shoaling + if (Sxxprev(itheta)<=0.0) then + ! + srcig_local(itheta, k) = 0.0 !Avoid big jumps in dSxx that can happen if a upwind point is a boundary point with Hinc=0 + ! + else + ! + if (ig_opt == 1) 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) * dSxx / ds(itheta, k) + ! + endif + ! + else ! TL: option to add future parameterisations here for e.g. coral reef type coasts + ! + srcig_local(itheta, k) = 0.0 + ! + endif + ! + srcig_local(itheta, k) = max(srcig_local(itheta, k), 0.0) + ! + endif + ! + enddo + ! + end subroutine determine_infragravity_source_sink_term + subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) real*4, intent(in) :: beta real*4, intent(in) :: gam @@ -951,7 +1018,7 @@ subroutine estimate_shoaling_parameter_alphaig(beta, gam, alphaig) ! endif ! - ! Limit alphaig betwwen [0, 1] to prevent large overshoots in case of low gamma and very small beta + ! Limit alphaig between [0, 1] to prevent large overshoots in case of low gamma and very small beta ! alphaig = max(alphaig, 0.0) alphaig = min(alphaig, 1.0) @@ -974,7 +1041,7 @@ subroutine disper_approx(h,T,k,n,C,Cg,no_nodes) Cg = n*C ! end subroutine disper_approx - + ! ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino