Non-hydrostatic solver rewrite, velocity-form momentum scheme, and moving-bed source

.gitignore

@@ -65,6 +65,7 @@ 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.log
 /source/build_nvfortran_gpu_h7.sh
 /source/build_nvfortran_gpu.sh
 /source/Singularityfile-gpu.def

source/sfincs_lib/sfincs_lib.vfproj

@@ -30,16 +30,12 @@
 		<Filter Name="Resource Files" Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe"/>
 		<Filter Name="Source Files" Filter="f90;for;f;fpp;ftn;def;odl;idl">
 			<Filter Name="snapwave">
-				<File RelativePath="..\src\snapwave\interp.F90">
-				</File>
-				<File RelativePath="..\src\snapwave\snapwave_boundaries.f90">
-				</File>
+				<File RelativePath="..\src\snapwave\interp.F90"/>
+				<File RelativePath="..\src\snapwave\snapwave_boundaries.f90"/>
 				<File RelativePath="..\src\snapwave\snapwave_data.f90"/>
 				<File RelativePath="..\src\snapwave\snapwave_date.f90"/>
-				<File RelativePath="..\src\snapwave\snapwave_domain.f90">
-				</File>
-				<File RelativePath="..\src\snapwave\snapwave_infragravity.f90">
-				</File>
+				<File RelativePath="..\src\snapwave\snapwave_domain.f90"/>
+				<File RelativePath="..\src\snapwave\snapwave_infragravity.f90"/>
 				<File RelativePath="..\src\snapwave\snapwave_ncinput.F90">
 					<FileConfiguration Name="Release|x64">
 						<Tool Name="VFFortranCompilerTool" Preprocess="preprocessYes"/>
@@ -57,8 +53,7 @@
 					</FileConfiguration>
 				</File>
 				<File RelativePath="..\src\snapwave\snapwave_RFtable.f90"/>
-				<File RelativePath="..\src\snapwave\snapwave_solver.f90">
-				</File>
+				<File RelativePath="..\src\snapwave\snapwave_solver.f90"/>
 				<File RelativePath="..\src\snapwave\snapwave_windsource.f90"/>
 			</Filter>
 			<Filter Name="toml-f">
@@ -99,19 +94,15 @@
 				<File RelativePath="..\third_party_open\utils\toml-f\src\tomlf\version.f90"/>
 			</Filter>
 			<File RelativePath="..\third_party_open\Delft3D\astro.f90"/>
-			<File RelativePath="..\third_party_open\bicgstab\bicgstab_solver_ilu.f90"/>
 			<File RelativePath="..\third_party_open\utils\deg2utm.f90"/>
 			<File RelativePath="..\third_party_open\utils\geometry.f90"/>
 			<File RelativePath="..\src\sfincs_bathtub.f90"/>
 			<File RelativePath="..\src\sfincs_boundaries.f90"/>
-			<File RelativePath="..\src\sfincs_continuity.f90">
-			</File>
+			<File RelativePath="..\src\sfincs_continuity.f90"/>
 			<File RelativePath="..\src\sfincs_crosssections.f90"/>
 			<File RelativePath="..\src\sfincs_data.f90"/>
-			<File RelativePath="..\src\sfincs_data.f90"/>
 			<File RelativePath="..\src\sfincs_date.f90"/>
-			<File RelativePath="..\src\sfincs_discharges.f90">
-			</File>
+			<File RelativePath="..\src\sfincs_discharges.f90"/>
 			<File RelativePath="..\src\sfincs_domain.f90"/>
 			<File RelativePath="..\src\sfincs_error.f90"/>
 			<File RelativePath="..\src\sfincs_infiltration.f90"/>
@@ -121,6 +112,7 @@
 			<File RelativePath="..\src\sfincs_log.f90"/>
 			<File RelativePath="..\src\sfincs_meteo.f90"/>
 			<File RelativePath="..\src\sfincs_momentum.f90"/>
+			<File RelativePath="..\src\sfincs_momentum_velocity.f90"/>
 			<File RelativePath="..\src\sfincs_ncinput.F90">
 				<FileConfiguration Name="Release|x64">
 					<Tool Name="VFFortranCompilerTool" Preprocess="preprocessYes"/>
@@ -159,8 +151,7 @@
 			</File>
 			<File RelativePath="..\src\sfincs_read.f90"/>
 			<File RelativePath="..\src\sfincs_runup_gauges.f90"/>
-			<File RelativePath="..\src\sfincs_snapwave.f90">
-			</File>
+			<File RelativePath="..\src\sfincs_snapwave.f90"/>
 			<File RelativePath="..\src\sfincs_spiderweb.f90"/>
 			<File RelativePath="..\src\sfincs_structures.f90"/>
 			<File RelativePath="..\src\sfincs_subgrid.F90">
@@ -172,8 +163,7 @@
 				</FileConfiguration>
 			</File>
 			<File RelativePath="..\src\sfincs_timestep_analysis.f90"/>
-			<File RelativePath="..\src\sfincs_vegetation.f90">
-			</File>
+			<File RelativePath="..\src\sfincs_vegetation.f90"/>
 			<File RelativePath="..\src\sfincs_wave_enhanced_roughness.f90"/>
 			<File RelativePath="..\src\sfincs_wavemaker.f90"/>
 		</Filter>

source/src/Makefile.am

@@ -88,12 +88,12 @@ libsfincs_la_SOURCES = \
 	sfincs_snapwave.f90 \
 	../third_party_open/utils/deg2utm.f90 \
 	sfincs_meteo.f90 \
-        ../third_party_open/bicgstab/bicgstab_solver_ilu.f90 \
 	sfincs_nonhydrostatic.f90 \
 	sfincs_ncoutput_helpers.F90 \
 	sfincs_ncoutput.F90 \
 	sfincs_output.f90 \
 	sfincs_momentum.f90 \
+	sfincs_momentum_velocity.f90 \
 	sfincs_wavemaker.f90 \
 	sfincs_bathtub.f90 \
 	sfincs_openacc.f90 \

source/src/sfincs.f90

@@ -15,8 +15,9 @@ program sfincs
       !
       ! Set BMI flags to false 
       ! 
-      bmi = .false. 
+      bmi = .false.
       use_qext = .false.
+      use_dzbext = .false.
       ! 
       ierr = sfincs_initialize()
       ! 

source/src/sfincs_bmi.f90

@@ -23,7 +23,7 @@ module sfincs_bmi
    public :: get_end_time
    public :: get_time_step
    public :: get_current_time
-   public :: update_zbuv
+   public :: bmi_update_bed_level
    public :: update_apparent_roughness
    public :: get_sfincs_cell_index
    public :: get_sfincs_cell_indices
@@ -119,6 +119,12 @@ function get_value_ptr(c_var_name, c_data) result(ierr) &
          c_data = c_loc(subgrid_z_zmin)
       case("qext")
          c_data = c_loc(qext)
+      case("dzbext")
+         if (.not. allocated(dzbext)) then
+            allocate(dzbext(np))
+            dzbext = 0.0
+         endif
+         c_data = c_loc(dzbext)
       case("uorb")
          c_data = c_loc(uorb)
       case default
@@ -150,6 +156,12 @@ function get_var_shape(c_var_name, var_shape) result(ierr) &
          var_shape(1) = size(z_index_z_n)
       case("qext")
          var_shape(1) = size(qext)
+      case("dzbext")
+         if (.not. allocated(dzbext)) then
+            allocate(dzbext(np))
+            dzbext = 0.0
+         endif
+         var_shape(1) = size(dzbext)
       case default
          ierr = -1
       end select
@@ -170,7 +182,7 @@ function get_var_type(c_var_name, c_type) result(ierr) &
       var_name = char_array_to_string(c_var_name, strlen(c_var_name, BMI_LENVARADDRESS))
 
       select case(var_name)
-      case("z_xz", "z_yz", "zb", "subgrid_z_zmin", "qext", "uorb")
+      case("z_xz", "z_yz", "zb", "subgrid_z_zmin", "qext", "dzbext", "uorb")
          type_name = "float"
       case("zs")
          type_name = "double"
@@ -198,7 +210,7 @@ function get_var_rank(c_var_name, rank) result(ierr) &
       var_name = char_array_to_string(c_var_name, strlen(c_var_name, BMI_LENVARADDRESS))
 
       select case(var_name)
-      case("z_xz", "z_yz", "zs", "zb", "subgrid_z_zmin", "qext", "uorb")
+      case("z_xz", "z_yz", "zs", "zb", "subgrid_z_zmin", "qext", "dzbext", "uorb")
          rank = 1
       case default
          ierr = -1
@@ -227,7 +239,23 @@ function set_logical(c_flag_name, ival) result(ierr) bind(C, name="set_logical")
       select case(flag_name)
       case("qext")
          use_qext = bval
-         !write(*,*)'use_qext = ', use_qext 
+         !write(*,*)'use_qext = ', use_qext
+      case("dzbext")
+         !
+         ! Lazily allocate the external delta-bed-level array on first enable.
+         ! Once allocated we keep it around; toggling the flag off later just
+         ! disables the update path without freeing memory (same pattern as
+         ! qext is handled elsewhere).
+         !
+         if (bval .and. .not. allocated(dzbext)) then
+            !
+            allocate(dzbext(np))
+            dzbext = 0.0
+            !
+         endif
+         !
+         use_dzbext = bval
+         !
       case default
          ierr = -1
       end select
@@ -274,14 +302,15 @@ function get_current_time(tcurrent) result(ierr) bind(C, name="get_current_time"
 
    end function get_current_time
 
-   function update_zbuv() result(ierr) bind(C, name="update_zbuv")
-   ! Update bed level at uv points
-   !DEC$ ATTRIBUTES DLLEXPORT :: update_zbuv
+   function bmi_update_bed_level() result(ierr) bind(C, name="update_bed_level")
+   ! Apply dzbext to the bed-level arrays (zb or subgrid_z_zmin/zmax and
+   ! subgrid_uv_zmin/zmax) and rebuild zbuvmx for non-subgrid runs.
+   !DEC$ ATTRIBUTES DLLEXPORT :: bmi_update_bed_level
       integer(kind=c_int) :: ierr
-      call compute_zbuvmx()
+      call update_bed_level()
       ierr = 0
-   
-   end function update_zbuv
+
+   end function bmi_update_bed_level
 
    function update_apparent_roughness() result(ierr) bind(C, name="update_apparent_roughness")
    ! Update apparent roughness at uv points

source/src/sfincs_continuity.f90

@@ -569,6 +569,15 @@ subroutine compute_water_levels_subgrid(dt,t)
             !
          endif
          !
+         ! Effective bed for subgrid models with velocity-form advection and/or
+         ! non-hydrostatics: zb = zs - V/A, the bed consistent with the volume
+         ! continuity conserves (the file zb is not a valid conveyance bed).
+         ! Bed SLOPES elsewhere stay frozen at their initialization values.
+         !
+         if (zb_effective) then
+            zb(nm) = zs(nm) - max(z_volume(nm) / a, 0.0)
+         endif
+         !
          !
          if (wiggle_suppression) then 
             ! 
@@ -618,9 +627,9 @@ subroutine compute_water_levels_subgrid(dt,t)
    enddo
    !$omp end do
    !$omp end parallel
-   !         
+   !
    !$acc end parallel
-   !         
+   !
    end subroutine
 
    subroutine compute_store_variables(dt)

source/src/sfincs_data.f90

@@ -16,6 +16,7 @@ module sfincs_data
       !!!
       logical       :: bmi
       logical       :: use_qext
+      logical       :: use_dzbext
       !!!
       !!! Constants
       !!!
@@ -97,10 +98,7 @@ module sfincs_data
       !real*4 dzdsbnd
       !real*4 manningbnd
       real*4 nuviscfac ! Factor on viscosity for 'difficult' points. Used in sfincs_momentum.f90.
-      real*4 nh_fnudge
-      real*4 nh_tstop
-      integer nh_itermax
-      real*4 nh_tol
+      logical zb_effective ! subgrid + (velocity scheme or nonh): zb is recomputed every step in continuity as the effective bed zs - z_volume/area
       real*4 runup_gauge_depth
       real*4 factor_wind
       real*4 factor_pres
@@ -114,6 +112,7 @@ module sfincs_data
       integer cd_nr
       integer baro
       integer advection_scheme
+      integer momentum_scheme   ! 0 = Bates flux form (default), 1 = velocity form
       !
       character*256 :: depfile
       character*256 :: mskfile
@@ -336,7 +335,6 @@ module sfincs_data
       integer*1,          dimension(:),   allocatable :: kfuv
       integer*1,          dimension(:),   allocatable :: mask_adv
       integer*1,          dimension(:),   allocatable :: scs_rain   ! logic if previous time step was raining
-      integer*1,          dimension(:),   allocatable :: mask_nonh
       !
       ! Quadtree
       !
@@ -597,6 +595,7 @@ module sfincs_data
       real*4, dimension(:),   allocatable :: zs0
       real*4, dimension(:),   allocatable :: zsderv
       real*4, dimension(:),   allocatable, target :: qext
+      real*4, dimension(:),   allocatable, target :: dzbext
       real*4, dimension(:),   allocatable, target :: uorb
       real*4, dimension(:),   allocatable :: gnapp2
       !
@@ -1018,6 +1017,7 @@ subroutine finalize_parameters()
     if(allocated(uv0)) deallocate(uv0)
     if(allocated(twet)) deallocate(twet)
     if(allocated(qext)) deallocate(qext)
+    if(allocated(dzbext)) deallocate(dzbext)
     !
 !    if(allocated(huu)) deallocate(huu)
 !    if(allocated(hvv)) deallocate(hvv)