MAPL_SphericalGeometry.F90 Source File


This file depends on

sourcefile~~mapl_sphericalgeometry.f90~~EfferentGraph sourcefile~mapl_sphericalgeometry.f90 MAPL_SphericalGeometry.F90 sourcefile~constants.f90 Constants.F90 sourcefile~mapl_sphericalgeometry.f90->sourcefile~constants.f90 sourcefile~mapl_errorhandling.f90 MAPL_ErrorHandling.F90 sourcefile~mapl_sphericalgeometry.f90->sourcefile~mapl_errorhandling.f90 sourcefile~mapl_keywordenforcer.f90 MAPL_KeywordEnforcer.F90 sourcefile~mapl_sphericalgeometry.f90->sourcefile~mapl_keywordenforcer.f90 sourcefile~internalconstants.f90 InternalConstants.F90 sourcefile~constants.f90->sourcefile~internalconstants.f90 sourcefile~mathconstants.f90 MathConstants.F90 sourcefile~constants.f90->sourcefile~mathconstants.f90 sourcefile~physicalconstants.f90 PhysicalConstants.F90 sourcefile~constants.f90->sourcefile~physicalconstants.f90 sourcefile~mapl_throw.f90 MAPL_Throw.F90 sourcefile~mapl_errorhandling.f90->sourcefile~mapl_throw.f90 sourcefile~physicalconstants.f90->sourcefile~mathconstants.f90

Files dependent on this one

sourcefile~~mapl_sphericalgeometry.f90~~AfferentGraph sourcefile~mapl_sphericalgeometry.f90 MAPL_SphericalGeometry.F90 sourcefile~base.f90 Base.F90 sourcefile~base.f90->sourcefile~mapl_sphericalgeometry.f90 sourcefile~base_base_implementation.f90 Base_Base_implementation.F90 sourcefile~base_base_implementation.f90->sourcefile~mapl_sphericalgeometry.f90 sourcefile~mapl.f90 MAPL.F90 sourcefile~mapl.f90->sourcefile~base.f90 sourcefile~mapl_nuopcwrappermod.f90 MAPL_NUOPCWrapperMod.F90 sourcefile~mapl_nuopcwrappermod.f90->sourcefile~base.f90 sourcefile~test_cfio_bundle.pf Test_CFIO_Bundle.pf sourcefile~test_cfio_bundle.pf->sourcefile~base.f90 sourcefile~tstqsat.f90 tstqsat.F90 sourcefile~tstqsat.f90->sourcefile~base.f90 sourcefile~ut_extdata.f90 ut_ExtData.F90 sourcefile~ut_extdata.f90->sourcefile~base.f90 sourcefile~utcfio_bundle.f90 utCFIO_Bundle.F90 sourcefile~utcfio_bundle.f90->sourcefile~base.f90 sourcefile~capdriver.f90 CapDriver.F90 sourcefile~capdriver.f90->sourcefile~mapl.f90 sourcefile~extdataroot_gridcomp.f90 ExtDataRoot_GridComp.F90 sourcefile~capdriver.f90->sourcefile~extdataroot_gridcomp.f90 sourcefile~comp_testing_driver.f90 Comp_Testing_Driver.F90 sourcefile~comp_testing_driver.f90->sourcefile~mapl.f90 sourcefile~extdatadriver.f90 ExtDataDriver.F90 sourcefile~extdatadriver.f90->sourcefile~mapl.f90 sourcefile~extdatadrivergridcomp.f90 ExtDataDriverGridComp.F90 sourcefile~extdatadriver.f90->sourcefile~extdatadrivergridcomp.f90 sourcefile~extdatadrivermod.f90 ExtDataDriverMod.F90 sourcefile~extdatadriver.f90->sourcefile~extdatadrivermod.f90 sourcefile~extdatadriver.f90->sourcefile~extdataroot_gridcomp.f90 sourcefile~extdatadrivergridcomp.f90->sourcefile~mapl.f90 sourcefile~extdatadrivermod.f90->sourcefile~mapl.f90 sourcefile~extdatadrivermod.f90->sourcefile~extdatadrivergridcomp.f90 sourcefile~extdatadrivermod.f90->sourcefile~extdataroot_gridcomp.f90 sourcefile~extdataroot_gridcomp.f90->sourcefile~mapl.f90 sourcefile~varspecdescription.f90 VarspecDescription.F90 sourcefile~extdataroot_gridcomp.f90->sourcefile~varspecdescription.f90 sourcefile~mapl_demo_fargparse.f90 MAPL_demo_fargparse.F90 sourcefile~mapl_demo_fargparse.f90->sourcefile~mapl.f90 sourcefile~pfio_mapl_demo.f90 pfio_MAPL_demo.F90 sourcefile~pfio_mapl_demo.f90->sourcefile~mapl.f90 sourcefile~regrid_util.f90 Regrid_Util.F90 sourcefile~regrid_util.f90->sourcefile~mapl.f90 sourcefile~time_ave_util.f90 time_ave_util.F90 sourcefile~time_ave_util.f90->sourcefile~mapl.f90 sourcefile~varspecdescription.f90->sourcefile~mapl.f90

Source Code

#include "MAPL_ErrLog.h"
module MAPL_SphericalGeometry
   use MAPL_KeywordEnforcerMod
   use mapl_ErrorHandlingMod
   use ESMF
   use MAPL_Constants
   use, intrinsic :: iso_fortran_env, only: REAL64,REAL32

implicit none
private
public get_points_in_spherical_domain
public get_area_spherical_polygon
contains

 ! get area of spherical rectangle given the four corners
 ! p4 ------ p3
 !    |    |
 !    |    |
 !    |    |
 ! p1 ------ p2
 function get_area_spherical_polygon(p1,p4,p2,p3) result(area)
    real(real64) :: area
    real(real64), intent(in) :: p1(2),p2(2),p3(2),p4(2)
    
    real(real64) :: e1(3),e2(3),e3(3)
    real(real64) :: ang1,ang2,ang3,ang4

    e1 = convert_to_cart(p1)
    e2 = convert_to_cart(p2)
    e3 = convert_to_cart(p4)
    ang1 = spherical_angles(e1, e2, e3)

    e1 = convert_to_cart(p2)
    e2 = convert_to_cart(p3)
    e3 = convert_to_cart(p1)
    ang2 = spherical_angles(e1, e2, e3)

    e1 = convert_to_cart(p3)
    e2 = convert_to_cart(p4)
    e3 = convert_to_cart(p2)
    ang3 = spherical_angles(e1, e2, e3)

    e1 = convert_to_cart(p4)
    e2 = convert_to_cart(p3)
    e3 = convert_to_cart(p1)
    ang4 = spherical_angles(e1, e2, e3)

    area = ang1 + ang2 + ang3 + ang4 - 2.0d0*MAPL_PI_R8

 end function get_area_spherical_polygon

 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 

 function point_in_polygon(p0,pinside,a1,a2,a3,a4) result(in_poly)
    real(real64), intent(in) :: p0(2),pinside(2),a1(2),a2(2),a3(2),a4(2)
    logical :: in_poly
 
    real(real64) :: p1c(3),p2c(3),a1c(3),a2c(3),a3c(3),a4c(3)
    logical :: intersect(4)
    p1c=convert_to_cart(p0)
    p2c=convert_to_cart(pinside)
    a1c=convert_to_cart(a1)
    a2c=convert_to_cart(a2)
    a3c=convert_to_cart(a3)
    a4c=convert_to_cart(a4)

    intersect(1) = lines_intersect(p1c,p2c,a1c,a2c)
    intersect(2) = lines_intersect(p1c,p2c,a2c,a3c)
    intersect(3) = lines_intersect(p1c,p2c,a3c,a4c)
    intersect(4) = lines_intersect(p1c,p2c,a4c,a1c)
    if (mod(count(intersect),2)==0) then
       in_poly=.true.
    else
       in_poly=.false.
    end if


 end function point_in_polygon

! it is claimed this should work but doesn't
 !function point_in_polygon_crosprod(p1,a1,a2,a3,a4) result(in_poly)
    !real(real64), intent(in) :: p1(2),a1(2),a2(2),a3(2),a4(2)
    !logical :: in_poly
 
    !real(real64) :: p1c(3),a1c(3),a2c(3),a3c(3),a4c(3)
    !real(real64) :: crs12(3),crs23(3),crs34(3),crs41(3)
    !real(real64) :: d12,d23,d34,d41
    !logical :: signs(4)
    ! a1 -> a2 -> a3 -> a4 so a4 connects to a1

    !p1c=convert_to_cart(p1)
    !a1c=convert_to_cart(a1)
    !a2c=convert_to_cart(a2)
    !a3c=convert_to_cart(a3)
    !a4c=convert_to_cart(a4)

    !crs12 = cross_prod(a1c,a2c)
    !crs23 = cross_prod(a2c,a3c)
    !crs34 = cross_prod(a3c,a4c)
    !crs41 = cross_prod(a4c,a1c)
    !d12=dot_product(p1c,crs12)
    !d23=dot_product(p1c,crs23)
    !d34=dot_product(p1c,crs34)
    !d41=dot_product(p1c,crs41)
    !signs(1)= (d12<0.0)
    !signs(2)= (d23<0.0)
    !signs(3)= (d34<0.0)
    !signs(4)= (d41<0.0)
    !in_poly=( (count(signs)==0) .or. (count(signs)==4) )

 !end function point_in_polygon_crossprod

 function lines_intersect(b0,b1,a0,a1)  result(intersect)
    real(real64), intent(in) :: b0(3),b1(3),a0(3),a1(3)
    logical :: intersect
    real(real64) :: p(3),q(3),t(3)
    real(real64) :: s1,s2,s3,s4
    logical :: signs(4)

    intersect=.false.
    q=cross_prod(b0,b1)
    p=cross_prod(a0,a1)
    t=normal_vect(cross_prod(p,q))

    s1=dot_product(cross_prod(a0,p),t)
    s2=dot_product(cross_prod(a1,p),t)
    s3=dot_product(cross_prod(b0,q),t)
    s4=dot_product(cross_prod(b1,q),t)

    signs(1) = -s1 <0.d0
    signs(2) = s2 <0.d0
    signs(3) = -s3 < 0.d0
    signs(4) = s4 < 0.d0

    intersect = ((count(signs)==0) .or. (count(signs)==4))

 end function lines_intersect

 function normal_vect(vin) result(vout)
    real(real64), intent(in) :: vin(3)
    real(real64) :: vout(3)
    vout=vin/sqrt(vin(1)*vin(1)+vin(2)*vin(2)+vin(3)*vin(3))

 end function normal_vect

 function cross_prod(v1,v2) result(vout)
    real(real64), intent(in) :: v1(3),v2(3)
    real(real64) :: vout(3)
    vout(1)=v1(2)*v2(3)-v1(3)*v2(2)
    vout(2)=v1(3)*v2(1)-v1(1)*v2(3)
    vout(3)=v1(1)*v2(2)-v1(2)*v2(1)
 end function cross_prod

 function convert_to_cart(v) result(xyz)
    real(real64), intent(in) :: v(2)
    real(real64) :: xyz(3)

    xyz(1)=cos(v(2))*cos(v(1))
    xyz(2)=cos(v(2))*sin(v(1))
    xyz(3)=sin(v(2))

 end function convert_to_cart

function vect_mag(v) result(mag)
   real(real64), intent(in) :: v(3)
   real :: mag
   mag = sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
end function vect_mag

function spherical_angles(p1,p2,p3) result(spherical_angle)
   real(real64) :: spherical_angle
   real(real64), intent(in) :: p1(3),p2(3),p3(3)

   real (real64):: e1(3), e2(3), e3(3)
   real (real64):: px, py, pz
   real (real64):: qx, qy, qz
   real (real64):: angle, ddd
   integer n

   do n=1,3
      e1(n) = p1(n)
      e2(n) = p2(n)
      e3(n) = p3(n)
   enddo

   !-------------------------------------------------------------------
   ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry
   !-------------------------------------------------------------------
   ! Vector P:
   px = e1(2)*e2(3) - e1(3)*e2(2)
   py = e1(3)*e2(1) - e1(1)*e2(3)
   pz = e1(1)*e2(2) - e1(2)*e2(1)
   ! Vector Q:
   qx = e1(2)*e3(3) - e1(3)*e3(2)
   qy = e1(3)*e3(1) - e1(1)*e3(3)
   qz = e1(1)*e3(2) - e1(2)*e3(1)

   ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz)

   if ( ddd <= 0.0d0 ) then
      angle = 0.d0
   else
      ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd)
      if ( abs(ddd)>1.d0) then
         angle = 0.5d0 * MAPL_PI_R8
         if (ddd < 0.d0) then
            angle = MAPL_PI_R8
         else
            angle = 0.d0
         end if
      else
         angle = acos( ddd )
      endif
   endif

   spherical_angle = angle
end function

end module MAPL_SphericalGeometry