FE-Project
Loading...
Searching...
No Matches
scale_cubedsphere_coord_cnv.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
20#include "scaleFElib.h"
21module scale_cubedsphere_coord_cnv
22 !-----------------------------------------------------------------------------
23 !
24 !++ used modules
25 !
26 use scale_const, only: &
27 pi => const_pi, &
28 eps => const_eps, &
29 rplanet => const_radius
30 use scale_precision
31 use scale_prc
32 use scale_io
33 use scale_prc
34
35 !-----------------------------------------------------------------------------
36 implicit none
37 private
38
39 !-----------------------------------------------------------------------------
40 !
41 !++ Public type & procedure
42 !
43 public :: cubedspherecoordcnv_cs2lonlatpos
44 public :: cubedspherecoordcnv_cs2lonlatvec
45 public :: cubedspherecoordcnv_lonlat2cspos
46 public :: cubedspherecoordcnv_lonlat2csvec
47 public :: cubedspherecoordcnv_cs2cartpos
48 public :: cubedspherecoordcnv_cart2csvec
49 public :: cubedspherecoordcnv_getmetric
50
51 interface cubedspherecoordcnv_cs2localorthvec_alpha
52 module procedure :: cubedspherecoordcnv_cs2localorthvec_alpha_0
53 module procedure :: cubedspherecoordcnv_cs2localorthvec_alpha_1
54 end interface
55 public :: cubedspherecoordcnv_cs2localorthvec_alpha
56
57 interface cubedspherecoordcnv_localorth2csvec_alpha
58 module procedure :: cubedspherecoordcnv_localorth2csvec_alpha_0
59 module procedure :: cubedspherecoordcnv_localorth2csvec_alpha_1
60 end interface
61 public :: cubedspherecoordcnv_localorth2csvec_alpha
62
63 !-----------------------------------------------------------------------------
64 !
65 !++ Public parameters & variables
66 !
67 !-----------------------------------------------------------------------------
68
69 !
70 !++ Private type & procedure
71 !
72 !-----------------------------------------------------------------------------
73
74 !-----------------------------------------------------------------------------
75 !
76 !++ Private parameters & variables
77 !
78 !-----------------------------------------------------------------------------
79
80
81contains
82!OCL SERIAL
83 subroutine cubedspherecoordcnv_cs2lonlatpos( &
84 panelID, alpha, beta, gam, Np, & ! (in)
85 lon, lat ) ! (out)
86 use scale_geographic_coord_cnv, only: &
87 geographiccoordcnv_orth_to_geo_pos
88 implicit none
89
90 integer, intent(in) :: panelid
91 integer, intent(in) :: np
92 real(rp), intent(in) :: alpha(np)
93 real(rp), intent(in) :: beta (np)
94 real(rp), intent(in) :: gam(np)
95 real(rp), intent(out) :: lon(np)
96 real(rp), intent(out) :: lat(np)
97
98 real(rp) :: cartpos(np,3)
99 real(rp) :: geogpos(np,3)
100
101 integer :: p
102 !-----------------------------------------------------------------------------
103
104 select case( panelid )
105 case( 1, 2, 3, 4 )
106 !$omp parallel
107 !$omp do
108 do p=1, np
109 lon(p) = alpha(p) + 0.5_rp * pi * dble(panelid - 1)
110 lat(p) = atan( tan( beta(p) ) * cos( alpha(p) ) )
111 end do
112 !$omp workshare
113 where( lon(:) < 0.0_rp )
114 lon(:) = lon(:) + 2.0_rp * pi
115 end where
116 !$omp end workshare
117 !$omp end parallel
118 case(5, 6)
119 call cubedspherecoordcnv_cs2cartpos( panelid, alpha, beta, gam, np, & ! (in)
120 cartpos(:,1), cartpos(:,2), cartpos(:,3) ) ! (out)
121
122 call geographiccoordcnv_orth_to_geo_pos( cartpos(:,:), np, & ! (in)
123 geogpos(:,:) ) ! (out)
124
125 !$omp parallel
126 !$omp do
127 do p=1, np
128 lon(p) = geogpos(p,1)
129 lat(p) = geogpos(p,2)
130 end do
131 !$omp workshare
132 where( alpha < 0.0_rp )
133 lon(:) = lon(:) + 2.0_rp * pi
134 end where
135 !$omp end workshare
136 !$omp end parallel
137 case default
138 log_error("CubedSphereCoordCnv_CS2LonLatPos",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
139 call prc_abort
140 end select
141
142 return
143 end subroutine cubedspherecoordcnv_cs2lonlatpos
144
145!OCL SERIAL
146 subroutine cubedspherecoordcnv_cs2lonlatvec( &
147 panelID, alpha, beta, gam, Np, & ! (in)
148 vecalpha, vecbeta, & ! (in)
149 veclon, veclat, & ! (out)
150 lat ) ! (in, optional)
151
152 use scale_const, only: &
153 eps => const_eps
154 implicit none
155
156 integer, intent(in) :: panelid
157 integer, intent(in) :: np
158 real(rp), intent(in) :: alpha(np)
159 real(rp), intent(in) :: beta (np)
160 real(rp), intent(in) :: gam(np) ! RPlanet / r
161 real(dp), intent(in) :: vecalpha(np)
162 real(dp), intent(in) :: vecbeta (np)
163 real(rp), intent(out) :: veclon(np)
164 real(rp), intent(out) :: veclat(np)
165 real(rp), intent(in), optional :: lat(np)
166
167 integer :: p
168 real(rp) :: x ,y, del2
169 real(rp) :: s
170
171 real(rp) :: radius
172 real(rp) :: cos_lat(np)
173 !-----------------------------------------------------------------------------
174
175 if (present(lat)) then
176 !$omp parallel do
177 do p=1, np
178 cos_lat(p) = cos(lat(p))
179 end do
180 end if
181
182 select case( panelid )
183 case( 1, 2, 3, 4 )
184 if (.not. present(lat)) then
185 !$omp parallel do
186 do p=1, np
187 cos_lat(p) = cos( atan( tan( beta(p) ) * cos( alpha(p) ) ) )
188 end do
189 end if
190
191 !$omp parallel do private( X, Y, del2, radius )
192 do p=1, np
193 x = tan( alpha(p) )
194 y = tan( beta(p) )
195 del2 = 1.0_rp + x**2 + y**2
196 radius = rplanet * gam(p)
197
198 veclon(p) = vecalpha(p) * cos_lat(p) * radius
199 veclat(p) = ( - x * y * vecalpha(p) + ( 1.0_rp + y**2 ) * vecbeta(p) ) &
200 * radius * sqrt( 1.0_rp + x**2 ) / del2
201 end do
202 case ( 5 )
203 s = 1.0_rp
204 case ( 6 )
205 s = - 1.0_rp
206 case default
207 log_error("CubedSphereCoordCnv_CS2LonLatVec",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
208 call prc_abort
209 end select
210
211 select case( panelid )
212 case( 5, 6 )
213 !$omp parallel do private( X, Y, del2, radius )
214 do p=1, np
215 x = tan( alpha(p) )
216 y = tan( beta(p) )
217 del2 = 1.0_rp + x**2 + y**2
218 radius = s * rplanet * gam(p)
219
220 if (.not. present(lat)) then
221 cos_lat(p) = cos( atan( sign(1.0_rp, s) / max( sqrt( x**2 + y**2 ), eps ) ) )
222 end if
223
224 veclon(p) = (- y * ( 1.0_rp + x**2 ) * vecalpha(p) + x * ( 1.0_rp + y**2 ) * vecbeta(p) ) &
225 * radius / max( x**2 + y**2, eps ) * cos_lat(p)
226 veclat(p) = (- x * ( 1.0_rp + x**2 ) * vecalpha(p) - y * ( 1.0_rp + y**2 ) * vecbeta(p) ) &
227 * radius / ( del2 * ( max( sqrt( x**2 + y**2 ), eps ) ) )
228 end do
229 end select
230
231 return
232 end subroutine cubedspherecoordcnv_cs2lonlatvec
233
234!OCL SERIAL
235 subroutine cubedspherecoordcnv_lonlat2cspos( &
236 panelID, lon, lat, Np, & ! (in)
237 alpha, beta ) ! (out)
238
239 implicit none
240
241 integer, intent(in) :: panelid
242 integer, intent(in) :: np
243 real(rp), intent(in) :: lon(np)
244 real(rp), intent(in) :: lat(np)
245 real(rp), intent(out) ::alpha(np)
246 real(rp), intent(out) ::beta (np)
247
248 integer :: p
249 real(rp) :: tan_lat
250 real(rp) :: lon_(np)
251 !-----------------------------------------------------------------------------
252
253 select case(panelid)
254 case ( 1, 2, 3, 4 )
255 !$omp parallel
256 if ( panelid == 1 ) then
257 !$omp workshare
258 where (lon(:) >= 2.0_rp * pi - 0.25_rp * pi )
259 lon_(:) = lon(:) - 2.0_rp * pi
260 elsewhere
261 lon_(:) = lon(:)
262 end where
263 !$omp end workshare
264 else
265 !$omp workshare
266 where (lon(:) < 0.0_rp )
267 lon_(:) = lon(:) + 2.0_rp * pi
268 elsewhere
269 lon_(:) = lon(:)
270 end where
271 !$omp end workshare
272 end if
273 !$omp do
274 do p=1, np
275 alpha(p) = lon_(p) - 0.5_rp * pi * ( dble(panelid) - 1.0_rp )
276 beta(p) = atan( tan(lat(p)) / cos(alpha(p)) )
277 end do
278 !$omp end parallel
279 case ( 5 )
280 !$omp parallel private(tan_lat)
281 !$omp do
282 do p=1, np
283 tan_lat = tan(lat(p))
284 alpha(p) = + atan( sin(lon(p)) / tan_lat )
285 beta(p) = - atan( cos(lon(p)) / tan_lat )
286 end do
287 !$omp end parallel
288 case ( 6 )
289 !$omp parallel private(tan_lat)
290 !$omp do
291 do p=1, np
292 tan_lat = tan(lat(p))
293 alpha(p) = - atan( sin(lon(p)) / tan_lat )
294 beta(p) = - atan( cos(lon(p)) / tan_lat )
295 end do
296 !$omp end parallel
297 case default
298 log_error("CubedSphereCoordCnv_LonLat2CSPos",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
299 call prc_abort
300 end select
301
302 return
303 end subroutine cubedspherecoordcnv_lonlat2cspos
304
305 !
306 !
307!OCL SERIAL
308 subroutine cubedspherecoordcnv_lonlat2csvec( &
309 panelID, alpha, beta, gam, Np, & ! (in)
310 veclon, veclat, & ! (in)
311 vecalpha, vecbeta, & ! (out)
312 lat ) ! (in, optional)
313
314 implicit none
315
316 integer, intent(in) :: panelid
317 integer, intent(in) :: np
318 real(rp), intent(in) :: alpha(np)
319 real(rp), intent(in) :: beta (np)
320 real(rp), intent(in) :: gam(np)
321 real(dp), intent(in) :: veclon(np)
322 real(dp), intent(in) :: veclat(np)
323 real(rp), intent(out) :: vecalpha(np)
324 real(rp), intent(out) :: vecbeta (np)
325 real(rp), intent(in), optional :: lat(np)
326
327 integer :: p
328 real(rp) :: x ,y, del2
329 real(rp) :: s
330
331 real(rp) :: radius
332 real(rp) :: cos_lat(np)
333 real(rp) :: veclon_ov_coslat
334 !-----------------------------------------------------------------------------
335
336 if (present(lat)) then
337 !$omp parallel do
338 do p=1, np
339 cos_lat(p) = cos(lat(p))
340 end do
341 end if
342
343 select case( panelid )
344 case( 1, 2, 3, 4 )
345 if (.not. present(lat)) then
346 !$omp parallel do
347 do p=1, np
348 cos_lat(p) = cos( atan( tan( beta(p) ) * cos( alpha(p) ) ) )
349 end do
350 end if
351
352 !$omp parallel do private( X, Y, del2, radius, VecLon_ov_cosLat )
353 do p=1, np
354 x = tan( alpha(p) )
355 y = tan( beta(p) )
356 del2 = 1.0_rp + x**2 + y**2
357 radius = rplanet * gam(p)
358 veclon_ov_coslat = veclon(p) / cos_lat(p)
359
360 vecalpha(p) = veclon_ov_coslat / radius
361 vecbeta(p) = ( x * y * veclon_ov_coslat + del2 / sqrt( 1.0_rp + x**2 ) * veclat(p) ) &
362 / ( radius * (1.0_rp + y**2) )
363 end do
364 case ( 5 )
365 s = 1.0_rp
366 case ( 6 )
367 s = -1.0_rp
368 case default
369 log_error("CubedSphereCoordCnv_LonLat2CSVec",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
370 call prc_abort
371 end select
372
373 select case( panelid )
374 case( 5, 6 )
375 !$omp parallel do private( X, Y, del2, radius, VecLon_ov_cosLat )
376 do p=1, np
377 x = tan( alpha(p) )
378 y = tan( beta(p) )
379 del2 = 1.0_rp + x**2 + y**2
380 radius = s * rplanet * gam(p)
381
382 if (.not. present(lat)) then
383 cos_lat(p) = cos( atan( s / sqrt(max(x**2 + y**2, eps)) ) )
384 end if
385 veclon_ov_coslat = veclon(p) / cos_lat(p)
386
387 vecalpha(p) = (- y * veclon_ov_coslat - del2 * x / sqrt( max(del2 - 1.0_rp,eps)) * veclat(p)) &
388 / ( radius * ( 1.0_rp + x**2 ) )
389 vecbeta(p) = ( x * veclon_ov_coslat - del2 * y / sqrt( max(del2 - 1.0_rp,eps)) * veclat(p)) &
390 / ( radius * ( 1.0_rp + y**2 ) )
391 end do
392 end select
393
394 return
395 end subroutine cubedspherecoordcnv_lonlat2csvec
396
397!OCL SERIAL
398 subroutine cubedspherecoordcnv_cs2cartpos( &
399 panelID, alpha, beta, gam, Np, & ! (in)
400 x, y, z ) ! (out)
401
402 implicit none
403 integer, intent(in) :: panelid
404 integer, intent(in) :: np
405 real(rp), intent(in) :: alpha(np)
406 real(rp), intent(in) :: beta (np)
407 real(rp), intent(in) :: gam(np)
408 real(rp), intent(out) :: x(np)
409 real(rp), intent(out) :: y(np)
410 real(rp), intent(out) :: z(np)
411
412 integer :: p
413 real(rp) :: x1, x2, fac
414
415 !-----------------------------------------------------------------------------
416
417 select case(panelid)
418 case(1)
419 !$omp parallel do private(x1, x2, fac)
420 do p=1, np
421 x1 = tan( alpha(p) )
422 x2 = tan( beta(p) )
423 fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
424 x(p) = fac
425 y(p) = fac * x1
426 z(p) = fac * x2
427 end do
428 case(2)
429 !$omp parallel do private(x1, x2, fac)
430 do p=1, np
431 x1 = tan( alpha(p) )
432 x2 = tan( beta(p) )
433 fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
434 x(p) = - fac * x1
435 y(p) = fac
436 z(p) = fac * x2
437 end do
438 case(3)
439 !$omp parallel do private(x1, x2, fac)
440 do p=1, np
441 x1 = tan( alpha(p) )
442 x2 = tan( beta(p) )
443 fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
444 x(p) = - fac
445 y(p) = - fac * x1
446 z(p) = fac * x2
447 end do
448 case(4)
449 !$omp parallel do private(x1, x2, fac)
450 do p=1, np
451 x1 = tan( alpha(p) )
452 x2 = tan( beta(p) )
453 fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
454 x(p) = fac * x1
455 y(p) = - fac
456 z(p) = fac * x2
457 end do
458 case(5)
459 !$omp parallel do private(x1, x2, fac)
460 do p=1, np
461 x1 = tan( alpha(p) )
462 x2 = tan( beta(p) )
463 fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
464 x(p) = - fac * x2
465 y(p) = fac * x1
466 z(p) = fac
467 end do
468 case(6)
469 !$omp parallel do private(x1, x2, fac)
470 do p=1, np
471 x1 = tan( alpha(p) )
472 x2 = tan( beta(p) )
473 fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
474 x(p) = fac * x2
475 y(p) = fac * x1
476 z(p) = - fac
477 end do
478 case default
479 log_error("CubedSphereCoordCnv_CS2CartPos",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
480 call prc_abort
481 end select
482
483 return
484 end subroutine cubedspherecoordcnv_cs2cartpos
485
486 !
487 !
488!OCL SERIAL
489 subroutine cubedspherecoordcnv_cart2csvec( &
490 panelID, alpha, beta, gam, Np, & ! (in)
491 vec_x, vec_y, vec_z, & ! (in)
492 vecalpha, vecbeta ) ! (out)
493
494 implicit none
495
496 integer, intent(in) :: panelid
497 integer, intent(in) :: np
498 real(rp), intent(in) :: alpha(np)
499 real(rp), intent(in) :: beta (np)
500 real(rp), intent(in) :: gam(np)
501 real(dp), intent(in) :: vec_x(np)
502 real(dp), intent(in) :: vec_y(np)
503 real(dp), intent(in) :: vec_z(np)
504 real(rp), intent(out) :: vecalpha(np)
505 real(rp), intent(out) :: vecbeta (np)
506
507 integer :: p
508 real(rp) :: x1, x2, fac
509 real(rp) :: r_sec2_alpha, r_sec2_beta
510 !-----------------------------------------------------------------------------
511
512 select case( panelid )
513 case(1)
514 !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
515 do p=1, np
516 x1 = tan( alpha(p) )
517 x2 = tan( beta(p) )
518 r_sec2_alpha = cos(alpha(p))**2
519 r_sec2_beta = cos(beta(p))**2
520 fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
521
522 vecalpha(p) = fac * r_sec2_alpha * ( - x1 * vec_x(p) + vec_y(p) )
523 vecbeta(p) = fac * r_sec2_beta * ( - x2 * vec_x(p) + vec_z(p) )
524 end do
525 case(2)
526 !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
527 do p=1, np
528 x1 = tan( alpha(p) )
529 x2 = tan( beta(p) )
530 r_sec2_alpha = cos(alpha(p))**2
531 r_sec2_beta = cos(beta(p))**2
532 fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
533
534 vecalpha(p) = - fac * r_sec2_alpha * ( vec_x(p) + x1 * vec_y(p) )
535 vecbeta(p) = fac * r_sec2_beta * ( - x2 * vec_y(p) + vec_z(p) )
536 end do
537 case(3)
538 !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
539 do p=1, np
540 x1 = tan( alpha(p) )
541 x2 = tan( beta(p) )
542 r_sec2_alpha = cos(alpha(p))**2
543 r_sec2_beta = cos(beta(p))**2
544 fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
545
546 vecalpha(p) = fac * r_sec2_alpha * ( x1 * vec_x(p) - vec_y(p) )
547 vecbeta(p) = fac * r_sec2_beta * ( x2 * vec_x(p) + vec_z(p) )
548 end do
549 case(4)
550 !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
551 do p=1, np
552 x1 = tan( alpha(p) )
553 x2 = tan( beta(p) )
554 r_sec2_alpha = cos(alpha(p))**2
555 r_sec2_beta = cos(beta(p))**2
556 fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
557
558 vecalpha(p) = fac * r_sec2_alpha * ( vec_x(p) + x1 * vec_y(p) )
559 vecbeta(p) = fac * r_sec2_beta * ( x2 * vec_y(p) + vec_z(p) )
560 end do
561 case ( 5 )
562 !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
563 do p=1, np
564 x1 = tan( alpha(p) )
565 x2 = tan( beta(p) )
566 r_sec2_alpha = cos(alpha(p))**2
567 r_sec2_beta = cos(beta(p))**2
568 fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
569
570 vecalpha(p) = fac * r_sec2_alpha * ( vec_y(p) - x1 * vec_z(p) )
571 vecbeta(p) = - fac * r_sec2_beta * ( vec_x(p) + x2 * vec_z(p) )
572 end do
573 case ( 6 )
574 !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
575 do p=1, np
576 x1 = tan( alpha(p) )
577 x2 = tan( beta(p) )
578 r_sec2_alpha = cos(alpha(p))**2
579 r_sec2_beta = cos(beta(p))**2
580 fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
581
582 vecalpha(p) = fac * r_sec2_alpha * ( vec_y(p) + x1 * vec_z(p) )
583 vecbeta(p) = fac * r_sec2_beta * ( vec_x(p) + x2 * vec_z(p) )
584 end do
585 case default
586 log_error("CubedSphereCoordCnv_Cart2CSVec",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
587 call prc_abort
588 end select
589
590 return
591 end subroutine cubedspherecoordcnv_cart2csvec
592
593!OCL SERIAL
594 subroutine cubedspherecoordcnv_cs2localorthvec_alpha_0( &
595 alpha, beta, radius, Np, & ! (in)
596 vecorth1, vecorth2 ) ! (inout)
597
598 implicit none
599
600 integer, intent(in) :: Np
601 real(RP), intent(in) :: alpha(Np)
602 real(RP), intent(in) :: beta (Np)
603 real(RP), intent(in) :: radius(Np)
604 real(RP), intent(inout) :: VecOrth1(Np)
605 real(RP), intent(inout) :: VecOrth2(Np)
606
607 integer :: p
608 real(RP) :: x1, x2, del
609 real(RP) :: fac
610 real(RP) :: tmp
611 !-----------------------------------------------------------------------------
612
613 !$omp parallel do private(x1, x2, del, fac, tmp)
614 do p=1, np
615 x1 = tan( alpha(p) )
616 x2 = tan( beta(p) )
617 del = sqrt(1.0_rp + x1**2 + x2**2)
618 fac = radius(p) * sqrt(1.0_rp + x1**2) / del**2
619
620 tmp = vecorth1(p)
621 vecorth1(p) = fac * del * tmp
622 vecorth2(p) = fac * ( - x1 * x2 * tmp + (1.0_rp + x2**2) * vecorth2(p) )
623 end do
624
625 return
626 end subroutine cubedspherecoordcnv_cs2localorthvec_alpha_0
627
628!OCL SERIAL
629 subroutine cubedspherecoordcnv_cs2localorthvec_alpha_1( &
630 alpha, beta, radius, Np, & ! (in)
631 vecalpha, vecbeta, & ! (in)
632 vecorth1, vecorth2 ) ! (out)
633
634 implicit none
635
636 integer, intent(in) :: Np
637 real(RP), intent(in) :: alpha(Np)
638 real(RP), intent(in) :: beta (Np)
639 real(RP), intent(in) :: radius(Np)
640 real(RP), intent(in) :: VecAlpha(Np)
641 real(RP), intent(in) :: VecBeta (Np)
642 real(RP), intent(out) :: VecOrth1(Np)
643 real(RP), intent(out) :: VecOrth2(Np)
644
645 !-----------------------------------------------------------------------------
646
647 !$omp parallel workshare
648 vecorth1(:) = vecalpha(:)
649 vecorth2(:) = vecbeta(:)
650 !$omp end parallel workshare
651
652 call cubedspherecoordcnv_cs2localorthvec_alpha_0( &
653 alpha, beta, radius, np, & ! (in)
654 vecorth1, vecorth2 ) ! (inout)
655
656 return
657 end subroutine cubedspherecoordcnv_cs2localorthvec_alpha_1
658
659!OCL SERIAL
660 subroutine cubedspherecoordcnv_localorth2csvec_alpha_0( &
661 alpha, beta, radius, Np, & ! (in)
662 vecalpha, vecbeta ) ! (inout)
663
664 implicit none
665
666 integer, intent(in) :: Np
667 real(RP), intent(in) :: alpha(Np)
668 real(RP), intent(in) :: beta (Np)
669 real(RP), intent(in) :: radius(Np)
670 real(RP), intent(inout) :: VecAlpha(Np)
671 real(RP), intent(inout) :: VecBeta (Np)
672
673 integer :: p
674 real(RP) :: x1, x2, del
675 real(RP) :: fac
676 real(RP) :: tmp
677 !-----------------------------------------------------------------------------
678
679 !$omp parallel do private(x1, x2, del, fac, tmp)
680 do p=1, np
681 x1 = tan( alpha(p) )
682 x2 = tan( beta(p) )
683 del = sqrt(1.0_rp + x1**2 + x2**2)
684 fac = del / ( radius(p) * ( 1.0_rp + x2**2 ) * sqrt( 1.0_rp + x1**2 ) )
685
686 tmp = vecalpha(p)
687 vecalpha(p) = fac * ( 1.0_rp + x2**2 ) * tmp
688 vecbeta(p) = fac * ( x1 * x2 * tmp + del * vecbeta(p) )
689 end do
690
691 return
692 end subroutine cubedspherecoordcnv_localorth2csvec_alpha_0
693
694!OCL SERIAL
695 subroutine cubedspherecoordcnv_localorth2csvec_alpha_1( &
696 alpha, beta, radius, Np, & ! (in)
697 vecorth1, vecorth2, & ! (in)
698 vecalpha, vecbeta ) ! (out)
699
700 implicit none
701
702 integer, intent(in) :: Np
703 real(RP), intent(in) :: alpha(Np)
704 real(RP), intent(in) :: beta (Np)
705 real(RP), intent(in) :: radius(Np)
706 real(RP), intent(in) :: VecOrth1(Np)
707 real(RP), intent(in) :: VecOrth2(Np)
708 real(RP), intent(out) :: VecAlpha(Np)
709 real(RP), intent(out) :: VecBeta (Np)
710 !-----------------------------------------------------------------------------
711
712 !$omp parallel workshare
713 vecalpha(:) = vecorth1(:)
714 vecbeta(:) = vecorth2(:)
715 !$omp end parallel workshare
716
717 call cubedspherecoordcnv_localorth2csvec_alpha_0( &
718 alpha, beta, radius, np, & ! (in)
719 vecalpha, vecbeta ) ! (inout)
720
721 return
722 end subroutine cubedspherecoordcnv_localorth2csvec_alpha_1
723
724!OCL SERIAL
725 subroutine cubedspherecoordcnv_getmetric( &
726 alpha, beta, Np, radius, & ! (in)
727 g_ij, gij, gsqrt ) ! (out)
728
729 implicit none
730
731 integer, intent(in) :: np
732 real(rp), intent(in) :: alpha(np)
733 real(rp), intent(in) :: beta(np)
734 real(rp), intent(in) :: radius
735 real(rp), intent(out) :: g_ij(np,2,2)
736 real(rp), intent(out) :: gij (np,2,2)
737 real(rp), intent(out) :: gsqrt(np)
738
739 real(rp) :: x, y
740 real(rp) :: r2
741 real(rp) :: fac
742 real(rp) :: g_ij_(2,2)
743 real(rp) :: gij_ (2,2)
744
745 integer :: p
746 real(rp) :: oneplusx2, oneplusy2
747 !-----------------------------------------------------------------------------
748
749 !$omp parallel do private( &
750 !$omp X, Y, r2, OnePlusX2, OnePlusY2, fac, &
751 !$omp G_ij_, GIJ_ )
752 do p=1, np
753 x = tan(alpha(p))
754 y = tan(beta(p))
755 r2 = 1.0_rp + x**2 + y**2
756 oneplusx2 = 1.0_rp + x**2
757 oneplusy2 = 1.0_rp + y**2
758
759 fac = oneplusx2 * oneplusy2 * ( radius / r2 )**2
760 g_ij_(1,1) = fac * oneplusx2
761 g_ij_(1,2) = - fac * (x * y)
762 g_ij_(2,1) = - fac * (x * y)
763 g_ij_(2,2) = fac * oneplusy2
764 g_ij(p,:,:) = g_ij_(:,:)
765
766 gsqrt(p) = radius**2 * oneplusx2 * oneplusy2 / ( r2 * sqrt(r2) )
767
768 fac = 1.0_rp / gsqrt(p)**2
769 gij_(1,1) = fac * g_ij_(2,2)
770 gij_(1,2) = - fac * g_ij_(1,2)
771 gij_(2,1) = - fac * g_ij_(2,1)
772 gij_(2,2) = fac * g_ij_(1,1)
773 gij(p,:,:) = gij_(:,:)
774 end do
775
776 return
777 end subroutine cubedspherecoordcnv_getmetric
778
779end module scale_cubedsphere_coord_cnv
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cs2localorthvec_alpha
Definition scale_cubedsphere_coord_cnv.F90:51
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_localorth2csvec_alpha
Definition scale_cubedsphere_coord_cnv.F90:57
scale_cubedsphere_coord_cnv
module common / Coordinate conversion with a cubed-sphere
Definition scale_cubedsphere_coord_cnv.F90:21
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cs2lonlatvec
subroutine, public cubedspherecoordcnv_cs2lonlatvec(panelid, alpha, beta, gam, np, vecalpha, vecbeta, veclon, veclat, lat)
Definition scale_cubedsphere_coord_cnv.F90:151
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cs2lonlatpos
subroutine, public cubedspherecoordcnv_cs2lonlatpos(panelid, alpha, beta, gam, np, lon, lat)
Definition scale_cubedsphere_coord_cnv.F90:86
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_lonlat2csvec
subroutine, public cubedspherecoordcnv_lonlat2csvec(panelid, alpha, beta, gam, np, veclon, veclat, vecalpha, vecbeta, lat)
Definition scale_cubedsphere_coord_cnv.F90:313
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_getmetric
subroutine, public cubedspherecoordcnv_getmetric(alpha, beta, np, radius, g_ij, gij, gsqrt)
Definition scale_cubedsphere_coord_cnv.F90:728
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_lonlat2cspos
subroutine, public cubedspherecoordcnv_lonlat2cspos(panelid, lon, lat, np, alpha, beta)
Definition scale_cubedsphere_coord_cnv.F90:238
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cart2csvec
subroutine, public cubedspherecoordcnv_cart2csvec(panelid, alpha, beta, gam, np, vec_x, vec_y, vec_z, vecalpha, vecbeta)
Definition scale_cubedsphere_coord_cnv.F90:493
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cs2cartpos
subroutine, public cubedspherecoordcnv_cs2cartpos(panelid, alpha, beta, gam, np, x, y, z)
Definition scale_cubedsphere_coord_cnv.F90:401
scale_geographic_coord_cnv
module common / Coordinate conversion with a geographic coordinate
Definition scale_geographic_coord_cnv.F90:9
scale_geographic_coord_cnv::geographiccoordcnv_orth_to_geo_pos
subroutine, public geographiccoordcnv_orth_to_geo_pos(orth_p, np, geo_p)
Definition scale_geographic_coord_cnv.F90:44
Generated by doxygen 1.10.0