From ccfc3ccacd1659f1a2d23e88ad40069908ba79df Mon Sep 17 00:00:00 2001 From: maartenvanormondt Date: Wed, 14 Jun 2023 09:57:20 -0400 Subject: [PATCH 1/4] added quadtree netcdf file and storage volume --- source/sfincs_lib/sfincs_lib.vfproj | 4 +- source/src/quadtree.f90 | 347 ++++++++++++++++++------ source/src/sfincs_continuity.f90 | 83 ++++-- source/src/sfincs_data.f90 | 6 + source/src/sfincs_domain.f90 | 44 ++- source/src/sfincs_input.f90 | 6 + source/src/sfincs_lib.f90 | 2 +- source/src/snapwave/snapwave_domain.f90 | 19 +- 8 files changed, 386 insertions(+), 125 deletions(-) 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..a5e19a7a 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,17 @@ 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_netcdf = .true. + quadtree_nmax = 0 + quadtree_mmax = 0 ! - close(500) + if (quadtree_netcdf) then + quadtree_netcdf = .true. + call quadtree_read_file_netcdf(qtrfile) + else + quadtree_netcdf = .false. + call quadtree_read_file_binary(qtrfile) + endif ! quadtree_rotation = quadtree_rotation*pi/180 cosrot = cos(quadtree_rotation) @@ -142,18 +104,12 @@ 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) ! - 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) - ! 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 +166,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 +414,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 +1022,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 f39ee98f..829d2c3e 100644 --- a/source/src/sfincs_continuity.f90 +++ b/source/src/sfincs_continuity.f90 @@ -443,12 +443,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 @@ -466,7 +468,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 @@ -484,14 +487,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 @@ -511,24 +517,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 @@ -536,37 +540,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 @@ -578,17 +582,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 @@ -596,28 +599,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 bf1c79c7..7e85831d 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 @@ -200,6 +201,7 @@ module sfincs_data logical :: interpolate_zst logical :: advection logical :: fixed_output_intervals + logical :: use_storage_volume !!! !!! sfincs_input.f90 switches integer storevelmax @@ -345,6 +347,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 020cf455..b42ee827 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 @@ -2632,6 +2658,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 c28cdc77..21622eb4 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 ! @@ -451,6 +452,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 3d9a2908..c7bd285d 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/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 From c6d85937e5ab975c6b2f365c1516cfa975e01f34 Mon Sep 17 00:00:00 2001 From: Maarten van Ormondt Date: Fri, 30 Jun 2023 10:22:16 +0200 Subject: [PATCH 2/4] Update netcdf4.vcxproj --- source/third_party_open/netcdf/netcdf4.vcxproj | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 2d3f6121ebe49da6cf23e70918a7369a764d2126 Mon Sep 17 00:00:00 2001 From: Maarten van Ormondt Date: Sat, 1 Jul 2023 10:28:19 +0200 Subject: [PATCH 3/4] fixed check on quadtree netcdf (if .nc in qtrfile, assume netcdf) --- source/src/quadtree.f90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/source/src/quadtree.f90 b/source/src/quadtree.f90 index a5e19a7a..13390409 100644 --- a/source/src/quadtree.f90 +++ b/source/src/quadtree.f90 @@ -73,15 +73,17 @@ subroutine quadtree_read_file(qtrfile) ! pi = 3.141592653589793 ! - quadtree_netcdf = .true. quadtree_nmax = 0 quadtree_mmax = 0 + ! + quadtree_netcdf = .false. + if (index(qtrfile, '.nc') > 0) then + quadtree_netcdf = .true. + endif ! if (quadtree_netcdf) then - quadtree_netcdf = .true. call quadtree_read_file_netcdf(qtrfile) else - quadtree_netcdf = .false. call quadtree_read_file_binary(qtrfile) endif ! From d2f68202ff44575293f0663b2b0d3191b829670c Mon Sep 17 00:00:00 2001 From: maartenvanormondt Date: Mon, 17 Jul 2023 14:17:59 -0400 Subject: [PATCH 4/4] Bug fix for theta smoothing (smoothing would take place in quadtree for points with refinement neighbors where qx_nmu and qx_nmd were not defined). Also, smoothing now no longer occurs when neighbors are dry (the way it used to be in old versions). Also fixed bug where quadtree_nmax was not determined, which resulted in WRONG obs points for quadtree grids! This error was introduced in recent version with quadtree netcdf file option. --- source/src/quadtree.f90 | 7 +++++++ source/src/sfincs_momentum.f90 | 20 +++++++++++++++++--- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/source/src/quadtree.f90 b/source/src/quadtree.f90 index 13390409..a1fcb5a3 100644 --- a/source/src/quadtree.f90 +++ b/source/src/quadtree.f90 @@ -112,6 +112,13 @@ subroutine quadtree_read_file(qtrfile) 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) ! 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)