eggdir.F90

Go to the documentation of this file.

SUBROUTINE eggdir (PRPI, PRA, PDELX, PDELY, KPROF,&
 & KBEG, KEND, KULOUT, PGX, PGY)  
!****
!---------------------------------------------------------------------

!     GEOGRAPHY OF GRID-POINTS, ROTATION/PROJECTION FROM GEOGRAPHICAL
!     SPHERE TO ARPEGE-ALADIN
!     --------------------------------------------------------------------

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

!      MUST BE USED IN CONNECTION WITH SUBROUTINE EGGX
!      AFTER IT HAS INITIALIZED /YEMGGCM/

!     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
!      PGX(KPROF) : INPUT: TRUE GEOGRAPHICAL LONGITUDE IN RADIANS
!                          (SEE EGGX FOR CONVENTIONS)
!      PGY(KPROF) : INPUT : TRUE GEOGRAPHICAL LATITUDE IN RADIANS

!            IN SHORT, THE OUTPUT IS SUITABLE FOR USE IN EGGRVS

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

!     OUTPUT PARAMETERS
!     -----------------
!      PGX (KPROF):
!                   OUTPUT: X LOCATION OF POINTS, DISTANCE UNDER PROJECTION,
!                            RELATIVE ROTATED LONGITUDE UNDER ROTATION,
!                            ALWAYS DEFINED AS (JLON-KDLUN)*PDELX
!      PGY (KPROF):
!                   OUTPUT : Y LOCATION OF POINTS, DISTANCE UNDER PROJECTION,
!                            RELATIVE ROTATED LATITUDE UNDER ROTATION,
!                            ALWAYS DEFINED AS (JLAT-KDGUN)*PDELY
!            UNDER ROTATION, THE POSITION OF THE ORIGIN (XLAT1R,XLON1U)
!            IS HANDLED BY THIS SUBROUTINE : ONLY RELATIVE LOCATION
!            ARE GIVEN ON OUTPUT.

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

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

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

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

USE yemggcm  , ONLY : nymggi   ,nymggr   ,nymggwh  ,xlatr    ,&
 & xlonr    ,xggpk    ,xlat0r   ,xlon0u   ,xipore   ,&
 & xjpore   ,xlon1r   ,xlon1u   ,xlat1r   ,xlon2r   ,&
 & 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) :: PGX(kprof) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PGY(kprof) 

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

INTEGER(KIND=JPIM) :: JJ

LOGICAL :: LLGWH

REAL(KIND=JPRB) :: Z2PI, Z2PIPK, ZCOBETA, ZCOS, ZCOSO, ZGAM,&
 & ZLAT, ZLIMIT, ZLON, ZPIS2, ZPIS4, ZRC, ZRR, &
 & ZSECAN, ZSECUR, ZSIBETA, ZSIN, ZSINO, ZXE, &
 & ZXP, ZYE, ZYP  
REAL(KIND=JPRB) :: ZHOOK_HANDLE

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

#include "abor1.intfb.h"

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

zrc=0.0_jprb

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

IF ( hsud < 0.0_jprb ) THEN
  DO jj = kbeg, kend
    pgy(jj) = abs( pgy(jj) )
  ENDDO
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.
!     CHECK LONGITUDES
DO jj = kbeg, kend
  IF ( pgx(jj) < 0.0_jprb ) THEN
    pgx(jj) = pgx(jj) + 2.0_jprb*prpi
  ENDIF
ENDDO

!*
!----------------------------------------------------------------------
!     1.- DIRECT ROTATION
!     --------------------
IF ( nymggr /= 0 ) THEN
  DO jj = kbeg, kend
    zsin = cos( xlatr )*sin( pgy(jj) ) -&
     & sin( xlatr )*cos( pgy(jj) )*cos( pgx(jj)-xlonr )  
    zlat = asin( zsin )
    IF ( abs(zlat) >= zpis2 ) THEN
      zlon = 0.0_jprb
    ELSE
      zcos = cos( zlat )
      zcoso = ( sin( xlatr )*sin( pgy(jj) ) +&
       & cos( xlatr )*cos( pgy(jj) )*cos( pgx(jj)-xlonr ) )/zcos  
      zcoso = min(1.0_jprb,max(-1.0_jprb,zcoso))
      zsino = ( cos( pgy(jj) )*sin( pgx(jj)-xlonr ) )/zcos
      zsino = min(1.0_jprb,max(-1.0_jprb,zsino))
      zlon = acos( zcoso )
      IF ( asin( zsino ) < 0.0_jprb ) zlon = 2.0_jprb*prpi - zlon
    ENDIF
    pgx(jj) = zlon
    pgy(jj) = zlat
  ENDDO
ENDIF

!     CORRECTION OF POSITION IF DOMAIN ASTRIDE GREENWICH
!     --------------------------------------------------
IF ( llgwh ) THEN
  DO jj = kbeg, kend
    zlimit=0.5_jprb*(xlon2r+xlon1r-zpis2)
    IF ( pgx(jj) < xlon1r .AND. pgx(jj) < zlimit )&
     & pgx(jj) = pgx(jj) + 2.0_jprb*prpi  
  ENDDO
ENDIF

!     FINAL RESULT IN THE ABSENCE OF PROJECTION
!     -----------------------------------------
IF ( xggpk < 0.0_jprb ) THEN
  DO jj = kbeg, kend
    pgx(jj) = pgx(jj) - xlon1u
    pgy(jj) = pgy(jj) - xlat1r
  ENDDO
ENDIF

!*
!------------------------------------------------------------------
!     2.- DIRECT STEREO-LAMBERT PROJECTION
!     ------------------------------------
IF ( xggpk > 0.0_jprb ) THEN
  IF ( xggpk < 1.0_jprb ) THEN
    zrc = pra*( cos( xlat0r )**(1.0_jprb-xggpk) )*&
     & ( (1.0_jprb+sin( xlat0r ))**xggpk )/xggpk  
    z2pi = 2.0_jprb*prpi
    z2pipk = 2.0_jprb*prpi*xggpk
    DO jj = kbeg, kend
      IF ( (pgx(jj)-xlon0u)  >  z2pipk ) THEN
        pgx(jj) = pgx(jj) - z2pi
      ENDIF
    ENDDO
  ELSEIF ( xggpk == 1.0_jprb ) THEN
    zrc = pra*( 1.0_jprb + sin( xlat0r ) )
  ELSE
    WRITE (kulout,*) ' *** EGGDIR *** UNKNOWN PROJECTION '
  ENDIF
  zxp = xipore*pdelx
  zyp = xjpore*pdely

  DO jj = kbeg, kend
    zrr = zrc*( cos(pgy(jj)) /&
     & ( 1.0_jprb+ sqrt(1.0_jprb-cos(pgy(jj))**2 ))  )**xggpk      
    zgam = xggpk*( pgx(jj)-xlon0u ) - xbeta
    pgx(jj) = zxp + zrr*sin( zgam )
    pgy(jj) = zyp - hsud*zrr*cos( zgam )
  ENDDO
ENDIF

!*
!----------------------------------------------------------------------
!     3.- DIRECT MERCATOR PROJECTION
!     ------------------------------
IF ( xggpk == 0.0_jprb ) THEN
  zrc = pra*cos( xlat0r )
  zxe = xipore*pdelx
  zye = xjpore*pdely
  zsibeta = sin( xbeta )
  zcobeta = cos( xbeta )
  DO jj = kbeg, kend
    zxp = zxe + zrc*( pgx(jj)-xlon0u )
    zyp = zye - zrc*log( tan(zpis4 - 0.5_jprb*pgy(jj)) )
    pgx(jj) = zxp*zcobeta + zyp*zsibeta
    pgy(jj) = -zxp*zsibeta + zyp*zcobeta
  ENDDO
ENDIF

!-------------------------------------------------------------------
IF (lhook) CALL dr_hook('EGGDIR',1,zhook_handle)
END SUBROUTINE eggdir