eggx.F90

Go to the documentation of this file.

SUBROUTINE eggx (PRPI, PRA, KROTEQ, PLONR, PLATR, PBETA,&
 & PLON1, PLAT1, PLON2, PLAT2, PLON0, PLAT0, PRPK, KULOUT,&
 & KSOTRP, KGIV0,&
 & PGELAM, PGELAT, PGM, PGNORX, PGNORY,&
 & KDLSA, KDLSUR, KDGSA, KDGEN, KDLUN, KDLUX, KDGUN, KDGUX,&
 & PDELX, PDELY)  
!****
!----------------------------------------------------------------------

!     GEOGRAPHY OF GRID-POINTS
!     ARPEGE-ALADIN
!     -------------------------

!       -------------------------------------------------------
!     PURPOSE
!     -------
!      PROVIDES THE BASIC GEOGRAPHICAL PARAMETERS OF EACH GRID-POINT
!      IN THE WINDOW OF INTERNAL+COUPLING ZONE, DEFINING THE
!      GEOGRAPHICAL DOMAIN OF INTEREST.
!      SIMULTANEOUSLY, INITIALISES THE INTERNAL GEOGRAPHICAL COMMON
!      /YEMGGCM/ FOR POSSIBLE RE-USE OF POSITION INVERSION SUBROUTINE
!      EGGRVS.

!      THE GEOGRAPHICAL PARAMETERS ARE :
!       - GEOGRAPHIC LONGITUDE
!       - GEOGRAPHIC LATITUDE
!       - MAP FACTOR
!       - COMPONENTS OF VECTOR DIRECTED TOWARDS THE GEOGRAPHIC NORTH
!       POLE FOR PROJECTION OF VECTORS
!       - X GRID-SIZE (DISTANCE OR LONGITUDE INCREMENT)
!       - Y GRID-SIZE (DISTANCE OR LATITUDE INCREMENT)
!      THE LAST TWO DEPEND ON THE GEOMETRY. THE LATTER CAN BE
!       - SPHERICAL GEOMETRY, WITH ROTATION OF DOMAIN TO THE EQUATOR
!       - CARTESIAN GEOMETRY, WITH PROJECTION OF DOMAIN ON A PLANE

!      ROTATION AND PROJECTION CAN BE COMBINED. HOWEVER, PROJECTIONS
!      OTHER THAN MERCATOR SHOULD BE USED WITHOUT ROTATION.

!      EGGX CAN RUN UNDER FULL USER CONTROL OR IN MORE OR LESS
!      AUTOMATIC MODES

!      FROM THESE BASIC PARAMETERS, OTHER CAN EASILY BE DEDUCED, SUCH
!      AS THE CORIOLIS PARAMETER, ETC.

!      (GENERAL WARNING ABOUT AUTOMATIC MODE)
!      --------------------------------------
!       NOTE THAT AUTOMATIC MODE IS NOT SURE TO CONVERGE, AND MOREOVER,
!      IT MAY CONVERGE TOWARD A SOLUTION DIFFERENT FROM WHAT YOU EXPECTED
!      THIS IS SO IN PARTICULAR FOR A DOMAIN OVER THE POLE OR NEAR THE
!      EQUATOR. IT IS RECOMMENDED TO CHECK CAREFULLY A CHOICE OF
!      PARAMETERS BEFORE RUNNING A COMPLETE SUITE

!      (CONVENTIONS FOR LATITUDE, LONGITUDE)
!      -------------------------------------
!       THEY ARE IN RADIANS
!       LATITUDE VARY FROM  -PI/2  (SOUTH POLE) TO  PI/2  (NORTH POLE)
!       LONGITUDE VARY FROM  0 (GREENWICH OR EQUIVALENT) TO 2*PI
!         THE CONVENTION FOR A POLE IS THAT IT HAS LONGITUDE 0

!      (HANDLING OF PROJECTION IN SOUTHERN HEMISPHERE)
!      ------------------------------------------------
!       IT IS ASSUMED THAT DOMAINS NOT ROTATED AND ENTERELY CONTAINED
!       IN THE SOUTHERN HEMISPHERE WILL BE PROJECTED USING STEREO/LAMBERT

!       IN THAT CASE, AND IN THAT CASE ONLY,
!         HSUD IS SET TO -1. IN THE INTERNAL COMMON

!       THE SIGN OF LATITUDES IS CHANGED FOR ALL INTERNAL CALCULATIONS,
!       SO THAT ALL THE TRIGONOMETRY IS UNCHANGED, HOWEVER, THE RESPECTIVE
!       POSITIONS OF POLE AND GRID ARE REVERSED

!       WHEREVER THIS DIFFERENCE MATTERS, THAT IS:
!       - MAPPING POLAR COORDINATES TO (X,Y)
!       - VECTOR ROTATION
!       ALTERNATIVE LINES OR MULTIPLICATIONS BY HSUD ARE INCLUDED

!       THE SIGNS OF OUTGOING OR INCOMING LATITUDES ARE CHANGED AT THE LAST
!        MINUTE OR PRIOR TO ANY COMPUTATION

!      (CONVENTIONS FOR ARRAY ORGANISATION)
!      ------------------------------------
!       THE POINT (KDLUN,KDGUN) HAS GEOGRAPHIC COORDINATES (PLON1,PLAT1)
!       THE POINT (KDLUX,KDGUX)                            (PLON2,PLAT2)
!        THE FIRST INDEX (X DIRECTION) VARIES ROUGHLY
!        WITH INCREASING LONGITUDE
!        THE SECOND INDEX (Y DIRECTION) VARIES
!        ROUGHLY WITH INCREASING LATITUDE
!        THEY ARE ACTUAL LONGITUDE, LATITUDE ONLY IN THE CASE WITHOUT PROJ.

!     ONLY A PART OF THE ARRAYS IS INITIALIZED : THAT PART CORRESPONDING
!     TO INTERNAL + COUPLING ZONES. IT MIGHT BE NECESSARY TO COMPLETE
!     THE ARRAYS (SUCH AS GM) IN ORDER TO MAKE THEM DOUBLY-PERIODIC.

!     ONE TYPE OF ROTATION AND THREE TYPES OF PROJECTIONS ARE RECOGNISED

!      (CONVENTIONS FOR WIND AND OTHER VECTORS ROTATION)
!      --------------------------------------------------
!       GEOGRAPHIC U = PGNORY * U PROJ - PGNORX * V PROJ
!       GEOGRAPHIC V = PGNORX * U PROJ + PGNORY * V PROJ

!       WHEN USING THIS TRANSFORM, DO NOT FORGET THE MAP FACTOR EFFECT

!       THE INVERSE TRANSFORM IS EASY TO DEDUCE (IT IS A ROTATION)

!       (OVERVIEW OF ALGORITHM)
!       -----------------------
!     PARAMETERS ARE EDITED  AND POSSIBLY CHECKED
!     THEY ARE OPTIONNALLY ROTATED
!     SOME AUTOMATIC ADJUSTMENTS OF PROJECTION PARAMETER ARE OPTIONALLY
!     PERFORMED
!     PROJECTION IS PERFORMED (OPTIONNALLY) ON THE ROTATED SPHERE
!     AT THIS STAGE, INCREMENTS OR DISTANCE ON THE MAP ARE DEDUCED
!     COMPONENTS OF ROTATION UNDER PROJECTION ARE COMPUTED
!     INVERSE ROTATION BACK TO GEOGRAPHICAL VALUES IS OPTIONNALLY DONE
!     ON THE RESULTS, COMPONENTS OF ROTATION ARE COMPOSED WITH THOSE
!     OF PROJECTION

!     INPUT PARAMETERS
!     ----------------
!      PRPI : PI (3.14ETC AS GENERALLY DEFINED ELSEWHERE : YOMCST IN ARPEGE)
!      PRA  : A, RADIUS OF SPHERICAL PLANET (M)
!       WARNING :
!       ---------  DEPENDING ON OPTIONS, INPUT PARAMETERS MAY BECOME OUTPUT
!                  PARAMETERS. EVEN IF UNITIALISED, THE VARIABLES MUST BE
!                  DECLARE INDEPENDENTLY.

!      KROTEQ = 0 : NO ROTATION
!             = 1 : POINT (PLONR,PLATR) IS ROTATED TO EQUATOR, THE NORTH POLE
!                   IS ON THE NEW GREENWICH MERIDIAN

!      PLONR : GEOGRAPHIC LONGITUDE OF REFERENCE POINT OF ROTATION
!      PLATR : GEOGRAPHIC LATITUDE OF REFERENCE POINT OF ROTATION

!      PBETA : ANGLE (IN RD) BETWEEN X-AXIS AND ROTATED LATITUDE CIRCLES
!              AT THE REFERENCE LONGITUDE
!              (USUALLY, PBETA = 0. : GIVES PURE PROJECTIONS)

!      PLON1, PLAT1 : GEOGRAPHIC LONGITUDE, LATITUDE OF THE SOUTH-WEST
!                     CORNER OF USEFUL DOMAIN
!      PLON2, PLAT2 : GEOGRAPHIC LONGITUDE, LATITUDE OF THE NORTH-EAST
!                     CORNER OF USEFUL DOMAIN
!       IF ROTATION IS REQUIRED, (PLON1,PLAT1) IS SELF DETERMINED
!       FROM (PLONR,PLATR) AND (PLON2,PLAT2) WHEN KSOTRP = 2
!       IF PROJECTION IS REQUIRED, PLAT1 IS SEFL DETERMINED WHEN
!       KSOTRP = 2. PLON1 MUST BE GIVEN.

!       IF ROTATION IS REQUIRED, (PLON2,PLAT2) IS SELF DETERMINED
!       FROM (PLONR,PLATR) AND (PLON1,PLAT1) WHEN KSOTRP = 1
!       IF PROJECTION IS REQUIRED, PLAT2 IS SEFL DETERMINED WHEN
!       KSOTRP = 1. PLON2 MUST BE GIVEN.

!      PLON0 : GEOGRAPHIC LONGITUDE OF REFERENCE FOR THE PROJECTION
!              (THE VERTICAL ONE IN STEREO/LAMBERT)
!        WARNING :
!        ---------  STRANGE RESULTS MAIN OCCUR IF PLON0 IS OUT OF THE
!                   DOMAIN
!        PLON0 IS NOT REQUIRED WHEN KGIV0 =1 OR 3
!      PLAT0 : GEOGRAPHIC LATITUDE OF REFERENCE FOR THE PROJECTION
!              (WHERE M = 1)
!        PLAT0 IS NOT REQUIRED WHEN KGIV0 = 2 OR 3
!      PRPK  : PROJECTION PARAMETER AND DEFINITION
!              PRPK = 10. PROJECTION TYPE SELF DETERMINED
!                         BY MINIMIZING THE VARIATION OF THE MAP FACTOR
!              PRPK = 1.  POLAR STEREOGRAPHIC PROJECTION
!              0. < PRPK < 1.  LAMBERT CONFORMAL PROJECTION WITH
!                              CONE PARAMETER PRPK
!              PRPK = 0.  MERCATOR CONFORMAL PROJECTION
!              PRPK < 0.  NO PROJECTION
!             ON OUTPUT, PRPK CONTAINS THE EFFECTIVE PROJECTION
!             PARAMETER THAT HAS BEEN USED
!      KSOTRP : ISOTROPY PARAMETER UNDER PROJECTION
!               = 0, PLAT1, PLAT2 USED AS GIVEN. GRID SPACING ISOTROPY
!                    GENERALLY NOT AVAILABLE ( PDELX DIFFERENT FORM PDELY )
!               = 1, GRID SPACING ISOTROPIC, PLAT1 USED AS GIVEN,
!                    PLAT2 CHANGED (PLON2 SELF DETERMINED IF KROTEQ = 1)
!               = 2, GRID SPACING ISOTROPIC, PLAT2 USED AS GIVEN,
!                    PLAT1 CHANGED (PLON1 SELF DETERMINED IF KROTEQ = 1)

!      KGIV0 : CHOICE OF REFERENCE POINT FOR PROJECTION
!               = 0, PLAT0 AND PLON0 USED AS GIVEN
!               = 1, PLAT0 REQUIRED, PLON0 SELF DETERMINED
!               = 2, PLON0 REQUIRED, PLAT0 SELF DETERMINED
!               = 3, PLAT0 AND PLON0 SELF DETERMINED

!      KDLSA:KDLSUR : LOWER AND UPPER FIRST DIMENSIONS OF ARRAYS (X)
!      KDGSA:KDGEN  : LOWER AND UPPER SECOND DIMENSIONS OF ARRAYS (Y)

!      KDLUN:KDLUX  : LOWER AND UPPER FIRST DIMENSIONS OF
!                     THE DOMAIN OF INTEREST, WHERE ARRAYS ARE
!                     INITIALIZED.
!      KDGUN:KDGUX  : LOWER AND UPPER SECOND DIMENSIONS OF
!                     THE DOMAIN OF INTEREST, WHERE ARRAYS ARE
!                     INITIALIZED.
!                   TOGETHER WITH THE CORNERS OF THE DOMAIN, THESE
!                   DEFINE THE GRID RESOLUTION

!      KULOUT : UNIT OF OUTPUT FILE

!     IMPLICIT INPUT
!     --------------
!       NONE
!       THIS IS AN ENVIRONMENT INDEPENDENT SUBROUTINE

!     OTHER INPUT  OR EXTERNALS
!     -------------------------
!       ALL MANNER OF TRIGONOMETRIC FUNCTIONS AND LOGARITHMS
!       CALLS EGGMLT TO COMPUTE MISSING LATITUDE WHEN KSOTRP > 0
!       CALLS EGGRVS TO DEDUCE GEOGRAPHIC PARAMETERS FROM FINAL
!       PROJECTION AND ROTATION PARAMETERS.

!       AFTER A FIRST CALL TO EGGX, ROUTINE EGGRVS CAN BE RE-CALLED
!       EXTERNALLY WITH DIFFERENT (PROJECTED OR) ROTATED POSITIONS

!     OUTPUT PARAMETERS
!     -----------------
!       PGELAM(KDLSA:KDLSUR,KDGSA:KDGEN) :
!                          GEOGRAPHIC LONGITUDE
!       PGELAT(KDLSA:KDLSUR,KDGSA:KDGEN) :
!                          GEOGRAPHIC LATITUDE
!       PGM(KDLSA:KDLSUR,KDGSA:KDGSUR) : MAP FACTOR
!       PGNORX(KDLSA:KDLSUR,KDGSA:KDGEN) :
!              PROJECTION OF GEOGRAPHICAL NORTH UNIT VECTOR ON X AXIS
!       PGNORY(KDLSA:KDLSUR,KDGSA:KDGEN) :
!              PROJECTION OF GEOGRAPHICAL NORTH UNIT VECTOR ON Y AXIS

!       PDELX : GRID SIZE IN M ALONG X IF PROJECTION
!               LONGITUDE INCREMENT IN RD IF SPHERICAL GEOMETRY
!               NECESSARY TO COMPUTE PERIOD AND DERIVATIVES
!       PDELY : GRID SIZE IN M ALONG Y IF PROJECTION
!               LATITUDE INCREMENT IN RD IF SPEHRICAL GEOMETRY

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

!     NEW NORTHERN HEMISPHERE VERSION : 27/2/92
!     SOUTHERN HEM ISPHERE VERSION : 27/1/93

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

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    ,xlon1r   ,&
 & xlon1u   ,xlat1r   ,xlon2r   ,hsud     ,xbeta  

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

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KDLSA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KDLSUR 
INTEGER(KIND=JPIM),INTENT(IN)    :: KDGSA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KDGEN 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRPI
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRA
INTEGER(KIND=JPIM),INTENT(IN)    :: KROTEQ 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLONR 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLATR 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PBETA 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLON1 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLAT1 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLON2 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLAT2 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLON0 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLAT0 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRPK
INTEGER(KIND=JPIM),INTENT(IN)    :: KULOUT
INTEGER(KIND=JPIM),INTENT(INOUT) :: KSOTRP
INTEGER(KIND=JPIM),INTENT(IN)    :: KGIV0 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGELAM(kdlsa:kdlsur,kdgsa:kdgen)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGELAT(kdlsa:kdlsur,kdgsa:kdgen)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGM(kdlsa:kdlsur,kdgsa:kdgen)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGNORX(kdlsa:kdlsur,kdgsa:kdgen)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PGNORY(kdlsa:kdlsur,kdgsa:kdgen)
INTEGER(KIND=JPIM),INTENT(IN)    :: KDLUN
INTEGER(KIND=JPIM),INTENT(IN)    :: KDLUX
INTEGER(KIND=JPIM),INTENT(IN)    :: KDGUN
INTEGER(KIND=JPIM),INTENT(IN)    :: KDGUX
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PDELX
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PDELY

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

INTEGER(KIND=JPIM) :: INBESS, INNEGA, ISPECA, ITERK, ITERKX, JLAT, JLON
INTEGER(KIND=JPIM) :: ISOTRP

LOGICAL :: LLGWH, LL510, LL520
LOGICAL :: LLPLANEX
LOGICAL :: LLPLANEY

REAL(KIND=JPRB) :: ZCONDEG, ZCONRAD, ZCOSA, ZCOSO, ZDCLA0, ZDCLA1,&
 & ZDCLA2, ZDLON, ZDM, ZDMMAX, ZDMMIN, ZDRK, &
 & ZDTLAT, ZFACE, ZIPV, ZJPV, ZLAT, ZLAT2R, &
 & ZLATLIM, ZLON, ZLON2U, ZPIS2, ZPIS4, ZRKAX, &
 & ZRKI, ZRKII, ZRKIN, ZRKOLD, ZRKT, ZRPK, ZSECAN, &
 & ZSECUR, ZSINA, ZSINO, ZUSKP  
REAL(KIND=JPRB) :: ZHOOK_HANDLE

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

#include "eggmlt.h"
#include "eggrvs.h"

#include "abor1.intfb.h"

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

zpis2 = prpi*0.5_jprb
zpis4 = prpi*0.25_jprb
zsecur = 1.e-12_jprb
zsecan = 1.e-05_jprb

!*
!---------------------------------------------------------------------
!     1.- PRINTING INPUT PARAMETERS
zconrad = prpi/180._jprb
zcondeg = 180._jprb/prpi

WRITE (kulout,*) ' '
WRITE (kulout,*) ' ---------------------------------- '
WRITE (kulout,*) ' '
WRITE (kulout,*) '      ARPEGE-ALADIN '
WRITE (kulout,*) ' '
WRITE (kulout,*) '    GEOGRAPHY OF GRID-POINTS '
WRITE (kulout,*) ' '
WRITE (kulout,*) ' ---------------------------------- '
WRITE (kulout,*) ' '
WRITE (kulout,*) ' INPUT PARAMETERS '
WRITE (kulout,*) ' '
WRITE (kulout,*) ' PI = ',prpi
WRITE (kulout,*) ' RADIUS OF PLANET A = ',pra*1.e-03_jprb,' KM '
WRITE (kulout,*) ' '
WRITE (kulout,*) ' X-SIZE OF ARRAYS ',kdlsur-kdlsa+1
WRITE (kulout,*) ' Y-SIZE OF ARRAYS ',kdgen-kdgsa+1
WRITE (kulout,*) ' '
WRITE (kulout,*) ' X WINDOW KDLUN = ',kdlun,' KDLUX = ',kdlux
WRITE (kulout,*) '          SIZE = ',kdlux-kdlun+1
WRITE (kulout,*) ' Y WINDOW KDGUN = ',kdgun,' KDGUX = ',kdgux
WRITE (kulout,*) '          SIZE = ',kdgux-kdgun+1
WRITE (kulout,*) ' '
WRITE (kulout,*) ' ROTATION PARAMETER KROTEQ = ',kroteq
WRITE (kulout,*) '          PLONR = ',plonr
WRITE (kulout,*) '          PLATR = ',platr
WRITE (kulout,*) ' '
WRITE (kulout,*) ' ANGLE WITH X/LATITUDE AT PLON0 = ',pbeta
WRITE (kulout,*) ' '
WRITE (kulout,*) ' SW CORNER      PLON1 = ',plon1
WRITE (kulout,*) '                PLAT1 = ',plat1
WRITE (kulout,*) ' '
WRITE (kulout,*) ' NE CORNER      PLON2 = ',plon2
WRITE (kulout,*) '                PLAT2 = ',plat2
WRITE (kulout,*) ' '
WRITE (kulout,*) ' PROJECTION PARAMETER PRPK = ',prpk
WRITE (kulout,*) ' REF POINT      PLON0 = ',plon0
WRITE (kulout,*) '                PLAT0 = ',plat0
WRITE (kulout,*) ' '
WRITE (kulout,*) ' ISOTROPY PARAMETER KSOTRP = ',ksotrp
WRITE (kulout,*) ' '
WRITE (kulout,*) ' PROJECTION REF POINT KGIV0 = ',kgiv0
WRITE (kulout,*) ' '
WRITE (kulout,*) ' ---------------------------------- '
WRITE (kulout,*) ' '

!* have a look whether we calculate plane model version
llplanex = kdlun == kdlux
llplaney = kdgun == kdgux

IF( llplanex.OR.llplaney )THEN
  isotrp = ksotrp
  ksotrp = 0
  WRITE (kulout,*) '!'
  WRITE (kulout,*) 'YOU RUN PLANE MODEL'
  WRITE (kulout,*) 'KSOTRP LOCALLY RESET TO 0 (ISOTROPIC PARAMETER) '
  WRITE (kulout,*) 'ISOTROPY PARAMETER KSOTRP = ',ksotrp
  WRITE (kulout,*) '!'
ENDIF

!     SECURITY CORRECTIONS
IF ( kroteq == 0 ) THEN
  plonr = 0.0_jprb
  platr = 0.0_jprb
ENDIF
IF ( prpk < 0.0_jprb ) THEN
  plon0 = plonr
  plat0 = platr
ENDIF
!     CHECK LONGITUDES
IF ( kroteq == 1 ) THEN
  IF ( plonr < 0.0_jprb ) THEN
    plonr = plonr + 2.0_jprb*prpi
    WRITE (kulout,*) ' *** EGGX ERROR *** WRONG CONVENTION',&
     & ' USED FOR LONGITUDES '  
    WRITE (kulout,*) ' *** NEW PLONR = ',plonr
  ENDIF
ENDIF
IF ( plon1 < 0.0_jprb ) THEN
  plon1 = plon1 + 2.0_jprb*prpi
  WRITE (kulout,*) ' *** EGGX ERROR *** WRONG CONVENTION',&
   & ' USED FOR LONGITUDES '  
  WRITE (kulout,*) ' *** NEW PLON1 = ',plon1
ENDIF
IF ( plon2 < 0.0_jprb ) THEN
  plon2 = plon2 + 2.0_jprb*prpi
  WRITE (kulout,*) ' *** EGGX ERROR *** WRONG CONVENTION',&
   & ' USED FOR LONGITUDES '  
  WRITE (kulout,*) ' *** NEW PLON2 = ',plon2
ENDIF
IF ( kgiv0 == 0.OR. kgiv0 == 2 ) THEN
  IF ( plon0 < 0.0_jprb ) THEN
    plon0 = plon0 + 2.0_jprb*prpi
    WRITE (kulout,*) ' *** EGGX ERROR *** WRONG CONVENTION',&
     & ' USED FOR LONGITUDES '  
    WRITE (kulout,*) ' *** NEW PLON0 = ',plon0
  ENDIF
ENDIF

nymggr = kroteq
xlatr = platr
xlonr = plonr
xbeta = pbeta
!*
!--------------------------------------------------------------------
!     2.- ROTATION TO EQUATOR

hsud = 1.0_jprb
IF ( kroteq == 0 ) THEN
  ! SOUTH HEMISPHERE DETECTION
  IF ( prpk == 10._jprb .OR. (prpk > 0.0_jprb .AND. prpk <= 1.0_jprb) ) THEN
    ! WHEN PROJECTION IS UNKNOWN, IF ONE OF THE REQUESTED LATITUDE IS
    ! NEGATIVE, A GUESS FOR THE OTHER MUST BE PROVIDED : SO BOTH PLAT1
    ! AND PLAT2 ARE ASSUMED KNOWN
    IF ( ksotrp == 0 ) THEN
      IF ( plat1 < 0.0_jprb .AND. plat2 < 0.0_jprb ) THEN
        hsud = -1.0_jprb
        plat1 = abs(plat1)
        plat2 = abs(plat2)
      ENDIF
    ELSEIF (ksotrp == 1 ) THEN
      IF ( plat1 < 0.0_jprb ) THEN
        IF ( zpis2 >= plat2.AND. plat2 >= -zpis2 ) THEN
          IF ( plat2 < 0.0_jprb ) THEN
            hsud = -1.0_jprb
            plat1 = abs(plat1)
            plat2 = abs(plat2)
          ELSEIF ( prpk == 10._jprb .AND. plat2 >= 0.0_jprb ) THEN
            prpk = 0.0_jprb
          ELSE
            WRITE (kulout,*) ' *** EGGX WARNING *** ',&
             & ' YOU SHOULD USE MERCATOR PROJECTION '  
          ENDIF
        ELSEIF ( prpk /= 10._jprb ) THEN
          hsud = -1.0_jprb
          plat1 = abs(plat1)
        ELSE
          WRITE (kulout,*) ' *** EGGX ERROR *** ',&
           & ' REFERENCE POLE (HSUD) CANNOT BE DECIDED '  
          WRITE (kulout,*) ' RERUN WITH A REASONABLE GUESS ',' FOR PLAT2 '
          CALL abor1(' EGGX: abor1 2.1')
        ENDIF
      ENDIF

    ELSEIF ( ksotrp == 2 ) THEN
      IF ( plat2 < 0.0_jprb ) THEN
        IF ( zpis2 >= plat1.AND. plat1 >= -zpis2 ) THEN
          IF ( plat1 < 0.0_jprb ) THEN
            hsud = -1.0_jprb
            plat1 = abs(plat1)
            plat2 = abs(plat2)
          ELSEIF ( plat1 >= 0.0_jprb ) THEN
            WRITE (kulout,*) ' *** EGGX ERROR *** ',&
             & ' PLAT1 CANNOT BE GREATER THAN PLAT2 '  
            WRITE (kulout,*) ' RERUN WITH A REASONABLE GUESS ',' FOR PLAT1 '
            CALL abor1(' EGGX: abor1 2.2')
          ENDIF
        ELSEIF ( prpk /= 10._jprb ) THEN
          hsud = -1.0_jprb
          plat2 = abs(plat2)
        ELSE
          WRITE (kulout,*) ' *** EGGX ERROR *** ',&
           & ' REFERENCE POLE (HSUD) CANNOT BE DECIDED '  
          WRITE (kulout,*) ' RERUN WITH A REASONABLE GUESS ',' FOR PLAT1 '
          CALL abor1(' EGGX: abor1 2.3')
        ENDIF
      ENDIF
    ENDIF

  ENDIF

  IF ( hsud < 0.0_jprb ) THEN
    plat0 = abs( plat0 )
  ENDIF

  ! NO ROTATION
  xlon1r =mod(plon1,2.0_jprb*prpi)
  IF ( ksotrp /= 2 ) xlat1r = plat1
  xlon2r = mod(plon2,2.0_jprb*prpi)
  IF ( ksotrp /= 1 ) zlat2r = plat2
  IF ( kgiv0 == 0.OR. kgiv0 == 2 ) xlon0r = plon0
  IF ( kgiv0 == 0.OR. kgiv0 == 1 ) xlat0r = plat0

ELSE
  ! ROTATION
  IF ( ksotrp == 0.OR. ksotrp == 1 ) THEN
    ! ROTATION OF SW CORNER
    ! ----------------------
    zsina = cos( platr )*sin( plat1 )&
     & - sin( platr )*cos( plat1 )*cos( plon1-plonr )  
    xlat1r = asin( zsina )
    IF ( abs( xlat1r ) >= zpis2 ) THEN
      xlon1r = 0.0_jprb
    ELSE
      zcosa = cos( xlat1r )
      zcoso = ( sin( platr )*sin( plat1 ) +&
       & cos( platr )*cos( plat1 )*cos( plon1-plonr ) )/zcosa  
      zcoso = min(1.0_jprb,max(-1.0_jprb,zcoso))
      zsino = ( cos( plat1 )*sin( plon1-plonr ) )/zcosa
      zsino = min(1.0_jprb,max(-1.0_jprb,zsino))
      xlon1r = acos( zcoso )
      IF ( asin( zsino ) < 0.0_jprb ) xlon1r = 2.0_jprb*prpi - xlon1r
    ENDIF
  ENDIF

  IF ( ksotrp == 0.OR. ksotrp == 2 ) THEN
    ! ROTATION OF NE CORNER
    ! ----------------------
    zsina = cos( platr )*sin( plat2 )&
     & - sin( platr )*cos( plat2 )*cos( plon2-plonr )  
    zlat2r = asin( zsina )
    IF ( abs( zlat2r ) >= zpis2 ) THEN
      xlon2r = 0.0_jprb
    ELSE
      zcosa = cos( zlat2r )
      zcoso = ( sin( platr )*sin( plat2 ) +&
       & cos( platr )*cos( plat2 )*cos( plon2-plonr ) )/zcosa  
      zcoso = min(1.0_jprb,max(-1.0_jprb,zcoso))
      zsino = ( cos( plat2 )*sin( plon2-plonr ) )/zcosa
      zsino = min(1.0_jprb,max(-1.0_jprb,zsino))
      xlon2r = acos( zcoso )
      IF ( asin( zsino ) < 0.0_jprb ) xlon2r = 2.0_jprb*prpi - xlon2r
    ENDIF
  ENDIF

  !     ROTATION OF PROJECTION REFERENCE POINT
  !     --------------------------------------
  IF ( kgiv0 == 2 ) THEN
    WRITE (kulout,*) ' *** EGGX ERROR '
    WRITE (kulout,*) ' KGIV0 = 2 IMPOSSIBLE WITH KROTEQ = 1 '
    CALL abor1(' EGGX: abor1 2.4')
  ENDIF
  IF ( kgiv0 == 0.OR. kgiv0 == 1 ) THEN
    zsina = cos( platr )*sin( plat0 )&
     & - sin( platr )*cos( plat0 )*cos( plon0-plonr )  
    xlat0r = asin( zsina )
  ENDIF
  IF ( kgiv0 == 0 ) THEN
    IF ( abs( xlat0r ) >= zpis2 ) THEN
      xlon0r = 0.0_jprb
    ELSE
      zcosa = max( cos( xlat0r ) , zsecur )
      zcoso = ( sin( platr )*sin( plat0 ) +&
       & cos( platr )*cos( plat0 )*cos( plon0-plonr ) )/zcosa  
      zcoso = min(1.0_jprb,max(-1.0_jprb,zcoso))
      zsino = ( cos( plat0 )*sin( plon0-plonr ) )/zcosa
      zsino = min(1.0_jprb,max(-1.0_jprb,zsino))
      xlon0r = acos( zcoso )
      IF ( asin( zsino ) < 0.0_jprb ) xlon0r = 2.0_jprb*prpi - xlon0r
    ENDIF
  ENDIF
ENDIF

IF ( kroteq /= 0 ) THEN
  IF ( prpk > 0.0_jprb ) THEN
    WRITE (kulout,*) ' *** EGGX WARNING ',&
     & ' USE OF ROTATION + NON-MERCATOR PROJECTION WILL LEAD ',&
     & ' TO UNPREDICTABLE RESULTS, ESP. IN SOUTH HEMISPHERE '  
    IF ( xlat1r < 0.0_jprb .AND. zlat2r < 0.0_jprb ) THEN
      hsud = -1.0_jprb
      xlat1r = abs( xlat1r )
      zlat2r = abs( zlat2r )
      xlat0r = abs( xlat0r )
    ENDIF
  ENDIF
ENDIF

WRITE (kulout,*) ' '
WRITE (kulout,*) ' HEMISPHERE INDICATOR HSUD = ',hsud
WRITE (kulout,*) ' '
WRITE (kulout,*) ' ROTATED COORDINATES '
WRITE (kulout,*) ' SW CORNER      XLON1R = ',xlon1r
WRITE (kulout,*) '                XLAT1R = ',xlat1r
WRITE (kulout,*) ' '
WRITE (kulout,*) ' NE CORNER      XLON2R = ',xlon2r
WRITE (kulout,*) '                ZLAT2R = ',zlat2r
WRITE (kulout,*) ' '
WRITE (kulout,*) ' REF POINT      XLON0R = ',xlon0r
WRITE (kulout,*) '                XLAT0R = ',xlat0r
WRITE (kulout,*) ' '

!*
!--------------------------------------------------------------------
!     3.- SPHERICAL GEOMETRY

IF ( prpk < 0.0_jprb ) THEN

  !   THE GRID IS MADE UP OF REGULARLY SPACED POINTS IN LONGITUDE (X)
  !   AND LATITUDE (Y) CORRDINATES

  !   PDELX AND PDELY ARE NON DIMENSIONNAL IN THAT CASE

  IF ( ksotrp == 0 ) THEN
    xlon1u = xlon1r
    zlon2u = xlon2r
    IF ( zlon2u < xlon1u ) zlon2u = 2.0_jprb*prpi + zlon2u
  ENDIF
  !   DEFINES OTHER CORNER
  IF ( ksotrp == 1 ) THEN
    xlon1u = xlon1r
    zlat2r = - xlat1r
    zlon2u = xlon1u + ( zlat2r-xlat1r )*&
     & REAL(kdlux-kdlun,jprb)/REAL(KDGUX-KDGUN,JPRB)  
    xlon2r = zlon2u
    IF ( zlon2u >= 2.0_jprb*prpi ) xlon2r = zlon2u - 2.0_jprb*prpi
  ENDIF
  IF ( ksotrp == 2 ) THEN
    zlon2u = xlon2r
    xlat1r = - zlat2r
    xlon1u = zlon2u - ( zlat2r-xlat1r )*&
     & REAL(kdlux-kdlun,jprb)/REAL(KDGUX-KDGUN,JPRB)  
    xlon1r = xlon1u
    IF ( xlon1u >= 2.0_jprb*prpi ) xlon1r = xlon1u - 2.0_jprb*prpi
  ENDIF

  IF( .NOT.llplanex.AND..NOT.llplaney )THEN
    pdelx = ( zlon2u-xlon1u )/REAL( kdlux-kdlun ,jprb)
    pdely = ( zlat2r-xlat1r )/REAL( kdgux-kdgun ,jprb)
  ELSEIF( llplanex )THEN
    pdely = ( zlat2r-xlat1r )/REAL( kdgux-kdgun )
    pdelx = pdely
  ELSEIF( llplaney )THEN
    pdelx = ( zlon2u-xlon1u )/REAL( kdlux-kdlun )
    pdely = pdelx
  ENDIF

  xipore=0.0_jprb
  xjpore=0.0_jprb
  xrpksm=prpk
  xggpk = prpk
  nymggi = 10

  DO jlat = kdgun, kdgux
    zlat = REAL(jlat-kdgun,jprb)*PDELY

    DO jlon = kdlun, kdlux
      zlon = REAL(jlon-kdlun,jprb)*PDELX
      pgelam(jlon,jlat) = zlon
      pgelat(jlon,jlat) = zlat
    ENDDO

    CALL eggrvs (prpi, pra, pdelx, pdely, kdlsur-kdlsa+1,&
     & 1, kdlux-kdlun+1, kulout,&
     & pgelam(kdlun,jlat), pgelat(kdlun,jlat), pgm(kdlun,jlat),&
     & pgnorx(kdlun,jlat), pgnory(kdlun,jlat))  
  ENDDO
ENDIF

!*
!--------------------------------------------------------------------
!     4.- MAP FACTOR ON SPHERICAL GEOMETRY

!*
!--------------------------------------------------------------------
!     5.- PREPARATION OF PROJECTION PARAMETERS

IF ( prpk >= 0.0_jprb ) THEN

  !   PROJECTION LONGITUDES ARE ALWAYS KNOWN

  !   SPECIAL LOGICAL TRUE WHEN DOMAIN ASTRIDE GREENWICH MERIDIAN
  !    LLGWH = .F. : GREENWICH OUT OF DOMAIN
  !    LLGWH = .T. : GREENWICH WITHIN DOMAIN
  llgwh = ( xlon2r < xlon1r )
  nymggwh = 0
  IF ( llgwh ) nymggwh = 1
  !   SHIFT LONGITUDES WHEN GREENWICH MERIDIAN IS WITHIN THE DOMAIN
  xlon1u = xlon1r
  zlon2u = xlon2r
  IF ( llgwh ) THEN
    zlon2u = xlon2r + 2.0_jprb*prpi
  ENDIF

  IF ( kgiv0 == 1.OR. kgiv0 == 3 ) THEN
    IF ( abs(abs(zlon2u-xlon1u)-prpi) > zsecan ) THEN
      WRITE (kulout,*) ' NORMAL LONGITUDE DIFFERENCE '
      xlon0u = 0.5_jprb*( xlon1u + zlon2u )
      xlon0r = xlon0u
    ELSE
      WRITE (kulout,*) ' LONGITUDE DIFFERENCE = PI '
      xlon0u = xlon1u + zpis4
      xlon0r = xlon0u
    ENDIF
    IF ( xlon0u >= 2.0_jprb*prpi ) xlon0r = xlon0r - 2.0_jprb*prpi
    WRITE (kulout,*) ' '
    WRITE (kulout,*) ' PROJECTION REFERENCE LONGITUDE '
    WRITE (kulout,*) ' (ON ROTATED SPHERE) LON0R = ',xlon0r
    WRITE (kulout,*) ' GREENWICH LOGICAL = ',llgwh
    WRITE (kulout,*) ' '
  ENDIF

  xlon0u = xlon0r
  IF ( llgwh .AND. xlon0r < xlon1r ) THEN
    xlon0u = xlon0r + 2.0_jprb*prpi
  ENDIF

  !   PROJECTION TYPE IS GIVEN BY THE USER
  !   ------------------------------------

  IF ( prpk <= 1.0_jprb ) THEN

    ! DETERMINES MISSING LATITUDE IF ANY
    IF ( ksotrp >= 1 ) THEN

      ! ADJUSTMENT OF ONE EXTREME LATITUDE IN ORDER TO HAVE PDELX = PDELY
      CALL eggmlt (prpi,kdlux,kdlun,kdgux,kdgun,kulout,1,&
       & prpk,xlon0u,xlon1u,zlon2u,ksotrp,xlat1r,zlat2r,&
       & hsud,xbeta)  

    ENDIF

    ! DETERMINES REFERENCE LATITUDE

    IF ( kgiv0 == 2.OR. kgiv0 == 3 ) THEN
      xlat0r = 0.5_jprb*( zlat2r + xlat1r )
      xlat0r = min(zpis2,max(-zpis2,xlat0r))
      WRITE (kulout,*) ' '
      WRITE (kulout,*) ' PROJECTION REFERENCE LATITUDE ',&
       & ' (ON ROTATED SPHERE) '  
      WRITE (kulout,*) ' LAT0R = ',xlat0r
      WRITE (kulout,*) ' '
    ENDIF

  ENDIF

  !     MAP PROJECTION TYPE AND POSSIBLY ONE CORNER MUST BE DETERMINED
  !     --------------------------------------------------------------
  IF ( prpk == 10._jprb ) THEN

    ! EXPLORE VARIATIONS OF MAP FACTOR
    ! --------------------------------

    zlatlim = prpi/6._jprb
    ! FIRST GUESS : EITHER STEREO OR MERCATOR
    IF ( ksotrp == 0 ) THEN
      zrpk = 1.0_jprb
      IF ( xlat1r < zlatlim ) zrpk = 0.0_jprb
    ELSEIF ( ksotrp == 1 ) THEN
      zrpk = 1.0_jprb
      IF ( xlat1r < zlatlim ) zrpk = 0.0_jprb
    ELSEIF ( ksotrp == 2 ) THEN
      zrpk = 1.0_jprb
      IF ( zlat2r < zlatlim ) zrpk = 0.0_jprb
    ENDIF

    ! FIRST GUESS POSSIBLE OTHER LATITUDE
    IF ( ksotrp >= 1 ) THEN
      CALL eggmlt (prpi,kdlux,kdlun,kdgux,kdgun,kulout,1,&
       & zrpk,xlon0u,xlon1u,zlon2u,ksotrp,xlat1r,zlat2r,&
       & hsud,xbeta)  
    ENDIF
    ! FIRST GUESS REFERENCE LATITUDE
    IF ( kgiv0 == 2.OR. kgiv0 == 3 ) THEN
      xlat0r = 0.5_jprb*( zlat2r + xlat1r )
      xlat0r = min(zpis2,max(-zpis2,xlat0r))
    ENDIF

    ! PREPARING FOR OUTER ITERATION LOOP (RK + LATITUDE)
    iterkx = 10
    iterk = 0
    zrkold = zrpk

    ll510=.true.

    DO WHILE(ll510)

      zrkii = 0.01_jprb
      zdrk = 0.01_jprb
      zrki = 0.0_jprb
      ispeca = 0

      ! IDENTIFYING SPECIAL CASES THAT MUST GO INTO THE INNER LOOP
      IF ( xlat1r*zlat2r  <=  0.0_jprb ) THEN
        ! SPECIAL CASE : DOMAIN ASTRIDE EQUATOR
        zrpk = 0.0_jprb
        zrkold = zrpk
        ispeca = 1
      ELSEIF ( hsud > 0.0_jprb .AND.abs(zlat2r) <= abs(xlat1r) ) THEN
        ! SPECIAL CASE : DOMAIN INCLUDING NORTH POLE (THE TEST IS A WEAK ONE)
        zrpk = 1.0_jprb
        zrkold = zrpk
        ispeca = 1
      ELSEIF ( hsud < 0.0_jprb .AND.abs(xlat1r) <= abs(zlat2r) ) THEN
        ! SPECIAL CASE : DOMAIN INCLUDING SOUTH POLE (THE TEST IS A WEAK ONE)
        zrpk = 1.0_jprb
        zrkold = zrpk
        ispeca = 1
      ELSEIF ( abs(abs(zlat2r)-zpis2) < zsecan ) THEN
        ! VARIOUS OTHER SPECIAL CASES THAT MUST NOT COVER THE WHOLE LOOP
        zrki = zrkii
      ELSEIF ( abs(abs(xlat1r)-zpis2) < zsecan ) THEN
        ! VARIOUS OTHER SPECIAL CASES THAT MUST NOT COVER THE WHOLE LOOP
        zrki = zrkii
      ENDIF

      IF ( ispeca /= 1 ) THEN

        ! PREPARE FOR INNER LOOP ON RK
        zdmmax = -1.0_jprb
        zdmmin = 1.e+05_jprb
        zrkax = 1.0_jprb
        zrkin = 1.0_jprb
        innega = 0
        inbess = 0
        zrkt = zrki

        ! INNER LOOP ON RK

        ll520=.true.

        DO WHILE(ll520)
          ! COMPUTES MAP FACTOR VARIATION FOR VARIOUS PROJECTION
          IF ( zrkt == 0.0_jprb ) THEN
            zdm = cos(xlat0r)/cos(zlat2r) - cos(xlat0r)/cos(xlat1r)
          ELSEIF ( zrkt /= 1.0_jprb ) THEN
            zdm = (cos(xlat0r)**(1.0_jprb-zrkt))*&
             & ((1.0_jprb+sin(xlat0r))**zrkt)*( (cos(xlat1r)**(zrkt-1.0_jprb))*&
             & ((1.0_jprb+sin(xlat1r))**(-zrkt)) - (cos(zlat2r)**(zrkt-1.0_jprb))*&
             & ((1.0_jprb+sin(zlat2r))**(-zrkt)) )  
          ELSEIF ( zrkt == 1.0_jprb ) THEN
            zdm = (1.0_jprb+sin(xlat0r))*( 1.0_jprb/(1.0_jprb+sin(xlat1r)) -&
             & 1.0_jprb/(1.0_jprb+sin(zlat2r)) )  
          ENDIF
          inbess = inbess + 1
          IF ( zdm <= 0.0_jprb ) THEN
            innega = innega + 1
          ENDIF
          IF ( zdm >= zdmmax .AND. zdm > 0.0_jprb ) THEN
            zdmmax = zdm
            zrkax = zrkt
          ENDIF
          IF ( zdm <= zdmmin .AND. zdm > 0.0_jprb ) THEN
            zdmmin = zdm
            zrkin = zrkt
          ENDIF
          zrkt = zrkt + zdrk
          ll520=(zrkt <= 1.0_jprb)
        ENDDO ! DO WHILE(LL520)

        WRITE (kulout,*) ' '
        WRITE (kulout,*) ' CHOICE OF OPTIMAL RK ITERATION ',iterk
        WRITE (kulout,*) ' PREVIOUS RK ',zrkold
        WRITE (kulout,*) ' TEST OVER ',inbess,' VALUES '
        WRITE (kulout,*) ' INCLUDING ',innega,' NEGATIVE VALUES'
        WRITE (kulout,*) ' RK MINI = ',zrki,' INCRMENT = ',zdrk
        WRITE (kulout,*) ' '
        WRITE (kulout,*) ' DELTA(M) MAXI = ',zdmmax,' AT RK = ',zrkax
        WRITE (kulout,*) ' DELTA(M) MINI = ',zdmmin,' AT RK = ',zrkin

        ! UPDATES VALUE OF GUESS ZRPK
        zrkold = zrpk
        zrpk = zrkin

      ENDIF ! ISPECA /= 1

      IF ( ksotrp == 0 ) THEN
        ! END THE OUTER LOOP IF BOTH LATITUDES WERE KNOWN
        zrkold = zrpk
      ELSEIF ( ksotrp >= 1 ) THEN
        ! DETRMINES NEW OTHER LATITUDE
        CALL eggmlt (prpi,kdlux,kdlun,kdgux,kdgun,kulout,1,&
         & zrpk,xlon0u,xlon1u,zlon2u,ksotrp,xlat1r,zlat2r,&
         & hsud,xbeta)  
      ENDIF
      ! NEW REFERENCE LATITUDE
      IF ( kgiv0 == 2.OR. kgiv0 == 3 ) THEN
        xlat0r = 0.5_jprb*( zlat2r + xlat1r )
        xlat0r = min(zpis2,max(-zpis2,xlat0r))
      ENDIF

      ! COMPLETES INTERATION OF OUTER LOOP (RK AND LATITUDE)
      iterk = iterk + 1
      IF ( iterk > iterkx ) THEN
        WRITE (kulout,*) ' *** EGGX **** TROUBLE '
        WRITE (kulout,*) ' NO CONVERGENCE OF AUTOMATIC CHOICE '
        CALL abor1(' EGGX: abor1 5.1')
      ENDIF

      ll510=(abs(zrpk-zrkold) > zdrk)

    ENDDO ! DO WHILE(LL510)

    prpk = zrpk
    WRITE (kulout,*) ' '
    WRITE (kulout,*) ' --- EGGX AUTOMATIC CHOICE '
    WRITE (kulout,*) ' '
    WRITE (kulout,*) ' FINAL VALUE OF PRPK = ',prpk
    WRITE (kulout,*) ' '
    WRITE (kulout,*) ' PROJECTION REFERENCE LATITUDE ',&
     & ' (ON ROTATED SPHERE) '  
    WRITE (kulout,*) ' FINAL LATITUDE LAT1 R = ',xlat1r
    WRITE (kulout,*) ' FINAL LATITUDE LAT2 R = ',zlat2r
    WRITE (kulout,*) ' LAT0R = ',xlat0r
    WRITE (kulout,*) ' '
  ENDIF

ENDIF
!     MEMORIZES THE FINAL VALUE OF PRPK
xggpk = prpk

!*
!---------------------------------------------------------------------
!     6.- PROJECTION ON CARTESIAN PLAN

!     REMARK : THE CODE IS GENERAL, BUT STEREO/LAMBERT SHOULD
!              NOT BE COMBINED TO ROTATION

IF ( prpk > 1.0_jprb ) THEN
  WRITE (kulout,*) ' *** EGGX ERROR : NON-EXISTING PROJ.'
  CALL abor1(' EGGX: abor1 6.1')
ENDIF

IF ( prpk > 0.0_jprb ) THEN
  WRITE (kulout,*) ' STEREO OR LAMBERT PROJECTION '
  IF ( prpk == 1.0_jprb ) WRITE (kulout,*)&
   & ' EFFECTIVELY STREOGRAPHIC PROJECTION '  

  !   COMPUTES BASIC PARAMETERS
  !   -------------------------

  !    HALF COLATITUDE
  zdcla1 = zpis4 - 0.5_jprb*xlat1r
  zdcla2 = zpis4 - 0.5_jprb*zlat2r
  zdcla0 = zpis4 - 0.5_jprb*xlat0r
  !   PROJECTION CONSTANTS
  IF ( prpk < 1.0_jprb ) THEN
    xggm0 = ( cos( xlat0r )**(1.0_jprb-prpk) )*&
     & ( ( 1.0_jprb + sin( xlat0r ) )**prpk )  
    zuskp = xggm0/prpk
    xrpksm = 1.0_jprb/zuskp
  ELSE
    xggm0 = 1.0_jprb + sin( xlat0r )
    zuskp = xggm0
    xrpksm = 1.0_jprb/zuskp
  ENDIF

  !   COMPUTES RESOLUTION
  !   -------------------
  IF( .NOT.llplanex.AND..NOT.llplaney )THEN
    pdelx = pra*zuskp*( (tan(zdcla2)**prpk)*sin( prpk*(&
     & zlon2u-xlon0u)-xbeta ) - (tan(zdcla1)**prpk)*&
     & sin( prpk*(xlon1u-xlon0u)-xbeta ) )/REAL(KDLUX-KDLUN,JPRB)  
    pdely = hsud*pra*zuskp*( (tan(zdcla1)**prpk)*cos( prpk*(&
     & xlon1u-xlon0u)-xbeta ) - (tan(zdcla2)**prpk)*&
     & cos( prpk*(zlon2u-xlon0u)-xbeta ) )/REAL(KDGUX-KDGUN,JPRB)  
  ELSEIF( llplanex )THEN
    pdely = hsud*pra*zuskp*( (tan(zdcla1)**prpk)*cos( prpk*(&
     & xlon1u-xlon0u)-xbeta ) - (tan(zdcla2)**prpk)*&
     & cos( prpk*(zlon2u-xlon0u)-xbeta ) )/REAL(KDGUX-KDGUN)  
    pdelx=pdely
  ELSEIF( llplaney )THEN
    pdelx = pra*zuskp*( (tan(zdcla2)**prpk)*sin( prpk*(&
     & zlon2u-xlon0u)-xbeta ) - (tan(zdcla1)**prpk)*&
     & sin( prpk*(xlon1u-xlon0u)-xbeta ) )/REAL(KDLUX-KDLUN)  
    pdely = pdelx
  ENDIF
      
  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' MAP FACTOR BASE XGGM0 = ',xggm0
  WRITE (kulout,*) ' ZUSKP = ',zuskp
  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' X GRID SIZE (KM) = ',pdelx*1.e-03_jprb
  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' Y GRID SIZE (KM) = ',pdely*1.e-03_jprb
  WRITE (kulout,*) ' '

  !   COMPUTES POLE LOCATION ON GRID
  !   ------------------------------

  xipore =  - pra*zuskp*( tan( zdcla1 )**prpk )*&
   & sin( prpk*(xlon1u-xlon0u)-xbeta )/pdelx  
  xjpore =  hsud*pra*zuskp*( tan( zdcla1 )**prpk )*&
   & cos( prpk*(xlon1u-xlon0u)-xbeta )/pdely  
  nymggi = 10

  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' POLE LOCATION ON GRID IP = ',xipore,' JP = ',xjpore

  zipv = REAL( KDLUX ,JPRB) - PRA*ZUSKP*( tan( zdcla2 )**prpk )*&
   & SIN( PRPK*(ZLON2U-XLON0U)-XBETA )/PDELX  - REAL(KDLUN,JPRB)  
  ZJPV = real( kdgux ,jprb) + pra*zuskp*( tan( zdcla2 )**prpk )*&
   & cos( prpk*(zlon2u-xlon0u)-xbeta )*hsud/pdely -&
   & REAL(KDGUN,JPRB)  

  WRITE (kulout,*) ' VRF POLE LOCATION ON GRID IP = ',zipv,' JP = ',zjpv

  !   GRID POINTS LOCATION
  !   --------------------

  DO jlat = kdgun, kdgux

    DO jlon = kdlun, kdlux
      pgelam(jlon,jlat) = REAL(jlon-kdlun,jprb)*PDELX
      pgelat(jlon,jlat) = REAL(jlat-kdgun,jprb)*PDELY
    ENDDO

    CALL eggrvs (prpi, pra, pdelx, pdely, kdlsur-kdlsa+1,&
     & 1, kdlux-kdlun+1, kulout,&
     & pgelam(kdlun,jlat), pgelat(kdlun,jlat), pgm(kdlun,jlat),&
     & pgnorx(kdlun,jlat), pgnory(kdlun,jlat))  
  ENDDO

ENDIF

IF ( prpk == 0.0_jprb ) THEN
  WRITE (kulout,*) ' MERCATOR PROJECTION '

  !   COMPUTES BASIC PARAMETERS
  !   -------------------------

  !    HALF COLATITUDE
  zdcla1 = zpis4 - 0.5_jprb*xlat1r
  zdcla2 = zpis4 - 0.5_jprb*zlat2r
  zdcla0 = zpis4 - 0.5_jprb*xlat0r

  !   COMPUTES RESOLUTION
  !   -------------------

  zface = pra*cos( xlat0r )
  zdlon =  zlon2u - xlon1u
  zdtlat = log( tan(zdcla1)/tan(zdcla2) )
  IF( .NOT.llplanex.AND..NOT.llplaney )THEN
    pdelx = zface*( zdlon*cos(xbeta) + zdtlat*sin(xbeta) )&
     & /REAL( KDLUX-KDLUN ,JPRB)  
    pdely = zface*( -zdlon*sin(xbeta) + zdtlat*cos(xbeta) )&
     & /REAL( KDGUX-KDGUN ,JPRB)  
  ELSEIF( llplanex )THEN
    pdely = zface*( -zdlon*sin(xbeta) + zdtlat*cos(xbeta) )&
     & /REAL( KDGUX-KDGUN )  
    pdelx = pdely
  ELSEIF( llplaney )THEN
    pdelx = zface*( zdlon*cos(xbeta) + zdtlat*sin(xbeta) )&
     & /REAL( KDLUX-KDLUN )  
    pdely = pdelx
  ENDIF

  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' MAP FACTOR BASE COS( LAT0 ) = ',cos( xlat0r )
  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' X GRID SIZE (KM) = ',pdelx*1.e-03_jprb
  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' Y GRID SIZE (KM) = ',pdely*1.e-03_jprb
  WRITE (kulout,*) ' '

  !   COMPUTES EQUATOR LOCATION ON GRID
  !   ------------------------------

  xipore = - pra*cos( xlat0r )*( xlon1u-xlon0u )/pdelx
  xjpore = + pra*cos( xlat0r )*log(tan( zdcla1 ) )/pdely
  xrpksm=1.0_jprb
  xggpk = prpk
  nymggi = 10

  WRITE (kulout,*) ' '
  WRITE (kulout,*) ' EQUATOR LOCATION ON GRID IE = ',xipore,' JE = ',xjpore

  !     GRID POINTS LOCATION
  !     --------------------

  DO jlat = kdgun, kdgux

    DO jlon = kdlun, kdlux
      pgelam(jlon,jlat) = REAL(jlon-kdlun,jprb)*PDELX
      pgelat(jlon,jlat) = REAL(jlat-kdgun,jprb)*PDELY
    ENDDO

    CALL eggrvs (prpi, pra, pdelx, pdely, kdlsur-kdlsa+1,&
     & 1, kdlux-kdlun+1, kulout,&
     & pgelam(kdlun,jlat), pgelat(kdlun,jlat), pgm(kdlun,jlat),&
     & pgnorx(kdlun,jlat), pgnory(kdlun,jlat))  
  ENDDO

ENDIF

!*
!---------------------------------------------------------------------
!     7.- INVERSE ROTATION BACK TO GEOGRAPHICAL COORDINATES

!    THIS OPERATION IS PERFORMED BY EGGRVS

!    FINAL UPDATE OF ACTUAL CORNERS USED

IF ( ksotrp == 1 ) THEN
  plat2 = pgelat(kdlux,kdgux)
  plon2 = pgelam(kdlux,kdgux)
ENDIF
IF ( ksotrp == 2 ) THEN
  plat1 = pgelat(kdlun,kdgun)
  plon1 = pgelam(kdlun,kdgun)
ENDIF

!    PLON1 AND PLON2 MUST BE BETWEEN 0 AND 2*RPI

plon1=mod(plon1,2*prpi)
plon2=mod(plon2,2*prpi)

IF ( hsud < 0.0_jprb ) THEN
  IF ( ksotrp /= 2 ) THEN
    plat1 = hsud*plat1
  ENDIF
  IF ( ksotrp /= 1 ) THEN
    plat2 = hsud*plat2
  ENDIF
  plat0 = hsud*plat0
ENDIF

! KSOTRP reset back to original value
IF( llplanex.OR.llplaney )THEN
  ksotrp = isotrp
ENDIF

WRITE (kulout,*) ' '
WRITE (kulout,*) ' ---------- '
WRITE (kulout,*) ' '
WRITE (kulout,*) '  EGGX IS OVER '
WRITE (kulout,*) ' '
WRITE (kulout,*) ' ---------- '
WRITE (kulout,*) ' '
WRITE (kulout,*) ' '

!---------------------------------------------------------------------
IF (lhook) CALL dr_hook('EGGX',1,zhook_handle)
END SUBROUTINE eggx