feature/regional_grid This commit references #4.

Add Gerard's method to compute the global equivalent resolution
for a regional grid to filter_topo.F90.  Added as new routine
"global_equiv_resol".

sorc/fre-nctools.fd/tools/filter_topo/filter_topo.F90

@@ -63,6 +63,7 @@ program filter_topo
 
   !--- compute filter constants for the regional resolution
   if(regional)then
+    call global_equiv_resol
     call compute_filter_constants
   endif
 
@@ -1146,7 +1147,7 @@ subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles
     real::  a1(is-1:ie+2)
     real::  a2(is:ie,js-1:je+2)
     real::  a3(is:ie,js:je,ntiles)
-    real:: smax, smin, m_slope, fac
+    real:: smax, m_slope, fac
     integer:: i,j, nt, t
     integer:: is1, ie2, js1, je2
 
@@ -1763,6 +1764,173 @@ 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
+!
+  if(status /= nf90_noerr) then
+    write(0,*) ' check netcdf status = ', status
+    write(0,'("error ", a)') trim(nf90_strerror(status))
+    write(0,*) "Stopped"
+    stop 4
+  endif
+  end subroutine check
+
   !#######################################################################
   ! compute resolution-dependent values for the filtering.
 
@@ -1808,6 +1976,9 @@ subroutine compute_filter_constants
     enddo
   endif
 
+  print*
+  print*,'global cres equival, toms method ',res_regional
+
   n_del2_weak = nint(n_del2_weak_vals(index1)+factor*(n_del2_weak_vals(index2)-n_del2_weak_vals(index1)))
   cd4 = cd4_vals(index1)+factor*(cd4_vals(index2)-cd4_vals(index1))
   max_slope = max_slope_vals(index1)+factor*(max_slope_vals(index2)-max_slope_vals(index1))