Get dx and dy within SHOC from get_area_p

[PeterCaldwell]

Currently SHOC assumes dx and dy are input variables to SHOC interface. They should instead be calculated inside SHOC's interface. I think there are a couple steps needed:

1. create a function associated with the grid which returns the area of a given grid cell. This will be a C++ version of the existing "get_area_p" function in E3SM.
2. compute dx=dy = sqrt(area)

Note that I was originally advocating for dx=dy being replaced by a single "cell_len" variable since dx and dy are identical... but I realize now that if SHOC is ported to other models, those other models may have different dx and dy. So we should keep dx and dy separate as futureproofing.

[tcclevenger]

@PeterCaldwell @bartgol @AaronDonahue

In shoc_intr.F90 dx/dy are set in a function called grid_size():

subroutine grid_size(state, grid_dx, grid_dy)

do i=1,state%ncol
  ! determine the column area in radians
  column_area = get_area_p(state%lchnk,i)
  ! convert to degrees
  degree = sqrt(column_area)*(180._r8/shr_const_pi)

  ! Now find meters per degree latitude
  ! Below equation finds distance between two points on an ellipsoid, derived from expansion
  !  taking into account ellipsoid using World Geodetic System (WGS84) reference 
  mpdeglat = earth_ellipsoid1 - earth_ellipsoid2 * cos(2._r8*state%lat(i)) + earth_ellipsoid3 * cos(4._r8*state%lat(i))
  grid_dx(i) = mpdeglat * degree
  grid_dy(i) = grid_dx(i) ! Assume these are the same
enddo 

end subroutine grid_size

Will shoc need to perform these calculations? Or just set dx=dy=sqrt(column_area)? If we need them, will state%lat need to be added to the FM?

[PeterCaldwell]

Good question. It appears that get_area_p returns the column area in radians, while we need dx and dy in meters. So while we can still use dist=sqrt(area) - and in fact that's what the above code does in its calculation of "degree" - we need to take extra steps to convert to meters. One question is whether the dycore knows its cell area in square meters instead of just radians. If so, we should just use sqrt of that. If not, we'll probably and up using the code you've quoted above.

One somewhat tangential question in my mind is what the impact of non-sphericalness of the planet has on our solution. I think the surface area of a sphere in square radians is 4pi, so the fraction of the globe covered by a particular cell is its area in radians divided by 4pi. We could therefore probably multiply this ratio by the actual surface area of the Earth to get a good approximation to the cell area in meters without using latitudes?

I'm also a bit confused about how an elliptical earth would affect dx and dy. In particular, it seems like if latitude is stretched by non-sphericalness, then longitude would have gotten squashed (which isn't happening in the code above). In short, it might be more accurate to just ignore the nonsphericalness?

[bartgol]

A quick look to dyn_grid.F90 and Homme suggests that dynamics computes area as squared radians. At least that's what I understand from how the elem%D cube-field contravariant-on-sphere-field map is computed here. But some Homme guru may be able to confirm/refute this (@mt5555, @oksanaguba maybe).

[tcclevenger]

Finally back to looking at this.. I see the get/get_geometry_data() functionality for the grid. Currently the shoc_standalone test uses the PointGrid class which doesn't yet store anything for area or latitude (althought the impl of set_geometry_data() includes the check EKAT_REQUIRE_MSG (name=="lat" || name=="lon" || name=="area") showing it is expecting that data).

My question is what a grid like PointGrid should store? In shoc_intr.F90, there are 3 different ways they define dx/dy

if (single_column .and. .not. iop_mode) then
     host_dx(:) = 100000._r8
     host_dy(:) = 100000._r8
   else if (iop_mode) then
     call grid_size_uniform(host_dx, host_dy)
   else
     call grid_size(state1, host_dx_in, host_dy_in)
   endif

For now, could we just define PointGrid::geo_info["area"] to be some constant and take dx=dy=sqrt(area) (if-clause 1 above)? This would be us storing area in meters.

Else, not sure what a physically relevant non-uniform area for the point grid would look like.

[bartgol]

If we're running coupled with dynamics, then the DynamicsDrivenGridsManager already takes care of setting this data in the grids. The issue is when running a physics standalone test, which uses the PhysicsOnlyGridManager, which creates a dummy physiccs grid with the requested number of columns. For this case, I think we should do a back-of-the-envelope math calculation, and figure out a reasonable value for the area geo data, and set it in the grid (e.g., 4\pi*R^2/ncols). The lat/lon data is a bit more complicated. I would vote for setting views full of NaN in the grid, so that the view exist, and can be retrieved, avoiding crashes if not used, and likely causing crashes later on if one attempts to use them. And if we ever see that standalone tests do require lat/lon views, we can think about how to provide valid data (e.g., random, or constant, or actually compute it...).