Skip to content

implementing multiple tracers #598

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

Open
nmizukami wants to merge 28 commits into
base: main
Choose a base branch
from
Open

implementing multiple tracers #598

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
9e90da0
start adding tracer dimension in arrays
nmizukami Dec 16, 2025
5ec7638
allowing to specify multiple tracers in control file and reading trac…
nmizukami Jan 21, 2026
be9d232
fix tracer arrays in reach flux data structure and allcation
nmizukami Apr 4, 2026
0494d7c
finalized routing of multiple tracers
nmizukami Apr 4, 2026
a98c83e
additional dimension, and pio-domain decomp for multiple tracers
nmizukami Apr 4, 2026
06a79f5
including multiple tracers in history file writing
nmizukami Apr 4, 2026
ac2666d
mpi communication for mulitple tracers
nmizukami Apr 4, 2026
e68be55
restart writing for multiple tracers
nmizukami Apr 4, 2026
f5dbb63
finalize restart read including multiple tracers
nmizukami Apr 5, 2026
511b3b3
I think this is more efficient than using source in allocation
nmizukami Apr 5, 2026
7a0bc9a
consistent data initialization for both coldstart and restart
nmizukami Apr 5, 2026
9e3b55f
fix allocation bugs so that restart works properly
nmizukami Apr 5, 2026
cf5433b
small cleanup
nmizukami Apr 6, 2026
1983211
if tracer is off, need to ignore this
nmizukami Apr 6, 2026
2a67a5c
remove comm and pio variables not used in subroutines
nmizukami Apr 7, 2026
5f22297
basinSolute array allocation and initialization bugfix. Still not per…
nmizukami Apr 17, 2026
0729960
solute concentration array needs to be 2D
nmizukami May 4, 2026
ba41bdd
add sub-reach concentration in tracer state variables
nmizukami May 6, 2026
1643a22
add a capability to track solute movement within a reach to avoid ins…
nmizukami May 6, 2026
03364ce
add writing sub-reach tracer concentration in tracer states
nmizukami May 6, 2026
a800da2
adding the entire domain mass balance check
nmizukami May 6, 2026
23fe3f3
adding 2d double precision array scatter
nmizukami May 6, 2026
b2fcaec
add reading of sub-reach concentration as tracer state variable to be…
nmizukami May 6, 2026
c2d53a5
adding tracer sub-reach concentration scatter
nmizukami May 6, 2026
e405e8f
bugfix
nmizukami May 17, 2026
80a0a8b
typo
nmizukami May 17, 2026
7bfd9ca
for headwater reach, discharge from reach need to be just REACH_Q
nmizukami May 19, 2026
35e70c3
Merge branch 'main' into v3.0.1-26-g85bbe71f_tracer
nmizukami May 19, 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
4 changes: 2 additions & 2 deletions route/build/cpl/RtmMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,7 @@ SUBROUTINE route_ini()
! 6. initialize state variables
!-------------------------------------------------------

call init_state_data(iam, npes, mpicom_rof, ierr, cmessage)
call init_state_data(iam, npes, ierr, cmessage)
if(ierr/=0)then; call shr_sys_abort(trim(subname)//trim(cmessage)); endif

! put reach flux variables to associated HRUs
Expand Down Expand Up @@ -578,7 +578,7 @@ SUBROUTINE route_run(rstwr)
do ns = 1,nsub
call t_startf('mizuRoute_subcycling')

call mpi_route(iam, npes, mpicom_rof, ierr, cmessage, scatter_ro=.false.)
call mpi_route(iam, npes, ierr, cmessage, scatter_ro=.false.)
if(ierr/=0)then; call shr_sys_abort(trim(subname)//trim(cmessage)); endif

call t_stopf('mizuRoute_subcycling')
Expand Down
28 changes: 17 additions & 11 deletions route/build/src/basinUH.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ SUBROUTINE hru_irf(iSeg, & ! input: index of segment to be processed
ierr, message) ! output: error control

USE globalData, ONLY: FRAC_FUTURE !
USE globalData, ONLY: nTracer ! number of tracers
USE public_var, ONLY: is_lake_sim ! logical whether or not lake should be simulated
USE public_var, ONLY: tracer ! logical whether or not tracer is on

Expand All @@ -84,6 +85,7 @@ SUBROUTINE hru_irf(iSeg, & ! input: index of segment to be processed
integer(i4b), intent(out) :: ierr ! error code
character(*), intent(out) :: message ! error message
! Local variables
integer(i4b) :: iTrace ! loop index
real(dp), allocatable :: FRAC_FUTURE_local(:) ! local FRAC_FUTURE so that it can be changed for lakes to impulse
integer(i4b) :: ntdh ! number of time steps in IRF
character(len=strLen) :: cmessage ! error message from subroutine
Expand All @@ -100,10 +102,12 @@ SUBROUTINE hru_irf(iSeg, & ! input: index of segment to be processed
RCHFLX_out(iSeg)%QFUTURE(:) = 0._dp
end if

if (.not.allocated(RCHFLX_out(iSeg)%solute_future) .and. tracer)then
ntdh = size(FRAC_FUTURE)
allocate(RCHFLX_out(iSeg)%solute_future(ntdh), source=0._dp, stat=ierr)
if(ierr/=0)then; message=trim(message)//'unable to allocate space for RCHFLX_out(segIndex)%solute_future'; return; endif
if (tracer)then
if (.not.allocated(RCHFLX_out(iSeg)%solute_future))then
ntdh = size(FRAC_FUTURE)
allocate(RCHFLX_out(iSeg)%solute_future(ntdh,nTracer), source=0._dp, stat=ierr)
if(ierr/=0)then; message=trim(message)//'unable to allocate RCHFLX_out(segIndex)%solute_future'; return; endif
end if
end if

allocate(FRAC_FUTURE_local, source=FRAC_FUTURE, stat=ierr)
Expand All @@ -125,13 +129,15 @@ SUBROUTINE hru_irf(iSeg, & ! input: index of segment to be processed
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif

if(tracer) then
! perform river network UH routing
call irf_conv(FRAC_FUTURE_local, & ! input: unit hydrograph
RCHFLX_out(iSeg)%BASIN_solute_inst, & ! input: upstream fluxes
RCHFLX_out(iSeg)%solute_future, & ! inout: updated solute future time series
RCHFLX_out(iSeg)%BASIN_solute, & ! inout: updated fluxes at reach
ierr, message) ! output: error control
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
do iTrace=1,nTracer
! perform river network UH routing
call irf_conv(FRAC_FUTURE_local, & ! input: unit hydrograph
RCHFLX_out(iSeg)%BASIN_solute_inst(iTrace), & ! input: upstream fluxes
RCHFLX_out(iSeg)%solute_future(:,iTrace), & ! inout: updated solute future time series
RCHFLX_out(iSeg)%BASIN_solute(iTrace), & ! inout: updated fluxes at reach
ierr, message) ! output: error control
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
end do
end if

END SUBROUTINE hru_irf
Expand Down
14 changes: 7 additions & 7 deletions route/build/src/dataTypes.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@ MODULE dataTypes
real(dp) , allocatable :: basinRunoff(:) ! remapped river network catchment runoff [depth/time] (size: number of nHRU)
real(dp) , allocatable :: basinEvapo(:) ! remapped river network catchment evaporation [depth/time] (size: number of nHRU)
real(dp) , allocatable :: basinPrecip(:) ! remapped river network catchment precipitation [depth/time] (size: number of nHRU)
real(dp) , allocatable :: basinSolute(:) ! remapped river network catchment solute in water [mass/time] (size: number of nHRU)
real(dp) , allocatable :: basinSolute(:,:)! remapped river network catchment solute in water [mass/time] (size: nHRU, nTracer)
end type runoff

type, public, extends(inputData) :: wm ! water-management
Expand Down Expand Up @@ -318,7 +318,7 @@ MODULE dataTypes
real(dp), allocatable :: Q(:) ! Discharge at sub-reaches at current step (m3/s)
real(dp), allocatable :: A(:) ! Flow area at sub-reach at current step (m2)
real(dp), allocatable :: H(:) ! Flow height at sub-reach at current step (m)
real(dp), allocatable :: Solute_mass(:) ! solute mass at sub-reach at current step (mg)
real(dp), allocatable :: c_solute(:,:) ! solute concentration at sub-reach at current step (mg/m3)
end type SUBRCH

type, public :: kwRch
Expand Down Expand Up @@ -352,8 +352,8 @@ MODULE dataTypes
real(dp) :: REACH_WM_FLUX_actual ! water management fluxes to and from each reach [m3/s]
real(dp) :: WB ! reach water balance error [m3]
real(dp) :: Qerror ! simulated discharge error compared to obs [m3/s] -- only for data assimilation
real(dp) :: reach_solute_mass(0:1) ! constituent mass in channel [mg]
real(dp) :: reach_solute_flux ! constituent mass flux from reach outlet [mg/s]
real(dp),allocatable :: reach_solute_mass(:,:) ! mass of constituent(s) in channel [mg]
real(dp),allocatable :: reach_solute_flux(:) ! mass flux of constituent(s) from reach outlet [mg/s]
end type hydraulic

! fluxes and states in each reach
Expand All @@ -362,9 +362,9 @@ MODULE dataTypes
real(dp), allocatable :: QFUTURE_IRF(:) ! runoff volume in future time steps for IRF routing [m3/s]
real(dp), allocatable :: QPASTUP_IRF(:,:) ! runoff volume in the past time steps for lake upstream [m3/s]
real(dp), allocatable :: DEMANDPAST_IRF(:,:) ! demand volume for lake [m3/s]
real(dp), allocatable :: solute_future(:) ! lateral solute mass in future time steps [mg/s]
real(dp) :: BASIN_solute ! instantaneous constituent mass from the local basin [mg/s]
real(dp) :: BASIN_solute_inst ! instantaneous constituent mass from the local basin [mg/s]
real(dp), allocatable :: solute_future(:,:) ! lateral solute mass in future time steps [mg/s]
real(dp), allocatable :: BASIN_solute(:) ! mass of instantaneous constituent(s) from the local basin [mg/s]
real(dp), allocatable :: BASIN_solute_inst(:) ! mass of instantaneous constituent(s) from the local basin [mg/s]
real(dp) :: BASIN_QI ! instantaneous runoff volume from the local basin [m3/s]
real(dp) :: BASIN_QR(0:1) ! routed runoff volume from the local basin [m3/s]
type(hydraulic), allocatable :: ROUTE(:) ! reach fluxes and states for each routing method
Expand Down
62 changes: 34 additions & 28 deletions route/build/src/globalData.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,34 +140,37 @@ MODULE globalData

! ---------- MPI/OMP/PIO variables ----------------------------------------------------------------

integer(i4b), public :: mpicom_route ! communicator for this program
integer(i4b), public :: pid ! process id
integer(i4b), public :: nNodes ! number of MPI processors
integer(i4b), public :: nThreads ! number of threads
logical(lgt), public :: masterproc ! root logical. root processor => true, other => false
logical(lgt), public :: multiProcs ! MPI multi-processors logical. => number of processors>1 true, other => false
integer(i4b), public :: pio_numiotasks = -99 ! Number of iotasks (ntasks/stride) - see PIO documentation for more information
integer(i4b), public :: pio_rearranger = 2 ! 0=>PIO_rearr_none 1=> PIO_rearr_box 2=> PIO_rearr_subset
integer(i4b), public :: pio_root = 1 ! PIO root
integer(i4b), public :: pio_stride = 1 ! PIO stride - see PIO documentation for more information
type(iosystem_desc_t), pointer, public :: pioSystem ! PIO I/O system data
integer(i4b), public :: mpicom_route ! communicator for this program
integer(i4b), public :: pid ! process id
integer(i4b), public :: nNodes ! number of MPI processors
integer(i4b), public :: nThreads ! number of threads
logical(lgt), public :: masterproc ! root logical. root processor => true, other => false
logical(lgt), public :: multiProcs ! MPI multi-processors logical. => number of processors>1 true, other => false
integer(i4b), public :: pio_numiotasks = -99 ! Number of iotasks (ntasks/stride) - see PIO documentation for more information
integer(i4b), public :: pio_rearranger = 2 ! 0=>PIO_rearr_none 1=> PIO_rearr_box 2=> PIO_rearr_subset
integer(i4b), public :: pio_root = 1 ! PIO root
integer(i4b), public :: pio_stride = 1 ! PIO stride - see PIO documentation for more information
type(iosystem_desc_t), pointer, public :: pioSystem ! PIO I/O system data
! pio decomposition used for history file variables
type(io_desc_t), public :: ioDesc_hru_float ! [hru] (float)
type(io_desc_t), public :: ioDesc_rch_float ! [reach] (float)
type(io_desc_t), public :: ioDesc_gauge_float ! [gauge points] (float)
type(io_desc_t), public :: ioDesc_hru_float ! [hru] (float)
type(io_desc_t), public :: ioDesc_rch_float ! [reach] (float)
type(io_desc_t), public :: ioDesc_gauge_float ! [gauge points] (float)
! decomposition used for restart variables
type(io_desc_t), public :: ioDesc_rch_int ! [reach] (integer)
type(io_desc_t), public :: ioDesc_rch_double ! [reach] (double precision)
type(io_desc_t), public :: ioDesc_hru_double ! [hru] (double precision)
type(io_desc_t), public :: ioDesc_wave_int ! [reach x max. number of waves] (integer)
type(io_desc_t), public :: ioDesc_wave_double ! [reach x max. number of waves] (double precision)
type(io_desc_t), public :: ioDesc_mesh_kw_double ! [reach x Euler kinematic wave computational meshes] (double precision)
type(io_desc_t), public :: ioDesc_mesh_mc_double ! [reach x Muskingum-Cunge computational meshes] (double precision)
type(io_desc_t), public :: ioDesc_mesh_dw_double ! [reach x Duffusive wave computational meshes] (double precision)
type(io_desc_t), public :: ioDesc_irf_double ! [reach x IRF future timer-steps] (double precision)
type(io_desc_t), public :: ioDesc_vol_double ! [reach x time-bounds] (double precision)
type(io_desc_t), public :: ioDesc_irf_bas_double ! [reach x basin UH future time-steps] (double precision)
integer(i4b), allocatable, public :: index_write_gage(:) ! reach indices to gauge points w.r.t. distributed domains
type(io_desc_t), public :: ioDesc_rch_int ! [reach] (integer)
type(io_desc_t), public :: ioDesc_rch_double ! [reach] (double precision)
type(io_desc_t), public :: ioDesc_hru_double ! [hru] (double precision)
type(io_desc_t), public :: ioDesc_wave_int ! [reach, max number of waves] (integer)
type(io_desc_t), public :: ioDesc_wave_double ! [reach, max number of waves] (double precision)
type(io_desc_t), public :: ioDesc_mesh_kw_double ! [reach, Euler kinematic wave computational meshes] (double precision)
type(io_desc_t), public :: ioDesc_mesh_mc_double ! [reach, Muskingum-Cunge computational meshes] (double precision)
type(io_desc_t), public :: ioDesc_mesh_dw_double ! [reach, Duffusive wave computational meshes] (double precision)
type(io_desc_t), public :: ioDesc_tracer_bas_double ! [reach, basin UH future time-steps, tracers] (double precision)
type(io_desc_t), public :: ioDesc_tracer_double ! [reach, tracers] (double precision)
type(io_desc_t), public :: ioDesc_tracer_mesh_double ! [reach, sub_reach, tracers] (double precision)
type(io_desc_t), public :: ioDesc_irf_double ! [reach, IRF future timer-steps] (double precision)
type(io_desc_t), public :: ioDesc_vol_double ! [reach, time-bounds] (double precision)
type(io_desc_t), public :: ioDesc_irf_bas_double ! [reach, basin UH future time-steps] (double precision)
integer(i4b), allocatable, public :: index_write_gage(:) ! reach indices to gauge points w.r.t. distributed domains

! ---------- conversion factors -------------------------------------------------------------------
! runoff and solute mass flux unit conversion. Time conversion factor is also needed for solute because the time unit of solute flux MAY be different than runoff's
Expand Down Expand Up @@ -252,8 +255,11 @@ MODULE globalData
real(dp), allocatable, public :: basinEvapo_main(:) ! HRU evaporation array (m/s) for mainstem
real(dp), allocatable, public :: basinPrecip_trib(:) ! HRU precipitation array (m/s) for tributaries
real(dp), allocatable, public :: basinPrecip_main(:) ! HRU precipitation array (m/s) for mainstem
real(dp), allocatable, public :: basinSolute_trib(:) ! HRU constituent array (mg/s) for tributaries
real(dp), allocatable, public :: basinSolute_main(:) ! HRU constituent array (mg/s) for mainstem

! solute (tracer) data
integer(i4b) :: nTracer ! number of tracers
real(dp), allocatable, public :: basinSolute_trib(:,:) ! HRU constituent array [nhru, nTracer] (mg/s) for tributaries
real(dp), allocatable, public :: basinSolute_main(:,:) ! HRU constituent array [nhru, nTracer] (mg/s) for mainstem

! seg water management fluxes and target volume
type(wm), public :: wm_data ! SEG flux and target vol data structure for one time step for river network
Expand Down
Loading