Remove cs indg from tile_coord structure

[weiyuan-jiang]

Address issue #111

[gmao-rreichle]

@weiyuan-jiang : I don't think this is quite right yet:

1. In line 604 of GEOS_LdasGridComp.F90, the call to get_ij_ind_from_latlon() appears to reset %i_indg and %j_indg. These fields should NOT be touched after they are initialized. Shouldn't the final two arguments in this call be the new simple 1d arrays?

2. The name of the simple 1d arrays is still confusing because the underlying grid is a lat/lon grid, not a cs grid. I suggest renaming the simple arrays as follows:
   cs_i_indg --> latlonhash4cs_i_indg
   cs_j_indg --> latlonhash4cs_j_indg

3. There might still be some confusion about when "[latlonhash4]cs_i_indg" and "%i_indg" is used (likewise for "j").

Perhaps take another look at the write-up in #111 .

Also, we need to wait for @saraqzhang to get back and make sure we are testing her assimilation case on the cs grid.

[weiyuan-jiang]

Sara has tested it before I pulled request.

1. Indeed, I reset %i_indg and %j_indg. But that is a meaningful reset. It belongs to a lat/lon grid's i_indg and j_indg
2. the cs_i_indg and cs_j_indg is relative to cs grid not lat/lon grid. i_indg and j_indg is related to lat/lon or EASE grid.

[weiyuan-jiang]

Cut and past from Rolf's doc:

Reichle, 10 July 2020

Somewhere early in the code (or in preprocessing?), tile_coord_f is initialized using information from the *.til file (LocationStream) and/or an existing tile_coord file that is consistent with the LocationStream.

At this point, the fields %i_indg and %j_indg refer to the grid associated with the tile space. For an EASE-based tile space, this is an EASE grid. For a LatLon-based tile space, this is a LatLon grid. For a cs-based tile space, this is a cs grid.

%i_indg and %j_indg should not be re-set or re-defined after initialization. Otherwise, there is potential for misinterpretation later on and s/w maintenance and development becomes unnecessarily complicated.

The modifications outlined below suggest simple 1d arrays hash_[x]_indg[_f].
Alternatively, the same info could also be included as two new fields in the tile_coord[_f] structure under the names %hash_i_indg and %hash_j_indg.
It’s not obvious which alternative is ultimately easier.
If we go with simple arrays, we probably need to change the interface to several of the helper subroutines in TileCoordRoutines.
If we go with additional tile_coord fields, we would need to edit the guts of these helper routines.

Perhaps it would be best to first create a new branch that removes obsolete subroutines and “use” statements (search for “obsolete” below), without making any of the other changes. This should be very straightforward, easy to test and helpful in any case. Once we have cleaned up, the rest of the modifications may be a little bit easier.

---

src/Components/GEOSldas_GridComp/GEOSmetforce_GridComp/LDAS_Forcing.F90:

1. Define simple 1d arrays:

   public :: hash_i_indg, hash_j_indg ! initialized in GEOS_LdasGridComp

(that is, same as in git branch “remove_cs_indg”, but change the name because these arrays will be used for all tile spaces, not just cs-based tile spaces)

---

src/Components/GEOSldas_GridComp/GEOS_LdasGridComp.F90

1. re-name cs_i_indg_f to hash_i_indg_f
   re-name cs_j_indg_f to hash_j_indg_f

2. Define hash_i_indg_f and tile_grid_f for LatLon and EASE tile spaces (~line 584)

   hash_i_indg_f = tile_coord_f%i_indg
   hash_j_indg_f = tile_coord_f%j_indg

3. (Re-)Define hash_i_indg_f and tile_grid_f for cs tile space:

~595-596: remove the two lines that “save [the] original index”

~line 604, set hash_i_indg_f and hash_j_indg_f via the call to get_ij_ind_from_latlon:

call get_ij_ind_from_latlon(latlon_tmp_g,tile_coord_f(i)%com_lat,tile_coord_f(i)%com_lon, &
             hash_i_indg_f, hash_j_indg_f )	                          

~line 607, the call to get_tile_grid() in case of the c3 grid can presumably be replaced with

tile_grid_f = tile_grid_g      

(It looks like we always use tile_grid_f with global extent for cs tile spaces. In future, this will require more development, otherwise small cs domains may run too slowly.)

4. ~line 635, MPI_BCAST hash_i_indg_f and hash_j_indg_f for all tile spaces (incl LatLon, EASE and cs). That is, always execute this block (for the renamed hash_[x]_indg), not just the cs tile space.

---

At this point, we should have initialized the “[x]_indg” variables as follows:

tile_coord[_f?]%[x]_indg : original meaning relative to grid associated with tile space
relative to LatLon grid for LatLon-based tile space
relative to EASE grid for EASE-based tile space
relative to cs grid for cs-based tile space

hash_[x]_indg[_f] : index relative to “hash” grid:
identical to tile_coord%[x]_indg for EASE & LatLon tile space
relative to LatLon grid for cs tile space

---

The BIG challenge is to identify all instances where “tile_coord[_f]%[x]indg” is used but “hash[x]indg[f]” should now be used, and to change the call statements accordingly. This probably requires changing the interfaces of some of the helper subroutines.

---

src/Components/GEOSldas_GridComp/GEOSlandpert_GridComp/GEOS_LandPertGridComp.F90

~line 1063-1064
internal%i_indgs(:)=tile_coord(:)%i_indg
internal%j_indgs(:)=tile_coord(:)%j_indg

I’m guessing this should change to:

internal%i_indgs(:)=hash_i_indg
internal%j_indgs(:)=hash_j_indg

Also, the four calls to “tile2grid_simple()” may need to use “hash_[x]_indg” rather than “%[x]_indg”.

Likewise for the eight calls to “grid2tile()”.

Finally, the eight calls to “tile_mask_grid()” with arguments “internal%[x]_indgs” may need change.

But I’m not sure at all. My hesitation is that I don’t know on what grid the perturbations are computed. For EASE and LatLon, it should be the grid associated with the tile space, but it could also be a “coarsened” version of that grid. (See field %coarsen in “pert_param” variables.)
For the cs tile space, the perturbations grid should be some LatLon grid, but I’m not sure where and how exactly this is defined.
Also, I’m not sure what [x]_indgs is (as opposed to [x]_indg).

---

src/Components/GEOSldas_GridComp/Shared/LDAS_TileCoordRoutines.F90

subroutine LDAS_create_grid_g -- not relevant

subroutine LDAS_read_land_tile -- not relevant (preprocessing only)
subroutine read_catchment_def -- not relevant (preprocessing only)
subroutine read_grid_elev -- not relevant (preprocessing only)
subroutine fix_dateline_bug_in_tilecoord -- not relevant (preprocessing only)

subroutine read_til_file -- not relevant (obsolete, remove?)
subroutine get_land_tile_info -- not relevant (obsolete, remove?)
subroutine extract_land_tiles -- not relevant (obsolete, remove?)
subroutine reorder_tiles -- not relevant (obsolete, remove?)
subroutine tile2grid tiles -- not relevant (obsolete, remove?)
subroutine grid2tile_old -- not relevant (obsolete, remove?)
subroutine get_is_land -- not relevant (obsolete, remove?)

subroutine get_tile_grid -- bypass for cs tiles in LdasGridComp init as suggested above

I think the following subroutines would need to use “hash_[x]_indg” instead of “tile_coord%[x]_indg”:

subroutine get_number_of_tiles_in_cell_ij
subroutine get_tile_num_in_cell_ij

Called for mapping from tiles to obs space:

called from get_obs_pred()        in clsm_ensupd_upd_routines.F90 
called from get_enkf_increments() in clsm_ensupd_enkf_update.F90  

subroutine get_tile_num_in_ellipse

called from get_obs_pred() in clsm_ensupd_upd_routines.F90

subroutine get_tile_num_from_latlon

Called to find a tile to which the obs is assigned:

called from get_tile_num_for_obs() in clsm_ensupd_read_obs.F90

The following helper subroutine already takes simple 1d arrays for [x]_indg:

subroutine get_ij_ind_from_latlon

called from get_tile_num_from_latlon() in LDAS_TileCoordRoutines.F90
called from get_tile_num_in_ellipse()  in LDAS_TileCoordRoutines.F90
called from GEOS_LdasGridComp.F90      (see above)
called from get_obs_pert()             in clsm_ensupd_upd_routines.F90 (three times!)

--> probably need to pass through “hash_[x]_indg” via the calling subroutines

Note obsolete “use” statements for get_ij_ind_from_latlon in:
GEOS_LandPertGridComp.F90
LDAS_PertRoutines.F90
Remove?

The following subroutines are used only in GEOS_LandPertGridComp.F90 (see above):

subroutine tile2grid_simple
subroutine tile_mask_grid

The following subroutines

subroutine grid2tile_new
subroutine grid2tile_real8

called from write_smapL4SMaup() in clsm_ensupd_enkf_update.F90,
should use %[x]indg (not “hash[x]_indg”)

see above for the calls from GEOS_LandPertGridComp.F90

[tclune]

CMake changes ok.

[weiyuan-jiang]

closed and replaced by the PR #295