Skip to content
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

update bin_to_cube tool #6294

Merged
merged 4 commits into from
May 31, 2024
Merged

update bin_to_cube tool #6294

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
af6cbb5
update bin_to_cube according to https://github.com/NCAR/Topo.git
jsbamboo Mar 13, 2024
80817ee
make input parameters to command arguments and clean up signatures
jsbamboo Mar 22, 2024
064de81
remove landfrac in bin_to_cube and cube_to_target; correct some print…
jsbamboo Mar 22, 2024
90d4e2d
remove comment in bin_to_cube
whannah1 Apr 15, 2024
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
88 changes: 70 additions & 18 deletions components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ program convterr
!
integer :: im, jm

integer, parameter :: ncube = 3000 !dimension of cubed-sphere grid
integer :: ncube !dimension of cubed-sphere grid
! integer, parameter :: ncube = 540 !dimension of cubed-sphere grid
! integer, parameter :: ncube = 361 ! for debugging

Expand Down Expand Up @@ -72,11 +72,25 @@ program convterr
real(r8) :: vol,dx_rad,vol_cube,area_latlon,darea_latlon ! latitude array
real(r8), allocatable, dimension(:,:) :: darea_cube

INTEGER :: UNIT

character(len=1024) :: raw_latlon_data_file,output_file

namelist /binparams/ &
raw_latlon_data_file,output_file,ncube

UNIT=221
OPEN( UNIT=UNIT, FILE="bin_to_cube.nl" ) !, NML = cntrls )
READ( UNIT=UNIT, NML=binparams)
CLOSE(UNIT=UNIT)
write(*,*) "Intermediate cubed-sphere resolution ",ncube

!
! read in USGS data from netCDF file
!
! status = nf_open('topo-lowres.nc', 0, ncid) !for debugging
status = nf_open('usgs-rawdata.nc', 0, ncid)
status = nf_open(raw_latlon_data_file, 0, ncid)
write(*,*) "Opening: ",TRIM(raw_latlon_data_file)
IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR(STATUS)

status = NF_INQ_DIMID(ncid, 'lat', dimlatid)
Expand Down Expand Up @@ -105,13 +119,13 @@ program convterr

allocate ( lon(im),stat=alloc_error )
if( alloc_error /= 0 ) then
print*,'Program could not allocate space for landfrac'
print*,'Program could not allocate space for lon'
stop
end if

allocate ( lat(jm),stat=alloc_error )
if( alloc_error /= 0 ) then
print*,'Program could not allocate space for landfrac'
print*,'Program could not allocate space for lat'
stop
end if

Expand Down Expand Up @@ -153,15 +167,15 @@ program convterr

WRITE(*,*) 'done reading in USGS data from netCDF file'

WRITE(*,*) "Adjustments to land fraction: Extend land fraction for Ross Ice shelf by"
WRITE(*,*) "setting all landfractions south of 79S to 1"
DO j=1,jm
IF (lat(j)<-79.0) THEN
DO i=1,im
landfrac(i,j) = 1
END DO
END IF
END DO
!WRITE(*,*) "Adjustments to land fraction: Extend land fraction for Ross Ice shelf by"
!WRITE(*,*) "setting all landfractions south of 79S to 1"
!DO j=1,jm
! IF (lat(j)<-79.0) THEN
! DO i=1,im
! landfrac(i,j) = 1
! END DO
! END IF
!END DO
Copy link
Contributor

Choose a reason for hiding this comment

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

This needs an explanation for why it's commented out

Copy link
Contributor

Choose a reason for hiding this comment

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

Or just remove it if we are not going to keep landfrac.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

in the updates all landfrac related content are removed

Copy link
Contributor

Choose a reason for hiding this comment

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

Then let's remove this code entirely.


WRITE(*,*) "compute volume for USGS raw data"
vol = 0.0
Expand Down Expand Up @@ -254,9 +268,29 @@ program convterr
idx(i,j) = icube
idy(i,j) = jcube
idp(i,j) = ipanel
!WRITE(*,*) "terr_cube(icube,jcube,ipanel) at j,i = ",icube,jcube,ipanel,j,i,terr_cube(icube,jcube,ipanel)
END DO
END DO



DO k=1,6
DO j=1,ncube
DO i=1,ncube
IF (ABS(weight(i,j,k))<1.0E-9) THEN
WRITE(*,*) "there is no lat-lon grid point in cubed sphere cell ",i,j,k
WRITE(*,*) "fatal error"
write(*,*) "weight ",i,j,k,weight(i,j,k)
STOP
ELSE
terr_cube (i,j,k) = terr_cube (i,j,k)/weight(i,j,k)
!+++ARH
Copy link
Contributor

Choose a reason for hiding this comment

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

I'd drop these signatures, they clutter up the code and don't add value.

Copy link
Contributor

Choose a reason for hiding this comment

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

definitely drop these code tags, I can't stand clutter like this!

landfrac_cube (i,j,k) = landfrac_cube (i,j,k)/weight(i,j,k)
Copy link
Contributor

Choose a reason for hiding this comment

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

Do we need to keep landfrac around?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

in the updates, the second commit removed all landfrac related content

for my own experience so far landfrac is not needed so i vote to remove it, but it can be reverted if others want it?

Copy link
Contributor

Choose a reason for hiding this comment

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

let's drop it and then add it back if something breaks.

!---ARH
Copy link
Contributor

Choose a reason for hiding this comment

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

let's remove these "ARH" tags - they don't serve any purpose

END IF
END DO
END DO
END DO
WRITE(*,*) "min/max value of terr_cube:", MINVAL(terr_cube), MAXVAL(terr_cube)
!
! compute volume of topography on cubed-sphere
Expand Down Expand Up @@ -306,7 +340,7 @@ program convterr
!
! write data to NetCDF file
!
CALL wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube)
CALL wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube,raw_latlon_data_file,output_file)
DEALLOCATE(weight,terr,landfrac,idx,idy,idp,lat,lon)
WRITE(*,*) "done writing cubed sphere data"
end program convterr
Expand Down Expand Up @@ -431,7 +465,7 @@ END SUBROUTINE CubedSphereABPFromRLL
!
! write netCDF file
!
subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube)
subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube,raw_latlon_data_file,output_file)
use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
# include <netcdf.inc>
Expand All @@ -441,6 +475,7 @@ subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube)
!
integer, intent(in) :: ncube
real (r8), dimension(6*ncube*ncube), intent(in) :: terr_cube,landfrac_cube,sgh30_cube
character(len=1024) :: raw_latlon_data_file, git_http, tmp_string, output_file
!
! Local variables
!
Expand Down Expand Up @@ -468,7 +503,6 @@ subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube)
integer, dimension(2) :: nc_dims2_id ! netCDF dim id array for 2-d arrays
integer :: grid_dims

character(18), parameter :: grid_file_out = 'USGS-topo-cube.nc'
character(90), parameter :: grid_name = 'equi-angular gnomonic cubed sphere grid'

character (len=32) :: fout ! NetCDF output file
Expand Down Expand Up @@ -510,13 +544,31 @@ subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube)
end do
end do

WRITE(*,*) "Create NetCDF file for output"
ncstat = nf_create (grid_file_out, NF_64BIT_OFFSET,nc_grid_id)
WRITE(*,*) "Create NetCDF file for output: ", TRIM(output_file)
!ncstat = nf_create (TRIM(output_file), NF_64BIT_OFFSET,nc_grid_id)
whannah1 marked this conversation as resolved.
Show resolved Hide resolved
ncstat = nf_create (TRIM(output_file), NF_64BIT_DATA,nc_grid_id)
call handle_err(ncstat)

ncstat = nf_put_att_text (nc_grid_id, NF_GLOBAL, 'title',len_trim(grid_name), grid_name)
call handle_err(ncstat)

ncstat = nf_put_att_text (nc_grid_id, NF_GLOBAL, 'raw_topo',len_trim(raw_latlon_data_file), TRIM(raw_latlon_data_file))
call handle_err(ncstat)

git_http='https://github.com/E3SM.git (modified according to https://github.com/NCAR/Topo.git)'

ncstat = nf_put_att_text (nc_grid_id, NF_GLOBAL, 'source_code',len_trim(git_http), TRIM(git_http))
call handle_err(ncstat)

call DATE_AND_TIME(DATE=datestring)
tmp_string = 'Written on date: ' // datestring
status = nf_put_att_text (nc_grid_id,NF_GLOBAL,'history',len_trim(tmp_string), TRIM(tmp_string))
call handle_err(ncstat)

tmp_string='Peter Hjort Lauritzen (NCAR)'
ncstat = nf_put_att_text (nc_grid_id, NF_GLOBAL, 'author',len_trim(tmp_string), TRIM(tmp_string))
call handle_err(ncstat)

WRITE(*,*) "define grid size dimension"
ncstat = nf_def_dim (nc_grid_id, 'grid_size', 6*ncube*ncube, nc_gridsize_id)
call handle_err(ncstat)
Expand Down
5 changes: 5 additions & 0 deletions components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
&binparams
raw_latlon_data_file = '/p/lustre2/zhang73/DATA/DEM/S5P_OPER_REF_DEM_15_00000000T000000_99999999T999999_20160111T104226.NCL_24-3.nc'
output_file = 'S5P_OPER_REF_DEM_15_24-3_cube12000.nc'
ncube=12000
Copy link
Contributor

Choose a reason for hiding this comment

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

are we sure we want this to be the "default"?
Maybe it would be better to create a few different namelist files with the ncube value in the file name so it forces the user to make a decision on which resolution they need rather than blindly using whatever is set up in the single namelist file?

Copy link
Contributor

Choose a reason for hiding this comment

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

These namelists are small enough we could just make them input (command line) arguments as well. I'd prefer that to namelists personally.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

yes it shouldn't be the default, i just picked the one for my test

please check the updates which use the command arguments

Copy link
Contributor

Choose a reason for hiding this comment

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

Maybe I'm missing something, but I don't see any changes where command line argument support was added?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

this nl had been deleted since 03/21 commits, and the command line arguments are put in bin_to_cube.F90:

  iargc = command_argument_count()
  if (iargc /= 3) then
      print *, "Usage: <executable> raw_latlon_data_file output_file ncube"
      stop
  end if

  call get_command_argument(1, raw_latlon_data_file)
  call get_command_argument(2, output_file)
  call get_command_argument(3, string_ncube)
  read(string_ncube, *) ncube

/