From af6cbb547ce0c16ff8c330da1d671d5e4f07b8d2 Mon Sep 17 00:00:00 2001 From: jsbamboo Date: Wed, 13 Mar 2024 13:54:57 -0700 Subject: [PATCH 1/4] update bin_to_cube according to https://github.com/NCAR/Topo.git --- .../topo_tool/bin_to_cube/bin_to_cube.F90 | 88 +++++++++++++++---- .../topo_tool/bin_to_cube/bin_to_cube.nl | 5 ++ 2 files changed, 75 insertions(+), 18 deletions(-) mode change 100644 => 100755 components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 create mode 100755 components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl diff --git a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 old mode 100644 new mode 100755 index 2c48727b7ca5..e3f1b9388b19 --- a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 +++ b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 @@ -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 @@ -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) @@ -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 @@ -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 WRITE(*,*) "compute volume for USGS raw data" vol = 0.0 @@ -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 + landfrac_cube (i,j,k) = landfrac_cube (i,j,k)/weight(i,j,k) +!---ARH + 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 @@ -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 @@ -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 @@ -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 ! @@ -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 @@ -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) + 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) diff --git a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl new file mode 100755 index 000000000000..6ad7669a9e88 --- /dev/null +++ b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl @@ -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 +/ From 80817ee5c6f7588fee35e5885d94b3a19c2285ec Mon Sep 17 00:00:00 2001 From: jsbamboo Date: Thu, 21 Mar 2024 21:23:41 -0700 Subject: [PATCH 2/4] make input parameters to command arguments and clean up signatures --- .../topo_tool/bin_to_cube/bin_to_cube.F90 | 29 ++++++++++--------- .../topo_tool/bin_to_cube/bin_to_cube.nl | 5 ---- 2 files changed, 16 insertions(+), 18 deletions(-) delete mode 100755 components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl diff --git a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 index e3f1b9388b19..5b677b032b24 100755 --- a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 +++ b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 @@ -20,7 +20,6 @@ program convterr integer :: im, jm integer :: ncube !dimension of cubed-sphere grid -! integer, parameter :: ncube = 540 !dimension of cubed-sphere grid ! integer, parameter :: ncube = 361 ! for debugging integer*2, allocatable, dimension(:,:) :: terr ! global 30-sec terrain data @@ -72,18 +71,24 @@ program convterr real(r8) :: vol,dx_rad,vol_cube,area_latlon,darea_latlon ! latitude array real(r8), allocatable, dimension(:,:) :: darea_cube - INTEGER :: UNIT + INTEGER :: iargc - character(len=1024) :: raw_latlon_data_file,output_file + character(len=1024) :: raw_latlon_data_file,output_file,string_ncube - 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 + iargc = command_argument_count() + if (iargc /= 3) then + print *, "Usage: 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 + + write(*,*) "Raw lat/lon data file: ", trim(raw_latlon_data_file) + write(*,*) "Output file: ", trim(output_file) + write(*,*) "ncube: ", ncube ! ! read in USGS data from netCDF file @@ -284,9 +289,7 @@ program convterr STOP ELSE terr_cube (i,j,k) = terr_cube (i,j,k)/weight(i,j,k) -!+++ARH landfrac_cube (i,j,k) = landfrac_cube (i,j,k)/weight(i,j,k) -!---ARH END IF END DO END DO diff --git a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl deleted file mode 100755 index 6ad7669a9e88..000000000000 --- a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.nl +++ /dev/null @@ -1,5 +0,0 @@ -&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 -/ From 064de8105e0638c63d15f0eb986e7f0b1a2024c8 Mon Sep 17 00:00:00 2001 From: jsbamboo Date: Thu, 21 Mar 2024 23:03:31 -0700 Subject: [PATCH 3/4] remove landfrac in bin_to_cube and cube_to_target; correct some print message --- .../topo_tool/bin_to_cube/bin_to_cube.F90 | 55 +------- .../cube_to_target/cube_to_target.F90 | 123 +++--------------- 2 files changed, 23 insertions(+), 155 deletions(-) diff --git a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 index 5b677b032b24..c81fa9bced48 100755 --- a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 +++ b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 @@ -23,7 +23,6 @@ program convterr ! integer, parameter :: ncube = 361 ! for debugging integer*2, allocatable, dimension(:,:) :: terr ! global 30-sec terrain data - integer*1, allocatable, dimension(:,:) :: landfrac ! global 30-sec land fraction integer :: alloc_error,dealloc_error integer :: i,j,n,k,index ! index @@ -43,7 +42,7 @@ program convterr real(r8) :: alpha, beta,da,wt,dlat integer :: ipanel,icube,jcube - real(r8), allocatable, dimension(:,:,:) :: weight,terr_cube,landfrac_cube,sgh30_cube + real(r8), allocatable, dimension(:,:,:) :: weight,terr_cube,sgh30_cube integer , allocatable, dimension(:,:) :: idx,idy,idp ! real(r8) :: dx,dy @@ -110,12 +109,6 @@ program convterr WRITE(*,*) "lon-lat dimensions: ",im,jm - allocate ( landfrac(im,jm),stat=alloc_error ) - if( alloc_error /= 0 ) then - print*,'Program could not allocate space for landfrac' - stop - end if - allocate ( terr(im,jm),stat=alloc_error ) if( alloc_error /= 0 ) then print*,'Program could not allocate space for terr' @@ -135,15 +128,6 @@ program convterr end if terr = -32768 ! integer*2 - landfrac = -99.0 - - status = NF_INQ_VARID(ncid, 'landfract', landid) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - - status = NF_GET_VAR_INT1(ncid, landid,landfrac) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - WRITE(*,*) "min/max of 30sec land fraction",MINVAL(landfrac),MAXVAL(landfrac) - status = NF_INQ_VARID(ncid, 'htopo', topoid) IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) @@ -172,16 +156,6 @@ 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(*,*) "compute volume for USGS raw data" vol = 0.0 dx = (lon(2)-lon(1)) @@ -216,12 +190,6 @@ program convterr stop end if terr_cube = 0.0 - allocate ( landfrac_cube(ncube,ncube,6),stat=alloc_error ) - if( alloc_error /= 0 ) then - print*,'Program could not allocate space for terr_cube' - stop - end if - landfrac_cube = 0.0 allocate ( idx(im,jm),stat=alloc_error ) if( alloc_error /= 0 ) then print*,'Program could not allocate space for idx' @@ -266,7 +234,6 @@ program convterr weight(icube,jcube,ipanel) = weight(icube,jcube,ipanel)+wt ! terr_cube (icube,jcube,ipanel) = terr_cube (icube,jcube,ipanel)+wt*DBLE(terr(i,j)) - landfrac_cube(icube,jcube,ipanel) = landfrac_cube(icube,jcube,ipanel)+wt*DBLE(landfrac(i,j)) ! ! save "index-association" for variance computation ! @@ -289,7 +256,6 @@ program convterr STOP ELSE terr_cube (i,j,k) = terr_cube (i,j,k)/weight(i,j,k) - landfrac_cube (i,j,k) = landfrac_cube (i,j,k)/weight(i,j,k) END IF END DO END DO @@ -343,8 +309,8 @@ program convterr ! ! write data to NetCDF file ! - 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) + CALL wrt_cube(ncube,terr_cube,sgh30_cube,raw_latlon_data_file,output_file) + DEALLOCATE(weight,terr,idx,idy,idp,lat,lon) WRITE(*,*) "done writing cubed sphere data" end program convterr @@ -468,7 +434,7 @@ END SUBROUTINE CubedSphereABPFromRLL ! ! write netCDF file ! -subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube,raw_latlon_data_file,output_file) +subroutine wrt_cube(ncube,terr_cube,sgh30_cube,raw_latlon_data_file,output_file) use shr_kind_mod, only: r8 => shr_kind_r8 implicit none # include @@ -477,7 +443,7 @@ subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube,raw_latlon_data_fil ! Dummy arguments ! integer, intent(in) :: ncube - real (r8), dimension(6*ncube*ncube), intent(in) :: terr_cube,landfrac_cube,sgh30_cube + real (r8), dimension(6*ncube*ncube), intent(in) :: terr_cube,sgh30_cube character(len=1024) :: raw_latlon_data_file, git_http, tmp_string, output_file ! ! Local variables @@ -499,7 +465,6 @@ subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube,raw_latlon_data_fil integer :: nc_grdcntrlat_id ! netCDF grid center lat id integer :: nc_grdcntrlon_id ! netCDF grid center lon id integer :: nc_terr_id - integer :: nc_landfrac_id integer :: nc_var_id @@ -602,13 +567,6 @@ subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube,raw_latlon_data_fil ncstat = nf_put_att_text (nc_grid_id, nc_terr_id, 'units',1, 'm') call handle_err(ncstat) - WRITE(*,*) "define landfrac_cube array" - ncstat = nf_def_var (nc_grid_id, 'LANDFRAC', NF_DOUBLE,1, nc_gridsize_id, nc_landfrac_id) - call handle_err(ncstat) - ncstat = nf_put_att_text (nc_grid_id, nc_landfrac_id, 'long_name',70,& - 'land ocean transition mask: ocean (0), continent (1), transition (0-1)') - call handle_err(ncstat) - WRITE(*,*) "define sgh30_cube array" ncstat = nf_def_var (nc_grid_id, 'SGH30', NF_DOUBLE,1, nc_gridsize_id, nc_var_id) call handle_err(ncstat) @@ -641,9 +599,6 @@ subroutine wrt_cube(ncube,terr_cube,landfrac_cube,sgh30_cube,raw_latlon_data_fil ncstat = nf_put_var_double(nc_grid_id, nc_terr_id, terr_cube) call handle_err(ncstat) - ncstat = nf_put_var_double(nc_grid_id, nc_landfrac_id, landfrac_cube) - call handle_err(ncstat) - ncstat = nf_put_var_double(nc_grid_id, nc_var_id, sgh30_cube) call handle_err(ncstat) diff --git a/components/eam/tools/topo_tool/cube_to_target/cube_to_target.F90 b/components/eam/tools/topo_tool/cube_to_target/cube_to_target.F90 index 97d23fa8476f..c16f6a6e4758 100644 --- a/components/eam/tools/topo_tool/cube_to_target/cube_to_target.F90 +++ b/components/eam/tools/topo_tool/cube_to_target/cube_to_target.F90 @@ -35,10 +35,6 @@ program convterr ! logical :: lexternal_smooth_terr = .TRUE. ! - ! set PHIS=0.0 if LANDFRAC<0.01 - ! - logical :: lzero_out_ocean_point_phis = .FALSE. - ! ! For internal smoothing (experimental at this point) ! =================================================== ! @@ -63,7 +59,7 @@ program convterr integer :: im, jm, ncoarse integer :: ncube !dimension of cubed-sphere grid - real(r8), allocatable, dimension(:) :: landfrac, terr, sgh30 + real(r8), allocatable, dimension(:) :: terr, sgh30 real(r8), allocatable, dimension(:) :: terr_coarse !for internal smoothing integer :: alloc_error,dealloc_error @@ -89,7 +85,7 @@ program convterr real(r8) :: wt,dlat integer :: ipanel,icube,jcube - real(r8), allocatable, dimension(:,:,:) :: weight,terr_cube,landfrac_cube,sgh30_cube + real(r8), allocatable, dimension(:,:,:) :: weight,terr_cube,sgh30_cube integer, allocatable, dimension(:,:) :: idx,idy,idp integer :: npatch, isub,jsub, itmp, iplm1,jmin,jmax real(r8) :: sum,dx,scale,dmax,arad,jof,term,s1,c1,clon,iof,dy,s2,c2,dist @@ -104,7 +100,7 @@ program convterr integer :: src_grid_dim ! for netCDF weight file integer :: n_a,n_b,n_s,n_aid,n_bid,n_sid integer :: count - real(r8), allocatable, dimension(:) :: landfrac_target, terr_target, sgh30_target, sgh_target + real(r8), allocatable, dimension(:) :: terr_target, sgh30_target, sgh_target real(r8), allocatable, dimension(:) :: area_target ! ! this is only used if target grid is a lat-lon grid @@ -304,27 +300,12 @@ program convterr ncube = INT(SQRT(DBLE(n/6))) WRITE(*,*) "cubed-sphere dimension, ncube: ",ncube - ! - ! read LANDFRAC - ! - allocate ( landfrac(n),stat=alloc_error ) - if( alloc_error /= 0 ) then - print*,'Program could not allocate space for landfrac' - stop - end if - - status = NF_INQ_VARID(ncid, 'LANDFRAC', landid) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - - status = NF_GET_VAR_DOUBLE(ncid, landid,landfrac) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - WRITE(*,*) "min/max of landfrac",MINVAL(landfrac),MAXVAL(landfrac) ! ! read terr ! allocate ( terr(n),stat=alloc_error ) if( alloc_error /= 0 ) then - print*,'Program could not allocate space for landfrac' + print*,'Program could not allocate space for terr' stop end if @@ -339,7 +320,7 @@ program convterr ! allocate ( sgh30(n),stat=alloc_error ) if( alloc_error /= 0 ) then - print*,'Program could not allocate space for landfrac' + print*,'Program could not allocate space for sgh30' stop end if @@ -366,11 +347,6 @@ program convterr print*,'Program could not allocate space for terr_target' stop end if - allocate (landfrac_target(ntarget),stat=alloc_error ) - if( alloc_error /= 0 ) then - print*,'Program could not allocate space for landfrac_target' - stop - end if allocate (sgh30_target(ntarget),stat=alloc_error ) if( alloc_error /= 0 ) then print*,'Program could not allocate space for sgh30_target' @@ -378,7 +354,6 @@ program convterr end if allocate (area_target(ntarget),stat=alloc_error ) terr_target = 0.0 - landfrac_target = 0.0 sgh30_target = 0.0 area_target = 0.0 @@ -403,7 +378,6 @@ program convterr wt = weights_all(count,1) terr_target (i) = terr_target (i) + wt*terr (ii)/area_target(i) - landfrac_target (i) = landfrac_target (i) + wt*landfrac (ii)/area_target(i) sgh30_target (i) = sgh30_target (i) + wt*sgh30 (ii)/area_target(i) tmp = tmp+wt*terr(ii) @@ -445,7 +419,6 @@ program convterr WRITE(*,*) WRITE(*,*) "min/max of unsmoothed terr_target : ",MINVAL(terr_target ),MAXVAL(terr_target ) - WRITE(*,*) "min/max of landfrac_target : ",MINVAL(landfrac_target),MAXVAL(landfrac_target) WRITE(*,*) "min/max of var30_target : ",MINVAL(sgh30_target ),MAXVAL(sgh30_target ) ! ! compute mean height (globally) of topography about sea-level for target grid unfiltered elevation @@ -572,7 +545,7 @@ program convterr WRITE(*,*) "resolution of coarse grid", 90.0/ncube_coarse allocate ( terr_coarse(ncoarse),stat=alloc_error ) if( alloc_error /= 0 ) then - print*,'Program could not allocate space for landfrac' + print*,'Program could not allocate space for terr_coarse' stop end if WRITE(*,*) "coarsening" @@ -837,10 +810,6 @@ program convterr ! WRITE(*,*) "min/max of terr_target : ",MINVAL(terr_target),MAXVAL(terr_target) - if (lzero_out_ocean_point_phis) then - WRITE(*,*) "if ocean mask PHIS=0.0" - end if - sgh_target=0.0 do count=1,jall @@ -856,9 +825,6 @@ program convterr wt = weights_all(count,1) - if (lzero_out_ocean_point_phis.AND.landfrac_target(i).lt.0.01_r8) then - terr_target(i) = 0.0_r8 !5*terr_target(i) - end if sgh_target(i) = sgh_target(i)+wt*((terr_target(i)-terr(ii))**2)/area_target(i) end do @@ -868,7 +834,6 @@ program convterr ! zero out small values ! DO i=1,ntarget - IF (landfrac_target(i)<.001_r8) landfrac_target(i) = 0.0 IF (sgh_target(i)<0.5) sgh_target(i) = 0.0 IF (sgh30_target(i)<0.5) sgh30_target(i) = 0.0 END DO @@ -877,17 +842,17 @@ program convterr WRITE(*,*) "min/max of sgh_target : ",MINVAL(sgh_target),MAXVAL(sgh_target) WRITE(*,*) "min/max of sgh30_target : ",MINVAL(sgh30_target),MAXVAL(sgh30_target) - DEALLOCATE(terr,weights_all,weights_eul_index_all,landfrac) + DEALLOCATE(terr,weights_all,weights_eul_index_all) IF (ltarget_latlon) THEN - CALL wrtncdf_rll(nlon,nlat,lpole,ntarget,terr_target,landfrac_target,sgh_target,sgh30_target,& + CALL wrtncdf_rll(nlon,nlat,lpole,ntarget,terr_target,sgh_target,sgh30_target,& target_center_lon,target_center_lat,.true.,output_topography_file) ELSE - CALL wrtncdf_unstructured(ntarget,terr_target,landfrac_target,sgh_target,sgh30_target,& + CALL wrtncdf_unstructured(ntarget,terr_target,sgh_target,sgh30_target,& target_center_lon,target_center_lat,output_topography_file) END IF - DEALLOCATE(terr_target,landfrac_target,sgh30_target,sgh_target) + DEALLOCATE(terr_target,sgh30_target,sgh_target) end program convterr @@ -1016,7 +981,7 @@ end subroutine usage ! ! ! -subroutine wrtncdf_unstructured(n,terr,landfrac,sgh,sgh30,lon,lat,fout) +subroutine wrtncdf_unstructured(n,terr,sgh,sgh30,lon,lat,fout) use shr_kind_mod, only: r8 => shr_kind_r8 implicit none @@ -1026,7 +991,7 @@ subroutine wrtncdf_unstructured(n,terr,landfrac,sgh,sgh30,lon,lat,fout) ! Dummy arguments ! integer, intent(in) :: n - real(r8),dimension(n) , intent(in) :: terr, landfrac,sgh,sgh30,lon, lat + real(r8),dimension(n) , intent(in) :: terr, sgh,sgh30,lon, lat ! ! Local variables ! @@ -1035,7 +1000,7 @@ subroutine wrtncdf_unstructured(n,terr,landfrac,sgh,sgh30,lon,lat,fout) integer :: lonid, lonvid integer :: latid, latvid integer :: terrid,nid - integer :: terrdim,landfracid,sghid,sgh30id + integer :: terrdim,sghid,sgh30id integer :: status ! return value for error control of netcdf routin integer :: i,j integer, dimension(2) :: nc_lat_vid,nc_lon_vid @@ -1062,9 +1027,6 @@ subroutine wrtncdf_unstructured(n,terr,landfrac,sgh,sgh30,lon,lat,fout) status = nf_def_var (foutid,'PHIS', NF_DOUBLE, 1, nid, terrid) if (status .ne. NF_NOERR) call handle_err(status) - status = nf_def_var (foutid,'LANDFRAC', NF_DOUBLE, 1, nid, landfracid) - if (status .ne. NF_NOERR) call handle_err(status) - status = nf_def_var (foutid,'SGH', NF_DOUBLE, 1, nid, sghid) if (status .ne. NF_NOERR) call handle_err(status) @@ -1100,11 +1062,6 @@ subroutine wrtncdf_unstructured(n,terr,landfrac,sgh,sgh30,lon,lat,fout) status = nf_put_att_text (foutid, sgh30id, 'units' , 1, 'm') ! status = nf_put_att_text (foutid, sgh30id, 'filter' , 4, 'none') - status = nf_put_att_double (foutid, landfracid, 'missing_value', nf_double, 1, fillvalue) - status = nf_put_att_double (foutid, landfracid, '_FillValue' , nf_double, 1, fillvalue) - status = nf_put_att_text (foutid, landfracid, 'long_name', 21, 'gridbox land fraction') - ! status = nf_put_att_text (foutid, landfracid, 'filter', 40, 'area averaged from 30-sec USGS raw data') - status = nf_put_att_text (foutid,latvid,'long_name', 8, 'latitude') if (status .ne. NF_NOERR) call handle_err(status) @@ -1141,11 +1098,6 @@ subroutine wrtncdf_unstructured(n,terr,landfrac,sgh,sgh30,lon,lat,fout) if (status .ne. NF_NOERR) call handle_err(status) print*,"done writing terrain data" - print*,"writing landfrac data",MINVAL(landfrac),MAXVAL(landfrac) - status = nf_put_var_double (foutid, landfracid, landfrac) - if (status .ne. NF_NOERR) call handle_err(status) - print*,"done writing landfrac data" - print*,"writing sgh data",MINVAL(sgh),MAXVAL(sgh) status = nf_put_var_double (foutid, sghid, sgh) if (status .ne. NF_NOERR) call handle_err(status) @@ -1179,7 +1131,7 @@ end subroutine wrtncdf_unstructured ! !************************************************************** ! -subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon,lat,lprepare_fv_smoothing_routine,fout) +subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,sgh_in,sgh30_in,lon,lat,lprepare_fv_smoothing_routine,fout) use shr_kind_mod, only: r8 => shr_kind_r8 implicit none @@ -1193,7 +1145,7 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon ! lprepare_fv_smoothing_routine is to make a NetCDF file that can be used with the CAM-FV smoothing software ! logical , intent(in) :: lpole,lprepare_fv_smoothing_routine - real(r8),dimension(n) , intent(in) :: terr_in, landfrac_in,sgh_in,sgh30_in,lon, lat + real(r8),dimension(n) , intent(in) :: terr_in, sgh_in,sgh30_in,lon, lat ! ! Local variables ! @@ -1202,7 +1154,7 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon integer :: lonid, lonvid integer :: latid, latvid integer :: terrid,nid - integer :: terrdim,landfracid,sghid,sgh30id + integer :: terrdim,sghid,sgh30id integer :: status ! return value for error control of netcdf routin integer :: i,j integer, dimension(2) :: nc_lat_vid,nc_lon_vid @@ -1214,8 +1166,8 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon real(r8),dimension(nlon) :: lonar ! longitude array real(r8),dimension(nlat) :: latar ! latitude array - integer, dimension(2) :: htopodim,landfdim,sghdim,sgh30dim - real(r8),dimension(n) :: terr, landfrac,sgh,sgh30 + integer, dimension(2) :: htopodim,sghdim,sgh30dim + real(r8),dimension(n) :: terr, sgh,sgh30 IF (nlon*nlat.NE.n) THEN WRITE(*,*) "inconsistent input for wrtncdf_rll" @@ -1239,7 +1191,6 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon terr = terr_in sgh=sgh_in sgh30 =sgh30_in - landfrac = landfrac_in if (lpole) then write(*,*) "average pole control volume" @@ -1294,23 +1245,6 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon end do sgh30(n-(nlon+1):n) = ave/DBLE(nlon) - ! - ! North pole - landfrac - ! - ave = 0.0 - do i=1,nlon - ave = ave + landfrac_in(i) - end do - landfrac(1:nlon) = ave/DBLE(nlon) - ! - ! South pole - ! - ave = 0.0 - do i=n-(nlon+1),n - ave = ave + landfrac_in(i) - end do - landfrac(n-(nlon+1):n) = ave/DBLE(nlon) - end if @@ -1343,17 +1277,6 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon end if if (status .ne. NF_NOERR) call handle_err(status) - landfdim(1)=lonid - landfdim(2)=latid - - if (lprepare_fv_smoothing_routine) then - status = nf_def_var (foutid,'ftopo', NF_DOUBLE, 2, landfdim, landfracid) - else - status = nf_def_var (foutid,'LANDFRAC', NF_DOUBLE, 2, landfdim, landfracid) - end if - - if (status .ne. NF_NOERR) call handle_err(status) - sghdim(1)=lonid sghdim(2)=latid @@ -1396,11 +1319,6 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon status = nf_put_att_text (foutid, sgh30id, 'units' , 1, 'm') status = nf_put_att_text (foutid, sgh30id, 'filter' , 4, 'none') - status = nf_put_att_double (foutid, landfracid, 'missing_value', nf_double, 1, fillvalue) - status = nf_put_att_double (foutid, landfracid, '_FillValue' , nf_double, 1, fillvalue) - status = nf_put_att_text (foutid, landfracid, 'long_name', 21, 'gridbox land fraction') - status = nf_put_att_text (foutid, landfracid, 'filter', 40, 'area averaged from 30-sec USGS raw data') - status = nf_put_att_text (foutid,latvid,'long_name', 8, 'latitude') if (status .ne. NF_NOERR) call handle_err(status) @@ -1441,11 +1359,6 @@ subroutine wrtncdf_rll(nlon,nlat,lpole,n,terr_in,landfrac_in,sgh_in,sgh30_in,lon if (status .ne. NF_NOERR) call handle_err(status) print*,"done writing terrain data" - print*,"writing landfrac data",MINVAL(landfrac),MAXVAL(landfrac) - status = nf_put_var_double (foutid, landfracid, landfrac) - if (status .ne. NF_NOERR) call handle_err(status) - print*,"done writing landfrac data" - print*,"writing sgh data",MINVAL(sgh),MAXVAL(sgh) status = nf_put_var_double (foutid, sghid, sgh) if (status .ne. NF_NOERR) call handle_err(status) From 90d4e2da8b9b12f69bfd4120f855374a5c27985e Mon Sep 17 00:00:00 2001 From: Walter Hannah Date: Mon, 15 Apr 2024 09:02:23 -0600 Subject: [PATCH 4/4] remove comment in bin_to_cube --- components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 index c81fa9bced48..cabea18d44ba 100755 --- a/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 +++ b/components/eam/tools/topo_tool/bin_to_cube/bin_to_cube.F90 @@ -513,7 +513,6 @@ subroutine wrt_cube(ncube,terr_cube,sgh30_cube,raw_latlon_data_file,output_file) end do WRITE(*,*) "Create NetCDF file for output: ", TRIM(output_file) - !ncstat = nf_create (TRIM(output_file), NF_64BIT_OFFSET,nc_grid_id) ncstat = nf_create (TRIM(output_file), NF_64BIT_DATA,nc_grid_id) call handle_err(ncstat)