Skip to content

fix ease index in ease pfaf tile file and remove pole tiles #1430

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

Open
weiyuan-jiang wants to merge 5 commits into
base: develop
Choose a base branch
from
+39 −6
Open

fix ease index in ease pfaf tile file and remove pole tiles #1430

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
ff1ada1
fix EASE index for EASE pfaf nc4 file
weiyuan-jiang May 21, 2026
9f47996
change the indexes in ocean too
weiyuan-jiang May 21, 2026
83c0226
fix codes to use the new indexing
weiyuan-jiang May 21, 2026
f50a92e
remove pole tiles
weiyuan-jiang May 22, 2026
51ccb06
Adjust yi_tile calculation for sub-catchment indexing
weiyuan-jiang May 22, 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
2 changes: 1 addition & 1 deletion ...dComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -692,7 +692,7 @@ subroutine create_mapping_handler(tilegrid, pfaf_tilegrid, rc)
kk=0
do ii=1,np
if(type_tile(ii)==100)then
tile_id=map_tile(xi_tile(ii)+1,ny-yi_tile(ii)) ! note: EASE indices are 0-based
tile_id=map_tile(xi_tile(ii)+1, yi_tile(ii)+1) ! note: EASE indices are 0-based
Copy link
Copy Markdown
Contributor

@gmao-rreichle gmao-rreichle May 21, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@zyj8881357 : Is there a place elsewhere in the code (e.g., in preprocessing) where the orientation of the latitude indices (j_indg) from the EASE-Pfaf file needs to be flipped?

Copy link
Copy Markdown

@zyj8881357 zyj8881357 May 22, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. for the file GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_num_sub_catchment.f90.

The line 54: yi_tile = yi_tile + 1
should be changed to
yi_tile = nlatE - yi_tile
Also, add nlatE in line 6 as:
use routing_model_constants, only : nc, nlon, nlat, np=>np_tot, nlatE

That is all.

if(1<=tile_id.and.tile_id<=route%nt_global)then
kk=kk+1
tid_patch( kk)=tile_id ! the EASE id for the patch (src)
Expand Down
39 changes: 36 additions & 3 deletions ...omp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -968,12 +968,17 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id)
real*4, allocatable, target :: q0 (:,:)
real(REAL64), allocatable :: rTable(:,:)
integer, allocatable :: iTable(:,:)
real(REAL64), allocatable :: rTable_keep(:,:)
integer, allocatable :: iTable_keep(:,:)
character(len=128) :: gName(2)
logical, allocatable :: IsOcean(:)
logical, allocatable :: keep_tile(:)
! This is used to adjust EASE grid from 1-based to 0-based indexes
! The tile file with only one EASE grid is already 0-based and may not go through this subroutine
! This is a special case for river-routing. The ocean grid is also EASE just for convenience
logical :: two_EASE
integer :: ip_keep, k
real(REAL64), parameter :: EASE_LAT_MAX = 85.04459_REAL64
Copy link
Copy Markdown
Contributor

@gmao-rreichle gmao-rreichle May 22, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@weiyuan-jiang : The min/max latitudes for the EASE grid depend on the resolution. For consistency, we should probably get this from MAPL_ease_extent()


! -----------------------------------------------------
!
Expand Down Expand Up @@ -1222,10 +1227,12 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id)

if (two_EASE) then
iTable(:,2) = iTable(:,2) - 1
iTable(:,3) = iTable(:,3) - 1
!flip to make it consistent with conventional EASE indexing
iTable(:,3) = JM(1) - iTable(:,3)

where (iTable(:,0) == 0)
iTable(:,4) = iTable(:,4) -1
iTable(:,5) = iTable(:,5) -1
iTable(:,5) = JM(1) - iTable(:,5)
endwhere
j = index(gName(2), "-Pfafstetter")
gName(2) = gName(2)(1:j-1)
Expand All @@ -1247,7 +1254,33 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id)
endwhere

fname=trim(fnameTil)//'.nc4'
call MAPL_WriteTilingNC4(fname, gName(1:n_grid), im(1:n_grid), jm(1:n_grid), nx, ny, iTable, rTable, N_PfafCat=SRTM_maxcat, rc=status)
if (two_EASE) then
! remove tiles outside EASE grid domain
allocate(keep_tile(ip))
keep_tile = (rTable(1:ip,2) <= EASE_LAT_MAX) .and. (rTable(1:ip,2) >= -EASE_LAT_MAX)
ip_keep = count(keep_tile)
ASSERT_(ip_keep > 0)

allocate(iTable_keep(ip_keep,0:7))
allocate(rTable_keep(ip_keep,10))

k = 0
do n = 1, ip
if (keep_tile(n)) then
k = k + 1
iTable_keep(k,0:7) = iTable(n,0:7)
rTable_keep(k,1:10) = rTable(n,1:10)
endif
enddo

call MAPL_WriteTilingNC4(fname, gName(1:n_grid), im(1:n_grid), jm(1:n_grid), &
nx, ny, iTable_keep, rTable_keep, N_PfafCat=SRTM_maxcat, rc=status)

deallocate(iTable_keep, rTable_keep, keep_tile)
else
call MAPL_WriteTilingNC4(fname, gName(1:n_grid), im(1:n_grid), jm(1:n_grid), &
nx, ny, iTable, rTable, N_PfafCat=SRTM_maxcat, rc=status)
endif

deallocate (rTable, iTable)
deallocate (limits)
Expand Down
4 changes: 2 additions & 2 deletions ...ridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_num_sub_catchment.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ program main

use MAPL
use river_ncfile_helper
use routing_model_constants, only : nc, nlon, nlat, np=>np_tot
use routing_model_constants, only : nc, nlon, nlat, np=>np_tot, nlatE

implicit none

Expand Down Expand Up @@ -51,7 +51,7 @@ program main
call read_ncfile_int1d(trim(file_path1),"i_indg",xi_tile,np)
xi_tile = xi_tile + 1
call read_ncfile_int1d(trim(file_path1),"j_indg",yi_tile,np)
yi_tile = yi_tile + 1
yi_tile = nlatE - yi_tile
Copy link
Copy Markdown
Contributor

@gmao-rreichle gmao-rreichle May 22, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@zyj8881357 : I tried to understand why the lat index is flipped here. It looks like yi_tile is needed to get the climatological runoff from the file SMAPL4_NRv11p4_M09_runoff_mean_2001_2023.nc, which is read in get_Qr_clmt.f90
It looks like the lat index orientation in the nc file above is consistent with the EASE convention.
But after reading the data, the lat index is flipped here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_Qr_clmt.f90

Line 32 in ff111ad

temp = runoff(:, nlat:1:-1)

It looks like if we don't flip the lat index at all, the processing should work out. But I didn't inspect every little bit of the code and didn't test it. There may be other processing involved that requires the lat index.
If my hunch is right, though, we should avoid flipping the lat index twice.

call read_ncfile_int1d(trim(file_path1),"pfaf_index",catid_tile,np)
call read_ncfile_real1d(trim(file_path1),"area",area_tile,np)
area_tile = area_tile * MAPL_RADIUS**2 / 1.e6 !km^2
Expand Down