scale_mapprojection.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_mapprojection
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: mapprojection_setup
  public :: mapprojection_xy2lonlat
  public :: mapprojection_lonlat2xy
  public :: mapprojection_mapfactor
  public :: mapprojection_rotcoef
  public :: mapprojection_get_attributes

  interface mapprojection_xy2lonlat
     module procedure mapprojection_xy2lonlat_0d
     module procedure mapprojection_xy2lonlat_2d
  end interface
  interface mapprojection_lonlat2xy
     module procedure mapprojection_lonlat2xy_0d
     module procedure mapprojection_lonlat2xy_2d
  end interface
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public :: mapprojection_basepoint_lon = 135.221_rp ! position of base point (domain center) in real world [deg]
  real(RP), public :: mapprojection_basepoint_lat =  34.653_rp ! position of base point (domain center) in real world [deg]

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: mapprojection_none_setup
  private :: mapprojection_lambertconformal_setup
  private :: mapprojection_polarstereographic_setup
  private :: mapprojection_mercator_setup
  private :: mapprojection_equidistantcylindrical_setup

  private :: mapprojection_none_xy2lonlat
  private :: mapprojection_lambertconformal_xy2lonlat
  private :: mapprojection_polarstereographic_xy2lonlat
  private :: mapprojection_mercator_xy2lonlat
  private :: mapprojection_equidistantcylindrical_xy2lonlat

  private :: mapprojection_none_lonlat2xy
  private :: mapprojection_lambertconformal_lonlat2xy
  private :: mapprojection_polarstereographic_lonlat2xy
  private :: mapprojection_mercator_lonlat2xy
  private :: mapprojection_equidistantcylindrical_lonlat2xy

  private :: mapprojection_none_mapfactor
  private :: mapprojection_lambertconformal_mapfactor
  private :: mapprojection_polarstereographic_mapfactor
  private :: mapprojection_mercator_mapfactor
  private :: mapprojection_equidistantcylindrical_mapfactor

  private :: mapprojection_none_rotcoef
  private :: mapprojection_lambertconformal_rotcoef
  private :: mapprojection_polarstereographic_rotcoef
  private :: mapprojection_mercator_rotcoef
  private :: mapprojection_equidistantcylindrical_rotcoef

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  character(len=H_SHORT), private :: mapprojection_type = 'NONE'
                                               ! 'NONE'
                                               ! 'LC'
                                               ! 'PS'
                                               ! 'MER'
                                               ! 'EC'

  real(DP), private :: mapprojection_hemisphere         ! hemisphere flag: 1=north, -1=south

  real(DP), private :: mapprojection_basepoint_x        ! position of base point in the model [m]
  real(DP), private :: mapprojection_basepoint_y        ! position of base point in the model [m]

  real(DP), private :: mapprojection_pole_x             ! position of north/south pole in the model [m]
  real(DP), private :: mapprojection_pole_y             ! position of north/south pole in the model [m]
  real(DP), private :: mapprojection_eq_x               ! position of equator at the base lon. in the model [m]
  real(DP), private :: mapprojection_eq_y               ! position of equator at the base lon. in the model [m]

  real(DP), private :: mapprojection_rotation =  0.0_dp ! rotation factor
  real(DP), private :: mapprojection_rot_fact_sin
  real(DP), private :: mapprojection_rot_fact_cos

  real(DP), private :: mapprojection_lc_lat1  = 30.0_dp ! standard latitude1 for L.C. projection [deg]
  real(DP), private :: mapprojection_lc_lat2  = 60.0_dp ! standard latitude2 for L.C. projection [deg]
  real(DP), private :: mapprojection_lc_c               ! conformal factor
  real(DP), private :: mapprojection_lc_fact            ! pre-calc factor

  real(DP), private :: mapprojection_ps_lat             ! standard latitude1 for P.S. projection [deg]
  real(DP), private :: mapprojection_ps_fact            ! pre-calc factor

  real(DP), private :: mapprojection_m_lat    =  0.0_dp ! standard latitude1 for Mer. projection [deg]
  real(DP), private :: mapprojection_m_fact             ! pre-calc factor

  real(DP), private :: mapprojection_ec_lat   =  0.0_dp ! standard latitude1 for E.C. projection [deg]
  real(DP), private :: mapprojection_ec_fact            ! pre-calc factor

  real(DP), private :: pi
  real(DP), private :: d2r
  real(DP), private :: radius

  character(len=H_SHORT) :: mapprojection_mapping = ""
  real(DP), private :: mapprojection_false_easting
  real(DP), private :: mapprojection_false_northing
  real(DP), private :: mapprojection_longitude_of_central_meridian
  real(DP), private :: mapprojection_longitude_of_projection_origin
  real(DP), private :: mapprojection_latitude_of_projection_origin
  real(DP), private :: mapprojection_straight_vertical_longitude_from_pole
  real(DP), private :: mapprojection_standard_parallel(2)

  interface
     subroutine xy2lonlat_s( x, y, lon, lat )
       use scale_precision
       real(DP), intent(in)  :: x, y
       real(RP), intent(out) :: lon, lat
     end subroutine xy2lonlat_s
     subroutine lonlat2xy_s( lon, lat, x, y )
       use scale_precision
       real(RP), intent(in)  :: lon, lat
       real(DP), intent(out) :: x, y
     end subroutine lonlat2xy_s
     subroutine mapfactor_s( lat, m1, m2 )
       use scale_precision
       real(RP), intent(in)  :: lat
       real(RP), intent(out) :: m1, m2
     end subroutine mapfactor_s
     subroutine rotcoef_s( lon, lat, rotc_cos, rotc_sin )
       use scale_precision
       real(RP), intent(in)  :: lon, lat
       real(RP), intent(out) :: rotc_cos, rotc_sin
     end subroutine rotcoef_s
  end interface
  procedure(xy2lonlat_s), pointer :: xy2lonlat => null()
  procedure(lonlat2xy_s), pointer :: lonlat2xy => null()
  procedure(mapfactor_s), pointer :: mapfactor => null()
  procedure(rotcoef_s),   pointer :: rotcoef   => null()

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine mapprojection_setup( DOMAIN_CENTER_X, DOMAIN_CENTER_Y )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       undef     => const_undef,  &
       pi_rp     => const_pi,     &
       d2r_rp    => const_d2r,    &
       radius_rp => const_radius
    implicit none

    real(RP), intent(in) :: DOMAIN_CENTER_X
    real(RP), intent(in) :: DOMAIN_CENTER_Y

    namelist / param_mapprojection / &
       mapprojection_basepoint_lon, &
       mapprojection_basepoint_lat, &
       mapprojection_basepoint_x,   &
       mapprojection_basepoint_y,   &
       mapprojection_type,          &
       mapprojection_rotation,      &
       mapprojection_lc_lat1,       &
       mapprojection_lc_lat2,       &
       mapprojection_ps_lat,        &
       mapprojection_m_lat,         &
       mapprojection_ec_lat

    integer :: ierr
    !---------------------------------------------------------------------------

    log_newline
    log_info("MAPPROJECTION_setup",*) 'Setup'

    pi     = real(pi_rp,     kind=dp)
    d2r    = real(d2r_rp,    kind=dp)
    radius = real(radius_rp, kind=dp)

    mapprojection_basepoint_x = undef
    mapprojection_basepoint_y = undef
    mapprojection_ps_lat      = undef

    mapprojection_basepoint_x = domain_center_x
    mapprojection_basepoint_y = domain_center_y

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_mapprojection,iostat=ierr)

    if( ierr < 0 ) then !--- missing
       log_info("MAPPROJECTION_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("MAPPROJECTION_setup",*) 'Not appropriate names in namelist PARAM_MAPPROJECTION. Check!'
       call prc_abort
    endif
    log_nml(param_mapprojection)

    log_newline
    log_info("MAPPROJECTION_setup",*) 'Map projection information'
    log_info_cont('(1x,A,F15.3)') 'Basepoint(x)    [m] : ', mapprojection_basepoint_x
    log_info_cont('(1x,A,F15.3)') 'Basepoint(y)    [m] : ', mapprojection_basepoint_y
    log_info_cont(*)              'Map projection type : ', trim(mapprojection_type)

    mapprojection_mapping = ""
    mapprojection_false_easting = undef
    mapprojection_false_northing = undef
    mapprojection_longitude_of_central_meridian = undef
    mapprojection_longitude_of_projection_origin = undef
    mapprojection_latitude_of_projection_origin = undef
    mapprojection_straight_vertical_longitude_from_pole = undef
    mapprojection_standard_parallel(:) = undef

    if( mapprojection_ps_lat == undef ) mapprojection_ps_lat = mapprojection_basepoint_lat

    select case(trim(mapprojection_type))
    case('NONE')
       log_info_cont(*) '=> NO map projection'
       call mapprojection_none_setup
       xy2lonlat => mapprojection_none_xy2lonlat
       lonlat2xy => mapprojection_none_lonlat2xy
       mapfactor => mapprojection_none_mapfactor
       rotcoef   => mapprojection_none_rotcoef
    case('LC')
       log_info_cont(*) '=> Lambert Conformal projection'
       call mapprojection_lambertconformal_setup
       xy2lonlat => mapprojection_lambertconformal_xy2lonlat
       lonlat2xy => mapprojection_lambertconformal_lonlat2xy
       mapfactor => mapprojection_lambertconformal_mapfactor
       rotcoef   => mapprojection_lambertconformal_rotcoef
    case('PS')
       log_info_cont(*) '=> Polar Stereographic projection'
       call mapprojection_polarstereographic_setup
       xy2lonlat => mapprojection_polarstereographic_xy2lonlat
       lonlat2xy => mapprojection_polarstereographic_lonlat2xy
       mapfactor => mapprojection_polarstereographic_mapfactor
       rotcoef   => mapprojection_polarstereographic_rotcoef
    case('MER')
       log_info_cont(*) '=> Mercator projection'
       call mapprojection_mercator_setup
       xy2lonlat => mapprojection_mercator_xy2lonlat
       lonlat2xy => mapprojection_mercator_lonlat2xy
       mapfactor => mapprojection_mercator_mapfactor
       rotcoef   => mapprojection_mercator_rotcoef
    case('EC')
       log_info_cont(*) '=> Equidistant Cylindrical projection'
       call mapprojection_equidistantcylindrical_setup
       xy2lonlat => mapprojection_equidistantcylindrical_xy2lonlat
       lonlat2xy => mapprojection_equidistantcylindrical_lonlat2xy
       mapfactor => mapprojection_equidistantcylindrical_mapfactor
       rotcoef   => mapprojection_equidistantcylindrical_rotcoef
    case default
       log_error("MAPPROJECTION_setup",*) 'Unsupported MAPPROJECTION_type. STOP'
       call prc_abort
    endselect

    mapprojection_rotation = mapprojection_rotation * d2r
    mapprojection_rot_fact_sin = sin(mapprojection_rotation)
    mapprojection_rot_fact_cos = cos(mapprojection_rotation)

    return
  end subroutine mapprojection_setup

  !-----------------------------------------------------------------------------
  subroutine mapprojection_xy2lonlat_0d( &
       x, y,    &
       lon, lat )
    implicit none
    real(RP), intent(in)  :: x
    real(RP), intent(in)  :: y
    real(RP), intent(out) :: lon ! [rad]
    real(RP), intent(out) :: lat ! [rad]

    real(DP) :: xx, yy
    !---------------------------------------------------------------------------

    xx = mapprojection_basepoint_x &
       + ( x - mapprojection_basepoint_x ) * mapprojection_rot_fact_cos &
       - ( y - mapprojection_basepoint_y ) * mapprojection_rot_fact_sin
    yy = mapprojection_basepoint_y &
       + ( y - mapprojection_basepoint_y ) * mapprojection_rot_fact_cos &
       + ( x - mapprojection_basepoint_x ) * mapprojection_rot_fact_sin

    call xy2lonlat( xx, yy, lon, lat )

    return
  end subroutine mapprojection_xy2lonlat_0d

  subroutine mapprojection_xy2lonlat_2d( &
       IA, IS, IE, JA, JS, JE, &
       x, y,    &
       lon, lat )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(in)  :: x(ia,ja)
    real(RP), intent(in)  :: y(ia,ja)
    real(RP), intent(out) :: lon(ia,ja) ! [rad]
    real(RP), intent(out) :: lat(ia,ja) ! [rad]

    integer :: i, j
    !---------------------------------------------------------------------------

    !$omp parallel do
    do j = js, je
    do i = is, ie
       call mapprojection_xy2lonlat_0d( x(i,j), y(i,j), lon(i,j), lat(i,j) )
    end do
    end do

    return
  end subroutine mapprojection_xy2lonlat_2d

  !-----------------------------------------------------------------------------
  subroutine mapprojection_lonlat2xy_0d( &
       lon, lat, &
       x, y      )
    implicit none
    real(RP), intent(in)  :: lon ! [rad]
    real(RP), intent(in)  :: lat ! [rad]
    real(RP), intent(out) :: x
    real(RP), intent(out) :: y

    real(DP) :: xx, yy
    !---------------------------------------------------------------------------

    call lonlat2xy( lon, lat, xx, yy )

    xx = xx - mapprojection_basepoint_x
    yy = yy - mapprojection_basepoint_y

    x = mapprojection_basepoint_x + ( xx * mapprojection_rot_fact_cos &
                                    + yy * mapprojection_rot_fact_sin )
    y = mapprojection_basepoint_y + ( yy * mapprojection_rot_fact_cos &
                                    - xx * mapprojection_rot_fact_sin )

    return
  end subroutine mapprojection_lonlat2xy_0d

  subroutine mapprojection_lonlat2xy_2d( &
       IA, IS, IE, JA, JS, JE, &
       lon, lat, &
       x, y      )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(in)  :: lon(ia,ja) ! [rad]
    real(RP), intent(in)  :: lat(ia,ja) ! [rad]
    real(RP), intent(out) :: x(ia,ja)
    real(RP), intent(out) :: y(ia,ja)

    integer  :: i, j
    !---------------------------------------------------------------------------

    !$omp parallel do
    do j = js, je
    do i = is, ie
       call mapprojection_lonlat2xy_0d( lon(i,j), lat(i,j), x(i,j), y(i,j) )
    end do
    end do

    return
  end subroutine mapprojection_lonlat2xy_2d

  !-----------------------------------------------------------------------------
  subroutine mapprojection_mapfactor( &
       IA, IS, IE, JA, JS, JE, &
       lat,   &
       m1, m2 )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(in)  :: lat(ia,ja) ! [rad]
    real(RP), intent(out) :: m1 (ia,ja)
    real(RP), intent(out) :: m2 (ia,ja)

    !real(RP) :: mm1, mm2
    integer :: i, j
    !---------------------------------------------------------------------------

    !$omp parallel do
    do j = js, je
    do i = is, ie
       call mapfactor( lat(i,j), m1(i,j), m2(i,j) )
!!$       call mapfactor( lat(i,j), mm1, mm2 )
!!$       m1(i,j) = sqrt( (MAPPROJECTION_rot_fact_cos * mm1)**2 + (MAPPROJECTION_rot_fact_sin * mm2)**2 )
!!$       m2(i,j) = sqrt( (MAPPROJECTION_rot_fact_cos * mm2)**2 + (MAPPROJECTION_rot_fact_sin * mm1)**2 )
    end do
    end do

    return
  end subroutine mapprojection_mapfactor

  !-----------------------------------------------------------------------------
  subroutine mapprojection_rotcoef( &
       IA, IS, IE, JA, JS, JE, &
       lon, lat, &
       rotc      )
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE

    real(RP), intent(in)  :: lon(ia,ja)    ! [rad]
    real(RP), intent(in)  :: lat(ia,ja)    ! [rad]
    real(RP), intent(out) :: rotc(ia,ja,2)

    integer :: i, j
    !---------------------------------------------------------------------------

    !$omp parallel do
    do j = js, je
    do i = is, ie
       call rotcoef( lon(i,j), lat(i,j), rotc(i,j,1), rotc(i,j,2) )
    end do
    end do

    return
  end subroutine mapprojection_rotcoef

  !-----------------------------------------------------------------------------
  subroutine mapprojection_get_attributes( &
       mapping,                               &
       false_easting,                         &
       false_northing,                        &
       longitude_of_central_meridian,         &
       longitude_of_projection_origin,        &
       latitude_of_projection_origin,         &
       straight_vertical_longitude_from_pole, &
       standard_parallel                      )
    implicit none

    character(len=*), intent(out) :: mapping

    real(DP), intent(out), optional :: false_easting
    real(DP), intent(out), optional :: false_northing
    real(DP), intent(out), optional :: longitude_of_central_meridian
    real(DP), intent(out), optional :: longitude_of_projection_origin
    real(DP), intent(out), optional :: latitude_of_projection_origin
    real(DP), intent(out), optional :: straight_vertical_longitude_from_pole
    real(DP), intent(out), optional :: standard_parallel(2)
    !---------------------------------------------------------------------------

    mapping = mapprojection_mapping

    if ( present(false_easting)                         ) &
                 false_easting                          = mapprojection_false_easting
    if ( present(false_northing)                        ) &
                 false_northing                         = mapprojection_false_northing
    if ( present(longitude_of_central_meridian)         ) &
                 longitude_of_central_meridian          = mapprojection_longitude_of_central_meridian
    if ( present(longitude_of_projection_origin)        ) &
                 longitude_of_projection_origin         = mapprojection_longitude_of_projection_origin
    if ( present(latitude_of_projection_origin)         ) &
                 latitude_of_projection_origin          = mapprojection_latitude_of_projection_origin
    if ( present(straight_vertical_longitude_from_pole) ) &
                 straight_vertical_longitude_from_pole  = mapprojection_straight_vertical_longitude_from_pole
    if ( present(standard_parallel)                     ) &
                 standard_parallel(:)                   = mapprojection_standard_parallel(:)

    return
  end subroutine mapprojection_get_attributes

  !-----------------------------------------------------------------------------
  subroutine mapprojection_none_setup
    implicit none
    !---------------------------------------------------------------------------

    log_newline
    log_info("MAPPROJECTION_None_setup",*) 'MAPPROJECTION_rotation   = ', mapprojection_rotation

    mapprojection_mapping = ""

    return
  end subroutine mapprojection_none_setup

  !-----------------------------------------------------------------------------
  subroutine mapprojection_none_xy2lonlat( &
       x,   &
       y,   &
       lon, &
       lat  )
    implicit none

    real(DP), intent(in)  :: x
    real(DP), intent(in)  :: y
    real(RP), intent(out) :: lon ! [rad]
    real(RP), intent(out) :: lat ! [rad]

    real(DP) :: gno(2)
    real(DP) :: rho, gmm
    !---------------------------------------------------------------------------

    gno(1) = ( x - mapprojection_basepoint_x ) / radius
    gno(2) = ( y - mapprojection_basepoint_y ) / radius

    rho = sqrt( gno(1) * gno(1) + gno(2) * gno(2) )
    gmm = atan( rho )

    if ( rho == 0.0_dp ) then
       lon = mapprojection_basepoint_lon * d2r
       lat = mapprojection_basepoint_lat * d2r
    else
       lon = mapprojection_basepoint_lon * d2r &
           + atan( gno(1)*sin(gmm) / ( rho   *cos(mapprojection_basepoint_lat*d2r)*cos(gmm) &
                                     - gno(2)*sin(mapprojection_basepoint_lat*d2r)*sin(gmm) ) )
       lat = asin(        sin(mapprojection_basepoint_lat*d2r)*cos(gmm)       &
                 + gno(2)*cos(mapprojection_basepoint_lat*d2r)*sin(gmm) / rho )
    endif

    return
  end subroutine mapprojection_none_xy2lonlat

  !-----------------------------------------------------------------------------
  subroutine mapprojection_none_lonlat2xy( &
       lon, &
       lat, &
       x,   &
       y    )
    use scale_prc, only: &
       prc_abort
    implicit none

    real(RP), intent(in)  :: lon ! [rad]
    real(RP), intent(in)  :: lat ! [rad]
    real(DP), intent(out) :: x
    real(DP), intent(out) :: y

    real(DP) :: lat_d, lat0_d, dlon
    real(DP) :: gno(2)
    real(DP) :: cos_gmm
    !---------------------------------------------------------------------------
    ! http://mathworld.wolfram.com/GnomonicProjection.html

    lat_d  = real(lat,kind=dp)
    lat0_d = real(mapprojection_basepoint_lat*d2r,kind=dp)

    dlon = lon - mapprojection_basepoint_lon * d2r

    cos_gmm = sin(lat0_d) * sin(lat_d) &
            + cos(lat0_d) * cos(lat_d) * cos(dlon)

    gno(1) = ( cos(lat_d)  * sin(dlon)  ) / cos_gmm
    gno(2) = ( cos(lat0_d) * sin(lat_d) &
             - sin(lat0_d) * cos(lat_d) * cos(dlon) ) / cos_gmm

    x = gno(1) * radius + mapprojection_basepoint_x
    y = gno(2) * radius + mapprojection_basepoint_y

    return
  end subroutine mapprojection_none_lonlat2xy

  !-----------------------------------------------------------------------------
  subroutine mapprojection_none_mapfactor( &
       lat,   &
       m1, m2 )
    implicit none

    real(RP), intent(in)  :: lat ! [rad]
    real(RP), intent(out) :: m1
    real(RP), intent(out) :: m2
    !---------------------------------------------------------------------------

    m1 = 1.0_rp
    m2 = m1

    return
  end subroutine mapprojection_none_mapfactor

  !-----------------------------------------------------------------------------
  subroutine mapprojection_none_rotcoef( &
       lon, lat,          &
       rotc_cos, rotc_sin )
    implicit none
    real(RP), intent(in)  :: lon, lat
    real(RP), intent(out) :: rotc_cos, rotc_sin
    !---------------------------------------------------------------------------

    rotc_cos = mapprojection_rot_fact_cos
    rotc_sin = mapprojection_rot_fact_sin

    return
  end subroutine mapprojection_none_rotcoef

  !-----------------------------------------------------------------------------
  subroutine mapprojection_lambertconformal_setup
    use scale_prc, only: &
       prc_abort
    implicit none

    real(DP) :: lat1rot, lat2rot
    real(DP) :: dlon, latrot, dist
    !---------------------------------------------------------------------------

    if ( mapprojection_lc_lat1 >= mapprojection_lc_lat2 ) then
       log_error("MAPPROJECTION_LambertConformal_setup",*) 'Please set MAPPROJECTION_LC_lat1 < MAPPROJECTION_LC_lat2 in degree. STOP'
       call prc_abort
    endif

    ! check hemisphere: 1=north, -1=south
    mapprojection_hemisphere = sign(1.0_dp,mapprojection_lc_lat1+mapprojection_lc_lat2)

    lat1rot = 0.5_dp*pi - mapprojection_hemisphere * mapprojection_lc_lat1 * d2r
    lat2rot = 0.5_dp*pi - mapprojection_hemisphere * mapprojection_lc_lat2 * d2r

    ! calc conformal factor c
    mapprojection_lc_c = ( log( sin(lat1rot)        ) - log( sin(lat2rot)        ) ) &
              / ( log( tan(0.5_dp*lat1rot) ) - log( tan(0.5_dp*lat2rot) ) )

    ! pre-calc factor
    mapprojection_lc_fact = sin(lat1rot) / mapprojection_lc_c / tan(0.5_dp*lat1rot)**mapprojection_lc_c

    ! calc (x,y) at pole point
    dlon = 0.0_dp

    latrot = 0.5_dp*pi - mapprojection_hemisphere * mapprojection_basepoint_lat * d2r

    dist = mapprojection_lc_fact * radius * tan(0.5_dp*latrot)**mapprojection_lc_c

    mapprojection_pole_x = mapprojection_basepoint_x -                   dist * sin(mapprojection_lc_c*dlon)
    mapprojection_pole_y = mapprojection_basepoint_y + mapprojection_hemisphere * dist * cos(mapprojection_lc_c*dlon)

    log_newline
    log_info("MAPPROJECTION_LambertConformal_setup",*) 'Input parameters'
    log_info_cont(*) 'MAPPROJECTION_LC_lat1    = ', mapprojection_lc_lat1
    log_info_cont(*) 'MAPPROJECTION_LC_lat2    = ', mapprojection_lc_lat2
    log_info_cont(*) 'MAPPROJECTION_hemisphere = ', mapprojection_hemisphere
    log_info_cont(*) 'MAPPROJECTION_LC_c       = ', mapprojection_lc_c
    log_info_cont(*) 'MAPPROJECTION_LC_fact    = ', mapprojection_lc_fact
    log_info_cont(*) 'MAPPROJECTION_pole_x     = ', mapprojection_pole_x
    log_info_cont(*) 'MAPPROJECTION_pole_y     = ', mapprojection_pole_y
    log_info_cont(*) 'MAPPROJECTION_pole_y     = ', mapprojection_pole_y
    log_info_cont(*) 'MAPPROJECTION_rotation   = ', mapprojection_rotation

    mapprojection_mapping = "lambert_conformal_conic"
    mapprojection_standard_parallel(:) = (/ mapprojection_lc_lat1, mapprojection_lc_lat2 /)
    mapprojection_longitude_of_central_meridian = mapprojection_basepoint_lon
    mapprojection_latitude_of_projection_origin = mapprojection_basepoint_lat
    mapprojection_false_easting  = mapprojection_basepoint_x
    mapprojection_false_northing = mapprojection_basepoint_y

    return
  end subroutine mapprojection_lambertconformal_setup

  !-----------------------------------------------------------------------------
  subroutine mapprojection_lambertconformal_xy2lonlat( &
       x,   &
       y,   &
       lon, &
       lat  )
    implicit none

    real(DP), intent(in)  :: x
    real(DP), intent(in)  :: y
    real(RP), intent(out) :: lon ! [rad]
    real(RP), intent(out) :: lat ! [rad]

    real(DP) :: xx, yy, dist
    !---------------------------------------------------------------------------

    xx =                             ( x - mapprojection_pole_x ) / radius / mapprojection_lc_fact
    yy = -mapprojection_hemisphere * ( y - mapprojection_pole_y ) / radius / mapprojection_lc_fact

    dist = sqrt( xx*xx + yy*yy )

    lon = mapprojection_basepoint_lon * d2r + atan2(xx,yy) / mapprojection_lc_c

    ! check hemisphere: 1=north, -1=south
    lat = mapprojection_hemisphere * ( 0.5_dp*pi - 2.0_dp*atan( dist**(1.0_dp/mapprojection_lc_c) ) )

    return
  end subroutine mapprojection_lambertconformal_xy2lonlat

  !-----------------------------------------------------------------------------
  subroutine mapprojection_lambertconformal_lonlat2xy( &
       lon, &
       lat, &
       x,   &
       y    )
    implicit none

    real(RP), intent(in)  :: lon ! [rad]
    real(RP), intent(in)  :: lat ! [rad]
    real(DP), intent(out) :: x
    real(DP), intent(out) :: y

    real(DP) :: dlon, latrot, dist
    !---------------------------------------------------------------------------

    dlon = lon - mapprojection_basepoint_lon * d2r
    if ( dlon >  pi ) dlon = dlon - pi*2.0_dp
    if ( dlon < -pi ) dlon = dlon + pi*2.0_dp

    latrot = 0.5_dp*pi - mapprojection_hemisphere * lat

    dist = mapprojection_lc_fact * radius * tan(0.5_dp*latrot)**mapprojection_lc_c

    x = mapprojection_pole_x +                   dist * sin(mapprojection_lc_c*dlon)
    y = mapprojection_pole_y - mapprojection_hemisphere * dist * cos(mapprojection_lc_c*dlon)

    return
  end subroutine mapprojection_lambertconformal_lonlat2xy

  !-----------------------------------------------------------------------------
  subroutine mapprojection_lambertconformal_mapfactor( &
       lat,   &
       m1, m2 )
    implicit none

    real(RP), intent(in)  :: lat ! [rad]
    real(RP), intent(out) :: m1
    real(RP), intent(out) :: m2

    real(DP) :: latrot
    !---------------------------------------------------------------------------

    latrot = 0.5_dp*pi - mapprojection_hemisphere * lat

    m1 = mapprojection_lc_fact / sin(latrot) * mapprojection_lc_c * tan(0.5_dp*latrot)**mapprojection_lc_c
    m2 = m1

    return
  end subroutine mapprojection_lambertconformal_mapfactor

  !-----------------------------------------------------------------------------
  subroutine mapprojection_lambertconformal_rotcoef( &
       lon, lat,          &
       rotc_cos, rotc_sin )
    implicit none
    real(RP), intent(in)  :: lon, lat
    real(RP), intent(out) :: rotc_cos, rotc_sin

    real(DP) :: dlon
    real(DP) :: alpha
    !---------------------------------------------------------------------------

    dlon = lon - mapprojection_basepoint_lon * d2r
    if ( dlon >  pi ) dlon = dlon - pi*2.0_dp
    if ( dlon < -pi ) dlon = dlon + pi*2.0_dp

    alpha = - mapprojection_lc_c * dlon * mapprojection_hemisphere &
            + mapprojection_rotation

    rotc_cos = cos( alpha )
    rotc_sin = sin( alpha )

    return
  end subroutine mapprojection_lambertconformal_rotcoef

  !-----------------------------------------------------------------------------
  subroutine mapprojection_polarstereographic_setup
    implicit none

    real(DP) :: lat0
    real(DP) :: dlon, latrot, dist
    !---------------------------------------------------------------------------

    ! check hemisphere: 1=north, -1=south
    mapprojection_hemisphere = sign(1.0_dp,mapprojection_ps_lat)

    lat0 = mapprojection_hemisphere * mapprojection_ps_lat * d2r

    ! pre-calc factor
    mapprojection_ps_fact = 1.0_dp + sin(lat0)

    ! calc (x,y) at pole point
    dlon = 0.0_dp

    latrot = 0.5_dp*pi - mapprojection_hemisphere * mapprojection_basepoint_lat * d2r

    dist = mapprojection_ps_fact * radius * tan(0.5_dp*latrot)

    mapprojection_pole_x = mapprojection_basepoint_x -                   dist * sin(dlon)
    mapprojection_pole_y = mapprojection_basepoint_y + mapprojection_hemisphere * dist * cos(dlon)

    log_newline
    log_info("MAPPROJECTION_PolarStereographic_setup",*) 'MAPPROJECTION_PS_lat1    = ', mapprojection_ps_lat
    log_info("MAPPROJECTION_PolarStereographic_setup",*) 'MAPPROJECTION_hemisphere = ', mapprojection_hemisphere
    log_info("MAPPROJECTION_PolarStereographic_setup",*) 'MAPPROJECTION_PS_fact    = ', mapprojection_ps_fact
    log_info("MAPPROJECTION_PolarStereographic_setup",*) 'MAPPROJECTION_pole_x     = ', mapprojection_pole_x
    log_info("MAPPROJECTION_PolarStereographic_setup",*) 'MAPPROJECTION_pole_y     = ', mapprojection_pole_y

    mapprojection_mapping = "polar_stereographic"
    mapprojection_straight_vertical_longitude_from_pole = mapprojection_basepoint_lon
    mapprojection_latitude_of_projection_origin = mapprojection_basepoint_lat
    mapprojection_standard_parallel(1) = mapprojection_ps_lat
    mapprojection_false_easting = mapprojection_basepoint_x
    mapprojection_false_northing = mapprojection_basepoint_y

    return
  end subroutine mapprojection_polarstereographic_setup

  !-----------------------------------------------------------------------------
  subroutine mapprojection_polarstereographic_xy2lonlat( &
       x,   &
       y,   &
       lon, &
       lat  )
    implicit none

    real(DP), intent(in)  :: x
    real(DP), intent(in)  :: y
    real(RP), intent(out) :: lon ! [rad]
    real(RP), intent(out) :: lat ! [rad]

    real(DP) :: xx, yy, dist
    !---------------------------------------------------------------------------

    xx =                             ( x - mapprojection_pole_x ) / radius / mapprojection_ps_fact
    yy = -mapprojection_hemisphere * ( y - mapprojection_pole_y ) / radius / mapprojection_ps_fact

    dist = sqrt( xx*xx + yy*yy )

    lon = mapprojection_basepoint_lon * d2r + atan2(xx,yy)

    lat = mapprojection_hemisphere * ( 0.5_dp*pi - 2.0_dp*atan(dist) )

    return
  end subroutine mapprojection_polarstereographic_xy2lonlat

  !-----------------------------------------------------------------------------
  subroutine mapprojection_polarstereographic_lonlat2xy( &
       lon, &
       lat, &
       x,   &
       y    )
    implicit none

    real(RP), intent(in)  :: lon ! [rad]
    real(RP), intent(in)  :: lat ! [rad]
    real(DP), intent(out) :: x
    real(DP), intent(out) :: y

    real(DP) :: dlon, latrot, dist
    !---------------------------------------------------------------------------

    dlon = lon - mapprojection_basepoint_lon * d2r

    latrot = 0.5_dp*pi - mapprojection_hemisphere * lat

    dist = mapprojection_ps_fact * radius * tan(0.5_dp*latrot)

    x = mapprojection_pole_x +                   dist * sin(dlon)
    y = mapprojection_pole_y - mapprojection_hemisphere * dist * cos(dlon)

    return
  end subroutine mapprojection_polarstereographic_lonlat2xy

  !-----------------------------------------------------------------------------
  subroutine mapprojection_polarstereographic_mapfactor( &
       lat,   &
       m1, m2 )
    implicit none

    real(RP), intent(in)  :: lat ! [rad]
    real(RP), intent(out) :: m1
    real(RP), intent(out) :: m2
    !---------------------------------------------------------------------------

    m1 = mapprojection_ps_fact / ( 1.0_dp + sin(mapprojection_hemisphere*lat) )
    m2 = m1

    return
  end subroutine mapprojection_polarstereographic_mapfactor

  !-----------------------------------------------------------------------------
  subroutine mapprojection_polarstereographic_rotcoef( &
       lon, lat,          &
       rotc_cos, rotc_sin )
    implicit none
    real(RP), intent(in)  :: lon, lat
    real(RP), intent(out) :: rotc_cos, rotc_sin

    real(DP) :: dlon
    real(DP) :: alpha
    !---------------------------------------------------------------------------

    dlon = lon - mapprojection_basepoint_lon * d2r
    if ( dlon >  pi ) dlon = dlon - pi*2.0_dp
    if ( dlon < -pi ) dlon = dlon + pi*2.0_dp

    alpha = - dlon * mapprojection_hemisphere &
            + mapprojection_rotation

    rotc_cos = cos( alpha )
    rotc_sin = sin( alpha )

    return
  end subroutine mapprojection_polarstereographic_rotcoef

  !-----------------------------------------------------------------------------
  subroutine mapprojection_mercator_setup
    use scale_prc, only: &
       prc_abort
    implicit none

    real(DP) :: lat0
    real(DP) :: latrot, dist
    !---------------------------------------------------------------------------

    lat0 = mapprojection_m_lat * d2r

    ! pre-calc factor
    mapprojection_m_fact = cos(lat0)

    if ( mapprojection_m_fact == 0.0_dp ) then
       log_error("MAPPROJECTION_Mercator_setup",*) 'MAPPROJECTION_M_lat cannot be set to pole point! value=', mapprojection_m_lat
       call prc_abort
    endif

    ! calc (x,y) at (lon,lat) = (base,0)
    latrot = 0.5_dp*pi - mapprojection_basepoint_lat * d2r

    dist = 1.0_dp / tan(0.5_dp*latrot)

    mapprojection_eq_x = mapprojection_basepoint_x
    mapprojection_eq_y = mapprojection_basepoint_y - radius * mapprojection_m_fact * log(dist)

    log_newline
    log_info("MAPPROJECTION_Mercator_setup",*) 'MAPPROJECTION_M_lat      = ', mapprojection_m_lat
    log_info("MAPPROJECTION_Mercator_setup",*) 'MAPPROJECTION_M_fact     = ', mapprojection_m_fact
    log_info("MAPPROJECTION_Mercator_setup",*) 'MAPPROJECTION_eq_x       = ', mapprojection_eq_x
    log_info("MAPPROJECTION_Mercator_setup",*) 'MAPPROJECTION_eq_y       = ', mapprojection_eq_y

    mapprojection_mapping = "mercator"
    mapprojection_longitude_of_projection_origin = mapprojection_basepoint_lon
    mapprojection_standard_parallel(1) = mapprojection_m_lat
    mapprojection_false_easting = mapprojection_basepoint_x
    mapprojection_false_northing = mapprojection_basepoint_y

    return
  end subroutine mapprojection_mercator_setup

  !-----------------------------------------------------------------------------
  subroutine mapprojection_mercator_xy2lonlat( &
       x,   &
       y,   &
       lon, &
       lat  )
    implicit none

    real(DP), intent(in)  :: x
    real(DP), intent(in)  :: y
    real(RP), intent(out) :: lon ! [rad]
    real(RP), intent(out) :: lat ! [rad]

    real(DP) :: xx, yy
    !---------------------------------------------------------------------------

    xx = ( x - mapprojection_eq_x ) / radius / mapprojection_m_fact
    yy = ( y - mapprojection_eq_y ) / radius / mapprojection_m_fact

    lon = xx + mapprojection_basepoint_lon * d2r

    lat = 0.5_dp*pi - 2.0_dp*atan( 1.0_dp/exp(yy) )

    return
  end subroutine mapprojection_mercator_xy2lonlat

  !-----------------------------------------------------------------------------
  subroutine mapprojection_mercator_lonlat2xy( &
       lon, &
       lat, &
       x,   &
       y    )
    implicit none

    real(RP), intent(in)  :: lon ! [rad]
    real(RP), intent(in)  :: lat ! [rad]
    real(DP), intent(out) :: x
    real(DP), intent(out) :: y

    real(DP) :: dlon, latrot, dist
    !---------------------------------------------------------------------------

    dlon = lon - mapprojection_basepoint_lon * d2r

    latrot = 0.5_dp*pi - lat

    dist = 1.0_dp / tan(0.5_dp*latrot)

    x = mapprojection_eq_x + radius * mapprojection_m_fact * dlon
    y = mapprojection_eq_y + radius * mapprojection_m_fact * log(dist)

    return
  end subroutine mapprojection_mercator_lonlat2xy

  !-----------------------------------------------------------------------------
  subroutine mapprojection_mercator_mapfactor( &
       lat,   &
       m1, m2 )
    implicit none

    real(RP), intent(in)  :: lat ! [rad]
    real(RP), intent(out) :: m1
    real(RP), intent(out) :: m2
    !---------------------------------------------------------------------------

    m1 = mapprojection_m_fact / cos(real(lat,kind=dp))
    m2 = m1

    return
  end subroutine mapprojection_mercator_mapfactor

  !-----------------------------------------------------------------------------
  subroutine mapprojection_mercator_rotcoef( &
       lon, lat,          &
       rotc_cos, rotc_sin )
    implicit none
    real(RP), intent(in)  :: lon, lat
    real(RP), intent(out) :: rotc_cos, rotc_sin
    !---------------------------------------------------------------------------

    rotc_cos = mapprojection_rot_fact_cos
    rotc_sin = mapprojection_rot_fact_sin

    return
  end subroutine mapprojection_mercator_rotcoef

  !-----------------------------------------------------------------------------
  subroutine mapprojection_equidistantcylindrical_setup
    use scale_prc, only: &
       prc_abort
    implicit none

    real(DP) :: lat0
    !---------------------------------------------------------------------------

    lat0 = mapprojection_ec_lat * d2r

    ! pre-calc factor
    mapprojection_ec_fact = cos(lat0)

    if ( mapprojection_ec_fact == 0.0_dp ) then
       log_error("MAPPROJECTION_EquidistantCylindrical_setup",*) 'MAPPROJECTION_EC_lat cannot be set to pole point! value=', mapprojection_ec_lat
       call prc_abort
    endif

    mapprojection_eq_x = mapprojection_basepoint_x
    mapprojection_eq_y = mapprojection_basepoint_y - radius * mapprojection_basepoint_lat * d2r

    log_newline
    log_info("MAPPROJECTION_EquidistantCylindrical_setup",*) 'MAPPROJECTION_EC_lat     = ', mapprojection_ec_lat
    log_info("MAPPROJECTION_EquidistantCylindrical_setup",*) 'MAPPROJECTION_EC_fact    = ', mapprojection_ec_fact
    log_info("MAPPROJECTION_EquidistantCylindrical_setup",*) 'MAPPROJECTION_eq_x       = ', mapprojection_eq_x
    log_info("MAPPROJECTION_EquidistantCylindrical_setup",*) 'MAPPROJECTION_eq_y       = ', mapprojection_eq_y

    mapprojection_mapping = "equirectangular"
    mapprojection_standard_parallel(1) = mapprojection_ec_lat
    mapprojection_longitude_of_central_meridian = mapprojection_basepoint_lon
    mapprojection_false_easting  = mapprojection_basepoint_x
    mapprojection_false_northing = mapprojection_basepoint_y

    return
  end subroutine mapprojection_equidistantcylindrical_setup

  !-----------------------------------------------------------------------------
  subroutine mapprojection_equidistantcylindrical_xy2lonlat( &
       x,   &
       y,   &
       lon, &
       lat  )
    use scale_prc, only: &
       prc_abort
    implicit none

    real(DP), intent(in)  :: x
    real(DP), intent(in)  :: y
    real(RP), intent(out) :: lon ! [rad]
    real(RP), intent(out) :: lat ! [rad]

    real(DP) :: xx, yy
    !---------------------------------------------------------------------------

    xx = ( x - mapprojection_eq_x ) / radius / mapprojection_ec_fact
    yy = ( y - mapprojection_eq_y ) / radius

    lon = xx + mapprojection_basepoint_lon * d2r

    lat = yy

    if ( abs(lat) >  0.5_dp*pi ) then
       log_error("MAPPROJECTION_EquidistantCylindrical_xy2lonlat",*) 'Invalid latitude range! value=', lat
       call prc_abort
    endif

    return
  end subroutine mapprojection_equidistantcylindrical_xy2lonlat

  !-----------------------------------------------------------------------------
  subroutine mapprojection_equidistantcylindrical_lonlat2xy( &
       lon, &
       lat, &
       x,   &
       y    )
    implicit none

    real(RP), intent(in)  :: lon ! [rad]
    real(RP), intent(in)  :: lat ! [rad]
    real(DP), intent(out) :: x
    real(DP), intent(out) :: y

    real(DP) :: dlon
    !---------------------------------------------------------------------------

    dlon = lon - mapprojection_basepoint_lon * d2r

    x = mapprojection_eq_x + radius * mapprojection_ec_fact * dlon
    y = mapprojection_eq_y + radius * lat

    return
  end subroutine mapprojection_equidistantcylindrical_lonlat2xy

  !-----------------------------------------------------------------------------
  subroutine mapprojection_equidistantcylindrical_mapfactor( &
       lat,   &
       m1, m2 )
    implicit none

    real(RP), intent(in)  :: lat ! [rad]
    real(RP), intent(out) :: m1
    real(RP), intent(out) :: m2
    real(DP) :: mm1, mm2
    !---------------------------------------------------------------------------

    mm1 = mapprojection_ec_fact / cos(real(lat,kind=dp))
    mm2 = 1.0_rp

    m1 = sqrt( (mapprojection_rot_fact_cos * mm1)**2 + (mapprojection_rot_fact_sin * mm2)**2 )
    m2 = sqrt( (mapprojection_rot_fact_cos * mm2)**2 + (mapprojection_rot_fact_sin * mm1)**2 )

    return
  end subroutine mapprojection_equidistantcylindrical_mapfactor

  !-----------------------------------------------------------------------------
  subroutine mapprojection_equidistantcylindrical_rotcoef( &
       lon, lat,          &
       rotc_cos, rotc_sin )
    implicit none
    real(RP), intent(in)  :: lon, lat
    real(RP), intent(out) :: rotc_cos, rotc_sin
    !---------------------------------------------------------------------------

    rotc_cos = mapprojection_rot_fact_cos
    rotc_sin = mapprojection_rot_fact_sin

    return
  end subroutine mapprojection_equidistantcylindrical_rotcoef

end module scale_mapprojection