Add river routing module for offline (GEOSldas) simulations

[weiyuan-jiang]

Add river routing module for offline land modeling (GEOSldas). Replaces #1023

---

Left for future PRs:

• Run routing module within the AGCM and coupled model (AOGCM).
• Run only the routing module (but not Catchment) in GEOSldas.
• Include landice tiles in routing (GEOSldas and GCM).

Related PRs:

• Add river routing GEOSldas#850
• Add river routing grid comp GEOSldas_GridComp#145
• added pfaf index to MAPL_Locstream MAPL#4000 (merged 5 Dec 2025; included in MAPL v2.64)

---

Testing:

• TBD

---

List of tasks that remain to be addressed (in no particular order):

• Move Routing GC from Land GC to Surface GC (@weiyuan-jiang)
• Make sure Routing works when GEOSldas is run in cube-sphere tile space (@weiyuan-jiang, @zyj8881357)
• Add routing parameters into Route GC restart files (similar to, say, porosity in Catch restart) (@zyj8881357)
• Clean up files in Utils/Raster/preproc/routing_model/ (@zyj8881357)
• Remove GEOSroute_GridComp/offline/ and save files elsewhere (for reference when we address the option to run only the routing module but not Catchment in GEOSldas) (@zyj8881357)
• Add files from bcs preprocessing (river_inputs.bc, cell_area.nc?) into shared bcs directory and update location of these inputs in make_bcs scripts (@biljanaorescanin, @zyj8881357)
• What (if anything) is needed for river routing in Utils/remap_restart? (@weiyuan-jiang, @zyj8881357) Add river-routing restart assembeling GEOS_Util#170

---

@YujiN, I am trying to bring back some codes in develop branch. The develop branch is not working as it is but I believe it would be more efficient. The runoff and the other variables are distributed as needed, which is way more efficient than allgatherV.

[gmao-rreichle]

@zyj8881357 : I'm still struggling to understand what exactly is done in EASE_pfaf_frac.F90.

[zyj8881357]

I’ve changed pfaf_frac to subtile_area in the *_tile2pfaf.nc file. This means the model no longer needs to read the catchment area from the restart file; instead, it can be derived by summing the subtile areas during initialization. This works for non-EASE grids as well.

[gmao-rreichle]

@zyj8881357 : We're getting close to being done, but I still have a few questions and suggestions, see below. Please let me know if anything is unclear.

[gmao-rreichle]

@zyj8881357 : I tried some more to understand what this program is doing, and I think I got further than last time. But I'm still nowhere near understanding enough to meaningfully review this code. I'll need more time (and perhaps some help from you) to work through this.

[gmao-rreichle]

@weiyuan-jiang : Does the mapping from tile space to Pfaf catch space work when the simulation domain is not global?
And what if the tile space contains landice tiles?

[zyj8881357]

@weiyuan-jiang Could you please help with this? Exit the model program when the simulation domain is not global. Thank you!

[zyj8881357]

Land ice and lakes (as defined in the *.rst raster file) are excluded from the catchment area in the routing model.

[mathomp4]

Suggested change

	real(REAL64), parameter :: PI = 3.14159265358979323846 ! ( MAPL_PI=3.14159265358979323846d0 as of Feb 2026 )
	real(REAL64), parameter :: PI = 3.14159265358979323846d0 ! ( MAPL_PI=3.14159265358979323846d0 as of Feb 2026 )

@gmao-rreichle et al, I think we want the d0 here. I talked with @tclune and technically Fortran says without the d0 the long constant is truncated to REAL32 and then cast to REAL64. So this is the "right" way.

NOTE: @tclune also said that in his experience, no compilers actually did this when built optimized, but might when building Debug.

[gmao-rreichle]

@zyj8881357, @weiyuan-jiang : I added a commit (93fef6e) and two comments in the hope of clarifying what is done and what we may want to do re. the EASEtile-to-Pfafcatch mapping. Please take a look, thanks

[gmao-rreichle]

@weiyuan-jiang, @zyj8881357 : I added this documentation based on a conversation with Yujin today. Perhaps the most important bit is the last sentence, which got me thinking that there may be a better way to deal with the EASEtile_to_Pfafcatch mapping. @weiyuan-jiang, how hard would it be to modify make_bcs so that in addition to the standard EASE tile space, we also create a supplemental *til file that is processed just like the cube-sphere tile file, except that the grid is EASE instead of cube-sphere. (That is, no restriction on the number of EASE tiles per EASE grid cell.)
It's possible that running the existing routines for an overlying ("atm") EASE grid requires only minimal modification. If the routines only need to know the center and edge coordinates of cube-sphere grid cells to do their work, we should be able to the same for the EASE grid now that EASE is in MAPL. We'd also have to deal with raster grid cells near the poles that are outside of the "global" EASE grid. But the rest of the tile generation shouldn't depend on the specifics of the overlying (atm) grid.
If this works, routing should be able use the supplemental "EASE" tile file to do the mapping.

[weiyuan-jiang]

I don't think it is hard to move this to make_bcs package

[gmao-rreichle]

Thanks, @weiyuan-jiang, for taking a look. Just to clarify, I didn't mean to "move" this program to make_bcs. Rather, I was wondering if we can use the existing make_bcs routines to generate a supplemental tile file (in nc4 format) with EASE tiles that respect the Pfaf catchment boundaries. Does your response apply to this latter option?

[weiyuan-jiang]

I am looking into it....

[zyj8881357]

The added documentation is correct and great.

[gmao-rreichle]

The "subarea" in the module and file name here refers to the portion of a Pfaf catchment that is within an EASE grid cell. In this sense, the "subarea" is identical to what would be a "tile" if we didn't limit the number of EASE tiles to at most one per EASE grid cell.
Within the script, the term "subtile" (or "sub-tile") is used to refer to the same thing. The term "subtile" is a bit confusing in this context because the Catchment model already has three "subtile" areas that vary dynamically ("wilting", "transpiring", "saturated").
My first instinct was to use "true tiles" instead of "subareas" because they correspond to what would be tiles in the cube-sphere or lat/lon space (i.e., tiles that respect the boudaries of Pfafstetter catchments). But "true tiles" isn't great because it's not very descriptive.
In any case, we should think about what we want to call these "subareas", and we need to make sure to use whatever term we pick consistently.

[zyj8881357]

How about to call them sub-catchments? Because everyone of them contributes a fraction of area and runoff to the catchment.

[weiyuan-jiang]

Now we have a EASE tile file which looks like a cubed-sphere tile file. Do you want to add elevation to it? @gmao-rreichle @zyj8881357

[zyj8881357]

Now we have a EASE tile file which looks like a cubed-sphere tile file. Do you want to add elevation to it? @gmao-rreichle @zyj8881357

Thank you! I have changed the routing model's code to read the new tile file instead of the tile2pfaf.nc4.

[gmao-rreichle]

we have a EASE tile file which looks like a cubed-sphere tile file. Do you want to add elevation to it?

Ideally, the new EASE tile file that looks like a cs tile file should have the same format and contents. So if elevation is in the cs tile file, it should also be in the new, cs-like EASE tile file, if only for consistency.