Skip to content

produce tile nc4 file while skipping land params generation #1495

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

Draft
weiyuan-jiang wants to merge 7 commits into
base: develop
Choose a base branch
from
Draft

produce tile nc4 file while skipping land params generation #1495

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
5f97329
fix typo
weiyuan-jiang Aug 6, 2026
9f06f6e
init doshmidt to avoid crash
weiyuan-jiang Aug 7, 2026
aa315c1
fix potential dateline bug
weiyuan-jiang Aug 7, 2026
19a3ab9
Merge branch 'hotfix/wjiang/bug_fixes' into feature/wjiang/tilenc4_fr…
weiyuan-jiang Aug 7, 2026
524baff
skipland with TileFile_ASCII_to_nc4
weiyuan-jiang Aug 10, 2026
dba0bc7
Merge branch 'develop' into feature/wjiang/tilenc4_from_til
weiyuan-jiang Aug 10, 2026
17846b8
move lake_type writing to MAPL
Aug 14, 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
5 changes: 2 additions & 3 deletions ...ridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CombineRasters.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ program mkOverlaySimple
logical :: DoZip
logical :: Verb
logical :: found
logical :: Merge, regrid
logical :: Merge
logical :: s_flag=.false.

character*4 :: tildir, rstdir
Expand Down Expand Up @@ -470,8 +470,7 @@ program mkOverlaySimple
if (.not. merge) then
if ( index(Grid1, "EASE") /=0 .and. index(Grid2, "EASE") /=0) then
l = index(TilFile, '.til')
regrid = nx /= i_raster .or. ny /= j_raster
call supplemental_tile_attributes(nx, ny, regrid, 'DE',TilFile(1:l-1), Rst1)
call supplemental_tile_attributes(nx, ny, TilFile(1:l-1), RstFile)
endif
endif

Expand Down
8 changes: 4 additions & 4 deletions ...mp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CubedSphere_GridMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ subroutine Get_CubedSphere_Grid(npx, npy, xlocs, ylocs, grid_type, shift_west, s
real(r8) :: alocs(2)
real(r8) :: target_lon,target_lat

logical doShiftWest
logical :: doShiftWest
logical :: do_schmidt

isg = 1
Expand Down Expand Up @@ -91,10 +91,10 @@ subroutine Get_CubedSphere_Grid(npx, npy, xlocs, ylocs, grid_type, shift_west, s
call mirror_grid(grid_global, 0, npts, npts, 2, 6)
! Cell Vertices
doShiftWest = .false.
do_schmidt = .false.
if (present(shift_west)) doShiftWest = shift_west
if (present(stg)) then
do_schmidt = .true.
end if
if (present(stg)) do_schmidt = .true.

!------------------------
! Schmidt transformation:
!------------------------
Expand Down
294 changes: 19 additions & 275 deletions .../GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/TileFile_ASCII_to_nc4.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,300 +5,44 @@
!
! wjiang, rreichle, 29 Nov 2024

program TileFile_ASCII_to_nc4
program TileRst_to_nc4
use, intrinsic :: iso_fortran_env, only: REAL64
use MAPL, only: MAPL_WriteTilingNC4, MAPL_ease_extent
use LogRectRasterizeMod, only: MAPL_UNDEF_R8

use rmTinyCatchParaMod
implicit none

character(512) :: arg
integer :: i, unit, unit2
integer :: nc, nr

character(:), allocatable :: tile_file
character(:), allocatable :: catchmentdef_file
real(REAL64), allocatable :: rTable(:,:)
integer, allocatable :: iTable(:,:)
character(128) :: gName1, gName2
character(len=512) :: tmpline

character(:), allocatable :: array(:)
character(len=:), allocatable :: filenameNC4

real :: cell_area

integer :: n_tile, n_grid, n_lon1, n_lat1, n_cat, tmp_in1, tmp_in2
integer :: n_lon2, n_lat2, nx, ny, num, ll, maxcat
character(:), allocatable :: rst_file
character(:), allocatable :: gName
character(:), allocatable :: filenameNC4
integer :: n
logical :: file_exists

!
! usage: TileRst_to_nc4.x nx ny gName
! ----------------------------------------------------------------------
!
! process command-line arguments

CALL get_command_argument(1, arg)
tile_file = trim(arg)
read(arg,'(i6)') nc
CALL get_command_argument(2, arg)
catchmentdef_file = trim(arg)

! ----------------------------------------------------------------------
!
! open and read *.til ASCII file

open (newunit=unit, file=trim(tile_file), form='formatted', action='read')

read (unit,"(A)") tmpline ! header line 1: N_tile [maxcat] nx ny (see below)
read (unit,*) N_grid ! header line 2: N_grid [=1 for EASE, =2 otherwise]
read (unit,*) gName1 ! header line 3: name of atm grid
read (unit,*) n_lon1 ! header line 4: N_lon of atm grid
read (unit,*) n_lat1 ! header line 5: N_lat of atm grid

! special treatment needed for header line 1 because maxcat is not included in legacy bcs

call split(tmpline, array, " ")
read(array(1), *) n_tile
num = size(array)
ll = 0
if (num == 4) then
ll = 1
read(array(2), *) maxcat ! number of Pfafstetter catchments
else
maxcat = -1 ! maxcat not available in legacy bcs
endif

read(array(2+ll), *) nx ! N_lon of raster grid
read(array(3+ll), *) ny ! N_lat of raster grid

if (N_grid == 1) then

! EASE grid tile space

! in some legacy bcs, dummy ocean grid info is included in header (despite N_grid=1);
! read next line and decide if it is dummy header or info for first tile

read (unit,"(A)") tmpline
if (index(tmpline,'OCEAN')/=0) then
read (unit,*)
read (unit,*)
read (unit,"(A)") tmpline
endif

else

! lat/lon or cube-sphere tile space

read (unit,*) gName2
read (unit,*) n_lon2
read (unit,*) n_lat2
read (unit,"(A)") tmpline ! read info for first tile (to accommodate legacy EASE grid issues above)

endif

allocate(iTable(N_tile,0:7))
allocate(rTable(N_tile,10))

rTable = MAPL_UNDEF_r8

! read ASCII tile file (NOTE: Info for first tile is already in tmpline!)

if ( index(gName1, 'EASE') /=0 ) then ! EASE grid tile space

read (tmpline,*) iTable(1,0), iTable(1,4), rTable(1,1), rTable(1,2), &
iTable(1,2), iTable(1,3), rTable(1,4)

do i = 2, N_tile
read (unit,*) iTable(i,0), iTable(i,4), rTable(i,1), rTable(i,2), &
iTable(i,2), iTable(i,3), rTable(i,4)
enddo

! rTable(:,4) is tile area fraction within grid cell (fr), convert to area;
! get fr back in WriteTilingNC4

call MAPL_ease_extent(gName1, tmp_in1, tmp_in2, cell_area=cell_area) ! get EASE grid cell area

rTable(:,3) = rTable(:,4)*cell_area
rTable(:,4) = cell_area

else ! lat/lon or cube-sphere tile space

read (tmpline,*) iTable(1,0), rTable(1,3), rTable(1,1), rTable(1,2), &
iTable(1,2), iTable(1,3), rTable(1,4), iTable(1,6), &
iTable(1,4), iTable(1,5), rTable(1,5), iTable(1,7)

do i = 2, N_tile
read (unit,*) iTable(i,0), rTable(i,3), rTable(i,1), rTable(i,2), &
iTable(i,2), iTable(i,3), rTable(i,4), iTable(i,6), &
iTable(i,4), iTable(i,5), rTable(i,5), iTable(i,7)
enddo
read(arg,'(i6)') nr
CALL get_command_argument(3, arg)
gName = trim(arg)

! re-define rTable(:,4) and rTable(:,5).
! fr will be re-created in WriteTilingNC4
tile_file = 'til/'//gName
rst_file = 'rst/'//gName//'.rst'

where (rTable(:,4) /=0.0)
rTable(:,4) = rTable(:,3)/rTable(:,4)
endwhere

where (rTable(:,5) /=0.0)
rTable(:,5) = rTable(:,3)/rTable(:,5)
endwhere
filenameNC4 = tile_file //'.nc4'

inquire(file=filenameNC4, exist=file_exists)
if (.not. file_exists) then
call supplemental_tile_attributes(nc,nr, tile_file, rst_file, write_catch=.false.)
endif

close(unit)

! ----------------------------------------------------------------------
!
! open and read catchment.def ASCII file

inquire( file= trim(catchmentdef_file), exist=file_exists)

if (file_exists) then

open (newunit=unit, file=trim(catchmentdef_file), form='formatted', action='read')

read(unit, *) n_cat ! number of *land* tiles

do i = 1, n_cat
read(unit, *) &
tmp_in1, &
tmp_in2, &
rTable(i, 6), &
rTable(i, 7), &
rTable(i, 8), &
rTable(i, 9), &
rTable(i,10)
enddo

close(unit)

endif

! assemble name of nc4 file

ll = index(tile_file, '.til')
filenameNC4 = tile_file(1:ll)//'nc4'

! write nc4 file

if (N_grid == 1) then
call MAPL_WriteTilingNC4(filenameNc4, [gName1 ], [n_lon1 ], [n_lat1 ], nx, ny, iTable, rTable, N_PfafCat=maxcat)
else
call MAPL_WriteTilingNC4(filenameNc4, [gName1, gName2], [n_lon1, n_lon2], [n_lat1, n_lat2], nx, ny, iTable, rTable, N_PfafCat=maxcat)
endif

contains

subroutine split(input_line,array,delimiters,order,nulls)

character(len=*),intent(in) :: input_line
character(len=*),optional,intent(in) :: delimiters
character(len=*),optional,intent(in) :: order
character(len=*),optional,intent(in) :: nulls
character(len=:),allocatable,intent(out) :: array(:)

integer :: n
integer,allocatable :: ibegin(:)
integer,allocatable :: iterm(:)
character(len=:),allocatable :: dlim
character(len=:),allocatable :: ordr
character(len=:),allocatable :: nlls
integer :: ii,iiii
integer :: icount
integer :: ilen
integer :: i10,i20,i30
integer :: icol
integer :: idlim
integer :: ifound
integer :: inotnull
integer :: ireturn
integer :: imax

! decide on value for optional DELIMITERS parameter
if (present(delimiters)) then ! optional delimiter list was present
if(delimiters/='')then ! if DELIMITERS was specified and not null use it
dlim=delimiters
else ! DELIMITERS was specified on call as empty string
dlim=' '//char(9)//char(10)//char(11)//char(12)//char(13)//char(0) ! use default delimiter when not specified
endif
else ! no delimiter value was specified
dlim=' '//char(9)//char(10)//char(11)//char(12)//char(13)//char(0) ! use default delimiter when not specified
endif
idlim=len(dlim) ! dlim a lot of blanks on some machines if dlim is a big string

if(present(order))then; ordr=adjustl(order); else; ordr='sequential'; endif ! decide on value for optional ORDER parameter

if(present(nulls))then; nlls=adjustl(nulls); else; nlls='ignore' ; endif ! optional parameter

n=len(input_line)+1 ! max number of strings INPUT_LINE could split into if all delimiter
allocate(ibegin(n)) ! allocate enough space to hold starting location of tokens if string all tokens
allocate(iterm(n)) ! allocate enough space to hold ending location of tokens if string all tokens
ibegin(:)=1
iterm(:)=1

ilen=len(input_line) ! ILEN is the column position of the last non-blank character
icount=0 ! how many tokens found
inotnull=0 ! how many tokens found not composed of delimiters
imax=0 ! length of longest token found

select case (ilen)

case (0) ! command was totally blank

case default ! there is at least one non-delimiter in INPUT_LINE if get here
icol=1 ! initialize pointer into input line
INFINITE: do i30=1,ilen,1 ! store into each array element
ibegin(i30)=icol ! assume start new token on the character
if(index(dlim(1:idlim),input_line(icol:icol))==0)then ! if current character is not a delimiter
iterm(i30)=ilen ! initially assume no more tokens
do i10=1,idlim ! search for next delimiter
ifound=index(input_line(ibegin(i30):ilen),dlim(i10:i10))
IF(ifound>0)then
iterm(i30)=min(iterm(i30),ifound+ibegin(i30)-2)
endif
enddo
icol=iterm(i30)+2 ! next place to look as found end of this token
inotnull=inotnull+1 ! increment count of number of tokens not composed of delimiters
else ! character is a delimiter for a null string
iterm(i30)=icol-1 ! record assumed end of string. Will be less than beginning
icol=icol+1 ! advance pointer into input string
endif
imax=max(imax,iterm(i30)-ibegin(i30)+1)
icount=i30 ! increment count of number of tokens found
if(icol>ilen)then ! no text left
exit INFINITE
endif
enddo INFINITE

end select

select case (trim(adjustl(nlls)))
case ('ignore','','ignoreend')
ireturn=inotnull
case default
ireturn=icount
end select
allocate(character(len=imax) :: array(ireturn)) ! allocate the array to return
!allocate(array(ireturn)) ! allocate the array to turn

select case (trim(adjustl(ordr))) ! decide which order to store tokens
case ('reverse','right') ; ii=ireturn ; iiii=-1 ! last to first
case default ; ii=1 ; iiii=1 ! first to last
end select

do i20=1,icount ! fill the array with the tokens that were found
if(iterm(i20)<ibegin(i20))then
select case (trim(adjustl(nlls)))
case ('ignore','','ignoreend')
case default
array(ii)=' '
ii=ii+iiii
end select
else
array(ii)=input_line(ibegin(i20):iterm(i20))
ii=ii+iiii
endif
enddo
end subroutine split

end program

! ======================= EOF ====================================================
9 changes: 5 additions & 4 deletions ...p/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/TileFile_ASCII_to_nc4.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@

import os
import glob
import subprocess
import subprocess as sp

if __name__ == "__main__":

Expand All @@ -32,7 +32,8 @@
catch_file = til.replace('geometry', 'land')
fname = os.path.basename(til)
catch_file = catch_file.replace(fname, 'clsm/catchment.def')
result = sp.run["./TileFile_ASCII_to_nc4.x", til, catch_file, capture_output=True, text=True)
print(result.stdout)


#result = sp.run(["./TileFile_ASCII_to_nc4.x", til, catch_file], capture_output=True, text=True)
#print(result.stdout)
print("Need to revisit this program")
os.chdir(PWD)
Loading