scale_cubedsphere_coord_cnv.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> Module common / Coordinate conversion with cubed-sphere projection
!!
!! @par Description
!!      A module to provide coordinate conversions with an equiangular gnomonic cubed-sphere projection
!!
!!
!! @par Reference
!!  - Nair et al. 2015:
!!    A Discontinuous Galerkin Transport Scheme on the Cubed Sphere. 
!!    Monthly Weather Review, 133, 814–828. 
!!    (Appendix A)
!!  - Yin et al. 2017:
!!    Parallel numerical simulation of the thermal convection in the Earth’s outer core on the cubed-sphere. 
!!    Geophysical Journal International, 209, 1934–1954. 
!!    (Appendix A)
!!
!! @author Yuta Kawai, Team SCALE
!!
#include "scaleFElib.h"
module scale_cubedsphere_coord_cnv
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_const, only: &
    pi => const_pi,         &
    eps => const_eps,       &
    rplanet => const_radius
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prc
 
#ifdef _OPENACC
    use openacc, only: acc_async_sync
#endif
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
  public :: cubedspherecoordcnv_cs2lonlatpos
  public :: cubedspherecoordcnv_cs2lonlatvec
  public :: cubedspherecoordcnv_lonlat2cspos
  public :: cubedspherecoordcnv_lonlat2csvec
  public :: cubedspherecoordcnv_cs2cartpos
  public :: cubedspherecoordcnv_cart2csvec  
  public :: cubedspherecoordcnv_getmetric
 
  interface cubedspherecoordcnv_cs2localorthvec_alpha
    module procedure :: cubedspherecoordcnv_cs2localorthvec_alpha_0
    module procedure :: cubedspherecoordcnv_cs2localorthvec_alpha_1
  end interface
  public :: cubedspherecoordcnv_cs2localorthvec_alpha
  
  interface cubedspherecoordcnv_localorth2csvec_alpha
    module procedure :: cubedspherecoordcnv_localorth2csvec_alpha_0
    module procedure :: cubedspherecoordcnv_localorth2csvec_alpha_1
  end interface
  public :: cubedspherecoordcnv_localorth2csvec_alpha
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !  
  !-----------------------------------------------------------------------------
 
 
contains
  !> Calculate longitude and latitude coordinates from local coordinates using the central angles in an equiangular gnomonic cubed-sphere projection
  !!
!OCL SERIAL
  subroutine cubedspherecoordcnv_cs2lonlatpos( &
    panelID, alpha, beta, gam, Np,             & ! (in)
    lon, lat )                                   ! (out)
    use scale_geographic_coord_cnv, only: &
      geographiccoordcnv_orth_to_geo_pos
    implicit none
 
    integer, intent(in) :: panelid     !< Panel ID of cubed-sphere coordinates
    integer, intent(in) :: np          !< Array size
    real(rp), intent(in) :: alpha(np)  !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: beta (np)  !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: gam(np)    !< A factor of r/a
    real(rp), intent(out) :: lon(np)   !< Longitude coordinate [rad]
    real(rp), intent(out) :: lat(np)   !< Latitude coordinate [rad]
 
    real(rp) :: cartpos(np,3)
    real(rp) :: geogpos(np,3)
 
    integer :: p
    !-----------------------------------------------------------------------------
 
    select case( panelid )
    case( 1, 2, 3, 4 )
      !$omp parallel do
      !$acc parallel loop present(alpha, beta, lon, lat)
      do p=1, np
        lon(p) = alpha(p) + 0.5_rp * pi * dble(panelid - 1)
        lat(p) = atan( tan( beta(p) ) * cos( alpha(p) ) )
 
        if ( lon(p) < 0.0_rp ) lon(p) = lon(p) + 2.0_rp * pi
      end do
      !$omp end parallel do
    case(5, 6)
      !$acc data create(CartPos, GeogPos) present(alpha, beta, gam, lon, lat)
      call cubedspherecoordcnv_cs2cartpos( panelid, alpha, beta, gam, np, & ! (in)
        cartpos(:,1), cartpos(:,2), cartpos(:,3)                          ) ! (out)
      
      call geographiccoordcnv_orth_to_geo_pos( cartpos(:,:), np, & ! (in)
        geogpos(:,:) ) ! (out)
      
      !$omp parallel do
      !$acc parallel loop
      do p=1, np
        lon(p) = geogpos(p,1)
        lat(p) = geogpos(p,2)
 
        if ( alpha(p) < 0.0_rp ) lon(p) = lon(p) + 2.0_rp * pi
      end do
      !$omp end parallel do
      !$acc end data
    case default
      log_error("CubedSphereCoordCnv_CS2LonLatPos",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
      call prc_abort
    end select
 
    return
  end subroutine cubedspherecoordcnv_cs2lonlatpos
 
  !> Convert the components of a vector in local coordinates with an equiangular gnomonic cubed-sphere projection to those in longitude and latitude coordinates
  !!
!OCL SERIAL
  subroutine cubedspherecoordcnv_cs2lonlatvec( &
    panelID, alpha, beta, gam, Np,         & ! (in)
    vecalpha, vecbeta,                     & ! (in)
    veclon, veclat,                        & ! (out)
    lat, gpu_async_id                      ) ! (in, optional)
 
    use scale_const, only: &
      eps => const_eps
    implicit none
 
    integer, intent(in) :: panelid       !< Panel ID of cubed-sphere coordinates
    integer, intent(in) :: np            !< Array size
    real(rp), intent(in) :: alpha(np)    !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: beta (np)    !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: gam(np)      !< A factor of RPlanet / r
    real(dp), intent(in) :: vecalpha(np) !< A component of vector in the alpha-coordinate
    real(dp), intent(in) :: vecbeta (np) !< A component of vector in the beta-coordinate
    real(rp), intent(out) :: veclon(np)  !< A component of vector in the longitude-coordinate
    real(rp), intent(out) :: veclat(np)  !< A component of vector in the latitude-coordinate
    real(rp), intent(in), optional :: lat(np) !< latitude [rad]
    integer, intent(in), optional :: gpu_async_id !< GPU asynchronous ID
 
    integer :: p
    real(rp) :: x ,y, del2
    real(rp) :: s
 
    real(rp) :: radius
    real(rp) :: cos_lat(np)
 
    integer :: gpu_async_id_
    !-----------------------------------------------------------------------------
 
#ifdef _OPENACC
    if (present(gpu_async_id)) then
      gpu_async_id_ = gpu_async_id
    else
      gpu_async_id_ = acc_async_sync
    end if
#endif
 
    !$acc data create(cos_Lat) present(alpha, beta, gam, VecLon, VecLat, VecAlpha, VecBeta)
 
    if (present(lat)) then
      !$omp parallel do
      !$acc parallel loop
      do p=1, np
        cos_lat(p) = cos(lat(p))
      end do
    end if
 
    select case( panelid )
    case( 1, 2, 3, 4 )
      if (.not. present(lat)) then
        !$omp parallel do
        !$acc parallel loop async(gpu_async_id_)
        do p=1, np
          cos_lat(p) = cos( atan( tan( beta(p) ) * cos( alpha(p) ) ) )
        end do
      end if
 
      !$omp parallel do private( X, Y, del2, radius )
      !$acc parallel loop async(gpu_async_id_)
      do p=1, np
        x = tan( alpha(p) )
        y = tan( beta(p) )
        del2 = 1.0_rp + x**2 + y**2
        radius = rplanet * gam(p)
 
        veclon(p) = vecalpha(p) * cos_lat(p) * radius
        veclat(p) = ( - x * y * vecalpha(p) + ( 1.0_rp + y**2 ) * vecbeta(p) ) &
                  * radius * sqrt( 1.0_rp + x**2 ) / del2 
      end do
    case ( 5 )
      s = 1.0_rp
    case ( 6 )
      s = - 1.0_rp
    case default
      log_error("CubedSphereCoordCnv_CS2LonLatVec",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
      call prc_abort
    end select
 
    select case( panelid )
    case( 5, 6 )
      !$omp parallel do private( X, Y, del2, radius )
      !$acc parallel loop async(gpu_async_id_)
      do p=1, np
        x = tan( alpha(p) )
        y = tan( beta(p) )
        del2 = 1.0_rp + x**2 + y**2
        radius = s * rplanet * gam(p)
 
        if (.not. present(lat)) then
          cos_lat(p) =  cos( atan( sign(1.0_rp, s) / max( sqrt( x**2 + y**2 ), eps ) ) )
        end if
 
        veclon(p) = (- y * ( 1.0_rp + x**2 ) * vecalpha(p) + x * ( 1.0_rp + y**2 ) * vecbeta(p) ) &
                  * radius / max( x**2 + y**2, eps ) * cos_lat(p)
        veclat(p) = (- x * ( 1.0_rp + x**2 ) * vecalpha(p) - y * ( 1.0_rp + y**2 ) * vecbeta(p) ) &
                  * radius / ( del2 * ( max( sqrt( x**2 + y**2 ), eps ) ) )
      end do
    end select
 
    !$acc end data
 
    return
  end subroutine cubedspherecoordcnv_cs2lonlatvec
 
  !> Calculate local coordinates using the central angles in an equiangular gnomonic cubed-sphere projection from longitude and latitude coordinates
  !!
!OCL SERIAL
  subroutine cubedspherecoordcnv_lonlat2cspos(  &
    panelID, lon, lat, Np,                 & ! (in)
    alpha, beta                            ) ! (out)
 
    implicit none
    integer, intent(in) :: panelid      !< Panel ID of cubed-sphere coordinates
    integer, intent(in) :: np           !< Array size
    real(rp), intent(in) :: lon(np)     !< Longitude coordinate [rad]
    real(rp), intent(in) :: lat(np)     !< Latitude coordinate [rad]
    real(rp), intent(out) ::alpha(np)   !< Local coordinate using the central angles [rad]
    real(rp), intent(out) ::beta (np)   !< Local coordinate using the central angles [rad]
 
    integer :: p
    real(rp) :: tan_lat
    real(rp) :: lon_(np)
    !-----------------------------------------------------------------------------
 
    !$acc data create(lon_) present(lon, lat, alpha, beta)
 
    select case(panelid)
    case ( 1, 2, 3, 4 )
      !$omp parallel
      if ( panelid == 1 ) then
        !$omp do
        !$acc parallel loop
        do p=1, np
          if ( lon(p) >= 2.0_rp * pi - 0.25_rp * pi ) then
            lon_(p) = lon(p) - 2.0_rp * pi
          else
            lon_(p) = lon(p)
          end if
        end do
      else
        !$omp do
        !$acc parallel loop
        do p=1, np
          if ( lon(p) < 0.0_rp ) then
            lon_(p) = lon(p) + 2.0_rp * pi
          else
            lon_(p) = lon(p)
          end if
        end do
      end if
      !$omp do
      !$acc parallel loop
      do p=1, np
        alpha(p) = lon_(p) - 0.5_rp * pi * ( dble(panelid) - 1.0_rp )
        beta(p) = atan( tan(lat(p)) / cos(alpha(p)) )
      end do
      !$omp end parallel
    case ( 5 )
      !$omp parallel private(tan_lat)   
      !$omp do
      !$acc parallel loop
      do p=1, np
        tan_lat = tan(lat(p)) 
        alpha(p) = + atan( sin(lon(p)) / tan_lat )
        beta(p) = - atan( cos(lon(p)) / tan_lat )
      end do
      !$omp end parallel
    case ( 6 )
      !$omp parallel private(tan_lat)    
      !$omp do
      !$acc parallel loop
      do p=1, np
        tan_lat = tan(lat(p)) 
        alpha(p) = - atan( sin(lon(p)) / tan_lat )
        beta(p) = - atan( cos(lon(p)) / tan_lat )
      end do
      !$omp end parallel
    case default
      log_error("CubedSphereCoordCnv_LonLat2CSPos",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
      call prc_abort
    end select
 
    !$acc end data
    return
  end subroutine cubedspherecoordcnv_lonlat2cspos
 
  !> Convert the components of a vector in longitude and latitude coordinates to those in local coordinates with an equiangular gnomonic cubed-sphere projection
  !!
!OCL SERIAL
  subroutine cubedspherecoordcnv_lonlat2csvec( &
    panelID, alpha, beta, gam, Np,         & ! (in)
    veclon, veclat,                        & ! (in)
    vecalpha, vecbeta,                     & ! (out)
    lat, gpu_async_id )                      ! (in, optional)
 
    implicit none
 
    integer, intent(in) :: panelid        !< Panel ID of cubed-sphere coordinates
    integer, intent(in) :: np             !< Array size
    real(rp), intent(in) :: alpha(np)     !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: beta (np)     !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: gam(np)       !< A factor of RPlanet / r
    real(dp), intent(in) :: veclon(np)    !< A component of vector in the longitude-coordinate
    real(dp), intent(in) :: veclat(np)    !< A component of vector in the latitude-coordinate
    real(dp), intent(out) :: vecalpha(np) !< A component of vector in the alpha-coordinate
    real(dp), intent(out) :: vecbeta (np) !< A component of vector in the beta-coordinate
    real(rp), intent(in), optional :: lat(np)     !< latitude [rad]
    integer, intent(in), optional :: gpu_async_id !< GPU asynchronous ID for OpenACC operations
 
    integer :: p
    real(rp) :: x ,y, del2
    real(rp) :: s
 
    real(rp) :: radius
    real(rp) :: cos_lat(np)
    real(rp) :: veclon_ov_coslat
 
    integer :: gpu_async_id_
    !-----------------------------------------------------------------------------
    
#ifdef _OPENACC
    if (present(gpu_async_id)) then
      gpu_async_id_ = gpu_async_id
    else
      gpu_async_id_ = acc_async_sync
    end if
#endif
 
    !$acc data create(cos_Lat) present(alpha, beta, gam, VecLon, VecLat, VecAlpha, VecBeta)
 
    if (present(lat)) then
      !$omp parallel do
      !$acc parallel loop async(gpu_async_id_)
      do p=1, np
        cos_lat(p) = cos(lat(p))
      end do
    end if
 
    select case( panelid )
    case( 1, 2, 3, 4 )
      if (.not. present(lat)) then
        !$omp parallel do
        !$acc parallel loop async(gpu_async_id_)
        do p=1, np
          cos_lat(p) = cos( atan( tan( beta(p) ) * cos( alpha(p) ) ) )
        end do
      end if
  
      !$omp parallel do private( X, Y, del2, radius, VecLon_ov_cosLat )
      !$acc parallel loop async(gpu_async_id_)
      do p=1, np
        x = tan( alpha(p) )
        y = tan( beta(p) )
        del2 = 1.0_rp + x**2 + y**2
        radius = rplanet * gam(p)
        veclon_ov_coslat = veclon(p) / cos_lat(p)
 
        vecalpha(p) = veclon_ov_coslat / radius
        vecbeta(p)  = ( x * y * veclon_ov_coslat + del2 / sqrt( 1.0_rp + x**2 ) * veclat(p) ) &
                     / ( radius * (1.0_rp + y**2) )
      end do
    case ( 5 )
      s = 1.0_rp
    case ( 6 )
      s = -1.0_rp
    case default
      log_error("CubedSphereCoordCnv_LonLat2CSVec",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
      call prc_abort
    end select
 
    select case( panelid )
    case( 5, 6 )
      !$omp parallel do private( X, Y, del2, radius, VecLon_ov_cosLat )
      !$acc parallel loop async(gpu_async_id_)
      do p=1, np
        x = tan( alpha(p) )
        y = tan( beta(p) )
        del2 = 1.0_rp + x**2 + y**2
        radius = s * rplanet * gam(p)
 
        if (.not. present(lat)) then
          cos_lat(p) =  cos( atan( s / sqrt(max(x**2 + y**2, eps)) ) )
        end if
        veclon_ov_coslat = veclon(p) / cos_lat(p)
 
        vecalpha(p) = (- y * veclon_ov_coslat - del2 * x / sqrt( max(del2 - 1.0_rp,eps)) * veclat(p)) &
                    / ( radius * ( 1.0_rp + x**2 ) )
        vecbeta(p) = (  x * veclon_ov_coslat - del2 * y / sqrt( max(del2 - 1.0_rp,eps)) * veclat(p)) &
                    / ( radius * ( 1.0_rp + y**2 ) )
      end do
    end select
 
    !$acc end data
 
    return
  end subroutine cubedspherecoordcnv_lonlat2csvec
 
  !> Calculate the Cartesian coordinates from local coordinates using the central angles in an equiangular gnomonic cubed-sphere projection
  !!
!OCL SERIAL
  subroutine cubedspherecoordcnv_cs2cartpos( &
    panelID, alpha, beta, gam, Np,           & ! (in)
    x, y, z                                  ) ! (out)
 
    implicit none
    integer, intent(in) :: panelid     !< Panel ID of cubed-sphere coordinates
    integer, intent(in) :: np          !< Array size
    real(rp), intent(in) :: alpha(np)  !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: beta (np)  !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: gam(np)    !< A factor of r/a
    real(rp), intent(out) :: x(np)     !< x-coordinate in the Cartesian coordinate
    real(rp), intent(out) :: y(np)     !< y-coordinate in the Cartesian coordinate
    real(rp), intent(out) :: z(np)     !< z-coordinate in the Cartesian coordinate
 
    integer :: p
    real(rp) :: x1, x2, fac
 
    !-----------------------------------------------------------------------------
 
    select case(panelid)
    case(1)
      !$omp parallel do private(x1, x2, fac)
      !$acc parallel loop present(alpha, beta, gam, X, Y, Z)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
        x(p) = fac
        y(p) = fac * x1
        z(p) = fac * x2
      end do
    case(2)
      !$omp parallel do private(x1, x2, fac)
      !$acc parallel loop present(alpha, beta, gam, X, Y, Z)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
        x(p) = - fac * x1 
        y(p) =   fac
        z(p) =   fac * x2
      end do
    case(3)
      !$omp parallel do private(x1, x2, fac)
      !$acc parallel loop present(alpha, beta, gam, X, Y, Z)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
        x(p) = - fac
        y(p) = - fac * x1 
        z(p) =   fac * x2
      end do
    case(4)
      !$omp parallel do private(x1, x2, fac)
      !$acc parallel loop present(alpha, beta, gam, X, Y, Z)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
        x(p) =   fac * x1 
        y(p) = - fac
        z(p) =   fac * x2
      end do
    case(5)
      !$omp parallel do private(x1, x2, fac)
      !$acc parallel loop present(alpha, beta, gam, X, Y, Z)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
        x(p) = - fac * x2
        y(p) =   fac * x1
        z(p) =   fac
      end do
    case(6)
      !$omp parallel do private(x1, x2, fac)
      !$acc parallel loop present(alpha, beta, gam, X, Y, Z)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        fac = rplanet * gam(p) / sqrt( 1.0_rp + x1**2 + x2**2 )
        x(p) =   fac * x2
        y(p) =   fac * x1
        z(p) = - fac
      end do
    case default
      log_error("CubedSphereCoordCnv_CS2CartPos",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
      call prc_abort
    end select
 
    return
  end subroutine cubedspherecoordcnv_cs2cartpos
 
  !> Convert the components of a vector in local coordinates with an equiangular gnomonic cubed-sphere projection to those in the Cartesian coordinates
  !!
!OCL SERIAL
  subroutine cubedspherecoordcnv_cart2csvec( &
    panelID, alpha, beta, gam, Np,         & ! (in)
    vec_x, vec_y, vec_z,                   & ! (in)
    vecalpha, vecbeta                      ) ! (out)
 
    implicit none
 
    integer, intent(in) :: panelid        !< Panel ID of cubed-sphere coordinates
    integer, intent(in) :: np             !< Array size
    real(rp), intent(in) :: alpha(np)     !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: beta (np)     !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: gam(np)       !< A factor of RPlanet / r
    real(dp), intent(in) :: vec_x(np)     !< A component of vector in the x-coordinate with the Cartesian coordinate
    real(dp), intent(in) :: vec_y(np)     !< A component of vector in the y-coordinate with the Cartesian coordinate
    real(dp), intent(in) :: vec_z(np)     !< A component of vector in the z-coordinate with the Cartesian coordinate
    real(rp), intent(out) :: vecalpha(np) !< A component of vector in the alpha-coordinate
    real(rp), intent(out) :: vecbeta (np) !< A component of vector in the beta-coordinate
 
    integer :: p
    real(rp) :: x1, x2, fac
    real(rp) :: r_sec2_alpha, r_sec2_beta
    !-----------------------------------------------------------------------------
 
    select case( panelid )
    case(1)
      !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
      !$acc parallel loop present(alpha, beta, gam, Vec_x, Vec_y, Vec_z, VecAlpha, VecBeta)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        r_sec2_alpha = cos(alpha(p))**2
        r_sec2_beta  = cos(beta(p))**2
        fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
 
        vecalpha(p) = fac * r_sec2_alpha * ( - x1  * vec_x(p) + vec_y(p) )
        vecbeta(p) = fac * r_sec2_beta  * ( - x2  * vec_x(p) + vec_z(p) )
      end do      
    case(2)
      !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
      !$acc parallel loop present(alpha, beta, gam, Vec_x, Vec_y, Vec_z, VecAlpha, VecBeta)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        r_sec2_alpha = cos(alpha(p))**2
        r_sec2_beta  = cos(beta(p))**2
        fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
 
        vecalpha(p) = - fac * r_sec2_alpha * ( vec_x(p) + x1 * vec_y(p) )
        vecbeta(p) =   fac * r_sec2_beta  * ( - x2  * vec_y(p) + vec_z(p) )
      end do
    case(3)
      !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
      !$acc parallel loop present(alpha, beta, gam, Vec_x, Vec_y, Vec_z, VecAlpha, VecBeta)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        r_sec2_alpha = cos(alpha(p))**2
        r_sec2_beta  = cos(beta(p))**2
        fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
 
        vecalpha(p) = fac * r_sec2_alpha * ( x1  * vec_x(p) - vec_y(p) )
        vecbeta(p) = fac * r_sec2_beta  * ( x2  * vec_x(p) + vec_z(p) )
      end do
    case(4)
      !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
      !$acc parallel loop present(alpha, beta, gam, Vec_x, Vec_y, Vec_z, VecAlpha, VecBeta)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        r_sec2_alpha = cos(alpha(p))**2
        r_sec2_beta  = cos(beta(p))**2
        fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
 
        vecalpha(p) = fac * r_sec2_alpha * ( vec_x(p) + x1 * vec_y(p) )
        vecbeta(p) = fac * r_sec2_beta  * ( x2 * vec_y(p) + vec_z(p) )
      end do
    case ( 5 )
      !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
      !$acc parallel loop present(alpha, beta, gam, Vec_x, Vec_y, Vec_z, VecAlpha, VecBeta)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        r_sec2_alpha = cos(alpha(p))**2
        r_sec2_beta  = cos(beta(p))**2
        fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
 
        vecalpha(p) =   fac * r_sec2_alpha * ( vec_y(p) - x1 * vec_z(p) )
        vecbeta(p) = - fac * r_sec2_beta  * ( vec_x(p) + x2 * vec_z(p) )
      end do
    case ( 6 )
      !$omp parallel do private(x1, x2, r_sec2_alpha, r_sec2_beta, fac)
      !$acc parallel loop present(alpha, beta, gam, Vec_x, Vec_y, Vec_z, VecAlpha, VecBeta)
      do p=1, np
        x1 = tan( alpha(p) )
        x2 = tan( beta(p) )
        r_sec2_alpha = cos(alpha(p))**2
        r_sec2_beta  = cos(beta(p))**2
        fac = sqrt( 1.0_rp + x1**2 + x2**2 ) / ( rplanet * gam(p) )
 
        vecalpha(p) = fac * r_sec2_alpha * ( vec_y(p) + x1 * vec_z(p) )
        vecbeta(p) = fac * r_sec2_beta  * ( vec_x(p) + x2 * vec_z(p) )
      end do
    case default
      log_error("CubedSphereCoordCnv_Cart2CSVec",'(a,i2,a)') "panelID ", panelid, " is invalid. Check!"
      call prc_abort
    end select
 
    return
  end subroutine cubedspherecoordcnv_cart2csvec
 
!OCL SERIAL  
  subroutine cubedspherecoordcnv_cs2localorthvec_alpha_0( &
    alpha, beta, radius, Np,                      & ! (in)
    vecorth1, vecorth2                            ) ! (inout)
 
    implicit none
    
    integer, intent(in) :: Np
    real(RP), intent(in) :: alpha(Np)
    real(RP), intent(in) :: beta (Np)
    real(RP), intent(in) :: radius(Np)
    real(RP), intent(inout) :: VecOrth1(Np)
    real(RP), intent(inout) :: VecOrth2(Np)
 
    integer :: p
    real(RP) :: x1, x2, del
    real(RP) :: fac
    real(RP) :: tmp
    !-----------------------------------------------------------------------------
 
    !$omp parallel do private(x1, x2, del, fac, tmp)
    !$acc parallel loop present(alpha, beta, radius, VecOrth1, VecOrth2)
    do p=1, np
      x1 = tan( alpha(p) )
      x2 = tan( beta(p) )
      del = sqrt(1.0_rp + x1**2 + x2**2)
      fac = radius(p) * sqrt(1.0_rp + x1**2) / del**2
 
      tmp = vecorth1(p)
      vecorth1(p) = fac * del * tmp
      vecorth2(p) = fac * ( - x1 * x2 * tmp + (1.0_rp + x2**2) * vecorth2(p) )
    end do
 
    return
  end subroutine cubedspherecoordcnv_cs2localorthvec_alpha_0
 
!OCL SERIAL  
  subroutine cubedspherecoordcnv_cs2localorthvec_alpha_1( &
    alpha, beta, radius, Np,                      & ! (in)
    vecalpha, vecbeta,                            & ! (in)
    vecorth1, vecorth2                            ) ! (out)
 
    implicit none
    
    integer, intent(in) :: Np
    real(RP), intent(in) :: alpha(Np)
    real(RP), intent(in) :: beta (Np)
    real(RP), intent(in) :: radius(Np)
    real(RP), intent(in) :: VecAlpha(Np)
    real(RP), intent(in) :: VecBeta (Np)
    real(RP), intent(out) :: VecOrth1(Np)
    real(RP), intent(out) :: VecOrth2(Np)
 
    integer :: p
    !-----------------------------------------------------------------------------
 
    !$omp parallel do
    !$acc parallel loop present(VecAlpha, VecBeta, VecOrth1, VecOrth2)
    do p=1, np
      vecorth1(p) = vecalpha(p)
      vecorth2(p) = vecbeta(p)
    end do
    !$omp end parallel do
 
    call cubedspherecoordcnv_cs2localorthvec_alpha_0( &
      alpha, beta, radius, np,                      & ! (in)
      vecorth1, vecorth2                            ) ! (inout)    
 
    return
  end subroutine cubedspherecoordcnv_cs2localorthvec_alpha_1
 
!OCL SERIAL  
  subroutine cubedspherecoordcnv_localorth2csvec_alpha_0( &
    alpha, beta, radius, Np,                      & ! (in)
    vecalpha, vecbeta                             ) ! (inout)
 
    implicit none
    
    integer, intent(in) :: Np
    real(RP), intent(in) :: alpha(Np)
    real(RP), intent(in) :: beta (Np)
    real(RP), intent(in) :: radius(Np)
    real(RP), intent(inout) :: VecAlpha(Np)
    real(RP), intent(inout) :: VecBeta (Np)    
 
    integer :: p
    real(RP) :: x1, x2, del
    real(RP) :: fac
    real(RP) :: tmp
    !-----------------------------------------------------------------------------
 
    !$omp parallel do private(x1, x2, del, fac, tmp)
    !$acc parallel loop present(alpha, beta, radius, VecAlpha, VecBeta)
    do p=1, np
      x1 = tan( alpha(p) )
      x2 = tan( beta(p) )
      del = sqrt(1.0_rp + x1**2 + x2**2)
      fac = del / ( radius(p) * ( 1.0_rp + x2**2 ) * sqrt( 1.0_rp + x1**2 ) )
 
      tmp = vecalpha(p)
      vecalpha(p) = fac * ( 1.0_rp + x2**2 ) * tmp
      vecbeta(p) = fac * ( x1 * x2 * tmp + del * vecbeta(p) )
    end do
 
    return
  end subroutine cubedspherecoordcnv_localorth2csvec_alpha_0
 
!OCL SERIAL  
  subroutine cubedspherecoordcnv_localorth2csvec_alpha_1( &
    alpha, beta, radius, Np,                      & ! (in)
    vecorth1, vecorth2,                           & ! (in)
    vecalpha, vecbeta                             ) ! (out)
    implicit none
    integer, intent(in) :: Np
    real(RP), intent(in) :: alpha(Np)
    real(RP), intent(in) :: beta (Np)
    real(RP), intent(in) :: radius(Np)
    real(RP), intent(in) :: VecOrth1(Np)
    real(RP), intent(in) :: VecOrth2(Np)
    real(RP), intent(out) :: VecAlpha(Np)
    real(RP), intent(out) :: VecBeta (Np)    
 
    integer :: p
    !-----------------------------------------------------------------------------
 
    !$omp parallel do
    !$acc parallel loop present(VecOrth1, VecOrth2, VecAlpha, VecBeta)
    do p=1, np
      vecalpha(p) = vecorth1(p)
      vecbeta(p) = vecorth2(p)
    end do
    !$omp end parallel do
 
    call cubedspherecoordcnv_localorth2csvec_alpha_0( &
      alpha, beta, radius, np,                        & ! (in)
      vecalpha, vecbeta                               ) ! (inout)    
 
    return
  end subroutine cubedspherecoordcnv_localorth2csvec_alpha_1
 
  !> Calculate the metrics associated with an equiangular gnomonic cubed-sphere projection to those in longitude and latitude coordinates
  !!
!OCL SERIAL  
  subroutine cubedspherecoordcnv_getmetric( &
    alpha, beta, Np, radius,            & ! (in)
    g_ij, gij, gsqrt                    ) ! (out)
 
    implicit none
 
    integer, intent(in) :: np             !< Array size
    real(rp), intent(in) :: alpha(np)     !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: beta (np)     !< Local coordinate using the central angles [rad]
    real(rp), intent(in) :: radius        !< Planetary radius
    real(rp), intent(out) :: g_ij(np,2,2) !< Horizontal covariant metric tensor
    real(rp), intent(out) :: gij (np,2,2) !< Horizontal contravariant metric tensor
    real(rp), intent(out) :: gsqrt(np)    !< Horizontal Jacobian
 
    real(rp) :: x, y
    real(rp) :: r2
    real(rp) :: fac
    real(rp) :: g_ij_(2,2)
    real(rp) :: gij_ (2,2)
 
    integer :: p
    real(rp) :: oneplusx2, oneplusy2
    !-----------------------------------------------------------------------------
 
    !$omp parallel do private( &
    !$omp X, Y, r2, OnePlusX2, OnePlusY2, fac,    &
    !$omp G_ij_, GIJ_                             )
    !$acc parallel loop private(G_ij_, GIJ_) present(alpha, beta, G_ij, GIJ, Gsqrt)
    do p=1, np
      x = tan(alpha(p))
      y = tan(beta(p))
      r2 = 1.0_rp + x**2 + y**2
      oneplusx2 = 1.0_rp + x**2
      oneplusy2 = 1.0_rp + y**2
 
      fac = oneplusx2 * oneplusy2 * ( radius / r2 )**2
      g_ij_(1,1) =   fac * oneplusx2
      g_ij_(1,2) = - fac * (x * y)
      g_ij_(2,1) = - fac * (x * y)
      g_ij_(2,2) =   fac * oneplusy2
      g_ij(p,:,:) = g_ij_(:,:)
 
      gsqrt(p) = radius**2 * oneplusx2 * oneplusy2 / ( r2 * sqrt(r2) )
 
      fac = 1.0_rp / gsqrt(p)**2
      gij_(1,1) =   fac * g_ij_(2,2)
      gij_(1,2) = - fac * g_ij_(1,2)
      gij_(2,1) = - fac * g_ij_(2,1)
      gij_(2,2) =   fac * g_ij_(1,1)
      gij(p,:,:) = gij_(:,:)
    end do
 
    return
  end subroutine cubedspherecoordcnv_getmetric
 
end module scale_cubedsphere_coord_cnv