oi_latlon_conf_proj.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!      ###################################################
       SUBROUTINE oi_latlon_conf_proj(NDIM,PLAT0,PLON0,PRPK,PBETA,PLATOR,PLONOR, &
                                       px,py,plat,plon                       )  
!      ###################################################
!
!!****  *LATLON_CONF_PROJ * - Routine to compute geographical coordinates
!!
!!     PURPOSE
!!     -------
!        This routine computes the latitude and longitude of
!      an array given in cartesian conformal coordinates
!        Five map projections are available: 
!      - polar-stereographic from south pole  (PRPK=1),
!      - lambert conformal from south pole  (0<PRPK<1),
!      - mercator                             (PRPK=0),
!      - lambert conformal from north pole (-1<PRPK<0),
!      - polar-stereographic from north pole  (PRPK=-1).
!
!
!!**   METHOD
!!     ------
!!       Spherical earth approximation is used. Longitude origin is 
!!     set in Greenwich, and is positive eastwards. An anticlockwise 
!!     rotation of PBETA degrees is applied to the conformal frame 
!!     with respect to the geographical directions.
!!
!!       WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
!!
!!     EXTERNAL
!!     --------
!!       None
!!
!!     REFERENCE
!!     ---------
!!      Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
!!            commun CNRM-LA, specifications techniques", 
!!            Note CNRM/GMME, 26, 139p, (Chapter 2).
!!      Ducrocq V., 1994, "Generation de la grille dans le modele",
!!            Note interne MNH, 5 mai, 3p.
!!      Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
!!            Internal note ARPEGE/ALADIN, february 27,28p.
!!      Levallois J., 1970, "Geodesie generale", Tome 2, Collection
!!             de l'IGN, Eyrolles, Paris, 408p.
!!       
!!     AUTHOR
!!     ------
!!      P.M.       *LA*
!!
!!     MODIFICATION
!!     ------------
!!       Original  PM  24/05/94
!!       Updated   PM  27/07/94
!!       Updated   VD  23/08/94
!!       Updated   VM  24/10/95 projection from north pole (PRPK<0) and 
!!                              longitudes set between PLON0-180. and PLON0+180.
!!       Update    VM  11/2004  externalized version
!-------------------------------------------------------------------------------
!
!*     0.     DECLARATIONS
!             ------------
!
USE modd_csts,ONLY : xpi, xradius
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*     0.1    Declarations of arguments and results
!
REAL,                 INTENT(IN) :: plat0  ! Reference latitude
REAL,                 INTENT(IN) :: plon0  ! Reference longitude
REAL,                 INTENT(IN) :: prpk   ! projection parameter 
!                                          !   K=1 : stereographic north pole
!                                          ! 0<K<1 : Lambert, north hemisphere
!                                          !   K=0 : Mercator
!                                          !-1<K<0 : Lambert, south hemisphere
!                                          !   K=-1: stereographic south pole
REAL,                 INTENT(IN) :: pbeta  ! angle between grid and reference longitude
REAL,                 INTENT(IN) :: plator ! Latitude of the origine point
REAL,                 INTENT(IN) :: plonor ! Longitude of the origine point
REAL, DIMENSION(NDIM),   INTENT(IN) :: px,py
                                           ! given conformal coordinates of the 
                                           ! processed points (meters);
REAL, DIMENSION(NDIM),   INTENT(OUT):: plat,plon    
                                           ! returned geographic latitudes and 
                                           ! longitudes of the processed points 
                                           ! (degrees).
INTEGER,              INTENT(IN) :: ndim
!
!*     0.2    Declarations of local variables
! 
REAL, DIMENSION(NDIM)     :: zy
REAL                      :: zrpk,zbeta,zlat0,zlon0,zlator,zlonor
REAL                      :: zrdsdg,zclat0,zslat0,zclator,zslator
REAL                      :: zxbm0,zybm0,zro0,zga0 
REAL                      :: zxp,zyp,zepsi,zt1,zcgam,zsgam,zraclat0
!
REAL, DIMENSION(NDIM)     :: zata,zro2,zt2,zxmi0,zymi0
REAL(KIND=JPRB) :: zhook_handle
!
!--------------------------------------------------------------------------------
!
!*     1.     Preliminary calculations for all projections
!             --------------------------------------------
!  
IF (lhook) CALL dr_hook('OI_LATLON_CONF_PROJ',0,zhook_handle)
zrdsdg = xpi/180.             ! Degree to radian conversion factor
zepsi  = 10.*epsilon(1.)      ! A small number
!
! By definition, (PLONOR,PLATOR) are the geographical 
! coordinates, and (ZXBM0,ZYBM0) the conformal cartesian 
! point of coordinates (x=0,y=0).
!
zxbm0 = 0.
zybm0 = 0.
!
!-------------------------------------------------------------------------------
!
!*     2.     POLAR STEREOGRAPHIC AND LAMBERT CONFORMAL CASES
!             -----------------------------------------------
!                   (PRPK=1 P-stereo, 0<PRPK<1 Lambert)
!
IF(prpk /= 0.) THEN
!
  IF (prpk<0.) THEN     ! projection from north pole
    zrpk=-prpk
    zbeta=-pbeta
    zlat0=-plat0
    zlon0=plon0+180.
    zlator=-plator
    zlonor=plonor+180.
    zy(:)=-py(:)
    zybm0=-zybm0
  ELSE                  ! projection from south pole
    zrpk=prpk
    zbeta=pbeta
    zlat0=plat0
    zlon0=plon0
    zlator=plator
    zlonor=plonor
    zy(:)=py(:)
  ENDIF    
!
!*     2.1    Preliminary calculations
!
  zclat0  = cos(zrdsdg*zlat0)
  zslat0  = sin(zrdsdg*zlat0)
  zclator = cos(zrdsdg*zlator)
  zslator = sin(zrdsdg*zlator)
  zro0    = (xradius/zrpk)*(abs(zclat0))**(1.-zrpk)     &
            * ((1.+zslat0)*abs(zclator)/(1.+zslator))**zrpk  
  zga0    = (zrpk*(zlonor-zlon0)-zbeta)*zrdsdg
  zxp     = zxbm0-zro0*sin(zga0)
  zyp     = zybm0+zro0*cos(zga0)
!
!*    2.2    Longitude
!
  WHERE  (abs(zy(:)-zyp) < zepsi    &
       .AND.abs(px(:)-zxp) < zepsi)  
    zata(:) = 0.
  ELSEWHERE
    zata(:) = atan2(-(zxp-px(:)),(zyp-zy(:)))/zrdsdg
  END WHERE
  !
  plon(:) = (zbeta+zata(:))/zrpk+zlon0
!
!*   2.3     Latitude
!
  zro2(:) = (px(:)-zxp)**2+(zy(:)-zyp)**2
  zt1     = (xradius*(abs(zclat0))**(1.-zrpk))**(2./zrpk)   &
            * (1+zslat0)**2  
  zt2(:)  = (zrpk**2*zro2(:))**(1./zrpk)
  !
  plat(:) = (xpi/2.-acos((zt1-zt2(:))/(zt1+zt2(:))))/zrdsdg
!
  IF (prpk<0.) THEN     ! projection from north pole
    plat(:)=-plat(:)
    plon(:)=plon(:)+180.
  ENDIF
!
!-------------------------------------------------------------------------------
!
!*  3.        MERCATOR PROJECTION WITH ROTATION
!             ---------------------------------
!                       (PRPK=0)
!
ELSE
!
!*  3.1       Preliminary calculations
!
  zcgam    = cos(-zrdsdg*pbeta)
  zsgam    = sin(-zrdsdg*pbeta)
  zraclat0 = xradius*cos(zrdsdg*plat0)
!
!*  3.2       Longitude
!
  zxmi0(:) = px(:)-zxbm0
  zymi0(:) = py(:)-zybm0
!
  plon(:) = (zxmi0(:)*zcgam+zymi0(:)*zsgam)     &
            / (zraclat0*zrdsdg)+plonor  
!
!*  3.3       Latitude
!
  zt1     = alog(tan(xpi/4.+plator*zrdsdg/2.))
  zt2(:)  = (-zxmi0(:)*zsgam+zymi0(:)*zcgam)/zraclat0
  !
  plat(:) = (-xpi/2.+2.*atan(exp(zt1+zt2(:))))/zrdsdg
!
!---------------------------------------------------------------------------------
!
!*  4.        EXIT
!             ----
!
END IF
plon(:)=plon(:)+nint((plon0-plon(:))/360.)*360.
IF (lhook) CALL dr_hook('OI_LATLON_CONF_PROJ',1,zhook_handle)
!---------------------------------------------------------------------------------
END SUBROUTINE oi_latlon_conf_proj