! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! + +
! + glimmer_map_trans.f90 - part of the Glimmer-CISM ice model +
! + +
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009, 2010
! Glimmer-CISM contributors - see AUTHORS file for list of contributors
!
! This file is part of Glimmer-CISM.
!
! Glimmer-CISM is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 2 of the License, or (at
! your option) any later version.
!
! Glimmer-CISM is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with Glimmer-CISM. If not, see .
!
! Glimmer-CISM is hosted on BerliOS.de:
! https://developer.berlios.de/projects/glimmer-cism/
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#ifdef HAVE_CONFIG_H
#include "config.inc"
#endif
!> convert between projections
module glimmer_map_trans
use glimmer_map_types
implicit none
private
public :: glimmap_ll_to_xy, glimmap_xy_to_ll, loncorrect
contains
!> Convert lat-long coordinates to grid coordinates.
!!
!! The subroutine returns the x-y coordinates as real values,
!! non-integer values indicating a position between grid-points.
subroutine glimmap_ll_to_xy(lon,lat,x,y,proj,grid)
use glimmer_log
use glimmer_coordinates
implicit none
real(rk),intent(in) :: lon !< The location of the point in lat-lon space (Longitude)
real(rk),intent(in) :: lat !< The location of the point in lat-lon space (Latitude)
real(rk),intent(out) :: x !< The location of the point in x-y space (x coordinate)
real(rk),intent(out) :: y !< The location of the point in x-y space (y coordinate)
type(glimmap_proj), intent(in) :: proj !< The projection being used
type(coordsystem_type),intent(in) :: grid !< the grid definition
real(rk) :: xx,yy ! These are real-space distances in meters
if (associated(proj%laea)) then
call glimmap_laea(lon,lat,xx,yy,proj%laea)
else if (associated(proj%aea)) then
call glimmap_aea(lon,lat,xx,yy,proj%aea)
else if (associated(proj%lcc)) then
call glimmap_lcc(lon,lat,xx,yy,proj%lcc)
else if (associated(proj%stere)) then
call glimmap_stere(lon,lat,xx,yy,proj%stere)
else
call write_log('No known projection found!',GM_WARNING)
end if
! Now convert the real-space distances to grid-points using the grid type
call space2grid(xx,yy,x,y,grid)
end subroutine glimmap_ll_to_xy
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Convert grid coordinates to lat-lon coordinates.
!!
!! The subroutine returns the lat-lon coordinates as real values,
!! non-integer values indicating a position between grid-points.
subroutine glimmap_xy_to_ll(lon,lat,x,y,proj,grid)
use glimmer_log
use glimmer_coordinates
implicit none
real(rk),intent(out) :: lon !< The location of the point in lat-lon space (Longitude)
real(rk),intent(out) :: lat !< The location of the point in lat-lon space (Latitude)
real(rk),intent(in) :: x !< The location of the point in x-y space (x coordinate)
real(rk),intent(in) :: y !< The location of the point in x-y space (y coordinate)
type(glimmap_proj), intent(in) :: proj !< The projection being used
type(coordsystem_type),intent(in) :: grid !< the grid definition
real(rk) :: xx,yy ! These are real-space distances in meters
! First convert grid-point space to real space
call grid2space(xx,yy,x,y,grid)
if (associated(proj%laea)) then
call glimmap_ilaea(lon,lat,xx,yy,proj%laea)
else if (associated(proj%aea)) then
call glimmap_iaea(lon,lat,xx,yy,proj%aea)
else if (associated(proj%lcc)) then
call glimmap_ilcc(lon,lat,xx,yy,proj%lcc)
else if (associated(proj%stere)) then
call glimmap_istere(lon,lat,xx,yy,proj%stere)
else
call write_log('No known projection found!',GM_WARNING)
end if
lon=loncorrect(lon,0.0_rk)
end subroutine glimmap_xy_to_ll
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! PRIVATE subroutines follow
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Lambert azimuthal equal area projection
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Forward transformation: lat-lon -> x-y of Lambert azimuthal equal area projection
subroutine glimmap_laea(lon,lat,x,y,params)
use glimmer_log
real(rk),intent(in) :: lon !< longitude
real(rk),intent(in) :: lat !< latitude
real(rk),intent(out) :: x !< x
real(rk),intent(out) :: y !< y
type(proj_laea),intent(in) :: params !< projection parameters
real(rk) :: sin_lat,cos_lat,sin_lon,cos_lon,c,dlon,dlat,tmp,k
character(80) :: errtxt
dlon = lon-params%longitude_of_central_meridian
! Check domain of longitude
dlon = loncorrect(dlon,-180.0_rk)
! Convert to radians and calculate sine and cos
dlon = dlon*D2R ; dlat = lat*D2R
call sincos(dlon,sin_lon,cos_lon);
call sincos(dlat,sin_lat,cos_lat);
c = cos_lat * cos_lon
! Mapping transformation
tmp = 1.0 + params%sinp * sin_lat + params%cosp * c
if (tmp > 0.0) then
k = EQ_RAD * sqrt (2.0 / tmp)
x = k * cos_lat * sin_lon
y = k * (params%cosp * sin_lat - params%sinp * c)
else
write(errtxt,*)'LAEA projection error:',lon,lat,params%latitude_of_projection_origin
call write_log(trim(errtxt),GM_FATAL,__FILE__,__LINE__)
endif
! Apply false eastings and northings
x = x + params%false_easting
y = y + params%false_northing
end subroutine glimmap_laea
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Inverse transformation: lat-lon -> x-y of Lambert azimuthal equal area projection
subroutine glimmap_ilaea(lon,lat,x,y,params)
use glimmer_log
real(rk),intent(out) :: lon !< longitude
real(rk),intent(out) :: lat !< latitude
real(rk),intent(in) :: x !< x
real(rk),intent(in) :: y !< y
type(proj_laea),intent(in) :: params !< projection parameters
real(rk) :: rho,c,sin_c,cos_c,xx,yy
character(80) :: errtxt
xx=x ; yy=y
! Account for false eastings and northings
xx = xx - params%false_easting
yy = yy - params%false_northing
rho=hypot (xx,yy)
if (abs(rho) < CONV_LIMIT) then
! If very near the centre of the map...
lat = params%latitude_of_projection_origin
lon = params%longitude_of_central_meridian
else
c = 2.0 * asin(0.5 * rho * i_EQ_RAD)
call sincos (c, sin_c, cos_c)
lat = asin (cos_c * params%sinp + (yy * sin_c * params%cosp / rho)) * R2D
select case(params%pole)
case(1)
lon = params%longitude_of_central_meridian + R2D * atan2 (xx, -yy)
case(-1)
lon = params%longitude_of_central_meridian + R2D * atan2 (xx, yy)
case(0)
lon = params%longitude_of_central_meridian + &
R2D * atan2 (xx * sin_c, (rho * params%cosp * cos_c - yy * params%sinp * sin_c))
case default
write(errtxt,*)'Inverse LAEA projection error:',params%pole
call write_log(trim(errtxt),GM_FATAL,__FILE__,__LINE__)
end select
endif
end subroutine glimmap_ilaea
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Albers equal area conic projection
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Forward transformation: lat-lon -> x-y of Albers equal area conic projection
subroutine glimmap_aea(lon,lat,x,y,params)
real(rk),intent(in) :: lon !< longitude
real(rk),intent(in) :: lat !< latitude
real(rk),intent(out) :: x !< x
real(rk),intent(out) :: y !< y
type(proj_aea),intent(in) :: params !< projection parameters
real(rk) :: dlon,theta,sint,cost,rho
dlon = lon-params%longitude_of_central_meridian
! Check domain of longitude
dlon = loncorrect(dlon,-180.0_rk)
theta = params%n * dlon * D2R
call sincos(theta,sint,cost)
rho = params%i_n*sqrt(params%c - 2.0*params%n*sin(lat*D2R))
x = EQ_RAD * rho * sint
y = EQ_RAD * (params%rho0_R - rho * cost)
! Apply false eastings and northings
x = x + params%false_easting
y = y + params%false_northing
end subroutine glimmap_aea
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Inverse transformation: lat-lon -> x-y of Albers equal area conic projection
subroutine glimmap_iaea(lon,lat,x,y,params)
real(rk),intent(out) :: lon !< longitude
real(rk),intent(out) :: lat !< latitude
real(rk),intent(in) :: x !< x
real(rk),intent(in) :: y !< y
type(proj_aea),intent(in) :: params !< projection parameters
real(rk) :: xx,yy,rho,theta
xx=x ; yy=y
! Account for false eastings and northings
xx = xx - params%false_easting
yy = yy - params%false_northing
rho = sqrt(xx**2.0 + (params%rho0 - yy)**2.0)
if (params%n >0.0) then
theta = atan2(xx,(params%rho0-yy))
else
theta = atan2(-xx,(yy-params%rho0))
end if
lat = asin((params%c-(rho*params%n/EQ_RAD)**2.0)*0.5*params%i_n)*R2D
lon = params%longitude_of_central_meridian+R2D*theta*params%i_n
end subroutine glimmap_iaea
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Lambert conformal conic projection
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Forward transformation: lat-lon -> x-y of Lambert conformal conic projection
subroutine glimmap_lcc(lon,lat,x,y,params)
real(rk),intent(in) :: lon !< longitude
real(rk),intent(in) :: lat !< latitude
real(rk),intent(out) :: x !< x
real(rk),intent(out) :: y !< y
type(proj_lcc),intent(in) :: params !< projection parameters
real(rk) :: dlon,rho,theta,sint,cost
dlon = lon-params%longitude_of_central_meridian
! Check domain of longitude
dlon = loncorrect(dlon,-180.0_rk)
rho = EQ_RAD * params%f/(tan(M_PI_4+lat*D2R/2.0))**params%n
theta = params%n*dlon*D2R
call sincos(theta,sint,cost)
x = rho * sint
y = params%rho0 - rho * cost
! Apply false eastings and northings
x = x + params%false_easting
y = y + params%false_northing
end subroutine glimmap_lcc
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Inverse transformation: lat-lon -> x-y of Lambert conformal conic projection
subroutine glimmap_ilcc(lon,lat,x,y,params)
real(rk),intent(out) :: lon !< longitude
real(rk),intent(out) :: lat !< latitude
real(rk),intent(in) :: x !< x
real(rk),intent(in) :: y !< y
type(proj_lcc),intent(in) :: params !< projection parameters
real(rk) :: xx,yy,rho,theta
xx=x ; yy=y
! Account for false eastings and northings
xx = xx - params%false_easting
yy = yy - params%false_northing
rho = sign(sqrt(xx**2.0 + (params%rho0-yy)**2.0),params%n)
if (params%n > 0.0) then
theta = atan2(xx,(params%rho0-yy))
else
theta = atan2(-xx,(yy-params%rho0))
end if
if (abs(rho) < CONV_LIMIT) then
lat = sign(real(90.0,kind=rk),params%n)
else
lat = R2D * (2.0 * atan((EQ_RAD*params%f/rho)**params%i_n) - M_PI_2)
end if
lon = params%longitude_of_central_meridian+R2D*theta*params%i_n
end subroutine glimmap_ilcc
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Stereographic projection
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Forward transformation: lat-lon -> x-y of Stereographic projection
subroutine glimmap_stere(lon,lat,x,y,params)
use glimmer_log
real(rk),intent(in) :: lon !< longitude
real(rk),intent(in) :: lat !< latitude
real(rk),intent(out) :: x !< x
real(rk),intent(out) :: y !< y
type(proj_stere),intent(in) :: params !< projection parameters
real(rk) :: dlon,k,dlat,slat,clat,slon,clon
character(80) :: errtxt
dlon = lon-params%longitude_of_central_meridian
! Check domain of longitude
dlon = loncorrect(dlon,-180.0_rk)
dlon = dlon * D2R
dlat = lat * D2R
call sincos(dlon,slon,clon)
select case(params%pole)
case(1) ! North pole
x = 2.0 * params%k0 * tan(M_PI_4 - dlat/2.0)*slon
y = -2.0 * params%k0 * tan(M_PI_4 - dlat/2.0)*clon
case(-1) ! South pole
x = 2.0 * params%k0 * tan(M_PI_4 + dlat/2.0)*slon
y = 2.0 * params%k0 * tan(M_PI_4 + dlat/2.0)*clon
case(0) ! Oblique
call sincos(dlat,slat,clat)
if (params%equatorial) then
k = 2.0 * params%k0 / (1.0 + clat*clon)
y = k * slat
else
k = 2.0 * params%k0 / (1.0 + params%sinp*slat + params%cosp*clat*clon)
y = k * (params%cosp*slat - params%sinp*clat*clon)
end if
x = k * clat * slon
case default
write(errtxt,*)'Stereographic projection error:',params%pole
call write_log(trim(errtxt),GM_FATAL,__FILE__,__LINE__)
end select
! Apply false eastings and northings
x = x + params%false_easting
y = y + params%false_northing
end subroutine glimmap_stere
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Inverse transformation: lat-lon -> x-y of Stereographic projection
subroutine glimmap_istere(lon,lat,x,y,params)
real(rk),intent(out) :: lon !< longitude
real(rk),intent(out) :: lat !< latitude
real(rk),intent(in) :: x !< x
real(rk),intent(in) :: y !< y
type(proj_stere),intent(in) :: params !< projection parameters
real(rk) :: xx,yy,rho,c,sinc,cosc
xx=x ; yy=y
! Account for false eastings and northings
xx = xx - params%false_easting
yy = yy - params%false_northing
rho = hypot(xx,yy)
if (abs(rho) transform from grid to space
subroutine grid2space(x,y,gx,gy,coordsys)
use glimmer_coordinates
implicit none
real(rk),intent(out) :: x !< x-location in real space
real(rk),intent(out) :: y !< y-location in real space
real(rk),intent(in) :: gx !< x-location in grid space
real(rk),intent(in) :: gy !< y-location in grid space
type(coordsystem_type), intent(in) :: coordsys !< coordinate system
x=coordsys%origin%pt(1) + real(gx - 1)*coordsys%delta%pt(1)
y=coordsys%origin%pt(2) + real(gy - 1)*coordsys%delta%pt(2)
end subroutine grid2space
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> convert from space to grid
subroutine space2grid(x,y,gx,gy,coordsys)
use glimmer_coordinates
implicit none
real(rk),intent(in) :: x !< x-location in real space
real(rk),intent(in) :: y !< y-location in real space
real(rk),intent(out) :: gx !< x-location in grid space
real(rk),intent(out) :: gy !< y-location in grid space
type(coordsystem_type), intent(in) :: coordsys !< coordinate system
gx = 1.0 + (x - coordsys%origin%pt(1))/coordsys%delta%pt(1)
gy = 1.0 + (y - coordsys%origin%pt(2))/coordsys%delta%pt(2)
end subroutine space2grid
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Calculates the sin and cos of an angle.
subroutine sincos(a,s,c)
implicit none
real(rk),intent(in) :: a !< Input value (radians).
real(rk),intent(out) :: s !< sin(a)
real(rk),intent(out) :: c !< cos(a)
s = sin (a)
c = cos (a)
end subroutine sincos
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Normalises a value of longitude to the range starting at min degrees.
!! \return The normalised value of longitude.
real(rk) function loncorrect(lon,minimum)
real(rk),intent(in) :: lon !< The longitude under consideration (degrees east)
real(rk),intent(in) :: minimum !< The lower end of the output range (degrees east)
real(rk) :: maximum
loncorrect = lon
maximum = minimum + 360.0
do while (loncorrect >= maximum)
loncorrect=loncorrect-360.0
enddo
do while (loncorrect < minimum)
loncorrect=loncorrect+360.0
enddo
end function loncorrect
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> compute \f$\sqrt{x^2+y^2}\f$
real(rk) function hypot(x,y)
implicit none
real(rk),intent(in) :: x !< One input value
real(rk),intent(in) :: y !< Another input value
hypot=sqrt(x*x+y*y)
end function hypot
end module glimmer_map_trans