diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj
index c241360b..331e4d9a 100644
--- a/source/sfincs_lib/sfincs_lib.vfproj
+++ b/source/sfincs_lib/sfincs_lib.vfproj
@@ -59,7 +59,9 @@
-
+
+
+
diff --git a/source/src/quadtree.f90 b/source/src/quadtree.f90
index 6beae513..a1fcb5a3 100644
--- a/source/src/quadtree.f90
+++ b/source/src/quadtree.f90
@@ -1,7 +1,11 @@
+#define NF90(nf90call) call handle_err(nf90call,__FILE__,__LINE__)
module quadtree
-
+ !
+ use netcdf
+ !
integer*4 :: quadtree_nr_points
integer*1 :: quadtree_nr_levels
+ logical :: quadtree_netcdf
real*4 :: quadtree_x0
real*4 :: quadtree_y0
real*4 :: quadtree_dx
@@ -36,9 +40,22 @@ module quadtree
integer*4, dimension(:), allocatable :: quadtree_last_point_per_level
real*4, dimension(:), allocatable :: quadtree_dxr
real*4, dimension(:), allocatable :: quadtree_dyr
-
+ integer*1, dimension(:), allocatable :: quadtree_mask
+ integer*1, dimension(:), allocatable :: quadtree_snapwave_mask
+ !
+ type net_type_qtr
+ integer :: ncid
+ integer :: np_dimid
+ integer :: n_varid, m_varid
+ integer :: level_varid
+ integer :: nu_varid, mu_varid, nd_varid, md_varid
+ integer :: nu1_varid, mu1_varid, nd1_varid, md1_varid, nu2_varid, mu2_varid, nd2_varid, md2_varid
+ integer :: z_varid, mask_varid, snapwave_mask_varid
+ end type
+ type(net_type_qtr) :: net_file_qtr
+ !
contains
-
+ !
subroutine quadtree_read_file(qtrfile)
!
! Reads quadtree file
@@ -56,72 +73,19 @@ subroutine quadtree_read_file(qtrfile)
!
pi = 3.141592653589793
!
- ! Read quadtree file (first time, only read number of active points)
- !
- open(unit = 500, file = trim(qtrfile), form = 'unformatted', access = 'stream')
- read(500)iversion
- read(500)iepsg
- read(500)np
- read(500)quadtree_nr_levels
- close(500)
- !
- allocate(quadtree_level(np))
- allocate(quadtree_md(np))
- allocate(quadtree_md1(np))
- allocate(quadtree_md2(np))
- allocate(quadtree_mu(np))
- allocate(quadtree_mu1(np))
- allocate(quadtree_mu2(np))
- allocate(quadtree_nd(np))
- allocate(quadtree_nd1(np))
- allocate(quadtree_nd2(np))
- allocate(quadtree_nu(np))
- allocate(quadtree_nu1(np))
- allocate(quadtree_nu2(np))
- allocate(quadtree_n(np))
- allocate(quadtree_m(np))
- allocate(quadtree_n_oddeven(np))
- allocate(quadtree_m_oddeven(np))
- allocate(quadtree_xz(np))
- allocate(quadtree_yz(np))
- allocate(quadtree_zz(np))
- !
- ! Read whole quadtree file
- !
- open(unit = 500, file = trim(qtrfile), form = 'unformatted', access = 'stream')
- !
- read(500)iversion
- read(500)iepsg
- read(500)quadtree_nr_points
- read(500)quadtree_nr_levels
- quadtree_nr_levels = quadtree_nr_levels
- !
- read(500)quadtree_x0
- read(500)quadtree_y0
- read(500)quadtree_dx
- read(500)quadtree_dy
- read(500)quadtree_rotation
- !
- read(500)quadtree_level
- read(500)quadtree_n
- read(500)quadtree_m
- !
- read(500)quadtree_nu
- read(500)quadtree_nu1
- read(500)quadtree_nu2
- read(500)quadtree_mu
- read(500)quadtree_mu1
- read(500)quadtree_mu2
- read(500)quadtree_nd
- read(500)quadtree_nd1
- read(500)quadtree_nd2
- read(500)quadtree_md
- read(500)quadtree_md1
- read(500)quadtree_md2
- !
- read(500)quadtree_zz
+ quadtree_nmax = 0
+ quadtree_mmax = 0
+ !
+ quadtree_netcdf = .false.
+ if (index(qtrfile, '.nc') > 0) then
+ quadtree_netcdf = .true.
+ endif
!
- close(500)
+ if (quadtree_netcdf) then
+ call quadtree_read_file_netcdf(qtrfile)
+ else
+ call quadtree_read_file_binary(qtrfile)
+ endif
!
quadtree_rotation = quadtree_rotation*pi/180
cosrot = cos(quadtree_rotation)
@@ -142,18 +106,19 @@ subroutine quadtree_read_file(qtrfile)
quadtree_dxr = dxr
quadtree_dyr = dyr
!
- quadtree_nmax = 0
- quadtree_mmax = 0
- !
do nm = 1, quadtree_nr_points
- !
+ !
n = quadtree_n(nm)
m = quadtree_m(nm)
iref = quadtree_level(nm)
!
+ ! Determine quadtree_nmax and quadtree_mmax (these are needed for binary search of nm indices)
+ !
quadtree_nmax = max(quadtree_nmax, int( (1.0*(n - 1) + 0.01) / (2**(iref - 1))) + 2)
quadtree_mmax = max(quadtree_mmax, int( (1.0*(m - 1) + 0.01) / (2**(iref - 1))) + 2)
!
+ ! Coordinates of cell centres
+ !
quadtree_xz(nm) = quadtree_x0 + cosrot*(1.0*(m - 0.5))*dxr(iref) - sinrot*(1.0*(n - 0.5))*dyr(iref)
quadtree_yz(nm) = quadtree_y0 + sinrot*(1.0*(m - 0.5))*dxr(iref) + cosrot*(1.0*(n - 0.5))*dyr(iref)
!
@@ -210,19 +175,225 @@ subroutine quadtree_read_file(qtrfile)
enddo
!
! Some write statements of interpreted quadtree grid to output for user:
- write(*,*)'Quadtree grid info - nr_levels:',quadtree_nr_levels
- write(*,*)'Quadtree grid info - x0:',quadtree_x0
- write(*,*)'Quadtree grid info - y0:',quadtree_y0
- write(*,*)'Quadtree grid info - dx:',quadtree_dx
- 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 - nr_levels :', quadtree_nr_levels
+ write(*,*)'Quadtree grid info - x0 :', quadtree_x0
+ write(*,*)'Quadtree grid info - y0 :', quadtree_y0
+ write(*,*)'Quadtree grid info - dx :', quadtree_dx
+ 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
!
end subroutine
-
- subroutine make_quadtree_from_indices(np, indices, nmax, mmax, x0, y0, dx, dy, rotation)
+
+ subroutine quadtree_read_file_binary(qtrfile)
+ !
+ ! Reads quadtree file
+ !
+ implicit none
+ !
+ character*256, intent(in) :: qtrfile
+ !
+ integer*1 :: iversion
+ integer :: np, ip, iepsg
+ !
+ ! Read quadtree file (first time, only read number of active points)
+ !
+ write(*,*)'Reading QuadTree binary file ...'
+ !
+ open(unit = 500, file = trim(qtrfile), form = 'unformatted', access = 'stream')
+ read(500)iversion
+ read(500)iepsg
+ read(500)np
+ read(500)quadtree_nr_levels
+ close(500)
+ !
+ allocate(quadtree_level(np))
+ allocate(quadtree_md(np))
+ allocate(quadtree_md1(np))
+ allocate(quadtree_md2(np))
+ allocate(quadtree_mu(np))
+ allocate(quadtree_mu1(np))
+ allocate(quadtree_mu2(np))
+ allocate(quadtree_nd(np))
+ allocate(quadtree_nd1(np))
+ allocate(quadtree_nd2(np))
+ allocate(quadtree_nu(np))
+ allocate(quadtree_nu1(np))
+ allocate(quadtree_nu2(np))
+ allocate(quadtree_n(np))
+ allocate(quadtree_m(np))
+ allocate(quadtree_n_oddeven(np))
+ allocate(quadtree_m_oddeven(np))
+ allocate(quadtree_xz(np))
+ allocate(quadtree_yz(np))
+ allocate(quadtree_zz(np))
+ !
+ ! Read whole quadtree file
+ !
+ open(unit = 500, file = trim(qtrfile), form = 'unformatted', access = 'stream')
+ !
+ read(500)iversion
+ read(500)iepsg
+ read(500)quadtree_nr_points
+ read(500)quadtree_nr_levels
+ quadtree_nr_levels = quadtree_nr_levels
+ !
+ read(500)quadtree_x0
+ read(500)quadtree_y0
+ read(500)quadtree_dx
+ read(500)quadtree_dy
+ read(500)quadtree_rotation
+ !
+ read(500)quadtree_level
+ read(500)quadtree_n
+ read(500)quadtree_m
+ !
+ read(500)quadtree_nu
+ read(500)quadtree_nu1
+ read(500)quadtree_nu2
+ read(500)quadtree_mu
+ read(500)quadtree_mu1
+ read(500)quadtree_mu2
+ read(500)quadtree_nd
+ read(500)quadtree_nd1
+ read(500)quadtree_nd2
+ read(500)quadtree_md
+ read(500)quadtree_md1
+ read(500)quadtree_md2
+ !
+ read(500)quadtree_zz
+ !
+ close(500)
+ !
+ end subroutine
+
+
+ subroutine quadtree_read_file_netcdf(qtrfile)
+ !
+ ! Reads quadtree file from netcdf file
+ !
+ implicit none
+ !
+ character*256, intent(in) :: qtrfile
+ !
+ integer*1 :: iversion
+ integer :: np, ip, iepsg
+ !
+ write(*,*)'Reading QuadTree netCDF file ...'
+ !
+ NF90(nf90_open(trim(qtrfile), NF90_CLOBBER, net_file_qtr%ncid))
+ !
+ ! Get dimensions id's: nr points
+ !
+ NF90(nf90_inq_dimid(net_file_qtr%ncid, "mesh2d_nFaces", net_file_qtr%np_dimid))
+ !
+ ! Get dimensions sizes
+ !
+ NF90(nf90_inquire_dimension(net_file_qtr%ncid, net_file_qtr%np_dimid, len = np)) ! nr of cells
+ !
+ quadtree_nr_points = np
+ !
+ ! Get variable id's
+ !
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'n', net_file_qtr%n_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'm', net_file_qtr%m_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'level', net_file_qtr%level_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'md', net_file_qtr%md_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'md1', net_file_qtr%md1_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'md2', net_file_qtr%md2_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'mu', net_file_qtr%mu_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'mu1', net_file_qtr%mu1_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'mu2', net_file_qtr%mu2_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'nd', net_file_qtr%nd_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'nd1', net_file_qtr%nd1_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'nd2', net_file_qtr%nd2_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'nu', net_file_qtr%nu_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'nu1', net_file_qtr%nu1_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'nu2', net_file_qtr%nu2_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'z', net_file_qtr%z_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'mask', net_file_qtr%mask_varid))
+ NF90(nf90_inq_varid(net_file_qtr%ncid, 'snapwave_mask', net_file_qtr%snapwave_mask_varid))
+ !
+ ! Allocate variables
+ !
+ allocate(quadtree_level(np))
+ allocate(quadtree_md(np))
+ allocate(quadtree_md1(np))
+ allocate(quadtree_md2(np))
+ allocate(quadtree_mu(np))
+ allocate(quadtree_mu1(np))
+ allocate(quadtree_mu2(np))
+ allocate(quadtree_nd(np))
+ allocate(quadtree_nd1(np))
+ allocate(quadtree_nd2(np))
+ allocate(quadtree_nu(np))
+ allocate(quadtree_nu1(np))
+ allocate(quadtree_nu2(np))
+ allocate(quadtree_n(np))
+ allocate(quadtree_m(np))
+ allocate(quadtree_n_oddeven(np))
+ allocate(quadtree_m_oddeven(np))
+ allocate(quadtree_xz(np))
+ allocate(quadtree_yz(np))
+ allocate(quadtree_zz(np))
+ allocate(quadtree_mask(np))
+ allocate(quadtree_snapwave_mask(np))
+ !
+ ! Read values
+ !
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%n_varid, quadtree_n(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%m_varid, quadtree_m(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%level_varid, quadtree_level(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%md_varid, quadtree_md(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%md1_varid, quadtree_md1(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%md2_varid, quadtree_md2(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%mu_varid, quadtree_mu(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%mu1_varid, quadtree_mu1(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%mu2_varid, quadtree_mu2(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%nd_varid, quadtree_nd(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%nd1_varid, quadtree_nd1(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%nd2_varid, quadtree_nd2(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%nu_varid, quadtree_nu(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%nu1_varid, quadtree_nu1(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%nu2_varid, quadtree_nu2(:)))
+ NF90(nf90_get_var(net_file_qtr%ncid, net_file_qtr%z_varid, quadtree_zz(:)))
+ 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))
+ NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'y0', quadtree_y0))
+ NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'dx', quadtree_dx))
+ NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'dy', quadtree_dy))
+ NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'nmax', quadtree_nmax))
+ NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'mmax', quadtree_mmax))
+ NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'rotation', quadtree_rotation))
+ NF90(nf90_get_att(net_file_qtr%ncid, nf90_global, 'nr_levels', quadtree_nr_levels))
+ !
+ NF90(nf90_close(net_file_qtr%ncid))
+ !
+ end subroutine
+
+
+subroutine make_quadtree_from_indices(np, indices, nmax, mmax, x0, y0, dx, dy, rotation)
!
implicit none
!
@@ -252,6 +423,7 @@ subroutine make_quadtree_from_indices(np, indices, nmax, mmax, x0, y0, dx, dy, r
logical :: global
!
global = .false.
+ quadtree_netcdf = .false.
!
allocate(quadtree_level(np))
allocate(quadtree_md(np))
@@ -859,4 +1031,18 @@ subroutine quadtree_deallocate()
!
end subroutine
+
+
+ subroutine handle_err(status,file,line)
+ !
+ integer, intent ( in) :: status
+ character(*), intent(in) :: file
+ integer, intent ( in) :: line
+ integer :: status2
+ !
+ if(status /= nf90_noerr) then
+ write(0,'("NETCDF ERROR: ",a,i6,":",a)') file,line,trim(nf90_strerror(status))
+ end if
+ end subroutine handle_err
+ !
end module
diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90
index 6b178aa7..95021414 100644
--- a/source/src/sfincs_continuity.f90
+++ b/source/src/sfincs_continuity.f90
@@ -441,10 +441,14 @@ subroutine compute_water_levels_subgrid(dt)
!$acc netprcp, cumprcpt, prcp, q, z_flags_type, z_flags_iref, uv_flags_iref, &
!$acc z_index_uv_md1, z_index_uv_md2, z_index_uv_nd1, z_index_uv_nd2, z_index_uv_mu1, z_index_uv_mu2, z_index_uv_nu1, z_index_uv_nu2, &
!$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), async(1)
+ !$acc z_index_wavemaker, wavemaker_uvmean, wavemaker_nmd, wavemaker_nmu, wavemaker_ndm, wavemaker_num, storage_volume), async(1)
!$acc loop independent, private( nm )
do nm = 1, np
- !
+ !
+ ! And now water level changes due to horizontal fluxes
+ !
+ dvol = 0.0
+ !
if (kcs(nm)==1 .or. kcs(nm)==4) then
!
if (crsgeo) then
@@ -462,7 +466,8 @@ subroutine compute_water_levels_subgrid(dt)
!
if (cumprcpt(nm)>0.001 .or. cumprcpt(nm)<-0.001) then
!
- z_volume(nm) = z_volume(nm) + cumprcpt(nm)*a
+! z_volume(nm) = z_volume(nm) + cumprcpt(nm)*a
+ dvol = dvol + cumprcpt(nm)*a
cumprcpt(nm) = 0.0
!
endif
@@ -480,14 +485,17 @@ subroutine compute_water_levels_subgrid(dt)
!
if (crsgeo) then
!
- z_volume(nm) = z_volume(nm) + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxm(nm) ) * dt
+ dvol = dvol + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxm(nm) ) * dt
+! z_volume(nm) = z_volume(nm) + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxm(nm) ) * dt
!
else
!
if (use_quadtree) then
- z_volume(nm) = z_volume(nm) + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxrm(uv_flags_iref(nm)) ) * dt
+ dvol = dvol + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxrm(uv_flags_iref(nm)) ) * dt
+! z_volume(nm) = z_volume(nm) + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxrm(uv_flags_iref(nm)) ) * dt
else
- z_volume(nm) = z_volume(nm) + ( (q(nmd) - q(nmu))*dy + (q(ndm) - q(num))*dx ) * dt
+ dvol = dvol + ( (q(nmd) - q(nmu))*dy + (q(ndm) - q(num))*dx ) * dt
+! z_volume(nm) = z_volume(nm) + ( (q(nmd) - q(nmu))*dy + (q(ndm) - q(num))*dx ) * dt
endif
!
endif
@@ -507,24 +515,22 @@ subroutine compute_water_levels_subgrid(dt)
ndm2 = z_index_uv_nd2(nm)
num2 = z_index_uv_nu2(nm)
!
- dvol = 0.0
- !
! U
!
if (nmd1>0) then
- dvol = dvol + q(nmd1)*dyrm(uv_flags_iref(nmd1))
+ dvol = dvol + q(nmd1)*dyrm(uv_flags_iref(nmd1))*dt
endif
!
if (nmd2>0) then
- dvol = dvol + q(nmd2)*dyrm(uv_flags_iref(nmd2))
+ dvol = dvol + q(nmd2)*dyrm(uv_flags_iref(nmd2))*dt
endif
!
if (nmu1>0) then
- dvol = dvol - q(nmu1)*dyrm(uv_flags_iref(nmu1))
+ dvol = dvol - q(nmu1)*dyrm(uv_flags_iref(nmu1))*dt
endif
!
if (nmu2>0) then
- dvol = dvol - q(nmu2)*dyrm(uv_flags_iref(nmu2))
+ dvol = dvol - q(nmu2)*dyrm(uv_flags_iref(nmu2))*dt
endif
!
! V
@@ -532,37 +538,37 @@ subroutine compute_water_levels_subgrid(dt)
if (crsgeo) then
!
if (ndm1>0) then
- dvol = dvol + q(ndm1)*dxm(ndm1)
+ dvol = dvol + q(ndm1)*dxm(ndm1)*dt
endif
!
if (ndm2>0) then
- dvol = dvol + q(ndm2)*dxm(ndm2)
+ dvol = dvol + q(ndm2)*dxm(ndm2)*dt
endif
!
if (num1>0) then
- dvol = dvol - q(num1)*dxm(num1)
+ dvol = dvol - q(num1)*dxm(num1)*dt
endif
!
if (num2>0) then
- dvol = dvol - q(num2)*dxm(num2)
+ dvol = dvol - q(num2)*dxm(num2)*dt
endif
!
else
!
if (ndm1>0) then
- dvol = dvol + q(ndm1)*dxrm(uv_flags_iref(ndm1))
+ dvol = dvol + q(ndm1)*dxrm(uv_flags_iref(ndm1))*dt
endif
!
if (ndm2>0) then
- dvol = dvol + q(ndm2)*dxrm(uv_flags_iref(ndm2))
+ dvol = dvol + q(ndm2)*dxrm(uv_flags_iref(ndm2))*dt
endif
!
if (num1>0) then
- dvol = dvol - q(num1)*dxrm(uv_flags_iref(num1))
+ dvol = dvol - q(num1)*dxrm(uv_flags_iref(num1))*dt
endif
!
if (num2>0) then
- dvol = dvol - q(num2)*dxrm(uv_flags_iref(num2))
+ dvol = dvol - q(num2)*dxrm(uv_flags_iref(num2))*dt
endif
!
endif
@@ -574,16 +580,16 @@ subroutine compute_water_levels_subgrid(dt)
ndm = z_index_uv_nd1(nm)
num = z_index_uv_nu1(nm)
!
- dvol = 0.0
+! dvol = 0.0
!
! U
!
if (nmd>0) then
- dvol = dvol + q(nmd)*dyrm(1)
+ dvol = dvol + q(nmd)*dyrm(1)*dt
endif
!
if (nmu>0) then
- dvol = dvol - q(nmu)*dyrm(1)
+ dvol = dvol - q(nmu)*dyrm(1)*dt
endif
!
! V
@@ -591,28 +597,52 @@ subroutine compute_water_levels_subgrid(dt)
if (crsgeo) then
!
if (ndm>0) then
- dvol = dvol + q(ndm)*dxm(ndm)
+ dvol = dvol + q(ndm)*dxm(ndm)*dt
endif
!
if (num>0) then
- dvol = dvol - q(num)*dxm(num)
+ dvol = dvol - q(num)*dxm(num)*dt
endif
!
else
!
if (ndm>0) then
- dvol = dvol + q(ndm)*dxrm(1)
+ dvol = dvol + q(ndm)*dxrm(1)*dt
endif
!
if (num>0) then
- dvol = dvol - q(num)*dxrm(1)
+ dvol = dvol - q(num)*dxrm(1)*dt
endif
!
endif
!
endif
!
- z_volume(nm) = z_volume(nm) + dvol*dt
+ endif
+ !
+ if (use_storage_volume) then
+ !
+ if (storage_volume(nm)>1.0e-6) then
+ !
+ ! Still some storage left
+ !
+ z_volume(nm) = z_volume(nm) + dvol - storage_volume(nm)
+ !
+ ! Compute remaining storage volume
+ !
+ storage_volume(nm) = max(storage_volume(nm) - dvol, 0.0)
+ !
+ else
+ !
+ ! Storage is full
+ !
+ z_volume(nm) = z_volume(nm) + dvol
+ !
+ endif
+ !
+ else
+ !
+ z_volume(nm) = z_volume(nm) + dvol
!
endif
!
diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90
index ad222570..90a80f6d 100644
--- a/source/src/sfincs_data.f90
+++ b/source/src/sfincs_data.f90
@@ -137,6 +137,7 @@ module sfincs_data
character*256 :: z0lfile
character*256 :: wvmfile
character*256 :: qtrfile
+ character*256 :: volfile
!
character*256 :: trefstr_iso8601
character*41 :: treftimefews
@@ -201,6 +202,7 @@ module sfincs_data
logical :: interpolate_zst
logical :: advection
logical :: fixed_output_intervals
+ logical :: use_storage_volume
!!!
!!! sfincs_input.f90 switches
integer storevelmax
@@ -347,6 +349,10 @@ module sfincs_data
real*4, dimension(:), allocatable :: scs_T1 ! time that it is not raining
real*4, dimension(:), allocatable :: scs_Se ! S for this 'event'
!
+ ! Storage volume
+ !
+ real*4, dimension(:), allocatable :: storage_volume ! Storage volume green infra
+ !
! Wind reduction for spiderweb winds
!
real*4, dimension(:,:), allocatable :: z0land ! z0 values over land for spiderweb wind speed reduction
diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90
index 9c445e87..73bb35c0 100644
--- a/source/src/sfincs_domain.f90
+++ b/source/src/sfincs_domain.f90
@@ -132,8 +132,6 @@ subroutine initialize_domain()
!
! Read quadtree file
!
- write(*,*)'Reading quadtree file ', trim(qtrfile), ' ...'
- !
call quadtree_read_file(qtrfile)
!
else
@@ -237,9 +235,21 @@ subroutine initialize_domain()
!
else
!
- open(unit = 500, file = trim(mskfile), form = 'unformatted', access = 'stream')
- read(500)msk
- close(500)
+ if (use_quadtree .and. quadtree_netcdf) then
+ !
+ ! Mask is stored in netcdf file
+ !
+ msk = quadtree_mask
+ !
+ else
+ !
+ ! Mask in separate file
+ !
+ open(unit = 500, file = trim(mskfile), form = 'unformatted', access = 'stream')
+ read(500)msk
+ close(500)
+ !
+ endif
!
endif
!
@@ -2253,6 +2263,22 @@ subroutine initialize_domain()
store_cumulative_precipitation = .false.
!
endif
+ !
+ if (use_storage_volume) then
+ !
+ ! Spatially-varying storage volume for green infra
+ !
+ write(*,*)'Turning on process: Storage Green Infrastructure'
+ !
+ ! Read spatially-varying storage per cell (m3)
+ !
+ write(*,*)'Reading ', trim(volfile), ' ...'
+ allocate(storage_volume(np))
+ open(unit = 500, file = trim(volfile), form = 'unformatted', access = 'stream')
+ read(500)storage_volume
+ close(500)
+ !
+ endif
!
end subroutine
@@ -2637,6 +2663,14 @@ subroutine initialize_hydro()
!
endif
!
+ if (use_storage_volume) then
+ !
+ ! Make sure wet cells do not have initial storage
+ !
+ storage_volume(nm) = max(storage_volume(nm) - z_volume(nm), 0.0)
+ !
+ endif
+ !
enddo
endif
!
diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90
index e84a1392..d7c1ea1d 100644
--- a/source/src/sfincs_input.f90
+++ b/source/src/sfincs_input.f90
@@ -108,6 +108,7 @@ subroutine read_sfincs_input()
call read_char_input(500,'weirfile',weirfile,'none')
call read_char_input(500,'manningfile',manningfile,'none')
call read_char_input(500,'drnfile',drnfile,'none')
+ call read_char_input(500,'volfile',volfile,'none')
!
! Forcing
!
@@ -456,6 +457,11 @@ subroutine read_sfincs_input()
store_wave_direction = .true.
endif
!
+ use_storage_volume = .false.
+ if (volfile(1:4) /= 'none') then
+ use_storage_volume = .true.
+ endif
+ !
end subroutine
diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90
index 34830901..f5808a52 100644
--- a/source/src/sfincs_lib.f90
+++ b/source/src/sfincs_lib.f90
@@ -236,7 +236,7 @@ function sfincs_update(dtrange) result(ierr)
!$acc patm, patm0, patm1, patmb, nmindbnd, &
!$acc prcp, prcp0, prcp1, cumprcp, cumprcpt, netprcp, prcp, q, qinfmap, cuminf, &
!$acc dxminv, dxrinv, dyrinv, dxm2inv, dxr2inv, dyr2inv, dxrinvc, dxm, dxrm, dyrm, cell_area_m2, cell_area, &
- !$acc gn2uv, fcorio2d, min_dt )
+ !$acc gn2uv, fcorio2d, min_dt, storage_volume )
!
! Set target time: if dt range is negative, do not modify t1
!
diff --git a/source/src/sfincs_momentum.f90 b/source/src/sfincs_momentum.f90
index c8a5cec7..25621771 100644
--- a/source/src/sfincs_momentum.f90
+++ b/source/src/sfincs_momentum.f90
@@ -478,10 +478,24 @@ subroutine compute_fluxes(dt, min_dt, tloop)
! Apply some smoothing if theta < 1.0 (not recommended anymore!)
! Note, for reliability in terms of precision, is written as 0.9999
!
+ qsm = qx_nm
+ !
if (theta<0.9999) then
- qsm = theta*qx_nm + 0.5*(1.0 - theta)*(qx_nmu + qx_nmd)
- else
- qsm = qx_nm
+ !
+ ! Apply theta smoothing
+ !
+ if (uv_flags_adv(ip)==1) then
+ !
+ ! But only at regular points
+ !
+ if (abs(qx_nmu) > 1.0e-6 .and. abs(qx_nmd) > 1.0e-6) then
+ !
+ ! And if both uv neighbors are active
+ !
+ qsm = theta*qx_nm + 0.5*(1.0 - theta)*(qx_nmu + qx_nmd)
+ !
+ endif
+ endif
endif
!
! Compute new flux for this uv point (Bates et al., 2010)
diff --git a/source/src/snapwave/snapwave_domain.f90 b/source/src/snapwave/snapwave_domain.f90
index 69d0bbc7..896429bc 100644
--- a/source/src/snapwave/snapwave_domain.f90
+++ b/source/src/snapwave/snapwave_domain.f90
@@ -1030,13 +1030,22 @@ subroutine read_snapwave_quadtree_mesh(load_quadtree)
!
! STEP 2 - Read mask file
!
- if (mskfile /= 'none') then
+ if (quadtree_netcdf) then
!
- write(*,*)'Reading SnapWave mask file ',trim(mskfile),' ...'
- open(unit = 500, file = trim(mskfile), form = 'unformatted', access = 'stream')
- read(500)msk_tmp
- close(500)
+ ! Mask is stored in netcdf file
+ !
+ msk_tmp = quadtree_snapwave_mask
!
+ else
+ !
+ if (mskfile /= 'none') then
+ !
+ write(*,*)'Reading SnapWave mask file ',trim(mskfile),' ...'
+ open(unit = 500, file = trim(mskfile), form = 'unformatted', access = 'stream')
+ read(500)msk_tmp
+ close(500)
+ !
+ endif
endif
!
! STEP 3 - Read depth file
diff --git a/source/third_party_open/netcdf/netcdf4.vcxproj b/source/third_party_open/netcdf/netcdf4.vcxproj
index 9a454c79..235821b5 100644
--- a/source/third_party_open/netcdf/netcdf4.vcxproj
+++ b/source/third_party_open/netcdf/netcdf4.vcxproj
@@ -33,7 +33,7 @@
StaticLibrary
MultiByte
- v142
+ v143
StaticLibrary