From beb9f1ac24557e9b5da9a11d3b7068019c1fa02f Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 13:28:25 +0100 Subject: [PATCH 01/19] - Initial commit, working on make a generic read_netcdf_quadtree_generic function - Apply first to netcdf storage volume - Work in progress still --- source/src/sfincs_lib.f90 | 4 +- source/src/sfincs_ncinput.F90 | 107 ++++++++++++++++++++++++++++------ 2 files changed, 91 insertions(+), 20 deletions(-) diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index a08899641..8405bef7b 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -92,8 +92,8 @@ function sfincs_initialize() result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = "$Rev: v2.3.0 mt. Faber Release" - build_date = "$Date: 2025-11-18" + build_revision = "$Rev: v2.3.0 mt. Faber Release-branch:268" + build_date = "$Date: 2025-02-04" ! call write_log('', 1) call write_log('------------ Welcome to SFINCS ------------', 1) diff --git a/source/src/sfincs_ncinput.F90 b/source/src/sfincs_ncinput.F90 index a29e49356..dc9b9a8c0 100644 --- a/source/src/sfincs_ncinput.F90 +++ b/source/src/sfincs_ncinput.F90 @@ -59,6 +59,11 @@ module sfincs_ncinput integer :: np_dimid integer :: vol_varid end type + type net_type_generic + integer :: ncid + integer :: np_dimid + integer :: gen_varid + end type ! type(net_type_bndbzsbzi) :: net_file_bndbzsbzi type(net_type_srcdis) :: net_file_srcdis @@ -66,7 +71,8 @@ module sfincs_ncinput type(net_type_amp) :: net_file_amp type(net_type_ampr) :: net_file_ampr type(net_type_spw) :: net_file_spw - type(net_type_vol) :: net_file_vol + type(net_type_vol) :: net_file_vol + type(net_type_generic) :: net_file_generic contains @@ -227,52 +233,117 @@ subroutine read_netcdf_discharge_data() subroutine read_netcdf_storage_volume() ! - use netcdf + !use netcdf use sfincs_data use quadtree ! implicit none ! - real*8, dimension(:), allocatable :: vols - integer :: nrcells, nm, ip + !real*8, dimension(:), allocatable :: vols + integer :: precision ! character (len=256), parameter :: vol_varname = 'vol' ! - NF90(nf90_open(trim(volfile), NF90_CLOBBER, net_file_vol%ncid)) + !NF90(nf90_open(trim(volfile), NF90_CLOBBER, net_file_vol%ncid)) + !! + !! Get dimensions id's: nr points + !! + !NF90(nf90_inq_dimid(net_file_vol%ncid, "mesh2d_nFaces", net_file_vol%np_dimid)) + !! + !! Get dimensions sizes + !! + !NF90(nf90_inquire_dimension(net_file_vol%ncid, net_file_vol%np_dimid, len = nrcells)) ! nr of cells + !! + !! Check that number of values in the cell matches quadtree_nr_points + !! + !! TODO: if (nrcells /=quadtree_nr_points) GIVE ERROR and stop simulation + !! + !NF90(nf90_inq_varid(net_file_vol%ncid, vol_varname, net_file_vol%vol_varid)) + !! + !allocate(vols(nrcells)) + !! + !NF90(nf90_get_var(net_file_vol%ncid, net_file_vol%vol_varid, vols(:))) + ! + ! Call the generic function + precision = 4 + ! + call read_netcdf_quadtree_generic(trim(volfile), vol_varname, storage_volume, np, precision) !ncfile, varname, varout, size, precision) + + !! Map quadtree to sfincs --> Question: better in or out of function? + !! + !do ip = 1, quadtree_nr_points + ! ! + ! nm = index_sfincs_in_quadtree(ip) + ! ! + ! storage_volume(nm) = vols(ip) + ! ! + !enddo + ! + !NF90(nf90_close(net_file_vol%ncid)) + ! + end subroutine + + subroutine read_netcdf_quadtree_generic(ncfile, varname, var, npoints, precision) + ! For instance: storage_volume.nc, vol, storage_volume, np, 5 + ! + use netcdf + use sfincs_data + use quadtree + ! + implicit none + ! + integer :: nm, npoints, ip, nrcells, precision + ! + character*256 :: ncfile + character*256 :: varname + !character (len=256), parameter :: varname + ! + !real*precision, dimension(:), allocatable :: varout + real*4, dimension(npoints), intent(inout) :: var ! variable that we are mapping to + + if (precision .eqv. 4) then + real*4, dimension(:), allocatable :: vartmp + elseif (precision .eqv. 8) then + real*8, dimension(:), allocatable :: vartmp + else + call write_log('ERROR : precision should be 4 or 8 ...', 0) + endif + ! + ! Open netcdf file + NF90(nf90_open(trim(ncfile), NF90_CLOBBER, net_file_generic%ncid)) ! ! Get dimensions id's: nr points ! - NF90(nf90_inq_dimid(net_file_vol%ncid, "mesh2d_nFaces", net_file_vol%np_dimid)) + NF90(nf90_inq_dimid(net_file_generic%ncid, "mesh2d_nFaces", net_file_generic%np_dimid)) ! ! Get dimensions sizes ! - NF90(nf90_inquire_dimension(net_file_vol%ncid, net_file_vol%np_dimid, len = nrcells)) ! nr of cells + NF90(nf90_inquire_dimension(net_file_generic%ncid, net_file_generic%np_dimid, len = nrcells)) ! nr of cells ! ! Check that number of values in the cell matches quadtree_nr_points - ! + write(logstr,*)'Info - number of cells in file: ',nrcells, ' vs quadtree_nr_points: ',quadtree_nr_points + call write_log(logstr, 1) ! TODO: if (nrcells /=quadtree_nr_points) GIVE ERROR and stop simulation ! - NF90(nf90_inq_varid(net_file_vol%ncid, vol_varname, net_file_vol%vol_varid)) + NF90(nf90_inq_varid(net_file_generic%ncid, varname, net_file_generic%gen_varid)) ! - allocate(vols(nrcells)) + allocate(vartmp(nrcells)) ! - NF90(nf90_get_var(net_file_vol%ncid, net_file_vol%vol_varid, vols(:))) + NF90(nf90_get_var(net_file_generic%ncid, net_file_generic%gen_varid, vartmp(:))) ! - ! Map quadtree to sfincs + ! Map quadtree to sfincs > externally for now ! do ip = 1, quadtree_nr_points ! nm = index_sfincs_in_quadtree(ip) ! - storage_volume(nm) = vols(ip) + var(nm) = vartmp(ip) ! enddo - ! - NF90(nf90_close(net_file_vol%ncid)) + ! + NF90(nf90_close(net_file_generic%ncid)) ! - end subroutine - - + end subroutine subroutine read_netcdf_amuv_data() ! Output is made exactly the same as original read_amuv_dimensions & read_amuv_file subroutines but then with data given by netcdf file From 68e1493dd88d527407a82074b0402b3c85d747be Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 15:11:22 +0100 Subject: [PATCH 02/19] - Working version for storage volume nc with fixed mapping to SFINCS (can always add a snapwave one) and with real*4 (can always make a double version as well) - In sfincs_domain.f90, 'call read_netcdf_storage_volume' is now replaced by 'call read_netcdf_quadtree_to_sfincs' - Succesfully tested for testbed run 'storage_volume_qt_sbg_thd2' --- source/src/sfincs_domain.f90 | 7 ++-- source/src/sfincs_lib.f90 | 2 +- source/src/sfincs_ncinput.F90 | 77 ++++++++--------------------------- 3 files changed, 22 insertions(+), 64 deletions(-) diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index a1b532f2d..8b789fb50 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2403,6 +2403,7 @@ subroutine initialize_storage_volume() ! integer :: nchar logical :: ok + character*256 :: varname ! if (use_storage_volume) then ! @@ -2424,9 +2425,9 @@ subroutine initialize_storage_volume() ! if (volfile(nchar - 1 : nchar) == 'nc') then ! - ! Read netcdf file - ! - call read_netcdf_storage_volume() + ! Call the generic quadtree nc file reader function + varname = 'vol' + call read_netcdf_quadtree_to_sfincs(volfile, varname, storage_volume) !ncfile, varname, varout) ! else ! diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index 8405bef7b..898905cc2 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -92,7 +92,7 @@ function sfincs_initialize() result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = "$Rev: v2.3.0 mt. Faber Release-branch:268" + build_revision = "$Rev: v2.3.1 mt. Faber" build_date = "$Date: 2025-02-04" ! call write_log('', 1) diff --git a/source/src/sfincs_ncinput.F90 b/source/src/sfincs_ncinput.F90 index dc9b9a8c0..dfe2c5198 100644 --- a/source/src/sfincs_ncinput.F90 +++ b/source/src/sfincs_ncinput.F90 @@ -239,52 +239,12 @@ subroutine read_netcdf_storage_volume() ! implicit none ! - !real*8, dimension(:), allocatable :: vols - integer :: precision - ! - character (len=256), parameter :: vol_varname = 'vol' - ! - !NF90(nf90_open(trim(volfile), NF90_CLOBBER, net_file_vol%ncid)) - !! - !! Get dimensions id's: nr points - !! - !NF90(nf90_inq_dimid(net_file_vol%ncid, "mesh2d_nFaces", net_file_vol%np_dimid)) - !! - !! Get dimensions sizes - !! - !NF90(nf90_inquire_dimension(net_file_vol%ncid, net_file_vol%np_dimid, len = nrcells)) ! nr of cells - !! - !! Check that number of values in the cell matches quadtree_nr_points - !! - !! TODO: if (nrcells /=quadtree_nr_points) GIVE ERROR and stop simulation - !! - !NF90(nf90_inq_varid(net_file_vol%ncid, vol_varname, net_file_vol%vol_varid)) - !! - !allocate(vols(nrcells)) - !! - !NF90(nf90_get_var(net_file_vol%ncid, net_file_vol%vol_varid, vols(:))) - ! - ! Call the generic function - precision = 4 - ! - call read_netcdf_quadtree_generic(trim(volfile), vol_varname, storage_volume, np, precision) !ncfile, varname, varout, size, precision) - - !! Map quadtree to sfincs --> Question: better in or out of function? - !! - !do ip = 1, quadtree_nr_points - ! ! - ! nm = index_sfincs_in_quadtree(ip) - ! ! - ! storage_volume(nm) = vols(ip) - ! ! - !enddo - ! - !NF90(nf90_close(net_file_vol%ncid)) + ! end subroutine - subroutine read_netcdf_quadtree_generic(ncfile, varname, var, npoints, precision) - ! For instance: storage_volume.nc, vol, storage_volume, np, 5 + subroutine read_netcdf_quadtree_to_sfincs(ncfile, varname, var) + ! For instance: storage_volume.nc, vol, storage_volume ! use netcdf use sfincs_data @@ -292,24 +252,17 @@ subroutine read_netcdf_quadtree_generic(ncfile, varname, var, npoints, precision ! implicit none ! - integer :: nm, npoints, ip, nrcells, precision + integer :: nm, ip, nrcells, precision ! character*256 :: ncfile character*256 :: varname - !character (len=256), parameter :: varname ! - !real*precision, dimension(:), allocatable :: varout - real*4, dimension(npoints), intent(inout) :: var ! variable that we are mapping to - - if (precision .eqv. 4) then - real*4, dimension(:), allocatable :: vartmp - elseif (precision .eqv. 8) then - real*8, dimension(:), allocatable :: vartmp - else - call write_log('ERROR : precision should be 4 or 8 ...', 0) - endif + real*4, dimension(np), intent(inout) :: var ! variable that we are mapping to + ! + real*4, dimension(:), allocatable :: vartmp ! ! Open netcdf file + ! NF90(nf90_open(trim(ncfile), NF90_CLOBBER, net_file_generic%ncid)) ! ! Get dimensions id's: nr points @@ -320,10 +273,14 @@ subroutine read_netcdf_quadtree_generic(ncfile, varname, var, npoints, precision ! NF90(nf90_inquire_dimension(net_file_generic%ncid, net_file_generic%np_dimid, len = nrcells)) ! nr of cells ! - ! Check that number of values in the cell matches quadtree_nr_points - write(logstr,*)'Info - number of cells in file: ',nrcells, ' vs quadtree_nr_points: ',quadtree_nr_points - call write_log(logstr, 1) - ! TODO: if (nrcells /=quadtree_nr_points) GIVE ERROR and stop simulation + ! Check that number of values in the cell matches quadtree_nr_points + ! (=all quadtree cells, not just the active ones) + ! + if (nrcells /=quadtree_nr_points) then + write(logstr,*)'Info - file: ',trim(ncfile),' contains ',nrcells, 'input points vs quadtree_nr_points: ',quadtree_nr_points,' in sfincs.nc quadtree grid' + call write_log(logstr, 0) + call stop_sfincs('Error ! Number of grid points in netcdf input file does not match with that in sfincs.nc quadtree grid!', 1) + endif ! NF90(nf90_inq_varid(net_file_generic%ncid, varname, net_file_generic%gen_varid)) ! @@ -331,7 +288,7 @@ subroutine read_netcdf_quadtree_generic(ncfile, varname, var, npoints, precision ! NF90(nf90_get_var(net_file_generic%ncid, net_file_generic%gen_varid, vartmp(:))) ! - ! Map quadtree to sfincs > externally for now + ! Map quadtree to sfincs variable ! do ip = 1, quadtree_nr_points ! From 245c1a14b7e3147e39967db1997c270139ff19fa Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 15:55:10 +0100 Subject: [PATCH 03/19] - Add error when wanted variable was not found in file and directly stop - Tested by supplying dummy stoarge volume nc file for manning, now nicely stops --- source/src/sfincs_ncinput.F90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/source/src/sfincs_ncinput.F90 b/source/src/sfincs_ncinput.F90 index dfe2c5198..7f4ef894c 100644 --- a/source/src/sfincs_ncinput.F90 +++ b/source/src/sfincs_ncinput.F90 @@ -252,7 +252,7 @@ subroutine read_netcdf_quadtree_to_sfincs(ncfile, varname, var) ! implicit none ! - integer :: nm, ip, nrcells, precision + integer :: nm, ip, nrcells, status ! character*256 :: ncfile character*256 :: varname @@ -277,12 +277,17 @@ subroutine read_netcdf_quadtree_to_sfincs(ncfile, varname, var) ! (=all quadtree cells, not just the active ones) ! if (nrcells /=quadtree_nr_points) then - write(logstr,*)'Info - file: ',trim(ncfile),' contains ',nrcells, 'input points vs quadtree_nr_points: ',quadtree_nr_points,' in sfincs.nc quadtree grid' - call write_log(logstr, 0) - call stop_sfincs('Error ! Number of grid points in netcdf input file does not match with that in sfincs.nc quadtree grid!', 1) + write(logstr,*)'Error : netcdf input file ',trim(ncfile),' contains: ',nrcells, 'input points, while expected is: ',quadtree_nr_points,' as in sfincs.nc quadtree grid' + call stop_sfincs(trim(logstr), 1) endif ! - NF90(nf90_inq_varid(net_file_generic%ncid, varname, net_file_generic%gen_varid)) + status = nf90_inq_varid(net_file_generic%ncid, varname, net_file_generic%gen_varid) + ! + ! Stop SFINCS if wanted variable was not found + if (status /= nf90_noerr) then + write(logstr,'(a,a,a,a,a)')'Error : netcdf input file ',trim(ncfile),' does not contain needed variable: ',trim(varname),' !' + call stop_sfincs(trim(logstr), 1) + endif ! allocate(vartmp(nrcells)) ! From d3db40e74bf1ab2d8cf37a4abdc4ea7306fe15ee Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 17:07:57 +0100 Subject: [PATCH 04/19] - Add netinfiltrationfile, netinftype variables - If netinfiltrationfile exists, then netinftype should be used - Do precheck of netinfiltrationfile before reading in all data --- source/src/sfincs_data.f90 | 3 +++ source/src/sfincs_domain.f90 | 50 +++++++++++++++++++++++++++++++++--- source/src/sfincs_input.f90 | 5 +++- source/src/sfincs_lib.f90 | 2 +- 4 files changed, 54 insertions(+), 6 deletions(-) diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index c778c4925..4f3b74156 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -160,6 +160,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 @@ -186,6 +187,7 @@ module sfincs_data character*3 :: outputtype_his character*3 :: utmzone character*3 :: inftype + character*3 :: netinftype integer :: epsg character*15 :: epsg_code integer :: nc_deflate_level @@ -225,6 +227,7 @@ module sfincs_data logical :: write_time_output logical :: bziwaves logical :: infiltration + LOGICAL :: netcdf_infiltration logical :: debug logical :: radstr logical :: crsgeo diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 8b789fb50..c5a7a9db5 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2034,6 +2034,7 @@ subroutine initialize_infiltration() ! Note, infiltration methods not designed to be stacked ! infiltration = .false. + netcdf_infiltration = .false. ! ! Four options for infiltration: ! @@ -2045,6 +2046,10 @@ subroutine initialize_infiltration() ! 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) + ! 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 ! @@ -2053,11 +2058,27 @@ subroutine initialize_infiltration() ! or: ! b) inftype == 'cna' or inftype == 'cnb' ! - ! First we determine precipitation type + !!!!!!!!!!!!!!!!!!!!! + ! Initializing steps: + !!!!!!!!!!!!!!!!!!!!! + ! + ! 1) First we determine infiltration type ! if (precip) then ! - if (qinf > 0.0) then + if (netinfiltrationfile /= 'none') then + ! + ! inftype is user defined, keyword: 'netinftype' in sfincs.inp: + ! + inftype = netinftype + ! + ! inftype is either: c2d, cna, cnb, gai, hor + ! 'inftype = con' is not relevant for netcdf input + ! + infiltration = .true. + netcdf_infiltration = .true. + ! + elseif (qinf > 0.0) then ! ! Spatially-uniform constant infiltration (specified as +mm/hr) ! @@ -2102,7 +2123,7 @@ subroutine initialize_infiltration() ! endif ! - ! We need cumprcp and cuminf + ! 2) We need cumprcp and cuminf ! allocate(cumprcp(np)) cumprcp = 0.0 @@ -2110,7 +2131,7 @@ subroutine initialize_infiltration() allocate(cuminf(np)) cuminf = 0.0 ! - ! Now allocate and read spatially-varying inputs + ! 3) Now allocate and read spatially-varying inputs ! if (infiltration) then ! @@ -2119,10 +2140,31 @@ subroutine initialize_infiltration() ! endif ! + ! 4) Pre-check whether netcdf infiltration file exists - once + ! + if (netcdf_infiltration) then + ! + write(logstr,'(a)')'Info : turning on infiltration from netcdf input file' + call write_log(logstr, 0) + ! + write(logstr,'(a,a)')'Info : reading netcdf infiltration file ', trim(netinfiltrationfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(netinfiltrationfile, 'Infiltration netcdf file', .true.) + ! + write(logstr,'(a,a)')'Info : specified inftype is ', trim(inftype) + call write_log(logstr, 0) + ! + endif + ! + ! 5) Read in data per type, either from ascii or general netcdf file + ! if (inftype == 'con') then ! ! Spatially-uniform constant infiltration (specified as +mm/hr) ! + ! Note : Input directly in sfincs.inp, so no file needs to be read + ! write(logstr,'(a)')'Info : turning on spatially-uniform constant infiltration' call write_log(logstr, 0) ! diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index aba244eeb..a5b93ddb7 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -207,7 +207,10 @@ subroutine read_sfincs_input() call read_char_input(500,'netamuamvfile',netamuamvfile,'none') 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,'netspwfile',netspwfile,'none') + ! + call read_char_input(500,'netinfiltrationfile',netinfiltrationfile,'none') + call read_char_input(500,'netinftype',netinftype,'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 898905cc2..b9f303f3a 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -92,7 +92,7 @@ function sfincs_initialize() result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = "$Rev: v2.3.1 mt. Faber" + build_revision = "$Rev: v2.3.1 mt. Faber+" build_date = "$Date: 2025-02-04" ! call write_log('', 1) From 435327424067c62af8f5887a46aa83315f66b351 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 17:48:48 +0100 Subject: [PATCH 05/19] - Implement netcdf reader option for all infiltration options - Current implementation will check and retrieve multiple variables from the same file (e.g. cnb) individually - Added some documentation what part of the code is reading/allocating/generic needed conversions - Kept naming in netcdf files for now consistent with the equivalent netcdf files --- source/src/sfincs_domain.f90 | 314 +++++++++++++++++++++++++---------- 1 file changed, 227 insertions(+), 87 deletions(-) diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index c5a7a9db5..c7c0d5547 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2020,6 +2020,7 @@ subroutine initialize_roughness() subroutine initialize_infiltration() ! use sfincs_data + use sfincs_ncinput ! implicit none ! @@ -2027,6 +2028,8 @@ subroutine initialize_infiltration() ! logical :: ok ! + character*256 :: varname + ! ! INFILTRATION ! ! Infiltration only works when rainfall is activated ! If you want infiltration without rainfall, use a precip file with 0.0s @@ -2195,17 +2198,32 @@ subroutine initialize_infiltration() write(logstr,'(a)')'Info : turning on spatially-varying constant infiltration' call write_log(logstr, 0) ! - ! Read spatially-varying infiltration (only binary, specified in +mm/hr) + allocate(qinffield(np)) ! - write(logstr,'(a,a)')'Info : reading infiltration file ', trim(qinffile) - call write_log(logstr, 0) + qinffield = 0.0 ! - ok = check_file_exists(qinffile, 'Infiltration qinf file', .true.) + ! Read spatially-varying infiltration (specified in +mm/hr) ! - allocate(qinffield(np)) - open(unit = 500, file = trim(qinffile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'qinf' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, qinffield) !ncfile, varname, varout) + ! + else ! from separate qinffile - only binary: + ! + write(logstr,'(a,a)')'Info : reading infiltration file ', trim(qinffile) + call write_log(logstr, 0) + ! + ok = check_file_exists(qinffile, 'Infiltration qinf file', .true.) + ! + open(unit = 500, file = trim(qinffile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + endif + ! + ! Generic needed conversion: ! qinffield = qinffield / 3600 / 1000 ! convert to +m/s ! @@ -2220,18 +2238,29 @@ subroutine initialize_infiltration() ! qinffield = 0.0 ! - write(logstr,'(a,a)')'Info : reading scs file ',trim(scsfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(scsfile, 'Infiltration scs file', .true.) - ! - open(unit = 500, file = trim(scsfile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'scs' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, qinffield) + ! + else ! from separate scsfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading scs file ',trim(scsfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(scsfile, 'Infiltration scs file', .true.) + ! + open(unit = 500, file = trim(scsfile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + endif ! - ! already convert qinffield from inches to m here + ! Generic needed conversion: ! qinffield = qinffield * 0.0254 ! to m + ! already convert qinffield from inches to m here ! elseif (inftype == 'cnb') then ! @@ -2243,48 +2272,83 @@ subroutine initialize_infiltration() ! Allocate Smax allocate(qinffield(np)) qinffield = 0.0 - write(logstr,'(a,a)')'Info : reading smax file ',trim(smaxfile) - call write_log(logstr, 0) ! - ok = check_file_exists(smaxfile, 'Infiltration smax file', .true.) - ! - open(unit = 500, file = trim(smaxfile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'smax' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, qinffield) + ! + else ! from separate smaxfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading smax file ',trim(smaxfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(smaxfile, 'Infiltration smax file', .true.) + ! + open(unit = 500, file = trim(smaxfile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + endif ! ! Allocate Se allocate(scs_Se(np)) scs_Se = 0.0 - write(logstr,'(a,a)')'Info : reading seff file ',trim(sefffile) - call write_log(logstr, 0) ! - ok = check_file_exists(sefffile, 'Infiltration seff file', .true.) - ! - open(unit = 501, file = trim(sefffile), form = 'unformatted', access = 'stream') - read(501)scs_Se - close(501) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'seff' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, scs_Se) + ! + else ! from separate sefffile - only binary: + ! + write(logstr,'(a,a)')'Info : reading seff file ',trim(sefffile) + call write_log(logstr, 0) + ! + ok = check_file_exists(sefffile, 'Infiltration seff file', .true.) + ! + open(unit = 501, file = trim(sefffile), form = 'unformatted', access = 'stream') + read(501)scs_Se + close(501) + ! + endif ! - ! Compute recovery ! Equation 4-36 ! Allocate Ks ! allocate(ksfield(np)) ksfield = 0.0 - write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) - call write_log(logstr, 0) ! - ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) - ! - open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') - read(502)ksfield - close(502) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'ks' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, ksfield) + ! + else ! from separate ksfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) + ! + open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') + read(502)ksfield + close(502) + ! + endif ! + ! Generic needed conversion: + ! ! 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 support variables: + ! allocate(scs_P1(np)) scs_P1 = 0.0 allocate(scs_F1(np)) @@ -2306,41 +2370,77 @@ subroutine initialize_infiltration() ! allocate(GA_head(np)) GA_head = 0.0 - write(logstr,'(a,a)')'Info : reading psi file ',trim(psifile) - call write_log(logstr, 0) ! - ok = check_file_exists(psifile, 'Infiltration psi file', .true.) - ! - open(unit = 500, file = trim(psifile), form = 'unformatted', access = 'stream') - read(500)GA_head - close(500) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'psi' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, GA_head) + ! + else ! from separate psifile - only binary: + ! + write(logstr,'(a,a)')'Info : reading psi file ',trim(psifile) + call write_log(logstr, 0) + ! + ok = check_file_exists(psifile, 'Infiltration psi file', .true.) + ! + open(unit = 500, file = trim(psifile), form = 'unformatted', access = 'stream') + read(500)GA_head + close(500) + ! + endif ! ! Allocate maximum soil moisture deficit ! allocate(GA_sigma_max(np)) GA_sigma_max = 0.0 - write(logstr,'(a,a)')'Info : reading sigma file ',trim(sigmafile) - call write_log(logstr, 0) ! - ok = check_file_exists(sigmafile, 'Infiltration sigma file', .true.) - ! - open(unit = 501, file = trim(sigmafile), form = 'unformatted', access = 'stream') - read(501)GA_sigma_max - close(501) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'sigma' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, GA_sigma_max) + ! + else ! from separate sigmafile - only binary: + ! + write(logstr,'(a,a)')'Info : reading sigma file ',trim(sigmafile) + call write_log(logstr, 0) + ! + ok = check_file_exists(sigmafile, 'Infiltration sigma file', .true.) + ! + open(unit = 501, file = trim(sigmafile), form = 'unformatted', access = 'stream') + read(501)GA_sigma_max + close(501) + ! + endif ! ! Allocate saturated hydraulic conductivity ! allocate(ksfield(np)) ksfield = 0.0 - write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) ! - open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') - read(502)ksfield - close(502) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'ks' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, ksfield) + ! + else ! from separate ksfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) + ! + open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') + read(502)ksfield + close(502) + ! + endif + ! + ! Generic needed conversion: + ! ! Compute recovery ! Equation 4-36 ! allocate(inf_kr(np)) @@ -2367,6 +2467,8 @@ subroutine initialize_infiltration() ! ! First time step doesnt have an estimate yet ! + ! Allocate support variables: + ! allocate(qinffield(np)) qinffield(nm) = 0.0 ! @@ -2377,48 +2479,86 @@ subroutine initialize_infiltration() call write_log('Info : turning on process infiltration (via modified Horton)', 0) ! ! Horton: final infiltration capacity (fc) - ! Note that qinffield = horton_fc + ! Note that qinffield = horton_fc (/3600/1000, see below) + ! allocate(horton_fc(np)) horton_fc = 0.0 - write(logstr,'(a,a)')'Info : reading fc file ',trim(fcfile) - call write_log(logstr, 0) ! - ok = check_file_exists(fcfile, 'Infiltration fc file', .true.) - ! - open(unit = 500, file = trim(fcfile), form = 'unformatted', access = 'stream') - read(500)horton_fc - close(500) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'fc' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, horton_fc) + ! + else ! from separate fcfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading fc file ',trim(fcfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(fcfile, 'Infiltration fc file', .true.) + ! + open(unit = 500, file = trim(fcfile), form = 'unformatted', access = 'stream') + read(500)horton_fc + close(500) + ! + endif ! ! Horton: initial infiltration capacity (f0) allocate(horton_f0(np)) horton_f0 = 0.0 - write(logstr,'(a,a)')'Info : reading f0 file ',trim(f0file) - call write_log(logstr, 0) - ! - ok = check_file_exists(f0file, 'Infiltration f0 file', .true.) - ! - 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 + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'f0' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, horton_f0) + ! + else ! from separate f0file - only binary: + ! + write(logstr,'(a,a)')'Info : reading f0 file ',trim(f0file) + call write_log(logstr, 0) + ! + ok = check_file_exists(f0file, 'Infiltration f0 file', .true.) + ! + open(unit = 501, file = trim(f0file), form = 'unformatted', access = 'stream') + read(501)horton_f0 + close(501) + ! + endif ! ! 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(logstr,'(a,a)')'Info : reading kd file ',trim(kdfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(kdfile, 'Infiltration kd file', .true.) ! - open(unit = 502, file = trim(kdfile), form = 'unformatted', access = 'stream') - read(502)horton_kd - close(502) + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'kd' + call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, horton_kd) + ! + else ! from separate kdfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading kd file ',trim(kdfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(kdfile, 'Infiltration kd file', .true.) + ! + open(unit = 502, file = trim(kdfile), form = 'unformatted', access = 'stream') + read(502)horton_kd + close(502) + ! + endif ! - write(logstr,'(a,a)')'Using constant recovery rate that is based on constant factor relative to ',trim(kdfile) + write(logstr,'(a,a)')'Info : Using constant recovery rate that is based on constant factor relative to ',trim(kdfile) call write_log(logstr, 0) ! + ! Generic needed conversion: + ! + ! Prescribe the current estimate (for output only; initial capacity) + qinffield = horton_f0/3600/1000 + ! + ! Allocate support variables: + ! ! Estimate of time allocate(rain_T1(np)) rain_T1 = 0.0 From 378e27e65b86d71b8108eadeda95a8a803397280 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 17:56:14 +0100 Subject: [PATCH 06/19] - Add stop_sfincs check on combination check of infiltration input type (orignal vs netcdf) vs grid type (regular vs quadtree) - Idea for now is if quadtree model, only netcdf input is support, and for regular model only the original binary files --- source/src/sfincs_domain.f90 | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index c7c0d5547..2b5bf86a5 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2159,8 +2159,28 @@ subroutine initialize_infiltration() call write_log(logstr, 0) ! endif - ! - ! 5) Read in data per type, either from ascii or general netcdf file + ! + ! 5) Check whether infiltration input type (orignal vs netcdf) are correctly matched to grid type (regular vs quadtree) + ! + if (netcdf_infiltration) then + ! + if (use_quadtree .eqv. .false.) then + ! + call stop_sfincs('Error ! Netcdf infiltration input format can only be specified for quadtree mesh model !', 1) + ! + endif + ! + else ! Original + ! + if (use_quadtree .eqv. .true.) then + ! + call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the netinfiltrationfile Netcdf format! !', 1) + ! + endif + ! + endif + ! + ! 6) Read in data per type, either from ascii or general netcdf file ! if (inftype == 'con') then ! From a2f67b4a2997b3f11a8b3f5395fc5846267cb469 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 18:05:49 +0100 Subject: [PATCH 07/19] - Constant uniform infiltration is allowed for both grid types --- source/src/sfincs_domain.f90 | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 2b5bf86a5..4ecb3258d 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2162,23 +2162,26 @@ subroutine initialize_infiltration() ! ! 5) Check whether infiltration input type (orignal vs netcdf) are correctly matched to grid type (regular vs quadtree) ! - if (netcdf_infiltration) then - ! - if (use_quadtree .eqv. .false.) then - ! - call stop_sfincs('Error ! Netcdf infiltration input format can only be specified for quadtree mesh model !', 1) + if (inftype /= 'con') then !constant uniform works for both options + ! + if (netcdf_infiltration) then + ! + if (use_quadtree .eqv. .false.) then + ! + call stop_sfincs('Error ! Netcdf infiltration input format can only be specified for quadtree mesh model !', 1) + ! + endif ! + else ! Original + ! + if (use_quadtree .eqv. .true.) then + ! + call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the netinfiltrationfile Netcdf ormat! !', 1) + endif + ! endif ! - else ! Original - ! - if (use_quadtree .eqv. .true.) then - ! - call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the netinfiltrationfile Netcdf format! !', 1) - ! - endif - ! - endif + endif ! ! 6) Read in data per type, either from ascii or general netcdf file ! From c654c0057b6198b5bbb2403771acbea064a011e2 Mon Sep 17 00:00:00 2001 From: Leynse Date: Wed, 4 Feb 2026 18:10:20 +0100 Subject: [PATCH 08/19] - Edit check - still possible that precip is forced but qinf=0 --- source/src/sfincs_domain.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 4ecb3258d..ca4c8a1dd 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2162,7 +2162,7 @@ subroutine initialize_infiltration() ! ! 5) Check whether infiltration input type (orignal vs netcdf) are correctly matched to grid type (regular vs quadtree) ! - if (inftype /= 'con') then !constant uniform works for both options + if (infiltration .and. inftype /= 'con') then !constant uniform works for both options ! if (netcdf_infiltration) then ! From 36a5aa8e085c95512e15c75011edf166856eea3d Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 12 Feb 2026 09:38:13 +0100 Subject: [PATCH 09/19] - Change names to infiltrationfile and infiltrationtype --- source/src/sfincs_data.f90 | 3 +-- source/src/sfincs_domain.f90 | 34 ++++++++++++++++------------------ source/src/sfincs_input.f90 | 4 ++-- 3 files changed, 19 insertions(+), 22 deletions(-) diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 4f3b74156..c367d1b3e 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -160,7 +160,7 @@ module sfincs_data character*256 :: netampfile character*256 :: netamprfile character*256 :: netspwfile - character*256 :: netinfiltrationfile + character*256 :: infiltrationfile character*256 :: scsfile character*256 :: smaxfile character*256 :: sefffile @@ -187,7 +187,6 @@ module sfincs_data character*3 :: outputtype_his character*3 :: utmzone character*3 :: inftype - character*3 :: netinftype integer :: epsg character*15 :: epsg_code integer :: nc_deflate_level diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index ca4c8a1dd..547617bf8 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2069,11 +2069,9 @@ subroutine initialize_infiltration() ! if (precip) then ! - if (netinfiltrationfile /= 'none') then + if (infiltrationfile /= 'none') then ! - ! inftype is user defined, keyword: 'netinftype' in sfincs.inp: - ! - inftype = netinftype + ! inftype is user defined, keyword: 'inftype' in sfincs.inp: ! ! inftype is either: c2d, cna, cnb, gai, hor ! 'inftype = con' is not relevant for netcdf input @@ -2150,10 +2148,10 @@ subroutine initialize_infiltration() write(logstr,'(a)')'Info : turning on infiltration from netcdf input file' call write_log(logstr, 0) ! - write(logstr,'(a,a)')'Info : reading netcdf infiltration file ', trim(netinfiltrationfile) + write(logstr,'(a,a)')'Info : reading netcdf infiltration file ', trim(infiltrationfile) call write_log(logstr, 0) ! - ok = check_file_exists(netinfiltrationfile, 'Infiltration netcdf file', .true.) + ok = check_file_exists(infiltrationfile, 'Infiltration netcdf file', .true.) ! write(logstr,'(a,a)')'Info : specified inftype is ', trim(inftype) call write_log(logstr, 0) @@ -2176,7 +2174,7 @@ subroutine initialize_infiltration() ! if (use_quadtree .eqv. .true.) then ! - call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the netinfiltrationfile Netcdf ormat! !', 1) + call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) endif ! endif @@ -2231,7 +2229,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'qinf' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, qinffield) !ncfile, varname, varout) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) !ncfile, varname, varout) ! else ! from separate qinffile - only binary: ! @@ -2265,7 +2263,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'scs' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, qinffield) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) ! else ! from separate scsfile - only binary: ! @@ -2300,7 +2298,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'smax' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, qinffield) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) ! else ! from separate smaxfile - only binary: ! @@ -2323,7 +2321,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'seff' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, scs_Se) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, scs_Se) ! else ! from separate sefffile - only binary: ! @@ -2347,7 +2345,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'ks' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, ksfield) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, ksfield) ! else ! from separate ksfile - only binary: ! @@ -2398,7 +2396,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'psi' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, GA_head) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, GA_head) ! else ! from separate psifile - only binary: ! @@ -2422,7 +2420,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'sigma' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, GA_sigma_max) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, GA_sigma_max) ! else ! from separate sigmafile - only binary: ! @@ -2446,7 +2444,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'ks' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, ksfield) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, ksfield) ! else ! from separate ksfile - only binary: ! @@ -2511,7 +2509,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'fc' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, horton_fc) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_fc) ! else ! from separate fcfile - only binary: ! @@ -2534,7 +2532,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'f0' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, horton_f0) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_f0) ! else ! from separate f0file - only binary: ! @@ -2557,7 +2555,7 @@ subroutine initialize_infiltration() ! ! Call the generic quadtree nc file reader function varname = 'kd' - call read_netcdf_quadtree_to_sfincs(netinfiltrationfile, varname, horton_kd) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_kd) ! else ! from separate kdfile - only binary: ! diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index a5b93ddb7..31b22daa9 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -209,8 +209,8 @@ subroutine read_sfincs_input() call read_char_input(500,'netampfile',netampfile,'none') call read_char_input(500,'netspwfile',netspwfile,'none') ! - call read_char_input(500,'netinfiltrationfile',netinfiltrationfile,'none') - call read_char_input(500,'netinftype',netinftype,'none') + call read_char_input(500,'infiltrationfile',infiltrationfile,'none') + call read_char_input(500,'infiltrationtype',inftype,'none') ! ! Output call read_char_input(500,'obsfile',obsfile,'none') From 30a0cb670a8bf205b81595ff490df795c3f576fb Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 12 Feb 2026 09:41:42 +0100 Subject: [PATCH 10/19] - Move whole subroutine 'initialize_infiltration' from sfincs_domain.f90 to sfincs_infiltration.f90 --- source/src/sfincs_domain.f90 | 583 +---------------------------- source/src/sfincs_infiltration.f90 | 583 +++++++++++++++++++++++++++++ 2 files changed, 585 insertions(+), 581 deletions(-) diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 547617bf8..9dbd31219 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -9,6 +9,7 @@ subroutine initialize_domain() ! use sfincs_data use quadtree + use sfincs_infiltration ! implicit none ! @@ -22,7 +23,7 @@ subroutine initialize_domain() ! call initialize_roughness() ! - call initialize_infiltration() + call initialize_infiltration() ! see: sfincs_infiltration.f90 ! call initialize_storage_volume() ! @@ -2017,586 +2018,6 @@ subroutine initialize_roughness() end subroutine - subroutine initialize_infiltration() - ! - use sfincs_data - use sfincs_ncinput - ! - implicit none - ! - integer :: nm - ! - logical :: ok - ! - character*256 :: varname - ! - ! INFILTRATION - ! - ! Infiltration only works when rainfall is activated ! If you want infiltration without rainfall, use a precip file with 0.0s - ! - ! Note, infiltration methods not designed to be stacked - ! - infiltration = .false. - netcdf_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) - ! 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' - ! - !!!!!!!!!!!!!!!!!!!!! - ! Initializing steps: - !!!!!!!!!!!!!!!!!!!!! - ! - ! 1) First we determine infiltration type - ! - if (precip) then - ! - if (infiltrationfile /= 'none') then - ! - ! inftype is user defined, keyword: 'inftype' in sfincs.inp: - ! - ! inftype is either: c2d, cna, cnb, gai, hor - ! 'inftype = con' is not relevant for netcdf input - ! - infiltration = .true. - netcdf_infiltration = .true. - ! - elseif (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 - ! - ! 2) We need cumprcp and cuminf - ! - allocate(cumprcp(np)) - cumprcp = 0.0 - ! - allocate(cuminf(np)) - cuminf = 0.0 - ! - ! 3) Now allocate and read spatially-varying inputs - ! - if (infiltration) then - ! - allocate(qinfmap(np)) - qinfmap = 0.0 - ! - endif - ! - ! 4) Pre-check whether netcdf infiltration file exists - once - ! - if (netcdf_infiltration) then - ! - write(logstr,'(a)')'Info : turning on infiltration from netcdf input file' - call write_log(logstr, 0) - ! - write(logstr,'(a,a)')'Info : reading netcdf infiltration file ', trim(infiltrationfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(infiltrationfile, 'Infiltration netcdf file', .true.) - ! - write(logstr,'(a,a)')'Info : specified inftype is ', trim(inftype) - call write_log(logstr, 0) - ! - endif - ! - ! 5) Check whether infiltration input type (orignal vs netcdf) are correctly matched to grid type (regular vs quadtree) - ! - if (infiltration .and. inftype /= 'con') then !constant uniform works for both options - ! - if (netcdf_infiltration) then - ! - if (use_quadtree .eqv. .false.) then - ! - call stop_sfincs('Error ! Netcdf infiltration input format can only be specified for quadtree mesh model !', 1) - ! - endif - ! - else ! Original - ! - if (use_quadtree .eqv. .true.) then - ! - call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) - endif - ! - endif - ! - endif - ! - ! 6) Read in data per type, either from ascii or general netcdf file - ! - if (inftype == 'con') then - ! - ! Spatially-uniform constant infiltration (specified as +mm/hr) - ! - ! Note : Input directly in sfincs.inp, so no file needs to be read - ! - write(logstr,'(a)')'Info : turning on spatially-uniform constant infiltration' - call write_log(logstr, 0) - ! - allocate(qinffield(np)) - ! - ! Note : qinf has already been converted to m/s in sfincs_input.f90 ! - ! - 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(logstr,'(a)')'Info : turning on spatially-varying constant infiltration' - call write_log(logstr, 0) - ! - allocate(qinffield(np)) - ! - qinffield = 0.0 - ! - ! Read spatially-varying infiltration (specified in +mm/hr) - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'qinf' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) !ncfile, varname, varout) - ! - else ! from separate qinffile - only binary: - ! - write(logstr,'(a,a)')'Info : reading infiltration file ', trim(qinffile) - call write_log(logstr, 0) - ! - ok = check_file_exists(qinffile, 'Infiltration qinf file', .true.) - ! - open(unit = 500, file = trim(qinffile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) - ! - endif - ! - ! Generic needed conversion: - ! - qinffield = qinffield / 3600 / 1000 ! convert to +m/s - ! - elseif (inftype == 'cna') then - ! - ! Spatially-varying infiltration with CN numbers (old) - ! - write(logstr,'(a)')'Info : turning on infiltration (via Curve Number method - A)' - call write_log(logstr, 0) - ! - allocate(qinffield(np)) - ! - qinffield = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'scs' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) - ! - else ! from separate scsfile - only binary: - ! - write(logstr,'(a,a)')'Info : reading scs file ',trim(scsfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(scsfile, 'Infiltration scs file', .true.) - ! - open(unit = 500, file = trim(scsfile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) - ! - endif - ! - ! Generic needed conversion: - ! - qinffield = qinffield * 0.0254 ! to m - ! already convert qinffield from inches to m here - ! - elseif (inftype == 'cnb') then - ! - ! Spatially-varying infiltration with CN numbers (new) - ! - write(logstr,'(a)')'Info : turning on infiltration (via Curve Number method - B)' - call write_log(logstr, 0) - ! - ! Allocate Smax - allocate(qinffield(np)) - qinffield = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'smax' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) - ! - else ! from separate smaxfile - only binary: - ! - write(logstr,'(a,a)')'Info : reading smax file ',trim(smaxfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(smaxfile, 'Infiltration smax file', .true.) - ! - open(unit = 500, file = trim(smaxfile), form = 'unformatted', access = 'stream') - read(500)qinffield - close(500) - ! - endif - ! - ! Allocate Se - allocate(scs_Se(np)) - scs_Se = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'seff' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, scs_Se) - ! - else ! from separate sefffile - only binary: - ! - write(logstr,'(a,a)')'Info : reading seff file ',trim(sefffile) - call write_log(logstr, 0) - ! - ok = check_file_exists(sefffile, 'Infiltration seff file', .true.) - ! - open(unit = 501, file = trim(sefffile), form = 'unformatted', access = 'stream') - read(501)scs_Se - close(501) - ! - endif - ! - ! Allocate Ks - ! - allocate(ksfield(np)) - ksfield = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'ks' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, ksfield) - ! - else ! from separate ksfile - only binary: - ! - write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) - ! - open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') - read(502)ksfield - close(502) - ! - endif - ! - ! Generic needed conversion: - ! - ! 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 - ! - call write_log('Info : turning on process infiltration (via Green-Ampt)', 0) - ! - ! Allocate suction head at the wetting front - ! - allocate(GA_head(np)) - GA_head = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'psi' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, GA_head) - ! - else ! from separate psifile - only binary: - ! - write(logstr,'(a,a)')'Info : reading psi file ',trim(psifile) - call write_log(logstr, 0) - ! - ok = check_file_exists(psifile, 'Infiltration psi file', .true.) - ! - open(unit = 500, file = trim(psifile), form = 'unformatted', access = 'stream') - read(500)GA_head - close(500) - ! - endif - ! - ! Allocate maximum soil moisture deficit - ! - allocate(GA_sigma_max(np)) - GA_sigma_max = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'sigma' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, GA_sigma_max) - ! - else ! from separate sigmafile - only binary: - ! - write(logstr,'(a,a)')'Info : reading sigma file ',trim(sigmafile) - call write_log(logstr, 0) - ! - ok = check_file_exists(sigmafile, 'Infiltration sigma file', .true.) - ! - open(unit = 501, file = trim(sigmafile), form = 'unformatted', access = 'stream') - read(501)GA_sigma_max - close(501) - ! - endif - ! - ! Allocate saturated hydraulic conductivity - ! - allocate(ksfield(np)) - ksfield = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'ks' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, ksfield) - ! - else ! from separate ksfile - only binary: - ! - write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) - ! - open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') - read(502)ksfield - close(502) - ! - endif - - ! - ! Generic needed conversion: - ! - ! 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 support variables: - ! - allocate(qinffield(np)) - qinffield(nm) = 0.0 - ! - elseif (inftype == 'hor') then - ! - ! Spatially-varying infiltration with the modified Horton Equation - ! - call write_log('Info : turning on process infiltration (via modified Horton)', 0) - ! - ! Horton: final infiltration capacity (fc) - ! Note that qinffield = horton_fc (/3600/1000, see below) - ! - allocate(horton_fc(np)) - horton_fc = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'fc' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_fc) - ! - else ! from separate fcfile - only binary: - ! - write(logstr,'(a,a)')'Info : reading fc file ',trim(fcfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(fcfile, 'Infiltration fc file', .true.) - ! - open(unit = 500, file = trim(fcfile), form = 'unformatted', access = 'stream') - read(500)horton_fc - close(500) - ! - endif - ! - ! Horton: initial infiltration capacity (f0) - allocate(horton_f0(np)) - horton_f0 = 0.0 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'f0' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_f0) - ! - else ! from separate f0file - only binary: - ! - write(logstr,'(a,a)')'Info : reading f0 file ',trim(f0file) - call write_log(logstr, 0) - ! - ok = check_file_exists(f0file, 'Infiltration f0 file', .true.) - ! - open(unit = 501, file = trim(f0file), form = 'unformatted', access = 'stream') - read(501)horton_f0 - close(501) - ! - endif - ! - ! 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 - ! - if (netcdf_infiltration) then - ! - ! Call the generic quadtree nc file reader function - varname = 'kd' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_kd) - ! - else ! from separate kdfile - only binary: - ! - write(logstr,'(a,a)')'Info : reading kd file ',trim(kdfile) - call write_log(logstr, 0) - ! - ok = check_file_exists(kdfile, 'Infiltration kd file', .true.) - ! - open(unit = 502, file = trim(kdfile), form = 'unformatted', access = 'stream') - read(502)horton_kd - close(502) - ! - endif - ! - write(logstr,'(a,a)')'Info : Using constant recovery rate that is based on constant factor relative to ',trim(kdfile) - call write_log(logstr, 0) - ! - ! Generic needed conversion: - ! - ! Prescribe the current estimate (for output only; initial capacity) - qinffield = horton_f0/3600/1000 - ! - ! Allocate support variables: - ! - ! Estimate of time - allocate(rain_T1(np)) - rain_T1 = 0.0 - ! - endif - ! - else - ! - ! Overrule input - ! - store_cumulative_precipitation = .false. - ! - endif - ! - end subroutine - - subroutine initialize_storage_volume() ! use sfincs_data diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 index 57a4e2f20..4dd7acc58 100644 --- a/source/src/sfincs_infiltration.f90 +++ b/source/src/sfincs_infiltration.f90 @@ -1,7 +1,590 @@ module sfincs_infiltration + use sfincs_log + use sfincs_error + contains + subroutine initialize_infiltration() + ! + use sfincs_data + use sfincs_ncinput + ! + implicit none + ! + integer :: nm + ! + logical :: ok + ! + character*256 :: varname + ! + ! INFILTRATION + ! + ! Infiltration only works when rainfall is activated ! If you want infiltration without rainfall, use a precip file with 0.0s + ! + ! Note, infiltration methods not designed to be stacked + ! + infiltration = .false. + netcdf_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) + ! 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' + ! + !!!!!!!!!!!!!!!!!!!!! + ! Initializing steps: + !!!!!!!!!!!!!!!!!!!!! + ! + ! 1) First we determine infiltration type + ! + if (precip) then + ! + if (infiltrationfile /= 'none') then + ! + ! inftype is user defined, keyword: 'inftype' in sfincs.inp: + ! + ! inftype is either: c2d, cna, cnb, gai, hor + ! 'inftype = con' is not relevant for netcdf input + ! + infiltration = .true. + netcdf_infiltration = .true. + ! + elseif (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 + ! + ! 2) We need cumprcp and cuminf + ! + allocate(cumprcp(np)) + cumprcp = 0.0 + ! + allocate(cuminf(np)) + cuminf = 0.0 + ! + ! 3) Now allocate and read spatially-varying inputs + ! + if (infiltration) then + ! + allocate(qinfmap(np)) + qinfmap = 0.0 + ! + endif + ! + ! 4) Pre-check whether netcdf infiltration file exists - once + ! + if (netcdf_infiltration) then + ! + write(logstr,'(a)')'Info : turning on infiltration from netcdf input file' + call write_log(logstr, 0) + ! + write(logstr,'(a,a)')'Info : reading netcdf infiltration file ', trim(infiltrationfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(infiltrationfile, 'Infiltration netcdf file', .true.) + ! + write(logstr,'(a,a)')'Info : specified inftype is ', trim(inftype) + call write_log(logstr, 0) + ! + endif + ! + ! 5) Check whether infiltration input type (orignal vs netcdf) are correctly matched to grid type (regular vs quadtree) + ! + if (infiltration .and. inftype /= 'con') then !constant uniform works for both options + ! + if (netcdf_infiltration) then + ! + if (use_quadtree .eqv. .false.) then + ! + call stop_sfincs('Error ! Netcdf infiltration input format can only be specified for quadtree mesh model !', 1) + ! + endif + ! + else ! Original + ! + if (use_quadtree .eqv. .true.) then + ! + call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) + endif + ! + endif + ! + endif + ! + ! 6) Read in data per type, either from ascii or general netcdf file + ! + if (inftype == 'con') then + ! + ! Spatially-uniform constant infiltration (specified as +mm/hr) + ! + ! Note : Input directly in sfincs.inp, so no file needs to be read + ! + write(logstr,'(a)')'Info : turning on spatially-uniform constant infiltration' + call write_log(logstr, 0) + ! + allocate(qinffield(np)) + ! + ! Note : qinf has already been converted to m/s in sfincs_input.f90 ! + ! + 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(logstr,'(a)')'Info : turning on spatially-varying constant infiltration' + call write_log(logstr, 0) + ! + allocate(qinffield(np)) + ! + qinffield = 0.0 + ! + ! Read spatially-varying infiltration (specified in +mm/hr) + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'qinf' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) !ncfile, varname, varout) + ! + else ! from separate qinffile - only binary: + ! + write(logstr,'(a,a)')'Info : reading infiltration file ', trim(qinffile) + call write_log(logstr, 0) + ! + ok = check_file_exists(qinffile, 'Infiltration qinf file', .true.) + ! + open(unit = 500, file = trim(qinffile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + endif + ! + ! Generic needed conversion: + ! + qinffield = qinffield / 3600 / 1000 ! convert to +m/s + ! + elseif (inftype == 'cna') then + ! + ! Spatially-varying infiltration with CN numbers (old) + ! + write(logstr,'(a)')'Info : turning on infiltration (via Curve Number method - A)' + call write_log(logstr, 0) + ! + allocate(qinffield(np)) + ! + qinffield = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'scs' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) + ! + else ! from separate scsfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading scs file ',trim(scsfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(scsfile, 'Infiltration scs file', .true.) + ! + open(unit = 500, file = trim(scsfile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + endif + ! + ! Generic needed conversion: + ! + qinffield = qinffield * 0.0254 ! to m + ! already convert qinffield from inches to m here + ! + elseif (inftype == 'cnb') then + ! + ! Spatially-varying infiltration with CN numbers (new) + ! + write(logstr,'(a)')'Info : turning on infiltration (via Curve Number method - B)' + call write_log(logstr, 0) + ! + ! Allocate Smax + allocate(qinffield(np)) + qinffield = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'smax' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, qinffield) + ! + else ! from separate smaxfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading smax file ',trim(smaxfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(smaxfile, 'Infiltration smax file', .true.) + ! + open(unit = 500, file = trim(smaxfile), form = 'unformatted', access = 'stream') + read(500)qinffield + close(500) + ! + endif + ! + ! Allocate Se + allocate(scs_Se(np)) + scs_Se = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'seff' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, scs_Se) + ! + else ! from separate sefffile - only binary: + ! + write(logstr,'(a,a)')'Info : reading seff file ',trim(sefffile) + call write_log(logstr, 0) + ! + ok = check_file_exists(sefffile, 'Infiltration seff file', .true.) + ! + open(unit = 501, file = trim(sefffile), form = 'unformatted', access = 'stream') + read(501)scs_Se + close(501) + ! + endif + ! + ! Allocate Ks + ! + allocate(ksfield(np)) + ksfield = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'ks' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, ksfield) + ! + else ! from separate ksfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) + ! + open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') + read(502)ksfield + close(502) + ! + endif + ! + ! Generic needed conversion: + ! + ! 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 + ! + call write_log('Info : turning on process infiltration (via Green-Ampt)', 0) + ! + ! Allocate suction head at the wetting front + ! + allocate(GA_head(np)) + GA_head = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'psi' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, GA_head) + ! + else ! from separate psifile - only binary: + ! + write(logstr,'(a,a)')'Info : reading psi file ',trim(psifile) + call write_log(logstr, 0) + ! + ok = check_file_exists(psifile, 'Infiltration psi file', .true.) + ! + open(unit = 500, file = trim(psifile), form = 'unformatted', access = 'stream') + read(500)GA_head + close(500) + ! + endif + ! + ! Allocate maximum soil moisture deficit + ! + allocate(GA_sigma_max(np)) + GA_sigma_max = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'sigma' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, GA_sigma_max) + ! + else ! from separate sigmafile - only binary: + ! + write(logstr,'(a,a)')'Info : reading sigma file ',trim(sigmafile) + call write_log(logstr, 0) + ! + ok = check_file_exists(sigmafile, 'Infiltration sigma file', .true.) + ! + open(unit = 501, file = trim(sigmafile), form = 'unformatted', access = 'stream') + read(501)GA_sigma_max + close(501) + ! + endif + ! + ! Allocate saturated hydraulic conductivity + ! + allocate(ksfield(np)) + ksfield = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'ks' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, ksfield) + ! + else ! from separate ksfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading ks file ',trim(ksfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(ksfile, 'Infiltration ks file', .true.) + ! + open(unit = 502, file = trim(ksfile), form = 'unformatted', access = 'stream') + read(502)ksfield + close(502) + ! + endif + + ! + ! Generic needed conversion: + ! + ! 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 support variables: + ! + allocate(qinffield(np)) + qinffield(nm) = 0.0 + ! + elseif (inftype == 'hor') then + ! + ! Spatially-varying infiltration with the modified Horton Equation + ! + call write_log('Info : turning on process infiltration (via modified Horton)', 0) + ! + ! Horton: final infiltration capacity (fc) + ! Note that qinffield = horton_fc (/3600/1000, see below) + ! + allocate(horton_fc(np)) + horton_fc = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'fc' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_fc) + ! + else ! from separate fcfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading fc file ',trim(fcfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(fcfile, 'Infiltration fc file', .true.) + ! + open(unit = 500, file = trim(fcfile), form = 'unformatted', access = 'stream') + read(500)horton_fc + close(500) + ! + endif + ! + ! Horton: initial infiltration capacity (f0) + allocate(horton_f0(np)) + horton_f0 = 0.0 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'f0' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_f0) + ! + else ! from separate f0file - only binary: + ! + write(logstr,'(a,a)')'Info : reading f0 file ',trim(f0file) + call write_log(logstr, 0) + ! + ok = check_file_exists(f0file, 'Infiltration f0 file', .true.) + ! + open(unit = 501, file = trim(f0file), form = 'unformatted', access = 'stream') + read(501)horton_f0 + close(501) + ! + endif + ! + ! 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 + ! + if (netcdf_infiltration) then + ! + ! Call the generic quadtree nc file reader function + varname = 'kd' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, horton_kd) + ! + else ! from separate kdfile - only binary: + ! + write(logstr,'(a,a)')'Info : reading kd file ',trim(kdfile) + call write_log(logstr, 0) + ! + ok = check_file_exists(kdfile, 'Infiltration kd file', .true.) + ! + open(unit = 502, file = trim(kdfile), form = 'unformatted', access = 'stream') + read(502)horton_kd + close(502) + ! + endif + ! + write(logstr,'(a,a)')'Info : Using constant recovery rate that is based on constant factor relative to ',trim(kdfile) + call write_log(logstr, 0) + ! + ! Generic needed conversion: + ! + ! Prescribe the current estimate (for output only; initial capacity) + qinffield = horton_f0/3600/1000 + ! + ! Allocate support variables: + ! + ! Estimate of time + allocate(rain_T1(np)) + rain_T1 = 0.0 + ! + endif + ! + else + ! + ! Overrule input + ! + store_cumulative_precipitation = .false. + ! + endif + ! + end subroutine + + subroutine update_infiltration_map(dt, tloop) ! ! Update infiltration rates in each grid cell From 0dc07201b981cd109e68bc79cc7229e5c76ec882 Mon Sep 17 00:00:00 2001 From: Leynse Date: Thu, 12 Feb 2026 09:52:01 +0100 Subject: [PATCH 11/19] - Add check whether specified infiltrationtype is correct --- source/src/sfincs_infiltration.f90 | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 index 4dd7acc58..bce4af73f 100644 --- a/source/src/sfincs_infiltration.f90 +++ b/source/src/sfincs_infiltration.f90 @@ -18,6 +18,11 @@ subroutine initialize_infiltration() ! character*256 :: varname ! + character(len=3), parameter :: allowed_types(5) = & + ['c2d', 'cna', 'cnb', 'gai', 'hor'] + + logical :: inftype_exists + ! ! INFILTRATION ! ! Infiltration only works when rainfall is activated ! If you want infiltration without rainfall, use a precip file with 0.0s @@ -64,8 +69,25 @@ subroutine initialize_infiltration() ! inftype is either: c2d, cna, cnb, gai, hor ! 'inftype = con' is not relevant for netcdf input ! - infiltration = .true. - netcdf_infiltration = .true. + ! Check if specified type is correct + ! + inftype_exists = any(inftype == allowed_types) + ! + if (inftype_exists) then + ! + infiltration = .true. + netcdf_infiltration = .true. + ! + write(logstr,'(a,a)')'Info : specified infiltrationtype is ', trim(inftype) + call write_log(logstr, 0) + ! + else + ! + write(logstr,*)'Error : infiltration input type ',trim(inftype),' is not part of supported types c2d cna cnb gai hor !' + call stop_sfincs(trim(logstr), 1) + ! + end if + ! ! elseif (qinf > 0.0) then ! @@ -141,9 +163,6 @@ subroutine initialize_infiltration() ! ok = check_file_exists(infiltrationfile, 'Infiltration netcdf file', .true.) ! - write(logstr,'(a,a)')'Info : specified inftype is ', trim(inftype) - call write_log(logstr, 0) - ! endif ! ! 5) Check whether infiltration input type (orignal vs netcdf) are correctly matched to grid type (regular vs quadtree) From 2a1a97489877705606b9b1573a3779a512f4ede2 Mon Sep 17 00:00:00 2001 From: Tim Leijnse Date: Fri, 6 Mar 2026 11:13:38 +0100 Subject: [PATCH 12/19] - Added version for reading in netcdf quadtree manningfile (#271) * - Added version for reading in netcdf quadtree manningfile - For now: input argument is still manningfile, internally checked whether it is a netcdf file or binary - For now: expected variable in the netcdf file is 'manning' - For now: tested only with dummy of stoarge nc file - Full testing pending on having ncinput example * - Add non-stopping warning if a manningfile is prescribed while there is also a subgridfile, because then it will be unused! * - Bugfix in general read_netcdf_quadtree_to_sfincs routine * - Bugfix that the netcdf manningfile is actually read in * - Now define 'rghfield' in sfincs_data.f90 rather than sfincs_domain.f90, so we can use that variable for output to netcdf file * - Add 'manning' as netcdf output in case not a subgrid model, and if manningfile supplied (either binary for regular model, or as netcdf for quadtree - Tested and working for new quadtree netcdf manningfile input --- source/src/sfincs_data.f90 | 1 + source/src/sfincs_domain.f90 | 39 +++++++++++++++++--- source/src/sfincs_lib.f90 | 2 +- source/src/sfincs_ncinput.F90 | 4 ++- source/src/sfincs_ncoutput.F90 | 65 +++++++++++++++++++++++++++++++++- 5 files changed, 103 insertions(+), 8 deletions(-) diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index c778c4925..ac55df396 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -373,6 +373,7 @@ module sfincs_data real*4, dimension(:), allocatable, target :: z_yz real*4, dimension(:), allocatable :: cell_area_m2 real*4, dimension(:), allocatable :: nuvisc + real*4, dimension(:), allocatable :: rghfield ! ! UV-points ! diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 8b789fb50..b5bc0a447 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -1945,16 +1945,19 @@ subroutine initialize_boundaries() subroutine initialize_roughness() ! use sfincs_data + use sfincs_ncinput ! implicit none ! - real*4, dimension(:), allocatable :: rghfield ! integer :: ip integer :: nm integer :: nmu logical :: ok ! + integer :: nchar + character*256 :: varname + ! ! FRICTION COEFFICIENTS (only for regular bathymetry, as for subgrid the Manning's n values are stored in the tables) ! if (.not. subgrid) then @@ -1963,19 +1966,35 @@ subroutine initialize_roughness() ! gn2uv = 9.81*0.02*0.02 ! - if (manningfile(1:4) /= 'none') then + if (manningfile(1:4) /= 'none') then ! ! Read spatially-varying friction + ! File is either binary or netcdf ! allocate(rghfield(np)) + ! write(logstr,'(a,a)')'Info : reading roughness file ',trim(manningfile) call write_log(logstr, 0) ! + nchar = len_trim(manningfile) + ! ok = check_file_exists(manningfile, 'Roughness file', .true.) ! - open(unit = 500, file = trim(manningfile), form = 'unformatted', access = 'stream') - read(500)rghfield - close(500) + if (manningfile(nchar - 1 : nchar) == 'nc') then + ! + ! Call the generic quadtree nc file reader function + varname = 'manning' + call read_netcdf_quadtree_to_sfincs(manningfile, varname, rghfield) !ncfile, varname, varout) + ! + else + ! + ! Read from binary file + ! + open(unit = 500, file = trim(manningfile), form = 'unformatted', access = 'stream') + read(500)rghfield + close(500) + ! + endif ! do ip = 1, npuv nm = uv_index_z_nm(ip) @@ -2012,6 +2031,16 @@ subroutine initialize_roughness() enddo ! endif + else + ! + ! Give warning if manningfile is supplied, but also a subgrid file + ! + if (manningfile(1:4) /= 'none') then + ! + call write_log('Warning : manningfile input will be ignored because SFINCS will use the friction information from sbgfile!', 1) + ! + endif + ! endif ! end subroutine diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index 898905cc2..07aeb92fd 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -92,7 +92,7 @@ function sfincs_initialize() result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = "$Rev: v2.3.1 mt. Faber" + build_revision = "$Rev: v2.3.1 mt. Faber+branch-270" build_date = "$Date: 2025-02-04" ! call write_log('', 1) diff --git a/source/src/sfincs_ncinput.F90 b/source/src/sfincs_ncinput.F90 index 7f4ef894c..96687cd25 100644 --- a/source/src/sfincs_ncinput.F90 +++ b/source/src/sfincs_ncinput.F90 @@ -299,7 +299,9 @@ subroutine read_netcdf_quadtree_to_sfincs(ncfile, varname, var) ! nm = index_sfincs_in_quadtree(ip) ! - var(nm) = vartmp(ip) + if (nm>0) then + var(nm) = vartmp(ip) + endif ! enddo ! diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index 7042a0424..3fcfbe9d6 100644 --- a/source/src/sfincs_ncoutput.F90 +++ b/source/src/sfincs_ncoutput.F90 @@ -23,6 +23,7 @@ module sfincs_ncoutput integer :: fwx_varid, fwy_varid, beta_varid, snapwavedepth_varid integer :: zsm_varid, tsunami_arrival_time_varid integer :: inp_varid, total_runtime_varid, average_dt_varid, status_varid + integer :: manning_varid integer :: pnonh_varid integer :: subgridslope_varid ! @@ -241,6 +242,16 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, 'coordinates', 'x y')) endif ! + if (.not. subgrid) then + NF90(nf90_def_var(map_file%ncid, 'manning', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%manning_varid)) ! bed level in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%manning_varid, 1, 1, nc_deflate_level)) ! deflate + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, 'units', 's/m^1/3')) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, 'standard_name', 'manning')) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, 'long_name', 'manning_roughness')) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, 'coordinates', 'x y')) + endif + ! if (subgrid .and. store_hsubgrid .and. store_hmean) then ! ! The subgrid slope (zmax - zmin) / sqrt(A) is used for making high-res flood maps @@ -769,6 +780,25 @@ subroutine ncoutput_regular_map_init() ! endif ! + ! Write Manning (only non-subgrid model) + ! + if (.not. subgrid .and. manning2d) then + ! + zsg = FILL_VALUE + ! + do nm = 1, np + ! + n = z_index_z_n(nm) + m = z_index_z_m(nm) + ! + zsg(m, n) = rghfield(nm) ! gn2uv is on uv-points, but rghfield is in center + ! + enddo + ! + NF90(nf90_put_var(map_file%ncid, map_file%manning_varid, zsg, (/1, 1/))) + ! + endif + ! ! Write infiltration map ! if (infiltration) then @@ -815,7 +845,7 @@ subroutine ncoutput_quadtree_map_init() ! implicit none ! - integer :: nm, nmq, n, m, nn, ntmx, n_nodes, n_faces, iref + integer :: nm, nmq, nmu1, num1, n, m, nn, ntmx, n_nodes, n_faces, iref real*4 :: dxx, dyy ! real, dimension(:), allocatable :: nodes_x @@ -978,6 +1008,17 @@ subroutine ncoutput_quadtree_map_init() NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, 'standard_name', 'altitude')) NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, 'long_name', 'bed_level_above_reference_level')) ! + if (.not. subgrid) then + ! + NF90(nf90_def_var(map_file%ncid, 'manning', NF90_FLOAT, (/map_file%nmesh2d_face_dimid/), map_file%manning_varid)) ! bed level in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%manning_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, 'units', 's/m^1/3')) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, 'standard_name', 'manning')) + NF90(nf90_put_att(map_file%ncid, map_file%manning_varid, 'long_name', 'manning_roughness')) + ! + endif + ! if (subgrid .and. store_hsubgrid .and. store_hmean) then ! ! The subgrid slope (zmax - zmin) / sqrt(A) is used for making high-res flood maps @@ -1463,6 +1504,28 @@ subroutine ncoutput_quadtree_map_init() ! endif ! + ! Write Manning (only non-subgrid model) + ! + if (.not. subgrid .and. manning2d) then + ! + vtmp = FILL_VALUE + ! + do nmq = 1, quadtree_nr_points + ! + nm = index_sfincs_in_quadtree(nmq) + ! + if (nm>0) then + ! + vtmp(nmq) = rghfield(nm) ! gn2uv is on uv-points, but rghfield is in center + ! + endif + ! + enddo + ! + NF90(nf90_put_var(map_file%ncid, map_file%manning_varid, vtmp)) + ! + endif + ! ! Write infiltration map ! vtmp = FILL_VALUE From 167387b1f6b5dd83b0a70188392a975ca147b9be Mon Sep 17 00:00:00 2001 From: Leynse Date: Fri, 6 Mar 2026 15:51:13 +0100 Subject: [PATCH 13/19] - Bump date --- source/src/sfincs_lib.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index bdeb3e139..11318103e 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -93,7 +93,7 @@ function sfincs_initialize() result(ierr) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! build_revision = "$Rev: v2.3.1 mt. Faber+branch-272" - build_date = "$Date: 2025-02-04" + build_date = "$Date: 2026-03-06" ! call write_log('', 1) call write_log('------------ Welcome to SFINCS ------------', 1) From a0ed4a044a78cd8e522654237ed25ada8e423ac0 Mon Sep 17 00:00:00 2001 From: Kees Nederhoff Date: Wed, 18 Mar 2026 07:22:48 -0700 Subject: [PATCH 14/19] Progress with redoing continuity logic, adding new infiltration flavor, adding option to mimic drainage --- source/src/sfincs_continuity.f90 | 165 +++++++++++++----- source/src/sfincs_data.f90 | 21 +++ source/src/sfincs_domain.f90 | 80 ++++++++- source/src/sfincs_infiltration.f90 | 258 +++++++++++++++++++++++------ source/src/sfincs_input.f90 | 7 +- source/src/sfincs_lib.f90 | 18 +- source/src/sfincs_meteo.f90 | 38 ++--- source/src/sfincs_openacc.f90 | 10 +- 8 files changed, 457 insertions(+), 140 deletions(-) diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90 index 251951150..bcfc9699a 100644 --- a/source/src/sfincs_continuity.f90 +++ b/source/src/sfincs_continuity.f90 @@ -1,7 +1,62 @@ module sfincs_continuity contains - + + subroutine update_continuity(t, dt, tloopsrc, tloopinf, tloopcont) + ! + ! Unified continuity update: orchestrates all water balance terms + ! + ! A. Discharges, sources and sinks => computed in sfincs_discharges + ! B. Hydrodynamic fluxes (q) => already computed in sfincs_momentum + ! C. Main adjustments: + ! 1. Rainfall (+) => already computed in sfincs_meteo (prcp) + ! 2. Infiltration (-) => computed in sfincs_infiltration (qinfmap) + ! (includes: con, c2d, cna, cnb, gai, hor, bkt flavors) + ! 3. Drainage mimic (-) => simple constant rate (qdrain_rate) + ! 4. External source/sink qext (+/-) => set via BMI coupling + ! 5. Storage volume => depression storage (subgrid only) + ! + ! compute_water_levels then applies all terms to update zs/z_volume: + ! A. Discharges (+/-) => outside main loop + ! B. Hydrodynamic fluxes => div(q) * dt + ! C1. Rainfall (+) => prcp * dt + ! C2. Infiltration (-) => qinfmap * dt + ! C3. Drainage mimic (-) => qdrain_rate * dt + ! C4. External source/sink (+/-) => qext * dt + ! C5. Storage volume => absorbs excess volume + ! + use sfincs_data + use sfincs_infiltration + use sfincs_discharges + ! + implicit none + ! + real*8 :: t + real*4 :: dt + real :: tloopsrc + real :: tloopinf + real :: tloopcont + ! + ! A. Update discharges, sources and sinks + ! + call update_discharges(t, dt, tloopsrc) + ! + ! C2. Compute infiltration rates => qinfmap (all flavors including bucket) + ! + if (infiltration) then + call update_infiltration_map(dt, tloopinf) + endif + ! + ! C1, C3, C4, C5: rainfall, drainage mimic, qext, storage_volume + ! => nothing to compute, these are direct rates applied in compute_water_levels + ! + ! B + C: Update water levels (applies all terms) + ! + call compute_water_levels(t, dt, tloopcont) + ! + end subroutine + + subroutine compute_water_levels(t, dt, tloop) ! use sfincs_data @@ -75,7 +130,7 @@ subroutine compute_water_levels_regular(dt,t) ! endif ! - !$acc parallel present( kcs, zs, zb, netprcp, prcp, q, qext, zsmax, zsm, maxzsm, & + !$acc parallel present( kcs, zs, zb, prcp, q, qext, qinfmap, qdrain_rate, zsmax, zsm, maxzsm, & !$acc z_flags_iref, uv_flags_iref, & !$acc z_index_uv_md, z_index_uv_nd, z_index_uv_mu, z_index_uv_nu, & !$acc dxm, dxrm, dyrm, dxminv, dxrinv, dyrinv, cell_area_m2, cell_area, & @@ -85,24 +140,24 @@ subroutine compute_water_levels_regular(dt,t) ! First discharges (don't do this parallel, as it's probably not worth it) ! if (nsrcdrn > 0) then - ! + ! !$acc loop do isrc = 1, nsrcdrn - ! + ! nm = nmindsrc(isrc) - ! + ! if (crsgeo) then - ! + ! zs(nmindsrc(isrc)) = max(zs(nm) + qtsrc(isrc) * dt / cell_area_m2(nm), zb(nm)) - ! + ! else - ! + ! zs(nmindsrc(isrc)) = max(zs(nm) + qtsrc(isrc) * dt / cell_area(z_flags_iref(nm)), zb(nm)) - ! + ! endif - ! + ! enddo - ! + ! endif ! !$omp parallel & @@ -110,21 +165,31 @@ subroutine compute_water_levels_regular(dt,t) !$omp do schedule ( dynamic, 256 ) !$acc loop gang vector do nm = 1, np - ! + ! if (kcs(nm) == 1) then ! Regular point + ! + ! C1. Rainfall (+) ! if (precip) then - ! - zs(nm) = zs(nm) + netprcp(nm) * dt - ! + zs(nm) = zs(nm) + prcp(nm) * dt + endif + ! + ! C2. Infiltration (-) (includes all flavors: con, c2d, cna, cnb, gai, hor, bkt) + ! + if (infiltration) then + zs(nm) = zs(nm) - qinfmap(nm) * dt endif ! + ! C3. Drainage mimic (-) + ! + if (use_drainage_mimic) then + zs(nm) = zs(nm) - qdrain_rate(nm) * dt + endif + ! + ! C4. External source/sink (+/-) + ! if (use_qext) then - ! - ! Add external source (e.g. from XMI coupling) - ! - zs(nm) = zs(nm) + qext(nm) * dt - ! + zs(nm) = zs(nm) + qext(nm) * dt endif ! nmd = z_index_uv_md(nm) @@ -329,9 +394,9 @@ subroutine compute_water_levels_subgrid(dt,t) !$omp do schedule ( dynamic, 256 ) !$acc parallel present( kcs, zs, zs0, zb, z_volume, zsmax, zsm, maxzsm, zsderv, & !$acc subgrid_z_zmin, subgrid_z_zmax, subgrid_z_dep, subgrid_z_volmax, & - !$acc netprcp, prcp, q, qext, z_flags_iref, uv_flags_iref, & + !$acc prcp, q, qext, qinfmap, qdrain_rate, z_flags_iref, uv_flags_iref, & !$acc z_index_uv_md, z_index_uv_nd, z_index_uv_mu, z_index_uv_nu, & - !$acc dxm, dxrm, dyrm, dxminv, dxrinv, dyrinv, cell_area_m2, cell_area, & + !$acc dxm, dxrm, dyrm, dxminv, dxrinv, dyrinv, cell_area_m2, cell_area, & !$acc z_index_wavemaker, wavemaker_uvmean, wavemaker_nmd, wavemaker_nmu, wavemaker_ndm, wavemaker_num, storage_volume) !$acc loop gang vector do nm = 1, np @@ -470,32 +535,40 @@ subroutine compute_water_levels_subgrid(dt,t) ! endif ! - if (precip .or. use_qext) then - ! - dzsdt = 0.0 - ! - if (precip) then - ! - ! Add nett rainfall - ! - dzsdt = dzsdt + netprcp(nm) - ! - endif - ! - if (use_qext) then - ! - ! Add external source (e.g. from XMI coupling) - ! - dzsdt = dzsdt + qext(nm) - ! - endif - ! - ! dzsdt is still in m/s, so multiply with a * dt to get m^3 - ! - dvol = dvol + dzsdt * a * dt - ! + dzsdt = 0.0 + ! + ! C1. Rainfall (+) + ! + if (precip) then + dzsdt = dzsdt + prcp(nm) endif ! + ! C2. Infiltration (-) (includes all flavors: con, c2d, cna, cnb, gai, hor, bkt) + ! + if (infiltration) then + dzsdt = dzsdt - qinfmap(nm) + endif + ! + ! C3. Drainage mimic (-) + ! + if (use_drainage_mimic) then + dzsdt = dzsdt - qdrain_rate(nm) + endif + ! + ! C4. External source/sink (+/-) + ! + if (use_qext) then + dzsdt = dzsdt + qext(nm) + endif + ! + ! dzsdt is still in m/s, so multiply with a * dt to get m^3 + ! + if (dzsdt /= 0.0) then + dvol = dvol + dzsdt * a * dt + endif + ! + ! C5. Storage volume + ! if (use_storage_volume) then ! ! If water enters the cell through a point discharge, it will NOT end up in storage volume ! diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index b3108611f..111c4cc5c 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -170,6 +170,8 @@ module sfincs_data character*256 :: f0file character*256 :: fcfile character*256 :: kdfile + character*256 :: drainagefile + character*256 :: bucketfile character*256 :: z0lfile character*256 :: wvmfile character*256 :: qtrfile @@ -408,6 +410,20 @@ module sfincs_data ! real*4, dimension(:), allocatable :: storage_volume ! Storage volume green infra ! + ! Drainage mimic - constant removal rate representing subsurface drainage + ! + logical :: use_drainage_mimic = .false. + real*4 :: qdrain_uniform = 0.0 ! uniform drainage rate (mm/hr input, stored as m/s) + real*4, dimension(:), allocatable :: qdrain_rate ! drainage rate per cell (m/s) + ! + ! Bucket model - finite capacity reservoir with linear drainage + ! + logical :: use_bucket_model = .false. + real*4, dimension(:), allocatable :: bucket_volume ! current storage (m) + real*4, dimension(:), allocatable :: bucket_capacity ! max capacity S_max (m) + real*4, dimension(:), allocatable :: bucket_k ! drainage coefficient (1/s) + real*4, dimension(:), allocatable :: bucket_drain_rate ! net removal from surface this step (m/s) + ! ! Wind reduction for spiderweb winds ! real*4, dimension(:,:), allocatable :: z0land ! z0 values over land for spiderweb wind speed reduction @@ -924,6 +940,11 @@ subroutine finalize_parameters() if(allocated(qinffield)) deallocate(qinffield) if(allocated(ksfield)) deallocate(ksfield) if(allocated(scs_Se)) deallocate(scs_Se) + if(allocated(qdrain_rate)) deallocate(qdrain_rate) + if(allocated(bucket_volume)) deallocate(bucket_volume) + if(allocated(bucket_capacity)) deallocate(bucket_capacity) + if(allocated(bucket_k)) deallocate(bucket_k) + if(allocated(bucket_drain_rate)) deallocate(bucket_drain_rate) if(allocated(nuvisc)) deallocate(nuvisc) ! ! Boundary velocity points diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 6a1d2e229..7cfce450b 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -23,7 +23,9 @@ subroutine initialize_domain() ! call initialize_roughness() ! - call initialize_infiltration() ! see: sfincs_infiltration.f90 + call initialize_infiltration() ! see: sfincs_infiltration.f90 (includes bucket model if bucketfile specified) + ! + call initialize_drainage_mimic() ! call initialize_storage_volume() ! @@ -2047,6 +2049,82 @@ subroutine initialize_roughness() end subroutine + subroutine initialize_drainage_mimic() + ! + use sfincs_data + use sfincs_ncinput + ! + implicit none + ! + integer :: nm + integer :: nchar + logical :: ok + character*256 :: varname + ! + ! Check if drainage mimic is enabled + ! + if (qdrain_uniform > 0.0 .or. drainagefile /= 'none') then + ! + use_drainage_mimic = .true. + ! + allocate(qdrain_rate(np)) + ! + if (drainagefile /= 'none') then + ! + ! Spatially-varying drainage rate + ! + write(logstr,'(a)')'Info : turning on drainage mimic (spatially-varying)' + call write_log(logstr, 0) + ! + nchar = len_trim(drainagefile) + ok = check_file_exists(drainagefile, 'Drainage file', .true.) + ! + if (drainagefile(nchar - 1 : nchar) == 'nc') then + ! + varname = 'drainage_rate' + call read_netcdf_quadtree_to_sfincs(drainagefile, varname, qdrain_rate) + ! + ! Convert from mm/hr to m/s + ! + qdrain_rate = qdrain_rate / 3600.0 / 1000.0 + ! + else + ! + ! Read from binary file (assumed to be in mm/hr) + ! + open(unit = 500, file = trim(drainagefile), form = 'unformatted', access = 'stream') + read(500)qdrain_rate + close(500) + ! + ! Convert from mm/hr to m/s + ! + qdrain_rate = qdrain_rate / 3600.0 / 1000.0 + ! + endif + ! + else + ! + ! Uniform drainage rate (already converted to m/s in sfincs_input.f90) + ! + write(logstr,'(a,f10.4,a)')'Info : turning on drainage mimic (uniform, ', qdrain_uniform * 3600.0 * 1000.0, ' mm/hr)' + call write_log(logstr, 0) + ! + qdrain_rate = qdrain_uniform + ! + endif + ! + else + ! + ! Allocate minimal arrays for OpenACC compatibility + ! + allocate(qdrain_rate(1)) + qdrain_rate = 0.0 + ! + endif + ! + end subroutine + + subroutine initialize_storage_volume() ! use sfincs_data diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 index bce4af73f..7c6496155 100644 --- a/source/src/sfincs_infiltration.f90 +++ b/source/src/sfincs_infiltration.f90 @@ -18,8 +18,8 @@ subroutine initialize_infiltration() ! character*256 :: varname ! - character(len=3), parameter :: allowed_types(5) = & - ['c2d', 'cna', 'cnb', 'gai', 'hor'] + character(len=3), parameter :: allowed_types(6) = & + ['c2d', 'cna', 'cnb', 'gai', 'hor', 'bkt'] logical :: inftype_exists ! @@ -32,20 +32,22 @@ subroutine initialize_infiltration() infiltration = .false. netcdf_infiltration = .false. ! - ! Four options for infiltration: + ! Seven infiltration flavors (inftype): ! - ! 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) - ! 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 + ! 1) 'con' - Spatially-uniform constant infiltration + ! Requires: qinf (mm/hr in sfincs.inp) + ! 2) 'c2d' - Spatially-varying constant infiltration + ! Requires: qinffile or infiltrationfile + ! 3) 'cna' - SCS Curve Number (old, no recovery) + ! Requires: scsfile or infiltrationfile + ! 4) 'cnb' - SCS Curve Number (new, with recovery) + ! Requires: sefffile or infiltrationfile + ! 5) 'gai' - Green-Ampt infiltration + ! Requires: psifile or infiltrationfile + ! 6) 'hor' - Modified Horton equation + ! Requires: f0file or infiltrationfile + ! 7) 'bkt' - Bucket model (linear reservoir, HBV/wflow style) + ! Requires: bucketfile (netcdf with bucket_smax and bucket_k) ! ! cumprcp and cuminf are stored in the netcdf output if store_cumulative_precipitation == .true. which is the default ! @@ -83,7 +85,7 @@ subroutine initialize_infiltration() ! else ! - write(logstr,*)'Error : infiltration input type ',trim(inftype),' is not part of supported types c2d cna cnb gai hor !' + write(logstr,*)'Error : infiltration input type ',trim(inftype),' is not part of supported types c2d cna cnb gai hor bkt !' call stop_sfincs(trim(logstr), 1) ! end if @@ -124,7 +126,7 @@ subroutine initialize_infiltration() inftype = 'gai' infiltration = .true. ! - elseif (f0file /= 'none') then + elseif (f0file /= 'none') then ! ! The Horton Equation model for infiltration ! @@ -132,6 +134,13 @@ subroutine initialize_infiltration() infiltration = .true. store_meteo = .true. ! + elseif (bucketfile /= 'none') then + ! + ! Bucket model (linear reservoir) + ! + inftype = 'bkt' + infiltration = .true. + ! endif ! ! 2) We need cumprcp and cuminf @@ -498,7 +507,7 @@ subroutine initialize_infiltration() ! Allocate support variables: ! allocate(qinffield(np)) - qinffield(nm) = 0.0 + qinffield = 0.0 ! elseif (inftype == 'hor') then ! @@ -591,6 +600,14 @@ subroutine initialize_infiltration() allocate(rain_T1(np)) rain_T1 = 0.0 ! + elseif (inftype == 'bkt') then + ! + ! Bucket model (linear reservoir) - mimics hydrology models like wflow/HBV + ! + call write_log('Info : turning on process infiltration (via bucket model)', 0) + ! + call initialize_bucket_model() + ! endif ! else @@ -634,7 +651,7 @@ subroutine update_infiltration_map(dt, tloop) !$omp parallel & !$omp private ( nm ) !$omp do - !$acc parallel present( qinfmap, qinffield, z_volume, zs, zb, netprcp, cuminf ) + !$acc parallel present( qinfmap, qinffield, z_volume, zs, zb, cuminf ) !$acc loop independent gang vector do nm = 1, np ! @@ -656,17 +673,13 @@ subroutine update_infiltration_map(dt, tloop) ! endif ! - ! Compute nett precip - ! - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! if (store_cumulative_precipitation) then ! ! Compute cumulative infiltration ! cuminf(nm) = cuminf(nm) + qinfmap(nm) * dt ! - endif + endif ! enddo !$omp end do @@ -680,7 +693,7 @@ subroutine update_infiltration_map(dt, tloop) !$omp parallel & !$omp private ( Qq,I,nm ) !$omp do - !$acc parallel present( qinfmap, qinffield, prcp, netprcp, cumprcp, cuminf ) + !$acc parallel present( qinfmap, qinffield, prcp, cumprcp, cuminf ) !$acc loop independent gang vector do nm = 1, np ! @@ -702,10 +715,6 @@ subroutine update_infiltration_map(dt, tloop) ! endif ! - ! Compute nett precip - ! - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! if (store_cumulative_precipitation) then ! ! Compute cumulative infiltration @@ -726,7 +735,7 @@ subroutine update_infiltration_map(dt, tloop) !$omp parallel & !$omp private ( Qq,I,nm ) !$omp do - !$acc parallel present( qinfmap, prcp, netprcp, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, qinffield, inf_kr ) + !$acc parallel present( qinfmap, prcp, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, qinffield, inf_kr ) !$acc loop independent gang vector do nm = 1, np ! @@ -808,10 +817,6 @@ subroutine update_infiltration_map(dt, tloop) ! endif ! - ! Compute nett precip - ! - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! if (store_cumulative_precipitation) then ! ! Compute cumulative infiltration @@ -822,7 +827,7 @@ subroutine update_infiltration_map(dt, tloop) ! enddo !$omp end do - !$omp end parallel + !$omp end parallel !$acc end parallel ! elseif (inftype == 'gai') then @@ -832,7 +837,7 @@ subroutine update_infiltration_map(dt, tloop) !$omp parallel & !$omp private ( nm ) !$omp do - !$acc parallel present( qinfmap, prcp, netprcp, cuminf, rain_T1, & + !$acc parallel present( qinfmap, prcp, cuminf, rain_T1, & !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr ) !$acc loop independent gang vector do nm = 1, np @@ -853,8 +858,18 @@ subroutine update_infiltration_map(dt, tloop) ! ! Larger amounts of rainfall - Equation 4-27 from SWMM manual ! - qinfmap(nm) = (ksfield(nm) * (1.0 + (GA_head(np) * GA_sigma(np)) / GA_F(nm))) - qinfmap(nm) = max(min(qinfmap(nm), prcp(nm)), 0.0) ! never more than rainfall and and never negative + if (GA_F(nm) < 1.0e-10) then + ! + ! No cumulative infiltration yet (first timestep) - all rainfall infiltrates + ! + qinfmap(nm) = prcp(nm) + ! + else + ! + qinfmap(nm) = (ksfield(nm) * (1.0 + (GA_head(nm) * GA_sigma(nm)) / GA_F(nm))) + qinfmap(nm) = max(min(qinfmap(nm), prcp(nm)), 0.0) ! never more than rainfall and never negative + ! + endif ! endif ! @@ -891,11 +906,6 @@ subroutine update_infiltration_map(dt, tloop) endif endif ! - ! Compute nett precip - ! - !qinffield(nm) = qinfmap(nm) ! Really ? Why ? - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! if (store_cumulative_precipitation) then ! ! Compute cumulative infiltration @@ -906,7 +916,7 @@ subroutine update_infiltration_map(dt, tloop) ! enddo !$omp end do - !$omp end parallel + !$omp end parallel !$acc end parallel ! elseif (inftype == 'hor') then @@ -916,7 +926,7 @@ subroutine update_infiltration_map(dt, tloop) !$omp parallel & !$omp private ( nm, Qq, I, a, hh_local ) !$omp do - !$acc parallel present( qinfmap, prcp, netprcp, cuminf, cell_area_m2, cell_area, z_flags_iref, z_volume, zs, zb, rain_T1, & + !$acc parallel present( qinfmap, prcp, cuminf, cell_area_m2, cell_area, z_flags_iref, z_volume, zs, zb, rain_T1, & !$acc horton_kd, horton_fc, horton_f0 ) !$acc loop independent gang vector do nm = 1, np @@ -1007,10 +1017,6 @@ subroutine update_infiltration_map(dt, tloop) ! endif ! - ! Compute nett precip - ! - netprcp(nm) = netprcp(nm) - qinfmap(nm) - ! if (store_cumulative_precipitation) then ! ! Compute cumulative infiltration @@ -1021,14 +1027,164 @@ subroutine update_infiltration_map(dt, tloop) ! enddo !$omp end do - !$omp end parallel + !$omp end parallel !$acc end parallel ! + elseif (inftype == 'bkt') then + ! + ! Bucket model (linear reservoir) + ! + call compute_bucket_drainage(dt) + ! endif ! call system_clock(count1, count_rate, count_max) tloop = tloop + 1.0 * (count1 - count0) / count_rate ! - end subroutine + end subroutine + + + subroutine initialize_bucket_model() + ! + use sfincs_data + use sfincs_ncinput + ! + implicit none + ! + integer :: nchar + logical :: ok + character*256 :: varname + ! + if (bucketfile /= 'none') then + ! + use_bucket_model = .true. + ! + write(logstr,'(a)')'Info : turning on bucket model (linear reservoir)' + call write_log(logstr, 0) + ! + allocate(bucket_capacity(np)) + allocate(bucket_k(np)) + allocate(bucket_volume(np)) + allocate(bucket_drain_rate(np)) + ! + bucket_capacity = 0.0 + bucket_k = 0.0 + bucket_volume = 0.0 + bucket_drain_rate = 0.0 + ! + nchar = len_trim(bucketfile) + ok = check_file_exists(bucketfile, 'Bucket model file', .true.) + ! + if (bucketfile(nchar - 1 : nchar) == 'nc') then + ! + ! Read bucket capacity (S_max) in mm, convert to m + ! + varname = 'bucket_smax' + call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_capacity) + bucket_capacity = bucket_capacity / 1000.0 ! mm to m + ! + ! Read drainage coefficient (k) in 1/hr, convert to 1/s + ! + varname = 'bucket_k' + call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_k) + bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s + ! + else + ! + ! Read from binary files + ! + open(unit = 500, file = trim(bucketfile), form = 'unformatted', access = 'stream') + read(500)bucket_capacity + close(500) + bucket_capacity = bucket_capacity / 1000.0 ! mm to m + ! + ! For binary input, k needs a separate file - not supported yet + ! Default k = 0.1/hr + ! + bucket_k = 0.1 / 3600.0 + ! + endif + ! + write(logstr,'(a,f10.4,a)')'Info : bucket max capacity = ', maxval(bucket_capacity) * 1000.0, ' mm' + call write_log(logstr, 0) + ! + else + ! + ! Allocate minimal arrays for OpenACC compatibility + ! + allocate(bucket_capacity(1)) + allocate(bucket_k(1)) + allocate(bucket_volume(1)) + allocate(bucket_drain_rate(1)) + bucket_capacity = 0.0 + bucket_k = 0.0 + bucket_volume = 0.0 + bucket_drain_rate = 0.0 + ! + endif + ! + end subroutine + + + subroutine compute_bucket_drainage(dt) + ! + ! Bucket model: finite capacity reservoir with linear drainage (HBV/wflow style) + ! Recovery is inherent - bucket drains during dry periods via Q=k*S, restoring capacity + ! + ! Literature: Linear reservoir (Nash, 1957), HBV soil moisture bucket (Bergstrom, 1995) + ! + use sfincs_data + ! + implicit none + ! + real*4 :: dt + integer :: nm + real*4 :: exp_factor + real*4 :: drain_vol + real*4 :: available_water + real*4 :: available_cap + real*4 :: actual_inflow + ! + !$omp parallel do private(nm, exp_factor, drain_vol, available_water, available_cap, actual_inflow) + do nm = 1, np + ! + if (kcs(nm) == 1 .and. bucket_capacity(nm) > 0.0) then + ! + ! Step 1: Drain current storage (analytical linear reservoir solution) + ! S(t+dt) = S(t) * exp(-k*dt), drainage = S(t) - S(t+dt) + ! + exp_factor = exp(-bucket_k(nm) * dt) + drain_vol = bucket_volume(nm) * (1.0 - exp_factor) + bucket_volume(nm) = bucket_volume(nm) * exp_factor + ! + ! Step 2: Fill bucket from available rainfall + ! + available_water = max(prcp(nm), 0.0) * dt ! m of water available + available_cap = bucket_capacity(nm) - bucket_volume(nm) + actual_inflow = min(available_water, available_cap) + bucket_volume(nm) = bucket_volume(nm) + actual_inflow + ! + ! Step 3: Set qinfmap = what entered the bucket (removed from surface) + ! This is used by continuity as the infiltration loss term + ! + qinfmap(nm) = actual_inflow / dt + bucket_drain_rate(nm) = actual_inflow / dt + ! + if (store_cumulative_precipitation) then + cuminf(nm) = cuminf(nm) + qinfmap(nm) * dt + endif + ! + else + ! + qinfmap(nm) = 0.0 + bucket_drain_rate(nm) = 0.0 + ! + endif + ! + enddo + !$omp end parallel do + ! + end subroutine + end module diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index 31b22daa9..34163b0bd 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -209,8 +209,12 @@ subroutine read_sfincs_input() call read_char_input(500,'netampfile',netampfile,'none') call read_char_input(500,'netspwfile',netspwfile,'none') ! - call read_char_input(500,'infiltrationfile',infiltrationfile,'none') + ! Infiltration and losses + call read_char_input(500,'infiltrationfile',infiltrationfile,'none') call read_char_input(500,'infiltrationtype',inftype,'none') + call read_char_input(500,'bucketfile',bucketfile,'none') ! bucket model (infiltration flavor 'bkt') + call read_real_input(500,'qdrain',qdrain_uniform,0.0) ! drainage mimic (mm/hr) + call read_char_input(500,'drainagefile',drainagefile,'none') ! spatially-varying drainage rates ! ! Output call read_char_input(500,'obsfile',obsfile,'none') @@ -313,6 +317,7 @@ subroutine read_sfincs_input() gn2 = 9.81*0.02*0.02 ! Only to be used in subgrid ! qinf = qinf/(3600*1000) + qdrain_uniform = qdrain_uniform/(3600*1000) ! Convert mm/hr to m/s ! rotation = rotation*pi/180 cosrot = cos(rotation) diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index 11318103e..f40d12bf6 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -500,26 +500,12 @@ function sfincs_update(dtrange) result(ierr) ! call update_meteo_forcing(t, dt, tloopwnd2) ! - ! Update infiltration - ! - if (infiltration) then - ! - ! Compute infiltration rates - ! - call update_infiltration_map(dt, tloopinf) - ! - endif - ! endif ! ! Update boundary conditions ! call update_boundaries(t, dt, tloopbnd) ! - ! Update discharges - ! - call update_discharges(t, dt, tloopsrc) - ! if (snapwave .and. update_waves) then ! call timer(t3) @@ -570,9 +556,9 @@ function sfincs_update(dtrange) result(ierr) ! endif ! - ! Update water levels + ! Update continuity (discharges, infiltration, drainage, water levels) ! - call compute_water_levels(t, dt, tloopcont) + call update_continuity(t, dt, tloopsrc, tloopinf, tloopcont) ! ! OUTPUT ! diff --git a/source/src/sfincs_meteo.f90 b/source/src/sfincs_meteo.f90 index 733ae8467..62175caa3 100644 --- a/source/src/sfincs_meteo.f90 +++ b/source/src/sfincs_meteo.f90 @@ -1214,7 +1214,7 @@ subroutine update_meteo_forcing(t, dt, tloop) !$acc parallel, present( tauwu, tauwv, tauwu0, tauwv0, tauwu1, tauwv1, & !$acc windu, windv, windu0, windv0, windu1, windv1, windmax, & !$acc patm, patm0, patm1, & - !$acc prcp, prcp0, prcp1, cumprcp, netprcp, & + !$acc prcp, prcp0, prcp1, cumprcp, & !$acc zs, zb, z_volume ) !$acc loop gang vector do nm = 1, np @@ -1262,7 +1262,6 @@ subroutine update_meteo_forcing(t, dt, tloop) endif endif ! - netprcp(nm) = prcp(nm) cumprcp(nm) = cumprcp(nm) + prcp(nm) * dt ! endif @@ -1282,40 +1281,38 @@ subroutine update_meteo_forcing(t, dt, tloop) !$omp parallel & !$omp private ( nm ) !$omp do - !$acc parallel, present( tauwu, tauwv, patm, prcp, netprcp, zs, zb, z_volume ) + !$acc parallel, present( tauwu, tauwv, patm, prcp, zs, zb, z_volume ) !$acc loop gang vector do nm = 1, np ! if (wind) then tauwu(nm) = tauwu(nm) * smfac tauwv(nm) = tauwv(nm) * smfac - endif + endif ! if (patmos) then patm(nm) =patm(nm) * smfac + gapres * oneminsmfac - endif + endif ! if (precip) then - ! - netprcp(nm) = netprcp(nm) * smfac - ! - ! Don't allow negative netprcp during spinup (e.g. hardfixing infiltration/evaporation on model when forcing effective rainfall) when there's no water in the cell (same as check for constant infiltration) - ! - if (netprcp(nm) < 0.0) then - ! - ! No effective infiltration if there is no water - ! + ! + prcp(nm) = prcp(nm) * smfac + ! + ! Don't allow negative precip during spinup when there's no water in the cell + ! + if (prcp(nm) < 0.0) then + ! if (subgrid) then if (z_volume(nm) <= 0.0) then - netprcp(nm) = 0.0 + prcp(nm) = 0.0 endif else if (zs(nm) <= zb(nm)) then - netprcp(nm) = 0.0 + prcp(nm) = 0.0 endif - endif - ! - endif + endif + ! + endif endif ! enddo @@ -1465,12 +1462,11 @@ subroutine update_precipitation_from_timeseries(t, dt) !$omp parallel & !$omp private ( nm ) !$omp do - !$acc parallel present( prcp, cumprcp, netprcp ) + !$acc parallel present( prcp, cumprcp ) !$acc loop gang vector do nm = 1, np ! prcp(nm) = ptmp - netprcp(nm) = ptmp ! if (store_cumulative_precipitation) then cumprcp(nm) = cumprcp(nm) + ptmp * dt diff --git a/source/src/sfincs_openacc.f90 b/source/src/sfincs_openacc.f90 index d751d9d04..3d2a9cdcc 100644 --- a/source/src/sfincs_openacc.f90 +++ b/source/src/sfincs_openacc.f90 @@ -28,14 +28,15 @@ subroutine initialize_openacc() !$acc tauwu, tauwv, tauwu0, tauwv0, tauwu1, tauwv1, & !$acc windu, windv, windu0, windv0, windu1, windv1, windmax, & !$acc patm, patm0, patm1, patmb, nmindbnd, & - !$acc prcp, prcp0, prcp1, cumprcp, netprcp, prcp, qext, & + !$acc prcp, prcp0, prcp1, cumprcp, qext, & !$acc dxminv, dxrinv, dyrinv, dxm2inv, dxr2inv, dyr2inv, dxrinvc, dyrinvc, dxm, dxrm, dyrm, cell_area_m2, cell_area, & !$acc gn2uv, fcorio2d, storage_volume, nuvisc, & !$acc cuv_index_uv, cuv_index_uv1, cuv_index_uv2, & !$acc x73, & !$acc gnapp2, & !$acc qinffield, qinfmap, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, & - !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0 ) + !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0, & + !$acc qdrain_rate, bucket_volume, bucket_capacity, bucket_k, bucket_drain_rate ) ! end subroutine ! @@ -57,14 +58,15 @@ subroutine finalize_openacc() !$acc tauwu, tauwv, tauwu0, tauwv0, tauwu1, tauwv1, & !$acc windu, windv, windu0, windv0, windu1, windv1, windmax, & !$acc patm, patm0, patm1, patmb, nmindbnd, & - !$acc prcp, prcp0, prcp1, cumprcp, netprcp, prcp, qext, & + !$acc prcp, prcp0, prcp1, cumprcp, qext, & !$acc dxminv, dxrinv, dyrinv, dxm2inv, dxr2inv, dyr2inv, dxrinvc, dxm, dxrm, dyrm, cell_area_m2, cell_area, & !$acc gn2uv, fcorio2d, storage_volume, nuvisc, & !$acc cuv_index_uv, cuv_index_uv1, cuv_index_uv2, & !$acc x73, & !$acc gnapp2, & !$acc qinffield, qinfmap, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, & - !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0 ) + !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0, & + !$acc qdrain_rate, bucket_volume, bucket_capacity, bucket_k, bucket_drain_rate ) ! end ! From cf6a78d124d3ffae94ef9b8576c6950189e7c7c5 Mon Sep 17 00:00:00 2001 From: Kees Nederhoff Date: Thu, 19 Mar 2026 18:07:13 -0700 Subject: [PATCH 15/19] Included more output options --- source/sfincs/sfincs.log | 40 +++++ source/sfincs/sfincs.vfproj.keesn.user | 8 + .../sfincs_lib/sfincs_lib.vfproj.keesn.user | 8 + source/src/sfincs_continuity.f90 | 4 +- source/src/sfincs_data.f90 | 4 +- source/src/sfincs_domain.f90 | 4 +- source/src/sfincs_infiltration.f90 | 22 ++- source/src/sfincs_lib.f90 | 9 +- source/src/sfincs_ncoutput.F90 | 151 ++++++++++++++---- 9 files changed, 199 insertions(+), 51 deletions(-) create mode 100644 source/sfincs/sfincs.log create mode 100644 source/sfincs/sfincs.vfproj.keesn.user create mode 100644 source/sfincs_lib/sfincs_lib.vfproj.keesn.user diff --git a/source/sfincs/sfincs.log b/source/sfincs/sfincs.log new file mode 100644 index 000000000..2b94f1aae --- /dev/null +++ b/source/sfincs/sfincs.log @@ -0,0 +1,40 @@ + +------------ Welcome to SFINCS ------------ + + @@@@@ @@@@@@@ @@ @@ @@ @@@@ @@@@@ + @@@ @@@ @@@@@@@ @@ @@@ @@ @@@@@@@ @@@ @@@ + @@@ @@ @@ @@@ @@ @@ @@ @@@ + @@@@@ @@@@@@ @@ @@@@@@ @@ @@@@@ + @@@ @@ @@ @@ @@@ @@ @@ @@@ + @@@ @@@ @@ @@ @@ @@ @@@@@@ @@@ @@@ + @@@@@ @@ @@ @@ @ @@@@ @@@@@ + + .............. + ......:@@@@@@@@:...... + ..::::..@@........@@.:::::.. + ..:::::..@@..::..::..@@.::::::.. + .::::::..@@............@@.:::::::. + .::::::..@@..............@@.:::::::. + .::::::::..@@............@@..::::::::. + .:::::::::...@@.@..@@..@.@@..::::::::::. + .:::::::::...:@@@..@@..@@@:..:::::::::.. + ............@@.@@..@@..@@.@@............ + ^^^~~^^~~^^@@..............@@^^^~^^^~~^^ + .::::::::::@@..............@@.:::::::::. + .......:.@@.....@.....@....@@.:....... + .::....@@......@.@@@.@....@@.....::. + .:::~@@.:...:.@@...@@.:.:.@@~::::. + .::~@@@@@@@@@@.....@@@@@@@@@~::. + ..:~~~~~~~:.......:~~~~~~~:.. + ...................... + .............. + +------------------------------------------ + +Build-Revision: $Rev: v2.3.1 mt. Faber+branch-redo-infiltration +Build-Date: $Date: 2026-03-19 + +------ Preparing model simulation -------- + +Reading input file ... +Error : SFINCS input file "sfincs.inp" not found! SFINCS has stopped! diff --git a/source/sfincs/sfincs.vfproj.keesn.user b/source/sfincs/sfincs.vfproj.keesn.user new file mode 100644 index 000000000..6ae2aa06f --- /dev/null +++ b/source/sfincs/sfincs.vfproj.keesn.user @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/source/sfincs_lib/sfincs_lib.vfproj.keesn.user b/source/sfincs_lib/sfincs_lib.vfproj.keesn.user new file mode 100644 index 000000000..818b85e15 --- /dev/null +++ b/source/sfincs_lib/sfincs_lib.vfproj.keesn.user @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90 index bcfc9699a..d50437556 100644 --- a/source/src/sfincs_continuity.f90 +++ b/source/src/sfincs_continuity.f90 @@ -182,7 +182,7 @@ subroutine compute_water_levels_regular(dt,t) ! ! C3. Drainage mimic (-) ! - if (use_drainage_mimic) then + if (drainage) then zs(nm) = zs(nm) - qdrain_rate(nm) * dt endif ! @@ -551,7 +551,7 @@ subroutine compute_water_levels_subgrid(dt,t) ! ! C3. Drainage mimic (-) ! - if (use_drainage_mimic) then + if (drainage) then dzsdt = dzsdt - qdrain_rate(nm) endif ! diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 111c4cc5c..2eb81c5a0 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -410,9 +410,9 @@ module sfincs_data ! real*4, dimension(:), allocatable :: storage_volume ! Storage volume green infra ! - ! Drainage mimic - constant removal rate representing subsurface drainage + ! Drainage - constant removal rate representing subsurface drainage ! - logical :: use_drainage_mimic = .false. + logical :: drainage = .false. real*4 :: qdrain_uniform = 0.0 ! uniform drainage rate (mm/hr input, stored as m/s) real*4, dimension(:), allocatable :: qdrain_rate ! drainage rate per cell (m/s) ! diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 7cfce450b..9b0383092 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -2061,11 +2061,11 @@ subroutine initialize_drainage_mimic() logical :: ok character*256 :: varname ! - ! Check if drainage mimic is enabled + ! Check if drainage is enabled ! if (qdrain_uniform > 0.0 .or. drainagefile /= 'none') then ! - use_drainage_mimic = .true. + drainage = .true. ! allocate(qdrain_rate(np)) ! diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 index 7c6496155..8da1855d4 100644 --- a/source/src/sfincs_infiltration.f90 +++ b/source/src/sfincs_infiltration.f90 @@ -177,22 +177,18 @@ subroutine initialize_infiltration() ! 5) Check whether infiltration input type (orignal vs netcdf) are correctly matched to grid type (regular vs quadtree) ! if (infiltration .and. inftype /= 'con') then !constant uniform works for both options - ! - if (netcdf_infiltration) then - ! - if (use_quadtree .eqv. .false.) then - ! - call stop_sfincs('Error ! Netcdf infiltration input format can only be specified for quadtree mesh model !', 1) + ! + ! Netcdf infiltration works for both regular and quadtree grids + ! (regular grids populate quadtree_nr_points and index_sfincs_in_quadtree + ! via make_quadtree_from_indices) + ! + if (.not. netcdf_infiltration) then + ! + if (use_quadtree .eqv. .true.) then ! + call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) endif ! - else ! Original - ! - if (use_quadtree .eqv. .true.) then - ! - call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) - endif - ! endif ! endif diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index f40d12bf6..39e193ec2 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -92,8 +92,8 @@ function sfincs_initialize() result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = "$Rev: v2.3.1 mt. Faber+branch-272" - build_date = "$Date: 2026-03-06" + build_revision = "$Rev: v2.3.1 mt. Faber+branch-redo-infiltration" + build_date = "$Date: 2026-03-19" ! call write_log('', 1) call write_log('------------ Welcome to SFINCS ------------', 1) @@ -225,6 +225,11 @@ function sfincs_initialize() result(ierr) else call write_log('Infiltration : no', 1) endif + if (drainage) then + call write_log('Drainage : yes', 1) + else + call write_log('Drainage : no', 1) + endif if (snapwave) then call write_log('SnapWave : yes', 1) else diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index 3fcfbe9d6..6a9ecdf3d 100644 --- a/source/src/sfincs_ncoutput.F90 +++ b/source/src/sfincs_ncoutput.F90 @@ -214,12 +214,16 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'suction head at the wetting front - Green and Ampt')) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'm')) elseif (inftype == 'hor') then - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'f0')) - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'initial infiltration rate - Horton')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'f0')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'initial infiltration rate - Horton')) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'm')) + elseif (inftype == 'bkt') then + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'bucket_capacity')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'maximum bucket storage capacity')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'mm')) else - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'qinf')) - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'infiltration rate - constant in time')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'qinf')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'infiltration rate - constant in time')) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'mm h-1')) endif endif @@ -370,14 +374,25 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'coordinates', 'corner_x corner_y')) endif ! - ! Store current infiltration (only for Horton) + ! Store current infiltration capacity (only for Horton) ! if (inftype == 'hor') then - NF90(nf90_def_var(map_file%ncid, 'f', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%Seff_varid)) ! time-varying sigma - NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_def_var(map_file%ncid, 'f', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%Seff_varid)) ! time-varying f + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'units', 'mm h-1')) - NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'standard_name', 'sigma')) - NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'long_name', 'current infiltration capacity')) + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'standard_name', 'f')) + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'long_name', 'current infiltration capacity')) + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'coordinates', 'corner_x corner_y')) + endif + ! + ! Store current bucket storage (only for Bucket model) + ! + if (inftype == 'bkt') then + NF90(nf90_def_var(map_file%ncid, 'bucket_volume', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%Seff_varid)) ! time-varying bucket volume + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'units', 'm')) + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'standard_name', 'bucket_volume')) + NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'long_name', 'current bucket storage')) NF90(nf90_put_att(map_file%ncid, map_file%Seff_varid, 'coordinates', 'corner_x corner_y')) endif ! @@ -801,9 +816,9 @@ subroutine ncoutput_regular_map_init() ! ! Write infiltration map ! - if (infiltration) then + if (infiltration .and. allocated(qinffield)) then ! - zsg = FILL_VALUE + zsg = FILL_VALUE ! do nm = 1, np ! @@ -819,13 +834,29 @@ subroutine ncoutput_regular_map_init() zsg(m, n) = qinffield(nm) ! endif - ! + ! enddo ! NF90(nf90_put_var(map_file%ncid, map_file%qinf_varid, zsg, (/1, 1/))) ! write infiltration map ! endif ! + ! Write bucket capacity map (static) + ! + if (inftype == 'bkt' .and. allocated(bucket_capacity)) then + ! + zsg = FILL_VALUE + ! + do nm = 1, np + n = z_index_z_n(nm) + m = z_index_z_m(nm) + zsg(m, n) = bucket_capacity(nm) * 1000.0 ! m to mm + enddo + ! + NF90(nf90_put_var(map_file%ncid, map_file%qinf_varid, zsg, (/1, 1/))) ! write bucket capacity map + ! + endif + ! ! write away intermediate data ! NF90(nf90_sync(map_file%ncid)) !write away intermediate data @@ -1383,12 +1414,16 @@ subroutine ncoutput_quadtree_map_init() NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'suction head at the wetting front - Green and Ampt')) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'm')) elseif (inftype == 'hor') then - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'f0')) - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'initial infiltration rate - Horton')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'f0')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'initial infiltration rate - Horton')) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'mm h-1')) + elseif (inftype == 'bkt') then + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'bucket_capacity')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'maximum bucket storage capacity')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'mm')) else - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'qinf')) - NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'infiltration rate - constant in time')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'qinf')) + NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'long_name', 'infiltration rate - constant in time')) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'units', 'mm h-1')) endif endif @@ -1530,7 +1565,7 @@ subroutine ncoutput_quadtree_map_init() ! vtmp = FILL_VALUE ! - if (infiltration) then + if (infiltration .and. allocated(qinffield)) then ! if (inftype == 'con' .or. inftype == 'c2d') then do nmq = 1, quadtree_nr_points @@ -1538,20 +1573,37 @@ subroutine ncoutput_quadtree_map_init() if (nm>0) then vtmp(nmq) = qinffield(nm) * 3600 * 1000 endif - enddo + enddo else do nmq = 1, quadtree_nr_points nm = index_sfincs_in_quadtree(nmq) if (nm>0) then vtmp(nmq) = qinffield(nm) endif - enddo + enddo endif ! NF90(nf90_put_var(map_file%ncid, map_file%qinf_varid, vtmp)) ! write infiltration map ! endif ! + ! Write bucket capacity map (static) + ! + if (inftype == 'bkt' .and. allocated(bucket_capacity)) then + ! + vtmp = FILL_VALUE + ! + do nmq = 1, quadtree_nr_points + nm = index_sfincs_in_quadtree(nmq) + if (nm>0) then + vtmp(nmq) = bucket_capacity(nm) * 1000.0 ! m to mm + endif + enddo + ! + NF90(nf90_put_var(map_file%ncid, map_file%qinf_varid, vtmp)) ! write bucket capacity map + ! + endif + ! ! write away intermediate data ! NF90(nf90_sync(map_file%ncid)) !write away intermediate data @@ -1810,9 +1862,29 @@ subroutine ncoutput_his_init() ! if (inftype == 'gai') then NF90(nf90_def_var(his_file%ncid, 'point_S', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%S_varid)) ! time-varying S - NF90(nf90_put_att(his_file%ncid, his_file%S_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'units', 'm')) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'long_name', 'maximum soil moisture deficit')) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'coordinates', 'station_id station_name point_x point_y')) + endif + ! + ! More output for Horton method + ! + if (inftype == 'hor') then + NF90(nf90_def_var(his_file%ncid, 'point_S', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%S_varid)) ! time-varying f + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'units', 'mm hr-1')) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'long_name', 'current infiltration capacity')) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'coordinates', 'station_id station_name point_x point_y')) + endif + ! + ! More output for Bucket model + ! + if (inftype == 'bkt') then + NF90(nf90_def_var(his_file%ncid, 'point_S', NF90_FLOAT, (/his_file%points_dimid, his_file%time_dimid/), his_file%S_varid)) ! time-varying bucket volume + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'units', 'm')) - NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'long_name', 'maximum soil moisture deficit')) + NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'long_name', 'current bucket storage')) NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'coordinates', 'station_id station_name point_x point_y')) endif ! @@ -1987,10 +2059,10 @@ subroutine ncoutput_his_init() if (ndrn>0) then ! NF90(nf90_def_var(his_file%ncid, 'drainage_discharge', NF90_FLOAT, (/his_file%drain_dimid, his_file%time_dimid/), his_file%drain_varid)) ! time-varying discharge through drainage structure - NF90(nf90_put_att(his_file%ncid, his_file%discharge_varid, '_FillValue', FILL_VALUE)) - NF90(nf90_put_att(his_file%ncid, his_file%discharge_varid, 'units', 'm3 s-1')) - NF90(nf90_put_att(his_file%ncid, his_file%discharge_varid, 'long_name', 'discharge through drainage structure')) - NF90(nf90_put_att(his_file%ncid, his_file%discharge_varid, 'coordinates', 'drainage_name')) + NF90(nf90_put_att(his_file%ncid, his_file%drain_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(his_file%ncid, his_file%drain_varid, 'units', 'm3 s-1')) + NF90(nf90_put_att(his_file%ncid, his_file%drain_varid, 'long_name', 'discharge through drainage structure')) + NF90(nf90_put_att(his_file%ncid, his_file%drain_varid, 'coordinates', 'drainage_name')) ! endif ! @@ -2302,13 +2374,30 @@ subroutine ncoutput_update_regular_map(t,ntmapout) ! n = z_index_z_n(nm) m = z_index_z_m(nm) - ! + ! zsg(m, n) = qinfmap(nm) ! enddo ! NF90(nf90_put_var(map_file%ncid, map_file%Seff_varid, zsg, (/1, 1, ntmapout/))) ! + elseif (inftype == 'bkt') then + ! + ! Store current bucket volume + ! + zsg = FILL_VALUE + ! + do nm = 1, np + ! + n = z_index_z_n(nm) + m = z_index_z_m(nm) + ! + zsg(m, n) = bucket_volume(nm) + ! + enddo + ! + NF90(nf90_put_var(map_file%ncid, map_file%Seff_varid, zsg, (/1, 1, ntmapout/))) + ! endif ! if (store_meteo) then @@ -3047,8 +3136,12 @@ subroutine ncoutput_update_his(t,nthisout) tS_effective(iobs) = scs_Se(nm) elseif (inftype == 'gai') then tS_effective(iobs) = GA_sigma(nm) + elseif (inftype == 'hor') then + tS_effective(iobs) = qinfmap(nm)*3.6e3*1.0e3 ! current f in mm/hr + elseif (inftype == 'bkt') then + tS_effective(iobs) = bucket_volume(nm) ! current bucket storage in m endif - endif + endif ! if (store_meteo) then ! @@ -3127,9 +3220,7 @@ subroutine ncoutput_update_his(t,nthisout) ! NF90(nf90_put_var(his_file%ncid, his_file%qinf_varid, tqinf, (/1, nthisout/))) ! write qinf ! - if (inftype == 'cnb') then - NF90(nf90_put_var(his_file%ncid, his_file%S_varid, tS_effective, (/1, nthisout/))) ! write S - elseif (inftype == 'gai') then + if (inftype == 'cnb' .or. inftype == 'gai' .or. inftype == 'hor' .or. inftype == 'bkt') then NF90(nf90_put_var(his_file%ncid, his_file%S_varid, tS_effective, (/1, nthisout/))) ! write S endif ! From 3300914c1b82cb6c997132a916f707bfe890deb7 Mon Sep 17 00:00:00 2001 From: Leynse Date: Tue, 24 Mar 2026 15:18:36 +0100 Subject: [PATCH 16/19] - clean repo --- .gitignore | 8 +++++ source/Dockerfile.gpu | 31 ----------------- source/Dockerfile.gpu.25.5.ccall | 32 ------------------ source/Dockerfile.gpu.test | 15 --------- source/Dockerfile.gpu.update01 | 32 ------------------ source/Dockerfile.xpu | 21 ------------ source/Singularityfile-gpu.def | 58 -------------------------------- source/build_nvfortran_gpu.sh | 35 ------------------- source/build_nvfortran_gpu_h7.sh | 21 ------------ 9 files changed, 8 insertions(+), 245 deletions(-) delete mode 100644 source/Dockerfile.gpu delete mode 100644 source/Dockerfile.gpu.25.5.ccall delete mode 100644 source/Dockerfile.gpu.test delete mode 100644 source/Dockerfile.gpu.update01 delete mode 100644 source/Dockerfile.xpu delete mode 100644 source/Singularityfile-gpu.def delete mode 100644 source/build_nvfortran_gpu.sh delete mode 100644 source/build_nvfortran_gpu_h7.sh diff --git a/.gitignore b/.gitignore index 22752d026..6156e1289 100644 --- a/.gitignore +++ b/.gitignore @@ -65,3 +65,11 @@ source/third_party_open/netcdf/x64 source/sfincs/sfincs.opt.yaml /source/sfincs_lib/*.yaml /source/third_party_open/netcdf/netcdf-fortran-4.6.1/Debug +/source/build_nvfortran_gpu_h7.sh +/source/build_nvfortran_gpu.sh +/source/Singularityfile-gpu.def +/source/Dockerfile.xpu +/source/Dockerfile.gpu.update01 +/source/Dockerfile.gpu.test +/source/Dockerfile.gpu.25.5.ccall +/source/Dockerfile.gpu diff --git a/source/Dockerfile.gpu b/source/Dockerfile.gpu deleted file mode 100644 index 24da778bc..000000000 --- a/source/Dockerfile.gpu +++ /dev/null @@ -1,31 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:20.11-devel-cuda_multi-ubuntu20.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y dos2unix -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/20.11/REDIST/compilers && tar -czvhf /root/nvidia-20.11.tar.gz * && cd - -COPY . /usr/src/sfincs -WORKDIR /usr/src/sfincs - -RUN find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -RUN find . -name \*.f90|xargs dos2unix -RUN find . -name \*.F90|xargs dos2unix -RUN find . -name \*.am|xargs dos2unix -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=8" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -FROM nvidia/cuda:11.2.2-runtime-ubuntu20.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev tzdata -# copy nvidia SDK (REDIST folder contains symlinks) -COPY --from=0 /root/nvidia-20.11.tar.gz . -# copy software -COPY --from=0 /usr/local /usr/local -RUN mkdir -p /opt/nvidia/compilers -RUN tar -xzf nvidia-20.11.tar.gz -C /opt/nvidia/compilers -# add the REDIST libs to the environment -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] - diff --git a/source/Dockerfile.gpu.25.5.ccall b/source/Dockerfile.gpu.25.5.ccall deleted file mode 100644 index fb8e25234..000000000 --- a/source/Dockerfile.gpu.25.5.ccall +++ /dev/null @@ -1,32 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:25.5-devel-cuda_multi-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y dos2unix -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/2025/REDIST/compilers && tar -czvhf /root/nvidia-25.5.tar.gz * && cd - -COPY . /usr/src/sfincs -WORKDIR /usr/src/sfincs - -RUN find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -RUN find . -name \*.f90|xargs dos2unix -RUN find . -name \*.F90|xargs dos2unix -RUN find . -name \*.am|xargs dos2unix -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=8" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install -#RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -FROM nvidia/cuda:12.5.0-runtime-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev tzdata -# copy nvidia SDK (REDIST folder contains symlinks) -COPY --from=0 /root/nvidia-25.5.tar.gz . -# copy software -COPY --from=0 /usr/local /usr/local -RUN mkdir -p /opt/nvidia/compilers -RUN tar -xzf nvidia-25.5.tar.gz -C /opt/nvidia/compilers -# add the REDIST libs to the environment -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] - diff --git a/source/Dockerfile.gpu.test b/source/Dockerfile.gpu.test deleted file mode 100644 index 87870c2ac..000000000 --- a/source/Dockerfile.gpu.test +++ /dev/null @@ -1,15 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:20.11-devel-cuda_multi-ubuntu20.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/20.11/REDIST/compilers && tar -czvhf /root/nvidia-20.11.tar.gz * && cd - -WORKDIR /usr/src/sfincs -COPY . /usr/src/sfincs -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -# create a debug version -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -g -O0 -gpu=ccall" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] diff --git a/source/Dockerfile.gpu.update01 b/source/Dockerfile.gpu.update01 deleted file mode 100644 index bd54b6412..000000000 --- a/source/Dockerfile.gpu.update01 +++ /dev/null @@ -1,32 +0,0 @@ -FROM nvcr.io/nvidia/nvhpc:24.3-devel-cuda_multi-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y dos2unix -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -RUN cd /opt/nvidia/hpc_sdk/Linux_x86_64/2024/REDIST/compilers && tar -czvhf /root/nvidia-24.5.tar.gz * && cd - -COPY . /usr/src/sfincs -WORKDIR /usr/src/sfincs - -RUN find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -RUN find . -name \*.f90|xargs dos2unix -RUN find . -name \*.F90|xargs dos2unix -RUN find . -name \*.am|xargs dos2unix -# ccall -> Generate code for all compute capabilities supported by this platform and by the selected or default CUDA Toolkit. -RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=8" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install -#RUN autoreconf -vif && FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall" FC=nvfortran ./configure --disable-shared --disable-openmp && make && make install - -FROM nvidia/cuda:12.5.0-runtime-ubuntu22.04 -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev tzdata -# copy nvidia SDK (REDIST folder contains symlinks) -COPY --from=0 /root/nvidia-24.5.tar.gz . -# copy software -COPY --from=0 /usr/local /usr/local -RUN mkdir -p /opt/nvidia/compilers -RUN tar -xzf nvidia-24.5.tar.gz -C /opt/nvidia/compilers -# add the REDIST libs to the environment -ENV LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/opt/nvidia/compilers/lib -VOLUME /data -WORKDIR /data -CMD ["sfincs"] - diff --git a/source/Dockerfile.xpu b/source/Dockerfile.xpu deleted file mode 100644 index c311d2e67..000000000 --- a/source/Dockerfile.xpu +++ /dev/null @@ -1,21 +0,0 @@ -FROM ubuntu:latest -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -# copy the compilers -COPY aocc-compiler-4.0.0_1_amd64.deb . -RUN apt install -y ./aocc-compiler-4.0.0_1_amd64.deb -WORKDIR /usr/src/sfincs -COPY . /usr/src/sfincs -# Compile with the AMD fortran compiler -ENV PATH="/opt/AMD/aocc-compiler-4.0.0/bin:${PATH}" -RUN autoreconf -vif && FCFLAGS="-O3" FC=flang ./configure --disable-shared --disable-openmp && make && make install - -# Copy the files to empty docker container -FROM ubuntu:latest -ENV DEBIAN_FRONTEND=noninteractive -RUN apt update && apt install -y libnetcdf-dev tzdata -# copy nvidia SDK (REDIST folder contains symlinks) -COPY --from=0 /usr/local /usr/local -VOLUME /data -WORKDIR /data -CMD ["sfincs"] diff --git a/source/Singularityfile-gpu.def b/source/Singularityfile-gpu.def deleted file mode 100644 index 69ab0a258..000000000 --- a/source/Singularityfile-gpu.def +++ /dev/null @@ -1,58 +0,0 @@ -BootStrap: library -From: library://library/default/ubuntu:jammy -Stage: build - -# -# The source files of sfincs -# - -%files - . /usr/src/sfincs - -# -# Compiler flags -# -fallow-argument-mismatch needed for https://github.com/Unidata/netcdf-fortran/issues/212 -# - -%environment - FCFLAGS="-fopenmp -O3 -fallow-argument-mismatch -w" - FFLAGS="-fopenmp -O3 -fallow-argument-mismatch -w" - -# -# Installing compilers -# Compiling sfincs -# - -%post - apt clean && apt autoclean && apt -y update --fix-missing && apt -y upgrade -y && apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config gfortran gdb m4 - -# Install software into /usr/local by default - - cd /usr/src/sfincs - export CONFIG_SHELL=/bin/bash - autoreconf -vif - ./autogen.sh - ./configure FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=999" FC=nvfortran --disable-shared --disable-openmp --program-suffix="_async" - - make clean - make - make install - -# Install binary in a much smaller image - -BootStrap: library -From: library://library/default/ubuntu:jammy -Stage: final - -# install binary from stage one -%files from build - /usr/local/bin/sfincs /usr/local/bin/sfincs - -%files - run-sfincs.sh /usr/local/bin/run-sfincs.sh - -%post - apt clean && apt autoclean && apt update --fix-missing && apt upgrade -y - apt install -y libnetcdf19 libgfortran5 libgomp1 - apt clean && apt autoclean - chmod +x /usr/local/bin/run-sfincs.sh diff --git a/source/build_nvfortran_gpu.sh b/source/build_nvfortran_gpu.sh deleted file mode 100644 index a38ae8a7e..000000000 --- a/source/build_nvfortran_gpu.sh +++ /dev/null @@ -1,35 +0,0 @@ -#!/bin/bash - -LOGFILE=build_$(date +%Y%m%d_%H%M%S).log -exec > >(tee "$LOGFILE") 2>&1 - -echo "Starting build at $(date)" -echo "Saving log to $LOGFILE" - -find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -find . -name \*.f90|xargs dos2unix -find . -name \*.F90|xargs dos2unix -find . -name \*.am|xargs dos2unix -find . -name \*.sh|xargs dos2unix - -MANPATH=$MANPATH:/opt/nvidia/hpc_sdk/Linux_x86_64/24.5/compilers/man; export MANPATH -PATH=/opt/nvidia/hpc_sdk/Linux_x86_64/24.5/compilers/bin:$PATH; export PATH - -LD_LIBRARY_PATH=/usr/lib/wsl/lib:$LD_LIBRARY_PATH; export LD_LIBRARY_PATH - -apt install -y libnetcdf-dev build-essential autoconf automake libtool pkg-config tzdata -export CONFIG_SHELL=/bin/bash - -autoreconf -vif - -./autogen.sh - -./configure FCFLAGS="-acc=gpu -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=999" FC=nvfortran --disable-shared --disable-openmp --prefix=/usr/local/bin/sfincs/nvfortran_gpu_ccall - -make clean - -make - -make install - -echo "Build finished at $(date)" diff --git a/source/build_nvfortran_gpu_h7.sh b/source/build_nvfortran_gpu_h7.sh deleted file mode 100644 index c50fb36a6..000000000 --- a/source/build_nvfortran_gpu_h7.sh +++ /dev/null @@ -1,21 +0,0 @@ -#!/bin/sh - -module load nvidia/nvhpc/24.1 -module load netcdf - -export LDFLAGS="-L${NETCDF_C_LIBRARY}" - -find . -name \*.m4|xargs dos2unix && find . -name \*.ac|xargs dos2unix && find . -name \*.am|xargs dos2unix -find . -name \*.f90|xargs dos2unix -find . -name \*.F90|xargs dos2unix -find . -name \*.am|xargs dos2unix - -./autogen.sh - -./configure FCFLAGS="-acc -Minfo=accel -fast -O3 -gpu=ccall -DSIZEOF_PTRDIFF_T=999" FC=nvfortran --disable-shared --disable-openmp --prefix /u/${USER}/bin/sfincs_nvfortran_gpu - -make clean - -make - -make install From bb6166122ac78ca241b4cfda75eade374b7463ec Mon Sep 17 00:00:00 2001 From: Kees Nederhoff Date: Thu, 26 Mar 2026 11:18:37 -0700 Subject: [PATCH 17/19] submit progress --- source/src/sfincs_infiltration.f90 | 77 +++++++++++++++++++++--------- 1 file changed, 55 insertions(+), 22 deletions(-) diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 index 8da1855d4..d6784f47e 100644 --- a/source/src/sfincs_infiltration.f90 +++ b/source/src/sfincs_infiltration.f90 @@ -182,11 +182,28 @@ subroutine initialize_infiltration() ! (regular grids populate quadtree_nr_points and index_sfincs_in_quadtree ! via make_quadtree_from_indices) ! + ! Bucket model uses bucketfile (not infiltrationfile), but supports netcdf natively + ! if (.not. netcdf_infiltration) then ! if (use_quadtree .eqv. .true.) then ! - call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) + ! Allow bucket model with netcdf bucketfile on quadtree grids + ! + if (inftype == 'bkt' .and. bucketfile /= 'none') then + ! + if (bucketfile(len_trim(bucketfile) - 1 : len_trim(bucketfile)) /= 'nc') then + ! + call stop_sfincs('Error ! Bucket model on quadtree mesh requires a netcdf bucketfile (.nc) !', 1) + ! + endif + ! + else + ! + call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) + ! + endif + ! endif ! endif @@ -1051,7 +1068,7 @@ subroutine initialize_bucket_model() logical :: ok character*256 :: varname ! - if (bucketfile /= 'none') then + if (bucketfile /= 'none' .or. netcdf_infiltration) then ! use_bucket_model = .true. ! @@ -1068,36 +1085,52 @@ subroutine initialize_bucket_model() bucket_volume = 0.0 bucket_drain_rate = 0.0 ! - nchar = len_trim(bucketfile) - ok = check_file_exists(bucketfile, 'Bucket model file', .true.) - ! - if (bucketfile(nchar - 1 : nchar) == 'nc') then + if (netcdf_infiltration) then ! - ! Read bucket capacity (S_max) in mm, convert to m + ! Read from infiltrationfile (netcdf) - works for both regular and quadtree grids ! varname = 'bucket_smax' - call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_capacity) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_capacity) bucket_capacity = bucket_capacity / 1000.0 ! mm to m ! - ! Read drainage coefficient (k) in 1/hr, convert to 1/s - ! varname = 'bucket_k' - call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_k) + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_k) bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s ! - else - ! - ! Read from binary files - ! - open(unit = 500, file = trim(bucketfile), form = 'unformatted', access = 'stream') - read(500)bucket_capacity - close(500) - bucket_capacity = bucket_capacity / 1000.0 ! mm to m + elseif (bucketfile /= 'none') then ! - ! For binary input, k needs a separate file - not supported yet - ! Default k = 0.1/hr + nchar = len_trim(bucketfile) + ok = check_file_exists(bucketfile, 'Bucket model file', .true.) ! - bucket_k = 0.1 / 3600.0 + if (bucketfile(nchar - 1 : nchar) == 'nc') then + ! + ! Read bucket capacity (S_max) in mm, convert to m + ! + varname = 'bucket_smax' + call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_capacity) + bucket_capacity = bucket_capacity / 1000.0 ! mm to m + ! + ! Read drainage coefficient (k) in 1/hr, convert to 1/s + ! + varname = 'bucket_k' + call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_k) + bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s + ! + else + ! + ! Read from binary files + ! + open(unit = 500, file = trim(bucketfile), form = 'unformatted', access = 'stream') + read(500)bucket_capacity + close(500) + bucket_capacity = bucket_capacity / 1000.0 ! mm to m + ! + ! For binary input, k needs a separate file - not supported yet + ! Default k = 0.1/hr + ! + bucket_k = 0.1 / 3600.0 + ! + endif ! endif ! From f1c425b0be40a4b8482938b055380b1185308012 Mon Sep 17 00:00:00 2001 From: Kees Nederhoff Date: Tue, 31 Mar 2026 12:10:23 -0700 Subject: [PATCH 18/19] included loss function in bucket model --- source/src/sfincs_data.f90 | 7 +- source/src/sfincs_infiltration.f90 | 106 +++++++++++++++++++++++------ source/src/sfincs_input.f90 | 1 + source/src/sfincs_openacc.f90 | 4 +- 4 files changed, 95 insertions(+), 23 deletions(-) diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 2eb81c5a0..ef25bf796 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -419,10 +419,13 @@ module sfincs_data ! Bucket model - finite capacity reservoir with linear drainage ! logical :: use_bucket_model = .false. + real*4 :: bucket_loss_default = 0.0 ! uniform loss fraction from sfincs.inp (0-1) real*4, dimension(:), allocatable :: bucket_volume ! current storage (m) real*4, dimension(:), allocatable :: bucket_capacity ! max capacity S_max (m) real*4, dimension(:), allocatable :: bucket_k ! drainage coefficient (1/s) real*4, dimension(:), allocatable :: bucket_drain_rate ! net removal from surface this step (m/s) + real*4, dimension(:), allocatable :: bucket_loss ! loss fraction per cell (0-1), ET/deep percolation + real*4, dimension(:), allocatable :: bucket_runoff ! bucket drainage returned as surface runoff (m/s) ! ! Wind reduction for spiderweb winds ! @@ -945,7 +948,9 @@ subroutine finalize_parameters() if(allocated(bucket_capacity)) deallocate(bucket_capacity) if(allocated(bucket_k)) deallocate(bucket_k) if(allocated(bucket_drain_rate)) deallocate(bucket_drain_rate) - if(allocated(nuvisc)) deallocate(nuvisc) + if(allocated(bucket_loss)) deallocate(bucket_loss) + if(allocated(bucket_runoff)) deallocate(bucket_runoff) + if(allocated(nuvisc)) deallocate(nuvisc) ! ! Boundary velocity points ! diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 index d6784f47e..c3716496a 100644 --- a/source/src/sfincs_infiltration.f90 +++ b/source/src/sfincs_infiltration.f90 @@ -1059,12 +1059,13 @@ subroutine update_infiltration_map(dt, tloop) subroutine initialize_bucket_model() ! + use netcdf use sfincs_data use sfincs_ncinput ! implicit none ! - integer :: nchar + integer :: nchar, status, ncid, varid logical :: ok character*256 :: varname ! @@ -1079,11 +1080,15 @@ subroutine initialize_bucket_model() allocate(bucket_k(np)) allocate(bucket_volume(np)) allocate(bucket_drain_rate(np)) + allocate(bucket_loss(np)) + allocate(bucket_runoff(np)) ! - bucket_capacity = 0.0 - bucket_k = 0.0 - bucket_volume = 0.0 + bucket_capacity = 0.0 + bucket_k = 0.0 + bucket_volume = 0.0 bucket_drain_rate = 0.0 + bucket_loss = bucket_loss_default + bucket_runoff = 0.0 ! if (netcdf_infiltration) then ! @@ -1097,6 +1102,18 @@ subroutine initialize_bucket_model() call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_k) bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s ! + ! Try reading spatially-varying loss fraction (optional, falls back to uniform) + status = nf90_open(trim(infiltrationfile), NF90_NOWRITE, ncid) + if (status == nf90_noerr) then + status = nf90_inq_varid(ncid, 'bucket_loss', varid) + nchar = nf90_close(ncid) + if (status == nf90_noerr) then + varname = 'bucket_loss' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_loss) + call write_log('Info : read spatially-varying bucket_loss from infiltrationfile', 0) + endif + endif + ! elseif (bucketfile /= 'none') then ! nchar = len_trim(bucketfile) @@ -1116,6 +1133,18 @@ subroutine initialize_bucket_model() call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_k) bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s ! + ! Try reading spatially-varying loss fraction (optional, falls back to uniform) + status = nf90_open(trim(bucketfile), NF90_NOWRITE, ncid) + if (status == nf90_noerr) then + status = nf90_inq_varid(ncid, 'bucket_loss', varid) + nchar = nf90_close(ncid) + if (status == nf90_noerr) then + varname = 'bucket_loss' + call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_loss) + call write_log('Info : read spatially-varying bucket_loss from bucketfile', 0) + endif + endif + ! else ! ! Read from binary files @@ -1136,6 +1165,8 @@ subroutine initialize_bucket_model() ! write(logstr,'(a,f10.4,a)')'Info : bucket max capacity = ', maxval(bucket_capacity) * 1000.0, ' mm' call write_log(logstr, 0) + write(logstr,'(a,f6.3)')'Info : bucket loss fraction = ', maxval(bucket_loss) + call write_log(logstr, 0) ! else ! @@ -1145,10 +1176,14 @@ subroutine initialize_bucket_model() allocate(bucket_k(1)) allocate(bucket_volume(1)) allocate(bucket_drain_rate(1)) + allocate(bucket_loss(1)) + allocate(bucket_runoff(1)) bucket_capacity = 0.0 bucket_k = 0.0 bucket_volume = 0.0 bucket_drain_rate = 0.0 + bucket_loss = 0.0 + bucket_runoff = 0.0 ! endif ! @@ -1157,8 +1192,16 @@ subroutine initialize_bucket_model() subroutine compute_bucket_drainage(dt) ! - ! Bucket model: finite capacity reservoir with linear drainage (HBV/wflow style) - ! Recovery is inherent - bucket drains during dry periods via Q=k*S, restoring capacity + ! Bucket model with loss: linear reservoir + loss fraction (HBV/wflow style) + ! + ! Steps per cell: + ! 1. P_eff = P * (1 - loss) -- fraction lost to ET/deep percolation + ! 2. Fill bucket with P_eff (up to Smax capacity) + ! 3. Drain bucket: S(t+dt) = S(t)*exp(-k*dt), drainage returned as runoff + ! 4. qinfmap = P - runoff -- net removal from surface + ! + ! In continuity: zs += prcp*dt - qinfmap*dt = bucket_runoff*dt + ! => Only bucket drainage reaches the surface water level ! ! Literature: Linear reservoir (Nash, 1957), HBV soil moisture bucket (Bergstrom, 1995) ! @@ -1170,34 +1213,55 @@ subroutine compute_bucket_drainage(dt) integer :: nm real*4 :: exp_factor real*4 :: drain_vol - real*4 :: available_water + real*4 :: P_eff real*4 :: available_cap real*4 :: actual_inflow + real*4 :: precip_rate ! - !$omp parallel do private(nm, exp_factor, drain_vol, available_water, available_cap, actual_inflow) + !$omp parallel do private(nm, exp_factor, drain_vol, P_eff, available_cap, actual_inflow, precip_rate) + !$acc parallel present( kcs, prcp, qinfmap, cuminf, bucket_volume, bucket_capacity, bucket_k, & + !$acc bucket_drain_rate, bucket_loss, bucket_runoff ) + !$acc loop independent gang vector do nm = 1, np ! - if (kcs(nm) == 1 .and. bucket_capacity(nm) > 0.0) then + if (kcs(nm) == 1 .and. bucket_k(nm) > 0.0) then + ! + ! Step 1: Compute effective precipitation (after loss) + ! + precip_rate = max(prcp(nm), 0.0) + P_eff = precip_rate * (1.0 - bucket_loss(nm)) ! m/s after loss + ! + ! Step 2: Fill bucket with effective precip (up to capacity) + ! + if (bucket_capacity(nm) > 0.0) then + available_cap = bucket_capacity(nm) - bucket_volume(nm) + actual_inflow = min(P_eff * dt, available_cap) ! m + else + ! No capacity limit (Smax = 0 means infinite) + actual_inflow = P_eff * dt ! m + endif + bucket_volume(nm) = bucket_volume(nm) + actual_inflow ! - ! Step 1: Drain current storage (analytical linear reservoir solution) + ! Step 3: Drain bucket (analytical linear reservoir) ! S(t+dt) = S(t) * exp(-k*dt), drainage = S(t) - S(t+dt) ! exp_factor = exp(-bucket_k(nm) * dt) - drain_vol = bucket_volume(nm) * (1.0 - exp_factor) + drain_vol = bucket_volume(nm) * (1.0 - exp_factor) ! m drained this step bucket_volume(nm) = bucket_volume(nm) * exp_factor ! - ! Step 2: Fill bucket from available rainfall + ! Step 4: Bucket drainage becomes runoff returned to surface ! - available_water = max(prcp(nm), 0.0) * dt ! m of water available - available_cap = bucket_capacity(nm) - bucket_volume(nm) - actual_inflow = min(available_water, available_cap) - bucket_volume(nm) = bucket_volume(nm) + actual_inflow + bucket_runoff(nm) = drain_vol / dt ! m/s + ! + ! Step 5: Set qinfmap = loss + what entered bucket - what drained back + ! In continuity: zs += prcp*dt - qinfmap*dt + ! Water balance: qinfmap = prcp*loss + actual_inflow/dt - bucket_runoff + ! When bucket has room: actual_inflow = P_eff*dt => qinfmap = prcp - bucket_runoff + ! When bucket is full: actual_inflow = 0 => qinfmap can be negative (drainage > inflow) ! - ! Step 3: Set qinfmap = what entered the bucket (removed from surface) - ! This is used by continuity as the infiltration loss term + qinfmap(nm) = precip_rate * bucket_loss(nm) + actual_inflow / dt - bucket_runoff(nm) ! - qinfmap(nm) = actual_inflow / dt - bucket_drain_rate(nm) = actual_inflow / dt + bucket_drain_rate(nm) = bucket_runoff(nm) ! if (store_cumulative_precipitation) then cuminf(nm) = cuminf(nm) + qinfmap(nm) * dt @@ -1207,10 +1271,12 @@ subroutine compute_bucket_drainage(dt) ! qinfmap(nm) = 0.0 bucket_drain_rate(nm) = 0.0 + bucket_runoff(nm) = 0.0 ! endif ! enddo + !$acc end parallel !$omp end parallel do ! end subroutine diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index 34163b0bd..5f76d76b1 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -213,6 +213,7 @@ subroutine read_sfincs_input() call read_char_input(500,'infiltrationfile',infiltrationfile,'none') call read_char_input(500,'infiltrationtype',inftype,'none') call read_char_input(500,'bucketfile',bucketfile,'none') ! bucket model (infiltration flavor 'bkt') + call read_real_input(500,'bucket_loss_frac',bucket_loss_default,0.0) ! bucket loss fraction (0-1) call read_real_input(500,'qdrain',qdrain_uniform,0.0) ! drainage mimic (mm/hr) call read_char_input(500,'drainagefile',drainagefile,'none') ! spatially-varying drainage rates ! diff --git a/source/src/sfincs_openacc.f90 b/source/src/sfincs_openacc.f90 index 3d2a9cdcc..74d8a856a 100644 --- a/source/src/sfincs_openacc.f90 +++ b/source/src/sfincs_openacc.f90 @@ -36,7 +36,7 @@ subroutine initialize_openacc() !$acc gnapp2, & !$acc qinffield, qinfmap, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, & !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0, & - !$acc qdrain_rate, bucket_volume, bucket_capacity, bucket_k, bucket_drain_rate ) + !$acc qdrain_rate, bucket_volume, bucket_capacity, bucket_k, bucket_drain_rate, bucket_loss, bucket_runoff ) ! end subroutine ! @@ -66,7 +66,7 @@ subroutine finalize_openacc() !$acc gnapp2, & !$acc qinffield, qinfmap, cuminf, scs_rain, scs_Se, scs_P1, scs_F1, scs_S1, rain_T1, & !$acc ksfield, GA_head, GA_sigma, GA_sigma_max, GA_F, GA_Lu, inf_kr, horton_kd, horton_fc, horton_f0, & - !$acc qdrain_rate, bucket_volume, bucket_capacity, bucket_k, bucket_drain_rate ) + !$acc qdrain_rate, bucket_volume, bucket_capacity, bucket_k, bucket_drain_rate, bucket_loss, bucket_runoff ) ! end ! From 2a07fe88bdca891e71ae728324729104031b1fbb Mon Sep 17 00:00:00 2001 From: Kees Nederhoff Date: Wed, 15 Apr 2026 05:14:14 -0700 Subject: [PATCH 19/19] ongoing progress --- docs/input.rst | 66 ++++++++++++-- docs/parameters.rst | 40 ++++++--- source/src/sfincs_continuity.f90 | 2 +- source/src/sfincs_data.f90 | 3 - source/src/sfincs_domain.f90 | 58 +++++------- source/src/sfincs_infiltration.f90 | 138 ++++++++--------------------- source/src/sfincs_input.f90 | 50 +++++++---- source/src/sfincs_ncoutput.F90 | 6 ++ 8 files changed, 187 insertions(+), 176 deletions(-) diff --git a/docs/input.rst b/docs/input.rst index 0829c475c..f1cbd3167 100644 --- a/docs/input.rst +++ b/docs/input.rst @@ -357,14 +357,37 @@ SFINCS allows the specification of the following options for accounting for infi 3. The Curve Number method: empirical rainfall-runoff model 4. The Green-Ampt method: empirical rainfall-runoff model 5. The Horton infiltration method +6. The bucket model: linear reservoir with losses -Infiltration is specified with either constant in time values in mm/hr (both uniform and spatially varying), or using more detailed parameters for the Curve Number method, The Green-Ampt method or Horton method. +Spatially uniform infiltration is still specified directly in sfincs.inp with ``qinf``. All modern spatially varying infiltration and bucket-model input should be provided through ``infiltrationfile`` together with ``infiltrationtype``. The older binary keywords (``qinffile``, ``scsfile``, ``smaxfile``, ``sefffile``, ``ksfile``, ``psifile``, ``sigmafile``, ``f0file``, ``fcfile`` and ``kdfile``) remain available for backward compatibility only and should be removed in a future cleanup. **NOTE - Infiltration in SFINCS is only turned on when any rainfall is forced'** **NOTE - Infiltration methods in SFINCS are not designed to be stacked** +NetCDF infiltration input (recommended): +%%%%% + +For all spatially varying infiltration methods the recommended interface is: + +.. code-block:: text + + infiltrationfile = sfincs.infiltration.nc + infiltrationtype = c2d | cna | cnb | gai | hor | bkt + +The required variables in ``infiltrationfile`` depend on ``infiltrationtype``: + +* ``c2d``: ``qinf`` +* ``cna``: ``scs`` +* ``cnb``: ``smax``, ``seff``, ``ks`` +* ``gai``: ``psi``, ``sigma``, ``ks`` +* ``hor``: ``f0``, ``fc``, ``kd`` +* ``bkt``: ``bucket_smax``, ``bucket_k``, ``bucket_loss`` + +The older separate binary infiltration keywords are still supported for backward compatibility only. The former separate inputs ``bucketfile`` and ``bucket_loss_frac`` have been removed; for the bucket model, all required variables must now be present in ``infiltrationfile``. + + Spatially uniform constant in time: %%%%% @@ -381,7 +404,7 @@ Specify the keyword: Spatially varying constant in time: %%%%% -For spatially varying infiltration values per cell use the qinffile option, with the same grid based input as the depfile using a binary file. +For spatially varying infiltration values per cell use ``infiltrationfile`` with ``infiltrationtype = c2d``. The ``qinffile`` option below is kept for backward compatibility only and should be removed in a future cleanup. **qinffile = sfincs.qinf** @@ -424,7 +447,7 @@ where Smax = the soil's maximum moisture storage capacity. Smax typically derive **Without recovery** -For spatially varying infiltration values per cell using the Curve Number method without recovery use the scsfile option, with the same grid based input as the depfile using a binary file. Note here that in pre-processing the wanted CN values should be converted to S values following: +For spatially varying infiltration values per cell using the Curve Number method without recovery use ``infiltrationfile`` with ``infiltrationtype = cna``. The ``scsfile`` option below is kept for backward compatibility only and should be removed in a future cleanup. Note here that in pre-processing the wanted CN values should be converted to S values following: * scsfile: maximum soil moisture storage capacity in inches .. code-block:: text @@ -456,7 +479,7 @@ This option doesn't support restart functionality. **With recovery** -Within SFINCS, the Curve number method with recovery can be used as follows. The user needs to provide the following variables. For all variables, one needs to specify these values per cell with the same grid based input as the depfile using a binary file: +Within SFINCS, the Curve number method with recovery is preferably supplied through ``infiltrationfile`` with ``infiltrationtype = cnb``. The separate binary files listed below are kept for backward compatibility only. For all variables, one needs to specify these values per cell with the same grid based input as the depfile using a binary file: * smaxfile: maximum soil moisture storage capacity in m * sefffile: soil moisture storage capacity at the start in m @@ -496,7 +519,7 @@ The basic form of the Green-Ampt equation is expressed as follows: In which t is time, K is the saturated hydraulic conductivity, delta_theta is defined as the soil capacity (the difference between the saturated and initial moisture content) and sigma is the soil suction head. -Within SFINCS, the Green-Ampt method can be used as follows. The user needs to provide the following variables. For a range of typically values see Table 1. For all variables, one needs to specify these values per cell with the same grid based input as the depfile using a binary file: +Within SFINCS, the Green-Ampt method is preferably supplied through ``infiltrationfile`` with ``infiltrationtype = gai``. The separate binary files listed below are kept for backward compatibility only. For a range of typically values see Table 1. For all variables, one needs to specify these values per cell with the same grid based input as the depfile using a binary file: * ksfile: saturated hydraulic conductivity in mm/hr * sigmafile: soil moisture deficit in [-] @@ -520,7 +543,7 @@ The basic form of the Horton equation is expressed as follows: In which f_t is the infiltration rate at time, f_c is the final, constant infiltration rate, f_0 is the initial infiltration rate, k is a decay constant and t is the time since the start of infiltration. -Within SFINCS, the Horton method can be used as follows. The user needs to provide the following variables. For all variables, one needs to specify these values per cell with the same grid based input as the depfile using a binary file: +Within SFINCS, the Horton method is preferably supplied through ``infiltrationfile`` with ``infiltrationtype = hor``. The separate binary files listed below are kept for backward compatibility only. For all variables, one needs to specify these values per cell with the same grid based input as the depfile using a binary file: * f0file: maximum (Initial) Infiltration Capacity in mm/hr * fcfile: Minimum (Asymptotic) Infiltration Rate in mm/hr @@ -531,6 +554,37 @@ The recovery of the infiltration rate during dry weather (kr) is calculated as f This option also supports restart functionality. +The bucket model: +%%%%% + +The bucket model is a linear-reservoir representation of infiltration and losses. It is configured with: + +.. code-block:: text + + infiltrationfile = sfincs.infiltration.nc + infiltrationtype = bkt + +The ``infiltrationfile`` must contain the following variables: + +* ``bucket_smax``: maximum bucket storage in mm +* ``bucket_k``: drainage coefficient in 1/hr +* ``bucket_loss``: loss fraction in the range 0-1 + +The former separate inputs ``bucketfile`` and ``bucket_loss_frac`` are no longer supported. + + +Drainage mimic: +%%%%% + +Drainage mimic is configured separately from infiltration and now only supports ``drainagefile``: + +.. code-block:: text + + drainagefile = sfincs.drainage + +This file may be a binary map or a NetCDF file containing ``drainage_rate`` in mm/hr. The former uniform ``qdrain`` keyword has been removed. + + Storage volume ^^^^^ diff --git a/docs/parameters.rst b/docs/parameters.rst index 4ad213c52..eafb65937 100644 --- a/docs/parameters.rst +++ b/docs/parameters.rst @@ -465,54 +465,68 @@ Domain :units: s/m^(1/3) :required: no in case of regular mode, ignored in case of subgrid mode :format: bin + infiltrationfile = sfincs.infiltration.nc + :description: Recommended NetCDF input for spatially varying infiltration and bucket-model losses. Use together with infiltrationtype. + :units: depends on selected infiltrationtype and variables in the NetCDF file + :required: no + :format: net + infiltrationtype = c2d | cna | cnb | gai | hor | bkt + :description: Selects which infiltration method is read from infiltrationfile. Bucket mode requires bucket_smax, bucket_k and bucket_loss in infiltrationfile. + :units: - + :required: Only when infiltrationfile is used + :format: asc + drainagefile = sfincs.drainage + :description: Spatially varying drainage mimic input in mm/hr. Can be a binary map or a NetCDF file with variable drainage_rate. This replaces the removed qdrain keyword. + :units: mm/hr + :required: no + :format: bin or net qinffile = sfincs.qinf - :description: For spatially varying constant in time infiltration values per cell use the qinffile option, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For spatially varying constant in time infiltration values per cell prefer infiltrationfile with infiltrationtype = c2d. :units: mm/hr :required: no :format: bin scsfile = sfincs.scs - :description: For spatially varying infiltration values per cell using the Curve Number method A (without recovery) use the scsfile option, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For Curve Number method A (without recovery) prefer infiltrationfile with infiltrationtype = cna. :units: - :required: no :format: bin smaxfile = sfincs.smax - :description: For spatially varying infiltration values per cell using the Curve Number method B (with recovery) provide the smaxfile (as well as the sefffile and ksfile) as maximum soil moisture storage capacity in m, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For Curve Number method B (with recovery) prefer infiltrationfile with infiltrationtype = cnb. The smaxfile contains the maximum soil moisture storage capacity in m. :units: m :required: no :format: bin sefffile = sfincs.seff - :description: For spatially varying infiltration values per cell using the Curve Number method B (with recovery) provide the sefffile (as well as the smaxfile and ksfile) as soil moisture storage capacity at the start in m, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For Curve Number method B (with recovery) prefer infiltrationfile with infiltrationtype = cnb. The sefffile contains soil moisture storage capacity at the start in m. :units: m :required: no :format: bin ksfile = sfincs.ks - :description: For spatially varying infiltration values per cell using the Curve Number method B (with recovery) provide the ksfile (as well as the smaxfile and sefffile) as saturated hydraulic conductivity in mm/hr, with the same grid based input as the depfile using a binary file. - :description: For spatially varying infiltration values per cell using the Green & Ampt method (with recovery) provide the ksfile (as well as the sigmafile and psifile) as saturated hydraulic conductivity in mm/hr, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For Curve Number method B (with recovery) and Green & Ampt infiltration prefer infiltrationfile with infiltrationtype = cnb or gai. The ksfile contains saturated hydraulic conductivity in mm/hr. :units: mm/hr :required: no :format: bin sigmafile = sfincs.sigma - :description: For spatially varying infiltration values per cell using the Green & Ampt method (with recovery) provide the sigmafile (as well as the psifile and ksfile) as suction head at the wetting front in mm, with the same grid based input as the depfile using a binary file. - :units: mm + :description: Backward compatibility only. For Green & Ampt infiltration prefer infiltrationfile with infiltrationtype = gai. The sigmafile contains soil moisture deficit in [-]. + :units: - :required: no :format: bin psifile = sfincs.psi - :description: For spatially varying infiltration values per cell using the Green & Ampt method (with recovery) provide the psifile (as well as the sigmafile and ksfile) as soil moisture deficit in [-], with the same grid based input as the depfile using a binary file. - :units: - + :description: Backward compatibility only. For Green & Ampt infiltration prefer infiltrationfile with infiltrationtype = gai. The psifile contains suction head at the wetting front in mm. + :units: mm :required: no :format: bin f0file = sfincs.f0 - :description: For spatially varying infiltration values per cell using the Horton method (with recovery) provide the f0file (as well as the fcfile and kdfile) as maximum (Initial) Infiltration Capacity in mm/hr, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For Horton infiltration prefer infiltrationfile with infiltrationtype = hor. The f0file contains maximum (initial) infiltration capacity in mm/hr. :units: mm/hr :required: no :format: bin fcfile = sfincs.fc - :description: For spatially varying infiltration values per cell using the Horton method (with recovery) provide the fcfile (as well as the f0file and kdfile) as Minimum (Asymptotic) Infiltration Rate in mm/hr, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For Horton infiltration prefer infiltrationfile with infiltrationtype = hor. The fcfile contains the minimum (asymptotic) infiltration rate in mm/hr. :units: mm/hr :required: no :format: bin kdfile = sfincs.kd - :description: For spatially varying infiltration values per cell using the Horton method (with recovery) provide the kdfile (as well as the f0file and fcfile) as empirical constant (hr-1) of decay, with the same grid based input as the depfile using a binary file. + :description: Backward compatibility only. For Horton infiltration prefer infiltrationfile with infiltrationtype = hor. The kdfile contains the empirical decay constant in hr-1. :units: hr-1 :required: no :format: bin diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90 index d50437556..ac42e2c57 100644 --- a/source/src/sfincs_continuity.f90 +++ b/source/src/sfincs_continuity.f90 @@ -12,7 +12,7 @@ subroutine update_continuity(t, dt, tloopsrc, tloopinf, tloopcont) ! 1. Rainfall (+) => already computed in sfincs_meteo (prcp) ! 2. Infiltration (-) => computed in sfincs_infiltration (qinfmap) ! (includes: con, c2d, cna, cnb, gai, hor, bkt flavors) - ! 3. Drainage mimic (-) => simple constant rate (qdrain_rate) + ! 3. Drainage mimic (-) => drainage rate field (qdrain_rate) ! 4. External source/sink qext (+/-) => set via BMI coupling ! 5. Storage volume => depression storage (subgrid only) ! diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index ef25bf796..a2a4d4d0f 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -171,7 +171,6 @@ module sfincs_data character*256 :: fcfile character*256 :: kdfile character*256 :: drainagefile - character*256 :: bucketfile character*256 :: z0lfile character*256 :: wvmfile character*256 :: qtrfile @@ -413,13 +412,11 @@ module sfincs_data ! Drainage - constant removal rate representing subsurface drainage ! logical :: drainage = .false. - real*4 :: qdrain_uniform = 0.0 ! uniform drainage rate (mm/hr input, stored as m/s) real*4, dimension(:), allocatable :: qdrain_rate ! drainage rate per cell (m/s) ! ! Bucket model - finite capacity reservoir with linear drainage ! logical :: use_bucket_model = .false. - real*4 :: bucket_loss_default = 0.0 ! uniform loss fraction from sfincs.inp (0-1) real*4, dimension(:), allocatable :: bucket_volume ! current storage (m) real*4, dimension(:), allocatable :: bucket_capacity ! max capacity S_max (m) real*4, dimension(:), allocatable :: bucket_k ! drainage coefficient (1/s) diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 9b0383092..a2df13ae5 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -23,7 +23,7 @@ subroutine initialize_domain() ! call initialize_roughness() ! - call initialize_infiltration() ! see: sfincs_infiltration.f90 (includes bucket model if bucketfile specified) + call initialize_infiltration() ! see: sfincs_infiltration.f90 (includes bucket model if infiltrationtype='bkt') ! call initialize_drainage_mimic() ! @@ -2063,53 +2063,41 @@ subroutine initialize_drainage_mimic() ! ! Check if drainage is enabled ! - if (qdrain_uniform > 0.0 .or. drainagefile /= 'none') then + if (drainagefile /= 'none') then ! drainage = .true. ! allocate(qdrain_rate(np)) ! - if (drainagefile /= 'none') then + ! + ! Spatially-varying drainage rate + ! + write(logstr,'(a)')'Info : turning on drainage mimic (spatially-varying)' + call write_log(logstr, 0) + ! + nchar = len_trim(drainagefile) + ok = check_file_exists(drainagefile, 'Drainage file', .true.) + ! + if (drainagefile(nchar - 1 : nchar) == 'nc') then ! - ! Spatially-varying drainage rate + varname = 'drainage_rate' + call read_netcdf_quadtree_to_sfincs(drainagefile, varname, qdrain_rate) ! - write(logstr,'(a)')'Info : turning on drainage mimic (spatially-varying)' - call write_log(logstr, 0) + ! Convert from mm/hr to m/s ! - nchar = len_trim(drainagefile) - ok = check_file_exists(drainagefile, 'Drainage file', .true.) - ! - if (drainagefile(nchar - 1 : nchar) == 'nc') then - ! - varname = 'drainage_rate' - call read_netcdf_quadtree_to_sfincs(drainagefile, varname, qdrain_rate) - ! - ! Convert from mm/hr to m/s - ! - qdrain_rate = qdrain_rate / 3600.0 / 1000.0 - ! - else - ! - ! Read from binary file (assumed to be in mm/hr) - ! - open(unit = 500, file = trim(drainagefile), form = 'unformatted', access = 'stream') - read(500)qdrain_rate - close(500) - ! - ! Convert from mm/hr to m/s - ! - qdrain_rate = qdrain_rate / 3600.0 / 1000.0 - ! - endif + qdrain_rate = qdrain_rate / 3600.0 / 1000.0 ! else ! - ! Uniform drainage rate (already converted to m/s in sfincs_input.f90) + ! Read from binary file (assumed to be in mm/hr) ! - write(logstr,'(a,f10.4,a)')'Info : turning on drainage mimic (uniform, ', qdrain_uniform * 3600.0 * 1000.0, ' mm/hr)' - call write_log(logstr, 0) + open(unit = 500, file = trim(drainagefile), form = 'unformatted', access = 'stream') + read(500)qdrain_rate + close(500) + ! + ! Convert from mm/hr to m/s ! - qdrain_rate = qdrain_uniform + qdrain_rate = qdrain_rate / 3600.0 / 1000.0 ! endif ! diff --git a/source/src/sfincs_infiltration.f90 b/source/src/sfincs_infiltration.f90 index c3716496a..f020477e4 100644 --- a/source/src/sfincs_infiltration.f90 +++ b/source/src/sfincs_infiltration.f90 @@ -47,7 +47,7 @@ subroutine initialize_infiltration() ! 6) 'hor' - Modified Horton equation ! Requires: f0file or infiltrationfile ! 7) 'bkt' - Bucket model (linear reservoir, HBV/wflow style) - ! Requires: bucketfile (netcdf with bucket_smax and bucket_k) + ! Requires: infiltrationfile with bucket_smax, bucket_k and bucket_loss ! ! cumprcp and cuminf are stored in the netcdf output if store_cumulative_precipitation == .true. which is the default ! @@ -63,6 +63,12 @@ subroutine initialize_infiltration() ! 1) First we determine infiltration type ! if (precip) then + ! + if (inftype == 'bkt' .and. infiltrationfile == 'none') then + ! + call stop_sfincs('Error ! Bucket model requires infiltrationfile together with infiltrationtype = bkt !', 1) + ! + endif ! if (infiltrationfile /= 'none') then ! @@ -134,13 +140,6 @@ subroutine initialize_infiltration() infiltration = .true. store_meteo = .true. ! - elseif (bucketfile /= 'none') then - ! - ! Bucket model (linear reservoir) - ! - inftype = 'bkt' - infiltration = .true. - ! endif ! ! 2) We need cumprcp and cuminf @@ -182,27 +181,11 @@ subroutine initialize_infiltration() ! (regular grids populate quadtree_nr_points and index_sfincs_in_quadtree ! via make_quadtree_from_indices) ! - ! Bucket model uses bucketfile (not infiltrationfile), but supports netcdf natively - ! if (.not. netcdf_infiltration) then ! if (use_quadtree .eqv. .true.) then ! - ! Allow bucket model with netcdf bucketfile on quadtree grids - ! - if (inftype == 'bkt' .and. bucketfile /= 'none') then - ! - if (bucketfile(len_trim(bucketfile) - 1 : len_trim(bucketfile)) /= 'nc') then - ! - call stop_sfincs('Error ! Bucket model on quadtree mesh requires a netcdf bucketfile (.nc) !', 1) - ! - endif - ! - else - ! - call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) - ! - endif + call stop_sfincs('Error ! Infiltration input for quadtree mesh model can only be specified using the infiltrationfile Netcdf format! !', 1) ! endif ! @@ -1065,11 +1048,10 @@ subroutine initialize_bucket_model() ! implicit none ! - integer :: nchar, status, ncid, varid - logical :: ok + integer :: status, ncid, varid character*256 :: varname ! - if (bucketfile /= 'none' .or. netcdf_infiltration) then + if (netcdf_infiltration) then ! use_bucket_model = .true. ! @@ -1087,82 +1069,38 @@ subroutine initialize_bucket_model() bucket_k = 0.0 bucket_volume = 0.0 bucket_drain_rate = 0.0 - bucket_loss = bucket_loss_default + bucket_loss = 0.0 bucket_runoff = 0.0 ! - if (netcdf_infiltration) then - ! - ! Read from infiltrationfile (netcdf) - works for both regular and quadtree grids - ! - varname = 'bucket_smax' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_capacity) - bucket_capacity = bucket_capacity / 1000.0 ! mm to m - ! - varname = 'bucket_k' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_k) - bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s - ! - ! Try reading spatially-varying loss fraction (optional, falls back to uniform) - status = nf90_open(trim(infiltrationfile), NF90_NOWRITE, ncid) - if (status == nf90_noerr) then - status = nf90_inq_varid(ncid, 'bucket_loss', varid) - nchar = nf90_close(ncid) - if (status == nf90_noerr) then - varname = 'bucket_loss' - call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_loss) - call write_log('Info : read spatially-varying bucket_loss from infiltrationfile', 0) - endif - endif - ! - elseif (bucketfile /= 'none') then - ! - nchar = len_trim(bucketfile) - ok = check_file_exists(bucketfile, 'Bucket model file', .true.) - ! - if (bucketfile(nchar - 1 : nchar) == 'nc') then - ! - ! Read bucket capacity (S_max) in mm, convert to m - ! - varname = 'bucket_smax' - call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_capacity) - bucket_capacity = bucket_capacity / 1000.0 ! mm to m - ! - ! Read drainage coefficient (k) in 1/hr, convert to 1/s - ! - varname = 'bucket_k' - call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_k) - bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s - ! - ! Try reading spatially-varying loss fraction (optional, falls back to uniform) - status = nf90_open(trim(bucketfile), NF90_NOWRITE, ncid) - if (status == nf90_noerr) then - status = nf90_inq_varid(ncid, 'bucket_loss', varid) - nchar = nf90_close(ncid) - if (status == nf90_noerr) then - varname = 'bucket_loss' - call read_netcdf_quadtree_to_sfincs(bucketfile, varname, bucket_loss) - call write_log('Info : read spatially-varying bucket_loss from bucketfile', 0) - endif - endif - ! - else - ! - ! Read from binary files - ! - open(unit = 500, file = trim(bucketfile), form = 'unformatted', access = 'stream') - read(500)bucket_capacity - close(500) - bucket_capacity = bucket_capacity / 1000.0 ! mm to m - ! - ! For binary input, k needs a separate file - not supported yet - ! Default k = 0.1/hr - ! - bucket_k = 0.1 / 3600.0 - ! - endif - ! + ! + ! Read from infiltrationfile (netcdf) - works for both regular and quadtree grids + ! + varname = 'bucket_smax' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_capacity) + bucket_capacity = bucket_capacity / 1000.0 ! mm to m + ! + varname = 'bucket_k' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_k) + bucket_k = bucket_k / 3600.0 ! 1/hr to 1/s + ! + status = nf90_open(trim(infiltrationfile), NF90_NOWRITE, ncid) + if (status /= nf90_noerr) then + call stop_sfincs('Error ! Cannot open infiltrationfile for bucket model input !', 1) + endif + ! + status = nf90_inq_varid(ncid, 'bucket_loss', varid) + if (nf90_close(ncid) /= nf90_noerr) then + call stop_sfincs('Error ! Cannot close infiltrationfile after checking bucket model variables !', 1) endif ! + if (status /= nf90_noerr) then + call stop_sfincs('Error ! Bucket model requires variable bucket_loss in infiltrationfile !', 1) + endif + ! + varname = 'bucket_loss' + call read_netcdf_quadtree_to_sfincs(infiltrationfile, varname, bucket_loss) + call write_log('Info : read spatially-varying bucket_loss from infiltrationfile', 0) + ! write(logstr,'(a,f10.4,a)')'Info : bucket max capacity = ', maxval(bucket_capacity) * 1000.0, ' mm' call write_log(logstr, 0) write(logstr,'(a,f6.3)')'Info : bucket loss fraction = ', maxval(bucket_loss) diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index 5f76d76b1..be4a8bc66 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -41,6 +41,7 @@ subroutine read_sfincs_input() ! character*256 wmsigstr character*256 advstr + character*256 removed_input ! ok = check_file_exists('sfincs.inp', 'SFINCS input file', .true.) ! @@ -187,20 +188,6 @@ subroutine read_sfincs_input() call read_char_input(500,'amprfile',amprfile,'none') call read_char_input(500,'z0lfile',z0lfile,'none') call read_char_input(500,'wvmfile',wvmfile,'none') - call read_char_input(500,'qinffile',qinffile,'none') - ! Curve Number files - call read_char_input(500,'scsfile',scsfile,'none') - call read_char_input(500,'smaxfile',smaxfile,'none') - call read_char_input(500,'sefffile',sefffile,'none') - ! Green and Ampt files - call read_char_input(500,'psifile',psifile,'none') ! suction head [mm] - call read_char_input(500,'sigmafile',sigmafile,'none') ! maximum moisture deficit θdmax [-] - call read_char_input(500,'ksfile',ksfile,'none') ! saturated hydraulic conductivity [mm/hr] - ! Horton file - call read_char_input(500,'f0file',f0file,'none') ! Maximum (Initial) Infiltration Capacity, F0 - call read_char_input(500,'fcfile',fcfile,'none') ! Minimum (Asymptotic) Infiltration Rate, Fc - call read_char_input(500,'kdfile',kdfile,'none') ! k = empirical constant (hr-1) of decay - call read_real_input(500,'horton_kr_kd',horton_kr_kd,10.0) ! recovery goes 10 times as SLOW as decay ! Netcdf input call read_char_input(500,'netbndbzsbzifile',netbndbzsbzifile,'none') call read_char_input(500,'netsrcdisfile',netsrcdisfile,'none') @@ -212,10 +199,38 @@ subroutine read_sfincs_input() ! Infiltration and losses call read_char_input(500,'infiltrationfile',infiltrationfile,'none') call read_char_input(500,'infiltrationtype',inftype,'none') - call read_char_input(500,'bucketfile',bucketfile,'none') ! bucket model (infiltration flavor 'bkt') - call read_real_input(500,'bucket_loss_frac',bucket_loss_default,0.0) ! bucket loss fraction (0-1) - call read_real_input(500,'qdrain',qdrain_uniform,0.0) ! drainage mimic (mm/hr) call read_char_input(500,'drainagefile',drainagefile,'none') ! spatially-varying drainage rates + call read_char_input(500,'bucketfile',removed_input,'__removed_keyword_not_present__') + if (trim(removed_input) /= '__removed_keyword_not_present__') then + write(logstr,'(a)') 'Error : keyword bucketfile has been removed. Use infiltrationfile together with infiltrationtype = bkt.' + call stop_sfincs(trim(logstr), 1) + endif + call read_char_input(500,'bucket_loss_frac',removed_input,'__removed_keyword_not_present__') + if (trim(removed_input) /= '__removed_keyword_not_present__') then + write(logstr,'(a)') 'Error : keyword bucket_loss_frac has been removed. Add bucket_loss to infiltrationfile instead.' + call stop_sfincs(trim(logstr), 1) + endif + call read_char_input(500,'qdrain',removed_input,'__removed_keyword_not_present__') + if (trim(removed_input) /= '__removed_keyword_not_present__') then + write(logstr,'(a)') 'Error : keyword qdrain has been removed. Use drainagefile for drainage mimic input.' + call stop_sfincs(trim(logstr), 1) + endif + ! + ! Legacy binary infiltration input (backward compatibility only; remove in a future cleanup) + call read_char_input(500,'qinffile',qinffile,'none') + ! Curve Number files (legacy binary support) + call read_char_input(500,'scsfile',scsfile,'none') + call read_char_input(500,'smaxfile',smaxfile,'none') + call read_char_input(500,'sefffile',sefffile,'none') + ! Green and Ampt files (legacy binary support) + call read_char_input(500,'psifile',psifile,'none') ! suction head [mm] + call read_char_input(500,'sigmafile',sigmafile,'none') ! maximum moisture deficit theta_dmax [-] + call read_char_input(500,'ksfile',ksfile,'none') ! saturated hydraulic conductivity [mm/hr] + ! Horton files (legacy binary support) + call read_char_input(500,'f0file',f0file,'none') ! Maximum (Initial) Infiltration Capacity, F0 + call read_char_input(500,'fcfile',fcfile,'none') ! Minimum (Asymptotic) Infiltration Rate, Fc + call read_char_input(500,'kdfile',kdfile,'none') ! k = empirical constant (hr-1) of decay + call read_real_input(500,'horton_kr_kd',horton_kr_kd,10.0) ! recovery goes 10 times as SLOW as decay ! ! Output call read_char_input(500,'obsfile',obsfile,'none') @@ -318,7 +333,6 @@ subroutine read_sfincs_input() gn2 = 9.81*0.02*0.02 ! Only to be used in subgrid ! qinf = qinf/(3600*1000) - qdrain_uniform = qdrain_uniform/(3600*1000) ! Convert mm/hr to m/s ! rotation = rotation*pi/180 cosrot = cos(rotation) diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index 6a9ecdf3d..dc9e5cf99 100644 --- a/source/src/sfincs_ncoutput.F90 +++ b/source/src/sfincs_ncoutput.F90 @@ -3963,6 +3963,9 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'amvfile',amvfile)) NF90(nf90_put_att(ncid, varid, 'ampfile',ampfile)) NF90(nf90_put_att(ncid, varid, 'amprfile',amprfile)) + NF90(nf90_put_att(ncid, varid, 'infiltrationfile',infiltrationfile)) + NF90(nf90_put_att(ncid, varid, 'infiltrationtype',inftype)) + NF90(nf90_put_att(ncid, varid, 'drainagefile',drainagefile)) NF90(nf90_put_att(ncid, varid, 'qinffile',qinffile)) NF90(nf90_put_att(ncid, varid, 'scsfile',scsfile)) NF90(nf90_put_att(ncid, varid, 'smaxfile',smaxfile)) @@ -3970,6 +3973,9 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'ksfile',ksfile)) NF90(nf90_put_att(ncid, varid, 'psifile',psifile)) NF90(nf90_put_att(ncid, varid, 'sigmafile',sigmafile)) + NF90(nf90_put_att(ncid, varid, 'f0file',f0file)) + NF90(nf90_put_att(ncid, varid, 'fcfile',fcfile)) + NF90(nf90_put_att(ncid, varid, 'kdfile',kdfile)) NF90(nf90_put_att(ncid, varid, 'z0lfile',z0lfile)) NF90(nf90_put_att(ncid, varid, 'wvmfile',wvmfile)) !