LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - coordinate_systems_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 81 | 266 | 30.5 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 7 | 25 | 28.0 %

              Line data    Source code
           1             : module coordinate_systems_mod
           2             :   use shr_kind_mod,     only: r8=>shr_kind_r8
           3             :   use physconst,        only: pi
           4             :   use cam_abortutils,   only: endrun
           5             :
           6             : ! WARNING:  When using this class be sure that you know if the
           7             : ! cubic coordinates are on the unit cube or the [-\pi/4,\pi/4] cube
           8             : ! and if the spherical longitude is in [0,2\pi] or [-\pi,\pi]
           9             :   implicit none
          10             :   private
          11             :
          12             :   real(kind=r8), public, parameter :: DIST_THRESHOLD= 1.0e-9_r8
          13             :   real(kind=r8), parameter :: one=1.0_r8
          14             :   real(kind=r8), parameter :: two=2.0_r8
          15             :
          16             :   type, public :: cartesian2D_t
          17             :      real(r8) :: x             ! x coordinate
          18             :      real(r8) :: y             ! y coordinate
          19             :   end type cartesian2D_t
          20             :
          21             :   type, public :: cartesian3D_t
          22             :      real(r8) :: x             ! x coordinate
          23             :      real(r8) :: y             ! y coordinate
          24             :      real(r8) :: z             ! z coordinate
          25             :   end type cartesian3D_t
          26             :
          27             :   type, public :: spherical_polar_t
          28             :      real(r8) :: r             ! radius
          29             :      real(r8) :: lon           ! longitude
          30             :      real(r8) :: lat           ! latitude
          31             :   end type spherical_polar_t
          32             :
          33             :
          34             :   interface assignment ( = )
          35             :      module procedure copy_cart2d
          36             :      module procedure copy_spherical_polar
          37             :   end interface
          38             :
          39             :   interface operator( == )
          40             :      module procedure eq_cart2d
          41             :   end interface
          42             :
          43             :   interface distance
          44             :      module procedure distance_cart2D
          45             :      module procedure distance_cart2D_v
          46             :      module procedure distance_cart3D
          47             :      module procedure distance_cart3D_v
          48             :   end interface
          49             :
          50             :   interface change_coordinates
          51             :      module procedure spherical_to_cart_v
          52             :      module procedure spherical_to_cart
          53             :      module procedure cart_to_spherical_v
          54             :      module procedure cart_to_spherical
          55             :      module procedure aray_to_spherical
          56             :   end interface
          57             :
          58             :
          59             :   ! ==========================================
          60             :   ! Public Interfaces
          61             :   ! ==========================================
          62             :
          63             :   public :: sphere_tri_area
          64             :   public :: surfareaxy
          65             :   public :: distance
          66             :   public :: change_coordinates
          67             :   public :: cart2cubedsphere    ! (x,y,z)           -> equal-angle (x,y)
          68             :   public :: cart2cubedsphere_failsafe
          69             :   public :: spherical_to_cart   ! (lat,lon)         ->  (x,y,z)
          70             :   public :: projectpoint        ! equal-angle (x,y) ->  (lat,lon)
          71             :                                 ! should be called cubedsphere2spherical
          72             :   public :: cubedsphere2cart    ! equal-angle (x,y) ->  (x,y,z)
          73             :   public :: sphere2cubedsphere  ! (lat,lon)         ->  equal-angle (x,y)
          74             :   public :: cube_face_number_from_cart
          75             :   public :: cube_face_number_from_sphere
          76             :
          77             : ! CE
          78             :   public :: cart2cubedspherexy  !  (x,y,z)          -> gnomonic (x,y)
          79             :   public :: cart2spherical      !  gnominic (x,y)   -> (lat,lon)
          80             :
          81             :   private :: copy_cart2d
          82             :   private :: copy_spherical_polar
          83             :   private :: eq_cart2d
          84             :   private :: distance_cart2D
          85             :   private :: distance_cart2D_v
          86             :   private :: distance_cart3D
          87             :   private :: distance_cart3D_v
          88             :   private :: spherical_to_cart_v
          89             :   !private :: spherical_to_cart
          90             :   private :: cart_to_spherical_v
          91             :   private :: cart_to_spherical
          92             :   private :: aray_to_spherical
          93             :
          94             : contains
          95             :
          96             :   ! ============================================
          97             :   ! copy_cart2d:
          98             :   !
          99             :   ! Overload assignment operator for cartesian2D_t
         100             :   ! ============================================
         101             :
         102             :   subroutine copy_cart2d(cart2,cart1)
         103             :
         104             :     type(cartesian2D_t), intent(out) :: cart2
         105             :     type(cartesian2D_t), intent(in)  :: cart1
         106             :     cart2%x=cart1%x
         107             :     cart2%y=cart1%y
         108             :   end subroutine copy_cart2d
         109             :
         110             :   ! ============================================
         111             :   ! copy_spherical_polar:
         112             :   !
         113             :   ! Overload assignment operator for spherical_polar_t
         114             :   ! ============================================
         115             :
         116      541344 :   pure subroutine copy_spherical_polar(sph2, sph1)
         117             :
         118             :     type(spherical_polar_t), intent(out) :: sph2
         119             :     type(spherical_polar_t), intent(in)  :: sph1
         120      541344 :     sph2%r   = sph1%r
         121      541344 :     sph2%lat = sph1%lat
         122      541344 :     sph2%lon = sph1%lon
         123      541344 :   end subroutine copy_spherical_polar
         124             :
         125             :   ! ============================================
         126             :   ! eq_cart2d:
         127             :   !
         128             :   ! Overload == operator for cartesian2D_t
         129             :   ! ============================================
         130             :
         131             :   pure function eq_cart2d(cart2,cart1) result(is_same)
         132             :
         133             :     type(cartesian2D_t), intent(in)  :: cart2
         134             :     type(cartesian2D_t), intent(in)  :: cart1
         135             :
         136             :     logical :: is_same
         137             :
         138             :     if (distance(cart1,cart2)<DIST_THRESHOLD) then
         139             :        is_same=.true.
         140             :     else
         141             :        is_same=.false.
         142             :     end if
         143             :
         144             :   end function eq_cart2d
         145             :
         146             :   ! ===================================================
         147             :   ! distance_cart2D  : scalar version
         148             :   ! distance_cart2D_v: vector version
         149             :   !
         150             :   ! computes distance between cartesian 2D coordinates
         151             :   ! ===================================================
         152             :
         153           0 :   pure function distance_cart2D(cart1,cart2) result(dist)
         154             :     implicit none
         155             :     type(cartesian2D_t), intent(in)           :: cart1
         156             :     type(cartesian2D_t), intent(in), optional :: cart2
         157             :     real(r8)   :: dist
         158             :
         159           0 :     if (present(cart2)) then
         160             :        dist = SQRT((cart1%x-cart2%x)**2 + &
         161           0 :                    (cart1%y-cart2%y)**2   )
         162             :     else
         163             :        dist = SQRT(cart1%x*cart1%x + &
         164           0 :                    cart1%y*cart1%y   )
         165             :     end if
         166             :
         167           0 :   end function distance_cart2D
         168             :
         169           0 :   pure function distance_cart2D_v(cart1,cart2) result(dist)
         170             :     implicit none
         171             :     type(cartesian2D_t), intent(in)           :: cart1(:)
         172             :     type(cartesian2D_t), intent(in), optional :: cart2(:)
         173             :     real(r8)                           :: dist(SIZE(cart1))
         174             :
         175             :     integer             :: i
         176             :
         177           0 :     if (present(cart2)) then
         178           0 :        forall (i=1:SIZE(cart1)) dist(i) = distance_cart2D(cart1(i),cart2(i))
         179             :     else
         180           0 :        forall (i=1:SIZE(cart1)) dist(i) = distance_cart2D(cart1(i))
         181             :     end if
         182             :
         183           0 :   end function distance_cart2D_v
         184             :
         185             :
         186             :   ! ===================================================
         187             :   ! distance_cart3D  : scalar version
         188             :   ! distance_cart3D_v: vector version
         189             :   ! ===================================================
         190             :
         191           0 :   pure function distance_cart3D(cart1,cart2) result(dist)
         192             :     implicit none
         193             :     type(cartesian3D_t), intent(in)          :: cart1
         194             :     type(cartesian3D_t), intent(in),optional :: cart2
         195             :     real(r8)                          :: dist
         196             :
         197           0 :     if (present(cart2)) then
         198             :        dist = SQRT((cart1%x-cart2%x)**2 + &
         199             :                    (cart1%y-cart2%y)**2 + &
         200           0 :                    (cart1%z-cart2%z)**2   )
         201             :     else
         202             :        dist = SQRT(cart1%x*cart1%x + &
         203             :                    cart1%y*cart1%y + &
         204           0 :                    cart1%z*cart1%z   )
         205             :     end if
         206           0 :   end function distance_cart3D
         207             :
         208           0 :   pure function distance_cart3D_v(cart1,cart2) result(dist)
         209             :     implicit none
         210             :     type(cartesian3D_t), intent(in)          :: cart1(:)
         211             :     type(cartesian3D_t), intent(in),optional :: cart2(:)
         212             :     real(r8)                          :: dist(SIZE(cart1))
         213             :
         214             :     integer             :: i
         215             :
         216           0 :     if (present(cart2)) then
         217           0 :        forall (i=1:SIZE(cart1)) dist(i) = distance_cart3D(cart1(i),cart2(i))
         218             :     else
         219           0 :        forall (i=1:SIZE(cart1)) dist(i) = distance_cart3D(cart1(i))
         220             :     end if
         221           0 :   end function distance_cart3D_v
         222             :
         223             :   ! ===================================================================
         224             :   ! spherical_to_cart:
         225             :   ! converts spherical polar {lon,lat}  to 3D cartesian {x,y,z}
         226             :   ! on unit sphere.  Note: spherical longitude is [0,2\pi]
         227             :   ! ===================================================================
         228             :
         229      195744 :   pure function spherical_to_cart(sphere) result (cart)
         230             :     implicit none
         231             :     type(spherical_polar_t), intent(in) :: sphere
         232             :     type(cartesian3D_t)                 :: cart
         233             :
         234      195744 :     cart%x=sphere%r*COS(sphere%lat)*COS(sphere%lon)
         235      195744 :     cart%y=sphere%r*COS(sphere%lat)*SIN(sphere%lon)
         236      195744 :     cart%z=sphere%r*SIN(sphere%lat)
         237             :
         238      195744 :   end function spherical_to_cart
         239             :
         240             :   ! ===================================================================
         241             :   ! spherical_to_cart_v:
         242             :   ! converts spherical polar {lon,lat}  to 3D cartesian {x,y,z}
         243             :   ! on unit sphere.  Note: spherical longitude is [0,2\pi]
         244             :   ! ===================================================================
         245             :
         246           0 :   pure function spherical_to_cart_v(sphere) result (cart)
         247             :     implicit none
         248             :     type(spherical_polar_t), intent(in) :: sphere(:)
         249             :     type(cartesian3D_t)                 :: cart(SIZE(sphere))
         250             :
         251             :     integer                 :: i
         252             :
         253           0 :     forall (i=1:SIZE(sphere)) cart(i) = spherical_to_cart(sphere(i))
         254           0 :   end function spherical_to_cart_v
         255             :
         256             :   ! ==========================================================================
         257             :   ! cart_to_spherical:
         258             :   !
         259             :   ! converts 3D cartesian {x,y,z} to spherical polar {lon,lat}
         260             :   ! on unit sphere. Note: spherical longitude is [0,2\pi]
         261             :   ! ==========================================================================
         262             :
         263             :   ! scalar version
         264             :
         265           0 :   pure function cart_to_spherical(cart) result (sphere)
         266             :
         267             :     type(cartesian3D_t), intent(in) :: cart
         268             :     type(spherical_polar_t)         :: sphere
         269             :
         270           0 :     sphere%r=distance(cart)
         271           0 :     sphere%lat=ASIN(cart%z/sphere%r)
         272           0 :     sphere%lon=0
         273             :
         274             :     ! ==========================================================
         275             :     ! enforce three facts:
         276             :     !
         277             :     ! 1) lon at poles is defined to be zero
         278             :     !
         279             :     ! 2) Grid points must be separated by about .01 Meter (on earth)
         280             :     !    from pole to be considered "not the pole".
         281             :     !
         282             :     ! 3) range of lon is { 0<= lon < 2*pi }
         283             :     !
         284             :     ! ==========================================================
         285             :
         286             : !   if point is away from the POLE.  distance(cart) = distance from center of earth,
         287             : !   so this was a bug:
         288             : !    if (distance(cart) >= DIST_THRESHOLD) then
         289             :
         290           0 :     if ( abs(abs(sphere%lat)-PI/2)  >= DIST_THRESHOLD ) then
         291           0 :        sphere%lon=ATAN2(cart%y,cart%x)
         292           0 :        if (sphere%lon<0) then
         293           0 :           sphere%lon=sphere%lon + 2*PI
         294             :        end if
         295             :     end if
         296             :
         297           0 :   end function cart_to_spherical
         298             :
         299           0 :   pure function aray_to_spherical(coordinates) result (sphere)
         300             :     implicit none
         301             :     real(kind=r8),  intent(in)    :: coordinates(3)
         302             :     type(spherical_polar_t)              :: sphere
         303             :     type(cartesian3D_t)                  :: cart
         304           0 :     cart%x = coordinates(1)
         305           0 :     cart%y = coordinates(2)
         306           0 :     cart%z = coordinates(3)
         307           0 :     sphere = cart_to_spherical(cart)
         308           0 :   end function aray_to_spherical
         309             :
         310             :
         311           0 :   pure function cart_to_spherical_v(cart) result (sphere)
         312             :
         313             :     type(cartesian3D_t), intent(in) :: cart(:)
         314             :     type(spherical_polar_t)         :: sphere(SIZE(cart))
         315             :
         316             :     integer                 :: i
         317           0 :     forall (i=1:SIZE(cart)) sphere(i) = cart_to_spherical(cart(i))
         318           0 :   end function cart_to_spherical_v
         319             :
         320             :
         321             :
         322             :
         323      541344 :   function unit_face_based_cube_to_unit_sphere(cart, face_no) result(sphere)
         324             :
         325             : ! Note: Output spherical longitude is [-pi,pi]
         326             :
         327             : ! Project from a UNIT cube to a UNIT sphere.  ie, the lenght of the cube edge is 2.
         328             : ! Face 1 of the cube touches the sphere at longitude, latitude (0,0). The negative
         329             : ! x axis is negative longitude (ie. going west is negative), the positive x axis
         330             : ! is increasing longitude.  Face 1 maps the Face 1 to the lat,lon on the sphere:
         331             : !    [-1,1] x [-1,1] => [-\pi/4,\pi/4] x [-\pi/4, \pi/4]
         332             :
         333             : ! Face 2 continues with increasing longitude (ie to the east of Face 1).
         334             : ! The left edge of Face 2 (negative x) is the right edge of Face 1 (positive x)
         335             : ! The latitude is the same as Face 1, but the longitude increases:
         336             : !    [-1,1] x [-1,1] => [\pi/4, 3\pi/4] x [-\pi/4, \pi/4]
         337             :
         338             : ! Face 3 continues with increasing longitude (ie to the east of Face 2).
         339             : ! Face 3 is like Face 1, but the x coordinates are reversed, ie. decreasing x
         340             : ! is increasing longitude:
         341             : !    [-1,1] x [-1,1]  =    [-1,0] x [-1,1] U  [0,1] x [-1,1] =>
         342             : !            [3\pi/4,\pi] x [-\pi, -3\pi/4]
         343             :
         344             : ! Face 4 finally connects Face 3 to Face 1.  Like Face 2, but wtih opposite x
         345             : !    [-1,1] x [-1,1] => [-3\pi/4, -\pi/4] x [-\pi/4, \pi/4]
         346             :
         347             : ! Face 5 is along the bottom edges of Faces 1,2,3,and 4 so the latitude goes from
         348             : ! -\pi/4 to -\pi/2.  The tricky part is lining up the longitude.  The zero longitude
         349             : ! must line up with the center of Face 1. ATAN2(x,1) = 0 => x = 0.
         350             : ! So the (0,1) point on Face 5 is the zero longitude on the sphere.  The top edge of
         351             : ! Face 5 is the bottom edge of Face 1.
         352             : ! ATAN(x,0) = \pi/2 => x = 1, so the right edge of Face 5 is the bottom of Face 2.
         353             : ! Continueing, the bottom edge of 5 is the bottom of 3.  Left of 5 is bottom of 4.
         354             :
         355             : ! Face 6 is along the top edges of Faces 1,2,3 and 4 so the latitude goes from
         356             : ! \pi/4 to \pi/2.   The zero longitude must line up with the center of Face 1.
         357             : ! This is just like Face 5, but the y axis is reversed.  So the bottom edge of Face 6
         358             : ! is the top edge of Face 1.  The right edge of Face 6 is the top of Face 2.  The
         359             : ! top of 6 the top of 3 and the left of 6 the top of 4.
         360             :
         361             :     type (cartesian2d_t), intent(in)     :: cart   ! On face_no of a unit cube
         362             :     integer,              intent(in)     :: face_no
         363             :
         364             :     type (spherical_polar_t)             :: sphere
         365             :
         366             :     integer i,j
         367             :     real(kind=r8) :: r!, l_inf
         368             :
         369             : ! MNL: removing check that points are on the unit cube because we allow
         370             : ! spherical grids to map beyond the extent of the cube (though we probably
         371             : ! should still have an upper bound for how far past the edge the element lies)
         372             : !    l_inf = MAX(ABS(cart%x), ABS(cart%y))
         373             : !    if (1.01 < l_inf) then
         374             : !      call endrun('unit_face_based_cube_to_unit_sphere: Input not on unit cube.')
         375             : !    end if
         376             :
         377      541344 :     sphere%r=one
         378      541344 :     r = SQRT( one + (cart%x)**2 + (cart%y)**2)
         379      659328 :     select case (face_no)
         380             :     case (1)
         381      117984 :        sphere%lat=ASIN((cart%y)/r)
         382      117984 :        sphere%lon=ATAN2(cart%x,one)
         383             :     case (2)
         384       84672 :        sphere%lat=ASIN((cart%y)/r)
         385       84672 :        sphere%lon=ATAN2(one,-cart%x)
         386             :     case (3)
         387       84672 :        sphere%lat=ASIN((cart%y)/r)
         388       84672 :        sphere%lon=ATAN2(-cart%x,-one)
         389             :     case (4)
         390       84672 :        sphere%lat=ASIN((cart%y)/r)
         391       84672 :        sphere%lon=ATAN2(-one,cart%x)
         392             :     case (5)
         393       84672 :        if (ABS(cart%y) > DIST_THRESHOLD .or. ABS(cart%x) > DIST_THRESHOLD ) then
         394       84640 :           sphere%lon=ATAN2(cart%x, cart%y )
         395             :        else
         396          32 :           sphere%lon= 0.0_r8     ! longitude is meaningless at south pole set to 0.0
         397             :        end if
         398       84672 :        sphere%lat=ASIN(-one/r)
         399             :     case (6)
         400       84672 :        if (ABS(cart%y) > DIST_THRESHOLD .or. ABS(cart%x) > DIST_THRESHOLD ) then
         401       84640 :           sphere%lon = ATAN2(cart%x, -cart%y)
         402             :        else
         403          32 :           sphere%lon= 0.0_r8     ! longitude is meaningless at north pole set to 0.0
         404             :        end if
         405       84672 :        sphere%lat=ASIN(one/r)
         406             :     case default
         407      541344 :        call endrun('unit_face_based_cube_to_unit_sphere: Face number not 1 to 6.')
         408             :     end select
         409             :
         410      541344 :     if (sphere%lon < 0.0_r8) then
         411      270640 :        sphere%lon=sphere%lon + two*PI
         412             :     end if
         413             :
         414      541344 :   end function unit_face_based_cube_to_unit_sphere
         415             :
         416      194400 :   function cart2spherical(x,y, face_no) result(sphere)
         417             : ! IMPORTANT: INPUT ARE the REAL cartesian from the cube sphere
         418             : ! Note: Output spherical longitude is [-pi,pi]
         419             :
         420             : ! Project from a UNIT cube to a UNIT sphere.  ie, the lenght of the cube edge is 2.
         421             : ! Face 1 of the cube touches the sphere at longitude, latitude (0,0). The negative
         422             : ! x axis is negative longitude (ie. going west is negative), the positive x axis
         423             : ! is increasing longitude.  Face 1 maps the Face 1 to the lat,lon on the sphere:
         424             : !    [-1,1] x [-1,1] => [-\pi/4,\pi/4] x [-\pi/4, \pi/4]
         425             :
         426             : ! Face 2 continues with increasing longitude (ie to the east of Face 1).
         427             : ! The left edge of Face 2 (negative x) is the right edge of Face 1 (positive x)
         428             : ! The latitude is the same as Face 1, but the longitude increases:
         429             : !    [-1,1] x [-1,1] => [\pi/4, 3\pi/4] x [-\pi/4, \pi/4]
         430             :
         431             : ! Face 3 continues with increasing longitude (ie to the east of Face 2).
         432             : ! Face 3 is like Face 1, but the x coordinates are reversed, ie. decreasing x
         433             : ! is increasing longitude:
         434             : !    [-1,1] x [-1,1]  =    [-1,0] x [-1,1] U  [0,1] x [-1,1] =>
         435             : !            [3\pi/4,\pi] x [-\pi, -3\pi/4]
         436             :
         437             : ! Face 4 finally connects Face 3 to Face 1.  Like Face 2, but wtih opposite x
         438             : !    [-1,1] x [-1,1] => [-3\pi/4, -\pi/4] x [-\pi/4, \pi/4]
         439             :
         440             : ! Face 5 is along the bottom edges of Faces 1,2,3,and 4 so the latitude goes from
         441             : ! -\pi/4 to -\pi/2.  The tricky part is lining up the longitude.  The zero longitude
         442             : ! must line up with the center of Face 1. ATAN2(x,1) = 0 => x = 0.
         443             : ! So the (0,1) point on Face 5 is the zero longitude on the sphere.  The top edge of
         444             : ! Face 5 is the bottom edge of Face 1.
         445             : ! ATAN(x,0) = \pi/2 => x = 1, so the right edge of Face 5 is the bottom of Face 2.
         446             : ! Continueing, the bottom edge of 5 is the bottom of 3.  Left of 5 is bottom of 4.
         447             :
         448             : ! Face 6 is along the top edges of Faces 1,2,3 and 4 so the latitude goes from
         449             : ! \pi/4 to \pi/2.   The zero longitude must line up with the center of Face 1.
         450             : ! This is just like Face 5, but the y axis is reversed.  So the bottom edge of Face 6
         451             : ! is the top edge of Face 1.  The right edge of Face 6 is the top of Face 2.  The
         452             : ! top of 6 the top of 3 and the left of 6 the top of 4.
         453             :
         454             :     implicit none
         455             :     real(kind=r8), intent(in)     :: x,y   ! On face_no of a unit cube
         456             :     integer,              intent(in)     :: face_no
         457             :
         458             :     type (spherical_polar_t)             :: sphere
         459             :
         460             :     integer i,j
         461             :     real(kind=r8) :: r!, l_inf
         462             :
         463             : ! MNL: removing check that points are on the unit cube because we allow
         464             : ! spherical grids to map beyond the extent of the cube (though we probably
         465             : ! should still have an upper bound for how far past the edge the element lies)
         466             : !    l_inf = MAX(ABS(cart%x), ABS(cart%y))
         467             : !    if (1.01 < l_inf) then
         468             : !      call endrun('unit_face_based_cube_to_unit_sphere: Input not on unit cube.')
         469             : !    end if
         470             :
         471      194400 :     sphere%r=one
         472      194400 :     r = SQRT( one + x**2 + y**2)
         473      226800 :     select case (face_no)
         474             :     case (1)
         475       32400 :        sphere%lat=ASIN(y/r)
         476       32400 :        sphere%lon=ATAN2(x,one)
         477             :     case (2)
         478       32400 :        sphere%lat=ASIN(y/r)
         479       32400 :        sphere%lon=ATAN2(one,-x)
         480             :     case (3)
         481       32400 :        sphere%lat=ASIN(y/r)
         482       32400 :        sphere%lon=ATAN2(-x,-one)
         483             :     case (4)
         484       32400 :        sphere%lat=ASIN(y/r)
         485       32400 :        sphere%lon=ATAN2(-one,x)
         486             :     case (5)
         487       32400 :        if (ABS(y) > DIST_THRESHOLD .or. ABS(x) > DIST_THRESHOLD ) then
         488       32400 :           sphere%lon=ATAN2(x, y )
         489             :        else
         490           0 :           sphere%lon= 0.0_r8     ! longitude is meaningless at south pole set to 0.0
         491             :        end if
         492       32400 :        sphere%lat=ASIN(-one/r)
         493             :     case (6)
         494       32400 :        if (ABS(y) > DIST_THRESHOLD .or. ABS(x) > DIST_THRESHOLD ) then
         495       32400 :           sphere%lon = ATAN2(x, -y)
         496             :        else
         497           0 :           sphere%lon= 0.0_r8     ! longitude is meaningless at north pole set to 0.0
         498             :        end if
         499       32400 :        sphere%lat=ASIN(one/r)
         500             :     case default
         501      194400 :        call endrun('unit_face_based_cube_to_unit_sphere: Face number not 1 to 6.')
         502             :     end select
         503             :
         504      194400 :     if (sphere%lon < 0.0_r8) then
         505       97200 :        sphere%lon=sphere%lon + two*PI
         506             :     end if
         507             :
         508      194400 :   end function cart2spherical
         509             :
         510             :
         511             :
         512             :
         513             :
         514             :
         515             :
         516             :
         517             : ! Note: Output spherical longitude is [-pi,pi]
         518      541344 :   function projectpoint(cartin, face_no) result(sphere)
         519             :
         520             : ! Projection from a [-pi/4, \pi/4] sized cube.
         521             : ! This will be checked because unit_face_based_cube_to_unit_sphere checks the ranges.
         522             : ! See unit_face_based_cube_to_unit_sphere for documentation.
         523             :
         524             :     implicit none
         525             :     type (cartesian2d_t), intent(in)     :: cartin
         526             :     integer,              intent(in)     :: face_no
         527             :     type (spherical_polar_t)             :: sphere
         528             :     type (cartesian2d_t)                 :: cart
         529             :
         530             :     !ASC  This is X and Y and not xhi eta ...
         531             :
         532      541344 :     cart%x = TAN(cartin%x)
         533      541344 :     cart%y = TAN(cartin%y)
         534             :
         535      541344 :     sphere = unit_face_based_cube_to_unit_sphere(cart, face_no)
         536             :
         537      541344 :   end function projectpoint
         538             :
         539             :   ! takes a 2D point on a face of the cube of size [-\pi/4, \pi/4] and projects it
         540             :   ! onto a 3D point on a cube of size [-1,1] in R^3
         541      195744 :   function cubedsphere2cart(cartin, face_no) result(cart)
         542             :     implicit none
         543             :     type (cartesian2d_t), intent(in)    :: cartin   ! assumed to be cartesian coordinates of cube
         544             :     integer,              intent(in)    :: face_no
         545             :
         546             :     type(cartesian3D_t)                 :: cart
         547             :
         548      195744 :     cart = spherical_to_cart(projectpoint(cartin, face_no))
         549             :
         550      195744 :   end function cubedsphere2cart
         551             :
         552             :
         553             :   ! onto a cube of size [-\pi/2,\pi/2] in R^3
         554             :   ! the spherical longitude can be either in [0,2\pi] or [-\pi,\pi]
         555           0 :   pure function sphere2cubedsphere (sphere, face_no) result(cart)
         556             :     implicit none
         557             :     type(spherical_polar_t), intent(in) :: sphere
         558             :     integer,                 intent(in) :: face_no
         559             :
         560             :     type(cartesian2d_t)                 :: cart
         561             :     real(kind=r8)               :: xp,yp
         562             :     real(kind=r8)               :: lat,lon
         563             :     real(kind=r8)               :: twopi, pi2, pi3, pi4
         564             :
         565           0 :     lat = sphere%lat
         566           0 :     lon = sphere%lon
         567             :
         568           0 :     twopi = 2.0_r8 * pi
         569           0 :     pi2   = pi * 0.5_r8
         570           0 :     pi3   = pi * 1.5_r8
         571           0 :     pi4   = pi * 0.25_r8
         572             :
         573           0 :     select case (face_no)
         574             :     case  (1)
         575           0 :        xp = lon
         576           0 :        if (pi < lon) xp = lon - twopi !if lon in [0,2\pi]
         577           0 :        yp = atan(tan(lat)/cos(xp))
         578             :     case  (2)
         579           0 :        xp = lon - pi2
         580           0 :        yp = atan(tan(lat)/cos(xp))
         581             :     case  (3)
         582           0 :        xp = lon - pi
         583           0 :        if (lon < 0) xp = lon + pi  !if lon in [0,2\pi]
         584           0 :        yp = atan(tan(lat)/cos(xp))
         585             :     case  (4)
         586           0 :        xp = lon - pi3
         587           0 :        if (lon < 0) xp = lon + pi2  !if lon in [0,2\pi]
         588           0 :        yp = atan(tan(lat)/cos(xp))
         589             :     case  (5)
         590           0 :        xp = atan(-sin(lon)/tan(lat))
         591           0 :        yp = atan(-cos(lon)/tan(lat))
         592             :     case  (6)
         593           0 :        xp = atan( sin(lon)/tan(lat))
         594           0 :        yp = atan(-cos(lon)/tan(lat))
         595             :     end select
         596             :
         597             :     ! coordinates on the cube:
         598           0 :     cart%x = xp
         599           0 :     cart%y = yp
         600             :
         601           0 :   end function sphere2cubedsphere
         602             :
         603             : ! Go from an arbitrary sized cube in 3D
         604             : ! to a [-\pi/4,\pi/4] sized cube with (face,2d) coordinates.
         605             : !
         606             : !                        Z
         607             : !                        |
         608             : !                        |
         609             : !                        |
         610             : !                        |
         611             : !                        ---------------Y
         612             : !                       /
         613             : !                      /
         614             : !                     /
         615             : !                    /
         616             : !                   X
         617             : !
         618             : ! NOTE: Face 1 =>  X positive constant face of cube
         619             : !       Face 2 =>  Y positive constant face of cube
         620             : !       Face 3 =>  X negative constant face of cube
         621             : !       Face 4 =>  Y negative constant face of cube
         622             : !       Face 5 =>  Z negative constant face of cube
         623             : !       Face 6 =>  Z positive constant face of cube
         624      109344 :   pure function cart2cubedsphere(cart3D, face_no) result(cart)
         625             :
         626             :     implicit none
         627             :     type(cartesian3D_t),intent(in) :: cart3d
         628             :     integer,            intent(in) :: face_no
         629             :     type (cartesian2d_t)           :: cart
         630             :
         631             :     real(kind=r8) :: x,y
         632             :
         633      122016 :     select case (face_no)
         634             :     case (1)
         635       12672 :        x =  cart3D%y/cart3D%x
         636       12672 :        y =  cart3D%z/cart3D%x
         637             :     case (2)
         638       21264 :        x = -cart3D%x/cart3D%y
         639       21264 :        y =  cart3D%z/cart3D%y
         640             :     case (3)
         641       12672 :        x =  cart3D%y/cart3D%x
         642       12672 :        y = -cart3D%z/cart3D%x
         643             :     case (4)
         644       21264 :        x = -cart3D%x/cart3D%y
         645       21264 :        y = -cart3D%z/cart3D%y
         646             :     case (5)
         647       20736 :        x  = -cart3D%y/cart3D%z
         648       20736 :        y  = -cart3D%x/cart3D%z
         649             :     case (6)
         650       20736 :        x  =  cart3D%y/cart3D%z
         651      109344 :        y  = -cart3D%x/cart3D%z
         652             :     end select
         653      109344 :     cart%x = ATAN(x)
         654      109344 :     cart%y = ATAN(y)
         655      109344 :   end function cart2cubedsphere
         656             :
         657           0 :   function cart2cubedsphere_failsafe(cart3D, face_no) result(cart)
         658             :     implicit none
         659             :     type(cartesian3D_t),intent(in) :: cart3d
         660             :     integer,            intent(in) :: face_no
         661             :     type (cartesian2d_t)           :: cart
         662             :
         663             :     real(kind=r8) :: x,y
         664             :
         665           0 :     select case (face_no)
         666             :     case (1)
         667           0 :       if (abs(cart3D%x) < 1.E-13_r8) then
         668           0 :         cart%x=9.0E9_r8
         669           0 :         cart%y=9.0E9_r8
         670           0 :         return
         671             :       end if
         672           0 :        x =  cart3D%y/cart3D%x
         673           0 :        y =  cart3D%z/cart3D%x
         674             :     case (2)
         675           0 :       if (abs(cart3D%y)<1.0E-13_r8) then
         676           0 :         cart%x=9.0E9_r8
         677           0 :         cart%y=9.0E9_r8
         678           0 :         return
         679             :       end if
         680           0 :        x = -cart3D%x/cart3D%y
         681           0 :        y =  cart3D%z/cart3D%y
         682             :     case (3)
         683           0 :       if (abs(cart3D%x)<1.0E-13_r8) then
         684           0 :         cart%x=9.0E9_r8
         685           0 :         cart%y=9.0E9_r8
         686           0 :         return
         687             :       end if
         688           0 :        x =  cart3D%y/cart3D%x
         689           0 :        y = -cart3D%z/cart3D%x
         690             :     case (4)
         691           0 :       if (abs(cart3D%y)<1.0E-13_r8) then
         692           0 :         cart%x=9.0E9_r8
         693           0 :         cart%y=9.0E9_r8
         694           0 :         return
         695             :       end if
         696           0 :        x = -cart3D%x/cart3D%y
         697           0 :        y = -cart3D%z/cart3D%y
         698             :     case (5)
         699           0 :       if (abs(cart3D%z)<1.0E-13_r8) then
         700           0 :         cart%x=9.0E9_r8
         701           0 :         cart%y=9.0E9_r8
         702           0 :         return
         703             :       end if
         704           0 :        x  = -cart3D%y/cart3D%z
         705           0 :        y  = -cart3D%x/cart3D%z
         706             :     case (6)
         707           0 :       if (abs(cart3D%z)<1.0E-13_r8) then
         708           0 :         cart%x=9.0E9_r8
         709           0 :         cart%y=9.0E9_r8
         710           0 :         return
         711             :       end if
         712           0 :        x  =  cart3D%y/cart3D%z
         713           0 :        y  = -cart3D%x/cart3D%z
         714             :      case default
         715           0 :        write(*,*) "face_no out out range ",face_no
         716             :     end select
         717           0 :     cart%x = ATAN(x)
         718           0 :     cart%y = ATAN(y)
         719           0 :   end function cart2cubedsphere_failsafe
         720             :
         721             :
         722             :
         723             : ! This function divides three dimentional space up into
         724             : ! six sectors.  These sectors are then considered as the
         725             : ! faces of the cube.  It should work for any (x,y,z) coordinate
         726             : ! if on a sphere or on a cube.
         727           0 :   pure function cube_face_number_from_cart(cart) result(face_no)
         728             :
         729             :     implicit none
         730             :     type(cartesian3D_t),intent(in) :: cart
         731             :     integer :: face_no
         732             :
         733             :     real(r8) :: x,y,z
         734           0 :     x=cart%x
         735           0 :     y=cart%y
         736           0 :     z=cart%z
         737             :
         738             : ! Divide the X-Y plane into for quadrants of
         739             : ! [-\pi/2,\pi/2], [\pi/2,3\pi/2], .....
         740             : ! based on the lines X=Y and X=-Y.  This divides
         741             : ! 3D space up into four sections.  Doing the same
         742             : ! for the XZ and YZ planes divides space into six
         743             : ! sections.  Can also be thought of as conic sections
         744             : ! in the L_infinity norm.
         745             :
         746           0 :     if (y<x .and. y>-x) then      ! x>0, Face 1,5 or 6
         747           0 :       if (z>x) then
         748             :          face_no=6  ! north pole
         749           0 :       else if (z<-x) then
         750             :          face_no=5 ! south pole
         751             :       else
         752           0 :          face_no = 1
         753             :       endif
         754           0 :    else if (y>x .and. y<-x) then  ! x<0
         755           0 :       if (z>-x) then
         756             :          face_no=6 ! north pole
         757           0 :       else if (z<x) then
         758             :          face_no=5 ! south pole
         759             :       else
         760           0 :          face_no=3
         761             :       endif
         762           0 :    else if (y>x .and. y>-x) then  ! y>0
         763           0 :       if (z>y) then
         764             :          face_no=6 ! north pole
         765           0 :       else if (z<-y) then
         766             :         face_no = 5 ! south pole
         767             :       else
         768           0 :          face_no=2
         769             :       endif
         770           0 :    else if (y<x .and. y<-x) then  ! y<0
         771           0 :       if (z>-y) then
         772             :          face_no=6 ! north pole
         773           0 :       else if (z<y) then
         774             :          face_no=5 ! south pole
         775             :       else
         776           0 :          face_no=4
         777             :       endif
         778             :    else
         779             :       ! abs(y) = abs(x).  point is on cube edge, or on face 5 or 6:
         780           0 :       if (abs(x)<z) then
         781             :          face_no=6  ! north pole
         782           0 :       else if (z<-abs(x)) then
         783             :          face_no=5 ! south pole
         784           0 :       else if (0 < x .and. 0 < y) then
         785             :          face_no = 1
         786           0 :       else if (x < 0 .and. 0 < y) then
         787             :          face_no = 2
         788           0 :       else if (x < 0 .and. y < 0) then
         789             :          face_no = 3
         790             :       else
         791           0 :          face_no = 4
         792             :       endif
         793             :    endif
         794             :
         795           0 :    end function cube_face_number_from_cart
         796             :
         797             : ! This could be done directly by using the lon, lat coordinates,
         798             : ! but call cube_face_number_from_cart just so that there is one place
         799             : ! to do the conversions and they are all consistant.
         800           0 :   pure function cube_face_number_from_sphere(sphere) result(face_no)
         801             :     implicit none
         802             :     type (spherical_polar_t), intent(in) :: sphere
         803             :     type (cartesian3d_t)                 :: cart
         804             :     integer                              :: face_no
         805             :
         806           0 :     cart =  spherical_to_cart(sphere)
         807           0 :     face_no = cube_face_number_from_cart(cart)
         808           0 :   end function cube_face_number_from_sphere
         809             :
         810             : ! CE, need real (cartesian) xy coordinates on the cubed sphere
         811           0 : subroutine cart2cubedspherexy(cart3d,face_no,cartxy)
         812             :
         813             :   type(cartesian3D_t),intent(in)      :: cart3d
         814             :   integer, intent(in)                 :: face_no
         815             :   type (cartesian2d_t), intent(out)   :: cartxy
         816             :
         817             :   ! a (half length of a cube side) is supposed to be 1
         818           0 :   select case (face_no)
         819             :   case (1)
         820           0 :      cartxy%x = cart3D%y/cart3D%x
         821           0 :      cartxy%y = cart3D%z/cart3D%x
         822             :   case (2)
         823           0 :      cartxy%x = -cart3D%x/cart3D%y
         824           0 :      cartxy%y = cart3D%z/cart3D%y
         825             :   case (3)
         826           0 :      cartxy%x = cart3D%y/cart3D%x
         827           0 :      cartxy%y = -cart3D%z/cart3D%x
         828             :   case (4)
         829           0 :      cartxy%x = -cart3D%x/cart3D%y
         830           0 :      cartxy%y = -cart3D%z/cart3D%y
         831             :   case (5)       !bottom face
         832           0 :      cartxy%x  = -cart3D%y/cart3D%z
         833           0 :      cartxy%y  = -cart3D%x/cart3D%z
         834             :   case (6)        !top face
         835           0 :      cartxy%x  = cart3D%y/cart3D%z
         836           0 :      cartxy%y  = -cart3D%x/cart3D%z
         837             :   end select
         838           0 : end subroutine cart2cubedspherexy
         839             : ! CE END
         840             :
         841             :
         842             :
         843             :
         844           0 :   subroutine sphere_tri_area( v1, v2, v3, area )
         845             :   !  input: v1(3),v2(3),v3(3)  cartesian coordinates of triangle
         846             :   !  output: area
         847             :   !  based on formulas in STRI_QUAD:
         848             :   !  http://people.sc.fsu.edu/~burkardt/f_src/stri_quad/stri_quad.html
         849             :   !
         850             :   ! MT:  note that the usual Girard formula used below is ill-conditioned
         851             :   ! for small nearly flat triangles (which is probably our case).
         852             :   ! I put in a crude fix below.
         853             :   ! We should instead be using l'Huiller's formula,
         854             :   ! see:  http://williams.best.vwh.net/avform.htm
         855             :
         856             :   real(kind=r8) area
         857             :   real(kind=r8) a,b,c,al,bl,cl,sina,sinb,sinc,sins,a1,b1,c1
         858             :   type (cartesian3D_t) v1,v2,v3
         859             :
         860             :   ! compute great circle lengths
         861           0 :   al = acos( v2%x * v3%x + v2%y * v3%y + v2%z * v3%z )
         862           0 :   bl = acos( v3%x * v1%x + v3%y * v1%y + v3%z * v1%z )
         863           0 :   cl = acos( v1%x * v2%x + v1%y * v2%y + v1%z * v2%z )
         864             :
         865             :   ! compute angles
         866           0 :   sina = sin( (bl+cl-al)/2 )  ! sin(sl-al)
         867           0 :   sinb = sin( (al+cl-bl)/2 )  ! sin(sl-bl)
         868           0 :   sinc = sin( (al+bl-cl)/2 )
         869           0 :   sins = sin( (al+bl+cl)/2 )
         870             :
         871             :   ! for small areas, formula above looses precision.
         872             :   ! 2atan(x) + 2atan(1/x) = pi
         873             :   ! 2atan(x) - pi = -2atan(1/x)
         874           0 :   a = sqrt( (sinb*sinc) / (sins*sina) )
         875           0 :   b = sqrt( (sina*sinc) / (sins*sinb) )
         876           0 :   c = sqrt( (sina*sinb) / (sins*sinc) )
         877             :
         878             :
         879           0 :   a1 = 2*atan(a)
         880           0 :   b1 = 2*atan(b)
         881           0 :   c1 = 2*atan(c)
         882             :
         883           0 :   if (a.gt.b.and.a.gt.c) then
         884           0 :      a1 = -2*atan(1/a)
         885           0 :   else if (b.gt.c) then
         886           0 :      b1 = -2*atan(1/b)
         887             :   else
         888           0 :      c1 = -2*atan(1/c)
         889             :   endif
         890             :   ! apply Girard's theorem
         891           0 :   area = a1+b1+c1
         892             :
         893             :
         894           0 :   end subroutine sphere_tri_area
         895             :
         896             : !INPUT: Points in xy cubed sphere coordinates, counterclockwise
         897             : !OUTPUT: corresponding area on the sphere
         898           0 : function surfareaxy(x1,x2,y1,y2) result(area)
         899             :   implicit none
         900             :   real (kind=r8), intent(in)    :: x1, x2, y1, y2
         901             :   real (kind=r8)   :: area
         902             :   real (kind=r8)   :: a1,a2,a3,a4
         903             :
         904             :   ! cubed-sphere cell area, from Lauritzen & Nair MWR 2008
         905             :   ! central angles:
         906             :   ! cube face: -pi/4,-pi/4 -> pi/4,pi/4
         907             :   ! this formula gives 2   so normalize by 4pi/6 / 2 = pi/3
         908             :   ! use implementation where the nodes a counterclockwise (not as in the paper)
         909           0 :   a1 = acos(-sin(atan(x1))*sin(atan(y1)))
         910           0 :   a2 =-acos(-sin(atan(x2))*sin(atan(y1)))
         911           0 :   a3 = acos(-sin(atan(x2))*sin(atan(y2)))
         912           0 :   a4 =-acos(-sin(atan(x1))*sin(atan(y2)))
         913           0 :   area = (a1+a2+a3+a4)
         914             :   return
         915             : end function surfareaxy
         916             :
         917             :
         918             :
         919           0 : end module coordinate_systems_mod

Generated by: LCOV version 1.14