diff --git a/CHANGELOG.md b/CHANGELOG.md
index cd0d451e9f43..7ac210167fde 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,10 +9,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- Corrected bug in HorizontalFluxRegridder.  Fluxes need to be
+  multiplied by edge length for correct treatment.
+
 ### Added
 
 ### Changed
 
+- Changed default decomposition algorithm in Base/base.
+  - More optimal for CS case
+  - Hopefully aligns with common choices for native decomp to reduce need for nontrivial regridding.
+
 ### Removed
 
 ### Deprecated
diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90
index 6d2cb9053397..d43dc716bba5 100644
--- a/base/Base/Base_Base_implementation.F90
+++ b/base/Base/Base_Base_implementation.F90
@@ -832,7 +832,7 @@ module subroutine MAPL_MakeDecomposition(nx, ny, unusable, reduceFactor, rc)
 
     type (ESMF_VM) :: vm
     integer :: pet_count
-
+    integer :: bias
     character(len=*), parameter :: Iam= __FILE__ // '::MAPL_MakeDecomposition()'
     integer :: status
 
@@ -842,11 +842,18 @@ module subroutine MAPL_MakeDecomposition(nx, ny, unusable, reduceFactor, rc)
     _VERIFY(status)
     call ESMF_VMGet(vm, petCount=pet_count, rc=status)
     _VERIFY(status)
-    if (present(reduceFactor)) pet_count=pet_count/reduceFactor
+    if (present(reduceFactor)) then
+       pet_count=pet_count/reduceFactor
+       ! Assume CS
+       bias = 1
+    else
+       ! Assume Lat-Lon
+       bias =2
+    end if
 
     ! count down from sqrt(n)
     ! Note: inal iteration (nx=1) is guaranteed to succeed.
-    do nx = floor(sqrt(real(2*pet_count))), 1, -1
+    do nx = nint(sqrt(real(bias*pet_count))), 1, -1
        if (mod(pet_count, nx) == 0) then ! found a decomposition
           ny = pet_count / nx
           exit
diff --git a/base/HorizontalFluxRegridder.F90 b/base/HorizontalFluxRegridder.F90
index b73b2bebdca4..4a2887b6663f 100644
--- a/base/HorizontalFluxRegridder.F90
+++ b/base/HorizontalFluxRegridder.F90
@@ -10,6 +10,7 @@ module mapl_HorizontalFluxRegridder
    use mapl_KeywordEnforcerMod
    use mapl_ErrorHandlingMod
    use mapl_BaseMod
+   use mapl_SphericalGeometry
    implicit none
    private
 
@@ -20,6 +21,8 @@ module mapl_HorizontalFluxRegridder
       integer :: resolution_ratio = -1
       integer :: im_in, jm_in
       integer :: im_out, jm_out
+      real, allocatable :: dx_in(:,:), dy_in(:,:)
+      real, allocatable :: dx_out(:,:), dy_out(:,:)
    contains
       procedure, nopass :: supports
       procedure :: initialize_subclass
@@ -70,6 +73,8 @@ subroutine initialize_subclass(this, unusable, rc)
 
      integer :: counts(5)
      integer :: status
+     integer :: units ! unused
+     real(kind=ESMF_KIND_R8), pointer :: lons(:,:), lats(:,:)
 
      _UNUSED_DUMMY(unusable)
      spec = this%get_spec()
@@ -91,8 +96,36 @@ subroutine initialize_subclass(this, unusable, rc)
          _ASSERT((IM_in / IM_out) == (JM_in / JM_out), 'inconsistent aspect ratio')
          
          this%resolution_ratio = (IM_in / IM_out)
+
+         call ESMF_GridGetCoord(grid_in, coordDim=1, farrayPtr=lons, &
+              localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__)
+         call ESMF_GridGetCoord(grid_in, coordDim=2, farrayPtr=lats, &
+              localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__)
+
+         this%dx_in = distance( &
+              lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), &
+              lons(2:IM_in+1,1:JM_in), lats(2:IM_in+1,1:JM_in))
+
+         this%dy_in = distance( &
+              lons(1:IM_in,1:JM_in), lats(1:IM_in,1:JM_in), &
+              lons(1:IM_in,2:JM_in+1), lats(1:IM_in,2:JM_in+1))
+
+         call ESMF_GridGetCoord(grid_out, coordDim=1, farrayPtr=lons, &
+              localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__)
+         call ESMF_GridGetCoord(grid_out, coordDim=2, farrayPtr=lats, &
+              localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, __RC__)
+
+         this%dx_out = distance( &
+              lons(1:IM_out,1:JM_out), lats(1:IM_out,1:JM_out), &
+              lons(2:IM_out+1,1:JM_out), lats(2:IM_out+1,1:JM_out))
+
+         this%dy_out = distance( &
+              lons(1:IM_out,1:JM_out), lats(1:IM_out,1:JM_out), &
+              lons(1:IM_out,2:JM_out+1), lats(1:IM_out,2:JM_out+1))
+
        end associate
      end associate
+
      
      _RETURN(_SUCCESS)
   end subroutine initialize_subclass
@@ -129,9 +162,14 @@ subroutine regrid_vector_2d_real32(this, u_in, v_in, u_out, v_out, rotate, rc)
              do i  = 1, IM
                 m_y = 0
                 do ii = 1 + (i-1)*N, i*N
-                   m_y = m_y + v_in(ii,jj)
+                   associate (d_in => this%dx_in(ii,jj))
+                     m_y = m_y + v_in(ii,jj) * d_in
+                   end associate
                 end do
-                v_out(i,j) = m_y
+
+                associate (d_out => this%dx_out(i,j))
+                  v_out(i,j) = m_y / d_out
+                end associate
              end do
           end do
           
@@ -141,9 +179,13 @@ subroutine regrid_vector_2d_real32(this, u_in, v_in, u_out, v_out, rotate, rc)
              do j  = 1, JM
                 m_x = 0
                 do jj = 1 + (j-1)*N, j*N
-                   m_x = m_x + u_in(ii,jj)
+                   associate (d_in => this%dy_in(ii,jj))
+                     m_x = m_x + u_in(ii,jj) * d_in
+                   end associate
                 end do
-                u_out(i,j) = m_x
+                associate (d_out => this%dy_out(i,j))
+                  u_out(i,j) = m_x / d_out
+                end associate
              end do
           end do
 
@@ -186,9 +228,13 @@ subroutine regrid_vector_2d_real64(this, u_in, v_in, u_out, v_out, rotate, rc)
              do i  = 1, IM
                 m_y = 0
                 do ii = 1 + (i-1)*N, i*N
-                   m_y = m_y + v_in(ii,jj)
+                   associate (d_in => this%dx_in(ii,jj))
+                     m_y = m_y + v_in(ii,jj) * d_in
+                   end associate
                 end do
-                v_out(i,j) = m_y
+                associate (d_out => this%dx_out(i,j))
+                  v_out(i,j) = m_y / d_out
+                end associate
              end do
           end do
           
@@ -198,9 +244,13 @@ subroutine regrid_vector_2d_real64(this, u_in, v_in, u_out, v_out, rotate, rc)
              do j  = 1, JM
                 m_x = 0
                 do jj = 1 + (j-1)*N, j*N
-                   m_x = m_x + u_in(ii,jj)
+                   associate (d_in => this%dy_in(ii,jj))
+                     m_x = m_x + u_in(ii,jj) * d_in
+                   end associate
                 end do
-                u_out(i,j) = m_x
+                associate (d_out => this%dy_out(i,j))
+                  u_out(i,j) = m_x / d_out
+                end associate
              end do
           end do
 
diff --git a/base/MAPL_SphericalGeometry.F90 b/base/MAPL_SphericalGeometry.F90
index 8cadc85c3d18..fe799b489faf 100644
--- a/base/MAPL_SphericalGeometry.F90
+++ b/base/MAPL_SphericalGeometry.F90
@@ -1,220 +1,252 @@
 #include "MAPL_ErrLog.h"
-module MAPL_SphericalGeometry
-   use MAPL_KeywordEnforcerMod
+module mapl_SphericalGeometry
+   use mapl_KeywordEnforcerMod
    use mapl_ErrorHandlingMod
    use ESMF
-   use MAPL_Constants
-   use, intrinsic :: iso_fortran_env, only: REAL64,REAL32
+   use mapl_Constants
+   use, intrinsic :: iso_fortran_env, only: REAL64, REAL32
+   implicit none
+   private
+
+   public :: get_points_in_spherical_domain
+   public :: distance
+
+
+   interface distance
+      module procedure distance_r32
+      module procedure distance_r64
+   end interface distance
 
-implicit none
-private
-public get_points_in_spherical_domain
 contains
 
- 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
-
-end module MAPL_SphericalGeometry
+   ! Computes distance between two points (in lat lon as radians),
+   ! and returns distance in radians (unit sphere)
+   ! Using formulae from: https://www.movable-type.co.uk/scripts/latlong.html
+
+   elemental real(kind=REAL64) function distance_r64(lon1, lat1, lon2, lat2) result(d)
+      real(kind=REAL64), intent(in) :: lon1, lat1
+      real(kind=REAL64), intent(in) :: lon2, lat2
+
+      associate(a => sin(lat2-lat1)**2 + cos(lat1)*cos(lat2)*sin((lon2-lon1)/2)**2)
+        d = 2*atan2(sqrt(a), sqrt(1-a))
+      end associate
+
+   end function distance_r64
+
+   elemental real(kind=REAL32) function distance_r32(lon1, lat1, lon2, lat2) result(d)
+      real(kind=REAL32), intent(in) :: lon1, lat1
+      real(kind=REAL32), intent(in) :: lon2, lat2
+
+      associate(a => sin(lat2-lat1)**2 + cos(lat1)*cos(lat2)*sin((lon2-lon1)/2)**2)
+        d = 2*atan2(sqrt(a), sqrt(1-a))
+      end associate
+
+   end function distance_r32
+
+   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
+
+end module mapl_SphericalGeometry