From b01858d8887ebdb44f295e813a6fde7608d9955b Mon Sep 17 00:00:00 2001 From: Yujia Cai Date: Tue, 16 Jun 2026 11:36:54 -0700 Subject: [PATCH 1/9] Factor out advective flux calculation DSM2-1419 --- dsm2/src/transport/advection.f90 | 170 ++++++++++++++++++++++++++----- 1 file changed, 147 insertions(+), 23 deletions(-) diff --git a/dsm2/src/transport/advection.f90 b/dsm2/src/transport/advection.f90 index c4526fd1..464b7d84 100644 --- a/dsm2/src/transport/advection.f90 +++ b/dsm2/src/transport/advection.f90 @@ -153,6 +153,152 @@ subroutine advect(mass, & ncell, & nvar) + call advective_flux(flow, & + flow_prev, & + flow_lo, & + flow_hi, & + area, & + area_prev, & + area_lo, & + area_hi, & + explicit_diffuse_op, & + ncell, & + nvar, & + time, & + dt, & + dx, & + use_limiter, & + width, & + width_prev, & + depth, & + depth_prev, & + wet_p, & + wet_p_prev, & + constraint, & + name, & + LL, & + sed_percent, & + div_flux, & + source_prev) + + !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 + end subroutine + + subroutine advective_flux(flow, & + flow_prev, & + flow_lo, & + flow_hi, & + area, & + area_prev, & + area_lo, & + area_hi, & + explicit_diffuse_op, & + ncell, & + nvar, & + time, & + dt, & + dx, & + use_limiter, & + width, & + width_prev, & + depth, & + depth_prev, & + wet_p, & + wet_p_prev, & + constraint, & + name, & + LL, & + sed_percent, & + div_flux, & + source_prev) + + use constants + use primitive_variable_conversion + use gradient + 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 + integer, intent(in) :: ncell !< Number of cells + integer, intent(in) :: nvar !< Number of variables + integer, intent(in) :: LL !< Time step + real(gtm_real),intent(out) :: source_prev(ncell,nvar) + real(gtm_real),intent(out) :: div_flux(ncell,nvar) !< cell centered flux divergence, time centered + real(gtm_real),intent(out) :: sed_percent(n_node,n_qext,nvar)! difference call conc_gradient(grad_lo, & @@ -265,29 +411,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 From f5fe4d770f84520ac57513eec45b867cf4d56809 Mon Sep 17 00:00:00 2001 From: Yujia Cai Date: Tue, 16 Jun 2026 11:37:59 -0700 Subject: [PATCH 2/9] Calculate "expected" diffusive flux and write to hdf5 DSM2-1419 --- dsm2/src/gtm/gtm.f90 | 30 +++++- dsm2/src/process_io/gtm_hdf_ts_write.f90 | 18 ++++ dsm2/src/transport/advection.f90 | 19 +++- dsm2/src/transport/mass_balance.f90 | 112 +++++++++++++++++++++++ 4 files changed, 172 insertions(+), 7 deletions(-) create mode 100644 dsm2/src/transport/mass_balance.f90 diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index f218f98f..04855b36 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,9 @@ 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 :: 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(:,:,:)! +! 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, & + conc, & + conc_prev, & + area, & + area_prev, & + ncell, & + nvar, & + advective_div_flux, & + 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 + real(gtm_real),intent(in) :: area(ncell) !< cell-centered area at new time + real(gtm_real),intent(in) :: area_prev(ncell) !< cell-centered area at old time + real(gtm_real),intent(in) :: conc(ncell,nvar) !< concentration at current time + real(gtm_real),intent(in) :: conc_prev(ncell,nvar) !< concentration at previous 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 + !-----locals + real(gtm_real):: calc_div_flux(ncell,nvar) !< cell centered flux divergence, time centered + real(gtm_real):: mass_change(ncell,nvar) !< actual mass change calculated from conc*area change between new and old time step + real(gtm_real):: source_prev(ncell,nvar) !< source term at previous time step + real(gtm_real):: sed_percent(n_node,n_qext,nvar)! Date: Mon, 29 Jun 2026 14:21:59 -0700 Subject: [PATCH 3/9] Write out advective divergence to hdf5 --- dsm2/src/gtm/gtm.f90 | 9 +++++++++ dsm2/src/process_io/gtm_hdf_ts_write.f90 | 18 ++++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index 04855b36..4d5ca7ea 100644 --- a/dsm2/src/gtm/gtm.f90 +++ b/dsm2/src/gtm/gtm.f90 @@ -88,6 +88,7 @@ module dsm2gtm ! 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(:,:,:)! Date: Mon, 29 Jun 2026 14:22:52 -0700 Subject: [PATCH 4/9] Fix diffusive divergence calculated from mass balance --- dsm2/src/gtm/gtm.f90 | 4 --- dsm2/src/transport/mass_balance.f90 | 53 +++++------------------------ 2 files changed, 8 insertions(+), 49 deletions(-) diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index 4d5ca7ea..140c46cf 100644 --- a/dsm2/src/gtm/gtm.f90 +++ b/dsm2/src/gtm/gtm.f90 @@ -682,10 +682,6 @@ subroutine gtm_loop() call expected_net_diffusive_flux(mass, & mass_prev, & - conc, & - conc_prev, & - area, & - area_prev, & n_cell, & n_var, & sub_gtm_time_step*sixty, & diff --git a/dsm2/src/transport/mass_balance.f90 b/dsm2/src/transport/mass_balance.f90 index 61e8023d..88ec859d 100644 --- a/dsm2/src/transport/mass_balance.f90 +++ b/dsm2/src/transport/mass_balance.f90 @@ -29,12 +29,10 @@ module mass_balance subroutine expected_net_diffusive_flux(mass, & mass_prev, & - conc, & - conc_prev, & - area, & - area_prev, & - ncell, & - nvar, & + ncell, & + nvar, & + dt, & + dx, & advective_div_flux, & calc_net_diffusive_flux) @@ -54,59 +52,24 @@ subroutine expected_net_diffusive_flux(mass, & !--- args real(gtm_real),intent(out) :: calc_net_diffusive_flux(ncell,nvar) !< net diffusive flux at new time - real(gtm_real),intent(in) :: area(ncell) !< cell-centered area at new time - real(gtm_real),intent(in) :: area_prev(ncell) !< cell-centered area at old time - real(gtm_real),intent(in) :: conc(ncell,nvar) !< concentration at current time - real(gtm_real),intent(in) :: conc_prev(ncell,nvar) !< concentration at previous 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 !-----locals - real(gtm_real):: calc_div_flux(ncell,nvar) !< cell centered flux divergence, time centered real(gtm_real):: mass_change(ncell,nvar) !< actual mass change calculated from conc*area change between new and old time step - real(gtm_real):: source_prev(ncell,nvar) !< source term at previous time step - real(gtm_real):: sed_percent(n_node,n_qext,nvar)! Date: Mon, 6 Jul 2026 13:41:38 -0700 Subject: [PATCH 5/9] Remove mass adjustment in advection step for calculating diffusive divergence --- dsm2/src/gtm/gtm.f90 | 3 +++ dsm2/src/transport/mass_balance.f90 | 4 +++- dsm2/src/transport/state_variables.f90 | 6 ++++++ 3 files changed, 12 insertions(+), 1 deletion(-) diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index 140c46cf..c4f81422 100644 --- a/dsm2/src/gtm/gtm.f90 +++ b/dsm2/src/gtm/gtm.f90 @@ -638,7 +638,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 ---------- @@ -687,6 +689,7 @@ subroutine gtm_loop() 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) diff --git a/dsm2/src/transport/mass_balance.f90 b/dsm2/src/transport/mass_balance.f90 index 88ec859d..387bf0b4 100644 --- a/dsm2/src/transport/mass_balance.f90 +++ b/dsm2/src/transport/mass_balance.f90 @@ -34,6 +34,7 @@ subroutine expected_net_diffusive_flux(mass, & dt, & dx, & advective_div_flux, & + mass_adv_adjust, & calc_net_diffusive_flux) use constants @@ -59,6 +60,7 @@ subroutine expected_net_diffusive_flux(mass, & 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 @@ -66,7 +68,7 @@ subroutine expected_net_diffusive_flux(mass, & 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) + calc_net_diffusive_flux(icell,ivar) = - mass_change(icell,ivar)/(dt/dx(icell)) - advective_div_flux(icell,ivar) + mass_adv_adjust(icell,ivar) end do end do return 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) From fdee23903133d26daffbffcfec62131bdc393c19 Mon Sep 17 00:00:00 2001 From: Yujia Cai Date: Mon, 6 Jul 2026 13:42:06 -0700 Subject: [PATCH 6/9] Fix the wrong message --- dsm2/src/gtm/gtm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index c4f81422..d1901bdf 100644 --- a/dsm2/src/gtm/gtm.f90 +++ b/dsm2/src/gtm/gtm.f90 @@ -476,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 From ea1225c34967d510bd31b1cefa6cb947b5b5eba0 Mon Sep 17 00:00:00 2001 From: Yujia Cai Date: Wed, 8 Jul 2026 13:00:15 -0700 Subject: [PATCH 7/9] Add missing dt/dx --- dsm2/src/transport/mass_balance.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsm2/src/transport/mass_balance.f90 b/dsm2/src/transport/mass_balance.f90 index 387bf0b4..36e06961 100644 --- a/dsm2/src/transport/mass_balance.f90 +++ b/dsm2/src/transport/mass_balance.f90 @@ -68,7 +68,7 @@ subroutine expected_net_diffusive_flux(mass, & 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) + 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 From bfe2cec18ec141928963923bb55efd52536915d5 Mon Sep 17 00:00:00 2001 From: Yujia Cai Date: Thu, 9 Jul 2026 12:02:49 -0700 Subject: [PATCH 8/9] Wrtie number of subcycles to hdf5 when debug_print is true DM-238 --- dsm2/src/gtm/gtm.f90 | 2 +- dsm2/src/process_io/gtm_hdf_ts_write.f90 | 82 +++++++++++++++++++++++- 2 files changed, 82 insertions(+), 2 deletions(-) diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index d1901bdf..817fb126 100644 --- a/dsm2/src/gtm/gtm.f90 +++ b/dsm2/src/gtm/gtm.f90 @@ -506,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 @@ -804,6 +803,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, & diff --git a/dsm2/src/process_io/gtm_hdf_ts_write.f90 b/dsm2/src/process_io/gtm_hdf_ts_write.f90 index 398aff0e..99163744 100644 --- a/dsm2/src/process_io/gtm_hdf_ts_write.f90 +++ b/dsm2/src/process_io/gtm_hdf_ts_write.f90 @@ -42,6 +42,7 @@ 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 @@ -181,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 @@ -766,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, & @@ -1223,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) From 530795729de7f266647b57cc48787951f77e4e12 Mon Sep 17 00:00:00 2001 From: Kijin Nam Date: Mon, 13 Jul 2026 11:50:31 -0700 Subject: [PATCH 9/9] Guard non-advective flux calculation Non-advective flux calculation is only for debug information recording, and it is guarded to make it run only when needed. --- dsm2/src/gtm/gtm.f90 | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/dsm2/src/gtm/gtm.f90 b/dsm2/src/gtm/gtm.f90 index 817fb126..1b07c3bb 100644 --- a/dsm2/src/gtm/gtm.f90 +++ b/dsm2/src/gtm/gtm.f90 @@ -681,17 +681,19 @@ subroutine gtm_loop() call prim2cons(mass,conc,area,n_cell,n_var) end if - 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) + 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