Skip to content

Commit

Permalink
feature/regional_grid This commit references NOAA-EMC#4
Browse files Browse the repository at this point in the history 
Remove computation of global_equiv_resol from the filter topo
program.  Retain the stand-alone version of global_equiv_resol.

Set CRES value for both esg and gfdl regional grids from
global_equiv_resol.
  • Loading branch information
@GeorgeGayno-NOAA
GeorgeGayno-NOAA committed Jul 1, 2020
1 parent e9cacf2 commit 8d01219
Show file tree
Hide file tree
Showing 3 changed files with 28 additions and 169 deletions.
156 changes: 0 additions & 156 deletions sorc/fre-nctools.fd/tools/filter_topo/filter_topo.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,6 @@ program filter_topo

!--- compute filter constants for the regional resolution
if(regional)then
call global_equiv_resol
call compute_filter_constants
endif

Expand Down Expand Up @@ -1764,161 +1763,6 @@ subroutine read_namelist

end subroutine read_namelist

!=======================================================================
! Determine the global equivalent resolution for a jim purser
! regional grid.
!=======================================================================

subroutine global_equiv_resol

use netcdf

implicit none

integer, parameter :: dp = kind(1.0d0)
real(dp), parameter :: pi_geom = 4.0*atan(1.0), &
radius_Earth = 6371000.0

character(len=50) :: tile_file
integer :: ncid, nxSG_dimid, nySG_dimid, dASG_varid
integer :: nxSG, nySG, nx, ny, RES_equiv, id_var
real(dp) :: avg_cell_size, min_cell_size, max_cell_size
real(dp), dimension(:,:), allocatable :: &
quarter_dA_ll, quarter_dA_lr, quarter_dA_ur, quarter_dA_ul, &
dASG, dA, sqrt_dA

WRITE(*,500)
WRITE(*,500) "Compute global equivalent resolution."
WRITE(*,500) "Opening NetCDF mosaic file for reading:"
WRITE(*,500) " file=", trim(grid_file)

!=======================================================================
! Obtain the grid file name from the mosaic file.
!=======================================================================

call check( nf90_open(trim(grid_file), NF90_NOWRITE, ncid) )

call check( nf90_inq_varid(ncid, 'gridfiles', id_var) )

call check ( nf90_get_var(ncid, id_var, tile_file ) )

call check ( nf90_close(ncid) )
!
!=======================================================================
!
! Open the grid file and read in the dimensions of the supergrid. The
! supergrid is a grid that has twice the resolution of the actual/compu-
! tational grid. In the file, the names of the supergrid dimensions are
! nx and ny. Here, however, we reserve those names for the dimensions
! of the actual grid (since in the FV3 code and in other data files, nx
! and ny are used to denote the dimensions of the actual grid) and in-
! stead use the variables nxSG and nySG to denote the dimensions of the
! supergrid.
!
!=======================================================================
!
WRITE(*,500)
WRITE(*,500) "Opening NetCDF grid file for reading:"
WRITE(*,500) " file=", trim(tile_file)

call check( nf90_open(trim(tile_file), NF90_NOWRITE, ncid) )

call check( nf90_inq_dimid(ncid, "nx", nxSG_dimid) )
call check( nf90_inquire_dimension(ncid, nxSG_dimid, len=nxSG) )

call check( nf90_inq_dimid(ncid, "ny", nySG_dimid) )
call check( nf90_inquire_dimension(ncid, nySG_dimid, len=nySG) )

WRITE(*,500)
WRITE(*,500) "Dimensions of supergrid are:"
WRITE(*,520) " nxSG = ", nxSG
WRITE(*,520) " nySG = ", nySG
!
!=======================================================================
!
! Read in the cell areas on the supergrid. Then add the areas of the
! four supergrid cells that make up one grid cell to obtain the cell
! areas on the actual grid.
!
!=======================================================================
!
allocate(dASG(0:nxSG-1, 0:nySG-1))
call check( nf90_inq_varid(ncid, "area", dASG_varid) )
call check( nf90_get_var(ncid, dASG_varid, dASG) )

call check ( nf90_close(ncid) )

nx = nxSG/2
ny = nySG/2

WRITE(*,500)
WRITE(*,500) "Dimensions of (actual, i.e. computational) grid are:"
WRITE(*,520) " nx = ", nx
WRITE(*,520) " ny = ", ny

allocate(quarter_dA_ll(0:nx-1, 0:ny-1))
allocate(quarter_dA_lr(0:nx-1, 0:ny-1))
allocate(quarter_dA_ul(0:nx-1, 0:ny-1))
allocate(quarter_dA_ur(0:nx-1, 0:ny-1))

quarter_dA_ll = dASG(0:nxSG-1:2, 0:nySG-1:2)
quarter_dA_lr = dASG(0:nxSG-1:2, 1:nySG-1:2)
quarter_dA_ur = dASG(1:nxSG-1:2, 1:nySG-1:2)
quarter_dA_ul = dASG(1:nxSG-1:2, 0:nySG-1:2)

deallocate(dASG)

allocate(dA(0:nx-1, 0:ny-1))
allocate(sqrt_dA(0:nx-1, 0:ny-1))

dA = quarter_dA_ll + quarter_dA_lr + quarter_dA_ur + quarter_dA_ul

deallocate(quarter_dA_ll, quarter_dA_lr, quarter_dA_ur, quarter_dA_ul)

!=======================================================================
!
! Calculate a typical/representative cell size for each cell by taking
! the square root of the area of the cell. Then calculate the minimum,
! maximum, and average cell sizes over the whole grid.
!
!=======================================================================
!
sqrt_dA = sqrt(dA)
deallocate(dA)
min_cell_size = minval(sqrt_dA)
max_cell_size = maxval(sqrt_dA)
avg_cell_size = sum(sqrt_dA)/(nx*ny)
deallocate(sqrt_dA)

WRITE(*,500)
WRITE(*,500) "Minimum, maximum, and average cell sizes are (based on square"
WRITE(*,500) "root of cell area):"
WRITE(*,530) " min_cell_size = ", min_cell_size
WRITE(*,530) " max_cell_size = ", max_cell_size
WRITE(*,530) " avg_cell_size = ", avg_cell_size
!
!=======================================================================
!
! Use the average cell size to calculate an equivalent global uniform
! cubed-sphere resolution (in units of number of cells) for the regional
! grid. This is the RES that a global uniform (i.e. stretch factor of
! 1) cubed-sphere grid would need to have in order to have the same no-
! minal cell size as the average cell size of the regional grid.
!
!=======================================================================
!
RES_equiv = nint( (2.0*pi_geom*radius_Earth)/(4.0*avg_cell_size) )

WRITE(*,500)
WRITE(*,500) "Equivalent global uniform cubed-sphere resolution is:"
WRITE(*,530) " RES_equiv = ", RES_equiv

500 FORMAT(A)
520 FORMAT(A, I7)
530 FORMAT(A, G)

end subroutine global_equiv_resol

subroutine check(status)
use netcdf
integer,intent(in) :: status
Expand Down
14 changes: 9 additions & 5 deletions ush/fv3gfs_driver_grid.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,10 +315,7 @@ elif [ $gtype = regional ]; then

export ntiles=1
tile=7
rn=$( echo "$stretch_fac * 10" | bc | cut -c1-2 )
name=C${res}r${rn}n${refine_ratio}_${title}
out_dir=$out_dir/C${res}
mkdir -p $out_dir
name=gfdl
grid_dir=$TMPDIR/${name}/grid
orog_dir=$TMPDIR/$name/orog
filter_dir=$orog_dir # nested grid topography will be filtered online
Expand All @@ -336,6 +333,13 @@ elif [ $gtype = regional ]; then
exit $err
fi

# redefine res for gfdl regional grids.

res=$( ncdump -h ${grid_dir}/C*_grid.tile7.nc | grep -o ":RES_equiv = [0-9]\+" | grep -o "[0-9]" )
res=${res//$'\n'/}
out_dir=$out_dir/C${res}
mkdir -p $out_dir

echo "Begin regional orography generation at `date`"

#----------------------------------------------------------------------------------
Expand Down Expand Up @@ -434,7 +438,7 @@ elif [ $gtype = regional_esg ]; then

halop1=$(( halo + 1 ))
tile=7
name=regional
name=regional_esg
grid_dir=$TMPDIR/${name}/grid
orog_dir=$TMPDIR/${name}/orog
filter_dir=$TMPDIR/${name}/filter_topo
Expand Down
27 changes: 19 additions & 8 deletions ush/fv3gfs_make_grid.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,16 +121,21 @@ fi
if [ $gtype = regional_esg ]; then

$APRUN $exec_dir/global_equiv_resol regional_grid.nc
if [ $? -ne 0 ]; then
set +x
echo
echo "FATAL ERROR running global_equiv_resol."
echo
set -x
exit 2
fi

elif [ $gtype = regional ]; then

$APRUN $exec_dir/global_equiv_resol C${res}_grid.tile7.nc

fi

if [ $? -ne 0 ]; then
set +x
echo
echo "FATAL ERROR running global_equiv_resol."
echo
set -x
exit 2
fi

#---------------------------------------------------------------------------------------
# Create mosaic file.
Expand Down Expand Up @@ -169,6 +174,12 @@ elif [ $gtype = nest ]; then

elif [ $gtype = regional ];then

# For gfdl regional grids, redefine the CRES value.

res_save=$res
res=$( ncdump -h C${res}_grid.tile7.nc | grep -o ":RES_equiv = [0-9]\+" | grep -o "[0-9]" )
res=${res//$'\n'/}
mv C${res_save}_grid.tile7.nc C${res}_grid.tile7.nc
$APRUN $executable --num_tiles $ntiles --dir $outdir --mosaic C${res}_mosaic --tile_file C${res}_grid.tile7.nc

elif [ $gtype = regional_esg ]; then
Expand Down

0 comments on commit 8d01219

Please sign in to comment.