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/mk_restart? (@weiyuan-jiang, @zyj8881357)

---

@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 added a few commits and in-line comments (see below).

Here's a link to the aggregated changes from my commits (with white-space changes hidden):
https://github.com/GEOS-ESM/GEOSgcm_GridComp/compare/49790fa..f452878?w=1
The changes are focused on cleaning up the code and improving the documentation/comments (including units).
Please review my changes and let me know if anything looks wrong.

The in-line comments below include a couple of minor changes but also bigger things like changing which parameters should be in the bcs and routing restart files. It's probably best to address the changes that are zero-diff first, and verify that the revised branch produces zero-diff results. Once that is done, we can tackle the change in the routing (and reservoir) parameters, and verify that the resulting simulations produce scientifically identical results (should be within "roundoff")

Please let me know if you have any questions

cc: @weiyuan-jiang

[zyj8881357]

@gmao-rreichle Thank you for all the comments and changes. I will review them and edit the branch accordingly.

[gmao-rreichle]

@zyj8881357 : Thanks for making the changes. I "resolved" the comments that were fully addressed. I also added another commit 64775ef. Please double-check my changes. You may want to "hide white-space changes" when viewing.

---

@weiyuan-jiang, @biljanaorescanin:

While I was going through the routing PR I noticed that in the atm topo package we hardwire MAPL_Radius here (4 times):

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90

Line 750 in 64775ef

great_circle_dist(p1,p2,6371.d0), &

I assume it would be 0-diff if we replaced the hardcoded value with MAPL_Radius via
"use MAPL_ConstantsMod, only: MAPL_radius"
I'd be ok with making this change on the routing PR.

Same goes for using MAPL_PI here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90

Line 55 in 64775ef

real(REAL64), parameter :: PI = 3.14159265358979323846

[weiyuan-jiang]

@zyj8881357 : Thanks for making the changes. I "resolved" the comments that were fully addressed. I also added another commit 64775ef. Please double-check my changes. You may want to "hide white-space changes" when viewing.

@weiyuan-jiang, @biljanaorescanin:

While I was going through the routing PR I noticed that in the atm topo package we hardwire MAPL_Radius here (4 times):

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90

Line 750 in 64775ef

great_circle_dist(p1,p2,6371.d0), &

I assume it would be 0-diff if we replaced the hardcoded value with MAPL_Radius via
"use MAPL_ConstantsMod, only: MAPL_radius"
I'd be ok with making this change on the routing PR.
Same goes for using MAPL_PI here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90

Line 55 in 64775ef

real(REAL64), parameter :: PI = 3.14159265358979323846

That topo package has no dependence on MAPL. It may be used by others without MAPL.

[weiyuan-jiang]

I have not turned on the switch here. I checked the area calculated and the values in the file. About 4 significant digits are the same. Please very it and flip the switch. @zyj8881357

[gmao-rreichle]

@weiyuan-jiang, many thanks for helping address this. I'm confident your solution works -- and I'm sure @zyj8881357 can verify. But the solution seems more complicated than I had anticipated. To recap, what is needed here is the area of each grid cell on the "raster" grid, which is a regular 21,600-by-10,800 lat/lon grid (i.e., 1 arc-minute). The grid cell area should thus be area=R*R*delta_lat*delta_lon*cos(lat), where area is in km2, R=6371 km (MAPL_Radius) and delta_lat, delta_lon, and lat are in radians. Note that n_lat and n_lon are defined (hard-coded) at the top of the file.

This simple calculation has several benefits:

• The area depends only on lat. That is, it can be computed outside of the loop through longitudes.
• We can probably avoid allocating a 21,600-by-10,800 array for cellarea.
• We can avoid the dependency on subroutines in the makebcs/topo package. If any dependencies are added, they should be on MAPL or ESMF.

[weiyuan-jiang]

Ahh, your formula is not just a better one, it is a correct one. The formula I used is based on big cricle, your is based on lat-lon. I will correct it.

[gmao-rreichle]

@weiyuan-jiang, @zyj8881357 : I spoke with Weiyuan about further simplifying the calculation of the area of each raster grid cell, and Weiyuan will address it. Since speaking with Weiyuan, I noticed that further changes and cleanup are necessary in the preproc routines:

Utils/Raster/preproc/routing_model/get_num_sub_catchment.f90
In this routine, we call subroutine EASE_Find_subs(), and after Weiyuan's changes the cellarea will be calculated inside this subroutine. That is, code related to cellarea in get_num_sub_catchment.f90 needs to be cleaned up.

Utils/Raster/preproc/routing_model/get_area.f90
Do we need this program? It appears to be specific to the M09 grid, and the preprocessing should work for other EASE grid resolutions as well. Could this be a leftover?

Utils/Raster/preproc/routing_model/get_latloni_cellarea.py
What does this py script do? If all it does is calculate the area of each grid cell on the 1 arc minute raster grid and then write the. cellarea.nc file, it can be removed.

Utils/Raster/preproc/routing_model/run_routing_preproc.py
Remove call to get_latloni_cellarea.py

I might be wrong, and there might be other changes and cleanup needed.

[zyj8881357]

All have been cleaned up.

[gmao-rreichle]

@zyj8881357 : Thanks for making the changes. I "resolved" the comments that were fully addressed. I also added another commit 64775ef. Please double-check my changes. You may want to "hide white-space changes" when viewing.
@weiyuan-jiang, @biljanaorescanin:
While I was going through the routing PR I noticed that in the atm topo package we hardwire MAPL_Radius here (4 times):

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90

Line 750 in 64775ef

great_circle_dist(p1,p2,6371.d0), &

I assume it would be 0-diff if we replaced the hardcoded value with MAPL_Radius via
"use MAPL_ConstantsMod, only: MAPL_radius"
I'd be ok with making this change on the routing PR.
Same goes for using MAPL_PI here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90

Line 55 in 64775ef

real(REAL64), parameter :: PI = 3.14159265358979323846

That topo package has no dependence on MAPL. It may be used by others without MAPL.

@mathomp4, @weiyuan-jiang : I didn't realize that we are avoiding a dependency on MAPL in the new makebcs topo package, which is fine. But that shouldn't stop us from using constants values that are consistent with those in MAPL when we run the package for GEOS. I added a commit to make this a bit clearer: 6efc70f (please double-check my work; PS: I fixed a silly syntax error in the subsequent commit). Fortunately, the value for the Earth radius is already consistent with MAPL_RADIUS. But there are 5 places where the value of PI is hard-coded, and it's hard-coded to two distinct values, which is ugly. Do we really want to keep the distinct values of PI? I realize changing this would not be zero-diff, although it may still be possible to make the change before the soon-to-be-released v12 GCM (which I think is meant to use newly generated topography files) goes any further. At least I think we should make a conscious decision either way.

cc: @biljanaorescanin

[gmao-rreichle]

@zyj8881357 : Thanks for implementing the change to using the alpha parameters. There is probably a good amount of cleanup to do, see comments below. The comments aren't necessarily complete or 100% correct, I didn't go through the code in great detail. Maybe I missed something.

[gmao-rreichle]

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

[gmao-rreichle]

@zyj8881357 : I'm still struggling to understand what exactly is done in EASE_pfaf_frac.F90. In the above lines, is the objective to find the EASE grid cell i and j indices for a given raster grid cell? If so, the simpler way to do this would be to use MAPL_ease_convert(). This subroutine will give you exactly that information. Note the somewhat odd ordering of that subroutine's argument list.