scale_geographic_coord_cnv.F90

Go to the documentation of this file.

 
#include "scaleFElib.h"
module scale_geographic_coord_cnv
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_const, only: &
    pi => const_pi,      &
    eps => const_eps
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prc
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
  public :: geographiccoordcnv_orth_to_geo_pos
  public :: geographiccoordcnv_geo_to_orth_pos
 
  public :: geographiccoordcnv_orth_to_geo_vec
  public :: geographiccoordcnv_geo_to_orth_vec
 
  public :: geographiccoordcnv_rotatex
  public :: geographiccoordcnv_rotatey
  public :: geographiccoordcnv_rotatez
 
contains
!OCL SERIAL
  subroutine geographiccoordcnv_orth_to_geo_pos( orth_p, Np, &
    geo_p )
    
    implicit none
    integer, intent(in) :: np
    real(rp), intent(in) :: orth_p(np,3)
    real(rp), intent(out) :: geo_p(np,3)
 
    integer :: p
    !---------------------------------------------------------  
 
    !$omp parallel private(p)
    !$omp do
    do p=1, np
      geo_p(p,3) = sqrt( orth_p(p,1)**2 + orth_p(p,2)**2 + orth_p(p,3)**2 )
      geo_p(p,2) = asin( orth_p(p,3) / geo_p(p,3) )
      geo_p(p,1) = atan( orth_p(p,2) / orth_p(p,1) )
    end do
    !$omp do
    do p=1, np
      if ( geo_p(p,1) <= 0.0_rp .and. orth_p(p,1) < 0.0_rp ) then
        geo_p(p,1) = geo_p(p,1) + pi
      else if ( geo_p(p,1) >= 0.0_rp .and. orth_p(p,1) < 0.0_rp ) then
        geo_p(p,1) = geo_p(p,1) - pi
      end if
      if ( orth_p(p,1) == 0.0_rp .and. orth_p(p,2) == 0.0_rp ) then
        geo_p(p,1) = 0.0_rp
      else if ( orth_p(p,1) == 0.0_rp ) then
        geo_p(p,1) = sign(1.0_rp, orth_p(p,2)) * 0.5_rp * pi
      end if   
    end do
    !$omp end parallel
 
    return
  end subroutine geographiccoordcnv_orth_to_geo_pos
 
!OCL SERIAL
  subroutine geographiccoordcnv_geo_to_orth_pos( geo_p, Np, &
    orth_p )
 
    implicit none
    integer, intent(in) :: np
    real(rp), intent(in) :: geo_p(np,3)
    real(rp), intent(out) :: orth_p(np,3)
 
    integer :: p
    !---------------------------------------------------------  
 
    !$omp parallel do private(p)
    do p=1, np
      orth_p(p,1) = geo_p(p,3) * cos(geo_p(p,2)) * cos(geo_p(p,1))
      orth_p(p,2) = geo_p(p,3) * cos(geo_p(p,2)) * sin(geo_p(p,1))
      orth_p(p,3) = geo_p(p,3) * sin(geo_p(p,2))
    end do
 
    return
  end subroutine geographiccoordcnv_geo_to_orth_pos
 
!OCL SERIAL
  subroutine geographiccoordcnv_orth_to_geo_vec( orth_v, geo_p, Np, &
    geo_v ) 
    implicit none
 
    integer, intent(in) :: np
    real(rp), intent(in) :: orth_v(np,3)
    real(rp), intent(in) :: geo_p(np,3)
    real(rp), intent(out) :: geo_v(np,3)
 
    integer :: p
    real(rp) :: sin_geo_p1
    real(rp) :: cos_geo_p1
    real(rp) :: sin_geo_p2
    real(rp) :: cos_geo_p2    
    !---------------------------------------------------------  
 
    !$omp parallel do private(p, cos_geo_p2)
    do p=1, np
      sin_geo_p1 = sin(geo_p(p,1))
      cos_geo_p1 = cos(geo_p(p,1))
      sin_geo_p2 = sin(geo_p(p,2))
      cos_geo_p2 = cos(geo_p(p,2))
 
      geo_v(p,3) =   orth_v(p,1) * cos_geo_p2 * cos_geo_p1  &
                   + orth_v(p,2) * cos_geo_p2 * sin_geo_p1 &
                   + orth_v(p,3) * sin_geo_p2
  
      geo_v(p,2) = - orth_v(p,1) * sin_geo_p2 * cos_geo_p1 &
                   - orth_v(p,2) * sin_geo_p2 * sin_geo_p1 &
                   + orth_v(p,3) * cos_geo_p2
    
      geo_v(p,1) = - orth_v(p,1) * sin_geo_p1 &
                   + orth_v(p,2) * cos_geo_p1
    end do
 
    return
  end subroutine geographiccoordcnv_orth_to_geo_vec
 
!OCL SERIAL
  subroutine geographiccoordcnv_geo_to_orth_vec( geo_v, geo_p, Np, &
    orth_v ) 
    implicit none
 
    integer, intent(in) :: np
    real(rp), intent(in) :: geo_v(np,3)
    real(rp), intent(in) :: geo_p(np,3)
    real(rp), intent(out) :: orth_v(np,3)
 
    integer :: p
    real(rp) :: sin_geo_p1
    real(rp) :: cos_geo_p1
    real(rp) :: sin_geo_p2
    real(rp) :: cos_geo_p2    
    !---------------------------------------------------------  
 
    !$omp parallel do private(p, sin_geo_p1, cos_geo_p1, sin_geo_p2, cos_geo_p2)
    do p=1, np
      sin_geo_p1 = sin(geo_p(p,1))
      cos_geo_p1 = cos(geo_p(p,1))
      sin_geo_p2 = sin(geo_p(p,2))
      cos_geo_p2 = cos(geo_p(p,2))
 
      orth_v(p,1) =   geo_v(p,3) * cos_geo_p1 * cos_geo_p2 &
                    - geo_v(p,2) * cos_geo_p1 * sin_geo_p2 &
                    - geo_v(p,1) * sin_geo_p1
 
      orth_v(p,2) =   geo_v(p,3) * sin_geo_p1 * cos_geo_p2 &
                    - geo_v(p,2) * sin_geo_p1 * sin_geo_p2 &
                    + geo_v(p,1) * cos_geo_p1
 
      orth_v(p,3) = geo_v(p,3) * sin_geo_p2 + geo_v(p,2) * cos_geo_p2                  
    end do
 
    return
  end subroutine geographiccoordcnv_geo_to_orth_vec
 
!OCL SERIAL
  subroutine geographiccoordcnv_rotatex( pos_vec, angle, &
    rotated_pos_vec )
 
    implicit none
    real(rp), intent(in) :: pos_vec(3)
    real(rp), intent(in) :: angle
    real(rp), intent(out) :: rotated_pos_vec(3)
 
    real(rp) :: sin_angle, cos_angle
    !---------------------------------------------------------  
 
    sin_angle = sin(angle)
    cos_angle = cos(angle)
 
    rotated_pos_vec(1) = pos_vec(1)
    rotated_pos_vec(2) = cos_angle * pos_vec(2) - sin_angle * pos_vec(3)
    rotated_pos_vec(3) = sin_angle * pos_vec(2) + cos_angle * pos_vec(3)
 
    return
  end subroutine geographiccoordcnv_rotatex
 
!OCL SERIAL
  subroutine geographiccoordcnv_rotatey( pos_vec, angle, &
    rotated_pos_vec )
    
    implicit none
    real(rp), intent(in) :: pos_vec(3)
    real(rp), intent(in) :: angle
    real(rp), intent(out) :: rotated_pos_vec(3)
 
    real(rp) :: sin_angle, cos_angle
    !---------------------------------------------------------  
 
    sin_angle = sin(angle)
    cos_angle = cos(angle)
 
    rotated_pos_vec(1) = cos_angle * pos_vec(1) + sin_angle * pos_vec(3)
    rotated_pos_vec(2) = pos_vec(2)
    rotated_pos_vec(3) = - sin_angle * pos_vec(1) + cos_angle * pos_vec(3)
 
    return
  end subroutine geographiccoordcnv_rotatey
 
!OCL SERIAL
  subroutine geographiccoordcnv_rotatez( pos_vec, angle, &
    rotated_pos_vec )
 
    implicit none
    real(rp), intent(in) :: pos_vec(3)
    real(rp), intent(in) :: angle
    real(rp), intent(out) :: rotated_pos_vec(3)
 
    real(rp) :: sin_angle, cos_angle
    !---------------------------------------------------------  
 
    sin_angle = sin(angle)
    cos_angle = cos(angle)
 
    rotated_pos_vec(1) = cos_angle * pos_vec(1) - sin_angle * pos_vec(2)
    rotated_pos_vec(2) = sin_angle * pos_vec(1) + cos_angle * pos_vec(2)
    rotated_pos_vec(3) = pos_vec(3)
 
    return
  end subroutine geographiccoordcnv_rotatez
  
end  module scale_geographic_coord_cnv