Add river routing module for offline (GEOSldas) simulations#1143
Add river routing module for offline (GEOSldas) simulations#1143weiyuan-jiang wants to merge 125 commits intodevelopfrom
Conversation
…t in GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp
…ent (GEOS_RouteGridComp.F90)
That topo package has no dependence on MAPL. It may be used by others without MAPL. |
...m_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_pfaf_frac.F90
Outdated
Show resolved
Hide resolved
@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. |
gmao-rreichle
left a comment
gmao-rreichle
left a comment
There was a problem hiding this comment.
@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.
...GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90
Outdated
Show resolved
Hide resolved
GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/reservoir.F90
Outdated
Show resolved
Hide resolved
GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/reservoir.F90
Outdated
Show resolved
Hide resolved
GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/reservoir.F90
Outdated
Show resolved
Hide resolved
...GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90
Outdated
Show resolved
Hide resolved
...ysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/create_river_input.py
Outdated
Show resolved
Hide resolved
gmao-rreichle
left a comment
gmao-rreichle
left a comment
There was a problem hiding this comment.
@zyj8881357 : I'm still struggling to understand what exactly is done in EASE_pfaf_frac.F90.
...m_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_pfaf_frac.F90
Outdated
Show resolved
Hide resolved
…also change variable pfaf_frac to subtile_area in the *_tile2pfaf.nc4 file
|
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. |
…S_RouteGridComp.F90, create_river_input.py)
…UT'; added comments (GEOS_RouteGridComp.F90)
…; moved mk_RouteRestarts.F90 to dir of obsolete files (routing_model.F90, CMakeLists.txt, mk_RouteRestarts.F90)
gmao-rreichle
left a comment
gmao-rreichle
left a comment
There was a problem hiding this comment.
@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.
| v = Variable(type=PFIO_INT32, dimensions='tile') | ||
| call v%add_attribute('units', '1') | ||
| call v%add_attribute('long_name', 'tile_id') | ||
| call meta%add_variable('tile_id', v) | ||
|
|
||
| v = Variable(type=pFio_INT32, dimensions='tile') | ||
| call v%add_attribute('units', '1') | ||
| call v%add_attribute('long_name', 'pfaf_index') | ||
| call meta%add_variable('pfaf_index', v) | ||
|
|
||
| v = Variable(type=pFio_REAL32, dimensions='tile') | ||
| call v%add_attribute('units', 'm+2') | ||
| call v%add_attribute('long_name', 'area_of_subtile') | ||
| call meta%add_variable('subtile_area', v) |
There was a problem hiding this comment.
@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.
|
|
||
| ! ---------------------------------------------------- | ||
|
|
||
| subroutine create_mapping_handler(tilegrid, pfaf_tilegrid, rc) |
There was a problem hiding this comment.
@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?
| /bin/mv {GRIDNAME}.j geometry/{GRIDNAME}/. | ||
| /bin/cp til/{GRIDNAME}{RS}.til geometry/{GRIDNAME}/. | ||
| /bin/cp til/{GRIDNAME}{RS}.nc4 geometry/{GRIDNAME}/. | ||
| /bin/cp til/{GRIDNAME}_tile_pfaf.nc4 geometry/{GRIDNAME}/. |
There was a problem hiding this comment.
Should this read
/bin/cp til/{GRIDNAME}_tile2pfaf.nc4 geometry/{GRIDNAME}/.
(i.e., tile2pfaf instead of tile_pfaf)
| lati09_output_file = "temp/lati_1m_M09.txt" | ||
| loni09_output_file = "temp/loni_1m_M09.txt" |
There was a problem hiding this comment.
This py script is specific to the M09 EASE grid. I think its product is only needed by "get_Qr_clmt.f90", and the M09 dependency stems from the fact that "get_Qr_clmt.f90" creates a runoff climatology from a SMAP L4 Open Loop simulation.
Is the above correct?
If so, we need to document this dependency in this file and also in readme.txt
| open(10, file="temp/lati_1m_M09.txt") | ||
| read(10, *) lati | ||
| open(11, file="temp/loni_1m_M09.txt") | ||
| read(11, *) loni |
There was a problem hiding this comment.
This program is specific to the M09 EASE grid. I think what's produced by this program is only needed by "get_Qr_clmt.f90", and the M09 dependency stems from the fact that "get_Qr_clmt.f90" creates a runoff climatology from a SMAP L4 Open Loop simulation.
Is the above correct?
If so, we need to document this dependency in this file and also in readme.txt
| WRIVER = WRIVER + QINFLOW_LOCAL*real(route_dt) | ||
|
|
||
| if (associated(QSFLOW)) QSFLOW = QSFLOW_OUT ! outflow from local streams | ||
| if (associated(QOUTFLOW)) QOUTFLOW = QOUTFLOW_OUT ! outflow from main river |
There was a problem hiding this comment.
@zyj8881357 : I'm a bit confused about the meaning of QOUT_CAT vs. QOUTFLOW_OUT. I think that QOUT_CAT is the outflow from the Catchment, which for active reservoirs is the same as the outflow from the reservoir, whereas QOUTFLOW_OUT is the discharge from the main river before taking into account any reservoir.
If that's indeed the case, it's not a problem as such.
But I'm wondering how users would interpret our export variables (QSFLOW, QOUTFLOW, QRES). I don't think we provide the user with sufficiently clear information in the LONG_NAME attributes to assist them with the correct interpretation. Specifically, I think we're asking the user to apply the "active reservoir flag" on their end if what they want it the outflow from the catchment after the reservoir. I don't think that should be left to the user. It's probably better to also add QOUT_CAT as an export variable.
| call MAPL_GetPointer(EXPORT, QSFLOW, 'QSFLOW', _RC) | ||
| call MAPL_GetPointer(EXPORT, QINFLOW, 'QINFLOW' , _RC) | ||
| call MAPL_GetPointer(EXPORT, QOUTFLOW, 'QOUTFLOW', _RC) | ||
| call MAPL_GetPointer(EXPORT, QRES, 'QRES', _RC) |
There was a problem hiding this comment.
@weiyuan-jiang, @zyj8881357, I think we should add ALLOC=.true. here, as in:
call MAPL_GetPointer(EXPORT, QSFLOW, 'QSFLOW', ALLOC=.true., _RC)
With this change, these export variables are always allocated and can be used regardless of whether the user writes them out via HISTORY.
Then we wouldn't need the corresponding local variables *_IN and *_OUT. It's a bit confusing to essentially have two copies of each of these.
I'm not 100% sure I'm getting this right. @weiyuan-jiang, please advise.
| @@ -0,0 +1,193 @@ | |||
| module river_read | |||
There was a problem hiding this comment.
This file isn't listed in readme.txt. It looks like it's just a wrapper for netcdf calls, and as such doesn't really require documentation in readme.txt. But the file name river_read.f90 is misleading. I suggest renaming this to something like read_ncfile_helper.F90. It doesn't look like there's anything specific to the river routing model in this helper.
| @@ -0,0 +1,1029 @@ | |||
| module k_module | |||
There was a problem hiding this comment.
I suggest adding the file name k_module_cali.f90 into readme.txt, under what is currently item 6 ("get_K_model_calik.f90"). Currently, this file isn't mentioned in readme.txt
| integer,parameter :: nlath = 33600 | ||
|
|
||
| integer,parameter :: nl_USGS = 3352492 ! Total number of USGS data records | ||
| integer,parameter :: nt09 = 1684725 ! Total number of catchment gridcell in M09 |
There was a problem hiding this comment.
Same as "nt_global09" above. We shouldn't need this extra parameter
…ld error from earlier commit (routing_model.F90)
5 participants
Add river routing module for offline land modeling (GEOSldas). Replaces #1023
Left for future PRs:
Related PRs:
Testing:
List of tasks that remain to be addressed (in no particular order):
@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.