Renamed bins to levels

source/sfincs/sfincs.vfproj.leijnse.user

@@ -3,5 +3,5 @@
 	<Configurations>
 		<Configuration Name="Debug|Win32"/>
 		<Configuration Name="Release|Win32"/>
-		<Configuration Name="Debug|x64" CommandArguments="sfincs.inp -autostartstop -t 1" WorkingDirectory="d:\repos\SFINCS_testbed_branch\sfincs_testbed_v2.0.3_alpha_03112023\models\03_Application\Hernani2D_advection\Input\"/>
-		<Configuration Name="Release|x64" WorkingDirectory="d:\repos\SFINCS_testbed_branch\sfincs_testbed_v2.0.3_alpha_03112023\models\03_Application\Puget_Sound_2020\Input\"/></Configurations></VisualStudioUserFile>
+		<Configuration Name="Debug|x64" CommandArguments="sfincs.inp -autostartstop -t 1" WorkingDirectory="d:\tmp\Majuro_Koen\"/>
+		<Configuration Name="Release|x64" WorkingDirectory="d:\repos\phd_waves\paper_IG_parametrisation\05_application_case\Duck_NC\01_modelsetup\ext_delilah_02_ig_final_withoutreduction_herbers\"/></Configurations></VisualStudioUserFile>

source/sfincs_lib/sfincs_lib.vfproj

@@ -50,6 +50,7 @@
 		<File RelativePath="..\src\snapwave\snapwave_data.f90"/>
 		<File RelativePath="..\src\snapwave\snapwave_date.f90"/>
 		<File RelativePath="..\src\snapwave\snapwave_domain.f90"/>
+		<File RelativePath="..\src\snapwave\snapwave_infragravity.f90"/>
 		<File RelativePath="..\src\snapwave\snapwave_ncoutput.F90">
 			<FileConfiguration Name="Release|x64">
 				<Tool Name="VFFortranCompilerTool" Preprocess="preprocessYes"/></FileConfiguration>

source/src/Makefile.am

@@ -24,6 +24,7 @@ libsfincs_la_SOURCES = \
 	snapwave/snapwave_results.f90 \
 	snapwave/snapwave_data.f90 \
 	snapwave/snapwave_boundaries.f90 \
+	snapwave/snapwave_infragravity.f90 \
 	snapwave/snapwave_date.f90 \
 	snapwave/snapwave_ncoutput.F90 \
 	snapwave/snapwave_domain.f90 \

source/src/quadtree.F90

@@ -183,7 +183,7 @@ subroutine quadtree_read_file(qtrfile)
    write(*,*)'Quadtree grid info -        dy :', quadtree_dy
    write(*,*)'Quadtree grid info -      mmax :', quadtree_mmax
    write(*,*)'Quadtree grid info -      nmax :', quadtree_nmax
-   write(*,*)'Quadtree grid info -  rotation :', quadtree_rotation
+   write(*,*)'Quadtree grid info -  rotation :', quadtree_rotation / pi * 180
    !
    end subroutine
 
@@ -363,20 +363,6 @@ subroutine quadtree_read_file_netcdf(qtrfile)
    NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%mask_varid,  quadtree_mask(:)))
    NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%snapwave_mask_varid,  quadtree_snapwave_mask(:)))
    !
-   ! Indices in netcdf file are zero-based so add 1
-   !
-   quadtree_level = quadtree_level + 1
-   quadtree_n     = quadtree_n + 1
-   quadtree_m     = quadtree_m + 1
-   quadtree_md1   = quadtree_md1 + 1
-   quadtree_md2   = quadtree_md2 + 1
-   quadtree_mu1   = quadtree_mu1 + 1
-   quadtree_mu2   = quadtree_mu2 + 1
-   quadtree_nd1   = quadtree_nd1 + 1
-   quadtree_nd2   = quadtree_nd2 + 1
-   quadtree_nu1   = quadtree_nu1 + 1
-   quadtree_nu2   = quadtree_nu2 + 1
-   !   
    ! Read attibute (should read EPSG code here ?)
    !
    NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'x0', quadtree_x0))

source/src/sfincs.f90

@@ -18,9 +18,13 @@ program sfincs
    if (ierr==0) then
       ierr = sfincs_update(deltat)
    endif
-   !
-   if (ierr==0) then
-      ierr = sfincs_finalize()
-   endif
+   ! Output:
+   !ierr = 0 ! simulation ran succesfully
+   !ierr = 1 ! simulation was stopped because of instabilities   
+   !
+   ! Always finalize, especially in case of error:
+   !if (ierr==0) then
+   ierr = sfincs_finalize()
+   !endif
    !
 end program

source/src/sfincs_boundaries.f90

@@ -418,7 +418,7 @@ subroutine find_boundary_indices()
                !
                ind1_bnd_gbp(nb)  = ib1
                ind2_bnd_gbp(nb)  = ib2
-               fac_bnd_gbp(nb) = dst2/(dst1 + dst2)
+               fac_bnd_gbp(nb) = dst2/max(dst1 + dst2, 1.0e-9)
                !
             else
                !
@@ -719,7 +719,7 @@ subroutine update_boundary_fluxes(dt)
                !
                ! Interpolation required
                !
-               dzuv    = (subgrid_uv_zmax(ip) - subgrid_uv_zmin(ip)) / (subgrid_nbins - 1)
+               dzuv    = (subgrid_uv_zmax(ip) - subgrid_uv_zmin(ip)) / (subgrid_nlevels - 1)
                iuv     = int((zsuv - subgrid_uv_zmin(ip))/dzuv) + 1
                facint  = (zsuv - (subgrid_uv_zmin(ip) + (iuv - 1)*dzuv) ) / dzuv
                depthuv = subgrid_uv_havg(iuv, ip) + (subgrid_uv_havg(iuv + 1, ip) - subgrid_uv_havg(iuv, ip))*facint

source/src/sfincs_continuity.f90

@@ -29,12 +29,12 @@ subroutine compute_water_levels(dt, tloop)
    endif  
    !
    ! Put non-default store options in a separate subroutine (all but zsmax) to save computation time when not used (both regular and subgrid):
-   if (store_maximum_velocity .eqv. .true. .or. store_maximum_flux .eqv. .true. .or. store_twet .eqv. .true.) then
+   if ((store_maximum_velocity .eqv. .true.) .or. (store_maximum_flux .eqv. .true.) .or. (store_twet .eqv. .true.)) then    
       ! 
       call compute_store_variables(dt)       
       !    
    endif
-   !
+   !   
    call system_clock(count1, count_rate, count_max)
    tloop = tloop + 1.0*(count1 - count0)/count_rate
    !
@@ -494,8 +494,8 @@ subroutine compute_water_levels_subgrid(dt)
             !
             ! Interpolation from subgrid tables needed.
             !
-            dzvol    = subgrid_z_volmax(nm) / (subgrid_nbins - 1)
-            iuv      = min(int(z_volume(nm)/dzvol) + 1, subgrid_nbins - 1)
+            dzvol    = subgrid_z_volmax(nm) / (subgrid_nlevels - 1)
+            iuv      = min(int(z_volume(nm)/dzvol) + 1, subgrid_nlevels - 1)
             facint   = (z_volume(nm) - (iuv - 1)*dzvol ) / dzvol
             zs(nm)   = subgrid_z_dep(iuv, nm) + (subgrid_z_dep(iuv + 1, nm) - subgrid_z_dep(iuv, nm))*facint                    
             !

source/src/sfincs_data.f90

@@ -7,6 +7,9 @@ module sfincs_data
       real    :: tstart_all, tfinish_all
       real*4  :: dtavg
       !!!
+      !!! Error code
+      integer :: error
+      !!!
       !!! Constants
       !!!
       real*4 g                                   ! gravitational constant g
@@ -70,7 +73,7 @@ module sfincs_data
       real*4 dym2
       real*4 nuvisc
       real*4 nuviscdim
-      real*4 nuviscinp      
+      real*4 nuviscinp   
       real*4 spw_merge_frac
       real*4 tsunami_arrival_threshold
       real*4 dtwave
@@ -216,6 +219,7 @@ module sfincs_data
       logical       :: use_spw_precip
       logical       :: friction2d
       logical       :: advection_limiter
+      logical       :: advection_mask
       !!!
       !!! sfincs_input.f90 switches
       integer storevelmax
@@ -288,6 +292,7 @@ module sfincs_data
       integer*1,          dimension(:),   allocatable :: kcs
       integer*1,          dimension(:),   allocatable :: kcuv
       ! integer*1,          dimension(:),   allocatable :: kfu
+      integer*1,          dimension(:),   allocatable :: mask_adv
       integer*1,          dimension(:),   allocatable :: scs_rain   ! logic if previous time step was raining
       !
       ! Quadtree
@@ -455,7 +460,7 @@ module sfincs_data
       !
       ! Sub-grid
       !
-      integer                             :: subgrid_nbins
+      integer                             :: subgrid_nlevels
       !
       real*4, dimension(:),   allocatable :: subgrid_z_zmin
       real*4, dimension(:),   allocatable :: subgrid_z_zmax
@@ -519,7 +524,17 @@ module sfincs_data
       real*4, dimension(:),   allocatable :: fwuv
       real*4, dimension(:),   allocatable :: mean_wave_direction
       real*4, dimension(:),   allocatable :: wave_directional_spreading
-!      real*4, dimension(:),   allocatable :: tauwavv
+      real*4, dimension(:),   allocatable :: dw
+      real*4, dimension(:),   allocatable :: df      
+      real*4, dimension(:),   allocatable :: dwig
+      real*4, dimension(:),   allocatable :: dfig
+      real*4, dimension(:),   allocatable :: cg    
+      real*4, dimension(:),   allocatable :: qb      
+      real*4, dimension(:),   allocatable :: betamean
+      real*4, dimension(:),   allocatable :: srcsh      
+      real*4, dimension(:),   allocatable :: alphaig      
+
+      !      real*4, dimension(:),   allocatable :: tauwavv
       !
       !!!
       !!! Boundary data