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 cube_to_target to calculate orographic shape parameter #6876

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

update cube_to_target to calculate orographic shape parameter #6876

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
7978edb
add orographic_shape_methods.F90
whannah1 Jan 6, 2025
688f0b9
update makefile
whannah1 Jan 6, 2025
36170b0
initial cube_to_target update for orographic shape params
whannah1 Jan 6, 2025
279d669
update file output methods for orographic shape data
whannah1 Jan 7, 2025
8653588
update help text for cube_to_target
whannah1 Jan 7, 2025
d7651ec
misc updates and fixes for oro shape mods
whannah1 Jan 8, 2025
d6089e0
cube_to_target fixes - tested and working on chrysalis
whannah1 Jan 9, 2025
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
4 changes: 2 additions & 2 deletions components/eam/tools/topo_tool/cube_to_target/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,14 +34,14 @@ UNAMEM := $(findstring CRAY,$(shell uname -m))
# Default rules and macros
#------------------------------------------------------------------------

OBJS := reconstruct.o remap.o cube_to_target.o shr_kind_mod.o
OBJS := reconstruct.o remap.o cube_to_target.o shr_kind_mod.o orographic_shape_methods.o

$(EXEDIR)/$(EXENAME): $(OBJS)
$(FC) -o $@ $(OBJS) $(LDFLAGS)

clean:
$(RM) -f $(OBJS) *.mod $(EXEDIR)/$(EXENAME)

cube_to_target.o: shr_kind_mod.o remap.o reconstruct.o
cube_to_target.o: shr_kind_mod.o remap.o reconstruct.o orographic_shape_methods.o
remap.o:
reconstruct.o: remap.o
250 changes: 230 additions & 20 deletions components/eam/tools/topo_tool/cube_to_target/cube_to_target.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
program convterr
use shr_kind_mod, only: r8 => shr_kind_r8
use reconstruct
use orographic_shape_methods
implicit none
#include <netcdf.inc>

Expand All @@ -18,7 +19,9 @@ program convterr
! USER SETTINGS BELOW
!
!**************************************
!
!
! flag to enable calculation of orographic shape parameters
logical :: calc_orographic_shape_params = .FALSE.
!
! If smoothed PHIS is available SGH needs to be recomputed to account for the
! sub-grid-scale variability introduced by the smoothing; This will be set
Expand Down Expand Up @@ -155,6 +158,18 @@ program convterr
character(len=512) :: input_topography_file
character(len=512) :: output_topography_file
character(len=512) :: smoothed_topography_file

! variables for orographic shape parameters
real(r8), allocatable, dimension(:,:) :: oa_target
real(r8), allocatable, dimension(:) :: oc_target
real(r8), allocatable, dimension(:,:) :: ol_target
integer, allocatable, dimension(:) :: indexb !max indice dimension
real(r8), allocatable, dimension(:,:,:) :: terrout
real(r8), allocatable, dimension(:,:) :: dxy
real(r8), allocatable, dimension(:) :: lat_terr
real(r8), allocatable, dimension(:) :: lon_terr
integer :: nvar_dirOA
integer :: nvar_dirOL

!
! turn extra debugging on/off
Expand All @@ -165,7 +180,29 @@ program convterr

call parse_arguments(target_grid_file , input_topography_file , &
output_topography_file, smoothed_topography_file, &
lsmooth_terr )
lsmooth_terr, calc_orographic_shape_params )


!
!*********************************************************
!
! check that options are compatible
!
!*********************************************************
!

if (ltarget_latlon.and.calc_orographic_shape_params) then
print *, 'ERROR: calc_orographic_shape_params not supported for lat/lon grids.'
stop
end if

!
!*********************************************************
!
! check input files
!
!*********************************************************
!

if (lsmooth_terr) then
status = nf_open(trim(smoothed_topography_file), 0, ncid)
Expand Down Expand Up @@ -315,8 +352,41 @@ program convterr
status = NF_GET_VAR_DOUBLE(ncid, landid,terr)
IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status)
WRITE(*,*) "min/max of terr",MINVAL(terr),MAXVAL(terr)

!
! read lat/lon coordinates from topo file if shape parameters are requested
!
if (calc_orographic_shape_params) then

! read latitude coordinate
allocate ( lat_terr(n),stat=alloc_error )
if( alloc_error /= 0 ) then
print*,'Program could not allocate space for lat_terr'
stop
end if

status = NF_INQ_VARID(ncid, 'lat', landid)
IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status)

status = NF_GET_VAR_DOUBLE(ncid, landid,lat_terr)
IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status)

! read longitude coordinate
allocate ( lon_terr(n),stat=alloc_error )
if( alloc_error /= 0 ) then
print*,'Program could not allocate space for lon_terr'
stop
end if

status = NF_INQ_VARID(ncid, 'lon', landid)
IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status)

status = NF_GET_VAR_DOUBLE(ncid, landid,lon_terr)
IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status)

end if
!
! read SGH30
!
allocate ( sgh30(n),stat=alloc_error )
if( alloc_error /= 0 ) then
Expand Down Expand Up @@ -841,18 +911,79 @@ program convterr
sgh30_target = SQRT(sgh30_target)
WRITE(*,*) "min/max of sgh_target : ",MINVAL(sgh_target),MAXVAL(sgh_target)
WRITE(*,*) "min/max of sgh30_target : ",MINVAL(sgh30_target),MAXVAL(sgh30_target)

!-----------------------------------------------------------------------------
! Orographic shape parameters
if (calc_orographic_shape_params) then

nvar_dirOA = 2 ! only 2 directions needed for asymmetry (i.e. lat/lon)
nvar_dirOL = 180 ! 180 => 2 degree angular resolution for effective length

! allocate variable for orographic shape calcualtions
allocate(oa_target(ntarget,nvar_dirOA), stat=alloc_error)
allocate(oc_target(ntarget), stat=alloc_error)
allocate(ol_target(ntarget,nvar_dirOL), stat=alloc_error)
allocate(indexb(ntarget), stat=alloc_error)

! initialize allocated vairables
oa_target = 0
oc_target = 0
ol_target = 0
indexb = 0

! Orographic asymmetry
print*,"calculating orographic asymmetry..."
call orographic_asymmetry_xie2020(terr, ntarget, ncube, n, nvar_dirOA, jall, &
weights_lgr_index_all, weights_eul_index_all(:,1), &
weights_eul_index_all(:,2), weights_eul_index_all(:,3), &
weights_all, target_center_lon, target_center_lat, &
lon_terr, lat_terr, area_target, oa_target)

! Orographic convexity
print*,"calculating orographic convexity..."
call orographic_convexity_kim2005(terr, ntarget, ncube, n, jall, &
weights_lgr_index_all, weights_eul_index_all(:,1), &
weights_eul_index_all(:,2), weights_eul_index_all(:,3), &
weights_all, area_target, sgh_target, terr_target, oc_target)

! Orographic effective length
print*,"calculating orographic effective length..."
do count=1,jall
i = weights_lgr_index_all(count)
indexb(i) = indexb(i)+1
enddo
allocate(terrout(4,ntarget,maxval(indexb)), stat=alloc_error)
allocate(dxy(ntarget,nvar_dirOL), stat=alloc_error)

call orographic_efflength_xie2020(terr, ntarget, ncube, n, jall, nlon, nlat, maxval(indexb), &
nvar_dirOL, weights_lgr_index_all, weights_eul_index_all(:,1), &
weights_eul_index_all(:,2), weights_eul_index_all(:,3), &
weights_all, target_center_lon, target_center_lat, &
target_corner_lon, target_corner_lat, &
lon_terr, lat_terr, sgh_target, area_target, ol_target, &
terrout, dxy)

end if
!-----------------------------------------------------------------------------

DEALLOCATE(terr,weights_all,weights_eul_index_all)


IF (ltarget_latlon) THEN
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,sgh_target,sgh30_target,&
target_center_lon,target_center_lat,output_topography_file)
END IF
if (ltarget_latlon) then
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, sgh_target, sgh30_target, &
calc_orographic_shape_params, &
nvar_dirOA, nvar_dirOL, &
oc_target, oa_target, ol_target, &
target_center_lon, target_center_lat, output_topography_file)
end if

DEALLOCATE(terr_target,sgh30_target,sgh_target)

if (calc_orographic_shape_params) deallocate(oa_target,oc_target,ol_target,indexb,terrout,dxy)

end program convterr

Expand All @@ -862,13 +993,14 @@ end program convterr
!
subroutine parse_arguments(target_grid_file , input_topography_file , &
output_topography_file, smoothed_topography_file, &
lsmooth_terr )
lsmooth_terr, calc_orographic_shape_params )
implicit none
character(len=*), intent(inout) :: target_grid_file
character(len=*), intent(inout) :: input_topography_file
character(len=*), intent(inout) :: output_topography_file
character(len=*), intent(inout) :: smoothed_topography_file
logical, intent(inout) :: lsmooth_terr
logical, intent(inout) :: lsmooth_terr
logical, intent(inout) :: calc_orographic_shape_params

integer :: n, nargs
character(len=512) :: arg
Expand All @@ -879,6 +1011,7 @@ subroutine parse_arguments(target_grid_file , input_topography_file , &
output_topography_file = ''
smoothed_topography_file = ''
lsmooth_terr = .false.
calc_orographic_shape_params = .false.

! Get number of arguments and make sure at least some arguments were passed
nargs = iargc()
Expand Down Expand Up @@ -917,6 +1050,10 @@ subroutine parse_arguments(target_grid_file , input_topography_file , &
smoothed_topography_file = trim(arg)
lsmooth_terr = .true.
n = n + 1
case ('--add-oro-shape')
! if flag is present then calculate the orographic shape parameters
calc_orographic_shape_params = .true.
n = n + 1
case ('--help')
call usage()
stop
Expand Down Expand Up @@ -966,6 +1103,12 @@ subroutine usage()
print *, ' 3km grid (SGH30) is also downscaled, '
print *, ' but does not depend on the smoothing. '
print *, ' '
print *, ' --add-oro-shape Enable the calculation of orographic '
print *, ' shape parameters needed for certain '
print *, ' orographic drag schemes. The parameters'
print *, ' are convexivity, asymmetry, and eff. '
print *, ' length. '
print *, ' '
print *, 'DESCRIPTION: '
print *, 'This code performs rigorous remapping of topography variables on a cubed- '
print *, 'sphere grid to any target grid. The code is documented in: '
Expand All @@ -981,17 +1124,31 @@ end subroutine usage
!
!
!
subroutine wrtncdf_unstructured(n,terr,sgh,sgh30,lon,lat,fout)
subroutine wrtncdf_unstructured(n, terr, sgh, sgh30, &
calc_orographic_shape_params, &
nvar_dirOA, nvar_dirOL, &
orographic_convexity, &
Copy link
Contributor

Choose a reason for hiding this comment

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

These changes require passing the shape parameter arrays even if we are not calculating/writing them. One alternative might be to add a separate function to write these to the output file, but maybe that's not worth the effort. Just an idea.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

yea, I struggled with this. I think this approach is better than adding a new routine. This is partly because I worry about appending to a file, I've never done that and I'm not sure how to ensure the file is closed before opening it again. Alternatively, I could open/close the file outside the routines that define variables and write the data, but that sounds like a lot more work.

The current approach just passes the flag used to trigger the shape parameter calculation, so it will ignore those extra (unallocated) variables when the shape parameters are not requested.

orographic_asymmetry, &
orographic_efflength, &
lon, lat, fout)
use shr_kind_mod, only: r8 => shr_kind_r8
implicit none

# include <netcdf.inc>

!
! Dummy arguments
!
integer, intent(in) :: n
real(r8),dimension(n) , intent(in) :: terr, sgh,sgh30,lon, lat
#include <netcdf.inc>
!
! Arguments
!
integer, intent(in) :: n
real(r8),dimension(n), intent(in) :: terr
real(r8),dimension(n), intent(in) :: sgh
real(r8),dimension(n), intent(in) :: sgh30
logical, intent(in) :: calc_orographic_shape_params
integer, intent(in) :: nvar_dirOA
integer, intent(in) :: nvar_dirOL
real(r8),dimension(n), intent(in) :: orographic_convexity
real(r8),dimension(n,nvar_dirOA), intent(in) :: orographic_asymmetry
real(r8),dimension(n,nvar_dirOL), intent(in) :: orographic_efflength
real(r8),dimension(n), intent(in) :: lon
real(r8),dimension(n), intent(in) :: lat
!
! Local variables
!
Expand All @@ -1008,6 +1165,14 @@ subroutine wrtncdf_unstructured(n,terr,sgh,sgh30,lon,lat,fout)
integer :: nc_gridcorn_id, lat_vid, lon_vid

real(r8), parameter :: fillvalue = 1.d36

integer :: nvar_dirOA_id
integer :: nvar_dirOL_id
integer :: ocid
integer :: oaid
integer :: olid
integer, dimension(3) :: oadim
integer, dimension(3) :: oldim

!
! Create NetCDF file for output
Expand All @@ -1020,6 +1185,10 @@ subroutine wrtncdf_unstructured(n,terr,sgh,sgh30,lon,lat,fout)
!
status = nf_def_dim (foutid, 'ncol', n, nid)
if (status .ne. NF_NOERR) call handle_err(status)
status = nf_def_dim (foutid, 'nvar_dirOA', nvar_dirOA, nvar_dirOA_id)
if (status .ne. NF_NOERR) call handle_err(status)
status = nf_def_dim (foutid, 'nvar_dirOL', nvar_dirOL, nvar_dirOL_id)
if (status .ne. NF_NOERR) call handle_err(status)
!
! Create variable for output
!
Expand All @@ -1038,6 +1207,23 @@ subroutine wrtncdf_unstructured(n,terr,sgh,sgh30,lon,lat,fout)

status = nf_def_var (foutid,'lon', NF_DOUBLE, 1, nid, lonvid)
if (status .ne. NF_NOERR) call handle_err(status)

if (calc_orographic_shape_params) then

status = nf_def_var (foutid,'OC', NF_DOUBLE, 1, nid, ocid)
if (status .ne. NF_NOERR) call handle_err(status)

oadim(1) = nid
oadim(2) = nvar_dirOA_id
status = nf_def_var (foutid,'OA', NF_DOUBLE, 2, oadim, oaid)
if (status .ne. NF_NOERR) call handle_err(status)

oldim(1) = nid
oldim(2) = nvar_dirOL_id
status = nf_def_var (foutid,'OL', NF_DOUBLE, 2, oldim, olid)
if (status .ne. NF_NOERR) call handle_err(status)

end if

!
! Create attributes for output variables
Expand Down Expand Up @@ -1084,6 +1270,13 @@ subroutine wrtncdf_unstructured(n,terr,sgh,sgh30,lon,lat,fout)
call DATE_AND_TIME(DATE=datestring)
status = nf_put_att_text (foutid,NF_GLOBAL,'history',25, 'Written on date: ' // datestring )
if (status .ne. NF_NOERR) call handle_err(status)

if (calc_orographic_shape_params) then

status = nf_put_att_text (foutid,oaid,'note', 40, '(2)+1 in nvar_dirOA to avoid bug in io')
if (status .ne. NF_NOERR) call handle_err(status)

end if

!
! End define mode for output file
Expand Down Expand Up @@ -1117,6 +1310,23 @@ subroutine wrtncdf_unstructured(n,terr,sgh,sgh30,lon,lat,fout)
status = nf_put_var_double (foutid, lonvid, lon)
if (status .ne. NF_NOERR) call handle_err(status)
print*,"done writing lon data"

if (calc_orographic_shape_params) then

status = nf_put_var_double (foutid, ocid, orographic_convexity)
if (status .ne. NF_NOERR) call handle_err(status)
print*,"done writing orographic convexity data"

status = nf_put_var_double (foutid, oaid, orographic_asymmetry)
if (status .ne. NF_NOERR) call handle_err(status)
print*,"done writing orographic asymmetry data"

status = nf_put_var_double (foutid, olid, orographic_efflength)
if (status .ne. NF_NOERR) call handle_err(status)
print*,"done writing orographic eff length data"

end if

!
! Close output file
!
Expand Down
Loading
Loading