eggrvs.F90

Go to the documentation of this file.

SUBROUTINE eggrvs (PRPI, PRA, PDELX, PDELY, KPROF,&
 & KBEG, KEND, KULOUT, PGELAM, PGELAT, PGM, PGNORX, PGNORY)  
!****
!---------------------------------------------------------------------

!     GEOGRAPHY OF GRID-POINTS, INVERSION FROM GRID TO GEOGRAPHICAL SPHERE
!     ARPEGE-ALADIN
!     --------------------------------------------------------------------

!       ---------------------------------------------------
!     PURPOSE
!     -------
!      KNOWING THE PRECISE GEOGRAPHICAL TRANSFORMATION FROM
!      ARGUMENTS AND COMMON /YEMGGCM/, COMPUTES THE LOCATION
!      ON THEE GEOGRAPHICAL POINTS GIVEN IN INPUT BY THEIR LOCATION
!      ON THE ARPEGE-ALADIN GRID

!      MUST BE USED IN CONNECTION WITH SUBROUTINE EGGX
!      EITHER WITHIN IT OR AFTER IT

!     INPUT PARAMETERS
!     ----------------
!      PRPI : PI (3.14ETC)
!      PRA  : A, RADIUS OF PLANET
!      PDELX, PDELY : GRID SIZE IN M IF PROJECTION, OR IN RADIANS
!      KPROF : SIZE OF (1D) ARRAYS
!      KBEG, KEND : BEGINNING AND END POINTS OF CALCULATIONS
!      KULOUT : LOGICAL UNIT OF LISTING
!      PGELAM(KPROF) : X LOCATION OF POINTS, DISTANCE UNDER PROJECTION,
!                            RELATIVE ROTATED LONGITUDE UNDER ROTATION,
!                      DEFINED AS (JLON-KDLUN)*PDELX
!      PGELAT(KPROF) : Y LOCATION OF POINTS, DISTANCE UNDER PROJECTION,
!                            RELATIVE ROTATED LATITUDE UNDER ROTATION,
!                      DEFINED AS (JLAT-KDGUN)*PDELY
!            UNDER ROTATION, THE POSITION OF THE ORIGIN (XLAT1R,XLON1U)
!            IS HANDLED BY THIS SUBROUTINE : ONLY RELATIVE LOCATION
!            NEED TO BE SPECIFIED

!     IMPLICIT INPUT
!     --------------
!      COMMON /YEMGGCM/ MUST HAVE BEEN INITIALIZED

!     OUTPUT PARAMETERS
!     -----------------
!      PGELAM (KPROF): GEOGRAPHICAL LONGITUDE
!      PGELAT (KPROF): GEOGRAPHICAL LATITUDE
!      PGM    (KRPOF): MAP FACTOR
!      PGNORX (KPROF): PROJECTION OF GEOGRAPHICAL NORTH ON X
!      PGNORY (KPROF): PROJECTION OF GEOGRAPHICAL NORTH ON Y

!     WRITTEN BY
!     ---------- ALAIN JOLY

!      ORIGINAL NORTHERN HEMISPHERE VERSION : 27/2/92
!      SOUTHERN HEMISPHERE VERSION : 27/1/93

!     Modified:
!     --------

!            98-05-07: P. Le Moigne :vectorization of eggrvs (LLSTOP)

!-------------------------------------------------------------------

USE parkind1  ,ONLY : jpim     ,jprb
USE yomhook   ,ONLY : lhook,   dr_hook

USE yemggcm  , ONLY : nymggi   ,nymggr   ,nymggwh  ,xlatr    ,&
 & xlonr    ,xggpk    ,xrpksm   ,xlat0r   ,xlon0r   ,&
 & xlon0u   ,xipore   ,xjpore   ,xggm0    ,xlon1u   ,&
 & xlat1r   ,hsud     ,xbeta  

!-------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KPROF 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRPI 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRA 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PDELX 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PDELY 
INTEGER(KIND=JPIM),INTENT(IN)    :: KBEG 
INTEGER(KIND=JPIM),INTENT(IN)    :: KEND 
INTEGER(KIND=JPIM),INTENT(IN)    :: KULOUT 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PGELAM(kprof) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PGELAT(kprof) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGM(kprof) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGNORX(kprof) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGNORY(kprof) 

!-------------------------------------------------------------------

INTEGER(KIND=JPIM) :: JJ

LOGICAL :: LLGWH, LLSTOP

REAL(KIND=JPRB) :: Z2PIPK, ZCOBETA, ZCOLA, ZCOSA, ZCOSOG, ZDIST,&
 & ZFUN, ZGAM, ZGM, ZKDL, ZLAT, ZLATG, ZLON, &
 & ZLONG, ZNORX, ZNORXP, ZNORY, ZNORYP, ZPIS2, &
 & ZPIS4, ZRPKSM2, ZSECAN, ZSECUR, ZSIBETA, &
 & ZSINAG, ZSINOG, ZURA2, ZX, ZY  
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------

#include "abor1.intfb.h"

!-------------------------------------------------------------------
IF (lhook) CALL dr_hook('EGGRVS',0,zhook_handle)
!-------------------------------------------------------------------

llstop=.false.
IF ( nymggi /= 10 ) THEN
  WRITE (kulout,*) '*** EGGRVS *** UNINITIALISED MODULE '
  CALL abor1(' EGGRVS: NYMGGI /= 10 ')
ENDIF

!     INITIALISE ROTATION ANGLE AND OTHER CONSTANTS
!     ---------------------------------------------
zpis2 = prpi*0.5_jprb
zpis4 = prpi*0.25_jprb
zsecur = 1.e-12_jprb
zsecan = 1.e-05_jprb
llgwh = .false.
IF ( nymggwh == 1 ) llgwh = .true.

DO jj = kbeg, kend
  pgm(jj) = 1.0_jprb
  pgnorx(jj) = 0.0_jprb
  pgnory(jj) = 1.0_jprb
ENDDO

!     CORRECTION OF POSITION UNDER ROTATION ONLY
!     -------------------------------------------

IF ( xggpk < 0.0_jprb ) THEN
  DO jj = kbeg, kend
    pgelam(jj) = xlon1u + pgelam(jj)
    IF ( pgelam(jj) >= 2.0_jprb*prpi )pgelam(jj) = pgelam(jj) - 2.0_jprb*prpi
    pgelat(jj) = xlat1r + pgelat(jj)
  ENDDO
ENDIF

!*
!--------------------------------------------------------------------
!     1.- REVERSE STEREO-LAMBERT PROJECTION
!     -------------------------------------
IF ( xggpk > 0.0_jprb ) THEN
  zrpksm2 = xrpksm*xrpksm
  zura2 = 1.0_jprb/( pra*pra )
  z2pipk = 2.0_jprb*prpi*xggpk

  DO jj = kbeg, kend
    zx = pgelam(jj) - xipore*pdelx
    zy = pgelat(jj) - xjpore*pdely
    zdist = (zrpksm2*( zx*zx + zy*zy )*zura2)**(1.0_jprb/(2.0_jprb*xggpk))

    zlat = zpis2 - 2.0_jprb*atan( zdist )

    IF ( zdist  < zsecur ) THEN
      ! THE POLE IS TOO NEAR TO DEFINE LONGITUDE
      zlon = 0.0_jprb
      IF ( xggpk /= 1.0_jprb ) llstop=.true.
      ! MAP FACTOR AT POLE FOR STEREOGRAPHIC PROJECTION
      zgm = ( 1.0_jprb + sin( xlat0r ) )/( 1.0_jprb + sin( zlat ) )
    ELSE
      IF ( hsud > 0.0_jprb ) THEN
        zkdl = atan2( zx,-zy )
      ELSE
        zkdl = atan2( zx,zy )
      ENDIF
      zlon = xlon0r + (zkdl + xbeta)/xggpk
      IF ( zlon < 0.0_jprb ) zlon = 2.0_jprb*prpi + zlon
      zgm = xggm0*( cos( zlat )**(xggpk-1.0_jprb) )*&
       & ( ( 1.0_jprb + sin( zlat ) )**(-xggpk) )  
    ENDIF

    pgelam(jj) = zlon
    pgelat(jj) = hsud*zlat
    pgm(jj) = zgm
    ! COMPONENTS OF VECTOR ROTATION MATRIX
    IF ( llgwh ) zlon = zlon + 2.0_jprb*prpi
    IF ( xggpk < 1.0_jprb .AND. (zlon-xlon0u)  >  z2pipk )&
     & zlon = zlon - 2.0_jprb*prpi  
    zgam = hsud*(xggpk*( zlon - xlon0u ) - xbeta)
    pgnorx(jj) = -sin( zgam )
    pgnory(jj) = cos( zgam )
  ENDDO
  IF ( llstop ) THEN
    WRITE (kulout,*) ' POLE WITHIN LAMBERT DOMAIN '
    CALL abor1(' EGGRVS: POLE WITHIN LAMBERT DOMAIN ')
  ENDIF
ENDIF

!*
!--------------------------------------------------------------------
!     2.- REVERSE MERCATOR PROJECTION
!     -------------------------------
IF ( xggpk == 0.0_jprb ) THEN
  zsibeta = sin( xbeta )
  zcobeta = cos( xbeta )
  DO jj = kbeg, kend
    zy = pgelam(jj)*zsibeta + pgelat(jj)*zcobeta - xjpore*pdely
    zdist = exp( -zy/( pra*cos(xlat0r) ) )
    zlat = zpis2 - 2.0_jprb*atan( zdist )

    zgm = cos( xlat0r )/cos( zlat )

    zx = pgelam(jj)*zcobeta - pgelat(jj)*zsibeta - xipore*pdelx
    zlon = xlon0u + zx/( pra*cos( xlat0r ) )
    IF ( zlon >= 2.0_jprb*prpi ) zlon = zlon - 2.0_jprb*prpi

    pgelam(jj) = zlon
    pgelat(jj) = zlat
    pgm(jj) = zgm
    ! COMPONENTS OF VECTOR ROTATION MATRIX
    pgnorx(jj) = zsibeta
    pgnory(jj) = zcobeta
  ENDDO
ENDIF

!*
!--------------------------------------------------------------------
!     3.- REVERSE ROTATION
!     --------------------
IF ( nymggr /= 0 ) THEN
  DO jj = kbeg, kend
    zlon = pgelam(jj)
    zlat = pgelat(jj)
    zsinag = sin( xlatr )*cos( zlat )*cos( zlon ) +cos( xlatr )*sin( zlat )
    zsinag = min(1.0_jprb,max(-1.0_jprb,zsinag))
    zlatg = asin( zsinag )
    IF ( abs( zlatg ) >= zpis2 ) THEN
      zlong = 0.0_jprb
    ELSE
      zcosa = cos( zlatg )
      zfun = cos( xlatr )*cos( zlat )*cos( zlon ) -sin( xlatr )*sin( zlat )
      zcosog = ( cos( xlonr )*zfun - sin( xlonr )*cos( zlat )*&
       & sin( zlon ) )/zcosa  
      zcosog = min(1.0_jprb,max(-1.0_jprb,zcosog))
      zsinog = ( sin( xlonr )*zfun + cos( xlonr )*cos( zlat )*&
       & sin( zlon ) )/zcosa  
      zsinog = min(1.0_jprb,max(-1.0_jprb,zsinog))
      zlong = acos( zcosog )
      IF ( asin(zsinog) < 0.0_jprb ) zlong = 2.0_jprb*prpi - zlong
    ENDIF
    pgelam(jj) = zlong
    pgelat(jj) = zlatg
    ! COMPONENTS OF ROTATION MATRIX DUE TO ROTATION
    zcola = sqrt( abs( 1.0_jprb - (cos( xlatr )*sin( zlatg )-&
     & sin( xlatr )*cos( zlatg )*cos( zlong-xlonr ))**2 ) )  
    IF ( zcola < zsecur ) THEN
      WRITE (kulout,*) ' *** EGGX QUASI ERROR ***',&
       & ' DOMAIN EXTENDS UP TO NEW POLE : IT IS PROBABLY TOO LARGE'  
      znory = 1.0_jprb
      znorx = 0.0_jprb
    ELSE
      znory = ( cos( xlatr )*cos( zlatg ) + sin( xlatr )*&
       & sin( zlatg )*cos( zlong-xlonr ) )/zcola  
      znorx = - sin( xlatr )*sin( zlong-xlonr )/zcola
    ENDIF
    ! COMPOSITION OF THIS ROTATION WITH THE ONE RESULTING FROM PROJECTION
    znoryp = pgnory(jj)
    znorxp = pgnorx(jj)
    pgnory(jj) = znory*znoryp - znorx*znorxp
    pgnorx(jj) = znorx*znoryp + znory*znorxp
  ENDDO
ENDIF

!-------------------------------------------------------------------
IF (lhook) CALL dr_hook('EGGRVS',1,zhook_handle)
END SUBROUTINE eggrvs