geo_coord.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  g e o _ c o o r d
!
!> @file
!!
!! Computation of longitude lambda and latitude phi for position (x,y) in the
!! numerical domain.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Ralf Greve, Reinhard Calov, Alex Robinson
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS 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 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS 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 SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Computation of longitude lambda and latitude phi for position (x,y) in the
!! numerical domain.
!<------------------------------------------------------------------------------
subroutine geo_coord(phi_val, lambda_val, x_val, y_val)

use sico_types
use sico_variables

implicit none

real(dp), intent(in) :: x_val, y_val

real(dp), intent(out) :: lambda_val, phi_val

real(dp) :: k

#if (GRID==0 || GRID==1)

!-------- Inverse stereographic projection --------

call stereo_inv_ellipsoid(x_val, y_val, a, b, &
                          lambda0, phi0, lambda_val, phi_val)

#elif GRID==2

!-------- Identity --------

lambda_val = x_val
phi_val    = y_val

#endif

end subroutine geo_coord

!-------------------------------------------------------------------------------
!> Inverse stereographic projection for an ellipsoidal planet.
!<------------------------------------------------------------------------------
subroutine stereo_inv_ellipsoid(x_val, y_val, A, B, &
                                lambda0, phi0, lambda_val, phi_val)

use sico_types

implicit none
real(dp), intent(in)  :: x_val, y_val, a, b, lambda0, phi0
real(dp), intent(out) :: lambda_val, phi_val

integer(i4b) :: l
integer(i4b) :: sign_phi0
real(dp) :: phi_aux, phi0_aux
real(dp) :: e, mc, t, tc, kp, rho, phi_p, residual
real(dp) :: sinphi0, sinlambda0, cosphi0, coslambda0

real(dp), parameter :: pi = 3.141592653589793_dp, &
                       eps = 1.0e-05_dp, &
                       eps_residual = 1.0e-09_dp

if (phi0 > eps) then   ! for northern hemisphere
   sign_phi0 =  1
else if (phi0 < (-eps)) then   ! for southern hemisphere
   sign_phi0 = -1
else
   stop ' geo_coord: phi0 must be different from zero!'
end if

phi0_aux = phi0 * sign_phi0

e=sqrt((a**2-b**2)/(a**2))

sinphi0    = sin(phi0_aux)
sinlambda0 = sin(lambda0)
cosphi0    = cos(phi0_aux)
coslambda0 = cos(lambda0)
  
tc=sqrt(((1.0_dp-sinphi0)/(1.0_dp+sinphi0))* &
       ((1.0_dp+e*sinphi0)/(1.0_dp-e*sinphi0))**e)
mc=cosphi0/sqrt(1.0_dp-e*e*sinphi0*sinphi0)
rho=sqrt(x_val*x_val+y_val*y_val)
t=rho*tc/(a*mc)

lambda_val = lambda0 + sign_phi0*atan2(y_val,x_val) + 0.5_dp*pi

!  fix point iteration

phi_p=0.5_dp*pi-2.0_dp*atan(t)
l=0
residual=3600.0_dp
do while(residual >= eps_residual)
   l=l+1
   phi_aux=0.5_dp*pi-2.0_dp*atan(t*((1.0_dp-e*sin(phi_p))/ &
           (1.0_dp+e*sin(phi_p)))**(0.5_dp*e))
   residual=abs(phi_aux-phi_p)
   phi_p=phi_aux
end do

phi_val = phi_aux * sign_phi0

if (lambda_val < 0.0_dp) then
   lambda_val = lambda_val + 2.0_dp*pi
else if (lambda_val >= (2.0_dp*pi)) then
   lambda_val = lambda_val - 2.0_dp*pi
end if

end subroutine stereo_inv_ellipsoid
!