subroutine get_gnomonic_ed_coords(im, start, lambda, theta)
integer, intent(in) :: im
integer, intent(in) :: start(:)
real(ESMF_KIND_R8), intent(out) :: lambda(start(1):, start(2):)
real(ESMF_KIND_R8), intent(out) :: theta(start(1):, start(2):)
! Local:
real(ESMF_KIND_R8) :: pp(3,lbound(lambda,1):ubound(lambda,1), lbound(lambda,2):ubound(lambda,2))
real(ESMF_KIND_R8) :: rsq3, alpha, beta, dely
real(ESMF_KIND_R8) :: b(1:im+1)
integer i, j
rsq3 = 1./sqrt(3.)
alpha = asin( rsq3 )
dely = alpha / im
! Compute distances along cube edge (segment) at equispaced
! angles. The formula is same for i and for j, but for local
! region these may or may not overlap. Worst case we
! do this twice for a "diagonal subdomain"
do i = lbound(lambda,1), ubound(lambda,1)
beta = -dely * (im + 2 - 2*i)
b(i) = tan(beta)*cos(alpha)
end do
do j = lbound(lambda,2), ubound(lambda,2)
beta = -dely * (im + 2 - 2*j)
b(j) = tan(beta)*cos(alpha)
end do
! Use the edge position to construct an array of
! cartesian points in the local domain.
do j = lbound(lambda,2), ubound(lambda,2)
do i = lbound(lambda,1), ubound(lambda,1)
pp(:,i,j) = [-rsq3, -b(i), b(j)]
end do
end do
call cart_to_latlon_new(pp, lambda, theta)
end subroutine get_gnomonic_ed_coords