Skip to content

Dsm2 1419 #40

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
kjnam merged 9 commits into from
Jul 13, 2026
Merged

Dsm2 1419 #40

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
b01858d
Factor out advective flux calculation
YujiaCai23 Jun 16, 2026
f5fe4d7
Calculate "expected" diffusive flux and write to hdf5
YujiaCai23 Jun 16, 2026
596d762
Write out advective divergence to hdf5
YujiaCai23 Jun 29, 2026
1997e85
Fix diffusive divergence calculated from mass balance
YujiaCai23 Jun 29, 2026
ca7644d
Remove mass adjustment in advection step for calculating diffusive di…
YujiaCai23 Jul 6, 2026
fdee239
Fix the wrong message
YujiaCai23 Jul 6, 2026
ea1225c
Add missing dt/dx
YujiaCai23 Jul 8, 2026
bfe2cec
Wrtie number of subcycles to hdf5 when debug_print is true
YujiaCai23 Jul 9, 2026
5307957
Guard non-advective flux calculation
kjnam Jul 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 40 additions & 4 deletions dsm2/src/gtm/gtm.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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(:,:,:)!<percentages of compositions at boundaries & 10 is the maximum number of
!external flows !<TODO: make array dimensions effective
! local variables to obtain time series data from HDF5 file
Expand Down Expand Up @@ -256,10 +261,16 @@ subroutine gtm_prepare_loop()
allocate(explicit_diffuse_op(n_cell,n_var))
allocate(linear_decay(n_var))
allocate(cfl(n_cell))
allocate(advective_div_flux(n_cell,n_var))
allocate(advective_div_flux_ts(n_cell))
allocate(calc_diffusive_div_flux_ts(n_cell))
allocate(disp_coef_lo(n_cell), disp_coef_hi(n_cell))
allocate(disp_coef_lo_prev(n_cell), disp_coef_hi_prev(n_cell))
allocate(calc_diffusive_div_flux(n_cell,n_var))
allocate(mass_closure(n_cell))

calc_diffusive_div_flux_ts = zero

write(*,*) "You need to have ",n_cell," number of cells in initial file."

call assign_node_ts
Expand Down Expand Up @@ -465,7 +476,7 @@ subroutine gtm_loop()
if ((max_cfl .gt. one).and.(sub_time_step)) then
if (ceil_max_cfl .gt. max_num_sub_ts) then
ceil_max_cfl = max_num_sub_ts
write(*,*) "exceed max number of sub timestep. consider to decrease the cell size."
write(*,*) "exceed max number of sub timestep. consider to increase the cell size."
end if
sub_gtm_time_step = gtm_time_interval/dble(ceil_max_cfl)
call deallocate_network_tmp
Expand Down Expand Up @@ -495,7 +506,6 @@ subroutine gtm_loop()
sub_st = 1
sub_gtm_time_step = gtm_time_interval
end if

do st = 1, sub_st
t_index = int(t_in_slice/sub_gtm_time_step) + st
new_current_time = current_time + dble(st-1)*sub_gtm_time_step
Expand Down Expand Up @@ -612,7 +622,8 @@ subroutine gtm_loop()
constraint, &
constituents(:)%use_module, &
LL, &
sed_percent)
sed_percent, &
advective_div_flux)
if (nonnegative) then
where (mass.lt.zero) mass = zero
end if
Expand All @@ -626,7 +637,9 @@ subroutine gtm_loop()
! call cons2prim(conc, mass, area, n_cell, n_var)
!elseif:
call cons2prim(conc, mass, area-mass_closure, n_cell, n_var) ! calculate conc based on the corrected mass.
mass_before_adjust = mass(:,:)
call prim2cons(mass, conc, area, n_cell, n_var) ! recalculate the mass based on the original area and the new conc: this mass imblance is caused by hydro and sub time step.
mass_adjust_advet = mass - mass_before_adjust
!endif

!--------- Diffusion ----------
Expand Down Expand Up @@ -668,6 +681,20 @@ subroutine gtm_loop()
call prim2cons(mass,conc,area,n_cell,n_var)
end if

if (debug_print == .true.) then
call expected_net_diffusive_flux(mass, &
mass_prev, &
n_cell, &
n_var, &
sub_gtm_time_step*sixty, &
dx_arr, &
advective_div_flux, &
mass_adjust_advet, &
calc_diffusive_div_flux )
advective_div_flux_ts = advective_div_flux(:,1)
calc_diffusive_div_flux_ts = calc_diffusive_div_flux(:,1)
end if

mass_prev = mass
conc_prev = conc
conc_resv_prev = conc_resv
Expand All @@ -688,7 +715,7 @@ subroutine gtm_loop()
prev_node_conc = node_conc
prev_conc_stip = conc_stip
node_conc = LARGEREAL
end do ! end of loop of sub time step
end do ! end of loop of sub time step
! at time step 1, do the above calculation to get old time and half time variables, but still keep the initial concentrations
if (run_pdaf) then
if (current_time .eq. gtm_start_jmin) then
Expand Down Expand Up @@ -778,6 +805,7 @@ subroutine gtm_loop()
time_index_in_gtm_hdf)
end if
if (debug_print==.true.) then
call write_gtm_hdf_ts_int(gtm_hdf%sub_ts_id, sub_st, time_index_in_gtm_hdf)
call write_gtm_hdf_ts(gtm_hdf%cell_flow_id, &
flow, &
n_cell, &
Expand All @@ -790,6 +818,14 @@ subroutine gtm_loop()
cfl, &
n_cell, &
time_index_in_gtm_hdf)
call write_gtm_hdf_ts(gtm_hdf%calc_net_diffusive_flux_id, &
calc_diffusive_div_flux_ts, &
n_cell, &
time_index_in_gtm_hdf)
call write_gtm_hdf_ts(gtm_hdf%net_advective_flux_id, &
advective_div_flux_ts, &
n_cell, &
time_index_in_gtm_hdf)
end if
end if
end if
Expand Down
118 changes: 117 additions & 1 deletion dsm2/src/process_io/gtm_hdf_ts_write.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,8 +42,11 @@ module gtm_hdf_ts_write
integer(HID_T) :: qext_conc_id
integer(HID_T) :: chan_budget_id
integer(HID_T) :: cell_flow_id
integer(HID_T) :: sub_ts_id
integer(HID_T) :: cell_area_id
integer(HID_T) :: cell_cfl_id
integer(HID_T) :: net_advective_flux_id
integer(HID_T) :: calc_net_diffusive_flux_id
integer(HSIZE_T) :: conc_dim
integer(HSIZE_T) :: cell_dim
integer(HSIZE_T) :: chan_dim
Expand Down Expand Up @@ -179,6 +182,7 @@ subroutine init_gtm_hdf(hdf_file, &
end if
! for debug print only
if (debug_print .eq. .true.) then
call init_sub_ts_gtm_hdf5(hdf_file, ntime)
call init_cell_gtm_hdf5_debug(hdf_file, ncell, ntime)
end if

Expand Down Expand Up @@ -299,6 +303,8 @@ subroutine init_cell_gtm_hdf5_debug(hdf_file, ncell, ntime)
integer(HID_T) :: fspace_flow_id ! File space identifier
integer(HID_T) :: fspace_area_id ! File space identifier
integer(HID_T) :: fspace_cfl_id ! File space identifier
integer(HID_T) :: fspace_net_advective_flux_id ! File space identifier
integer(HID_T) :: fspace_calc_net_diffusive_flux_id ! File space identifier

integer(HID_T) :: cparms ! dataset creatation property identifier
integer :: error ! HDF5 Error flag
Expand Down Expand Up @@ -376,6 +382,38 @@ subroutine init_cell_gtm_hdf5_debug(hdf_file, ncell, ntime)
call add_timeseries_attributes(hdf_file%cell_cfl_id, &
hdf_file%start_julmin, &
hdf_file%write_interval)
! initialize table for cell net advective flux
call h5screate_simple_f(cell_rank, &
cell_file_dims, &
fspace_net_advective_flux_id, &
error)
call h5dcreate_f(hdf_file%data_id, &
"cell net advective flux", &
H5T_NATIVE_DOUBLE, &
fspace_net_advective_flux_id, &
hdf_file%net_advective_flux_id, &
error, &
cparms)
call verify_error(error,"Cell net advective flux dataset creation")
call add_timeseries_attributes(hdf_file%net_advective_flux_id, &
hdf_file%start_julmin, &
hdf_file%write_interval)
! initialize table for cell net diffusive flux
call h5screate_simple_f(cell_rank, &
cell_file_dims, &
fspace_calc_net_diffusive_flux_id, &
error)
call h5dcreate_f(hdf_file%data_id, &
"cell net diffusive flux", &
H5T_NATIVE_DOUBLE, &
fspace_calc_net_diffusive_flux_id, &
hdf_file%calc_net_diffusive_flux_id, &
error, &
cparms)
call verify_error(error,"Cell net diffusive flux dataset creation")
call add_timeseries_attributes(hdf_file%calc_net_diffusive_flux_id, &
hdf_file%start_julmin, &
hdf_file%write_interval)
return
end subroutine

Expand Down Expand Up @@ -730,13 +768,61 @@ subroutine write_gtm_hdf_ts(data_id, &
error, &
memspace_id, &
fspace_id)
call verify_error(error,"Cell time series write (write_gtm_hdf_ts)")
call verify_error(error,"Cell time series write (write_gtm_hdf_ts)")
call h5sclose_f (fspace_id, error)
call h5sclose_f (memspace_id, error)
end if
return
end subroutine

!> 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, &
Expand Down Expand Up @@ -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)

Expand Down
Loading