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
|
