Skip to content

implementing multiple tracer #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

Draft
nmizukami wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

implementing multiple tracer #598

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d6ba771
start adding tracer dimension in arrays
nmizukami Dec 16, 2025
f4151be
Merge branch 'ESCOMP:main' into v3.0.1-26-g85bbe71f_tracer
nmizukami Jan 15, 2026
746db40
allowing to specify multiple tracers in control file and reading trac…
nmizukami Jan 21, 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
16 changes: 9 additions & 7 deletions route/build/src/basinUH.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,13 +125,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 iTracer=1,nTracer
! perform river network UH routing
call irf_conv(FRAC_FUTURE_local, & ! input: unit hydrograph
RCHFLX_out(iSeg)%BASIN_solute_inst(iTracer), & ! input: upstream fluxes
RCHFLX_out(iSeg)%solute_future, & ! inout: updated solute future time series
RCHFLX_out(iSeg)%BASIN_solute(iTracer), & ! 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
10 changes: 5 additions & 5 deletions route/build/src/dataTypes.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,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 @@ -357,8 +357,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(0:1,:) ! 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 @@ -368,8 +368,8 @@ MODULE dataTypes
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 :: 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
7 changes: 5 additions & 2 deletions route/build/src/globalData.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,8 +252,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
37 changes: 22 additions & 15 deletions route/build/src/main_route.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ MODULE main_route_module
SUBROUTINE main_route(basinRunoff_in, & ! input: basin (i.e.,HRU) runoff (m/s)
basinEvapo_in, & ! input: basin (i.e.,HRU) evaporation (m/s)
basinPrecip_in, & ! input: basin (i.e.,HRU) precipitation (m/s)
basinSolute_in, & ! input: basin (i.e.,HRU) constituent mass fluw (mg/s/m2)
basinSolute_in, & ! input: basin (i.e.,HRU) mass flux of constituent(s) (mg/s/m2)
reachflux_in, & ! input: reach (i.e.,reach) flux (m3/s)
reachvol_in, & ! input: reach (i.e.,reach) target volume for lakes (m3)
ixRchProcessed, & ! input: indices of reach to be routed
Expand All @@ -52,7 +52,8 @@ SUBROUTINE main_route(basinRunoff_in, & ! input: basin (i.e.,HRU) runoff (m/s)
USE globalData, ONLY: TSEC ! beginning/ending of simulation time step [sec]
USE globalData, ONLY: simDatetime ! current model datetime
USE globalData, ONLY: rch_routes !
USE public_var, ONLY: doesBasinRoute
USE globalData, ONLY: nTracer ! number of active tracers
USE public_var, ONLY: doesBasinRoute ! hillslope routing option
USE public_var, ONLY: is_flux_wm ! logical whether or not fluxes should be passed
USE public_var, ONLY: is_vol_wm ! logical whether or not target volume should be passed
USE public_var, ONLY: qmodOption ! options for streamflow modification (DA)
Expand All @@ -62,7 +63,7 @@ SUBROUTINE main_route(basinRunoff_in, & ! input: basin (i.e.,HRU) runoff (m/s)
real(dp), allocatable, intent(in) :: basinRunoff_in(:) ! basin (i.e.,HRU) runoff (m/s)
real(dp), allocatable, intent(in) :: basinEvapo_in(:) ! basin (i.e.,HRU) evaporation (m/s)
real(dp), allocatable, intent(in) :: basinPrecip_in(:) ! basin (i.e.,HRU) precipitation (m/s)
real(dp), allocatable, intent(in) :: basinSolute_in(:) ! basin (i.e.,HRU) constituent (mg/s/m2)
real(dp), allocatable, intent(in) :: basinSolute_in(:,:) ! basin (i.e.,HRU) constituent(s) (mg/s/m2)
real(dp), allocatable, intent(in) :: reachflux_in(:) ! reach (i.e.,reach) flux (m3/s)
real(dp), allocatable, intent(in) :: reachvol_in(:) ! reach (i.e.,reach) target volume for lakes (m3)
integer(i4b), allocatable, intent(in) :: ixRchProcessed(:) ! indices of reach to be routed
Expand All @@ -78,11 +79,12 @@ SUBROUTINE main_route(basinRunoff_in, & ! input: basin (i.e.,HRU) runoff (m/s)
character(len=strLen) :: cmessage ! error message of downwind routine
real(dp) :: qobs ! observed discharge at a time step and site
real(dp), allocatable :: reachRunoff_local(:) ! reach runoff (m3/s)
real(dp), allocatable :: reachSolute_local(:) ! reach constituent (mg/s)
real(dp), allocatable :: reachSolute_local(:,:) ! reach constituent(s) (mg/s)
real(dp), allocatable :: reachEvapo_local(:) ! reach evaporation (m3/s)
real(dp), allocatable :: reachPrecip_local(:) ! reach precipitation (m3/s)
integer(i4b) :: nSeg ! number of reach to be processed
integer(i4b) :: iSeg ! index of reach
integer(i4b) :: iTracer ! index of tracers
integer(i4b) :: ix,jx ! loop index
integer(i4b), allocatable :: reach_ix(:)
integer(i4b), parameter :: no_mod=0
Expand Down Expand Up @@ -149,16 +151,17 @@ SUBROUTINE main_route(basinRunoff_in, & ! input: basin (i.e.,HRU) runoff (m/s)
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif

if(tracer) then
allocate(reachSolute_local(nSeg), source=0._dp, stat=ierr)
allocate(reachSolute_local(nSeg,nTracer), source=0._dp, stat=ierr)
if(ierr/=0)then; message=trim(message)//'problem allocating arrays for [reachSolute_local]'; return; endif

call basin2reach_mass(basinSolute_in, & ! input: basin total constitunet (mg/s/m2)
NETOPO_in, & ! input: reach topology
RPARAM_in, & ! input: reach parameter
reachSolute_local, & ! output:total constituent going to reach (mg/s)
ierr, cmessage, & ! output: error control
ixRchProcessed) ! optional input: indices of reach to be routed
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
do iTracer=1,nTracer
call basin2reach_mass(basinSolute_in(:,iTracer), & ! input: basin total constitunet (mg/s/m2)
NETOPO_in, & ! input: reach topology
RPARAM_in, & ! input: reach parameter
reachSolute_local(:,iTracer), & ! output:total constituent going to reach (mg/s)
ierr, cmessage, & ! output: error control
ixRchProcessed) ! optional input: indices of reach to be routed
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
end do
end if

if (is_lake_sim) then
Expand Down Expand Up @@ -196,7 +199,9 @@ SUBROUTINE main_route(basinRunoff_in, & ! input: basin (i.e.,HRU) runoff (m/s)
RCHFLX_out(ixRchProcessed(iSeg))%BASIN_QI = reachRunoff_local(iSeg)
if (tracer) then
if (RCHFLX_out(ixRchProcessed(iSeg))%BASIN_QI>0) then ! this may cause mass inbalance between input and output. so may need to feed flag to land surface model
RCHFLX_out(ixRchProcessed(iSeg))%BASIN_solute_inst = reachSolute_local(iSeg)
do iTracer=1,nTracer
RCHFLX_out(ixRchProcessed(iSeg))%BASIN_solute_inst(iTracer) = reachSolute_local(iSeg,iTracer)
end do
else
RCHFLX_out(ixRchProcessed(iSeg))%BASIN_solute_inst = 0._dp
endif
Expand All @@ -218,7 +223,9 @@ SUBROUTINE main_route(basinRunoff_in, & ! input: basin (i.e.,HRU) runoff (m/s)
if (tracer) then
do iSeg = 1,nSeg
if (RCHFLX_out(ixRchProcessed(iSeg))%BASIN_QR(1)>0) then
RCHFLX_out(ixRchProcessed(iSeg))%BASIN_solute = reachSolute_local(iSeg) ! total constituent going to reach (mg/s)
do iTracer=1,nTracer
RCHFLX_out(ixRchProcessed(iSeg))%BASIN_solute(iTracer) = reachSolute_local(iSeg,iTracer) ! total constituent going to reach (mg/s)
end do
else
RCHFLX_out(ixRchProcessed(iSeg))%BASIN_solute = 0._dp
endif
Expand Down
32 changes: 18 additions & 14 deletions route/build/src/mpi_process.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,7 @@ SUBROUTINE comm_ntopo_data(pid, & ! input: proc id
USE globalData, ONLY: vol_wm_trib ! nRch target vol holder for mainstem
USE globalData, ONLY: nRch_trib ! number of tributary reaches
USE globalData, ONLY: nHRU_trib ! number of tributary HRUs
USE globalData, ONLY: nTracer ! number of tracer
USE globalData, ONLY: hru_per_proc ! number of hrus assigned to each proc (size = num of procs+1)
USE globalData, ONLY: rch_per_proc ! number of reaches assigned to each proc (size = num of procs+1)
USE globalData, ONLY: ixHRU_order ! global HRU index in the order of proc assignment (size = total number of HRUs contributing to any reaches, nContribHRU)
Expand Down Expand Up @@ -574,10 +575,10 @@ SUBROUTINE comm_ntopo_data(pid, & ! input: proc id
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif

if (tracer) then
allocate(basinSolute_main(nHRU_mainstem), source=0.0_dp, stat=ierr, errmsg=cmessage)
allocate(basinSolute_main(nHRU_mainstem,nTracer), source=0.0_dp, stat=ierr, errmsg=cmessage)
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif

allocate(basinSolute_trib(nHRU_trib), source=0.0_dp, stat=ierr, errmsg=cmessage)
allocate(basinSolute_trib(nHRU_trib,nTracer), source=0.0_dp, stat=ierr, errmsg=cmessage)
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
end if

Expand Down Expand Up @@ -629,7 +630,7 @@ SUBROUTINE comm_ntopo_data(pid, & ! input: proc id
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif

if (tracer) then
allocate(basinSolute_trib(nHRU_trib), source=0.0_dp, stat=ierr, errmsg=cmessage)
allocate(basinSolute_trib(nHRU_trib,nTracer), source=0.0_dp, stat=ierr, errmsg=cmessage)
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
end if

Expand Down Expand Up @@ -1394,6 +1395,7 @@ SUBROUTINE scatter_runoff(nNodes, & ! mpi variables: number nodes

USE globalData, ONLY: nHRU ! number of all HRUs
USE globalData, ONLY: nHRU_mainstem ! number of mainstem HRUs
USE globalData, ONLY: nTracer ! number of tracers
USE globalData, ONLY: runoff_data ! runoff data structure
USE globalData, ONLY: hru_per_proc ! number of hrus assigned to each proc (i.e., node)
USE globalData, ONLY: basinRunoff_main ! HRU runoff holder for mainstem
Expand All @@ -1413,10 +1415,11 @@ SUBROUTINE scatter_runoff(nNodes, & ! mpi variables: number nodes
integer(i4b), intent(out) :: ierr ! error code
character(len=strLen), intent(out) :: message ! error message
! local variables
integer(i4b) :: iTracer ! tracer loop index
real(dp) :: basinRunoff_local(nHRU) ! temporal basin runoff (m/s) for whole domain
real(dp) :: basinEvapo_local(nHRU) ! temporal basin evaporation (m/s) for whole domain
real(dp) :: basinPrecip_local(nHRU) ! temporal basin precipitation (m/s) for whole domain
real(dp) :: basinSolute_local(nHRU) ! temporal basin constituent (g) for whole domain
real(dp) :: basinSolute_local(nHRU,nTracer) ! temporal basin constituent (mg) for whole domain
character(len=strLen) :: cmessage ! error message from a subroutine

ierr=0; message='scatter_runoff/'
Expand All @@ -1433,7 +1436,7 @@ SUBROUTINE scatter_runoff(nNodes, & ! mpi variables: number nodes
end if

if (tracer) then
basinSolute_local (1:nHRU) = runoff_data%basinSolute(1:nHRU)
basinSolute_local (1:nHRU,1:nTracer) = runoff_data%basinSolute(1:nHRU,1:nTracer)
end if

if (nHRU_mainstem>0) then
Expand All @@ -1459,11 +1462,10 @@ SUBROUTINE scatter_runoff(nNodes, & ! mpi variables: number nodes

if (tracer) then
if (.not. allocated(basinSolute_main)) then
allocate(basinSolute_main(nHRU_mainstem), stat=ierr)
if(ierr/=0)then; message=trim(message)//'problem allocating array for [basinEvapo_main]'; return; endif
allocate(basinSolute_main(nHRU_mainstem,nTracer), source=0.0_dp, stat=ierr, errmsg=cmessage)
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
endif
basinSolute_main (1:nHRU_mainstem) = basinSolute_local (1:nHRU_mainstem)

basinSolute_main (1:nHRU_mainstem, 1:nTracer) = basinSolute_local(1:nHRU_mainstem,1:nTracer)
end if
end if

Expand Down Expand Up @@ -1491,11 +1493,13 @@ SUBROUTINE scatter_runoff(nNodes, & ! mpi variables: number nodes
end if

if (tracer) then
call shr_mpi_scatterV(basinSolute_local(nHRU_mainstem+1:nHRU), &
hru_per_proc(0:nNodes-1), &
basinSolute_trib, &
ierr, cmessage)
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
do iTracer=1,nTracer
call shr_mpi_scatterV(basinSolute_local(nHRU_mainstem+1:nHRU,iTracer), &
hru_per_proc(0:nNodes-1), &
basinSolute_trib(:,iTracer), &
ierr, cmessage)
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
end do
end if

END SUBROUTINE scatter_runoff
Expand Down
Loading