eggpack.F90

Go to the documentation of this file.

MODULE eggpack

! Version 2009.0317 by JD GRIL

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DOC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! This package contains many usefull functions about Projection in "Tangent 
!Case" for Lambert Conic Comformal, Polar Stereographic and Mercator. A Domain
!Maker subroutine for grid points domain defined by its center is present.
! * The functions are classed in three types :
!  - Independants functions
!  - Specifics functions -| theses functions work for single point data structure
!  - Generics functions  -| or for array of single point data structure (rank of
!                         | array is one)
! * One subroutine to make domain

! -1- Decription of Functions :

!  -1.1- Independants functions : 

!    Theses functions have not direct relations with projection but are usefull 
!    for this package.

!->logical function RETURN_PRINT(CODE_ERR,NUM_TEST,AUTO_STOP) 

!     Return a logical to know if in case of error the the program will stop
!     or return an error code to the caller depending AUTO_STOP logical flag.
!     In every cases it prints informations (ERROR,WARNING,INFO,OK) contained in
!     CODE_ERR structure of type ERROR and NUM_TEST. It's used in subroutine
!     MAKDO. See how inside. 

!->function TYPE_PROJ(REF_COORD)

!     Return 1 character (M,L,S) depending the coordinates of reference point
!     REF_COORD structure of type LOLA (in Degrees).

!->real function POLE_IS (REF_COORD) 

!     Return pole location (South is -1.0, North is 1.0) for Polar Stereographic or
!     Lambert projections (return 0.0 in Mercator : no sense) depending the coor-
!     dinates of reference point REF_COORD structure of type LOLA (in Degrees).

!  -1.2- Specifics functions : 

!    These functions are only used in Polar Stereographic or Lambert projection
!    type. Depending the type of arguments (arrays or not) the generic caller-
!    name defined by interface selects the good one (suffixed by _V for arrays
!     or _S for scalars).
!    All of these functions are prefixed by STPL_; they use P_PJ structure of
!    type PARAM_PROJ (the most important of this package) that defines the 
!    parameters of the projection depending the reference point, see below.

!->type (RTETA) function STPL_LATLON_TO_RTETA(PT_COORD,P_PJ,PI) result (PT_RTETA)

!->type (LOLA)  function STLP_RTETA_TO_LATLON(PT_RTETA,P_PJ,PI) result (PT_COORD)

!->type (XY)    function STLP_RTETA_TO_XY(PT_RTETA,P_PJ)        result (PT_XY)

!->type (RTETA) function STLP_XY_TO_RTETA(PT_XY,P_PJ,PI)        result (PT_RTETA)

!     These functions convert coordinates from LAT-LON (in Radians) to R-THETA
!     (in meters, radians) and from R-THETA to XY (in meters) and reverse. PT_???
!     are structures or arrays of structures. PT_XY coordinates are referenced by
!     STD (standard) origin that is projection of pole, the X axis is West to
!     East, the Y axis is South to North. PT_RTETA coordinates are referenced by
!     STD (standard) origin that is projection of pole, R is the distance to it,
!     TETA is angle with reference point longitude (negative in clockwise).
!            Remark -> PT_XY and PT_RTETA are on projection plane !

!  -1.3- Generics functions :

!    These functions are standards for all types of projections (M,L,S). the most
!    important is :

!->type (PARAM_PROJ) function REF_DATAS (REF_COORD,RA,TOZERO_COORD,LRT) 
!  result (P_P) 

!      It defines a P_P structure of type PARAM_PROJ that creates all the
!    parameters needs by all specifics and generics functions : this structure
!    defines the type of projection depending the reference point coordinates
!    (in Degrees, inside you have the reference point coordinates in Radians).
!      For all others functions, the generic caller-name selects the good one
!    (suffixed by _V for arrays or _S for scalars) depending the type of 
!    arguments (arrays or not). they are three classes of functions :

!    a) changing origin of XY type of coordinates :

!     In Mercator projection the YX origin is the reference point, because 
!     projection of pole is non sense. 

!->type (XY) function XY_NEW_TO_STD_ORIGIN(NEW_ORIGIN_COORD,PT_XY_IN_NEW_ORIGIN,
!  P_PJ,PI) result (PT_XY_IN_STD_ORIGIN) 

!     This function changes coordinates PT_XY_IN_NEW_ORIGIN (may be an array)
!     relative to NEW_ORIGIN_COORD (in Degrees) to STD origin (depending of 
!     projection).

!->type (XY) function XY_STD_TO_NEW_ORIGIN(NEW_ORIGIN_COORD,PT_XY_IN_STD_ORIGIN,
!  P_PJ,PI) result (PT_XY_IN_NEW_ORIGIN)

!     This function changes coordinates PT_XY_IN_STD_ORIGIN (may be an array)
!     relative to STD origin (depending of projection) to NEW_ORIGIN_COORD 
!     (in Degrees).

!            Remark -> To pass PT_XY_IN_OLD_ORIGIN to PT_XY_IN_NEW_ORIGIN, you
!                      can do (or create your own function) :
!       PT_XY_IN_NEW_ORIGIN =                                                   &
!           XY_STD_TO_NEW_ORIGIN(NEW_ORIGIN_COORD,                              &
!           XY_NEW_TO_STD_ORIGIN(OLD_ORIGIN_COORD,PT_XY_IN_OLD_ORIGIN,P_PJ,PI), &
!                                                 P_PJ,PI)  

!    b) converting LAT-LON coordinates in XY and reverse :

!->type (XY)   function LATLON_TO_XY(PT_COORD,P_PJ,PI) result (PT_XY)

!->type (LOLA) function XY_TO_LATLON(PT_XY,P_PJ,PI)    result (PT_COORD)

!     These functions convert coordinates from LAT-LON (in Radians) to XY (in
!     meters) and reverse. PT_??? are structures or arrays of structures. PT_XY
!     coordinates are referenced by STD (standard) origin that is projection of
!     pole, for Lambert or Polar Stereographic projection or reference point for
!     Mercator; the X axis is West to East, the Y axis is South to North.
!            Remark -> PT_XY is on projection plane !

!    c) Map Factor,  projection of Geographical North : 

!->real function MAP_FACTOR(PT_COORD,P_PJ,PI,RA)

!     Compute Map Factor for PT_COORD (in Radians) (array of) structure of type
!     LOLA depending type of projection.

!->type (PGN) function GN(PT_COORD,P_PJ)  

!     Compute Projection of Geographical North unit vector on XY for PT_COORD (in
!     Radians) (array of) structure of type LOLA depending type of projection.

! -2- Decription of Subroutine :       

!  Subroutine MAKDO (make domain) use CENTER point of the Domain as the Origin 
!  of XY points.

!->subroutine MAKDO(REF_COORD,CENTER_COORD,PDEL,NB_PTS,GRID_COORD,GRID_MF,
!  GRID_PGN,GRID_INFO,ERR_CODE,AUTO_STOP,PI,RA,LMRT)

!   This subroutine creates a grid domain usefull datas depending projection.
!  It makes validations tests, using functions of this package computes
!  outputs.
!   The X and Y points are computed from the center in a 2 dimensional array
!  GRID_XY_C structure of type XY
!   To pass 2 dimensional arrays to one dimensional in functions we must
!  use the pack F90 function and unpack for the reverse because all the
!  functions are defined in interfaces.

!  - Inputs-Outputs :
!    REF_COORD, CENTER_COORD : in Degrees in input, in Radians in output
!                              type LOLA
!                              defining Reference coordinates for the
!                              projection and center coordinates of domain.
!  - Inputs :
!    PDEL                    : resolution in meters of the grid in projection
!                              plane.
!                              type DELTA
!    NB_PTS                  : number of points on X and Y
!                              type NBPTS
!  - Optionals Inputs :
!    PI, RA                  : Pi and Raduis of Earth
!                              type real
!    AUTO_STOP               : specify if ERROR stops the program (default)
!                              or return ERR_CODE to caller
!                              type logical
!  - Outputs :
!    ERR_CODE                : error code returned
!                              type ERROR
!    GRID_INFO               : informations on grid domain (projections,
!                              reference point, center, corners, map factors)
!                              type DOMI
!    GRID_COORD              : grid coordinates LAT [-90.0,90.0] in Radian
!                                               LON [0.0,360.0[ in Radians
!                              array dim = 2 : NB_PTS%ONX,NB_PTS%ONY
!                              type LOLA
!    GRID_MF                 : grid map factor
!                              array dim = 2 : NB_PTS%ONX,NB_PTS%ONY
!                              type real
!    GRID_PGN                : projection of Geographical North unit vector
!                              on X and Y
!                              array dim = 2 : NB_PTS%ONX,NB_PTS%ONY
!                              type PGN

!            Remark -> the EBETA parameter in old EGGX is fully automatic
!                      for Lambert and Polar Stereographic : it depends only
!                      of reference longitude. In mercator, it has non sense
!                      because alls meridians are parallels to reference
!                      longitude.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Author : Jean-Daniel GRIL , CNRM/GMAP/EXT , March 28 2002
! Modified 04-08-19 R. El Khatib & J.-D. Gril : cleanups
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!  Modifications :
!  _______________

!  November 09, 2004 : JD GRIL
!  --------------------------
!      - supress of many "#define" values
!      - add 2 news optionals arguments in REF_DATAS :
!            - TOZERO_COORD coordinates of domain's center that is placed in (0,0)
!              coord. if Rotated-Tilted Mercator Projection is used
!            - LRT flag to .TRUE. if you use Rotated-Tilted Mercator Projection
!            The use of this flag defines a "W" value for P_P%TYPE_PJ (projection
!            type) and the use of EGGMRT module. The functions and MAKDO routine 
!            was modified according this new feature. 
!      - MAKDO has a new optional parameter to define Rotated-Tilted Mercator
!        Projection :
!            - LD_LMRT 

!      To use Rotated-Tilted Mercator projection, remember that reference coordinates
!      YD_REF_COORD%LAT = 0 (mercator) and YD_REF_COORD%LON is the tilting angle
!      of the domain (at the point of coordinate (0,0)),
!      after the rotation of the center of the domain (YD_CENTER_COORD coordinates)
!      to the point of coordinates (0,0). 
!      All use of functions is unchanged (the new description of projection is
!      done by the P_P argument).

!      - add some protections in Lambert (STLP_XY_TO_RTETA functions) to prevent
!        the nosense of XY values in the cutting part of Lambert projected plane.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!  Modified 11-07-2006 : JD Gril : use of LOLAD & LOLAR functions vs /DTR & *DTR
!  Modified 16-10-2006 : JD Gril : replace pack/unpack by reshape intrinsic func.
!  Modified 12-01-2007 : JD Gril : cleaning
!  Modified 04-07-2008 : JD Gril : latlon_to_xy and stpl_latlon_to_rteta edge
!                                : effect for huge domain : pb diff longitudes.
!                                : replace reshape by do loop
!  Modified 13-03-2009 : JD Gril : Optimization et cleanup
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

! ******************* Definition of parameters **********************************

!     Include Kinds
!     -------------

#ifndef DEBUG
! force no print info in makdo by default
#define _DEFP_ .FALSE.
#else
! force print info in makdo by default
#define _DEFP_ .TRUE.
#endif

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

! ******************* Loading module ********************************************

USE eggangles  ,ONLY : lola, val_coord, lolar, lolad, angle_domain, dist_2ref
USE eggmrt  ,ONLY : merotil, metilrot 

IMPLICIT NONE

! Physics constant
! ----------------

REAL(KIND=JPRD), PARAMETER, PUBLIC :: r_earth = 6371229._jprd

! ******************* Definition of type ****************************************

TYPE error
  INTEGER(KIND=JPIM) :: num
  CHARACTER(LEN=100) :: txt
END TYPE error

TYPE pgn
  REAL(KIND=JPRD) :: onx, ony
END TYPE pgn

TYPE nbpts
  INTEGER(KIND=JPIM) :: onx, ony
END TYPE nbpts

TYPE delta
  REAL(KIND=JPRD) :: onx, ony
END TYPE delta

TYPE rteta
  REAL(KIND=JPRD) :: r, teta
END TYPE rteta

TYPE xy
  REAL(KIND=JPRD) :: x, y
END TYPE xy

TYPE param_proj
  sequence
  type(lola) :: ref_pt,tzo_pt
  REAL(KIND=JPRD) :: kl, r_equateur, pole
  CHARACTER(LEN=1) :: type_pj
END TYPE param_proj

TYPE domi
  type(nbpts) :: g_size
  type(lola) :: ct_coord, rf_coord, sw_coord, se_coord, ne_coord, nw_coord
  REAL(KIND=JPRD) :: mf_ct, mf_rf, mf_sw, mf_se, mf_ne, mf_nw
  type(param_proj) :: info_proj
END TYPE domi

! ******************* Definition of Interface ***********************************

INTERFACE info_print
  MODULE PROCEDURE info_pp_print, info_domi_print
END INTERFACE

INTERFACE xy_new_to_std_origin
  MODULE PROCEDURE xy_new_to_std_origin_v, xy_new_to_std_origin_s
END INTERFACE

INTERFACE xy_std_to_new_origin
  MODULE PROCEDURE xy_std_to_new_origin_v, xy_std_to_new_origin_s
END INTERFACE

INTERFACE stpl_latlon_to_rteta
  MODULE PROCEDURE stpl_latlon_to_rteta_v, stpl_latlon_to_rteta_s
END INTERFACE

INTERFACE stlp_rteta_to_xy
  MODULE PROCEDURE stlp_rteta_to_xy_v, stlp_rteta_to_xy_s
END INTERFACE

INTERFACE stlp_xy_to_rteta
  MODULE PROCEDURE stlp_xy_to_rteta_v, stlp_xy_to_rteta_s
END INTERFACE

INTERFACE stlp_rteta_to_latlon
  MODULE PROCEDURE stlp_rteta_to_latlon_v, stlp_rteta_to_latlon_s
END INTERFACE

INTERFACE latlon_to_xy
  MODULE PROCEDURE latlon_to_xy_v, latlon_to_xy_s
END INTERFACE

INTERFACE xy_to_latlon
  MODULE PROCEDURE xy_to_latlon_v, xy_to_latlon_s
END INTERFACE

INTERFACE map_factor
  MODULE PROCEDURE map_factor_v, map_factor_s
END INTERFACE

INTERFACE gn
  MODULE PROCEDURE gn_v, gn_s
END INTERFACE

#include "abor1.intfb.h"

CONTAINS

! =================== FUNCTIONS =================================================

! ******************* Independants functions ************************************

! -------------------------------------------------------------------------------
SUBROUTINE info_domi_print(YD_G_INFO,KOUT,PI)
type(domi), INTENT(IN)                  :: yd_g_info
INTEGER(KIND=JPIM), INTENT(IN), OPTIONAL :: KOUT
REAL(KIND=JPRD), INTENT(IN), OPTIONAL    :: PI

REAL(KIND=JPRD)      :: TPI
INTEGER(KIND=JPIM)   :: TKOUT
REAL(KIND=JPRD)      :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:INFO_DOMI_PRINT',0,zhook_handle)
IF (PRESENT(kout))THEN
  tkout = kout
ELSE
  tkout = 6_jpim
ENDIF
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
CALL info_print(yd_g_info%INFO_PROJ,tkout,tpi)
WRITE(tkout,*) "============================================================="
WRITE(tkout,*) "===   Informations about Domain Information Structure    ===="
WRITE(tkout,*) "============================================================="
WRITE(tkout,*)
WRITE(tkout,*) "  -Size of Domain (in points) :"
WRITE(tkout,'(13X,A7,10X,A7)') " On X  "," On Y  "
WRITE(tkout,'(13X,I7,10X,I7)') yd_g_info%G_SIZE%ONX,yd_g_info%G_SIZE%ONY
WRITE(tkout,*) "  -Most important points informations :"
WRITE(tkout,'(1X,A7,1X,"|",1X,A9,1X,"|",1X,A8,1X,"|",1X,A16)') &
 & " Points ","Longitude","Latitude","   Map Factor   "
WRITE(tkout,'(A47)') "------------------------------------------------"
WRITE(tkout,'(1X,A7,1X,"|",2X,F7.2,2X,"|",2X,F7.2,1X,"|",1X,G16.10)') &
 & "Center ",yd_g_info%CT_COORD%LON,yd_g_info%CT_COORD%LAT,yd_g_info%MF_CT
WRITE(tkout,'(1X,A7,1X,"|",2X,F7.2,2X,"|",2X,F7.2,1X,"|",1X,G16.10)') &
 & "Refer. ",yd_g_info%RF_COORD%LON,yd_g_info%RF_COORD%LAT,yd_g_info%MF_RF
WRITE(tkout,'(1X,A7,1X,"|",2X,F7.2,2X,"|",2X,F7.2,1X,"|",1X,G16.10)') &
 & "S.West ",yd_g_info%SW_COORD%LON,yd_g_info%SW_COORD%LAT,yd_g_info%MF_SW
WRITE(tkout,'(1X,A7,1X,"|",2X,F7.2,2X,"|",2X,F7.2,1X,"|",1X,G16.10)') &
 & "S.East ",yd_g_info%SE_COORD%LON,yd_g_info%SE_COORD%LAT,yd_g_info%MF_SE
WRITE(tkout,'(1X,A7,1X,"|",2X,F7.2,2X,"|",2X,F7.2,1X,"|",1X,G16.10)') &
 & "N.East ",yd_g_info%NE_COORD%LON,yd_g_info%NE_COORD%LAT,yd_g_info%MF_NE
WRITE(tkout,'(1X,A7,1X,"|",2X,F7.2,2X,"|",2X,F7.2,1X,"|",1X,G16.10)') &
 & "N.West ",yd_g_info%NW_COORD%LON,yd_g_info%NW_COORD%LAT,yd_g_info%MF_NW
WRITE(tkout,*) "============================================================="
IF (lhook) CALL dr_hook('EGGPACK:INFO_DOMI_PRINT',1,zhook_handle)
END SUBROUTINE info_domi_print
! -------------------------------------------------------------------------------
SUBROUTINE info_pp_print(P_P,KOUT,PI)
type(param_proj), INTENT(IN)                :: p_p
INTEGER(KIND=JPIM), INTENT(IN), OPTIONAL     :: KOUT
REAL(KIND=JPRD), INTENT(IN), OPTIONAL        :: PI

REAL(KIND=JPRD)    :: TPI, DTR
INTEGER(KIND=JPIM) :: TKOUT
REAL(KIND=JPRD)    :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:INFO_PP_PRINT',0,zhook_handle)
IF (PRESENT(kout))THEN
  tkout = kout
ELSE
  tkout = 6_jpim
ENDIF
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
dtr = tpi/180.0_jprd
WRITE(tkout,*) "============================================================="
WRITE(tkout,*) "===  Informations about Parameters Projection Structure  ===="
WRITE(tkout,*) "============================================================="
WRITE(tkout,*)
WRITE(tkout,*) "  -Reference Point Coordinates :"
WRITE(tkout,'(13X,A7,10X,A7)') "Degrees","Radians"
WRITE(tkout,'(1X,"Longitude : ",F7.2,5X,G16.10)') p_p%REF_PT%LON/dtr,p_p%REF_PT%LON
WRITE(tkout,'(1X,"Latitude  : ",F7.2,5X,G16.10)') p_p%REF_PT%LAT/dtr,p_p%REF_PT%LAT
WRITE(tkout,*) "  -Projection Characteristics  :"
WRITE(tkout,'(13X,A16,5X,A18)') "   ERPK or KL   ","Type of Projection"
WRITE(tkout,'(13X,G16.10,13X,A1)') p_p%KL,p_p%TYPE_PJ
IF ((p_p%TYPE_PJ == "M").OR.(p_p%TYPE_PJ == "W")) THEN
  WRITE(tkout,*) "  -Rayon of Earth (in meters) :"
ELSE
  WRITE(tkout,*) "  -Distance between Equator and Pole of projection on"
  WRITE(tkout,*) "   projection plane (in meters) :"
ENDIF
WRITE(tkout,'(13X,G16.10)') p_p%R_EQUATEUR
WRITE(tkout,*) "  -Pole of projection (-1.0 for South, 1.0 for North, 0.0 for"
WRITE(tkout,*) "   Mercator projection : no sense) :"
WRITE(tkout,'(13X,F4.1)') p_p%POLE
WRITE(tkout,*) "============================================================="
IF (lhook) CALL dr_hook('EGGPACK:INFO_PP_PRINT',1,zhook_handle)
END SUBROUTINE info_pp_print
! -------------------------------------------------------------------------------
LOGICAL FUNCTION return_print(YD_CODE_ERR,K_NUM_TEST,AUTO_STOP,KOUT)
CHARACTER(LEN=36), PARAMETER                 :: FMT = '(A25,I4," ; test number = ",I3,A100)'
type(error), INTENT(IN)                     :: yd_code_err
INTEGER(KIND=JPIM), INTENT(IN)               :: K_NUM_TEST
LOGICAL, INTENT(IN), OPTIONAL                :: AUTO_STOP
INTEGER(KIND=JPIM), INTENT(IN), OPTIONAL     :: KOUT

CHARACTER(LEN=25)  :: CL_ADD_TEXT
LOGICAL            :: TAS
INTEGER(KIND=JPIM) :: TKOUT
REAL(KIND=JPRD)    :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:RETURN_PRINT',0,zhook_handle)
IF (PRESENT(kout))THEN
  tkout = kout
ELSE
  tkout = 6_jpim
ENDIF
IF (PRESENT(auto_stop))THEN
  tas = auto_stop
ELSE
  tas = .true.
ENDIF
SELECT CASE (yd_code_err%NUM)
CASE(:-1_jpim) ; cl_add_text = 'ERROR   : return value = '
CASE(0_jpim)   ; cl_add_text = 'OK      : return value = '
CASE(1_jpim)   ; cl_add_text = 'INFO    : return value = '
CASE(2_jpim:)  ; cl_add_text = 'WARNING : return value = '
END SELECT
IF (yd_code_err%NUM < 0_jpim) WRITE (tkout,*) "Subroutine last status when aborted : "
WRITE (tkout,fmt) cl_add_text,yd_code_err%NUM,k_num_test,yd_code_err%TXT
IF (yd_code_err%NUM < 0_jpim) THEN
  IF (tas) THEN
    CALL abor1("Abort by EGGPACK:RETURN_PRINT") 
  ELSE
    return_print = .true.
  ENDIF
ELSE
  return_print = .false.
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:RETURN_PRINT',1,zhook_handle)
END FUNCTION return_print
! -------------------------------------------------------------------------------
FUNCTION type_proj(REF_COORD) RESULT (TY_PJ)
! REF_COORD in Degrees
type(lola), INTENT(IN)       :: ref_coord
CHARACTER(LEN=1) :: TY_PJ
 
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:TYPE_PROJ',0,zhook_handle)    
IF (ref_coord%LAT == 0.0_jprd) THEN
  ty_pj = "M"
ELSEIF (abs(ref_coord%LAT) == 90.0_jprd) THEN
  ty_pj = "S"
ELSE
  ty_pj = "L"
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:TYPE_PROJ',1,zhook_handle)
END FUNCTION type_proj
! -------------------------------------------------------------------------------
REAL(KIND=JPRD) FUNCTION pole_is (REF_COORD) RESULT (POL)
! REF_COORD in Degrees
type(lola), INTENT(IN)         :: ref_coord

REAL(KIND=JPRD)    :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:POLE_IS',0,zhook_handle)
IF (ref_coord%LAT == 0.0_jprd) THEN
  pol = 0.0_jprd
ELSE
  pol = sign(1.0_jprd,ref_coord%LAT)
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:POLE_IS',1,zhook_handle)
END FUNCTION pole_is
! -------------------------------------------------------------------------------

! ******************* Specifics functions ***************************************

! STPL : STEREOGRAPHIQUE POLAIRE / LAMBERT MODE FONCTIONS 
! Transform LAT-LON to R-TETA, R-TETA to XY in STD Origin and reverse (-+R)
! -------------------------------------------------------------------------------
type(rteta) FUNCTION stpl_latlon_to_rteta_s(PT_COORD,P_PJ,PI) RESULT (PT_RTETA)
! PT_COORD in Radians
type(lola), INTENT(IN)                             :: pt_coord
type(param_proj), INTENT(IN)                       :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL               :: PI

REAL(KIND=JPRD) :: TPI
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STPL_LATLON_TO_RTETA_S',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
pt_rteta%R = p_pj%R_EQUATEUR*((tan((tpi/4.0_jprd)-((p_pj%POLE*pt_coord%LAT)/2.0_jprd)))**(p_pj%KL))
pt_rteta%TETA = p_pj%KL*dist_2ref(pt_coord,p_pj%REF_PT,tpi)
IF (lhook) CALL dr_hook('EGGPACK:STPL_LATLON_TO_RTETA_S',1,zhook_handle)
END FUNCTION stpl_latlon_to_rteta_s
! -------------------------------------------------------------------------------
type(xy) FUNCTION stlp_rteta_to_xy_s(PT_RTETA,P_PJ) RESULT (PT_XY)
! PT_XY in STD Origin that is pole ( Y positive to North, X positive to East )
type(rteta), INTENT(IN)            :: pt_rteta
type(param_proj), INTENT(IN)       :: p_pj

REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_XY_S',0,zhook_handle)
pt_xy%X = pt_rteta%R*sin(pt_rteta%TETA)
pt_xy%Y = -p_pj%POLE*pt_rteta%R*cos(pt_rteta%TETA)
IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_XY_S',1,zhook_handle)
END FUNCTION stlp_rteta_to_xy_s
! -------------------------------------------------------------------------------
type(rteta) FUNCTION stlp_xy_to_rteta_s(PT_XY,P_PJ,PI) RESULT (PT_RTETA)
! PT_XY in STD Origin that is pole ( Y positive to North, X positive to East )
type(xy), INTENT(IN)                     :: pt_xy
type(param_proj), INTENT(IN)             :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL     :: PI

REAL(KIND=JPRD)     :: TPI,TATNG
REAL(KIND=JPRD)     :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STLP_XY_TO_RTETA_S',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
pt_rteta%R = sqrt((pt_xy%X*pt_xy%X)+(pt_xy%Y*pt_xy%Y))
IF (pt_xy%Y == 0.0_jprd) THEN
  IF (pt_xy%X == 0.0_jprd) THEN
    tatng = tpi
  ELSE
    tatng = sign(tpi/2.0_jprd,-p_pj%POLE*pt_xy%X)
  ENDIF
ELSE
  tatng = atan(-p_pj%POLE*(pt_xy%X/pt_xy%Y))
ENDIF
pt_rteta%TETA = tpi*sign(1.0_jprd,pt_xy%X)*(sign(0.5_jprd,p_pj%POLE*pt_xy%Y)+0.5_jprd)+tatng
! This term : <-------------------------------------------------------------> modifies the value of atan() according
! the square defined by : y<0 => 0 ; x>0 and y>0 => pi ; x<0 and y>0 => -pi
!========>if PT_RTETA%TETA > TPI*P_PJ%KL then error (bad section of planed cone)
IF (abs(pt_rteta%TETA) > tpi*p_pj%KL) THEN
  print *,"Point at x = ",pt_xy%X," , y = ",pt_xy%Y
  print *,"Is out of the planed cone section ! Abort !!!"
  CALL abor1("Abort by EGGPACK:STLP_XY_TO_RTETA_S")
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:STLP_XY_TO_RTETA_S',1,zhook_handle)
END FUNCTION stlp_xy_to_rteta_s
! -------------------------------------------------------------------------------
type(lola) FUNCTION stlp_rteta_to_latlon_s(PT_RTETA,P_PJ,PI) RESULT (PT_COORD)
! PT_COORD in Radians
type(rteta), INTENT(IN)                    :: pt_rteta
type(param_proj), INTENT(IN)               :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL       :: PI

REAL(KIND=JPRD) :: TPI
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_LATLON_S',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
pt_coord%LON = p_pj%REF_PT%LON + pt_rteta%TETA/p_pj%KL
pt_coord%LAT = p_pj%POLE*((tpi/2.0_jprd)-2.0_jprd*atan((pt_rteta%R/p_pj%R_EQUATEUR)**(1.0_jprd/p_pj%KL)))
IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_LATLON_S',1,zhook_handle)
END FUNCTION stlp_rteta_to_latlon_s
! -------------------------------------------------------------------------------
FUNCTION stpl_latlon_to_rteta_v(PT_COORD,P_PJ,PI) RESULT (PT_RTETA)
! PT_COORD in Radians
type(lola), DIMENSION(:), INTENT(IN)                :: pt_coord
type(param_proj), INTENT(IN)                        :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL                :: PI
type(rteta), DIMENSION(SIZE(PT_COORD)) :: pt_rteta

REAL(KIND=JPRD) :: TPI
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STPL_LATLON_TO_RTETA_V',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
pt_rteta(:)%R = p_pj%R_EQUATEUR*((tan((tpi/4.0_jprd)-((p_pj%POLE*pt_coord(:)%LAT)/2.0_jprd)))**(p_pj%KL))
pt_rteta(:)%TETA = p_pj%KL*dist_2ref(pt_coord(:),p_pj%REF_PT,tpi)
IF (lhook) CALL dr_hook('EGGPACK:STPL_LATLON_TO_RTETA_V',1,zhook_handle)
END FUNCTION stpl_latlon_to_rteta_v
! -------------------------------------------------------------------------------
FUNCTION stlp_rteta_to_xy_v(PT_RTETA,P_PJ) RESULT (PT_XY)
! PT_XY in STD Origin that is pole ( Y positive to North, X positive to East )
type(rteta), DIMENSION(:), INTENT(IN)            :: pt_rteta
type(param_proj), INTENT(IN)                     :: p_pj
type(xy), DIMENSION(SIZE(PT_RTETA)) :: pt_xy

REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_XY_V',0,zhook_handle)
pt_xy(:)%X = pt_rteta(:)%R*sin(pt_rteta(:)%TETA)
pt_xy(:)%Y = -p_pj%POLE*pt_rteta(:)%R*cos(pt_rteta(:)%TETA)
IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_XY_V',1,zhook_handle)
END FUNCTION stlp_rteta_to_xy_v
! -------------------------------------------------------------------------------
FUNCTION stlp_xy_to_rteta_v(PT_XY,P_PJ,PI) RESULT (PT_RTETA)
! PT_XY in STD Origin that is pole ( Y positive to North, X positive to East )
type(xy), DIMENSION(:), INTENT(IN)               :: pt_xy
type(param_proj), INTENT(IN)                     :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL             :: PI
type(rteta), DIMENSION(SIZE(PT_XY)) :: pt_rteta

REAL(KIND=JPRD)                         :: TPI
REAL(KIND=JPRD), DIMENSION(SIZE(PT_XY)) :: TATNG
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STLP_XY_TO_RTETA_V',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
pt_rteta(:)%R = sqrt((pt_xy(:)%X*pt_xy(:)%X)+(pt_xy(:)%Y*pt_xy(:)%Y))
WHERE (pt_xy(:)%Y == 0.0_jprd)
  WHERE (pt_xy(:)%X == 0.0_jprd)
    tatng = tpi
  ELSEWHERE
    tatng = sign(tpi/2.0_jprd,-p_pj%POLE*pt_xy(:)%X)
  ENDWHERE
ELSEWHERE
  tatng = atan(-p_pj%POLE*(pt_xy(:)%X/pt_xy(:)%Y))
ENDWHERE
pt_rteta(:)%TETA = tpi*sign(1.0_jprd,pt_xy(:)%X)*(sign(0.5_jprd,p_pj%POLE*pt_xy(:)%Y)+0.5_jprd)+tatng(:)
! This term : <-------------------------------------------------------------------> modifies the value of atan()
! according the square defined by : y<0 => 0 ; x>0 and y>0 => pi ; x<0 and y>0 => -pi
!========>if one of PT_RTETA(:)%TETA > TPI*P_PJ%KL then error (bad section of planed cone)
IF (any(abs(pt_rteta(:)%TETA) > tpi*p_pj%KL)) THEN
  print *,"Some points are out of planed cone section ! Abort !!!"
  CALL abor1("Abort by EGGPACK:STLP_XY_TO_RTETA_V")       
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:STLP_XY_TO_RTETA_V',1,zhook_handle)
END FUNCTION stlp_xy_to_rteta_v
! -------------------------------------------------------------------------------
FUNCTION stlp_rteta_to_latlon_v(PT_RTETA,P_PJ,PI) RESULT (PT_COORD)
! PT_COORD in Radians
type(rteta), DIMENSION(:), INTENT(IN)              :: pt_rteta
type(param_proj), INTENT(IN)                       :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL               :: PI
type(lola), DIMENSION(SIZE(PT_RTETA)) :: pt_coord

REAL(KIND=JPRD) :: TPI
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_LATLON_V',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
pt_coord(:)%LON = p_pj%REF_PT%LON + pt_rteta(:)%TETA/p_pj%KL
pt_coord(:)%LAT = p_pj%POLE*((tpi/2.0_jprd)-2.0_jprd*atan((pt_rteta(:)%R/p_pj%R_EQUATEUR)**(1.0_jprd/p_pj%KL)))
IF (lhook) CALL dr_hook('EGGPACK:STLP_RTETA_TO_LATLON_V',1,zhook_handle)
END FUNCTION stlp_rteta_to_latlon_v
! -------------------------------------------------------------------------------

! ******************* Generics functions ****************************************

! Creation of Structure Definition projection type (-+D)
! -------------------------------------------------------------------------------
type(param_proj) FUNCTION ref_datas (REF_COORD,RA,TOZERO_COORD,LRT) RESULT (P_P)
! REF_COORD,TOZERO_COORD in -+Degrees
type(lola), INTENT(IN)                    :: ref_coord
REAL(KIND=JPRD), INTENT(IN), OPTIONAL      :: RA
type(lola), INTENT(IN), OPTIONAL          :: tozero_coord
LOGICAL, INTENT(IN), OPTIONAL              :: LRT

REAL(KIND=JPRD) :: RT
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:REF_DATAS',0,zhook_handle)
IF (PRESENT(ra))THEN
  rt = ra
ELSE
  rt = r_earth
ENDIF
! P_P%REF in -+Radians
p_p%REF_PT = lolar(angle_domain(ref_coord)) ! expect REF_COORD in Dg & put it in -+
IF (PRESENT(lrt)) THEN ! flag rot-tilt present
  IF ((lrt).AND.(p_p%REF_PT%LAT==0.0_jprd)) THEN ! flag rot-tilt true and mercator
    p_p%TYPE_PJ = "W"
    IF (PRESENT(tozero_coord)) THEN ! value of TOZERO_COORD
      p_p%TZO_PT = lolar(angle_domain(tozero_coord)) ! expect TOZERO_COORD in Dg & put it in -+
    ELSE ! default TOZERO_COORD value is (0,0)
      p_p%TZO_PT%LAT = 0.0_jprd
      p_p%TZO_PT%LON = 0.0_jprd
    ENDIF
  ELSE ! either flag false either not mercator so
    p_p%TYPE_PJ = type_proj(ref_coord) ! doing standard
    p_p%TZO_PT%LAT = -999.999_jprd ! no need of TOZERO_COORD
    p_p%TZO_PT%LON = -999.999_jprd
  ENDIF
ELSE ! not rot-tilt mode so
  p_p%TYPE_PJ = type_proj(ref_coord) ! doing standard
  p_p%TZO_PT%LAT = -999.999_jprd ! no need of TOZERO_COORD
  p_p%TZO_PT%LON = -999.999_jprd
ENDIF
p_p%POLE = pole_is(ref_coord)
p_p%KL = p_p%POLE*sin(p_p%REF_PT%LAT)
IF (p_p%KL /= 0.0_jprd) THEN
  ! Rho of projection of equator on polar plane (Rho-Theta)
  p_p%R_EQUATEUR = rt*((cos(p_p%REF_PT%LAT))**(1.0_jprd -p_p%KL))*(((1.0_jprd +p_p%KL)**p_p%KL)/p_p%KL)
ELSE
  p_p%R_EQUATEUR = rt
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:REF_DATAS',1,zhook_handle)
END FUNCTION ref_datas
! -------------------------------------------------------------------------------
! Change of origine between STD and NEW and reverse (-+D)
! -------------------------------------------------------------------------------
type(xy) FUNCTION xy_new_to_std_origin_s(NEW_ORIGIN_COORD,PT_XY_IN_NEW_ORIGIN,P_PJ,PI) RESULT (PT_XY_IN_STD_ORIGIN)
! NEW_ORIGIN_COORD in -+Degrees
type(lola), INTENT(IN)                               :: new_origin_coord
type(xy), INTENT(IN)                                 :: pt_xy_in_new_origin
type(param_proj), INTENT(IN)                         :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL                 :: PI

REAL(KIND=JPRD) :: TPI
type(xy)       :: n_o_pt_xy
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:XY_NEW_TO_STD_ORIGIN_S',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
n_o_pt_xy = latlon_to_xy(lolar(new_origin_coord),p_pj,tpi)
pt_xy_in_std_origin%X = pt_xy_in_new_origin%X+n_o_pt_xy%X
pt_xy_in_std_origin%Y = pt_xy_in_new_origin%Y+n_o_pt_xy%Y
IF (lhook) CALL dr_hook('EGGPACK:XY_NEW_TO_STD_ORIGIN_S',1,zhook_handle)
END FUNCTION xy_new_to_std_origin_s
! -------------------------------------------------------------------------------
type(xy) FUNCTION xy_std_to_new_origin_s(NEW_ORIGIN_COORD,PT_XY_IN_STD_ORIGIN,P_PJ,PI) RESULT (PT_XY_IN_NEW_ORIGIN)
! NEW_ORIGIN_COORD in -+Degrees
type(lola), INTENT(IN)                             :: new_origin_coord
type(xy), INTENT(IN)                               :: pt_xy_in_std_origin
type(param_proj), INTENT(IN)                       :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL               :: PI

REAL(KIND=JPRD) :: TPI
type(xy)       :: n_o_pt_xy
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:XY_STD_TO_NEW_ORIGIN_S',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
n_o_pt_xy = latlon_to_xy(lolar(new_origin_coord),p_pj,tpi)
pt_xy_in_new_origin%X = pt_xy_in_std_origin%X-n_o_pt_xy%X
pt_xy_in_new_origin%Y = pt_xy_in_std_origin%Y-n_o_pt_xy%Y
IF (lhook) CALL dr_hook('EGGPACK:XY_STD_TO_NEW_ORIGIN_S',1,zhook_handle)
END FUNCTION xy_std_to_new_origin_s
! -------------------------------------------------------------------------------
FUNCTION xy_new_to_std_origin_v(NEW_ORIGIN_COORD,PT_XY_IN_NEW_ORIGIN,P_PJ,PI) RESULT (PT_XY_IN_STD_ORIGIN)
! NEW_ORIGIN_COORD in -+Degrees
type(lola), INTENT(IN)                                      :: new_origin_coord
type(xy), DIMENSION(:), INTENT(IN)                          :: pt_xy_in_new_origin
type(param_proj), INTENT(IN)                                :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL                        :: PI
type(xy), DIMENSION(SIZE(PT_XY_IN_NEW_ORIGIN)) :: pt_xy_in_std_origin

REAL(KIND=JPRD) :: TPI
type(xy)       :: n_o_pt_xy
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:XY_NEW_TO_STD_ORIGIN_V',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
n_o_pt_xy = latlon_to_xy(lolar(new_origin_coord),p_pj,tpi)
pt_xy_in_std_origin(:)%X = pt_xy_in_new_origin(:)%X+n_o_pt_xy%X
pt_xy_in_std_origin(:)%Y = pt_xy_in_new_origin(:)%Y+n_o_pt_xy%Y
IF (lhook) CALL dr_hook('EGGPACK:XY_NEW_TO_STD_ORIGIN_V',1,zhook_handle)
END FUNCTION xy_new_to_std_origin_v
! -------------------------------------------------------------------------------
FUNCTION xy_std_to_new_origin_v(YL_NEW_ORIGIN_COORD,YL_PT_XY_IN_STD_ORIGIN,P_PJ,PI) RESULT (YD_PT_XY_IN_NEW_ORIGIN)
! YL_NEW_ORIGIN_COORD in -+Degrees
type(lola), INTENT(IN)                                         :: yl_new_origin_coord
type(xy), DIMENSION(:), INTENT(IN)                             :: yl_pt_xy_in_std_origin
type(param_proj), INTENT(IN)                                   :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL                           :: PI
type(xy), DIMENSION(SIZE(YL_PT_XY_IN_STD_ORIGIN)) :: yd_pt_xy_in_new_origin

REAL(KIND=JPRD) :: TPI
type(xy)       :: yl_n_o_pt_xy
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:XY_STD_TO_NEW_ORIGIN_V',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
yl_n_o_pt_xy = latlon_to_xy(lolar(yl_new_origin_coord),p_pj,tpi)
yd_pt_xy_in_new_origin(:)%X = yl_pt_xy_in_std_origin(:)%X-yl_n_o_pt_xy%X
yd_pt_xy_in_new_origin(:)%Y = yl_pt_xy_in_std_origin(:)%Y-yl_n_o_pt_xy%Y
IF (lhook) CALL dr_hook('EGGPACK:XY_STD_TO_NEW_ORIGIN_V',1,zhook_handle)
END FUNCTION xy_std_to_new_origin_v
! -------------------------------------------------------------------------------
! Coordinates transforms between XY in STD Origin and LAT-LON (-+R)
! -------------------------------------------------------------------------------
type(xy) FUNCTION latlon_to_xy_s(PT_COORD,P_PJ,PI) RESULT (PT_XY)
! PT_COORD in -+Radians
! PT_XY in STD Origin
type(lola), INTENT(IN)                           :: pt_coord
type(param_proj), INTENT(IN)                     :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL             :: PI

REAL(KIND=JPRD) :: TPI
type(lola)     :: pt_coord2
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:LATLON_TO_XY_S',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
IF ((p_pj%TYPE_PJ == "S").OR.(p_pj%TYPE_PJ == "L")) THEN
  pt_xy = stlp_rteta_to_xy(stpl_latlon_to_rteta(pt_coord,p_pj,tpi),p_pj)
ELSE
  pt_coord2 = pt_coord
  IF (p_pj%TYPE_PJ == "W") THEN
    pt_coord2 = metilrot(p_pj%REF_PT,p_pj%TZO_PT,pt_coord)
    pt_xy%X = pt_coord2%LON*p_pj%R_EQUATEUR
  ELSE
    pt_xy%X = p_pj%R_EQUATEUR*dist_2ref(pt_coord,p_pj%REF_PT,tpi)
  ENDIF
  pt_xy%Y = -p_pj%R_EQUATEUR*log(tan((tpi/4.0_jprd)-(pt_coord2%LAT/2.0_jprd)))
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:LATLON_TO_XY_S',1,zhook_handle)
END FUNCTION latlon_to_xy_s
! -------------------------------------------------------------------------------
type(lola) FUNCTION xy_to_latlon_s(PT_XY,P_PJ,PI) RESULT (PT_COORD)
! PT_XY in STD Origin
! PT_COORD in -+Radians??
type(xy), INTENT(IN)                           :: pt_xy
type(param_proj), INTENT(IN)                   :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL           :: PI

REAL(KIND=JPRD) :: TPI
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:XY_TO_LATLON_S',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
IF ((p_pj%TYPE_PJ == "S").OR.(p_pj%TYPE_PJ == "L")) THEN
  pt_coord = stlp_rteta_to_latlon(stlp_xy_to_rteta(pt_xy,p_pj,tpi),p_pj,tpi)
ELSE
  pt_coord%LON = (pt_xy%X/p_pj%R_EQUATEUR)
  pt_coord%LAT = (tpi/2.0_jprd)-2.0_jprd*atan(exp(-(pt_xy%Y/p_pj%R_EQUATEUR)))
  IF (p_pj%TYPE_PJ == "W") THEN
    pt_coord=merotil(p_pj%REF_PT,p_pj%TZO_PT,pt_coord)
  ELSE
    pt_coord%LON = p_pj%REF_PT%LON+pt_coord%LON
  ENDIF
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:XY_TO_LATLON_S',1,zhook_handle)
END FUNCTION xy_to_latlon_s
! -------------------------------------------------------------------------------
FUNCTION latlon_to_xy_v(PT_COORD,P_PJ,PI) RESULT (PT_XY)
! PT_COORD in -+Radians
! PT_XY in STD Origin
type(lola), DIMENSION(:), INTENT(IN)             :: pt_coord
type(param_proj), INTENT(IN)                     :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL             :: PI
type(xy), DIMENSION(SIZE(PT_COORD)) :: pt_xy

REAL(KIND=JPRD)                        :: TPI
type(lola), DIMENSION(SIZE(PT_COORD)) :: pt_coord2
REAL(KIND=JPRD)                        :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:LATLON_TO_XY_V',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
IF ((p_pj%TYPE_PJ == "S").OR.(p_pj%TYPE_PJ == "L")) THEN
  pt_xy = stlp_rteta_to_xy(stpl_latlon_to_rteta(pt_coord,p_pj,tpi),p_pj)
ELSE
  pt_coord2(:) = pt_coord(:)
  IF (p_pj%TYPE_PJ == "W") THEN
    pt_coord2(:) = metilrot(p_pj%REF_PT,p_pj%TZO_PT,pt_coord(:))
    pt_xy(:)%X = pt_coord2(:)%LON*p_pj%R_EQUATEUR
  ELSE
    pt_xy(:)%X = p_pj%R_EQUATEUR*dist_2ref(pt_coord(:),p_pj%REF_PT,tpi)
  ENDIF
  pt_xy(:)%Y = -p_pj%R_EQUATEUR*log(tan((tpi/4.0_jprd)-(pt_coord2(:)%LAT/2.0_jprd)))
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:LATLON_TO_XY_V',1,zhook_handle)
END FUNCTION latlon_to_xy_v
! -------------------------------------------------------------------------------
FUNCTION xy_to_latlon_v(YL_PT_XY,P_PJ,PI) RESULT (PT_COORD)
! YL_PT_XY in STD Origin
! PT_COORD in -+Radians??
type(xy), DIMENSION(:), INTENT(IN)                 :: yl_pt_xy
type(param_proj), INTENT(IN)                       :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL               :: PI
type(lola), DIMENSION(SIZE(YL_PT_XY)) :: pt_coord

REAL(KIND=JPRD) :: TPI
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:XY_TO_LATLON_V',0,zhook_handle)
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
IF ((p_pj%TYPE_PJ == "S").OR.(p_pj%TYPE_PJ == "L")) THEN
  pt_coord = stlp_rteta_to_latlon(stlp_xy_to_rteta(yl_pt_xy,p_pj,tpi),p_pj,tpi)
ELSE
  pt_coord(:)%LON = (yl_pt_xy(:)%X/p_pj%R_EQUATEUR)
  pt_coord(:)%LAT = (tpi/2.0_jprd)-2.0_jprd*atan(exp(-(yl_pt_xy(:)%Y/p_pj%R_EQUATEUR)))
  IF (p_pj%TYPE_PJ == "W") THEN
    pt_coord(:)=merotil(p_pj%REF_PT,p_pj%TZO_PT,pt_coord(:))
  ELSE
    pt_coord(:)%LON = p_pj%REF_PT%LON+pt_coord(:)%LON
  ENDIF
ENDIF
IF (lhook) CALL dr_hook('EGGPACK:XY_TO_LATLON_V',1,zhook_handle)
END FUNCTION xy_to_latlon_v
! -------------------------------------------------------------------------------
! Functions MAP_FACTOR and GN
! -------------------------------------------------------------------------------
REAL(KIND=JPRD) FUNCTION map_factor_s(PT_COORD,P_PJ,PI,RA)
! PT_COORD in -+Radians??
type(lola), INTENT(IN)                    :: pt_coord
type(param_proj), INTENT(IN)              :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL      :: RA, PI

REAL(KIND=JPRD) :: RT, TPI
type(rteta)    :: pt_rteta
type(lola)     :: pt_coord2
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:MAP_FACTOR_S',0,zhook_handle)
IF (PRESENT(ra))THEN
  rt = ra
ELSE
  rt = r_earth
ENDIF
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
SELECT CASE(p_pj%TYPE_PJ)
CASE('W') ; pt_coord2 = metilrot(p_pj%REF_PT,p_pj%TZO_PT,pt_coord)
  map_factor_s = 1.0_jprd/(cos(pt_coord2%LAT))
CASE('M') ; map_factor_s = 1.0_jprd/(cos(pt_coord%LAT))
CASE('S') ; map_factor_s = 2.0_jprd/(1.0_jprd +p_pj%POLE*sin(pt_coord%LAT))
CASE('L') ; pt_rteta = stpl_latlon_to_rteta(pt_coord,p_pj,tpi)
  map_factor_s = (p_pj%KL*pt_rteta%R)/(rt*cos(pt_coord%LAT))
END SELECT
IF (lhook) CALL dr_hook('EGGPACK:MAP_FACTOR_S',1,zhook_handle)
END FUNCTION map_factor_s
! -------------------------------------------------------------------------------
FUNCTION map_factor_v(PT_COORD,P_PJ,PI,RA)
! PT_COORD in -+Radians??
type(lola), DIMENSION(:), INTENT(IN)                    :: pt_coord
type(param_proj), INTENT(IN)                            :: p_pj
REAL(KIND=JPRD), INTENT(IN), OPTIONAL                    :: RA, PI
REAL(KIND=JPRD), DIMENSION(SIZE(PT_COORD))  :: MAP_FACTOR_V

REAL(KIND=JPRD)                         :: RT, TPI
type(rteta), DIMENSION(SIZE(PT_COORD)) :: yl_pt_rteta
type(lola), DIMENSION(SIZE(PT_COORD))  :: pt_coord2
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:MAP_FACTOR_V',0,zhook_handle)
IF (PRESENT(ra))THEN
  rt = ra
ELSE
  rt = r_earth
ENDIF
IF (PRESENT(pi)) THEN
  tpi = pi
ELSE
  tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
SELECT CASE(p_pj%TYPE_PJ)
CASE('W') ; pt_coord2(:) = metilrot(p_pj%REF_PT,p_pj%TZO_PT,pt_coord(:))
  map_factor_v(:) = 1.0_jprd/(cos(pt_coord2(:)%LAT))
CASE('M') ; map_factor_v(:) = 1.0_jprd/(cos(pt_coord(:)%LAT))
CASE('S') ; map_factor_v(:) = 2.0_jprd/(1.0_jprd +p_pj%POLE*sin(pt_coord(:)%LAT))
CASE('L') ; yl_pt_rteta     = stpl_latlon_to_rteta(pt_coord,p_pj,tpi)
  map_factor_v(:) = (p_pj%KL*yl_pt_rteta(:)%R)/(rt*cos(pt_coord(:)%LAT))
END SELECT
IF (lhook) CALL dr_hook('EGGPACK:MAP_FACTOR_V',1,zhook_handle)
END FUNCTION map_factor_v
! -------------------------------------------------------------------------------
type(pgn) FUNCTION gn_s(PT_COORD,P_PJ)
! PT_COORD in -+Radians??
type(lola), INTENT(IN)                           :: pt_coord
type(param_proj), INTENT(IN)                     :: p_pj

type(lola)     :: pt_coord2
REAL(KIND=JPRD) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:GN_S',0,zhook_handle)
SELECT CASE(p_pj%TYPE_PJ)
CASE('W') ; pt_coord2 = metilrot(p_pj%REF_PT,p_pj%TZO_PT,pt_coord)
  gn_s%ONY = (((cos(p_pj%REF_PT%LON)*((cos(p_pj%TZO_PT%LAT)*cos(pt_coord%LAT))+ &
   & (sin(p_pj%TZO_PT%LAT)*sin(pt_coord%LAT)*cos(pt_coord%LON-p_pj%TZO_PT%LON)))) &
   & -(sin(p_pj%REF_PT%LON)*sin(pt_coord%LAT)*sin(pt_coord%LON-p_pj%TZO_PT%LON))) &
   & /cos(pt_coord2%LAT))
  gn_s%ONX = -((cos(p_pj%REF_PT%LON)*sin(p_pj%TZO_PT%LAT)*sin(pt_coord%LON-p_pj%TZO_PT%LON)) &
   & +(sin(p_pj%REF_PT%LON)*cos(pt_coord%LON-p_pj%TZO_PT%LON)))/cos(pt_coord2%LAT)
CASE('M') ; gn_s%ONX = 0.0_jprd
  gn_s%ONY = 1.0_jprd
CASE DEFAULT ; gn_s%ONX = -p_pj%POLE*sin(p_pj%KL*(pt_coord%LON-p_pj%REF_PT%LON))
  gn_s%ONY = cos(p_pj%KL*(pt_coord%LON-p_pj%REF_PT%LON))
END SELECT
IF (lhook) CALL dr_hook('EGGPACK:GN_S',1,zhook_handle)
END FUNCTION gn_s
! -------------------------------------------------------------------------------
FUNCTION gn_v(YD_PT_COORD,YD_P_PJ)
! YD_PT_COORD in -+Radians??
type(lola), DIMENSION(:), INTENT(IN)                      :: yd_pt_coord
type(param_proj), INTENT(IN)                              :: yd_p_pj
type(pgn), DIMENSION(SIZE(YD_PT_COORD))      :: gn_v

type(lola), DIMENSION(SIZE(YD_PT_COORD)) :: yd_pt_coord2
REAL(KIND=JPRD)                           :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('EGGPACK:GN_V',0,zhook_handle)
SELECT CASE(yd_p_pj%TYPE_PJ)
CASE('W') ; yd_pt_coord2(:) = metilrot(yd_p_pj%REF_PT,yd_p_pj%TZO_PT,yd_pt_coord(:))
  gn_v(:)%ONY = (((cos(yd_p_pj%REF_PT%LON)*((cos(yd_p_pj%TZO_PT%LAT)*cos(yd_pt_coord(:)%LAT))+ &
   & (sin(yd_p_pj%TZO_PT%LAT)*sin(yd_pt_coord(:)%LAT)*cos(yd_pt_coord(:)%LON-yd_p_pj%TZO_PT%LON)))) &
   & -(sin(yd_p_pj%REF_PT%LON)*sin(yd_pt_coord(:)%LAT)*sin(yd_pt_coord(:)%LON-yd_p_pj%TZO_PT%LON))) &
   & /cos(yd_pt_coord2(:)%LAT))
  gn_v(:)%ONX = -((cos(yd_p_pj%REF_PT%LON)*sin(yd_p_pj%TZO_PT%LAT)*sin(yd_pt_coord(:)%LON-yd_p_pj%TZO_PT%LON)) &
   & +(sin(yd_p_pj%REF_PT%LON)*cos(yd_pt_coord(:)%LON-yd_p_pj%TZO_PT%LON)))/cos(yd_pt_coord2(:)%LAT)
CASE('M') ; gn_v(:)%ONX = 0.0_jprd
  gn_v(:)%ONY = 1.0_jprd
CASE DEFAULT ; gn_v(:)%ONX = -yd_p_pj%POLE*sin(yd_p_pj%KL*(yd_pt_coord(:)%LON-yd_p_pj%REF_PT%LON))
  gn_v(:)%ONY = cos(yd_p_pj%KL*(yd_pt_coord(:)%LON-yd_p_pj%REF_PT%LON))
END SELECT
IF (lhook) CALL dr_hook('EGGPACK:GN_V',1,zhook_handle)
END FUNCTION gn_v
! -------------------------------------------------------------------------------

! =================== SUBROUTINE ================================================

! ******************* subroutine to make grid domain ****************************

! -------------------------------------------------------------------------------
SUBROUTINE makdo(YD_REF_COORD,YD_CENTER_COORD,YD_PDEL,YD_NB_PTS,YD_GRID_COORD,P_GRID_MF, &
 & YD_GRID_PGN,YD_GRID_INFO,YD_ERR_CODE,LD_LIP,LD_AUTO_STOP,PI,P_RA,KOUT,LD_LMRT)
! Input Coordinates are in -+Degrees, output in 0+Radians excepted in YD_GRID_INFO : 0+D sauf Ref & Cen -+D
type(lola), INTENT(INOUT)                                            :: yd_ref_coord, yd_center_coord
type(delta), INTENT(IN)                                              :: yd_pdel
type(nbpts), INTENT(IN)                                              :: yd_nb_pts
REAL(KIND=JPRD), INTENT(IN), OPTIONAL                                 :: PI, P_RA
LOGICAL, INTENT(IN), OPTIONAL                                         :: LD_AUTO_STOP, LD_LIP, LD_LMRT
type(lola), DIMENSION(YD_NB_PTS%ONX,YD_NB_PTS%ONY), INTENT(OUT)      :: yd_grid_coord
REAL(KIND=JPRD), DIMENSION(YD_NB_PTS%ONX,YD_NB_PTS%ONY), INTENT(OUT)  :: P_GRID_MF
type(pgn), DIMENSION(YD_NB_PTS%ONX,YD_NB_PTS%ONY), INTENT(OUT)       :: yd_grid_pgn
type(domi), INTENT(OUT)                                              :: yd_grid_info
type(error), INTENT(OUT)                                             :: yd_err_code
INTEGER(KIND=JPIM), INTENT(IN), OPTIONAL                              :: KOUT
type(error), DIMENSION(-6:4)                     :: yl_tab_err_def = (/                                            &
 & error(-6_jpim, " : In Mercator, Deformations are too big out of [-85.0,85.0] (Rotated/Tilted mode)!"),           &
 & error(-5_jpim, " : In Lambert, the pole must be out of the domain !"),                                           &
 & error(-4_jpim, " : In Mercator, Deformations are too big out of [-85.0,85.0] !"),                                &
 & error(-3_jpim, " : Center point degrees coordinates out of bounds !"),                                           &
 & error(-2_jpim, " : Reference point degrees coordinates out of bounds !"),                                        &
 & error(-1_jpim, " : Subroutine aborted, sorry ..."),                                                              & 
 & error( 0_jpim, " : Subroutine finished successly."),                                                             &
 & error( 1_jpim, " : Test OK, subroutine go ahead !"),                                                             &
 & error( 2_jpim, " : In Mercator, It's better to use Lambert or St.Pol.  if ABS(CENTER_LAT) > 20.0 !"),            &
 & error( 3_jpim, " : In St.Pol. , It's better to use Lambert or Mercator if ABS(CENTER_LAT) < 70.0 !"),            &
 & error( 4_jpim, " : In Lambert , It's better to use St.Pol. or Mercator if ABS(REF_LAT) is out of [20.0,70.0] !") &
 & /)
REAL(KIND=JPRD)                                   :: Z_RT, Z_TPI
LOGICAL                                           :: LL_TAS, LL_TLIP, LL_TLMRT
INTEGER(KIND=JPIM)                                :: I, J, I_TKOUT
type(param_proj)                                 :: yl_p_p
type(xy), DIMENSION(YD_NB_PTS%ONX,YD_NB_PTS%ONY) :: yl_grid_xy_c, yl_grid_xy_p
REAL(KIND=JPRD)                                   :: ZHOOK_HANDLE
IF (lhook) CALL dr_hook('EGGPACK:MAKDO',0,zhook_handle)
IF (PRESENT(ld_lip))THEN
  ll_tlip = ld_lip
ELSE
  ll_tlip = _defp_
ENDIF
IF (PRESENT(ld_lmrt))THEN
  ll_tlmrt = ld_lmrt
ELSE
  ll_tlmrt = .false.
ENDIF
IF (PRESENT(kout))THEN
  i_tkout = kout
ELSE
  i_tkout = 6_jpim
ENDIF
IF (PRESENT(ld_auto_stop))THEN
  ll_tas = ld_auto_stop
ELSE
  ll_tas = .true.
ENDIF
IF (PRESENT(p_ra))THEN
  z_rt = p_ra
ELSE
  z_rt = r_earth
ENDIF
IF (PRESENT(pi)) THEN
  z_tpi = pi
ELSE
  z_tpi = asin(1.0_jprd)*2.0_jprd
ENDIF
WRITE(i_tkout,*) "Begining of Makin Domain MAKDO subroutine"
! Validation of inputs LAT_LON coordinates or init YD_ERR_CODE to no error
yd_err_code = yl_tab_err_def(val_coord(yd_ref_coord,-2_jpim,z_tpi))
test_1: IF (return_print(yd_err_code,1_jpim,ll_tas,i_tkout)) THEN
  IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
  RETURN
ENDIF test_1
yd_err_code = yl_tab_err_def(val_coord(yd_center_coord,-3_jpim,z_tpi))
test_2: IF (return_print(yd_err_code,2_jpim,ll_tas,i_tkout)) THEN
  IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
  RETURN
ENDIF test_2
! Compute parameters of projection and init of domain structures
yd_grid_info%CT_COORD  = yd_center_coord
yd_grid_info%RF_COORD  = yd_ref_coord
yl_p_p                 = ref_datas(yd_ref_coord,z_rt,yd_center_coord,ll_tlmrt)
yd_grid_info%INFO_PROJ = yl_p_p
! Degrees to Radians
yd_center_coord        = lolar(yd_center_coord)
yd_ref_coord           = yl_p_p%REF_PT
! Compute XY grid points under CENTER origin
DO i=1_jpim,yd_nb_pts%ONX
  DO j=1_jpim,yd_nb_pts%ONY
    yl_grid_xy_c(i,j)%X = (REAL(i,kind=jprd)-(REAL(yd_nb_pts%onx+1_jpim,kind=jprd)/2.0_jprd))*yd_pdel%ONX
    yl_grid_xy_c(i,j)%Y = (REAL(j,kind=jprd)-(REAL(yd_nb_pts%ony+1_jpim,kind=jprd)/2.0_jprd))*yd_pdel%ONY
  ENDDO
ENDDO
IF (yl_p_p%TYPE_PJ/='W') THEN
  ! Change XY coordinates in CENTER Origin in STD Origin
  DO j=1_jpim,yd_nb_pts%ONY
    yl_grid_xy_p(1:yd_nb_pts%ONX,j) = &
     & xy_new_to_std_origin(yd_grid_info%CT_COORD,yl_grid_xy_c(1:yd_nb_pts%ONX,j),yl_p_p,z_tpi)
  ENDDO
ELSE
  yl_grid_xy_p = yl_grid_xy_c
ENDIF
! Validation depending projection type
SELECT CASE (yl_p_p%TYPE_PJ)
CASE('M')
  ! Validation of Mercator projection
  IF (abs(yd_grid_info%CT_COORD%LAT) > 20.0_jprd) yd_err_code = yl_tab_err_def(2_jpim)
  test_3: IF (return_print(yd_err_code,3_jpim,ll_tas,i_tkout)) THEN
    IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
    RETURN
  ENDIF test_3
  IF (maxval(abs(yl_grid_xy_p(:,:)%Y)) > 3.0_jprd*z_rt) yd_err_code = yl_tab_err_def(-4_jpim)
  test_4: IF (return_print(yd_err_code,4_jpim,ll_tas,i_tkout)) THEN
    IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
    RETURN
  ENDIF test_4
CASE('S')
  ! Validation of Polar Stereographic projection
  IF (abs(yd_grid_info%CT_COORD%LAT) < 70.0_jprd) yd_err_code = yl_tab_err_def(3_jpim)
  test_5: IF (return_print(yd_err_code,5_jpim,ll_tas,i_tkout)) THEN
    IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
    RETURN
  ENDIF test_5
CASE('L')
  ! Validation of Lambert projection
  IF ((abs(yd_grid_info%RF_COORD%LAT) < 20.0_jprd).OR.(abs(yd_grid_info%RF_COORD%LAT) > 70.0_jprd)) &
   & yd_err_code = yl_tab_err_def(4_jpim)
  test_6: IF (return_print(yd_err_code,6_jpim,ll_tas,i_tkout)) THEN
    IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
    RETURN
  ENDIF test_6
  IF (((yl_grid_xy_p(1,1)%X*yl_grid_xy_p(yd_nb_pts%ONX,1)%X) < 0.0_jprd).AND. &
   & ((yl_grid_xy_p(1,1)%Y*yl_grid_xy_p(1,yd_nb_pts%ONY)%Y) < 0.0_jprd)) yd_err_code = yl_tab_err_def(-5_jpim)
  test_7: IF (return_print(yd_err_code,7_jpim,ll_tas,i_tkout)) THEN
    IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
    RETURN
  ENDIF test_7
CASE('W')
  ! Validation of Mercator projection in rotated/tilted case
  IF (maxval(abs(yl_grid_xy_p(:,:)%Y)) > 3.0_jprd*z_rt) yd_err_code = yl_tab_err_def(-6_jpim)
  test_8: IF (return_print(yd_err_code,8_jpim,ll_tas,i_tkout)) THEN
    IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
    RETURN
  ENDIF test_8
END SELECT
! Compute ouputs datas depending projection type
DO j=1_jpim,yd_nb_pts%ONY
  yd_grid_coord(1:yd_nb_pts%ONX,j) = xy_to_latlon(yl_grid_xy_p(1:yd_nb_pts%ONX,j),yl_p_p,z_tpi)
  p_grid_mf(1:yd_nb_pts%ONX,j)     = map_factor(yd_grid_coord(1:yd_nb_pts%ONX,j),yl_p_p,z_tpi,z_rt)
  yd_grid_pgn(1:yd_nb_pts%ONX,j)   = gn(yd_grid_coord(1:yd_nb_pts%ONX,j),yl_p_p)
  yd_grid_coord(1:yd_nb_pts%ONX,j) = angle_domain(yd_grid_coord(1:yd_nb_pts%ONX,j),z_tpi,'0+','R')
ENDDO
! Fill info structure
yd_grid_info%G_SIZE   = yd_nb_pts
yd_grid_info%SW_COORD = lolad(yd_grid_coord(1,1))
yd_grid_info%MF_SW    = p_grid_mf(1,1)
yd_grid_info%SE_COORD = lolad(yd_grid_coord(yd_nb_pts%ONX,1))
yd_grid_info%MF_SE    = p_grid_mf(yd_nb_pts%ONX,1)
yd_grid_info%NE_COORD = lolad(yd_grid_coord(yd_nb_pts%ONX,yd_nb_pts%ONY))
yd_grid_info%MF_NE    = p_grid_mf(yd_nb_pts%ONX,yd_nb_pts%ONY)
yd_grid_info%NW_COORD = lolad(yd_grid_coord(1,yd_nb_pts%ONY))
yd_grid_info%MF_NW    = p_grid_mf(1,yd_nb_pts%ONY)
yd_grid_info%MF_RF    = map_factor(yd_ref_coord,yl_p_p,z_tpi,z_rt)
yd_grid_info%MF_CT    = map_factor(yd_center_coord,yl_p_p,z_tpi,z_rt)
yd_ref_coord          = angle_domain(yd_ref_coord,z_tpi,'0+','R')
yd_center_coord       = angle_domain(yd_center_coord,z_tpi,'0+','R')
IF (ll_tlip) CALL info_print(yd_grid_info,i_tkout,z_tpi)
WRITE (i_tkout,*) "Subroutine last status when finished : "
IF (yd_err_code%NUM == 1_jpim) yd_err_code = yl_tab_err_def(0_jpim)
test_9: IF (return_print(yd_err_code,9_jpim,ll_tas,i_tkout)) THEN
  IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle) 
  RETURN
ENDIF test_9
IF (lhook) CALL dr_hook('EGGPACK:MAKDO',1,zhook_handle)
END SUBROUTINE makdo

#undef _DEFP_
! -------------------------------------------------------------------------------
END MODULE eggpack