From f4aa9fe07ce8383c11d8c050f527fa9554b22a30 Mon Sep 17 00:00:00 2001 From: Tim Leijnse Date: Thu, 9 Nov 2023 17:03:59 +0000 Subject: [PATCH] 56 make beta version for release (#58) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Prepare the v2.0.3-Cauberg release of SFINCS! Merged PR of Netcdf4, Netcdf spiderwebs and Quadtree updates, plus some small needed fixes to make it all working. In more detail: * First progress: netcdf data is read. Next step: interpolation in rest of SFINCS * working version of SFINCS with netcdf spiderweb (rainfall not tested yet) * made integers in sfincs_input.f90 temporary again * removed variable SFINCS doesnt need * major update of sfincs_domain and fluxes in case of quadtree * enabled combo of wave makers and subgrid * cosmetic updates * computing combined fluxes after reading restart file * added deallocation of quadtree arrays * cosmetic * Update netcdf4.vcxproj back to 142 (again ....) * storing logicals as ints in netcdf output attributes * Use NETCDF4 + nc_deflate_level * - add netspwfile to output var list * updated documentation for netcdf spiderwebs * - Fix quadtree map nc output file, 'NF90_64BIT_OFFSET' didn't work (/together with NF90_NETCDF4 compression), not sure why was needed * - We also need 'uv_index_u_nmd' and 'uv_index_u_nmu' in sfincs_momentum.f90 for exception (theta<0.9999 .and. .not.advection) for backward compatibility - For now done by always executing if-loops * - In case of inputtype 'asc', we need variable 'kcsg' from initialize_mesh() still in initialize_bathymetry(), because otherwise you loop over all nmax*mmax, which doesn't match with ip - Fixed by making kcsg a general available variable in sfincs_data.f90 * - Needed extra flag check for advection as discussed with Maarten * - V2.0.3-Cauberg! - Replace the (theta<0.9999) by a dedicated flag 'thetasmoothing' - That flag now used in sfincs_domain.f90 to determine whether extra uv indices should be determined, and in sfincs_momentum.f90 for calculating fluxes for smoothing - Redo commit #d2f6820 "Bug fix for theta smoothing [..] smoothing now no longer occurs when neighbors are dry (the way it used to be in old versions)." > was unintentianally removed in PR#47 > solved by adding within if thetasmoothing "if ( kcuv(uv_index_u_nmd(ip))==1 .and. kcuv(uv_index_u_nmu(ip))==1 ) then " * - Making horton infiltration a different option 6 in rstfile, as should have been done so in the first place in commit # 594218b * - Bump version to v2.0.3-Cauberg * - Make spw conversion to utm also available for netspwfile input --------- Co-authored-by: Kees Nederhoff Co-authored-by: maartenvanormondt Co-authored-by: David Gutiérrez-Barceló --- README.rst | 2 +- docs/conf.py | 4 +- docs/input_forcing.rst | 20 +- ...CENSING CONDITIONS for Free Trial Copy.txt | 2 +- source/sfincs/sfincs.vfproj.leijnse.user | 4 +- ...CENSING CONDITIONS for Free Trial Copy.txt | 2 +- source/src/sfincs_continuity.f90 | 509 ++---- source/src/sfincs_data.f90 | 67 +- source/src/sfincs_date.f90 | 36 + source/src/sfincs_domain.f90 | 1611 +++++++++-------- source/src/sfincs_input.f90 | 76 +- source/src/sfincs_lib.f90 | 16 +- source/src/sfincs_meteo.f90 | 46 +- source/src/sfincs_momentum.f90 | 173 +- source/src/sfincs_ncinput.F90 | 132 ++ source/src/sfincs_ncoutput.F90 | 381 ++-- source/src/sfincs_obspoints.f90 | 2 +- source/src/sfincs_output.f90 | 17 +- source/src/sfincs_spiderweb.f90 | 3 +- source/src/sfincs_wavemaker.f90 | 25 +- source/src/snapwave/snapwave_boundaries.f90 | 2 +- .../third_party_open/netcdf/netcdf4.vcxproj | 2 +- 22 files changed, 1660 insertions(+), 1472 deletions(-) diff --git a/README.rst b/README.rst index 43d145c3..ff919415 100644 --- a/README.rst +++ b/README.rst @@ -54,7 +54,7 @@ SFINCS seeks active contribution from the hydro modelling community, so feel fre :alt: Latest developers docs .. |docs_stable| image:: https://img.shields.io/badge/docs-stable-brightgreen.svg - :target: https://sfincs.readthedocs.io/en/v2.0.2_blockhaus_release/ + :target: https://sfincs.readthedocs.io/en/v2.0.3_cauberg_release/ :alt: Stable docs last release .. |license| image:: https://img.shields.io/github/license/Deltares/SFINCS diff --git a/docs/conf.py b/docs/conf.py index eaaf115f..52bfa3ca 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -71,9 +71,9 @@ def setup(app): # built documents. # # The short X.Y version. -version = '2.0.2_Blockhaus' +version = '2.0.3_Cauberg' # The full version, including alpha/beta/rc tags. -release = '2.0.2_Blockhaus' +release = '2.0.3_Cauberg' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/docs/input_forcing.rst b/docs/input_forcing.rst index 35c2b778..1722301d 100644 --- a/docs/input_forcing.rst +++ b/docs/input_forcing.rst @@ -249,11 +249,11 @@ Meteo There are a few different options to specify wind and rain input: -1) Use a spatially varying spiderweb input (as in Delft3D/Delft3D FM) for forcing tropical cyclones only the wind and pressure input, or for the wind as well as the rain input. +1) Use a spatially varying input in a polar coordinate system ('spiderweb'; as in Delft3D/Delft3D FM) for forcing tropical cyclones for wind and pressure input only, or also for rainfal. -2) Use a spatially varying grid input (as in Delft3D) for u- and v- wind velocities and/or the rain and/or pressure input. +2) Use a spatially varying grid input (as in Delft3D) for u- and v- wind velocities and/or the rain and/or pressure input. -3) Use a spatially varying grid input using a netcdf file based on a FEWS input type format for wind, rain and/or atmospheric pressure input. +3) Use a spatially varying grid input using a netcdf file based on a FEWS input type format for wind, rain, and/or atmospheric pressure input. 4) Use a spatially uniform input for wind and rain, which is faster but also more simplified. @@ -270,15 +270,23 @@ You can know how much rainfall / wind is added to the model in the output by spe Spatially varying spiderweb ^^^^^^^^^ -The option of forcing spiderweb files is only relevant for tropical cyclones, best is to put grid units in the same projected coordinate reference system (UTM zone) as the SFINCS grid. -For generation of these spiderweb files use Deltares' Wind Enhancement Scheme tool (WES, see https://content.oss.deltares.nl/delft3d/manuals/Delft3D-WES_User_Manual.pdf or OET Matlab equivalent) or get in touch. +The option of forcing spiderweb files is only relevant for tropical cyclones, as these time and spatially varying wind fields are constructed and written in a polar coordinate system. +For generation of these spiderweb files use Deltares' Wind Enhancement Scheme tool (WES, see https://content.oss.deltares.nl/delft3d/manuals/Delft3D-WES_User_Manual.pdf, OET Matlab equivalent, the Coastal Hazards Toolkit or other methods). If you have issues with the generation of a spiderweb, feel free to get in touch. +There are two options for spatially varying spiderweb. There is the 'traditional' ASCII format and the netcdf option. Both options support the possibility to include rainfall too. -**Spiderweb-input:** +**Spiderweb-input-ascii:** .. code-block:: text spwfile = tropical_cyclone.spw +**Spiderweb-input-netcdf:** + +.. code-block:: text + + netspwfile = tropical_cyclone.nc + + Spatially varying gridded ^^^^^^^^^ diff --git a/source/LICENSING CONDITIONS for Free Trial Copy.txt b/source/LICENSING CONDITIONS for Free Trial Copy.txt index f9ce7c71..43c36c16 100644 --- a/source/LICENSING CONDITIONS for Free Trial Copy.txt +++ b/source/LICENSING CONDITIONS for Free Trial Copy.txt @@ -3,7 +3,7 @@ For the Free Trial Copy Read this agreement -This is a legal agreement between the prospective user, either an individual or an entity, (hereafter: Licensee) and Stichting Deltares hereafter: Deltares. This Software [SFINCS v2.0.2 Blockhaus release 09-06-2023] is provided to you on the suspensive condition that you accept the following agreement. You are not allowed to copy, distribute, install or use the software without an agreement with Deltares. If you do not agree to the terms of this agreement, you must promptly stop installing this software and/or do not download the software from our site. Any further actions with respect to the software without an agreement with Deltares are in violation of copyright law. Downloading, installation and/or use of the Software implies acceptance of this agreement. +This is a legal agreement between the prospective user, either an individual or an entity, (hereafter: Licensee) and Stichting Deltares hereafter: Deltares. This Software [SFINCS v2.0.3 Cauberg release 15-11-2023] is provided to you on the suspensive condition that you accept the following agreement. You are not allowed to copy, distribute, install or use the software without an agreement with Deltares. If you do not agree to the terms of this agreement, you must promptly stop installing this software and/or do not download the software from our site. Any further actions with respect to the software without an agreement with Deltares are in violation of copyright law. Downloading, installation and/or use of the Software implies acceptance of this agreement. Article 1: Evaluation license diff --git a/source/sfincs/sfincs.vfproj.leijnse.user b/source/sfincs/sfincs.vfproj.leijnse.user index df720f15..730e6e33 100644 --- a/source/sfincs/sfincs.vfproj.leijnse.user +++ b/source/sfincs/sfincs.vfproj.leijnse.user @@ -3,5 +3,5 @@ - - + + diff --git a/source/sfincs/x64/Release/LICENSING CONDITIONS for Free Trial Copy.txt b/source/sfincs/x64/Release/LICENSING CONDITIONS for Free Trial Copy.txt index f9ce7c71..43c36c16 100644 --- a/source/sfincs/x64/Release/LICENSING CONDITIONS for Free Trial Copy.txt +++ b/source/sfincs/x64/Release/LICENSING CONDITIONS for Free Trial Copy.txt @@ -3,7 +3,7 @@ For the Free Trial Copy Read this agreement -This is a legal agreement between the prospective user, either an individual or an entity, (hereafter: Licensee) and Stichting Deltares hereafter: Deltares. This Software [SFINCS v2.0.2 Blockhaus release 09-06-2023] is provided to you on the suspensive condition that you accept the following agreement. You are not allowed to copy, distribute, install or use the software without an agreement with Deltares. If you do not agree to the terms of this agreement, you must promptly stop installing this software and/or do not download the software from our site. Any further actions with respect to the software without an agreement with Deltares are in violation of copyright law. Downloading, installation and/or use of the Software implies acceptance of this agreement. +This is a legal agreement between the prospective user, either an individual or an entity, (hereafter: Licensee) and Stichting Deltares hereafter: Deltares. This Software [SFINCS v2.0.3 Cauberg release 15-11-2023] is provided to you on the suspensive condition that you accept the following agreement. You are not allowed to copy, distribute, install or use the software without an agreement with Deltares. If you do not agree to the terms of this agreement, you must promptly stop installing this software and/or do not download the software from our site. Any further actions with respect to the software without an agreement with Deltares are in violation of copyright law. Downloading, installation and/or use of the Software implies acceptance of this agreement. Article 1: Evaluation license diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90 index 9c4a0eca..e19e180f 100644 --- a/source/src/sfincs_continuity.f90 +++ b/source/src/sfincs_continuity.f90 @@ -108,17 +108,17 @@ subroutine compute_water_levels_regular(dt) endif ! !$omp parallel & - !$omp private ( nm,dvol,nmd1,nmu1,ndm1,num1,nmd2,nmu2,ndm2,num2,nmd,nmu,ndm,num,qnmd,qnmu,qndm,qnum,iwm) + !$omp private ( nm,dvol,nmd,nmu,ndm,num,qnmd,qnmu,qndm,qnum,iwm) !$omp do schedule ( dynamic, 256 ) !$acc kernels present( kcs, zs, zb, netprcp, cumprcpt, prcp, q, zsmax, zsm, & - !$acc 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 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 z_index_wavemaker, wavemaker_uvmean, wavemaker_nmd, wavemaker_nmu, wavemaker_ndm, wavemaker_num), async(1) !$acc loop independent, private( nm ) do nm = 1, np ! - if (kcs(nm)==1) then + if (kcs(nm)==1) then ! Regular point ! if (precip) then ! @@ -137,156 +137,16 @@ subroutine compute_water_levels_regular(dt) ! endif ! - if (z_flags_type(nm) == 0) then - ! - ! Regular point with four surrounding cells of the same size - ! - nmd = z_index_uv_md1(nm) - nmu = z_index_uv_mu1(nm) - ndm = z_index_uv_nd1(nm) - num = z_index_uv_nu1(nm) - ! - if (crsgeo) then - zs(nm) = zs(nm) + (((q(nmd) - q(nmu))*dxminv(nmu) + (q(ndm) - q(num))*dyrinv(z_flags_iref(nm))))*dt - else - zs(nm) = zs(nm) + (((q(nmd) - q(nmu))*dxrinv(z_flags_iref(nm)) + (q(ndm) - q(num))*dyrinv(z_flags_iref(nm))))*dt - endif - ! - else - ! - if (use_quadtree) then - ! - ! More complicated - ! - nmd1 = z_index_uv_md1(nm) - nmu1 = z_index_uv_mu1(nm) - ndm1 = z_index_uv_nd1(nm) - num1 = z_index_uv_nu1(nm) - nmd2 = z_index_uv_md2(nm) - nmu2 = z_index_uv_mu2(nm) - 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)) - endif - ! - if (nmd2>0) then - dvol = dvol + q(nmd2)*dyrm(uv_flags_iref(nmd2)) - endif - ! - if (nmu1>0) then - dvol = dvol - q(nmu1)*dyrm(uv_flags_iref(nmu1)) - endif - ! - if (nmu2>0) then - dvol = dvol - q(nmu2)*dyrm(uv_flags_iref(nmu2)) - endif - ! - ! V - ! - if (crsgeo) then - ! - if (ndm1>0) then - dvol = dvol + q(ndm1)*dxm(ndm1) - endif - ! - if (ndm2>0) then - dvol = dvol + q(ndm2)*dxm(ndm2) - endif - ! - if (num1>0) then - dvol = dvol - q(num1)*dxm(num1) - endif - ! - if (num2>0) then - dvol = dvol - q(num2)*dxm(num2) - endif - ! - else - ! - if (ndm1>0) then - dvol = dvol + q(ndm1)*dxrm(uv_flags_iref(ndm1)) - endif - ! - if (ndm2>0) then - dvol = dvol + q(ndm2)*dxrm(uv_flags_iref(ndm2)) - endif - ! - if (num1>0) then - dvol = dvol - q(num1)*dxrm(uv_flags_iref(num1)) - endif - ! - if (num2>0) then - dvol = dvol - q(num2)*dxrm(uv_flags_iref(num2)) - endif - ! - endif - ! - else - ! - nmd = z_index_uv_md1(nm) - nmu = z_index_uv_mu1(nm) - ndm = z_index_uv_nd1(nm) - num = z_index_uv_nu1(nm) - ! - dvol = 0.0 - ! - ! U - ! - if (nmd>0) then - dvol = dvol + q(nmd)*dyrm(1) - endif - ! - if (nmu>0) then - dvol = dvol - q(nmu)*dyrm(1) - endif - ! - ! V - ! - if (crsgeo) then - ! - if (ndm>0) then - dvol = dvol + q(ndm)*dxm(ndm) - endif - ! - if (num>0) then - dvol = dvol - q(num)*dxm(num) - endif - ! - else - ! - if (ndm>0) then - dvol = dvol + q(ndm)*dxrm(1) - endif - ! - if (num>0) then - dvol = dvol - q(num)*dxrm(1) - endif - ! - endif - ! - endif - ! - if (crsgeo) then - zs(nm) = zs(nm) + dvol*dt/cell_area_m2(nm) - else - zs(nm) = zs(nm) + dvol*dt/cell_area(z_flags_iref(nm)) - endif - ! - endif + nmd = z_index_uv_md(nm) + nmu = z_index_uv_mu(nm) + ndm = z_index_uv_nd(nm) + num = z_index_uv_nu(nm) ! - ! No continuity update but keeping track of variables - ! zsmax used by default, therefore keep in standard continuity loop: - if (store_maximum_waterlevel) then - ! - zsmax(nm) = max(zsmax(nm), zs(nm)) - ! - endif + if (crsgeo) then + zs(nm) = zs(nm) + (((q(nmd) - q(nmu))*dxminv(nmu) + (q(ndm) - q(num))*dyrinv(z_flags_iref(nm))))*dt + else + zs(nm) = zs(nm) + (((q(nmd) - q(nmu))*dxrinv(z_flags_iref(nm)) + (q(ndm) - q(num))*dyrinv(z_flags_iref(nm))))*dt + endif ! endif ! @@ -294,7 +154,8 @@ subroutine compute_water_levels_regular(dt) ! if (kcs(nm)==4) then ! - ! Wave maker point + ! Wave maker point (seaward of wave maker) + ! Here we use the mean flux at the location of the wave maker ! iwm = z_index_wavemaker(nm) ! @@ -306,7 +167,7 @@ subroutine compute_water_levels_regular(dt) ! else ! - qnmd = q(z_index_uv_md1(nm)) + qnmd = q(z_index_uv_md(nm)) ! endif ! @@ -318,7 +179,7 @@ subroutine compute_water_levels_regular(dt) ! else ! - qnmu = q(z_index_uv_mu1(nm)) + qnmu = q(z_index_uv_mu(nm)) ! endif ! @@ -330,7 +191,7 @@ subroutine compute_water_levels_regular(dt) ! else ! - qndm = q(z_index_uv_nd1(nm)) + qndm = q(z_index_uv_nd(nm)) ! endif ! @@ -342,7 +203,7 @@ subroutine compute_water_levels_regular(dt) ! else ! - qnum = q(z_index_uv_nu1(nm)) + qnum = q(z_index_uv_nu(nm)) ! endif ! @@ -350,11 +211,26 @@ subroutine compute_water_levels_regular(dt) ! endif ! + endif + ! + if (snapwave) then + ! + ! Time-averaged water level used for SnapWave using exponential filter + ! ! Would double exponential filtering be better? ! zsm(nm) = factime*zs(nm) + (1.0 - factime)*zsm(nm) ! - endif + endif + ! + ! No continuity update but keeping track of variables + ! zsmax used by default, therefore keep in standard continuity loop: + ! + if (store_maximum_waterlevel) then + ! + zsmax(nm) = max(zsmax(nm), zs(nm)) + ! + endif ! enddo !$omp end do @@ -407,6 +283,11 @@ subroutine compute_water_levels_subgrid(dt) real*4 :: factime real*4 :: dvol ! + real*4 :: qnmu + real*4 :: qnmd + real*4 :: qnum + real*4 :: qndm + ! integer :: iuv real*4 :: dzvol real*4 :: facint @@ -434,12 +315,12 @@ subroutine compute_water_levels_subgrid(dt) endif ! !$omp parallel & - !$omp private ( dvol,nmd1,nmu1,ndm1,num1,nmd2,nmu2,ndm2,num2,nmd,nmu,ndm,num,a,iuv,facint,dzvol,ind) + !$omp private ( dvol,nmd,nmu,ndm,num,a,iuv,facint,dzvol,ind,iwm,qnmd,qnmu,qndm,qnum) !$omp do schedule ( dynamic, 256 ) !$acc kernels present( kcs, zs, zb, z_volume, zsmax, zsm, & !$acc subgrid_z_zmin, subgrid_z_zmax, subgrid_z_dep, subgrid_z_volmax, & - !$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 netprcp, cumprcpt, prcp, q, 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 z_index_wavemaker, wavemaker_uvmean, wavemaker_nmd, wavemaker_nmu, wavemaker_ndm, wavemaker_num, storage_volume), async(1) !$acc loop independent, private( nm ) @@ -449,7 +330,7 @@ subroutine compute_water_levels_subgrid(dt) ! dvol = 0.0 ! - if (kcs(nm)==1 .or. kcs(nm)==4) then + if (kcs(nm)==1) then ! if (crsgeo) then a = cell_area_m2(nm) @@ -466,7 +347,6 @@ 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 dvol = dvol + cumprcpt(nm)*a cumprcpt(nm) = 0.0 ! @@ -474,152 +354,94 @@ subroutine compute_water_levels_subgrid(dt) ! endif ! - if (z_flags_type(nm) == 0) then - ! - ! Regular point with four surrounding cells of the same size + nmd = z_index_uv_md(nm) + nmu = z_index_uv_mu(nm) + ndm = z_index_uv_nd(nm) + num = z_index_uv_nu(nm) + ! + if (crsgeo) then ! - nmd = z_index_uv_md1(nm) - nmu = z_index_uv_mu1(nm) - ndm = z_index_uv_nd1(nm) - num = z_index_uv_nu1(nm) + dvol = dvol + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxm(nm) ) * dt ! - if (crsgeo) then - ! - 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 - 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 - 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 - ! else ! if (use_quadtree) then - ! - ! More complicated - ! - nmd1 = z_index_uv_md1(nm) - nmu1 = z_index_uv_mu1(nm) - ndm1 = z_index_uv_nd1(nm) - num1 = z_index_uv_nu1(nm) - nmd2 = z_index_uv_md2(nm) - nmu2 = z_index_uv_mu2(nm) - ndm2 = z_index_uv_nd2(nm) - num2 = z_index_uv_nu2(nm) - ! - ! U - ! - if (nmd1>0) then - dvol = dvol + q(nmd1)*dyrm(uv_flags_iref(nmd1))*dt - endif - ! - if (nmd2>0) then - dvol = dvol + q(nmd2)*dyrm(uv_flags_iref(nmd2))*dt - endif - ! - if (nmu1>0) then - dvol = dvol - q(nmu1)*dyrm(uv_flags_iref(nmu1))*dt - endif - ! - if (nmu2>0) then - dvol = dvol - q(nmu2)*dyrm(uv_flags_iref(nmu2))*dt - endif - ! - ! V - ! - if (crsgeo) then - ! - if (ndm1>0) then - dvol = dvol + q(ndm1)*dxm(ndm1)*dt - endif - ! - if (ndm2>0) then - dvol = dvol + q(ndm2)*dxm(ndm2)*dt - endif - ! - if (num1>0) then - dvol = dvol - q(num1)*dxm(num1)*dt - endif - ! - if (num2>0) then - dvol = dvol - q(num2)*dxm(num2)*dt - endif - ! - else - ! - if (ndm1>0) then - dvol = dvol + q(ndm1)*dxrm(uv_flags_iref(ndm1))*dt - endif - ! - if (ndm2>0) then - dvol = dvol + q(ndm2)*dxrm(uv_flags_iref(ndm2))*dt - endif - ! - if (num1>0) then - dvol = dvol - q(num1)*dxrm(uv_flags_iref(num1))*dt - endif - ! - if (num2>0) then - dvol = dvol - q(num2)*dxrm(uv_flags_iref(num2))*dt - endif - ! - endif - ! + dvol = dvol + ( (q(nmd) - q(nmu))*dyrm(uv_flags_iref(nm)) + (q(ndm) - q(num))*dxrm(uv_flags_iref(nm)) ) * dt else - ! - nmd = z_index_uv_md1(nm) - nmu = z_index_uv_mu1(nm) - ndm = z_index_uv_nd1(nm) - num = z_index_uv_nu1(nm) - ! -! dvol = 0.0 - ! - ! U - ! - if (nmd>0) then - dvol = dvol + q(nmd)*dyrm(1)*dt - endif - ! - if (nmu>0) then - dvol = dvol - q(nmu)*dyrm(1)*dt - endif - ! - ! V - ! - if (crsgeo) then - ! - if (ndm>0) then - dvol = dvol + q(ndm)*dxm(ndm)*dt - endif - ! - if (num>0) then - dvol = dvol - q(num)*dxm(num)*dt - endif - ! - else - ! - if (ndm>0) then - dvol = dvol + q(ndm)*dxrm(1)*dt - endif - ! - if (num>0) then - dvol = dvol - q(num)*dxrm(1)*dt - endif - ! - endif - ! + dvol = dvol + ( (q(nmd) - q(nmu))*dy + (q(ndm) - q(num))*dx ) * dt endif ! - endif + endif + endif ! kcs==1 + ! + if (wavemaker .and. kcs(nm)==4) then + ! + ! Wave maker point (seaward of wave maker) + ! Here we use the mean flux at the location of the wave maker + ! + iwm = z_index_wavemaker(nm) + ! + if (wavemaker_nmd(iwm)>0) then + ! + ! Wave paddle on the left + ! + qnmd = wavemaker_uvmean(wavemaker_nmd(iwm)) + ! + else + ! + qnmd = q(z_index_uv_md(nm)) + ! + endif ! + if (wavemaker_nmu(iwm)>0) then + ! + ! Wave paddle on the right + ! + qnmu = wavemaker_uvmean(wavemaker_nmu(iwm)) + ! + else + ! + qnmu = q(z_index_uv_mu(nm)) + ! + endif + ! + if (wavemaker_ndm(iwm)>0) then + ! + ! Wave paddle below + ! + qndm = wavemaker_uvmean(wavemaker_ndm(iwm)) + ! + else + ! + qndm = q(z_index_uv_nd(nm)) + ! + endif + ! + if (wavemaker_num(iwm)>0) then + ! + ! Wave paddle above + ! + qnum = wavemaker_uvmean(wavemaker_num(iwm)) + ! + else + ! + qnum = q(z_index_uv_nu(nm)) + ! + endif + ! + if (use_quadtree) then + dvol = dvol + ( (qnmd - qnmu)*dyrm(uv_flags_iref(nm)) + (qndm - qnum)*dxrm(uv_flags_iref(nm)) ) * dt + else + dvol = dvol + ( (qnmd - qnmu)*dy + (qndm - qnum)*dx ) * dt + endif + ! + endif + ! + ! We got the volume change dvol in each active cell + ! Now update the volume and compute new water level + ! + if (kcs(nm) == 1 .or. kcs(nm) == 4) then + ! if (use_storage_volume) then ! if (storage_volume(nm)>1.0e-6) then @@ -645,6 +467,8 @@ subroutine compute_water_levels_subgrid(dt) z_volume(nm) = z_volume(nm) + dvol ! endif + ! + ! Obtain new water level from subgrid tables ! if (z_volume(nm)>=subgrid_z_volmax(nm) - 1.0e-6) then ! @@ -653,10 +477,14 @@ subroutine compute_water_levels_subgrid(dt) zs(nm) = max(subgrid_z_zmax(nm), -20.0) + (z_volume(nm) - subgrid_z_volmax(nm))/a ! elseif (z_volume(nm)<=1.0e-6) then + ! + ! No water in this cell. Set zs to z_zmin. ! zs(nm) = max(subgrid_z_zmin(nm), -20.0) ! else + ! + ! Interpolation from subgrid tables needed. ! dzvol = subgrid_z_volmax(nm) / (subgrid_nbins - 1) iuv = min(int(z_volume(nm)/dzvol) + 1, subgrid_nbins - 1) @@ -665,21 +493,25 @@ subroutine compute_water_levels_subgrid(dt) ! endif ! - ! No continuity update but keeping track of variables - ! zsmax used by default, therefore keep in standard continuity loop: - if (store_maximum_waterlevel) then - ! - zsmax(nm) = max(zsmax(nm), zs(nm)) - ! - endif - ! endif - ! - if (wavemaker) then + ! + if (snapwave) then + ! + ! Time-averaged water level used for SnapWave using exponential filter + ! + ! Would double exponential filtering be better? ! zsm(nm) = factime*zs(nm) + (1.0 - factime)*zsm(nm) ! - endif + endif + ! + ! No continuity update but keeping track of variables + ! zsmax used by default, therefore keep in standard continuity loop: + if (store_maximum_waterlevel) then + ! + zsmax(nm) = max(zsmax(nm), zs(nm)) + ! + endif ! enddo !$omp end do @@ -708,14 +540,12 @@ subroutine compute_store_variables(dt) real*4 :: qvz real*4 :: qz real*4 :: uvz - ! - ! when quadtree added add back in $omp private: nmd1,nmu1,ndm1,num1,nmd2,nmu2,ndm2,num2, ! !$omp parallel & - !$omp private (nmd,nmu,ndm,num,quz,qvz,qz,uvz ) + !$omp private ( nmd, nmu, ndm, num, quz, qvz, qz, uvz ) !$omp do schedule ( dynamic, 256 ) - !$acc kernels present( kcs, zs, zb, subgrid_z_zmin, q, z_flags_type, & - !$acc z_index_uv_md1, z_index_uv_nd1, z_index_uv_mu1, z_index_uv_nu1), async(1) + !$acc kernels present( kcs, zs, zb, subgrid_z_zmin, q, vmax, qmax, twet, & + !$acc z_index_uv_md1, z_index_uv_nd1, z_index_uv_mu1, z_index_uv_nu), async(1) !$acc loop independent, private( nm ) do nm = 1, np ! @@ -726,33 +556,24 @@ subroutine compute_store_variables(dt) ! if (kcs(nm)==1 .or. kcs(nm)==4) then ! TL: kcs(nm)==4 also correct for regular? ! - if (z_flags_type(nm) == 0) then - ! - ! Regular point with four surrounding cells of the same size - ! - nmd = z_index_uv_md1(nm) - nmu = z_index_uv_mu1(nm) - ndm = z_index_uv_nd1(nm) - num = z_index_uv_nu1(nm) - ! - if (store_maximum_velocity) then - quz = (uv(nmd) + uv(nmu)) / 2 - qvz = (uv(ndm) + uv(num)) / 2 - uvz = sqrt(quz**2 + qvz**2) - endif - ! - if (store_maximum_flux) then - quz = (q(nmd) + q(nmu)) / 2 - qvz = (q(ndm) + q(num)) / 2 - qz = sqrt(quz**2 + qvz**2) - endif - ! - else - ! - !if (use_quadtree) then - ! ! TODO - !endif - endif + ! Regular point with four surrounding cells of the same size + ! + nmd = z_index_uv_md(nm) + nmu = z_index_uv_mu(nm) + ndm = z_index_uv_nd(nm) + num = z_index_uv_nu(nm) + ! + if (store_maximum_velocity) then + quz = (uv(nmd) + uv(nmu)) / 2 + qvz = (uv(ndm) + uv(num)) / 2 + uvz = sqrt(quz**2 + qvz**2) + endif + ! + if (store_maximum_flux) then + quz = (q(nmd) + q(nmu)) / 2 + qvz = (q(ndm) + q(num)) / 2 + qz = sqrt(quz**2 + qvz**2) + endif ! ! No continuity update but keeping track of variables ! 1. store vmax diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 92cc4b20..f515d06b 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -129,6 +129,7 @@ module sfincs_data character*256 :: netamuamvfile character*256 :: netampfile character*256 :: netamprfile + character*256 :: netspwfile character*256 :: scsfile character*256 :: smaxfile character*256 :: sefffile @@ -157,6 +158,7 @@ module sfincs_data character*3 :: inftype integer :: epsg character*15 :: epsg_code + integer :: nc_deflate_level ! logical :: waves logical :: wind @@ -205,8 +207,11 @@ module sfincs_data logical :: use_quadtree logical :: interpolate_zst logical :: advection + logical :: thetasmoothing logical :: fixed_output_intervals logical :: use_storage_volume + logical :: output_irregular_grid + logical :: use_spw_precip !!! !!! sfincs_input.f90 switches integer storevelmax @@ -219,33 +224,18 @@ module sfincs_data integer storemeteo integer storehsubgrid integer wrttimeoutput - integer idebug - integer iradstr - integer igeo - integer icoriolis - integer iamprblock - integer iglobal - integer itsunamitime - integer ispinupmeteo - integer isnapwave - integer iwindmax - integer ioutfixed - integer iadvection - integer istorefw - integer istorewavdir - integer imanning2d - integer iviscosity - integer isubgrid - integer iwavemaker - integer iwavemaker_spectrum - integer ispwprecip !!! !!! Static data !!! integer*4 :: np integer*4 :: npuv + integer*4 :: ncuv integer*4 :: nkcuv2 ! + ! Temp for reading ascii depfile in initialize_bathymetry(): + ! + integer*1, dimension(:,:), allocatable :: kcsg + ! ! Internal wave maker ! integer*4 :: nkcs4 @@ -254,10 +244,8 @@ module sfincs_data ! Indices ! integer*4, dimension(:), allocatable :: nmindbnd -! integer*4, dimension(:,:), allocatable :: z_index - integer*4, dimension(:), allocatable :: z_index_z_n - integer*4, dimension(:), allocatable :: z_index_z_m -! integer*4, dimension(:), allocatable :: z_index_z_nm + integer*4, dimension(:), allocatable :: z_index_z_n + integer*4, dimension(:), allocatable :: z_index_z_m integer*4, dimension(:), allocatable :: z_index_uv_md1 integer*4, dimension(:), allocatable :: z_index_uv_md2 integer*4, dimension(:), allocatable :: z_index_uv_mu1 @@ -266,6 +254,10 @@ module sfincs_data integer*4, dimension(:), allocatable :: z_index_uv_nd2 integer*4, dimension(:), allocatable :: z_index_uv_nu1 integer*4, dimension(:), allocatable :: z_index_uv_nu2 + integer*4, dimension(:), allocatable :: z_index_uv_md + integer*4, dimension(:), allocatable :: z_index_uv_mu + integer*4, dimension(:), allocatable :: z_index_uv_nd + integer*4, dimension(:), allocatable :: z_index_uv_nu ! integer*4, dimension(:), allocatable :: uv_index_z_nm integer*4, dimension(:), allocatable :: uv_index_z_nmu @@ -277,21 +269,20 @@ module sfincs_data integer*4, dimension(:), allocatable :: uv_index_v_nm integer*4, dimension(:), allocatable :: uv_index_v_ndmu integer*4, dimension(:), allocatable :: uv_index_v_nmu + integer*4, dimension(:), allocatable :: cuv_index_uv + integer*4, dimension(:), allocatable :: cuv_index_uv1 + integer*4, dimension(:), allocatable :: cuv_index_uv2 ! ! Flags ! integer*1, dimension(:), allocatable :: z_flags_iref - integer*1, dimension(:), allocatable :: z_flags_type ! integer*1, dimension(:), allocatable :: uv_flags_iref integer*1, dimension(:), allocatable :: uv_flags_type integer*1, dimension(:), allocatable :: uv_flags_dir - integer*1, dimension(:), allocatable :: uv_flags_adv - integer*1, dimension(:), allocatable :: uv_flags_vis ! integer*1, dimension(:), allocatable :: kcs integer*1, dimension(:), allocatable :: kcuv - integer*1, dimension(:), allocatable :: kfuv integer*1, dimension(:), allocatable :: scs_rain ! logic if previous time step was raining ! ! Quadtree @@ -314,10 +305,6 @@ module sfincs_data real*4, dimension(:), allocatable :: dyrinvc real*4, dimension(:), allocatable :: cell_area ! -! integer*4, dimension(:), allocatable :: nr_points_per_level -! integer*4, dimension(:,:), allocatable :: nm_indices_per_level -! integer*4, dimension(:,:), allocatable :: cell_indices_per_level - ! ! Cell sizes ! real*4, dimension(:), allocatable :: dxm @@ -332,11 +319,6 @@ module sfincs_data ! ! UV-points ! -! integer*4, dimension(:,:), allocatable :: uv_index -! integer*4, dimension(:), allocatable :: uv_index_nm -! integer*4, dimension(:), allocatable :: uv_index_nmu -! integer*1, dimension(:,:), allocatable :: uv_flags - ! real*4, dimension(:), allocatable :: zb real*4, dimension(:), allocatable :: zbuv real*4, dimension(:), allocatable :: zbuvmx @@ -422,8 +404,7 @@ module sfincs_data real*4, dimension(:), allocatable :: wmf_hm0_ig_t real*4, dimension(:), allocatable :: wmf_tp_ig_t real*4, dimension(:), allocatable :: wmf_setup_t - ! - + ! ! integer*4 :: wavemaker_nr_cross ! integer*4 :: wavemaker_nr_along ! real*4 :: wavemaker_dx_cross @@ -466,10 +447,6 @@ module sfincs_data ! integer*4, dimension(:), allocatable :: wavemaker_idir_v ! real*4, dimension(:), allocatable :: wavemaker_qym ! real*4, dimension(:), allocatable :: wavemaker_angfac_v - ! - ! General grid - ! -! real*4, dimension(:,:), allocatable, target :: xg, yg, xz, yz ! ! Sub-grid ! @@ -614,9 +591,11 @@ module sfincs_data integer :: spw_ncols integer :: spw_nquant real*4 :: spw_radius + real*4, dimension(:), allocatable :: spw_radia real*4, dimension(:), allocatable :: spw_times real*4, dimension(:), allocatable :: spw_xe real*4, dimension(:), allocatable :: spw_ye + real*4, dimension(:), allocatable :: spw_pressure_eye real*4, dimension(:,:,:), allocatable :: spw_vmag real*4, dimension(:,:,:), allocatable :: spw_vdir real*4, dimension(:,:,:), allocatable :: spw_wu @@ -919,7 +898,7 @@ subroutine finalize_parameters() if(allocated(prcp0)) deallocate(prcp0) if(allocated(prcp1)) deallocate(prcp1) ! - if(allocated(kfuv)) deallocate(kfuv) +! if(allocated(kfuv)) deallocate(kfuv) ! ! Grid boundary points ! diff --git a/source/src/sfincs_date.f90 b/source/src/sfincs_date.f90 index 329d31eb..d712b365 100644 --- a/source/src/sfincs_date.f90 +++ b/source/src/sfincs_date.f90 @@ -259,6 +259,42 @@ function convert_fewsdate(timein, timent, treftimefews, trefstr) result (timeout timeout = real((int(timein) * 60) + dtsec) ! time from fews is in minutes, then correct for use wrt sfincs treftime ! end function + ! + function convert_spw_nc_date(timein, timent, treftimefews, trefstr) result (timeout) + ! + implicit none + ! + character*41 :: treftimefews + character*15 :: trefstr + ! + integer ijul1, ijul2, itb, timent + integer yyyy1,mm1,dd1,hh1,mn1,ss1,yyyy2,mm2,dd2,hh2,mn2,ss2, tmp + integer*8 sec1,sec2 + real*8 dtsec + ! + + real*4, dimension(:), allocatable :: timein + real*4, dimension(:), allocatable :: timeout + ! + read(trefstr,'(I4,2I2,1X,3I2)')yyyy1,mm1,dd1,hh1,mn1,ss1 + ! + ! netcdf time input: "days since 1970-01-01 00:00:00Z" + ! yyyy - mm - dd HH : MM : SS + read(treftimefews(12:), '(I4,A1,I2,A1,I2,A1,I2,A1,I2,A1,I2)') yyyy2, tmp, mm2, tmp, dd2, tmp, hh2, tmp, mn2, tmp, ss2 + ! + ijul1 = julian_date (yyyy1,mm1,dd1) + ijul2 = julian_date (yyyy2,mm2,dd2) + ! + sec1 = hh1*3600 + mn1*60 + ss1 + sec2 = hh2*3600 + mn2*60 + ss2 + dtsec = (ijul2 - ijul1)*86400.0 + sec2 - sec1 + ! + allocate(timeout(timent)) + ! + timeout = (timein * 86400) + dtsec ! time from spiderweb is in days + sec1 =0 + ! + end function ! function time_to_string(t_sec, tref_string) result (date_string) ! diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90 index 57b94005..c2dea2de 100644 --- a/source/src/sfincs_domain.f90 +++ b/source/src/sfincs_domain.f90 @@ -4,60 +4,46 @@ module sfincs_domain subroutine initialize_domain() ! - ! Initialize SFINCS domain (indices, flags and depths) - ! use sfincs_data use quadtree ! implicit none ! - integer ip, n, m, nm, nmu, num, ib, iref, npu, npv, ipuv, ipu, ipv, iud, iuu, ndm, nmd, j, iq, nmx, ibin, ind - integer ikcuv2 - integer npac, idummy - integer npuvq, npuvs, npzq + call initialize_processes() ! - ! Temporary arrays + call initialize_mesh() ! - integer*1, dimension(:,:), allocatable :: kcsg + call initialize_bathymetry() ! - real*4, dimension(:,:), allocatable :: zbg - real*4, dimension(:), allocatable :: rtmp - real*4, dimension(:), allocatable :: rtmpz - real*4, dimension(:), allocatable :: rtmpuv - real*4, dimension(:), allocatable :: rghfield + call initialize_boundaries() ! - integer*1, dimension(:), allocatable :: z_index_md - integer*4, dimension(:), allocatable :: z_index_md1 - integer*4, dimension(:), allocatable :: z_index_md2 - integer*1, dimension(:), allocatable :: z_index_mu - integer*4, dimension(:), allocatable :: z_index_mu1 - integer*4, dimension(:), allocatable :: z_index_mu2 - integer*1, dimension(:), allocatable :: z_index_nd - integer*4, dimension(:), allocatable :: z_index_nd1 - integer*4, dimension(:), allocatable :: z_index_nd2 - integer*1, dimension(:), allocatable :: z_index_nu - integer*4, dimension(:), allocatable :: z_index_nu1 - integer*4, dimension(:), allocatable :: z_index_nu2 + call initialize_roughness() ! - integer*1, dimension(:), allocatable :: msk - ! - integer*1, dimension(:), allocatable :: u_iref - integer*1, dimension(:), allocatable :: u_type - integer*4, dimension(:), allocatable :: u_index_nm - integer*4, dimension(:), allocatable :: u_index_nmu - integer*1, dimension(:), allocatable :: v_iref - integer*1, dimension(:), allocatable :: v_type - integer*4, dimension(:), allocatable :: v_index_nm - integer*4, dimension(:), allocatable :: v_index_num + call initialize_infiltration() ! - integer*4, dimension(:), allocatable :: uv_index_qt_in_sf + call initialize_storage_volume() ! - ! For 'old' input + call initialize_hydro() ! - integer*4, dimension(:), allocatable :: indices + if (quadtree_nr_levels==1 .and. use_quadtree) then + ! + ! Only one refinement level found in quadtree file. Reverting to use_quadtree is false. + ! This means netcdf output will be written to a regular grid. + ! + use_quadtree = .false. + ! + endif ! - real*4 :: ylat - real*4 :: dxymin + end subroutine + + + subroutine initialize_processes() + ! + ! Initialize physical processes + ! + use sfincs_data + ! + implicit none ! ! Set some flags ! @@ -79,7 +65,7 @@ subroutine initialize_domain() meteo3d = .false. include_boundaries = .false. ! - if (spwfile(1:4) /= 'none' .or. wndfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none') then + if (spwfile(1:4) /= 'none' .or. wndfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then ! wind = .true. write(*,*)'Turning on process: Wind' @@ -90,7 +76,7 @@ subroutine initialize_domain() endif endif ! - if ((ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. spwfile(1:4) /= 'none') .and. baro==1) then + if ((ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. spwfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') .and. baro==1) then patmos = .true. write(*,*)'Turning on process: Atmospheric pressure' endif @@ -102,13 +88,13 @@ subroutine initialize_domain() ! ! Check for time-spatially varying meteo ! - if (spwfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. amprfile(1:4) /= 'none' .or. ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. netamprfile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none') then + if (spwfile(1:4) /= 'none' .or. amufile(1:4) /= 'none' .or. amprfile(1:4) /= 'none' .or. ampfile(1:4) /= 'none' .or. netampfile(1:4) /= 'none' .or. netamprfile(1:4) /= 'none' .or. netamuamvfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then meteo3d = .true. write(*,*)'Turning on flag: meteo3d' endif ! wind_reduction = .false. - if (spwfile(1:4) /= 'none') then + if (spwfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then if (z0lfile(1:4) /= 'none') then wind_reduction = .true. write(*,*)'Turning on process: wind reduction over land' @@ -121,11 +107,57 @@ subroutine initialize_domain() store_meteo = .false. endif ! + end subroutine + + + subroutine initialize_mesh() + ! + ! Initialize SFINCS domain (indices, flags and neighbors) + ! + use sfincs_data + use quadtree + ! + implicit none + ! + integer ip, n, m, nm, nmu, num, ib, iref, npu, npv, ipuv, ipu, ipv, iud, iuu, ndm, nmd, j, iq, nmx, ind + integer ndmu, nmdu, numu + integer ikcuv2 + integer npac, idummy + integer nrefuv, irefuv, npuvtotal, icuv + ! + ! Temporary arrays + ! + integer*1, dimension(:), allocatable :: z_flags_md + integer*4, dimension(:), allocatable :: z_index_z_md1 + integer*4, dimension(:), allocatable :: z_index_z_md2 + integer*1, dimension(:), allocatable :: z_flags_mu + integer*4, dimension(:), allocatable :: z_index_z_mu1 + integer*4, dimension(:), allocatable :: z_index_z_mu2 + integer*1, dimension(:), allocatable :: z_flags_nd + integer*4, dimension(:), allocatable :: z_index_z_nd1 + integer*4, dimension(:), allocatable :: z_index_z_nd2 + integer*1, dimension(:), allocatable :: z_flags_nu + integer*4, dimension(:), allocatable :: z_index_z_nu1 + integer*4, dimension(:), allocatable :: z_index_z_nu2 + integer*4, dimension(:), allocatable :: z_index_cuv_md + integer*4, dimension(:), allocatable :: z_index_cuv_mu + integer*4, dimension(:), allocatable :: z_index_cuv_nd + integer*4, dimension(:), allocatable :: z_index_cuv_nu + ! + integer*1, dimension(:), allocatable :: msk + ! + ! For 'old' input + ! + integer*4, dimension(:), allocatable :: indices + ! + real*4 :: ylat + real*4 :: dxymin + ! ! READ MESH ! ! 2 options: ! - ! 1) quadtree grid + ! 1) quadtree grid (netcdf file) ! 2) 'old' inputs with index file ! if (use_quadtree) then @@ -183,7 +215,7 @@ subroutine initialize_domain() ! ! Read index file (first time, only read number of active points) ! - write(*,*)'Reading ', trim(indexfile), ' ...' + write(*,*)'Reading index file : ', trim(indexfile), ' ...' open(unit = 500, file = trim(indexfile), form = 'unformatted', access = 'stream') read(500)np allocate(indices(np)) @@ -217,10 +249,6 @@ subroutine initialize_domain() allocate(index_sfincs_in_quadtree(quadtree_nr_points)) index_sfincs_in_quadtree = 0 ! - ! Read mask file (must have same number of cells as quadtree file !) - ! - write(*,*)'Reading ', trim(mskfile), ' ...' - ! if (inputtype=='asc') then ! ip = 0 @@ -245,6 +273,10 @@ subroutine initialize_domain() ! ! Mask in separate file ! + ! Read mask file (must have same number of cells as quadtree file !) + ! + write(*,*)'Reading mask file : ', trim(mskfile), ' ...' + ! open(unit = 500, file = trim(mskfile), form = 'unformatted', access = 'stream') read(500)msk close(500) @@ -277,44 +309,41 @@ subroutine initialize_domain() index_sfincs_in_quadtree = 0 ! allocate(z_flags_iref(np)) - allocate(z_flags_type(np)) ! - z_flags_type = 0 z_flags_iref = 0 ! ! Temporary index arrays ! - allocate(z_index_md(np)) - allocate(z_index_mu(np)) - allocate(z_index_nd(np)) - allocate(z_index_nu(np)) - allocate(z_index_md1(np)) - allocate(z_index_md2(np)) - allocate(z_index_mu1(np)) - allocate(z_index_mu2(np)) - allocate(z_index_nd1(np)) - allocate(z_index_nd2(np)) - allocate(z_index_nu1(np)) - allocate(z_index_nu2(np)) + allocate(z_flags_md(np)) + allocate(z_flags_mu(np)) + allocate(z_flags_nd(np)) + allocate(z_flags_nu(np)) + allocate(z_index_z_md1(np)) + allocate(z_index_z_md2(np)) + allocate(z_index_z_mu1(np)) + allocate(z_index_z_mu2(np)) + allocate(z_index_z_nd1(np)) + allocate(z_index_z_nd2(np)) + allocate(z_index_z_nu1(np)) + allocate(z_index_z_nu2(np)) ! allocate(z_index_z_n(np)) allocate(z_index_z_m(np)) ! - z_index_md = 0 - z_index_mu = 0 - z_index_nd = 0 - z_index_nu = 0 - z_index_md1 = 0 - z_index_md2 = 0 - z_index_mu1 = 0 - z_index_mu2 = 0 - z_index_nd1 = 0 - z_index_nd2 = 0 - z_index_nu1 = 0 - z_index_nu2 = 0 + z_flags_md = 0 + z_flags_mu = 0 + z_flags_nd = 0 + z_flags_nu = 0 + z_index_z_md1 = 0 + z_index_z_md2 = 0 + z_index_z_mu1 = 0 + z_index_z_mu2 = 0 + z_index_z_nd1 = 0 + z_index_z_nd2 = 0 + z_index_z_nu1 = 0 + z_index_z_nu2 = 0 z_index_z_n = 0 z_index_z_m = 0 -! z_index_z_nm = 0 ! ! Permanent index arrays ! @@ -327,16 +356,30 @@ subroutine initialize_domain() allocate(z_index_uv_nu1(np)) allocate(z_index_uv_nu2(np)) ! - z_index_uv_md1 = 0 - z_index_uv_md2 = 0 - z_index_uv_mu1 = 0 - z_index_uv_mu2 = 0 - z_index_uv_nd1 = 0 - z_index_uv_nd2 = 0 - z_index_uv_nu1 = 0 - z_index_uv_nu2 = 0 - z_index_z_n = 0 - z_index_z_m = 0 + allocate(z_index_uv_md(np)) + allocate(z_index_uv_mu(np)) + allocate(z_index_uv_nd(np)) + allocate(z_index_uv_nu(np)) + ! + allocate(z_index_cuv_md(np)) + allocate(z_index_cuv_mu(np)) + allocate(z_index_cuv_nd(np)) + allocate(z_index_cuv_nu(np)) + ! + z_index_uv_md1 = 0 + z_index_uv_md2 = 0 + z_index_uv_mu1 = 0 + z_index_uv_mu2 = 0 + z_index_uv_nd1 = 0 + z_index_uv_nd2 = 0 + z_index_uv_nu1 = 0 + z_index_uv_nu2 = 0 + z_index_uv_md = 0 + z_index_uv_mu = 0 + z_index_uv_nd = 0 + z_index_uv_nu = 0 + z_index_z_n = 0 + z_index_z_m = 0 ! ! Re-mapping indices ! @@ -365,40 +408,40 @@ subroutine initialize_domain() z_index_z_n(ip) = quadtree_n(iq) z_index_z_m(ip) = quadtree_m(iq) ! - z_index_nu(ip) = quadtree_nu(iq) - z_index_nd(ip) = quadtree_nd(iq) - z_index_mu(ip) = quadtree_mu(iq) - z_index_md(ip) = quadtree_md(iq) + z_flags_nu(ip) = quadtree_nu(iq) + z_flags_nd(ip) = quadtree_nd(iq) + z_flags_mu(ip) = quadtree_mu(iq) + z_flags_md(ip) = quadtree_md(iq) ! ! Neighbors ! if (quadtree_nu1(iq)>0) then - z_index_nu1(ip) = index_sfincs_in_quadtree(quadtree_nu1(iq)) + z_index_z_nu1(ip) = index_sfincs_in_quadtree(quadtree_nu1(iq)) endif if (quadtree_nu2(iq)>0) then - z_index_nu2(ip) = index_sfincs_in_quadtree(quadtree_nu2(iq)) + z_index_z_nu2(ip) = index_sfincs_in_quadtree(quadtree_nu2(iq)) endif if (quadtree_nd1(iq)>0) then - z_index_nd1(ip) = index_sfincs_in_quadtree(quadtree_nd1(iq)) + z_index_z_nd1(ip) = index_sfincs_in_quadtree(quadtree_nd1(iq)) endif if (quadtree_nd2(iq)>0) then - z_index_nd2(ip) = index_sfincs_in_quadtree(quadtree_nd2(iq)) + z_index_z_nd2(ip) = index_sfincs_in_quadtree(quadtree_nd2(iq)) endif if (quadtree_mu1(iq)>0) then - z_index_mu1(ip) = index_sfincs_in_quadtree(quadtree_mu1(iq)) + z_index_z_mu1(ip) = index_sfincs_in_quadtree(quadtree_mu1(iq)) endif if (quadtree_mu2(iq)>0) then - z_index_mu2(ip) = index_sfincs_in_quadtree(quadtree_mu2(iq)) + z_index_z_mu2(ip) = index_sfincs_in_quadtree(quadtree_mu2(iq)) endif if (quadtree_md1(iq)>0) then - z_index_md1(ip) = index_sfincs_in_quadtree(quadtree_md1(iq)) + z_index_z_md1(ip) = index_sfincs_in_quadtree(quadtree_md1(iq)) endif if (quadtree_md2(iq)>0) then - z_index_md2(ip) = index_sfincs_in_quadtree(quadtree_md2(iq)) + z_index_z_md2(ip) = index_sfincs_in_quadtree(quadtree_md2(iq)) endif ! enddo - ! + ! ! Loop through quadtree to count the number of u and v points ! npu = 0 @@ -408,19 +451,19 @@ subroutine initialize_domain() ! ! Right ! - if (z_index_mu(nm)<1) then + if (z_flags_mu(nm)<1) then ! - if (z_index_mu1(nm)>0) then + if (z_index_z_mu1(nm)>0) then npu = npu + 1 endif ! else ! - if (z_index_mu1(nm)>0) then + if (z_index_z_mu1(nm)>0) then npu = npu + 1 endif ! - if (z_index_mu2(nm)>0) then + if (z_index_z_mu2(nm)>0) then npu = npu + 1 endif ! @@ -428,19 +471,19 @@ subroutine initialize_domain() ! ! Above ! - if (z_index_nu(nm)<1) then + if (z_flags_nu(nm)<1) then ! - if (z_index_nu1(nm)>0) then + if (z_index_z_nu1(nm)>0) then npv = npv + 1 endif ! else ! - if (z_index_nu1(nm)>0) then + if (z_index_z_nu1(nm)>0) then npv = npv + 1 endif ! - if (z_index_nu2(nm)>0) then + if (z_index_z_nu2(nm)>0) then npv = npv + 1 endif ! @@ -450,6 +493,129 @@ subroutine initialize_domain() ! npuv = npu + npv ! + ncuv = 0 + ! + if (quadtree_nr_levels>1) then + ! + ! Loop through all the points again and find 'combined' uv points at refinement boundaries + ! + do nm = 1, np + ! + ! Left + ! + if (z_flags_md(nm) == 1) then + ! + ! Finer on the left + ! + ncuv = ncuv + 1 + ! + endif + ! + ! Right + ! + if (z_flags_mu(nm) == 1) then + ! + ! Finer on the right + ! + ncuv = ncuv + 1 + ! + endif + ! + ! Below + ! + if (z_flags_nd(nm) == 1) then + ! + ! Finer below + ! + ncuv = ncuv + 1 + ! + endif + ! + ! Above + ! + if (z_flags_nu(nm) == 1) then + ! + ! Finer above + ! + ncuv = ncuv + 1 + ! + endif + ! + enddo + ! + write(*,*)'Number of refinement transitions : ', ncuv + ! + endif + ! + npuvtotal = npuv + ncuv ! total number of points in uv arrays (including combined uv points) + ! + ! Set neighbors of cells to dummy uv points + ! + z_index_uv_md = npuv + ncuv + 1 + z_index_uv_mu = npuv + ncuv + 1 + z_index_uv_nd = npuv + ncuv + 1 + z_index_uv_nu = npuv + ncuv + 1 + ! + ! Now allocate and fill arrays with combined uv points + ! + allocate(cuv_index_uv(ncuv)) ! index in uv array of combined point + allocate(cuv_index_uv1(ncuv)) ! index in uv array of first uv point + allocate(cuv_index_uv2(ncuv)) ! index in uv array of second uv point + ! + irefuv = 0 + ! + do nm = 1, np + ! + ! Left + ! + if (z_flags_md(nm) == 1) then + ! + ! Finer on the left + ! + irefuv = irefuv + 1 + cuv_index_uv(irefuv) = npuv + irefuv ! used to update uv and q in momentum + z_index_cuv_md(nm) = irefuv ! temporary + ! + endif + ! + ! Right + ! + if (z_flags_mu(nm) == 1) then + ! + ! Finer on the right + ! + irefuv = irefuv + 1 + cuv_index_uv(irefuv) = npuv + irefuv + z_index_cuv_mu(nm) = irefuv ! temporary + ! + endif + ! + ! Below + ! + if (z_flags_nd(nm) == 1) then + ! + ! Finer below + ! + irefuv = irefuv + 1 + cuv_index_uv(irefuv) = npuv + irefuv ! needed in momentum + z_index_cuv_nd(nm) = irefuv ! temporary + ! + endif + ! + ! Above + ! + if (z_flags_nu(nm) == 1) then + ! + ! Finer above + ! + irefuv = irefuv + 1 + cuv_index_uv(irefuv) = npuv + irefuv + z_index_cuv_nu(nm) = irefuv + ! + endif + ! + enddo + ! ! UV-points ! allocate(kcuv(npuv)) @@ -467,17 +633,10 @@ subroutine initialize_domain() uv_index_z_nm = 0 uv_index_z_nmu = 0 ! -! if (advection .or. coriolis) then - allocate(uv_flags_adv(npuv)) - uv_flags_adv = 0 -! endif + ! If advection, viscosity, coriolis or thetasmoothing, allocate extra arrays for 8 uv neighbor indices ! -! if (viscosity) then - allocate(uv_flags_vis(npuv)) - uv_flags_vis = 0 -! endif - ! -! if (advection) then + if (advection .or. coriolis .or. viscosity .or. thetasmoothing) then + ! allocate(uv_index_u_nmd(npuv)) allocate(uv_index_u_nmu(npuv)) allocate(uv_index_u_num(npuv)) @@ -486,6 +645,7 @@ subroutine initialize_domain() allocate(uv_index_v_nm(npuv)) allocate(uv_index_v_ndmu(npuv)) allocate(uv_index_v_nmu(npuv)) + ! uv_index_u_nmd = 0 uv_index_u_nmu = 0 uv_index_u_num = 0 @@ -494,619 +654,528 @@ subroutine initialize_domain() uv_index_v_nm = 0 uv_index_v_ndmu = 0 uv_index_v_nmu = 0 -! endif - ! -! if (viscosity .and. .not.advection) then -! allocate(uv_index_u_nmd(npuv)) -! allocate(uv_index_u_nmu(npuv)) -! allocate(uv_index_u_num(npuv)) -! allocate(uv_index_u_ndm(npuv)) -! uv_index_u_nmd = 0 -! uv_index_u_nmu = 0 -! uv_index_u_num = 0 -! uv_index_u_ndm = 0 -! endif - ! -! if (coriolis .and. .not.advection) then -! allocate(uv_index_v_nm(npuv)) -! uv_index_v_nm = 0 -! endif - ! - ip = 0 + ! + endif + ! + ip = 0 ! Number of uv points + ! + ! Loop through cells to set indices neighboring cells for each uv point + ! Also fill temporary arrays with indices of all 8 uv points for each cell + ! Also set some flags for each uv point ! do nm = 1, np ! ! Right ! - if (z_index_mu(nm)==0) then + if (z_flags_mu(nm)==0) then ! ! Same level ! - if (z_index_mu1(nm)>0) then + if (z_index_z_mu1(nm)>0) then ! and there is a neighbor + ! + ip = ip + 1 + ! + nmu = z_index_z_mu1(nm) + ! + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = nmu + z_index_uv_mu(nm) = ip + z_index_uv_md(nmu) = ip + z_index_uv_mu1(nm) = ip + z_index_uv_md1(nmu) = ip + ! + ! Flags to describe uv point ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_mu1(nm) - z_index_uv_mu1(nm) = ip - z_index_uv_md1(z_index_mu1(nm)) = ip uv_flags_iref(ip) = z_flags_iref(nm) uv_flags_dir(ip) = 0 ! u point uv_flags_type(ip) = 0 ! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection ! endif ! - elseif (z_index_mu(nm)==-1) then + elseif (z_flags_mu(nm)==-1) then ! ! Coarser to the right ! - if (z_index_mu1(nm)>0) then - ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_mu1(nm) - z_index_uv_mu1(nm) = ip - if (z_index_md1(z_index_mu1(nm)) == nm) then - z_index_uv_md1(z_index_mu1(nm)) = ip - else - z_index_uv_md2(z_index_mu1(nm)) = ip - endif - uv_flags_iref(ip) = z_flags_iref(nm) - uv_flags_dir(ip) = 0 ! u point - uv_flags_type(ip) = -1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection - ! - endif + ip = ip + 1 + nmu = z_index_z_mu1(nm) + ! + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = nmu + z_index_uv_mu(nm) = ip + z_index_uv_mu1(nm) = ip + ! + ! Set the combined uv point of the neighbor to the right (icuv is index of combined uv point in cuv array) + ! + icuv = z_index_cuv_md(nmu) + ! set indices + if (z_index_z_md1(z_index_z_mu1(nm)) == nm) then + ! Lower cell + cuv_index_uv1(icuv) = ip + z_index_uv_md1(nmu) = ip + else + ! Upper cell + cuv_index_uv2(icuv) = ip + z_index_uv_md2(nmu) = ip + endif + ! + ! Combined uv point for cell on the right + ! + z_index_uv_md(nmu) = cuv_index_uv(icuv) + ! + ! Flags to describe uv point + ! + uv_flags_iref(ip) = z_flags_iref(nm) + uv_flags_dir(ip) = 0 ! u point + uv_flags_type(ip) = -1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine ! else ! ! Finer to the right ! - if (z_index_mu1(nm)>0) then + if (z_index_z_mu1(nm)>0) then + ! + ! Lower cell + ! + ip = ip + 1 + nmu = z_index_z_mu1(nm) ! z index of lower fine cell on the right + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = nmu + z_index_uv_mu1(nm) = ip + z_index_uv_md(nmu) = ip + z_index_uv_md1(nmu) = ip + ! + icuv = z_index_cuv_mu(nm) ! (icuv is index of combined uv point in cuv array) + ! + ! Set the combined uv point of this cell + ! + z_index_uv_mu(nm) = cuv_index_uv(icuv) + ! + cuv_index_uv1(icuv) = ip + ! + ! Flags to describe uv point ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_mu1(nm) - z_index_uv_mu1(nm) = ip - z_index_uv_md1(z_index_mu1(nm)) = ip uv_flags_iref(ip) = z_flags_iref(nm) + 1 uv_flags_dir(ip) = 0 ! u point uv_flags_type(ip) = 1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection ! endif ! - if (z_index_mu2(nm)>0) then + if (z_index_z_mu2(nm)>0) then ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_mu2(nm) - z_index_uv_mu2(nm) = ip - z_index_uv_md1(z_index_mu2(nm)) = ip - uv_flags_iref(ip) = z_flags_iref(nm) + 1 - uv_flags_dir(ip) = 0 ! u point - uv_flags_type(ip) = 1 ! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection + ! Upper cell ! - endif - ! - endif - ! - ! Above - ! - if (z_index_nu(nm)==0) then - ! - ! Same level - ! - if (z_index_nu1(nm)>0) then + ip = ip + 1 + nmu = z_index_z_mu2(nm) ! z index of upper fine cell on the right + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = nmu + z_index_uv_mu2(nm) = ip + z_index_uv_md(nmu) = ip + z_index_uv_md1(nmu) = ip ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_nu1(nm) - z_index_uv_nu1(nm) = ip - z_index_uv_nd1(z_index_nu1(nm)) = ip - uv_flags_iref(ip) = z_flags_iref(nm) - uv_flags_dir(ip) = 1 ! v point - uv_flags_type(ip) = 0 ! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection + icuv = z_index_cuv_mu(nm) ! (icuv is index of combined uv point in cuv array) + ! + ! Set the combined uv point of this cell (this may already have been done in the lower cell) ! - endif - ! - elseif (z_index_nu(nm)==-1) then - ! - ! Coarser above - ! - if (z_index_nu1(nm)>0) then + z_index_uv_mu(nm) = cuv_index_uv(icuv) ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_nu1(nm) - z_index_uv_nu1(nm) = ip - if (z_index_nd1(z_index_nu1(nm)) == nm) then - z_index_uv_nd1(z_index_nu1(nm)) = ip - else - z_index_uv_nd2(z_index_nu1(nm)) = ip - endif - uv_flags_iref(ip) = z_flags_iref(nm) - uv_flags_dir(ip) = 1 ! u point - uv_flags_type(ip) = -1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection + cuv_index_uv2(icuv) = ip ! - endif - ! - else - ! - ! Finer above - ! - if (z_index_nu1(nm)>0) then + ! Flags to describe uv point ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_nu1(nm) - z_index_uv_nu1(nm) = ip - z_index_uv_nd1(z_index_nu1(nm)) = ip uv_flags_iref(ip) = z_flags_iref(nm) + 1 - uv_flags_dir(ip) = 1 ! v point + uv_flags_dir(ip) = 0 ! u point uv_flags_type(ip) = 1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection - - endif - ! - if (z_index_nu2(nm)>0) then - ! - ip = ip + 1 - uv_index_z_nm(ip) = nm - uv_index_z_nmu(ip) = z_index_nu2(nm) - z_index_uv_nu2(nm) = ip - z_index_uv_nd1(z_index_nu2(nm)) = ip - uv_flags_iref(ip) = z_flags_iref(nm) + 1 - uv_flags_dir(ip) = 1 ! v point - uv_flags_type(ip) = 1 ! -1 is fine too coarse, 0 is normal, 1 is coarse to fine - uv_flags_vis(ip) = 1 ! enable viscosity - if (advection .or. coriolis) uv_flags_adv(ip) = 1 ! enable advection ! endif ! endif ! - enddo - ! - ! And now for the uv neighbors of the uv points - ! - do ip = 1, npuv + ! Above ! - if (uv_flags_dir(ip)==0) then - ! - ! U point + if (z_flags_nu(nm)==0) then ! - nm = uv_index_z_nm(ip) - ! - ! left + ! Same level ! - if (z_index_md(nm)<1) then + if (z_index_z_nu1(nm)>0) then ! and there is a neighbor ! - ! Coarser or same level to the left + ip = ip + 1 ! - uv_index_u_nmd(ip) = z_index_uv_md1(nm) - ! - endif - ! - ! right - ! - if (z_index_mu(nm)==0) then + num = z_index_z_nu1(nm) ! - ! Same level to the right + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = num + z_index_uv_nu(nm) = ip + z_index_uv_nd(num) = ip + z_index_uv_nu1(nm) = ip + z_index_uv_nd1(num) = ip ! - nmu = z_index_mu1(nm) + ! Flags to describe uv point ! - if (nmu>0) then - if (z_index_mu(nmu)<1) then - ! - ! Not finer to the right of nmu - ! - uv_index_u_nmu(ip) = z_index_uv_mu1(nmu) - ! - endif - endif + uv_flags_iref(ip) = z_flags_iref(nm) + uv_flags_dir(ip) = 1 ! u point + uv_flags_type(ip) = 0 ! -1 is fine too coarse, 0 is normal, 1 is coarse to fine ! - endif - ! - ! below + endif ! - if (z_index_nd(nm)==0) then - ! - ! Same level below - ! - ndm = z_index_nd1(nm) - ! - if (ndm>0) then - if (z_index_mu(ndm)<1) then - ! - ! Not finer to the right of ndm - ! - uv_index_u_ndm(ip) = z_index_uv_mu1(ndm) - ! - endif - endif - ! - elseif (z_index_nd(nm)==-1) then - ! - ! Coarser below - ! - ndm = z_index_nd1(nm) - ! - if (z_index_nu2(ndm)==nm) then - if (z_index_mu(ndm)==1) then - ! - ! Finer to the right of ndm - ! - uv_index_u_ndm(ip) = z_index_uv_mu2(ndm) - ! - endif - endif - ! - endif + elseif (z_flags_nu(nm)==-1) then ! - ! 6 + ! Coarser above ! - if (z_index_nu(nm)==0) then - ! - ! Same level above - ! - num = z_index_nu1(nm) - ! - if (num>0) then - if (z_index_mu(num)<1) then - ! - ! Not finer to the right of num - ! - uv_index_u_num(ip) = z_index_uv_mu1(num) - ! - endif - endif - ! - elseif (z_index_nu(nm)==-1) then - ! - ! Coarser above - ! - num = z_index_nu1(nm) - ! - if (z_index_nd2(num)==nm) then - if (z_index_mu(num)==1) then - ! - ! Finer to the right of num - ! - uv_index_u_num(ip) = z_index_uv_mu1(num) - ! - endif - endif - ! + ip = ip + 1 + num = z_index_z_nu1(nm) + ! + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = num + z_index_uv_nu(nm) = ip + z_index_uv_nu1(nm) = ip + ! + ! Set the combined uv point of the neighbor above (icuv is index of combined uv point in cuv array) + ! + icuv = z_index_cuv_nd(num) + ! set indices + if (z_index_z_nd1(z_index_z_nu1(nm)) == nm) then + ! Lower cell + cuv_index_uv1(icuv) = ip + z_index_uv_nd1(num) = ip + else + ! Upper cell + cuv_index_uv2(icuv) = ip + z_index_uv_nd2(num) = ip endif ! - if (advection .or. coriolis) then - ! - ! 7 - ! - if (z_index_nd(nm)<1) then - ! - ! Same or coarser below - ! - uv_index_v_ndm(ip) = z_index_uv_nd1(nm) - ! - endif - ! - ! 8 - ! - if (z_index_mu(nm)==0) then - ! - ! Same on the right - ! - nmu = z_index_mu1(nm) - ! - if (nmu>0) then - ! - if (z_index_nd(nmu)<1) then - ! - ! Same or coarser below - ! - uv_index_v_ndmu(ip) = z_index_uv_nd1(nmu) - ! - endif - ! - endif - ! - endif - ! - ! 9 - ! - if (z_index_nu(nm)<1) then - ! - ! Same or coarser above - ! - uv_index_v_nm(ip) = z_index_uv_nu1(nm) - ! - endif - ! - ! 10 - ! - if (z_index_mu(nm)==0) then - ! - ! Same to the right - ! - nmu = z_index_mu1(nm) - ! - if (nmu>0) then - ! - if (z_index_nu(nmu)<1) then - ! - ! Same or coarser above - ! - uv_index_v_nmu(ip) = z_index_uv_nu1(nmu) - ! - endif - ! - endif - ! - endif - ! - endif + ! Combined uv point for cell above + ! + z_index_uv_nd(num) = cuv_index_uv(icuv) ! - else + ! Flags to describe uv point ! - ! V point + uv_flags_iref(ip) = z_flags_iref(nm) + uv_flags_dir(ip) = 1 ! v point + uv_flags_type(ip) = -1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine ! - nm = uv_index_z_nm(ip) + else ! - ! left + ! Finer above ! - if (z_index_nd(nm)<1) then + if (z_index_z_nu1(nm)>0) then ! - ! Coarser or same level to the left + ! Left cell ! - uv_index_u_nmd(ip) = z_index_uv_nd1(nm) + ip = ip + 1 + num = z_index_z_nu1(nm) ! z index of left fine cell above + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = num + z_index_uv_nd(num) = ip + z_index_uv_nu1(nm) = ip + z_index_uv_nd1(num) = ip ! - endif - ! - ! right - ! - if (z_index_nu(nm)==0) then + icuv = z_index_cuv_nu(nm) ! (icuv is index of combined uv point in cuv array) + ! + ! Set the combined uv point of this cell ! - ! Same level to the right + z_index_uv_nu(nm) = cuv_index_uv(icuv) ! - num = z_index_nu1(nm) + cuv_index_uv1(icuv) = ip ! - if (num>0) then - if (z_index_nu(num)<1) then - ! - ! Not finer to the right of nmu - ! - uv_index_u_nmu(ip) = z_index_uv_nu1(num) - ! - endif - endif + ! Flags to describe uv point ! - endif - ! - ! below + uv_flags_iref(ip) = z_flags_iref(nm) + 1 + uv_flags_dir(ip) = 1 ! v point + uv_flags_type(ip) = 1! -1 is fine too coarse, 0 is normal, 1 is coarse to fine + ! + endif ! - if (z_index_md(nm)==0) then + if (z_index_z_nu2(nm)>0) then ! - ! Same level below + ! Right cell ! - nmd = z_index_md1(nm) + ip = ip + 1 + num = z_index_z_nu2(nm) ! z index of left fine cell above + uv_index_z_nm(ip) = nm + uv_index_z_nmu(ip) = num + z_index_uv_nd(num) = ip + z_index_uv_nu2(nm) = ip + z_index_uv_nd1(num) = ip ! - if (nmd>0) then - if (z_index_nu(nmd)<1) then - ! - ! Not finer to the right of ndm - ! - uv_index_u_ndm(ip) = z_index_uv_nu1(nmd) - ! - endif - endif + icuv = z_index_cuv_nu(nm) ! (icuv is index of combined uv point in cuv array) + ! + ! Set the combined uv point of this cell ! - elseif (z_index_md(nm)==-1) then + z_index_uv_nu(nm) = cuv_index_uv(icuv) ! - ! Coarser below + cuv_index_uv2(icuv) = ip ! - nmd = z_index_md1(nm) + ! Flags to describe uv point ! - if (z_index_mu2(nmd)==nm) then - if (z_index_nu(nmd)==1) then - ! - ! Finer to the right of ndm - ! - uv_index_u_ndm(ip) = z_index_uv_nu2(nmd) - ! - endif - endif + uv_flags_iref(ip) = z_flags_iref(nm) + 1 + uv_flags_dir(ip) = 1 ! v point + uv_flags_type(ip) = 1 ! -1 is fine too coarse, 0 is normal, 1 is coarse to fine ! - endif - ! - ! 6 + endif + ! + endif + ! + enddo + ! + ! And now set the indices of the 8 uv neighbors of the uv points (u_nmu, u_nmd, u_num, u_ndm, + ! v_ndm, v_nm, v_ndmu, v_nmu) + if (advection .or. coriolis .or. viscosity .or. thetasmoothing) then + ! + do ip = 1, npuv ! - if (z_index_mu(nm)==0) then - ! - ! Same level above + if (uv_flags_dir(ip) == 0) then ! - nmu = z_index_mu1(nm) + ! U point ! - if (nmu>0) then - if (z_index_nu(nmu)<1) then - ! - ! Not finer to the right of num - ! - uv_index_u_num(ip) = z_index_uv_nu1(nmu) - ! - endif - endif + nm = uv_index_z_nm(ip) ! index of cell to the left of this uv point + nmu = uv_index_z_nmu(ip) ! index of cell to the right of this uv point ! - elseif (z_index_mu(nm)==-1) then + ! Left (u_nmd) + ! + uv_index_u_nmd(ip) = z_index_uv_md(nm) ! - ! Coarser above + ! Right (u_nmu) + ! + uv_index_u_nmu(ip) = z_index_uv_mu(nmu) ! - nmu = z_index_mu1(nm) + ! Below (u_ndm) ! - if (z_index_md2(nmu)==nm) then - if (z_index_nu(nmu)==1) then - ! - ! Finer to the right of num + if (z_flags_nd(nm)>-1) then ! The cell below is coarser or same level. Easy. + ! + ndm = z_index_z_nd1(nm) ! z-index of cell below + ! + if (ndm>0) then ! And it exists ! - uv_index_u_num(ip) = z_index_uv_nu1(nmu) + uv_index_u_ndm(ip) = z_index_uv_mu(ndm) ! - endif - endif - ! - endif - ! - if (advection .or. coriolis) then - ! - ! 7 - ! - if (z_index_md(nm)<1) then + endif ! - ! Same or coarser below - ! - uv_index_v_ndm(ip) = z_index_uv_md1(nm) + else + ! + ! Cell below is finer. Difficult. + ! + ! Check if cell ndmu is same level. If so, use z_index_uv_md of that cell. Otherwise, use z_index_uv_mu of cell below. + ! + ndm = z_index_z_nd2(nm) ! z-index of cell below ! + ! Check if cell below exists + ! + if (ndm>0) then + ! + if (z_flags_mu(ndm)==0) then ! cell ndmu is also finer level (could also be same level, but not coarser) + ! + uv_index_u_ndm(ip) = z_index_uv_mu(ndm) + ! + else ! cell ndmu is same level + ! + ndmu = z_index_z_mu1(ndm) ! z-index of cell below-right + ! + if (ndmu>0) then + ! + uv_index_u_ndm(ip) = z_index_uv_md(ndmu) + ! + endif + ! + endif + endif endif ! - ! 8 + ! Above (u_num) ! - if (z_index_nu(nm)==0) then + if (z_flags_nu(nm)>-1) then ! The cell above is coarser or same level. Easy. + ! + num = z_index_z_nu1(nm) ! z-index of cell above ! - ! Same on the right + if (num>0) then ! And it exists + ! + uv_index_u_num(ip) = z_index_uv_mu(num) + ! + endif ! - num = z_index_nu1(nm) + else + ! + ! Cell above is finer. Difficult. + ! + ! Check if cell numu is same level. If so, use z_index_uv_md of that cell. Otherwise, use z_index_uv_mu of cell above. + ! + num = z_index_z_nu2(nm) ! z-index of cell above ! + ! Check if cell above exists + ! if (num>0) then ! - if (z_index_nd(num)<1) then + if (z_flags_mu(num)==0) then ! cell numu is also finer level (could also same level, but not coarser) ! - ! Same or coarser below + uv_index_u_num(ip) = z_index_uv_mu(num) ! - uv_index_v_ndmu(ip) = z_index_uv_md1(num) + else ! cell numu is same level ! - endif - ! - endif - ! - endif - ! - ! 9 + numu = z_index_z_mu1(num) ! z-index of cell below-right + ! + if (numu>0) then + ! + uv_index_u_num(ip) = z_index_uv_md(numu) + ! + endif + ! + endif + endif + endif + ! + ! Left (v_nm) + ! + uv_index_v_nm(ip) = z_index_uv_nu(nm) + ! + ! Left-below (v_ndm) + ! + uv_index_v_ndm(ip) = z_index_uv_nd(nm) + ! + ! Right (v_nmu) ! - if (z_index_mu(nm)<1) then + uv_index_v_nmu(ip) = z_index_uv_nu(nmu) + ! + ! Right-below (v_ndmu) + ! + uv_index_v_ndmu(ip) = z_index_uv_nd(nmu) + ! + else + ! + ! V point (mirror-image of U point) + ! + nm = uv_index_z_nm(ip) ! index of cell below this uv point + num = uv_index_z_nmu(ip) ! index of cell above this uv point + ! + ! Left (u_nmd) + ! + uv_index_u_nmd(ip) = z_index_uv_nd(nm) + ! + ! Right (u_nmu) + ! + uv_index_u_nmu(ip) = z_index_uv_nu(num) + ! + ! Below (u_ndm) + ! + if (z_flags_md(nm)>-1) then ! The cell left is coarser or same level. Easy. + ! + nmd = z_index_z_md1(nm) ! z-index of cell left ! - ! Same or coarser above + if (nmd>0) then ! And it exists + ! + uv_index_u_ndm(ip) = z_index_uv_nu(nmd) + ! + endif ! - uv_index_v_nm(ip) = z_index_uv_mu1(nm) + else + ! + ! Cell below is finer. Difficult. + ! + ! Check if cell ndmu is same level. If so, use z_index_uv_md of that cell. Otherwise, use z_index_uv_mu of cell below. + ! + nmd = z_index_z_md2(nm) ! z-index of cell left ! + ! Check if cell below exists + ! + if (nmd>0) then + ! + if (z_flags_nu(nmd)==0) then ! cell ndmu is also finer level (could also be same level, but not coarser) + ! + uv_index_u_ndm(ip) = z_index_uv_nu(nmd) + ! + else ! cell nmdu is same level + ! + nmdu = z_index_z_nu1(nmd) ! z-index of cell below-right + ! + if (nmdu>0) then + ! + uv_index_u_ndm(ip) = z_index_uv_nd(nmdu) + ! + endif + ! + endif + endif endif ! - ! 10 + ! Above (u_num) ! - if (z_index_nu(nm)==0) then + if (z_flags_mu(nm)>-1) then ! The cell right coarser or same level. Easy. + ! + nmu = z_index_z_mu1(nm) ! z-index of cell right ! - ! Same to the right + if (nmu>0) then ! And it exists + ! + uv_index_u_num(ip) = z_index_uv_nu(nmu) + ! + endif ! - num = z_index_nu1(nm) + else + ! + ! Cell above is finer. Difficult. + ! + ! Check if cell numu is same level. If so, use z_index_uv_md of that cell. Otherwise, use z_index_uv_mu of cell above. + ! + nmu = z_index_z_mu2(nm) ! z-index of cell right ! - if (num>0) then + ! Check if cell above exists + ! + if (nmu>0) then ! - if (z_index_mu(num)<1) then + if (z_flags_nu(nmu)==0) then ! cell numu is also finer level (could also same level, but not coarser) ! - ! Same or coarser above + uv_index_u_num(ip) = z_index_uv_nu(nmu) ! - uv_index_v_nmu(ip) = z_index_uv_mu1(num) + else ! cell numu is same level ! - endif - ! - endif - ! - endif - ! - endif - ! - endif - ! - enddo - ! - ! Set flags for cells that do not have 4 normal neighbors (used in continuity) - ! - do nm = 1, np - ! - z_flags_type(nm) = 0 - ! - if (use_quadtree) then - ! - if (z_index_uv_md1(nm)==0 .or. z_index_uv_md2(nm)>0) then - z_flags_type(nm) = 1 - endif - ! - if (z_index_uv_mu1(nm)==0 .or. z_index_uv_mu2(nm)>0) then - z_flags_type(nm) = 1 - endif - ! - if (z_index_uv_nd1(nm)==0 .or. z_index_uv_nd2(nm)>0) then - z_flags_type(nm) = 1 - endif - ! - if (z_index_uv_nu1(nm)==0 .or. z_index_uv_nu2(nm)>0) then - z_flags_type(nm) = 1 - endif - ! - else - ! - if (z_index_uv_md1(nm)==0) then - z_flags_type(nm) = 1 - endif - ! - if (z_index_uv_mu1(nm)==0) then - z_flags_type(nm) = 1 - endif - ! - if (z_index_uv_nd1(nm)==0) then - z_flags_type(nm) = 1 - endif - ! - if (z_index_uv_nu1(nm)==0) then - z_flags_type(nm) = 1 + numu = z_index_z_nu1(nmu) ! z-index of cell below-right + ! + if (numu>0) then + ! + uv_index_u_num(ip) = z_index_uv_nd(numu) + ! + endif + ! + endif + endif + endif + ! + ! Left (v_nm) + ! + uv_index_v_nm(ip) = z_index_uv_mu(nm) + ! + ! Left-below (v_ndm) + ! + uv_index_v_ndm(ip) = z_index_uv_md(nm) + ! + ! Right (v_nmu) + ! + uv_index_v_nmu(ip) = z_index_uv_mu(num) + ! + ! Right-below (v_ndmu) + ! + uv_index_v_ndmu(ip) = z_index_uv_md(num) + ! endif ! - endif - ! - enddo -! ! -! ! Set flags to turn off viscosity in some cells -! ! -! do ip = 1, npuv -! do j = 3, 6 -! if (uv_index(j, ip)==0) then -! uv_flags(4, ip) = 0 -! endif -! enddo -! enddo -! ! -! ! Set flags to turn off advection in some cells -! ! -! if (advection .or. coriolis) then -! do ip = 1, npuv -! do j = 3, 10 -! if (uv_index(j, ip)==0) then -! uv_flags(5, ip) = 0 -! endif -! enddo -! enddo -! endif - ! - ! Okay, got all the quadtree cells including indices and flags + ! Missing uv neighbors + ! + ! For u dir points, set neighbor to this uv point itself + ! + if (uv_index_u_nmd(ip) == 0) uv_index_u_nmd(ip) = ip + if (uv_index_u_nmu(ip) == 0) uv_index_u_nmu(ip) = ip + if (uv_index_u_ndm(ip) == 0) uv_index_u_ndm(ip) = ip + if (uv_index_u_num(ip) == 0) uv_index_u_num(ip) = ip + ! + ! For v dir points, set neighbor to dummy index + ! + if (uv_index_v_ndm(ip) == 0) uv_index_v_ndm(ip) = npuvtotal + 1 + if (uv_index_v_ndmu(ip) == 0) uv_index_v_ndmu(ip) = npuvtotal + 1 + if (uv_index_v_nm(ip) == 0) uv_index_v_nm(ip) = npuvtotal + 1 + if (uv_index_v_nmu(ip) == 0) uv_index_v_nmu(ip) = npuvtotal + 1 + ! + enddo ! uv point + ! + endif ! advection, viscosity, coriolis or thetasmoothing + ! + ! Okay, got all the quadtree cells including indices, neighbors flags ! write(*,*)'Number of active z points : ', np write(*,*)'Number of active u/v points : ', npuv @@ -1117,6 +1186,27 @@ subroutine initialize_domain() enddo endif ! + ! Time to de-allocate some arrays + ! + if (allocated(z_flags_md)) deallocate(z_flags_md) + if (allocated(z_index_z_md1)) deallocate(z_index_z_md1) + if (allocated(z_index_z_md2)) deallocate(z_index_z_md2) + if (allocated(z_flags_mu)) deallocate(z_flags_mu) + if (allocated(z_index_z_mu1)) deallocate(z_index_z_mu1) + if (allocated(z_index_z_mu2)) deallocate(z_index_z_mu2) + if (allocated(z_flags_nd)) deallocate(z_flags_nd) + if (allocated(z_index_z_nd1)) deallocate(z_index_z_nd1) + if (allocated(z_index_z_nd2)) deallocate(z_index_z_nd2) + if (allocated(z_flags_nu)) deallocate(z_flags_nu) + if (allocated(z_index_z_nu1)) deallocate(z_index_z_nu1) + if (allocated(z_index_z_nu2)) deallocate(z_index_z_nu2) + if (allocated(z_index_cuv_md)) deallocate(z_index_cuv_md) + if (allocated(z_index_cuv_mu)) deallocate(z_index_cuv_mu) + if (allocated(z_index_cuv_nd)) deallocate(z_index_cuv_nd) + if (allocated(z_index_cuv_nu)) deallocate(z_index_cuv_nu) + if (allocated(msk)) deallocate(msk) + if (allocated(indices)) deallocate(indices) + ! ! MESH GEOMETRY ! ! Refinement levels @@ -1138,6 +1228,8 @@ subroutine initialize_domain() allocate(z_xz(np)) allocate(z_yz(np)) ! + write(*,*)'Computing cell centre coordinates ...' + ! do nm = 1, np ! n = z_index_z_n(nm) @@ -1271,12 +1363,38 @@ subroutine initialize_domain() ! enddo ! + ! Determine minimum and maximum time step + ! + dtmax = min(dtmax, dxymin/(sqrt(9.81*0.1))) + dtmin = alfa*dxymin/(sqrt(9.81*stopdepth)) ! If dt falls below this value, the simulation will stop + ! endif ! - ! Determine minimum and maximum time step + end subroutine + + subroutine initialize_bathymetry() + ! + use sfincs_data + use quadtree + ! + implicit none ! - dtmax = min(dtmax, dxymin/(sqrt(9.81*0.1))) - dtmin = alfa*dxymin/(sqrt(9.81*stopdepth)) ! If dt falls below this value, the simulation will stop + real*4, dimension(:,:), allocatable :: zbg + real*4, dimension(:), allocatable :: rtmp + real*4, dimension(:), allocatable :: rtmpz + real*4, dimension(:), allocatable :: rtmpuv + integer*4, dimension(:), allocatable :: uv_index_qt_in_sf + ! + integer :: idummy + integer :: ip + integer :: ibin + integer :: nm + integer :: nmu + integer :: n + integer :: m + integer :: npzq + integer :: npuvq + integer :: npuvs ! ! DEPTHS ! @@ -1296,7 +1414,7 @@ subroutine initialize_domain() ! if (depfile(1:4) /= 'none') then ! - write(*,*)'Reading ',trim(depfile) + write(*,*)'Reading dep file : ',trim(depfile) open(unit = 500, file = trim(depfile), form = 'unformatted', access = 'stream') read(500)zb close(500) @@ -1311,7 +1429,7 @@ subroutine initialize_domain() ! else ! - write(*,*)'Reading ',trim(depfile) + write(*,*)'Reading dep file : ',trim(depfile) ! if (inputtype=='asc') then ! @@ -1329,12 +1447,11 @@ subroutine initialize_domain() if (kcsg(n, m) > 0) then ip = ip + 1 zb(ip) = zbg(n, m) - endif + endif enddo enddo ! deallocate(zbg) - deallocate(kcsg) ! else ! @@ -1374,7 +1491,7 @@ subroutine initialize_domain() ! zbuvmx(ip) = max(zb(nm), zb(nmu)) + huthresh ! - enddo + enddo ! else ! @@ -1789,6 +1906,23 @@ subroutine initialize_domain() ! endif ! + end subroutine + + subroutine initialize_boundaries() + ! + use sfincs_data + ! + implicit none + ! + integer :: idummy + integer :: ip + integer :: ib + integer :: nm + integer :: nmu + integer :: n + integer :: m + integer :: ikcuv2 + ! ! BOUNDARIES ! ! First count number of boundary cells @@ -1951,68 +2085,24 @@ subroutine initialize_domain() endif enddo ! - ! Set flags to turn off viscosity in some cells + end subroutine + + subroutine initialize_roughness() ! - if (nmax/=1) then - do ip = 1, npuv - ! - if (uv_index_u_nmd(ip) == 0) uv_flags_vis(ip) = 0 - if (uv_index_u_nmu(ip) == 0) uv_flags_vis(ip) = 0 - if (uv_index_u_num(ip) == 0) uv_flags_vis(ip) = 0 - if (uv_index_u_ndm(ip) == 0) uv_flags_vis(ip) = 0 - ! - enddo - endif + use sfincs_data ! - ! Turn off advection in some cells (next to inactive points and boundaries) + implicit none ! - if ((advection .or. coriolis) .and. nmax/=1) then - ! - do ip = 1, npuv - ! - ! Turn off advection near inactive points - ! - if (uv_index_u_nmd(ip) == 0) uv_flags_adv(ip) = 0 - if (uv_index_u_nmu(ip) == 0) uv_flags_adv(ip) = 0 - if (uv_index_u_num(ip) == 0) uv_flags_adv(ip) = 0 - if (uv_index_u_ndm(ip) == 0) uv_flags_adv(ip) = 0 - if (uv_index_v_ndm(ip) == 0) uv_flags_adv(ip) = 0 - if (uv_index_v_nm(ip) == 0) uv_flags_adv(ip) = 0 - if (uv_index_v_nmu(ip) == 0) uv_flags_adv(ip) = 0 - if (uv_index_v_ndmu(ip) == 0) uv_flags_adv(ip) = 0 - ! - ! Turn off advection near the boundaries - ! - if (uv_index_u_nmd(ip)>0) then - if (kcuv(uv_index_u_nmd(ip)) == 2) uv_flags_adv(ip) = 0 - endif - if (uv_index_u_nmu(ip)>0) then - if (kcuv(uv_index_u_nmu(ip)) == 2) uv_flags_adv(ip) = 0 - endif - if (uv_index_u_num(ip)>0) then - if (kcuv(uv_index_u_num(ip)) == 2) uv_flags_adv(ip) = 0 - endif - if (uv_index_u_ndm(ip)>0) then - if (kcuv(uv_index_u_ndm(ip)) == 2) uv_flags_adv(ip) = 0 - endif - if (uv_index_v_ndm(ip)>0) then - if (kcuv(uv_index_v_ndm(ip)) == 2) uv_flags_adv(ip) = 0 - endif - if (uv_index_v_nm(ip)>0) then - if (kcuv(uv_index_v_nm(ip)) == 2) uv_flags_adv(ip) = 0 - endif - if (uv_index_v_nmu(ip)>0) then - if (kcuv(uv_index_v_nmu(ip)) == 2) uv_flags_adv(ip) = 0 - endif - if (uv_index_v_ndmu(ip)>0) then - if (kcuv(uv_index_v_ndmu(ip)) == 2) uv_flags_adv(ip) = 0 - endif - ! - enddo - ! - endif + real*4, dimension(:), allocatable :: rghfield ! - ! FRICTION COEFFICIENTS + integer :: idummy + integer :: ip + integer :: nm + integer :: nmu + integer :: n + integer :: m + ! + ! FRICTION COEFFICIENTS (only for regular bathymetry, as for subgrid the Manning's n values are stored in the tables) ! if (.not. subgrid) then ! @@ -2025,7 +2115,7 @@ subroutine initialize_domain() ! Read spatially-varying friction ! allocate(rghfield(np)) - write(*,*)'Reading ',trim(manningfile) + write(*,*)'Reading ',trim(manningfile), ' ...' open(unit = 500, file = trim(manningfile), form = 'unformatted', access = 'stream') read(500)rghfield close(500) @@ -2067,6 +2157,24 @@ subroutine initialize_domain() endif endif ! + end subroutine + + + subroutine initialize_infiltration() + ! + use sfincs_data + ! + implicit none + ! + real*4, dimension(:), allocatable :: rghfield + ! + integer :: idummy + integer :: ip + integer :: nm + integer :: nmu + integer :: n + integer :: m + ! ! INFILTRATION ! ! Infiltration only works when rainfall is activated ! If you want infiltration without rainfall, use a precip file with 0.0s @@ -2373,7 +2481,16 @@ subroutine initialize_domain() store_cumulative_precipitation = .false. ! endif - ! + ! + end subroutine + + + subroutine initialize_storage_volume() + ! + use sfincs_data + ! + implicit none + ! if (use_storage_volume) then ! ! Spatially-varying storage volume for green infra @@ -2405,7 +2522,7 @@ subroutine initialize_hydro() real*4, dimension(:), allocatable :: inizs real*4, dimension(:), allocatable :: iniq ! - integer :: nm, m, n, ivol, num, nmu, ibin, rsttype, ip, ind, iuv + integer :: nm, m, n, ivol, num, nmu, ibin, rsttype, ip, ind, iuv, icuv real*4 :: dzvol real*4 :: facint real*4 :: rdummy @@ -2419,18 +2536,19 @@ subroutine initialize_hydro() logical :: iok ! allocate(zs(np)) - allocate(q(npuv)) - allocate(q0(npuv)) - allocate(uv(npuv)) - allocate(uv0(npuv)) - allocate(kfuv(npuv)) ! Not needed anymore? + ! + ! q, q0, uv, and uv0 also contain the combined values at quadtree transition cells. Plus 1 for dummy values set to 0.0. + ! + allocate(q(npuv + ncuv + 1)) + allocate(q0(npuv + ncuv + 1)) + allocate(uv(npuv + ncuv + 1)) + allocate(uv0(npuv + ncuv + 1)) ! zs = 0.0 q = 0.0 q0 = 0.0 uv = 0.0 uv0 = 0.0 - kfuv = 0 ! if (snapwave) then ! @@ -2482,12 +2600,6 @@ subroutine initialize_hydro() allocate(tsunami_arrival_time(np)) endif ! - ! Make u and v points always active when they are boundary points - ! - do ip = 1, npuv - if (kcuv(ip)==2) kfuv(ip) = 1 ! Get rid of this - enddo - ! ! Initialize water level points ! if (rstfile(1:4) /= 'none') then @@ -2507,16 +2619,17 @@ subroutine initialize_hydro() ! 3: zs - ! 4: zs, q, uvmean and cnb infiltration (writing scs_Se) ! 5: zs, q, uvmean and gai infiltration (writing GA_sigma & GA_F) + ! 6: zs, q, uvmean and hor infiltration (writing rain_T1) ! read(500)rdummy read(500)rsttype read(500)rdummy ! - if (rsttype<1 .or. rsttype >5) then + if (rsttype<1 .or. rsttype >6) then ! ! Give warning, rstfile input rsttype not recognized ! - write(*,*)'WARNING! rstfile not recognized, skipping restartfile input! rsttype should be 1-5, but found rsttype= ', rsttype + write(*,*)'WARNING! rstfile not recognized, skipping restartfile input! rsttype should be 1-6, but found rsttype= ', rsttype ! close(500) ! @@ -2529,7 +2642,7 @@ subroutine initialize_hydro() ! ! Read fluxes q ! - if (rsttype==1 .or. rsttype==2 .or. rsttype==4 .or. rsttype==5) then + if (rsttype==1 .or. rsttype==2 .or. rsttype==4 .or. rsttype==5 .or. rsttype==6) then read(500)rdummy read(500)iniq read(500)rdummy @@ -2539,29 +2652,27 @@ subroutine initialize_hydro() read(500)rdummy endif ! - if (rsttype==4) then ! + if (rsttype==4) then ! Infiltration method cnb + ! + read(500)rdummy + read(500)scs_Se + write(*,*)'Reading scs_Se from rstfile, overwrites input values of: ',trim(sefffile) + ! + elseif (rsttype==5) then ! Infiltration method gai + ! read(500)rdummy - if (inftype == 'cnb') then - ! Infiltration method cnb - read(500)scs_Se - write(*,*)'Reading scs_Se from rstfile, overwrites input values of: ',trim(sefffile) - elseif (inftype == 'hor') then - ! Infiltration method horton - read(500)rain_T1 - write(*,*)'Reading rain_T1 from rstfile, complements input values of: ',trim(fcfile) - endif + read(500)GA_sigma + read(500)rdummy + read(500)GA_F + write(*,*)'Reading GA_sigma from rstfile, overwrites input values of: ',trim(sigmafile) ! - endif - ! - if (rsttype==5) then ! Infiltration method gai - read(500)rdummy - read(500)GA_sigma - read(500)rdummy - read(500)GA_F - ! - write(*,*)'Reading GA_sigma from rstfile, overwrites input values of: ',trim(sigmafile) - ! - endif + elseif (rsttype==6) then ! Infiltration method horton + ! + read(500)rdummy + read(500)rain_T1 + write(*,*)'Reading rain_T1 from rstfile, complements input values of: ',trim(fcfile) + ! + endif ! close(500) ! @@ -2688,6 +2799,17 @@ subroutine initialize_hydro() ! endif ! + ! Loop through combined uv points and determine average uv and q (this also happens at the end of sfincs_momentum.f90). + ! + do icuv = 1, ncuv + ! + ! Average of the two uv points + ! + q(cuv_index_uv(icuv)) = 0.5*(q(cuv_index_uv1(icuv)) + q(cuv_index_uv2(icuv))) + uv(cuv_index_uv(icuv)) = 0.5*(uv(cuv_index_uv1(icuv)) + uv(cuv_index_uv2(icuv))) + ! + enddo + ! if (wavemaker) then ! zsm = zs @@ -2826,5 +2948,56 @@ subroutine initialize_hydro() endif ! end subroutine + + subroutine deallocate_quadtree() + ! + ! De-allocate arrays that are no longer necessary + ! + use sfincs_data + use quadtree + ! + ! In sfincs_data + ! + if (allocated(z_index_uv_md1)) deallocate(z_index_uv_md1) + if (allocated(z_index_uv_md2)) deallocate(z_index_uv_md2) + if (allocated(z_index_uv_mu1)) deallocate(z_index_uv_mu1) + if (allocated(z_index_uv_mu2)) deallocate(z_index_uv_mu2) + if (allocated(z_index_uv_nd1)) deallocate(z_index_uv_nd1) + if (allocated(z_index_uv_nd2)) deallocate(z_index_uv_nd2) + if (allocated(z_index_uv_nu1)) deallocate(z_index_uv_nu1) + if (allocated(z_index_uv_nu2)) deallocate(z_index_uv_nu2) + ! + ! In quadtree + ! + if (allocated(quadtree_level)) deallocate(quadtree_level) + if (allocated(quadtree_md)) deallocate(quadtree_md) + if (allocated(quadtree_md1)) deallocate(quadtree_md1) + if (allocated(quadtree_md2)) deallocate(quadtree_md2) + if (allocated(quadtree_mu)) deallocate(quadtree_mu) + if (allocated(quadtree_mu1)) deallocate(quadtree_mu1) + if (allocated(quadtree_mu2)) deallocate(quadtree_mu2) + if (allocated(quadtree_nd)) deallocate(quadtree_nd) + if (allocated(quadtree_nd1)) deallocate(quadtree_nd1) + if (allocated(quadtree_nd2)) deallocate(quadtree_nd2) + if (allocated(quadtree_nu)) deallocate(quadtree_nu) + if (allocated(quadtree_nu1)) deallocate(quadtree_nu1) + if (allocated(quadtree_nu2)) deallocate(quadtree_nu2) + if (allocated(quadtree_n)) deallocate(quadtree_n) + if (allocated(quadtree_m)) deallocate(quadtree_m) + if (allocated(quadtree_n_oddeven)) deallocate(quadtree_n_oddeven) + if (allocated(quadtree_m_oddeven)) deallocate(quadtree_m_oddeven) + if (allocated(quadtree_xz)) deallocate(quadtree_xz) + if (allocated(quadtree_yz)) deallocate(quadtree_yz) + if (allocated(quadtree_zz)) deallocate(quadtree_zz) + if (allocated(quadtree_nm_indices)) deallocate(quadtree_nm_indices) + if (allocated(quadtree_first_point_per_level)) deallocate(quadtree_first_point_per_level) + if (allocated(quadtree_last_point_per_level)) deallocate(quadtree_last_point_per_level) + if (allocated(quadtree_dxr)) deallocate(quadtree_dxr) + if (allocated(quadtree_dyr)) deallocate(quadtree_dyr) + if (allocated(quadtree_mask)) deallocate(quadtree_mask) + if (allocated(quadtree_snapwave_mask)) deallocate(quadtree_snapwave_mask) ! + end subroutine + + end module diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index f3374ab0..4e0999f0 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -13,6 +13,28 @@ subroutine read_sfincs_input() ! integer*8 dtsec ! + ! Temporary variables + ! + integer iradstr + integer igeo + integer icoriolis + integer iamprblock + integer iglobal + integer itsunamitime + integer ispinupmeteo + integer isnapwave + integer iwindmax + integer ioutfixed + integer iadvection + integer istorefw + integer istorewavdir + integer imanning2d + integer iviscosity + integer isubgrid + integer iwavemaker + integer iwavemaker_spectrum + integer ispwprecip + ! character*256 wmsigstr ! write(*,*)'Reading input file ...' @@ -57,6 +79,7 @@ subroutine read_sfincs_input() call read_char_input(500,'outputformat',outputtype,'asc') call read_char_input(500,'outputtype_map',outputtype_map,'nil') call read_char_input(500,'outputtype_his',outputtype_his,'nil') + call read_int_input(500,'nc_deflate_level',nc_deflate_level,2) call read_int_input(500,'bndtype',bndtype,1) call read_int_input(500,'advection',iadvection,0) call read_int_input(500,'nfreqsig',nfreqsig,100) @@ -153,6 +176,7 @@ 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') ! ! Output call read_char_input(500,'obsfile',obsfile,'none') @@ -169,7 +193,7 @@ subroutine read_sfincs_input() call read_int_input(500,'storeqdrain',storeqdrain,1) call read_int_input(500,'storezvolume',storezvolume,0) call read_int_input(500,'writeruntime',wrttimeoutput,0) - call read_int_input(500,'debug',idebug,0) + call read_logical_input(500,'debug',debug,.false.) call read_int_input(500,'storemeteo',storemeteo,0) call read_int_input(500,'storemaxwind',iwindmax,0) call read_int_input(500,'storefw', istorefw, 0) @@ -337,11 +361,6 @@ subroutine read_sfincs_input() write_time_output = .true. endif ! - debug = .false. - if (idebug==1) then - debug = .true. - endif - ! radstr = .false. if (iradstr==1) then radstr = .true. @@ -427,6 +446,11 @@ subroutine read_sfincs_input() spinup_meteo = .false. endif ! + use_spw_precip = .true. + if (ispwprecip==0) then + use_spw_precip = .false. + endif + ! snapwave = .false. if (isnapwave==1) then snapwave = .true. @@ -442,6 +466,11 @@ subroutine read_sfincs_input() advection = .true. endif ! + thetasmoothing = .false. + if (theta<0.9999) then ! Note, for reliability in terms of precision, is written as 0.9999 + thetasmoothing = .true. + endif + ! store_wave_forces = .false. if (istorefw==1) then store_wave_forces = .true. @@ -593,7 +622,40 @@ subroutine read_char_input(fileid,keyword,value,default) ! end subroutine - + + subroutine read_logical_input(fileid,keyword,value,default) + ! + character(*), intent(in) :: keyword + character(len=256) :: keystr0 + character(len=256) :: keystr + character(len=256) :: valstr + character(len=256) :: line + integer, intent(in) :: fileid + logical, intent(in) :: default + logical, intent(out) :: value + integer j,stat,ilen + ! + value = default + rewind(fileid) + do while(.true.) + read(fileid,'(a)',iostat = stat)line + if (stat==-1) exit + j=index(line,'=') + keystr0 = trim(line(1:j-1)) + call notabs(keystr0,keystr,ilen) + if (trim(keystr)==trim(keyword)) then + valstr = adjustl(trim(line(j+1:256))) + if (valstr(1:1) == '1' .or. valstr(1:1) == 'y' .or. valstr(1:1) == 'Y' .or. valstr(1:1) == 't' .or. valstr(1:1) == 'T') then + value = .true. + else + value = .false. + endif + exit + endif + enddo + ! + end subroutine + subroutine notabs(INSTR,OUTSTR,ILEN) ! @(#) convert tabs in input to spaces in output while maintaining columns, assuming a tab is set every 8 characters ! diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index d108ae59..2a9fa1b7 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -84,8 +84,8 @@ function sfincs_initialize(config_file) result(ierr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - build_revision = '$Rev: v2.0.3-alpha' - build_date = '$Date: 2023-10-31' + build_revision = '$Rev: v2.0.3-Cauberg' + build_date = '$Date: 2023-11-15' ! write(*,'(a)')'' write(*,*)'----------- Welcome to SFINCS -----------' @@ -130,9 +130,7 @@ function sfincs_initialize(config_file) result(ierr) ! call read_meteo_data() ! Reads meteo data (amu, amv, spw file etc.) ! - call initialize_domain() ! Reads dep, msk, index files. Creates index, flag and depth arrays - ! - call initialize_hydro() ! Initializes water levels, fluxes, flags + call initialize_domain() ! Reads dep, msk, index files. Creates index, flag and depth arrays. Initializes water levels, fluxes, flags ! call read_structures() ! Reads thd files. Sets kcuv to zero where necessary ! @@ -204,6 +202,10 @@ function sfincs_initialize(config_file) result(ierr) ! call initialize_output(tmapout, tmaxout, thisout, trstout) ! + ! Quadtree no longer needed, so deallocate (this is done in sfincs_domain.f90) + ! + call deallocate_quadtree() + ! !call acc_init( acc_device_nvidia ) ! end function sfincs_initialize @@ -222,11 +224,11 @@ function sfincs_update(dtrange) result(ierr) ierr = -1 ! !$acc data, copyin( kcs, kcuv, zs, q, q0, uv, uv0, zb, zbuv, zbuvmx, zsmax, qmax, vmax, twet, zsm, z_volume, & - !$acc z_flags_iref, z_flags_type, uv_flags_iref, uv_flags_type, uv_flags_vis, uv_flags_adv, uv_flags_dir, & + !$acc z_flags_iref, uv_flags_iref, uv_flags_type, uv_flags_dir, & !$acc index_kcuv2, nmikcuv2, nmbkcuv2, ibkcuv2, zsb, zsb0, ibuvdir, uvmean, & !$acc subgrid_uv_zmin, subgrid_uv_zmax, subgrid_uv_hrep, subgrid_uv_navg, subgrid_uv_hrep_zmax, subgrid_uv_navg_zmax, & !$acc subgrid_z_zmin, subgrid_z_zmax, subgrid_z_dep, subgrid_z_volmax, & - !$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 z_index_uv_md, z_index_uv_nd, z_index_uv_mu, z_index_uv_nu, & !$acc uv_index_z_nm, uv_index_z_nmu, uv_index_u_nmd, uv_index_u_nmu, uv_index_u_ndm, uv_index_u_num, & !$acc uv_index_v_ndm, uv_index_v_ndmu, uv_index_v_nm, uv_index_v_nmu, & !$acc nmindsrc, qtsrc, drainage_type, drainage_params, & diff --git a/source/src/sfincs_meteo.f90 b/source/src/sfincs_meteo.f90 index c15991b4..764f973e 100644 --- a/source/src/sfincs_meteo.f90 +++ b/source/src/sfincs_meteo.f90 @@ -59,18 +59,34 @@ subroutine read_meteo_data() enddo ! if (spw_nquant==4) then - if (ispwprecip==0) then - spw_precip = .false. - write(*,*)'Info : Overruled to not use precipitation from spiderweb input ...' - else - spw_precip = .true. - endif + if (use_spw_precip) then + spw_precip = .true. + else + write(*,*)'Info : Overruled to not use precipitation from spiderweb input ...' + spw_precip = .false. + endif endif ! + endif + ! + if (netspwfile(1:4) /= 'none') then + ! + ! Write statement to log + write(*,*)'Reading netcdf spiderweb file ...' + ! + ! Read information + call read_netcdf_spw_data() + ! + endif + ! + if (spwfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then + ! if (utmzone/='nil') then ! ! Convert spiderweb coordinates to utm zone ! + write(*,*)'Converting spiderweb coordinates to utm zone ...' + ! do it = 1, spw_nt call deg2utm(spw_ye(it),spw_xe(it),xx,yy,utmzone) spw_xe(it) = xx @@ -79,8 +95,8 @@ subroutine read_meteo_data() ! endif ! - endif - ! + endif + ! if (amufile(1:4) /= 'none') then ! write(*,*)'Reading amu and amv file ...' @@ -1303,7 +1319,7 @@ subroutine update_meteo_fields(t, tloop) ! endif ! - if (spwfile(1:4) /= 'none') then + if (spwfile(1:4) /= 'none' .or. netspwfile(1:4) /= 'none') then ! call update_spiderweb_data() ! @@ -1317,18 +1333,6 @@ subroutine update_meteo_fields(t, tloop) ! endif ! -! if (precip .and. store_cumulative_precipitation) then -! ! -! call update_cumprcp_map() -! ! -! endif -! ! -! if (infiltration2d) then -! ! -! call update_infiltration_map() -! ! -! endif - ! call system_clock(count1, count_rate, count_max) tloop = tloop + 1.0*(count1 - count0)/count_rate ! diff --git a/source/src/sfincs_momentum.f90 b/source/src/sfincs_momentum.f90 index e5fb5b3c..1c965c4c 100644 --- a/source/src/sfincs_momentum.f90 +++ b/source/src/sfincs_momentum.f90 @@ -21,12 +21,17 @@ subroutine compute_fluxes(dt, min_dt, tloop) integer :: ip integer :: nm integer :: nmu + integer :: n + integer :: m integer :: idir integer :: iref integer :: itype integer :: iuv integer :: ind + integer :: icuv + integer :: iuv1 + integer :: iuv2 ! real*4 :: hu real*4 :: dxuvinv @@ -65,9 +70,9 @@ subroutine compute_fluxes(dt, min_dt, tloop) real*4 :: one_minus_facint ! real*4, parameter :: expo = 1.0/3.0 -! integer, parameter :: expo = 1 + ! integer, parameter :: expo = 1 ! - logical :: iadv, ivis, icorio, iok + logical :: iok ! call system_clock(count0, count_rate, count_max) ! @@ -80,7 +85,7 @@ subroutine compute_fluxes(dt, min_dt, tloop) !$omp do !$acc kernels, present(q, q0, uv, uv0), async(1) !$acc loop independent, private(nm) - do ip = 1, npuv + do ip = 1, npuv + ncuv ! q0(ip) = q(ip) uv0(ip) = uv(ip) @@ -91,13 +96,13 @@ subroutine compute_fluxes(dt, min_dt, tloop) !$omp end parallel ! !$omp parallel & - !$omp private ( ip,hu,qsm,qx_nm,nm,nmu,frc,adv,idir,itype,iadv,ivis,icorio,iref,dxuvinv,dxuv2inv,dyuvinv,dyuv2inv, & + !$omp private ( ip,hu,qsm,qx_nm,nm,nmu,frc,adv,idir,itype,iref,dxuvinv,dxuv2inv,dyuvinv,dyuv2inv, & !$omp qx_nmd,qx_nmu,uu_nm,uu_nmd,uu_nmu,uu_num,uu_ndm,vu, & !$omp fcoriouv,gnavg2,iok,zsu,dzuv,iuv,facint,fwmax,zmax,zmin,one_minus_facint,qvt ) & !$omp reduction ( min : min_dt ) !$omp do schedule ( dynamic, 256 ) !$acc kernels, present( kcuv, zs, q, q0, uv, uv0, min_dt, & - !$acc uv_flags_iref, uv_flags_type, uv_flags_vis, uv_flags_adv, uv_flags_dir, & + !$acc uv_flags_iref, uv_flags_type, uv_flags_dir, & !$acc subgrid_uv_zmin, subgrid_uv_zmax, subgrid_uv_hrep, subgrid_uv_navg, subgrid_uv_hrep_zmax, subgrid_uv_navg_zmax, & !$acc uv_index_z_nm, uv_index_z_nmu, uv_index_u_nmd, uv_index_u_nmu, uv_index_u_ndm, uv_index_u_num, & !$acc uv_index_v_ndm, uv_index_v_ndmu, uv_index_v_nm, uv_index_v_nmu, & @@ -145,46 +150,6 @@ subroutine compute_fluxes(dt, min_dt, tloop) ! idir = uv_flags_dir(ip) ! 0 is u, 1 is v ! - ! Check if viscosity term is computed for this point - ! - if (viscosity) then - if (uv_flags_vis(ip)==1) then - ivis = .true. - else - ivis = .false. - endif - else - ivis = .false. - endif - ! - ! Check if advection term is computed for this point - ! - if (advection) then - if (uv_flags_adv(ip)==1) then - iadv = .true. - else - iadv = .false. - endif - else - iadv = .false. - endif - ! - ! Check if Coriolis term is computed for this point - ! - if (coriolis) then - ! - if (uv_flags_adv(ip)==1) then - icorio = .true. - else - icorio = .false. - endif - ! - else - ! - icorio = .false. - ! - endif - ! ! Determine grid spacing (and coriolis factor fcoriouv) ! if (crsgeo) then @@ -278,7 +243,7 @@ subroutine compute_fluxes(dt, min_dt, tloop) ! qx_nm = q0(ip) ! - if (iadv) then + if (advection .or. coriolis .or. viscosity) then ! ! First term ! @@ -292,20 +257,7 @@ subroutine compute_fluxes(dt, min_dt, tloop) ! endif ! - if (ivis) then - ! - if (.not. iadv) then - uu_nm = uv0(ip) - uu_nmd = uv0(uv_index_u_nmd(ip)) - uu_nmu = uv0(uv_index_u_nmu(ip)) - endif - ! - uu_ndm = uv0(uv_index_u_ndm(ip)) - uu_num = uv0(uv_index_u_num(ip)) - ! - endif - ! - if (theta<0.9999 .and. .not.iadv) then ! for backward compatibility + if (thetasmoothing .and. .not.advection) then ! for backward compatibility ! Note, for reliability in terms of precision, is written as 0.9999 qx_nmd = q0(uv_index_u_nmd(ip)) qx_nmu = q0(uv_index_u_nmu(ip)) @@ -355,46 +307,53 @@ subroutine compute_fluxes(dt, min_dt, tloop) ! ! Advection term ! - if (iadv) then + if (advection) then ! ! 1D upwind advection slightly more robust than central scheme ! - if (qx_nm>1.0e-6) then - ! - adv = - ( qx_nm*uu_nm - qx_nmd*uu_nmd ) * dxuvinv - ! - elseif (qx_nm<-1.0e-6) then + if ( kcuv(uv_index_u_nmd(ip))==1 .and. kcuv(uv_index_u_nmu(ip))==1 ) then ! - adv = - ( qx_nmu*uu_nmu - qx_nm*uu_nm ) * dxuvinv - ! - else - ! - adv = 0.0 + ! But only at regular points + ! + if (qx_nm>1.0e-6) then + ! + adv = - ( qx_nm*uu_nm - qx_nmd*uu_nmd ) * dxuvinv + ! + elseif (qx_nm<-1.0e-6) then + ! + adv = - ( qx_nmu*uu_nmu - qx_nm*uu_nm ) * dxuvinv + ! + else + ! + adv = 0.0 + ! + endif ! - endif - ! - qvt = q0(uv_index_v_ndm(ip)) + q0(uv_index_v_ndmu(ip)) + q0(uv_index_v_nm(ip)) + q0(uv_index_v_nmu(ip)) - ! - if ( qvt > 1.0e-6 ) then + qvt = q0(uv_index_v_ndm(ip)) + q0(uv_index_v_ndmu(ip)) + q0(uv_index_v_nm(ip)) + q0(uv_index_v_nmu(ip)) ! - adv = adv - (q0(uv_index_v_ndm(ip)) + q0(uv_index_v_ndmu(ip))) * (uu_nm - uu_ndm) * dyuvinv / 2 + if ( qvt > 1.0e-6 ) then + ! + adv = adv - (q0(uv_index_v_ndm(ip)) + q0(uv_index_v_ndmu(ip))) * (uu_nm - uu_ndm) * dyuvinv / 2 + ! + elseif ( qvt < -1.0e-6 ) then + ! + adv = adv - (q0(uv_index_v_nm(ip)) + q0(uv_index_v_nmu(ip))) * (uu_num - uu_nm) * dyuvinv / 2 + ! + endif ! - elseif ( qvt < -1.0e-6 ) then + ! Let's try without the advection limiter ! - adv = adv - (q0(uv_index_v_nm(ip)) + q0(uv_index_v_nmu(ip))) * (uu_num - uu_nm) * dyuvinv / 2 + !adv = min(max(adv, -advlim), advlim) ! + frc = frc + adv + ! endif ! - ! Let's try without the advection limiter -! adv = min(max(adv, -advlim), advlim) - ! - frc = frc + adv - ! endif ! ! Viscosity term ! - if (ivis) then + if (viscosity) then ! frc = frc + nuvisc * hu * ( (uu_nmu - 2*uu_nm + uu_nmd )*dxuv2inv + (uu_num - 2*uu_nm + uu_ndm )*dyuv2inv ) ! @@ -402,7 +361,7 @@ subroutine compute_fluxes(dt, min_dt, tloop) ! ! Coriolis term ! - if (icorio) then + if (coriolis) then ! vu = (uv(uv_index_v_ndm(ip)) + uv(uv_index_v_ndmu(ip)) + uv(uv_index_v_nm(ip)) + uv(uv_index_v_nmu(ip))) / 4 ! @@ -476,26 +435,25 @@ subroutine compute_fluxes(dt, min_dt, tloop) endif ! ! 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 + if (thetasmoothing) then ! ! 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 + if ( kcuv(uv_index_u_nmd(ip))==1 .and. kcuv(uv_index_u_nmu(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) @@ -519,11 +477,24 @@ subroutine compute_fluxes(dt, min_dt, tloop) !$omp end parallel !$acc end kernels ! + ! Loop through combined uv points and determine average uv and q + ! + !$omp parallel & + !$omp private ( icuv ) + !$omp do + do icuv = 1, ncuv + ! Average of the two uv points + q(cuv_index_uv(icuv)) = 0.5*(q(cuv_index_uv1(icuv)) + q(cuv_index_uv2(icuv))) + uv(cuv_index_uv(icuv)) = 0.5*(uv(cuv_index_uv1(icuv)) + uv(cuv_index_uv2(icuv))) + enddo + !$omp end do + !$omp end parallel + ! !$acc update host(min_dt), async(1) ! call system_clock(count1, count_rate, count_max) tloop = tloop + 1.0*(count1 - count0)/count_rate ! end subroutine - + ! end module \ No newline at end of file diff --git a/source/src/sfincs_ncinput.F90 b/source/src/sfincs_ncinput.F90 index 53e27c84..1ba0d886 100644 --- a/source/src/sfincs_ncinput.F90 +++ b/source/src/sfincs_ncinput.F90 @@ -42,6 +42,15 @@ module sfincs_ncinput integer :: px_varid, py_varid integer :: time_varid integer :: prcp_varid + end type + type net_type_spw + integer :: ncid + integer :: ncols_dimid, nrows_dimid + integer :: time_dimid + integer :: time_varid + integer :: range_varid,azimuth_varid + integer :: xeye_varid, yeye_varid, peye_varid + integer :: wind_x_varid, wind_y_varid, pressure_varid, precip_varid end type ! type(net_type_bndbzsbzi) :: net_file_bndbzsbzi @@ -49,6 +58,7 @@ module sfincs_ncinput type(net_type_amuv) :: net_file_amuv type(net_type_amp) :: net_file_amp type(net_type_ampr) :: net_file_ampr + type(net_type_spw) :: net_file_spw contains @@ -551,6 +561,128 @@ subroutine read_netcdf_ampr_data() ! NF90(nf90_close(net_file_ampr%ncid)) ! + end subroutine + ! + ! + subroutine read_netcdf_spw_data() + ! + ! Read the orginal (Bert) implementation + ! + use sfincs_date + use netcdf + use sfincs_data + use sfincs_spiderweb + + implicit none + ! + integer nt, status + real*4, dimension(:,:,:), allocatable :: prtmp + real*4, dimension(:,:,:), allocatable :: ampr_prtmp + + write(*,*)'subroutine NetCDF spiderweb ...' + ! + ! Actual reading of data + NF90(nf90_open(trim(netspwfile), NF90_CLOBBER, net_file_spw%ncid)) + ! + ! Get dimensions id's: time, range and azimuth + NF90(nf90_inq_dimid(net_file_spw%ncid, "time", net_file_spw%time_dimid)) + NF90(nf90_inq_dimid(net_file_spw%ncid, "range", net_file_spw%nrows_dimid)) + NF90(nf90_inq_dimid(net_file_spw%ncid, "azimuth", net_file_spw%ncols_dimid)) + ! + ! Get dimensions sizes: time, cols, rows + NF90(nf90_inquire_dimension(net_file_spw%ncid, net_file_spw%time_dimid, len = spw_nt)) !nr of timesteps in file + NF90(nf90_inquire_dimension(net_file_spw%ncid, net_file_spw%nrows_dimid, len = spw_nrows)) !nr of rows + NF90(nf90_inquire_dimension(net_file_spw%ncid, net_file_spw%ncols_dimid, len = spw_ncols)) !nr of cols + + ! Get variable id's + NF90(nf90_inq_varid(net_file_spw%ncid, "time", net_file_spw%time_varid) ) + NF90(nf90_inq_varid(net_file_spw%ncid, "range", net_file_spw%range_varid) ) + NF90(nf90_inq_varid(net_file_spw%ncid, "azimuth", net_file_spw%azimuth_varid) ) + NF90(nf90_inq_varid(net_file_spw%ncid, "longitude_eye", net_file_spw%xeye_varid) ) + NF90(nf90_inq_varid(net_file_spw%ncid, "latitude_eye", net_file_spw%yeye_varid) ) + NF90(nf90_inq_varid(net_file_spw%ncid, "wind_x", net_file_spw%wind_x_varid) ) + NF90(nf90_inq_varid(net_file_spw%ncid, "wind_y", net_file_spw%wind_y_varid) ) + NF90(nf90_inq_varid(net_file_spw%ncid, "pressure", net_file_spw%pressure_varid) ) + ! + ! Attempt to get the variable ID for "precipitation" + status = NF90_INQ_VARID(net_file_spw%ncid, "precipitation", net_file_spw%precip_varid) + ! + ! Check the status + if (status /= NF90_NOERR) then + ! Handle error: variable does not exist + write(*,*) "Error: Variable 'precipitation' does not exist in the NetCDF file." + spw_precip = .false. + precip = .false. + else + NF90(nf90_inq_varid(net_file_spw%ncid, "precipitation", net_file_spw%precip_varid) ) + spw_precip = .true. + precip = .true. + write(*,*)'Turning on process: Precipitation from spwfile' + endif + ! + ! Allocate + allocate(spw_times(spw_nt)) + allocate(spw_xe(spw_nt)) + allocate(spw_ye(spw_nt)) + allocate(spw_pressure_eye(spw_nt)) + allocate(spw_radia(spw_nrows)) + allocate(spw_wu(spw_nt, spw_nrows, spw_ncols)) + allocate(spw_wv(spw_nt, spw_nrows, spw_ncols)) + allocate(spw_pdrp(spw_nt, spw_nrows, spw_ncols)) + allocate(spw_prcp(spw_nt, spw_nrows, spw_ncols)) + allocate(spw_pdrp01(spw_nrows, spw_ncols)) + allocate(spw_prcp01(spw_nrows, spw_ncols)) + allocate(spw_wu01(spw_nrows, spw_ncols)) + allocate(spw_wv01(spw_nrows, spw_ncols)) + + ! Support variables + allocate(prtmp(spw_nrows, spw_ncols,1)) + allocate(ampr_prtmp(1, spw_nrows, spw_ncols)) + ! + ! Read values from netcdf + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%time_varid, spw_times(:))) + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%xeye_varid, spw_xe(:))) + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%yeye_varid, spw_ye(:))) + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%range_varid, spw_radia(:))) + ! + ! We only need to know the maxima + spw_radius = spw_radia(spw_nrows) + dradspw = spw_radius/spw_nrows + dphispw = 2*pi/spw_ncols + ! + ! Read matrix values + do nt = 1, spw_nt + ! + ! Read wind_x + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%wind_x_varid, prtmp, start = (/ 1, 1, nt /), count = (/ spw_ncols, spw_nrows, 1 /))) ! be aware of start indices + ampr_prtmp = reshape( prtmp, (/ 1, spw_nrows, spw_ncols /), ORDER = (/ 3, 2, 1 /)) + spw_wu(nt,:,:) = ampr_prtmp(1,:,:) + ! + ! Read wind_y + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%wind_y_varid, prtmp, start = (/ 1, 1, nt /), count = (/ spw_ncols, spw_nrows, 1 /))) ! be aware of start indices + ampr_prtmp = reshape( prtmp, (/ 1, spw_nrows, spw_ncols /), ORDER = (/ 3, 2, 1 /)) + spw_wv(nt,:,:) = ampr_prtmp(1,:,:) + ! + ! Read pressure + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%pressure_varid, prtmp, start = (/ 1, 1, nt /), count = (/ spw_ncols, spw_nrows, 1 /))) ! be aware of start indices + ampr_prtmp = reshape( prtmp, (/ 1, spw_nrows, spw_ncols /), ORDER = (/ 3, 2, 1 /)) + spw_pdrp(nt,:,:) = ampr_prtmp(1,:,:) + ! + ! Read rainfall + if (spw_precip) then + NF90(nf90_get_var(net_file_spw%ncid, net_file_spw%precip_varid, prtmp, start = (/ 1, 1, nt /), count = (/ spw_ncols, spw_nrows, 1 /))) ! be aware of start indices + ampr_prtmp = reshape( prtmp, (/ 1, spw_nrows, spw_ncols /), ORDER = (/ 3, 2, 1 /)) + spw_prcp(nt,:,:) = ampr_prtmp(1,:,:) + endif + enddo + ! + ! Read time attibute and convert time + NF90(nf90_get_att(net_file_spw%ncid, net_file_spw%time_varid,'units', treftimefews)) + spw_times = convert_spw_nc_date(spw_times, spw_nt, treftimefews, trefstr) + ! + ! Close netcdf + NF90(nf90_close(net_file_spw%ncid)) + end subroutine ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90 index 9816dce7..03304f7f 100644 --- a/source/src/sfincs_ncoutput.F90 +++ b/source/src/sfincs_ncoutput.F90 @@ -80,7 +80,7 @@ subroutine ncoutput_regular_map_init() ! !write(*,*) trim(nf90_inq_libvers()) ! - NF90(nf90_create('sfincs_map.nc', NF90_CLOBBER, map_file%ncid)) + NF90(nf90_create('sfincs_map.nc', ior(NF90_CLOBBER,NF90_NETCDF4), map_file%ncid)) ! ! Create dimensions ! grid, time, points @@ -110,7 +110,8 @@ subroutine ncoutput_regular_map_init() ! ! Domain ! - NF90(nf90_def_var(map_file%ncid, 'x', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%face_x_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var(map_file%ncid, 'x', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%face_x_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%face_x_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%face_x_varid, '_FillValue', FILL_VALUE)) if (crsgeo) then NF90(nf90_put_att(map_file%ncid, map_file%face_x_varid, 'units', 'degrees')) @@ -123,7 +124,8 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%face_x_varid, 'grid_mapping', 'crs')) NF90(nf90_put_att(map_file%ncid, map_file%face_x_varid, 'grid', 'sfincsgrid')) ! - NF90(nf90_def_var(map_file%ncid, 'y', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%face_y_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var(map_file%ncid, 'y', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%face_y_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%face_y_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%face_y_varid, '_FillValue', FILL_VALUE)) if (crsgeo) then NF90(nf90_put_att(map_file%ncid, map_file%face_y_varid, 'units', 'degrees')) @@ -137,6 +139,7 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%face_y_varid, 'grid', 'sfincsgrid')) ! NF90(nf90_def_var(map_file%ncid, 'corner_x', NF90_FLOAT, (/map_file%corner_m_dimid, map_file%corner_n_dimid/), map_file%corner_x_varid)) ! location of u points in cell corner + NF90(nf90_def_var_deflate(map_file%ncid, map_file%corner_x_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%corner_x_varid, '_FillValue', FILL_VALUE)) if (crsgeo) then NF90(nf90_put_att(map_file%ncid, map_file%corner_x_varid, 'units', 'degrees')) @@ -150,6 +153,7 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%corner_x_varid, 'grid', 'sfincsgrid')) ! NF90(nf90_def_var(map_file%ncid, 'corner_y', NF90_FLOAT, (/map_file%corner_m_dimid, map_file%corner_n_dimid/), map_file%corner_y_varid)) ! location of v points in cell corner + NF90(nf90_def_var_deflate(map_file%ncid, map_file%corner_y_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%corner_y_varid, '_FillValue', FILL_VALUE)) if (crsgeo) then NF90(nf90_put_att(map_file%ncid, map_file%corner_y_varid, 'units', 'degrees')) @@ -162,11 +166,11 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%corner_y_varid, 'grid_mapping', 'crs')) NF90(nf90_put_att(map_file%ncid, map_file%corner_y_varid, 'grid', 'sfincsgrid')) ! - NF90(nf90_def_var(map_file%ncid, 'crs', NF90_INT, map_file%crs_varid)) ! For EPSG code + NF90(nf90_def_var(map_file%ncid, 'crs', NF90_INT, map_file%crs_varid)) ! For EPSG code NF90(nf90_put_att(map_file%ncid, map_file%crs_varid, 'EPSG', '-')) NF90(nf90_put_att(map_file%ncid, map_file%crs_varid, 'epsg_code', 'EPSG:' // trim(epsg_code) )) !--> add epsg_code like FEWS wants ! - NF90(nf90_def_var(map_file%ncid, 'sfincsgrid', NF90_INT, map_file%grid_varid)) ! For neat grid clarification + NF90(nf90_def_var(map_file%ncid, 'sfincsgrid', NF90_INT, map_file%grid_varid)) ! For neat grid clarification NF90(nf90_put_att(map_file%ncid, map_file%grid_varid, 'cf_role', 'grid_topology')) NF90(nf90_put_att(map_file%ncid, map_file%grid_varid, 'topology_dimension', 2)) NF90(nf90_put_att(map_file%ncid, map_file%grid_varid, 'node_dimensions', 'n m')) !or n m? @@ -175,18 +179,20 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%grid_varid, 'face_coordinates', 'x y')) NF90(nf90_put_att(map_file%ncid, map_file%grid_varid, 'corner_coordinates', 'corner_x corner_y')) ! - NF90(nf90_def_var(map_file%ncid, 'msk', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%msk_varid)) ! input msk value in cell centre + NF90(nf90_def_var(map_file%ncid, 'msk', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%msk_varid)) ! input msk value in cell centre NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'units', '-')) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'standard_name', 'land_binary_mask')) ! land_binary_mask but with added boundary=2 NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'long_name', 'msk_active_cells')) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'description', 'inactive=0, active=1, normal_boundary=2, outflow_boundary=3')) - NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'coordinates', 'x y')) + NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'coordinates', 'x y')) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%msk_varid, 1, 1, nc_deflate_level)) ! deflate ! ! Infiltration map ! if (infiltration) then - NF90(nf90_def_var(map_file%ncid, 'qinf', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%qinf_varid)) + NF90(nf90_def_var(map_file%ncid, 'qinf', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%qinf_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%qinf_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'coordinates', 'x y')) if (inftype == 'cna') then @@ -213,6 +219,7 @@ subroutine ncoutput_regular_map_init() endif ! NF90(nf90_def_var(map_file%ncid, 'zb', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%zb_varid)) ! bed level in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zb_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, 'standard_name', 'altitude')) @@ -230,6 +237,7 @@ subroutine ncoutput_regular_map_init() ! Time varying map output ! NF90(nf90_def_var(map_file%ncid, 'zs', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%zs_varid)) ! time-varying water level map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zs_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%zs_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zs_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%zs_varid, 'standard_name', 'sea_surface_height_above_reference_level')) @@ -238,6 +246,7 @@ subroutine ncoutput_regular_map_init() ! if (subgrid .eqv. .false. .or. store_hsubgrid .eqv. .true.) then NF90(nf90_def_var(map_file%ncid, 'h', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%h_varid)) ! time-varying water depth map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%h_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%h_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%h_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%h_varid, 'standard_name', 'depth')) @@ -249,7 +258,8 @@ subroutine ncoutput_regular_map_init() ! if (store_velocity) then ! - NF90(nf90_def_var(map_file%ncid, 'u', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%u_varid)) ! time-varying u map + NF90(nf90_def_var(map_file%ncid, 'u', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%u_varid)) ! time-varying u map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%u_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%u_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%u_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%u_varid, 'standard_name', 'eastward_sea_water_velocity')) ! not truly eastward when rotated, eastward_sea_water_velocity @@ -257,6 +267,7 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%u_varid, 'coordinates', 'x y')) ! NF90(nf90_def_var(map_file%ncid, 'v', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%v_varid)) ! time-varying u map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%v_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%v_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%v_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%v_varid, 'standard_name', 'northward_sea_water_velocity')) ! not truly eastward when rotated, eastward_sea_water_velocity @@ -265,8 +276,10 @@ subroutine ncoutput_regular_map_init() endif ! ! Store S_effective (only for CN method with recovery) + ! if (inftype == 'cnb') then NF90(nf90_def_var(map_file%ncid, 'Seff', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%Seff_varid)) ! time-varying S + NF90(nf90_def_var_deflate(map_file%ncid, map_file%Seff_varid, 1, 1, nc_deflate_level)) ! deflate 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', 'Se')) @@ -305,6 +318,7 @@ subroutine ncoutput_regular_map_init() ! if (store_maximum_waterlevel) then NF90(nf90_def_var(map_file%ncid, 'zsmax', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%timemax_dimid/), map_file%zsmax_varid)) ! time-varying maximum water level map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zsmax_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%zsmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zsmax_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%zsmax_varid, 'standard_name', 'maximum of sea_surface_height_above_reference_level')) @@ -314,6 +328,7 @@ subroutine ncoutput_regular_map_init() ! if (store_cumulative_precipitation) then NF90(nf90_def_var(map_file%ncid, 'cumprcp', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%timemax_dimid/), map_file%cumprcp_varid)) ! time-varying maximum water level map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%cumprcp_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%cumprcp_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%cumprcp_varid, 'units', 'mm')) NF90(nf90_put_att(map_file%ncid, map_file%cumprcp_varid, 'long_name', 'cumulative_precipitation_depth')) @@ -324,6 +339,7 @@ subroutine ncoutput_regular_map_init() ! if (store_twet) then NF90(nf90_def_var(map_file%ncid, 'tmax', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%timemax_dimid/), map_file%tmax_varid)) ! time-varying duration wet cell + NF90(nf90_def_var_deflate(map_file%ncid, map_file%tmax_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%tmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%tmax_varid, 'units', 'seconds')) NF90(nf90_put_att(map_file%ncid, map_file%tmax_varid, 'standard_name', 'duration cell is considered wet')) @@ -335,6 +351,7 @@ subroutine ncoutput_regular_map_init() if (store_maximum_waterlevel) then if (subgrid .eqv. .false. .or. store_hsubgrid .eqv. .true.) then NF90(nf90_def_var(map_file%ncid, 'hmax', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%timemax_dimid/), map_file%hmax_varid)) ! time-varying maximum water depth map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%hmax_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%hmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%hmax_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%hmax_varid, 'standard_name', 'sea_floor_depth_below_sea_surface')) @@ -346,6 +363,7 @@ subroutine ncoutput_regular_map_init() ! if (store_maximum_velocity) then NF90(nf90_def_var(map_file%ncid, 'vmax', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%timemax_dimid/), map_file%vmax_varid)) ! maximum flow velocity map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%vmax_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%vmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%vmax_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%vmax_varid, 'standard_name', 'maximum_flow_velocity')) ! no standard name available @@ -357,6 +375,7 @@ subroutine ncoutput_regular_map_init() ! if (store_maximum_flux) then NF90(nf90_def_var(map_file%ncid, 'qmax', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%timemax_dimid/), map_file%qmax_varid)) ! maximum flux map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%qmax_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%qmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%qmax_varid, 'units', 'm^2 s-1')) NF90(nf90_put_att(map_file%ncid, map_file%qmax_varid, 'standard_name', 'maximum_flux')) ! no standard name available @@ -367,6 +386,7 @@ subroutine ncoutput_regular_map_init() ! if (store_cumulative_precipitation) then NF90(nf90_def_var(map_file%ncid, 'cuminf', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%timemax_dimid/), map_file%cuminf_varid)) ! time-varying cumulative infiltration map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%cuminf_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%cuminf_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%cuminf_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%cuminf_varid, 'long_name', 'cumulative_infiltration_depth')) @@ -379,6 +399,7 @@ subroutine ncoutput_regular_map_init() if (wind) then ! NF90(nf90_def_var(map_file%ncid, 'wind_u', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%wind_u_varid)) ! cumulative precipitation map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%wind_u_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, 'standard_name', 'eastward_wind')) @@ -386,6 +407,7 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, 'coordinates', 'x y')) ! NF90(nf90_def_var(map_file%ncid, 'wind_v', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%wind_v_varid)) ! cumulative precipitation map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%wind_v_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%wind_v_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%wind_v_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%wind_v_varid, 'standard_name', 'northward_wind')) @@ -394,6 +416,7 @@ subroutine ncoutput_regular_map_init() ! if (meteo3d .and. store_wind_max) then NF90(nf90_def_var(map_file%ncid, 'windmax', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid/), map_file%windmax_varid)) ! maximum wind speed + NF90(nf90_def_var_deflate(map_file%ncid, map_file%windmax_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'long_name', 'maximum_wind_speed')) @@ -406,6 +429,7 @@ subroutine ncoutput_regular_map_init() if (patmos) then ! NF90(nf90_def_var(map_file%ncid, 'surface_air_pressure', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%patm_varid)) ! cumulative precipitation map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%patm_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%patm_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%patm_varid, 'units', 'N m-2')) NF90(nf90_put_att(map_file%ncid, map_file%patm_varid, 'standard_name', 'surface_air_pressure')) @@ -417,6 +441,7 @@ subroutine ncoutput_regular_map_init() if (precip) then ! NF90(nf90_def_var(map_file%ncid, 'precipitation_rate', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%precip_varid)) ! cumulative precipitation map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%precip_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%precip_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%precip_varid, 'units', 'mm h-1')) NF90(nf90_put_att(map_file%ncid, map_file%precip_varid, 'standard_name', 'precipitation_rate')) @@ -430,6 +455,7 @@ subroutine ncoutput_regular_map_init() if (snapwave) then ! NF90(nf90_def_var(map_file%ncid, 'hm0', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%hm0_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%hm0_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, 'standard_name', 'hm0_wave_height')) @@ -437,6 +463,7 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, 'coordinates', 'x y')) ! NF90(nf90_def_var(map_file%ncid, 'hm0ig', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%hm0ig_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%hm0ig_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%hm0ig_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%hm0ig_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%hm0ig_varid, 'standard_name', 'hm0_ig_wave_height')) @@ -446,6 +473,7 @@ subroutine ncoutput_regular_map_init() if (store_wave_forces) then ! NF90(nf90_def_var(map_file%ncid, 'fwx', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%fwx_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%fwx_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, 'standard_name', 'wave_force_x')) @@ -453,6 +481,7 @@ subroutine ncoutput_regular_map_init() NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, 'coordinates', 'x y')) ! NF90(nf90_def_var(map_file%ncid, 'fwy', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%fwy_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%fwy_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%fwy_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%fwy_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%fwy_varid, 'standard_name', 'wave_force_y')) @@ -464,6 +493,7 @@ subroutine ncoutput_regular_map_init() if (wavemaker) then ! NF90(nf90_def_var(map_file%ncid, 'zsm', NF90_FLOAT, (/map_file%m_dimid, map_file%n_dimid, map_file%time_dimid/), map_file%zsm_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zsm_varid, 1, 1, nc_deflate_level)) ! deflate NF90(nf90_put_att(map_file%ncid, map_file%zsm_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zsm_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%zsm_varid, 'standard_name', 'mean_water_level')) @@ -620,7 +650,6 @@ subroutine ncoutput_quadtree_map_init() face_nodes = 0 ! two = 2 - ! nn = 0 ! @@ -659,7 +688,7 @@ subroutine ncoutput_quadtree_map_init() ! endif enddo ! - NF90(nf90_create('sfincs_map.nc', ior(NF90_CLOBBER, NF90_64BIT_OFFSET), map_file%ncid)) + NF90(nf90_create('sfincs_map.nc', ior(NF90_CLOBBER,NF90_NETCDF4), map_file%ncid)) ! TL: removed 'NF90_64BIT_OFFSET' ! ! Create dimensions ! grid, time, points @@ -694,6 +723,7 @@ subroutine ncoutput_quadtree_map_init() ! Domain ! NF90(nf90_def_var(map_file%ncid, 'mesh2d', NF90_INT, (/map_file%runtime_dimid/), map_file%mesh2d_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%mesh2d_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_varid, 'cf_role', 'mesh_topology')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_varid, 'long_name', 'Topology data of 2D network')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_varid, 'topology_dimension', 2)) @@ -706,6 +736,7 @@ subroutine ncoutput_quadtree_map_init() if (crsgeo) then ! NF90(nf90_def_var(map_file%ncid, 'mesh2d_node_x', NF90_FLOAT, (/map_file%nmesh2d_node_dimid/), map_file%mesh2d_node_x_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%mesh2d_node_x_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'units', 'degrees')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'standard_name', 'longitude')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'long_name', 'longitude')) @@ -713,6 +744,7 @@ subroutine ncoutput_quadtree_map_init() NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'location', 'node')) ! NF90(nf90_def_var(map_file%ncid, 'mesh2d_node_y', NF90_FLOAT, (/map_file%nmesh2d_node_dimid/), map_file%mesh2d_node_y_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%mesh2d_node_y_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_y_varid, 'units', 'degrees')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_y_varid, 'standard_name', 'latitude')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_y_varid, 'long_name', 'latitude')) @@ -722,6 +754,7 @@ subroutine ncoutput_quadtree_map_init() else ! NF90(nf90_def_var(map_file%ncid, 'mesh2d_node_x', NF90_FLOAT, (/map_file%nmesh2d_node_dimid/), map_file%mesh2d_node_x_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%mesh2d_node_x_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'standard_name', 'projection_x_coordinate')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'long_name', 'x-coordinate of mesh nodes')) @@ -729,6 +762,7 @@ subroutine ncoutput_quadtree_map_init() NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_x_varid, 'location', 'node')) ! NF90(nf90_def_var(map_file%ncid, 'mesh2d_node_y', NF90_FLOAT, (/map_file%nmesh2d_node_dimid/), map_file%mesh2d_node_y_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%mesh2d_node_y_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_y_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_y_varid, 'standard_name', 'projection_y_coordinate')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_node_y_varid, 'long_name', 'y-coordinate of mesh nodes')) @@ -738,6 +772,7 @@ subroutine ncoutput_quadtree_map_init() endif ! NF90(nf90_def_var(map_file%ncid, 'mesh2d_face_nodes', NF90_INT, (/map_file%max_nmesh2d_face_nodes_dimid, map_file%nmesh2d_face_dimid/), map_file%mesh2d_face_nodes_varid)) ! location of zb, zs etc. in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%mesh2d_face_nodes_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_face_nodes_varid, 'cf_role', 'face_node_connectivity')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_face_nodes_varid, 'mesh', 'mesh2d')) NF90(nf90_put_att(map_file%ncid, map_file%mesh2d_face_nodes_varid, 'location', 'face')) @@ -749,12 +784,14 @@ subroutine ncoutput_quadtree_map_init() NF90(nf90_put_att(map_file%ncid, map_file%crs_varid, 'EPSG', '-')) ! NF90(nf90_def_var(map_file%ncid, 'zb', NF90_FLOAT, (/map_file%nmesh2d_face_dimid/), map_file%zb_varid)) ! bed level in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zb_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zb_varid, 'units', 'm')) 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')) ! NF90(nf90_def_var(map_file%ncid, 'msk', NF90_INT, (/map_file%nmesh2d_face_dimid/), map_file%msk_varid)) ! input msk value in cell centre + NF90(nf90_def_var_deflate(map_file%ncid, map_file%msk_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, '_FillValue', -999)) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'units', '-')) NF90(nf90_put_att(map_file%ncid, map_file%msk_varid, 'standard_name', 'mask')) @@ -772,6 +809,7 @@ subroutine ncoutput_quadtree_map_init() ! ! Time varying map output NF90(nf90_def_var(map_file%ncid, 'zs', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%zs_varid)) ! time-varying water level map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zs_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%zs_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zs_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%zs_varid, 'standard_name', 'sea_surface_height_above_reference_level')) @@ -779,22 +817,25 @@ subroutine ncoutput_quadtree_map_init() ! if (store_velocity) then ! - NF90(nf90_def_var(map_file%ncid, 'u', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%u_varid)) ! time-varying u map + NF90(nf90_def_var(map_file%ncid, 'u', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%u_varid)) ! time-varying u map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%u_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%u_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%u_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%u_varid, 'standard_name', 'sea_water_x_velocity')) ! not truly eastward when rotated, eastward_sea_water_velocity NF90(nf90_put_att(map_file%ncid, map_file%u_varid, 'long_name', 'flow_velocity_x_direction_in_cell_edge')) ! NF90(nf90_def_var(map_file%ncid, 'v', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%v_varid)) ! time-varying u map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%v_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%v_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%v_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%v_varid, 'standard_name', 'sea_water_y_velocity')) ! not truly eastward when rotated, eastward_sea_water_velocity NF90(nf90_put_att(map_file%ncid, map_file%v_varid, 'long_name', 'flow_velocity_y_direction_in_cell_edge')) + ! endif ! ! Time varying spatial output if (store_maximum_waterlevel) then - NF90(nf90_def_var(map_file%ncid, 'timemax', NF90_FLOAT, (/map_file%timemax_dimid/), map_file%timemax_varid)) ! time + NF90(nf90_def_var(map_file%ncid, 'timemax', NF90_FLOAT, (/map_file%timemax_dimid/), map_file%timemax_varid)) ! time NF90(nf90_put_att(map_file%ncid, map_file%timemax_varid, 'units', 'seconds since ' // trim(trefstr_iso8601) )) ! time stamp following ISO 8601 NF90(nf90_put_att(map_file%ncid, map_file%timemax_varid, 'standard_name', 'time')) NF90(nf90_put_att(map_file%ncid, map_file%timemax_varid, 'long_name', 'time_in_seconds_since_' // trim(trefstr_iso8601) )) @@ -802,6 +843,7 @@ subroutine ncoutput_quadtree_map_init() ! if (store_maximum_waterlevel) then NF90(nf90_def_var(map_file%ncid, 'zsmax', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%timemax_dimid/), map_file%zsmax_varid)) ! time-varying maximum water level map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zsmax_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%zsmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zsmax_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%zsmax_varid, 'standard_name', 'maximum of sea_surface_height_above_reference_level')) @@ -810,6 +852,7 @@ subroutine ncoutput_quadtree_map_init() ! if (store_twet) then NF90(nf90_def_var(map_file%ncid, 'tmax', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%timemax_dimid/), map_file%tmax_varid)) ! time-varying duration wet cell + NF90(nf90_def_var_deflate(map_file%ncid, map_file%tmax_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%tmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%tmax_varid, 'units', 'seconds')) NF90(nf90_put_att(map_file%ncid, map_file%tmax_varid, 'standard_name', 'duration cell is considered wet')) @@ -819,6 +862,7 @@ subroutine ncoutput_quadtree_map_init() ! if (store_maximum_velocity) then NF90(nf90_def_var(map_file%ncid, 'vmax', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%timemax_dimid/), map_file%vmax_varid)) ! maximum flow velocity map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%vmax_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%vmax_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%vmax_varid, 'units', 'm s-1')) NF90(nf90_put_att(map_file%ncid, map_file%vmax_varid, 'standard_name', 'maximum_flow_velocity')) ! no standard name available @@ -833,6 +877,7 @@ subroutine ncoutput_quadtree_map_init() ! Cumulative precipitation ! NF90(nf90_def_var(map_file%ncid, 'cumprcp', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%timemax_dimid/), map_file%cumprcp_varid)) ! cumulative precipitation map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%cumprcp_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%cumprcp_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%cumprcp_varid, 'units', 'mm')) NF90(nf90_put_att(map_file%ncid, map_file%cumprcp_varid, 'long_name', 'cumulative_precipitation_depth')) @@ -841,6 +886,7 @@ subroutine ncoutput_quadtree_map_init() ! Cumulative infiltration ! NF90(nf90_def_var(map_file%ncid, 'cuminf', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%timemax_dimid/), map_file%cuminf_varid)) ! cumulative infiltration map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%cuminf_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%cuminf_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%cuminf_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%cuminf_varid, 'long_name', 'cumulative_infiltration_depth')) @@ -848,25 +894,50 @@ subroutine ncoutput_quadtree_map_init() ! endif ! - ! Store maximum wind speed - ! - if (wind .and. store_wind_max .and. meteo3d) then - NF90(nf90_def_var(map_file%ncid, 'windmax', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%timemax_dimid/), map_file%windmax_varid)) ! maximum wind speed m/s - NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, '_FillValue', FILL_VALUE)) - NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'units', 'm s-1')) - NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'long_name', 'maximum_wind_speed')) - NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'cell_methods', 'time: maximum')) + if (store_meteo) then + ! + if (wind) then + ! + NF90(nf90_def_var(map_file%ncid, 'wind_u', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%wind_u_varid)) ! time-varying wind_u map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%wind_u_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, 'units', 'm s-1')) + NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, 'standard_name', 'eastward_wind')) ! not truly eastward when rotated, eastward_sea_water_velocity + NF90(nf90_put_att(map_file%ncid, map_file%wind_u_varid, 'long_name', 'wind_speed_u')) + ! + NF90(nf90_def_var(map_file%ncid, 'wind_v', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%wind_v_varid)) ! time-varying wind_u map + NF90(nf90_def_var_deflate(map_file%ncid, map_file%wind_v_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%wind_v_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%wind_v_varid, 'units', 'm s-1')) + NF90(nf90_put_att(map_file%ncid, map_file%wind_v_varid, 'standard_name', 'northward_wind')) ! not truly eastward when rotated, eastward_sea_water_velocity + NF90(nf90_put_att(map_file%ncid, map_file%wind_v_varid, 'long_name', 'wind_speed_v')) + ! + if (store_wind_max) then + ! + ! Store maximum wind speed + ! + NF90(nf90_def_var(map_file%ncid, 'windmax', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%timemax_dimid/), map_file%windmax_varid)) ! maximum wind speed m/s + NF90(nf90_def_var_deflate(map_file%ncid, map_file%windmax_varid, 1, 1, nc_deflate_level)) + NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, '_FillValue', FILL_VALUE)) + NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'units', 'm s-1')) + NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'long_name', 'maximum_wind_speed')) + NF90(nf90_put_att(map_file%ncid, map_file%windmax_varid, 'cell_methods', 'time: maximum')) + endif + ! + endif endif ! if (snapwave) then ! NF90(nf90_def_var(map_file%ncid, 'hm0', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%hm0_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%hm0_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, 'standard_name', 'hm0_wave_height')) NF90(nf90_put_att(map_file%ncid, map_file%hm0_varid, 'long_name', 'Hm0 wave height')) ! NF90(nf90_def_var(map_file%ncid, 'hm0ig', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%hm0ig_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%hm0ig_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%hm0ig_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%hm0ig_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%hm0ig_varid, 'standard_name', 'hm0_ig_wave_height')) @@ -875,12 +946,14 @@ subroutine ncoutput_quadtree_map_init() if (store_wave_forces) then ! NF90(nf90_def_var(map_file%ncid, 'fwx', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%fwx_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%fwx_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, 'standard_name', 'wave_force_x')) NF90(nf90_put_att(map_file%ncid, map_file%fwx_varid, 'long_name', 'Wave force in x-direction')) ! NF90(nf90_def_var(map_file%ncid, 'fwy', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%fwy_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%fwy_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%fwy_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%fwy_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%fwy_varid, 'standard_name', 'wave_force_y')) @@ -891,6 +964,7 @@ subroutine ncoutput_quadtree_map_init() if (wavemaker) then ! NF90(nf90_def_var(map_file%ncid, 'zsm', NF90_FLOAT, (/map_file%nmesh2d_face_dimid, map_file%time_dimid/), map_file%zsm_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%zsm_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%zsm_varid, '_FillValue', FILL_VALUE)) NF90(nf90_put_att(map_file%ncid, map_file%zsm_varid, 'units', 'm')) NF90(nf90_put_att(map_file%ncid, map_file%zsm_varid, 'standard_name', 'filtered_water_level')) @@ -901,7 +975,8 @@ subroutine ncoutput_quadtree_map_init() endif ! if (infiltration) then - NF90(nf90_def_var(map_file%ncid, 'qinf', NF90_FLOAT, (/map_file%nmesh2d_face_dimid/), map_file%qinf_varid)) + NF90(nf90_def_var(map_file%ncid, 'qinf', NF90_FLOAT, (/map_file%nmesh2d_face_dimid/), map_file%qinf_varid)) + NF90(nf90_def_var_deflate(map_file%ncid, map_file%qinf_varid, 1, 1, nc_deflate_level)) NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, '_FillValue', FILL_VALUE)) if (inftype == 'cna') then NF90(nf90_put_att(map_file%ncid, map_file%qinf_varid, 'standard_name', 'S')) @@ -1044,7 +1119,7 @@ subroutine ncoutput_his_init() return endif ! - NF90(nf90_create('sfincs_his.nc', NF90_CLOBBER, his_file%ncid)) + NF90(nf90_create('sfincs_his.nc', ior(NF90_CLOBBER,NF90_NETCDF4), his_file%ncid)) ! ! Create dimensions ! time, stations @@ -1237,7 +1312,7 @@ subroutine ncoutput_his_init() NF90(nf90_put_att(his_file%ncid, his_file%S_varid, 'long_name', 'current moisture storage (Se) 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 CN method with recovery ! if (inftype == 'gai') then @@ -1495,8 +1570,8 @@ subroutine ncoutput_update_regular_map(t,ntmapout) n = z_index_z_n(nm) m = z_index_z_m(nm) ! - nmd1 = z_index_uv_md1(nm) - nmu1 = z_index_uv_mu1(nm) + nmd1 = z_index_uv_md(nm) + nmu1 = z_index_uv_mu(nm) uz = 0.0 if (nmd1>0) then uz = uz + 0.5*uv(nmd1) @@ -1505,8 +1580,8 @@ subroutine ncoutput_update_regular_map(t,ntmapout) uz = uz + 0.5*uv(nmu1) endif ! - ndm1 = z_index_uv_nd1(nm) - num1 = z_index_uv_nu1(nm) + ndm1 = z_index_uv_nd(nm) + num1 = z_index_uv_nu(nm) vz = 0.0 if (ndm1>0) then vz = vz + 0.5*uv(ndm1) @@ -1805,78 +1880,16 @@ subroutine ncoutput_update_quadtree_map(t,ntmapout) ! if (nm>0) then ! - if (z_flags_type(nm) == 0) then - ! - ! Regular point with four surrounding cells of the same size - ! - nmd1 = z_index_uv_md1(nm) - nmu1 = z_index_uv_mu1(nm) - ndm1 = z_index_uv_nd1(nm) - num1 = z_index_uv_mu1(nm) - uz = 0.5*(uv(nmd1) + uv(nmu1)) - vz = 0.5*(uv(ndm1) + uv(num1)) - ! - else - ! - ! Bit more complicated - ! - nmd1 = z_index_uv_md1(nm) - nmu1 = z_index_uv_mu1(nm) - ndm1 = z_index_uv_nd1(nm) - num1 = z_index_uv_mu1(nm) - nmd2 = z_index_uv_md2(nm) - nmu2 = z_index_uv_mu2(nm) - ndm2 = z_index_uv_nd2(nm) - num2 = z_index_uv_nu2(nm) - ! - ! Left - ! - u1 = 0.0 - if (nmd1>0 .and. nmd2==0) then - u1 = uv(nmd1) - elseif (nmd1>0 .and. nmd2>0) then - u1 = 0.5*uv(nmd1) + 0.5*uv(nmd2) - elseif (nmd1==0 .and. nmd2>0) then - u1 = uv(nmd2) - endif - ! - ! Right - ! - u2 = 0.0 - if (nmu1>0 .and. nmu2==0) then - u2 = uv(nmu1) - elseif (nmu1>0 .and. nmu2>0) then - u2 = 0.5*uv(nmu1) + 0.5*uv(nmu2) - elseif (nmu2==0 .and. nmu2>0) then - u2 = uv(nmu2) - endif - ! - ! Bottom - ! - v1 = 0.0 - if (ndm1>0 .and. ndm2==0) then - v1 = uv(ndm1) - elseif (ndm1>0 .and. ndm2>0) then - v1 = 0.5*uv(ndm1) + 0.5*uv(ndm2) - elseif (ndm1==0 .and. ndm2>0) then - v1 = uv(ndm2) - endif - ! - ! Right - ! - v2 = 0.0 - if (num1>0 .and. num2==0) then - v2 = uv(num1) - elseif (num1>0 .and. num2>0) then - v2 = 0.5*uv(num1) + 0.5*uv(num2) - elseif (num2==0 .and. num2>0) then - v2 = uv(num2) - endif - ! - uz = 0.5*(u1 + u2) - vz = 0.5*(v1 + v2) - ! - endif + ! Regular point with four surrounding cells of the same size + ! + n = z_index_z_n(nm) + m = z_index_z_m(nm) + nmd1 = z_index_uv_md(nm) + nmu1 = z_index_uv_mu(nm) + ndm1 = z_index_uv_nd(nm) + num1 = z_index_uv_nu(nm) + uz = 0.5*(uv(nmd1) + uv(nmu1)) + vz = 0.5*(uv(ndm1) + uv(num1)) ! utmp(nmq) = cosrot*uz - sinrot*vz vtmp(nmq) = sinrot*uz + cosrot*vz @@ -1889,6 +1902,28 @@ subroutine ncoutput_update_quadtree_map(t,ntmapout) NF90(nf90_put_var(map_file%ncid, map_file%v_varid, vtmp, (/1, ntmapout/))) ! endif + ! + if (store_meteo .and. wind) then + ! + utmp = FILL_VALUE + vtmp = FILL_VALUE + ! + do nmq = 1, quadtree_nr_points + ! + nm = index_sfincs_in_quadtree(nmq) + ! + if (nm>0) then + ! + utmp(nmq) = windu(nm) + vtmp(nmq) = windv(nm) + ! + endif + enddo + ! + NF90(nf90_put_var(map_file%ncid, map_file%wind_u_varid, utmp, (/1, ntmapout/))) + NF90(nf90_put_var(map_file%ncid, map_file%wind_v_varid, vtmp, (/1, ntmapout/))) + ! + endif ! if (snapwave) then ! @@ -2037,78 +2072,14 @@ subroutine ncoutput_update_his(t,nthisout) ! if (store_velocity) then ! - if (z_flags_type(nm) == 0) then - ! - ! Regular point with four surrounding cells of the same size - ! - nmd1 = z_index_uv_md1(nm) - nmu1 = z_index_uv_mu1(nm) - ndm1 = z_index_uv_nd1(nm) - num1 = z_index_uv_mu1(nm) - uz = 0.5*(uv(nmd1) + uv(nmu1)) - vz = 0.5*(uv(ndm1) + uv(num1)) - ! - else - ! - ! Bit more complicated - ! - nmd1 = z_index_uv_md1(nm) - nmu1 = z_index_uv_mu1(nm) - ndm1 = z_index_uv_nd1(nm) - num1 = z_index_uv_mu1(nm) - nmd2 = z_index_uv_md2(nm) - nmu2 = z_index_uv_mu2(nm) - ndm2 = z_index_uv_nd2(nm) - num2 = z_index_uv_nu2(nm) - ! - ! Left - ! - u1 = 0.0 - if (nmd1>0 .and. nmd2==0) then - u1 = uv(nmd1) - elseif (nmd1>0 .and. nmd2>0) then - u1 = 0.5*uv(nmd1) + 0.5*uv(nmd2) - elseif (nmd1==0 .and. nmd2>0) then - u1 = uv(nmd2) - endif - ! - ! Right - ! - u2 = 0.0 - if (nmu1>0 .and. nmu2==0) then - u2 = uv(nmu1) - elseif (nmu1>0 .and. nmu2>0) then - u2 = 0.5*uv(nmu1) + 0.5*uv(nmu2) - elseif (nmu2==0 .and. nmu2>0) then - u2 = uv(nmu2) - endif - ! - ! Bottom - ! - v1 = 0.0 - if (ndm1>0 .and. ndm2==0) then - v1 = uv(ndm1) - elseif (ndm1>0 .and. ndm2>0) then - v1 = 0.5*uv(ndm1) + 0.5*uv(ndm2) - elseif (ndm1==0 .and. ndm2>0) then - v1 = uv(ndm2) - endif - ! - ! Right - ! - v2 = 0.0 - if (num1>0 .and. num2==0) then - v2 = uv(num1) - elseif (num1>0 .and. num2>0) then - v2 = 0.5*uv(num1) + 0.5*uv(num2) - elseif (num2==0 .and. num2>0) then - v2 = uv(num2) - endif - ! - uz = 0.5*(u1 + u2) - vz = 0.5*(v1 + v2) - ! - endif + ! Regular point with four surrounding cells of the same size + ! + nmd1 = z_index_uv_md(nm) + nmu1 = z_index_uv_mu(nm) + ndm1 = z_index_uv_nd(nm) + num1 = z_index_uv_mu(nm) + uz = 0.5*(uv(nmd1) + uv(nmu1)) + vz = 0.5*(uv(ndm1) + uv(num1)) ! uobs(iobs) = cosrot*uz - sinrot*vz vobs(iobs) = sinrot*uz + cosrot*vz @@ -2472,12 +2443,12 @@ subroutine ncoutput_update_quadtree_max(t,ntmaxout) ! if (kcs(nm)>0) then if (subgrid) then - if ( (zs(nm) - subgrid_z_zmin(nm)) > huthresh) then - zstmp(nmq) = zs(nm) + if ( (zsmax(nm) - subgrid_z_zmin(nm)) > huthresh) then + zstmp(nmq) = zsmax(nm) endif else - if ( (zs(nm) - zb(nm)) > huthresh) then - zstmp(nmq) = zs(nm) + if ( (zsmax(nm) - zb(nm)) > huthresh) then + zstmp(nmq) = zsmax(nm) endif endif endif @@ -2625,7 +2596,7 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'outputtype_map',outputtype_map)) NF90(nf90_put_att(ncid, varid, 'outputtype_his',outputtype_his)) NF90(nf90_put_att(ncid, varid, 'bndtype',bndtype)) - NF90(nf90_put_att(ncid, varid, 'advection',iadvection)) + NF90(nf90_put_att(ncid, varid, 'advection',logical2int(advection))) NF90(nf90_put_att(ncid, varid, 'nfreqsig',nfreqsig)) NF90(nf90_put_att(ncid, varid, 'freqminig',freqminig)) NF90(nf90_put_att(ncid, varid, 'freqmaxig',freqmaxig)) @@ -2641,17 +2612,17 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'qinf_zmin',qinf_zmin)) NF90(nf90_put_att(ncid, varid, 'btfilter', btfilter)) NF90(nf90_put_att(ncid, varid, 'sfacinf', sfacinf)) - NF90(nf90_put_att(ncid, varid, 'radstr', iradstr)) - NF90(nf90_put_att(ncid, varid, 'crsgeo',igeo)) - NF90(nf90_put_att(ncid, varid, 'coriolis',icoriolis)) - NF90(nf90_put_att(ncid, varid, 'amprblock',iamprblock)) + NF90(nf90_put_att(ncid, varid, 'radstr', logical2int(radstr))) + NF90(nf90_put_att(ncid, varid, 'crsgeo',logical2int(crsgeo))) + NF90(nf90_put_att(ncid, varid, 'coriolis',logical2int(coriolis))) + NF90(nf90_put_att(ncid, varid, 'amprblock',logical2int(ampr_block))) NF90(nf90_put_att(ncid, varid, 'spwmergefrac',spw_merge_frac)) - NF90(nf90_put_att(ncid, varid, 'usespwprecip',ispwprecip)) - NF90(nf90_put_att(ncid, varid, 'global',iglobal)) + NF90(nf90_put_att(ncid, varid, 'usespwprecip',logical2int(use_spw_precip))) + NF90(nf90_put_att(ncid, varid, 'global',logical2int(global))) NF90(nf90_put_att(ncid, varid, 'nuvisc',nuvisc)) - NF90(nf90_put_att(ncid, varid, 'spinup_meteo', ispinupmeteo)) + NF90(nf90_put_att(ncid, varid, 'spinup_meteo', logical2int(spinup_meteo))) NF90(nf90_put_att(ncid, varid, 'waveage',waveage)) - NF90(nf90_put_att(ncid, varid, 'snapwave', isnapwave)) + NF90(nf90_put_att(ncid, varid, 'snapwave', logical2int(snapwave))) NF90(nf90_put_att(ncid, varid, 'wmtfilter', wmtfilter)) NF90(nf90_put_att(ncid, varid, 'wmfred',wavemaker_freduv)) NF90(nf90_put_att(ncid, varid, 'horton_kr_kd',horton_kr_kd)) @@ -2700,7 +2671,7 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'sefffile',sefffile)) 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, 'sigmafile',sigmafile)) NF90(nf90_put_att(ncid, varid, 'z0lfile',z0lfile)) NF90(nf90_put_att(ncid, varid, 'wvmfile',wvmfile)) ! @@ -2710,6 +2681,7 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'netamuamvfile',netamuamvfile)) NF90(nf90_put_att(ncid, varid, 'netamprfile',netamprfile)) NF90(nf90_put_att(ncid, varid, 'netampfile',netampfile)) + NF90(nf90_put_att(ncid, varid, 'netspwfile',netspwfile)) ! ! Output ! @@ -2724,27 +2696,27 @@ subroutine ncoutput_add_params(ncid, varid) NF90(nf90_put_att(ncid, varid, 'storetwet',storetwet)) NF90(nf90_put_att(ncid, varid, 'storehsubgrid',storehsubgrid)) NF90(nf90_put_att(ncid, varid, 'twet_threshold',twet_threshold)) - NF90(nf90_put_att(ncid, varid, 'store_tsunami_arrival_time',itsunamitime)) + NF90(nf90_put_att(ncid, varid, 'store_tsunami_arrival_time',logical2int(store_tsunami_arrival_time))) NF90(nf90_put_att(ncid, varid, 'tsunami_arrival_threshold',tsunami_arrival_threshold)) NF90(nf90_put_att(ncid, varid, 'storeqdrain',storeqdrain)) NF90(nf90_put_att(ncid, varid, 'storezvolume',storezvolume)) NF90(nf90_put_att(ncid, varid, 'writeruntime',wrttimeoutput)) - NF90(nf90_put_att(ncid, varid, 'debug',idebug)) + NF90(nf90_put_att(ncid, varid, 'debug',logical2int(debug))) NF90(nf90_put_att(ncid, varid, 'storemeteo',storemeteo)) - NF90(nf90_put_att(ncid, varid, 'storemaxwind',iwindmax)) - NF90(nf90_put_att(ncid, varid, 'storefw',istorefw)) - NF90(nf90_put_att(ncid, varid, 'storewavdir', istorewavdir)) + NF90(nf90_put_att(ncid, varid, 'storemaxwind',logical2int(store_wind_max))) + NF90(nf90_put_att(ncid, varid, 'storefw',logical2int(store_wave_forces))) + NF90(nf90_put_att(ncid, varid, 'storewavdir', logical2int(store_wave_direction))) ! NF90(nf90_put_att(ncid, varid, 'cdnrb', cd_nr)) NF90(nf90_put_att(ncid, varid, 'cdwnd', cd_wnd)) NF90(nf90_put_att(ncid, varid, 'cdval', cd_val)) ! ! Internal code switches - note, you can't store logicals in netcdf, only integers for these type of switches - NF90(nf90_put_att(ncid, varid, 'manning2d', imanning2d)) - NF90(nf90_put_att(ncid, varid, 'subgrid', isubgrid)) - NF90(nf90_put_att(ncid, varid, 'viscosity', iviscosity)) - NF90(nf90_put_att(ncid, varid, 'wavemaker', iwavemaker)) - NF90(nf90_put_att(ncid, varid, 'wavemaker_spectrum', iwavemaker_spectrum)) + NF90(nf90_put_att(ncid, varid, 'manning2d', logical2int(manning2d))) + NF90(nf90_put_att(ncid, varid, 'subgrid', logical2int(subgrid))) + NF90(nf90_put_att(ncid, varid, 'viscosity', logical2int(viscosity))) + NF90(nf90_put_att(ncid, varid, 'wavemaker', logical2int(wavemaker))) + NF90(nf90_put_att(ncid, varid, 'wavemaker_spectrum', logical2int(wavemaker_spectrum))) end subroutine ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -2767,4 +2739,21 @@ subroutine handle_err(status,file,line) end if end subroutine handle_err ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! + function logical2int(lgc) result (i) + ! + implicit none + ! + logical :: lgc + integer :: i + ! + if (lgc) then + i = 1 + else + i = 0 + endif + ! + end function + ! end module diff --git a/source/src/sfincs_obspoints.f90 b/source/src/sfincs_obspoints.f90 index 04dcfa81..c9fc617d 100644 --- a/source/src/sfincs_obspoints.f90 +++ b/source/src/sfincs_obspoints.f90 @@ -115,7 +115,7 @@ subroutine read_obs_points() ! iref = z_flags_iref(nm) ! -! write(*,'(a,i0,a,a,a,i0,a,i0,a,i0,a,i0,a,f0.3)')' Observation point ',iobs,' : "',trim(nameobs(iobs)),'" nm=',nm,' n=',n,' m=',m,' iref=',iref,' z=',zbobs(iobs) + write(*,'(a,i0,a,a,a,i0,a,i0,a,i0,a,i0,a,f0.3)')' Observation point ',iobs,' : "',trim(nameobs(iobs)),'" nm=',nm,' n=',n,' m=',m,' iref=',iref,' z=',zbobs(iobs) ! else ! diff --git a/source/src/sfincs_output.f90 b/source/src/sfincs_output.f90 index 58b63bfc..13d88f85 100644 --- a/source/src/sfincs_output.f90 +++ b/source/src/sfincs_output.f90 @@ -631,8 +631,9 @@ subroutine write_rst_file(t) ! 3: zs - ! 4: zs, q, uvmean and cnb infiltration (writing scs_Se) ! 5: zs, q, uvmean and gai infiltration (writing GA_sigma & GA_F) + ! 6: zs, q, uvmean and hor infiltration (writing rain_T1) ! - ! Write for Infiltration methods (rsttype 4 or 5) + ! Write for Infiltration methods (rsttype 4 or 5 or 6) if (inftype == 'cnb' .or. inftype == 'gai' .or. inftype == 'hor') then ! if (inftype == 'cnb') then @@ -643,22 +644,22 @@ subroutine write_rst_file(t) write(911)uvmean write(911)scs_Se ! - elseif (inftype == 'hor') then + elseif (inftype == 'gai') then ! - write(911)4 + write(911)5 write(911)zs write(911)q write(911)uvmean - write(911)rain_T1 + write(911)GA_sigma + write(911)GA_F ! - elseif (inftype == 'gai') then + elseif (inftype == 'hor') then ! - write(911)5 + write(911)6 write(911)zs write(911)q write(911)uvmean - write(911)GA_sigma - write(911)GA_F + write(911)rain_T1 ! endif ! diff --git a/source/src/sfincs_spiderweb.f90 b/source/src/sfincs_spiderweb.f90 index 6d1881b0..d396851b 100644 --- a/source/src/sfincs_spiderweb.f90 +++ b/source/src/sfincs_spiderweb.f90 @@ -139,7 +139,8 @@ subroutine read_spw_file(filename,nt,nrows,ncols,spwrad,time,xe,ye,vmag,vdir,pdr close(888) ! end subroutine - + ! + ! Read amu subroutine read_amuv_file(filename,nt,nrows,ncols,time,uv,trefstr) ! implicit none diff --git a/source/src/sfincs_wavemaker.f90 b/source/src/sfincs_wavemaker.f90 index 4dcba6e3..84566be5 100644 --- a/source/src/sfincs_wavemaker.f90 +++ b/source/src/sfincs_wavemaker.f90 @@ -152,6 +152,7 @@ subroutine read_wavemaker_polylines() ! ! Check right and above ! +! nmu = z_index_uv_mu(ip) ! index of 1st uv neighbor to the right nmu = z_index_uv_mu1(ip) ! index of 1st uv neighbor to the right ! if (nmu>0) then @@ -753,6 +754,8 @@ subroutine read_wavemaker_polylines() iwm = 0 nok = 0 ! + write(*,*)'Setting wave makers ...' + ! do ip = 1, np ! if (kcs(ip)==4) then @@ -1124,6 +1127,8 @@ subroutine read_wavemaker_polylines() ! ! Set flags for kcuv points ! + write(*,*)'Setting wave maker flags ...' + ! do iwm = 1, wavemaker_nr_uv_points ! ip = wavemaker_index_uv(iwm) @@ -1137,14 +1142,16 @@ subroutine read_wavemaker_polylines() ! do ip = 1, npuv ! - if (kcuv(uv_index_u_nmd(ip)) == 4) uv_flags_adv(ip) = 0 !QUESTION - TL: goes wrong if you don't specify a quadtree grid !? > should make this clear in documentation - if (kcuv(uv_index_u_nmu(ip)) == 4) uv_flags_adv(ip) = 0 - if (kcuv(uv_index_u_num(ip)) == 4) uv_flags_adv(ip) = 0 - if (kcuv(uv_index_u_ndm(ip)) == 4) uv_flags_adv(ip) = 0 - if (kcuv(uv_index_v_ndm(ip)) == 4) uv_flags_adv(ip) = 0 - if (kcuv(uv_index_v_nm(ip)) == 4) uv_flags_adv(ip) = 0 - if (kcuv(uv_index_v_nmu(ip)) == 4) uv_flags_adv(ip) = 0 - if (kcuv(uv_index_v_ndmu(ip)) == 4) uv_flags_adv(ip) = 0 + ! TODO: re-implement uv_flags_adv ? + ! +! if (kcuv(uv_index_u_nmd(ip)) == 4) uv_flags_adv(ip) = 0 +! if (kcuv(uv_index_u_nmu(ip)) == 4) uv_flags_adv(ip) = 0 +! if (kcuv(uv_index_u_num(ip)) == 4) uv_flags_adv(ip) = 0 +! if (kcuv(uv_index_u_ndm(ip)) == 4) uv_flags_adv(ip) = 0 +! if (kcuv(uv_index_v_ndm(ip)) == 4) uv_flags_adv(ip) = 0 +! if (kcuv(uv_index_v_nm(ip)) == 4) uv_flags_adv(ip) = 0 +! if (kcuv(uv_index_v_nmu(ip)) == 4) uv_flags_adv(ip) = 0 +! if (kcuv(uv_index_v_ndmu(ip)) == 4) uv_flags_adv(ip) = 0 ! enddo endif @@ -1308,6 +1315,8 @@ subroutine read_wavemaker_polylines() ! call RANDOM_NUMBER(r) dphiig(ifreq) = 1.0e-6*2*3.1416/freqig(ifreq) enddo + ! + write(*,*)'Wavemakers initialized !' ! end subroutine diff --git a/source/src/snapwave/snapwave_boundaries.f90 b/source/src/snapwave/snapwave_boundaries.f90 index 9db6e2e7..21c0fa8b 100644 --- a/source/src/snapwave/snapwave_boundaries.f90 +++ b/source/src/snapwave/snapwave_boundaries.f90 @@ -55,7 +55,7 @@ subroutine read_boundary_enclosure () ! ! Read wave boundaries ! - write(*,*)'Reading wave boundary enclosure ...' + write(*,*)'Reading wave boundary enclosure ', trim(encfile), ' ...' open(500, file=trim(encfile)) !as in bwvfile of SFINCS do while(.true.) read(500,*,iostat = stat)dummy diff --git a/source/third_party_open/netcdf/netcdf4.vcxproj b/source/third_party_open/netcdf/netcdf4.vcxproj index 9a454c79..977f1623 100644 --- a/source/third_party_open/netcdf/netcdf4.vcxproj +++ b/source/third_party_open/netcdf/netcdf4.vcxproj @@ -43,7 +43,7 @@ StaticLibrary MultiByte - v143 + v142