get_points_in_spherical_domain Subroutine

 subroutine get_points_in_spherical_domain(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,rc)
    real(real64), intent(in) :: center_lats(:,:),center_lons(:,:)
    real(real64), intent(in) :: corner_lats(:,:),corner_lons(:,:)
    real(real64), intent(in) :: lons(:),lats(:)
    integer, intent(out) :: ii(:),jj(:)
    integer, intent(out), optional :: rc

    integer :: npts,i,n,niter,im,jm,ilb,jlb,iub,jub,ifound,jfound
    integer :: lold,uold,lnew,unew
    logical :: in_region,in_sub_region

    npts = size(lats)

    _ASSERT(size(lats)==size(lons),"lats and lons do not match")
    _ASSERT(npts==size(ii),"size of ii does not match")
    _ASSERT(npts==size(ii),"size of jj does not match")

    im=size(corner_lons,1)-1
    jm=size(corner_lons,2)-1
    niter = max(im,jm)

    do i=1,npts
       ifound=-1
       jfound=-1
       ilb=1
       iub=im
       jlb=1
       jub=jm
       in_region = point_in_polygon([lons(i),lats(i)],[center_lons(ilb,jlb),center_lats(ilb,jlb)],  &
           [corner_lons(ilb,jlb),corner_lats(ilb,jlb)], &
           [corner_lons(iub+1,jlb),corner_lats(iub+1,jlb)], &
           [corner_lons(iub+1,jub+1),corner_lats(iub+1,jub+1)], &
           [corner_lons(ilb,jub+1),corner_lats(ilb,jub+1)])
       if (in_region) then
          ! bisect first dimension
          lnew=ilb
          unew=iub 
          do n = 1,niter
             lold=lnew
             uold=unew
             unew=lold+(uold-lold)/2
             in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(lnew,jlb),center_lats(lnew,jlb)], &
                 [corner_lons(lnew,jlb),corner_lats(lnew,jlb)], &
                 [corner_lons(unew+1,jlb),corner_lats(unew+1,jlb)], &
                 [corner_lons(unew+1,jub+1),corner_lats(unew+1,jub+1)], &
                 [corner_lons(lnew,jub+1),corner_lats(lnew,jub+1)])
             if (in_sub_region) then
               lnew=lold
               unew=unew 
             else
               lnew=unew+1
               unew=uold
             end if
             if (unew==lnew) then
                ifound=unew
                exit
             end if
          enddo
          ! bisect 2nd dimension
          lnew=jlb
          unew=jub
          do n = 1,niter
             lold=lnew
             uold=unew
             unew=lold+(uold-lold)/2
             in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(ifound,lnew),center_lats(ifound,lnew)] , &
                 [corner_lons(ifound,lnew),corner_lats(ifound,lnew)], &
                 [corner_lons(ifound+1,lnew),corner_lats(ifound+1,lnew)], &
                 [corner_lons(ifound+1,unew+1),corner_lats(ifound+1,unew+1)], &
                 [corner_lons(ifound,unew+1),corner_lats(ifound,unew+1)])
             if (in_sub_region) then
               lnew=lold
               unew=unew 
             else
               lnew=unew+1
               unew=uold
             end if
             if (unew==lnew) then
                jfound=unew
                exit
             end if
          enddo
       end if
       ii(i)=ifound
       jj(i)=jfound
    enddo
    _RETURN(_SUCCESS)
         
 end subroutine get_points_in_spherical_domain