diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index f218f98f..1b07c3bb 100644 --- a/dsm2/src/gtm/gtm.f90 +++ b/dsm2/src/gtm/gtm.f90 @@ -48,6 +48,7 @@ module dsm2gtm use gradient_adjust use advection use diffusion + use mass_balance use dispersion_coefficient use source_sink use boundary_advection @@ -86,6 +87,10 @@ module dsm2gtm real(gtm_real), allocatable :: mass_closure(:) ! real(gtm_real) :: theta = half !< Crank-Nicolson implicitness coeficient real(gtm_real) :: constant_dispersion + real(gtm_real), allocatable :: advective_div_flux(:,:) + real(gtm_real), allocatable :: advective_div_flux_ts(:) !< time series of advective flux + real(gtm_real), allocatable :: calc_diffusive_div_flux(:,:) !< calculated net diffusive flux using mass change and advective flux + real(gtm_real), allocatable, target :: calc_diffusive_div_flux_ts(:) !< calculated net diffusive flux time series of the 1st water quality constituent real(gtm_real), allocatable :: sed_percent(:,:,:)! Initialize qual tide file for sub-timestep cycle time series + subroutine init_sub_ts_gtm_hdf5(hdf_file, ntime) + use hdf5 + implicit none + type(gtm_hdf_t), intent(inout) :: hdf_file + integer(HID_T) :: cparms ! dataset creation property identifier + integer :: error ! HDF5 Error flag + integer :: res_rank = 1 + integer(HSIZE_T), dimension(1) :: chunk_dims = 0 ! Dataset dimensions + integer :: ntime ! number of time points in tidefile + integer :: nqext + integer :: nconc + integer(HSIZE_T), dimension(1) :: file_dims = 0 ! Data size on file + integer(HID_T) :: fspace_id ! File space identifier + + !-------Create the datasets + res_rank = 1 + file_dims(1) = ntime + + chunk_dims(1) = min(TIME_CHUNK,ntime) + + ! Add chunking and compression + call h5pcreate_f(H5P_DATASET_CREATE_F, cparms, error) + if (ntime .gt. MIN_STEPS_FOR_CHUNKING) then + call h5pset_chunk_f(cparms, res_rank, chunk_dims, error) + call H5Pset_szip_f (cparms, H5_SZIP_NN_OM_F, & + HDF_SZIP_PIXELS_PER_BLOCK, error); + end if + + call h5screate_simple_f(res_rank, & + file_dims, & + fspace_id, & + error) + call h5dcreate_f(hdf_file%data_id, & + "number of sub-timesteps", & + H5T_NATIVE_INTEGER, & + fspace_id, & + hdf_file%sub_ts_id, & + error, & + cparms) + call add_timeseries_attributes(hdf_file%sub_ts_id, & + hdf_file%start_julmin, & + hdf_file%write_interval) + + return + end subroutine + + !> Write time series data to Qual tidefile (dimension cell) subroutine write_gtm_hdf(hdf_file, & @@ -1187,6 +1273,36 @@ subroutine write_1D_string_array(dest_id, name, arr, strlen, nstr) end subroutine + !> Write one integer value at a given time index into a pre-allocated 1D integer dataset + subroutine write_gtm_hdf_ts_int(dset_id, val, time_index) + use hdf5 + implicit none + integer(HID_T), intent(in) :: dset_id !< pre-allocated dataset identifier + integer, intent(in) :: val !< integer value to write + integer, intent(in) :: time_index !< 0-based time index + integer :: rank + integer(HID_T) :: fspace_id + integer(HID_T) :: memspace_id + integer(HSIZE_T), dimension(1) :: mdata_dims + integer(HSIZE_T), dimension(1) :: subset_dims + integer(HSIZE_T), dimension(1) :: h_offset + integer :: error + + h_offset(1) = time_index + subset_dims(1) = 1 + mdata_dims(1) = 1 + rank = 1 + call h5screate_simple_f(rank, mdata_dims, memspace_id, error) + call h5dget_space_f(dset_id, fspace_id, error) + call h5sselect_hyperslab_f(fspace_id, H5S_SELECT_SET_F, h_offset, subset_dims, error) + call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, val, mdata_dims, error, memspace_id, fspace_id) + call verify_error(error, "sub-timestep integer time series write") + call h5sclose_f(fspace_id, error) + call h5sclose_f(memspace_id, error) + return + end subroutine + + !> add time series attributes to tidefile subroutine add_timeseries_attributes(dset_id, ts_start, ts_interval) diff --git a/dsm2/src/transport/advection.f90 b/dsm2/src/transport/advection.f90 index c4526fd1..c7a9f83e 100644 --- a/dsm2/src/transport/advection.f90 +++ b/dsm2/src/transport/advection.f90 @@ -66,7 +66,8 @@ subroutine advect(mass, & constraint, & name, & LL, & - sed_percent) + sed_percent, & + div_flux) use constants use primitive_variable_conversion @@ -116,6 +117,7 @@ subroutine advect(mass, & real(gtm_real), intent(in) :: wet_p(ncell) real(gtm_real), intent(in) :: wet_p_prev(ncell) real(gtm_real), intent(in) :: constraint(ncell,nvar) + real(gtm_real), intent(out) :: div_flux(ncell,nvar) !< cell centered flux divergence, time centered character(len=32), intent(in) :: name(nvar) !-----locals @@ -131,7 +133,6 @@ subroutine advect(mass, & real(gtm_real) :: grad(ncell,nvar) !< cell centered difference adujsted for boundaries and hydraulic devices real(gtm_real) :: flux_lo(ncell,nvar) !< flux on lo side of cell, time centered real(gtm_real) :: flux_hi(ncell,nvar) !< flux on hi side of cell, time centered - real(gtm_real) :: div_flux(ncell,nvar) !< cell centered flux divergence, time centered real(gtm_real) :: mass_closure(ncell) ! real(gtm_real) :: sed_percent(n_node,n_qext,nvar)! difference call conc_gradient(grad_lo, & @@ -265,29 +420,7 @@ subroutine advect(mass, & flux_hi, & ncell, & nvar) - - !Conservative update including source. - call update_conservative(mass, & - mass_prev, & - div_flux, & - explicit_diffuse_op, & - source_prev, & - flow, & - area, & - area_prev, & - width, & - depth, & - wet_p, & - constraint, & - name, & - ncell, & - nvar, & - time, & - dt, & - dx) - !accumulative mass closure error - mass_closure =( (area - area_prev) + (dt/dx)*(flow_hi-flow_lo) ) - return + return end subroutine diff --git a/dsm2/src/transport/mass_balance.f90 b/dsm2/src/transport/mass_balance.f90 new file mode 100644 index 00000000..36e06961 --- /dev/null +++ b/dsm2/src/transport/mass_balance.f90 @@ -0,0 +1,77 @@ +! +! Copyright (C) 2017 State of California, +! Department of Water Resources. +! This file is part of DSM2-GTM. +! +! The Delta Simulation Model 2 (DSM2) - General Transport Model (GTM) +! is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! DSM2 is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with DSM2. If not, see . +! + +!> Module orchestrating the advection scheme. The main +!> routine in the module is advection(). +!>@ingroup transport +module mass_balance + + contains + !> Perform mass balance to calculate "predicted" diffusive flux at cell level + !> predicted diffusive flux = actual mass change - advective flux + + subroutine expected_net_diffusive_flux(mass, & + mass_prev, & + ncell, & + nvar, & + dt, & + dx, & + advective_div_flux, & + mass_adv_adjust, & + calc_net_diffusive_flux) + + use constants + use IO_Units + use primitive_variable_conversion + use gradient + use advection + use source_sink + use boundary_advection + use boundary_concentration + use gradient_adjust + use boundary_concentration + use common_vars, only: n_node, n_qext + + implicit none + + !--- args + real(gtm_real),intent(out) :: calc_net_diffusive_flux(ncell,nvar) !< net diffusive flux at new time + integer, intent(in) :: ncell !< Number of cells + integer, intent(in) :: nvar !< Number of variables + real(gtm_real),intent(in) :: mass(ncell,nvar) !< mass at new time + real(gtm_real),intent(in) :: mass_prev(ncell,nvar) !< mass at old time + real(gtm_real),intent(in) :: advective_div_flux(ncell,nvar) !< advective flux divergence + real(gtm_real),intent(in) :: dt !< current time step from old time to new time + real(gtm_real),intent(in) :: dx(ncell) !< spatial step + real(gtm_real),intent(in) :: mass_adv_adjust(ncell,nvar) !< mass adjustment at the end of the advection step + !-----locals + real(gtm_real):: mass_change(ncell,nvar) !< actual mass change calculated from conc*area change between new and old time step + integer :: icell, ivar + + mass_change = mass(:,:) - mass_prev(:,:) + do icell = 1,ncell + do ivar = 1,nvar + calc_net_diffusive_flux(icell,ivar) = - mass_change(icell,ivar)/(dt/dx(icell)) - advective_div_flux(icell,ivar) + mass_adv_adjust(icell,ivar)/(dt/dx(icell)) + end do + end do + return + end subroutine + +end module \ No newline at end of file diff --git a/dsm2/src/transport/state_variables.f90 b/dsm2/src/transport/state_variables.f90 index dcf1b4f6..afaff2ac 100644 --- a/dsm2/src/transport/state_variables.f90 +++ b/dsm2/src/transport/state_variables.f90 @@ -29,6 +29,12 @@ module state_variables !> Mass of constituent in the current/new time step, !> dimensions (ncell, nvar) real(gtm_real), save, allocatable :: mass(:,:) + !> Mass of constituent before the mass adjustment at the end of the advection step, + !> dimensions (ncell, nvar) + real(gtm_real), save, allocatable :: mass_before_adjust(:,:) + !> Change in mass adjust at the end of the advection step + !> dimensions (ncell, nvar) + real(gtm_real), save, allocatable :: mass_adjust_advet(:,:) !> Mass of constituent in the previous time step, !> dimensions (ncell, nvar)