!------------------------------------------------------------------------------
!        IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE         : Mohid Model
! PROJECT       : Mohid Base 2
! MODULE        : HorizontalGrid
! URL           : http://www.mohid.com
! AFFILIATION   : IST/MARETEC, Marine Modelling Group
! DATE          : May 2003
! REVISION      : Frank Braunschweig - v4.0
! DESCRIPTION   : Module to read and calculate a horizontal Grid
!
!------------------------------------------------------------------------------
Module ModuleHorizontalGrid

    use ModuleGlobalData
    use ModuleFunctions, only   : Rodaxy, FromCartesianToGrid, FromGridToCartesian
    use ModuleEnterData
    use ModuleDrawing
    use ModuleHDF5

    implicit none
    private

    !Subroutine----------------------------------------------------------------

    !Constructor
    public  :: ConstructHorizontalGrid
    private ::      AllocateInstance
    private ::      ConstructGlobalVariables
    private ::          AllocateVariables
    private ::          Mercator
    private ::              ConvertCoordenates     !Fromer CONVCO_BAT
    private ::          GridPointsLocation1D
    private ::          GridPointsLocation2D
    private ::          ComputeDistances
    private ::          ComputeRotation
    private ::          ComputeGridCellArea
    private ::          ComputeCoriolis
    public  :: ConstructFatherGridLocation
    private ::      ConstructNewFatherGrid1D
    private ::      ConstructNewFatherGrid2D
    private ::      Add_FatherGrid


    !Modifier
    public  :: WriteHorizontalGrid
    public  :: LocateCell

    !Selector
    public  :: GetHorizontalGridSize
    public  :: GetHorizontalGrid
    public  :: GetGridOrigin
    public  :: GetGridAngle
    public  :: GetGridZone
    public  :: GetLatitudeLongitude
    public  :: GetGridCoordType
    public  :: GetCoordTypeList
    public  :: GetGridMeanLatLong
    public  :: GetCoriolisFrequency
    public  :: GetGridLatitudeLongitude
    public  :: GetComputeZUV
    public  :: GetGridFileName
    public  :: GetCheckDistortion
    public  :: GetGridRotation
    public  :: GetDefineCellsMap
    public  :: GetFatherGridID
    public  :: GetGridCellArea
    private ::      Search_FatherGrid

    public  :: UnGetHorizontalGrid


    !Destructor
    public  :: KillHorizontalGrid
    private ::      DeallocateInstance
    private ::      KillFatherGridList


    !Management
    private ::      Ready

    !Auxilary
    private :: USCONVCO
    private :: UTGP83
    private :: GPUT83
    private :: IUTPC
    private :: IUTGP
    private :: DRGPUT
    private :: DRUTGP
    private :: TODMS
    private :: DATUMM
    private :: TMGRID
    private :: TCONST
    private :: TCONPC
    private :: TMGEOD

    private :: wgs84_to_rd
    private :: wgs842bessel_
    private :: bessel2rd_
    private :: rd_to_wgs84
    private :: rd2bessel_
    private :: bessel2wgs84_

    !Interfaces----------------------------------------------------------------

    private :: UnGetHorizontalGrid1d
    private :: UnGetHorizontalGrid2d
    private :: UngetHorizontalGrid2DInt
    interface  UnGetHorizontalGrid
        module procedure UnGetHorizontalGrid1D
        module procedure UnGetHorizontalGrid2D
        module procedure UnGetHorizontalGrid2DInt
    end interface  UnGetHorizontalGrid


    public :: InterpolRegularGrid
    interface  InterpolRegularGrid
        module procedure InterpolRegularGrid2D
        module procedure InterpolRegularGrid3D
        module procedure InterpolRegularGrid3D8
    end interface  InterpolRegularGrid

    private:: RotateVectorFieldToGrid2D
    private:: RotateVectorFieldToGrid3D
    public :: RotateVectorFieldToGrid
    interface  RotateVectorFieldToGrid
        module procedure RotateVectorFieldToGrid2D
        module procedure RotateVectorFieldToGrid3D
    end interface  RotateVectorFieldToGrid

    private:: RotateVectorGridToField2D
    private:: RotateVectorGridToField3D
    public :: RotateVectorGridToField
    interface  RotateVectorGridToField
        module procedure RotateVectorGridToField2D
        module procedure RotateVectorGridToField3D
    end interface  RotateVectorGridToField



    !Parameter-----------------------------------------------------------------

    !Types of coordinates
    integer, parameter :: GEOG_        = 1         !Coordenadas Geograficas  
    integer, parameter :: UTM_         = 2         !Coordenadas (UTM)
    integer, parameter :: MIL_PORT_    = 3         !Coordenadas Militares Portuguesas     
    integer, parameter :: SIMPLE_GEOG_ = 4         !Coordenadas Geograficas Simplificadas
                                                   !SIMPLE_GEOG_ - considera a terra como uma esfera
    integer, parameter :: GRID_COORD_  = 5         !Coordenadas relativas à origem da malha
    integer, parameter :: CIRCULAR_    = 6         !Coordenadas circulares (XX- raio, YY - anglo em graus)
    integer, parameter :: NLRD_        = 7         !Coordenadas Netherlands RD

    integer, parameter :: PORTUGUESE_UTM_ZONE_ = 29

    !Block information
    character(LEN = StringLength), parameter :: BeginXX         = '<BeginXX>'
    character(LEN = StringLength), parameter :: EndXX           = '<EndXX>'
    character(LEN = StringLength), parameter :: BeginYY         = '<BeginYY>'
    character(LEN = StringLength), parameter :: EndYY           = '<EndYY>'
    character(LEN = StringLength), parameter :: BeginCornersXY  = '<CornersXY>'
    character(LEN = StringLength), parameter :: EndCornersXY    = '<\CornersXY>'


    !Calculation Points
    integer, parameter :: ComputeZ_     = 1
    integer, parameter :: ComputeU_     = 2
    integer, parameter :: ComputeV_     = 3
    integer, parameter :: ComputeCross_ = 4
    integer, parameter :: ComputeZU_    = 5
    integer, parameter :: ComputeZV_    = 6


    integer, dimension(1:4), parameter :: diFace = (/0,1 ,0,-1/), djFace = (/-1, 0, 1,  0/)

    integer, parameter :: LargeScaleModel_ = 1
    integer, parameter :: Assimila_        = 2 
    integer, parameter :: SurfaceMM5_      = 3
    
    !Type----------------------------------------------------------------------
    type T_Compute
        real,    dimension(:),   pointer :: XX_Z, YY_Z, XX_U, YY_U, XX_V, YY_V
        real,    dimension(:),   pointer :: XX_Cross, YY_Cross
        real,    dimension(:,:), pointer :: XX2D_Z, YY2D_Z, XX2D_U, YY2D_U, XX2D_V, YY2D_V
    end type T_Compute

    type T_FatherGrid
        integer                          :: GridID
        real,    dimension(:,:), pointer :: XX_Z, YY_Z, XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross
        logical                          :: OkZ, OkU, OkV, OkCross
        logical                          :: CornersXYInput = .false.
        real                             :: Fhc 
        integer                          :: JX, IY
        integer, dimension(:,:), pointer :: IZ, JZ, IU, JU, IV, JV, ICross, JCross
        type (T_Size2D)                  :: MPI_Window
        type (T_FatherGrid),     pointer :: Next
        type (T_FatherGrid),     pointer :: Prev
    end type T_FatherGrid

    type T_HorizontalGrid
        integer                                 :: InstanceID

        !Former Bathymetry    
        real, pointer, dimension(:  )           :: XX, YY
        real                                    :: Xorig      = null_real
        real                                    :: Yorig      = null_real
        real                                    :: Latitude   = null_real
        real                                    :: Longitude  = null_real
        real                                    :: GRID_ANGLE = null_real
        integer                                 :: CoordType  = null_int
        integer                                 :: Zone       = null_int

        type(T_Compute)                         :: Compute

        type (T_FatherGrid),     pointer        :: FirstFatherGrid
        type (T_FatherGrid),     pointer        :: LastFatherGrid

        !Distances (DXX... DVY)
        real, dimension(:, :), pointer          :: DXX, DYY, DZX, DZY
        real, dimension(:, :), pointer          :: DUX, DUY, DVX, DVY
        real, dimension(:, :), pointer          :: XX_IE, YY_IE
        real, dimension(:, :), pointer          :: LatitudeConn, LongitudeConn

        real, dimension(:, :), pointer          :: RotationX, RotationY
        logical                                 :: CornersXYInput = .false., Distortion, RegularRotation

        integer, dimension(:,:), pointer        :: DefineCellsMap, DefineFacesUMap, &
                                                   DefineFacesVMap, DefineCrossMap

        !Latitude, Longitude
        real, dimension(:, :), pointer          :: LatitudeZ, LongitudeZ
        
        !Coriolis Factor
        real, dimension(:, :), pointer          :: F

        !Area           
        real, dimension(:, :), pointer          :: GridCellArea

        !Other
        type (T_Size2D)                         :: Size, WorkSize
        character(PathLength)                   :: FileName

        type (T_HorizontalGrid), pointer        :: Next
        
    end type T_HorizontalGrid

    !Global Module Variables
    type (T_HorizontalGrid), pointer            :: FirstHorizontalGrid
    type (T_HorizontalGrid), pointer            :: Me


    contains


    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    !CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR 

    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    subroutine ConstructHorizontalGrid(HorizontalGridID, DataFile, STAT)

        !Arguments-------------------------------------------------------------
        integer                             :: HorizontalGridID
        character(len=*), optional          :: DataFile
        integer, optional,  intent(OUT)     :: STAT    

        !Local-----------------------------------------------------------------
        integer                             :: STAT_, ready_
        integer                             :: STAT_CALL

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

        STAT_ = UNKNOWN_

        !Assures nullification of the global variable
        if (.not. ModuleIsRegistered(mHorizontalGrid_)) then
            nullify (FirstHorizontalGrid)
            call RegisterModule (mHorizontalGrid_) 
        endif

        call Ready(HorizontalGridID, ready_)    

cd2 :   if (ready_ .EQ. OFF_ERR_) then

            call AllocateInstance 
            
            nullify (Me%FirstFatherGrid)
            nullify (Me%LastFatherGrid )

            !Reads data file
            if (.not. present(DataFile)) then
                call ReadFileName('IN_BATIM', Me%FileName, Message = "Grid File", STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'ConstructHorizontalGrid - HorizontalGrid - ERR01'
            else
                Me%FileName = DataFile
            endif

            !Construct the variable common to all module
            call ConstructGlobalVariables


           !Constructs XX_IE, YY_IE, LatitudeZ and LongitudeZ
            call Mercator

            if (.not. Me%CornersXYInput) then

                !Constructs XX_Z, YY_Z, XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross
                call GridPointsLocation1D

            else
                !Constructs XX2D_Z, YY2D_Z, XX2D_U, YY2D_U, XX2D_V, YY2D_V, XX2D_Cross, YY2D_Cross
                call GridPointsLocation2D

            endif

            !Computes DXX, DYY, DZX, DZY, DUX, DUY, DVX, DVY
            call ComputeDistances

            !Computes RotationX, RotationY
            call ComputeRotation

            !Computes Area
            call ComputeGridCellArea

            !Computes Coriolis
            call ComputeCoriolis


            !Returns ID
            HorizontalGridID    = Me%InstanceID

            STAT_ = SUCCESS_
        else 
            
            stop 'HorizontalGrid - ConstructHorizontalGrid - ERR99' 

        end if cd2


        if (present(STAT)) STAT = STAT_

    end subroutine ConstructHorizontalGrid
 
    !--------------------------------------------------------------------------

    subroutine AllocateInstance 

        !Arguments-------------------------------------------------------------
    
        !Local-----------------------------------------------------------------
        type (T_HorizontalGrid), pointer            :: NewObjHorizontalGrid
        type (T_HorizontalGrid), pointer            :: PreviousObjHorizontalGrid


        !Allocates new instance
        allocate (NewObjHorizontalGrid)
        nullify  (NewObjHorizontalGrid%Next)

        !Insert New Instance into list and makes Current point to it
        if (.not. associated(FirstHorizontalGrid)) then
            FirstHorizontalGrid     => NewObjHorizontalGrid
            Me                      => NewObjHorizontalGrid
        else
            PreviousObjHorizontalGrid   => FirstHorizontalGrid
            Me                          => FirstHorizontalGrid%Next
            do while (associated(Me))
                PreviousObjHorizontalGrid  => Me
                Me                         => Me%Next
            enddo
            Me                             => NewObjHorizontalGrid
            PreviousObjHorizontalGrid%Next => NewObjHorizontalGrid
        endif

        Me%InstanceID = RegisterNewInstance (mHORIZONTALGRID_)

    end subroutine AllocateInstance

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

    subroutine ConstructFatherGridLocation(HorizontalGridID, HorizontalGridFatherID, & 
                                           GridID, OkCross, OkZ, OkU, OkV, Window,   &
                                           STAT)

        !Arguments-------------------------------------------------------------
        integer                                   :: HorizontalGridID
        integer                                   :: HorizontalGridFatherID
        integer, optional,           intent (IN)  :: GridID
        logical, optional,           intent (IN)  :: OkCross, OkZ, OkU, OkV
        type (T_Size2D), optional                 :: Window
        integer, optional,           intent (OUT) :: STAT    

        !Local-----------------------------------------------------------------
        type (T_HorizontalGrid), pointer          :: ObjHorizontalGridFather
        type (T_FatherGrid), pointer              :: NewFatherGrid
        integer                                   :: STAT_, ready_, GridID_
        logical                                   :: OkZ_, OkU_, OkV_, OkCross_

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

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

            if (present(GridID)) then 
                GridID_ = GridID
            else
                GridID_ = LargeScaleModel_
            endif


            if (present(OkCross)) then 
                OkCross_ = OkCross
            else
                OkCross_ = .false.
            endif

            if (present(OkZ)) then 
                OkZ_ = OkZ
            else
                OkZ_ = .true.
            endif

            if (present(OkU)) then 
                OkU_ = OkU
            else
                OkU_ = .true.
            endif


            if (present(OkV)) then 
                OkV_ = OkV
            else
                OkV_ = .true.
            endif

            nullify(ObjHorizontalGridFather)
            call LocateObjFather        (ObjHorizontalGridFather, HorizontalGridFatherID)

            if (.not. ObjHorizontalGridFather%CornersXYInput) then

                call ConstructNewFatherGrid1D (ObjHorizontalGridFather, NewFatherGrid, GridID_, &
                                             OkZ_, OkU_, OkV_, OkCross_)
            else
                call ConstructNewFatherGrid2D (ObjHorizontalGridFather, NewFatherGrid, GridID_, &
                                             OkZ_, OkU_, OkV_, OkCross_)
            endif

            call Add_FatherGrid         (NewFatherGrid)

            if (present(Window)) Window = NewFatherGrid%MPI_Window

            STAT_ = SUCCESS_

        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_


    end subroutine ConstructFatherGridLocation
 
    !--------------------------------------------------------------------------

    subroutine ConstructNewFatherGrid1D(ObjHorizontalGridFather, NewFatherGrid, GridID,    &
                                      OkZ, OkU, OkV, OkCross)

        !Arguments-------------------------------------------------------------
        type (T_HorizontalGrid), pointer          :: ObjHorizontalGridFather
        type (T_FatherGrid)    , pointer          :: NewFatherGrid
        integer,                intent (IN)       :: GridID
        logical,                intent (IN)       :: OkZ, OkU, OkV, OkCross

        !Local-----------------------------------------------------------------
        real,    dimension(:,:), pointer          :: XX_Z, YY_Z, XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross  
        integer, dimension(:,:), pointer          :: IZ,  JZ,  IU,  JU,  IV,  JV, ICross, JCross
        real                                      :: Xorig   , Yorig   , Rotation 
        integer                                   :: ILB, IUB, JLB, JUB, i, j
        integer                                   :: ILBwork, IUBwork, JLBwork, JUBwork

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


        if (Me%CornersXYInput) stop 'ConstructNewFatherGrid1D - ModuleHoriuzontalGrid - ERR01'

        ILB = Me%Size%ILB
        IUB = Me%Size%IUB

        JLB = Me%Size%JLB
        JUB = Me%Size%JUB

        ILBwork = Me%WorkSize%ILB
        IUBwork = Me%WorkSize%IUB

        JLBwork = Me%WorkSize%JLB
        JUBwork = Me%WorkSize%JUB

        allocate(NewFatherGrid)

        NewFatherGrid%GridID            = GridID
        NewFatherGrid%OkZ               = OkZ
        NewFatherGrid%OkU               = OkU
        NewFatherGrid%OkV               = OkV
        NewFatherGrid%OkCross           = OkCross
        NewFatherGrid%MPI_Window%ILB    = -null_int
        NewFatherGrid%MPI_Window%IUB    =  null_int
        NewFatherGrid%MPI_Window%JLB    = -null_int
        NewFatherGrid%MPI_Window%JUB    =  null_int

        NewFatherGrid%CornersXYInput    = .false.

        nullify (NewFatherGrid%XX_Z )
        nullify (NewFatherGrid%YY_Z )
        nullify (NewFatherGrid%XX_U )
        nullify (NewFatherGrid%YY_U )
        nullify (NewFatherGrid%XX_V )
        nullify (NewFatherGrid%YY_V )

        nullify (NewFatherGrid%Next )
        nullify (NewFatherGrid%Prev )


        !Compute points location ib the father grid
        if (NewFatherGrid%OkZ) then

            !Allocates Variables
            allocate (NewFatherGrid%XX_Z (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%YY_Z (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%IZ   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JZ   (ILB:IUB, JLB:JUB))

            !Initialize Values
            NewFatherGrid%XX_Z (:,:)  = FillValueReal
            NewFatherGrid%YY_Z (:,:)  = FillValueReal

            NewFatherGrid%IZ   (:,:)  = FillValueInt
            NewFatherGrid%JZ   (:,:)  = FillValueInt

            XX_Z  => NewFatherGrid%XX_Z
            YY_Z  => NewFatherGrid%YY_Z
            IZ    => NewFatherGrid%IZ
            JZ    => NewFatherGrid%JZ

        endif

        if (NewFatherGrid%OkU) then

            allocate (NewFatherGrid%XX_U (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%YY_U (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%IU   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JU   (ILB:IUB, JLB:JUB))

            NewFatherGrid%XX_U (:,:)  = FillValueReal
            NewFatherGrid%YY_U (:,:)  = FillValueReal
            NewFatherGrid%IU   (:,:)  = FillValueInt
            NewFatherGrid%JU   (:,:)  = FillValueInt

            XX_U  => NewFatherGrid%XX_U
            YY_U  => NewFatherGrid%YY_U
            IU    => NewFatherGrid%IU
            JU    => NewFatherGrid%JU

        endif

        if (NewFatherGrid%OkV) then

            allocate (NewFatherGrid%XX_V (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%YY_V (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%IV   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JV   (ILB:IUB, JLB:JUB))

            NewFatherGrid%XX_V (:,:)  = FillValueReal
            NewFatherGrid%YY_V (:,:)  = FillValueReal
            NewFatherGrid%IV   (:,:)  = FillValueInt
            NewFatherGrid%JV   (:,:)  = FillValueInt

            YY_V  => NewFatherGrid%YY_V
            XX_V  => NewFatherGrid%XX_V
            IV    => NewFatherGrid%IV
            JV    => NewFatherGrid%JV

        endif


        if (NewFatherGrid%OkCross) then

            allocate (NewFatherGrid%XX_Cross (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%YY_Cross (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%ICross   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JCross   (ILB:IUB, JLB:JUB))

            NewFatherGrid%XX_Cross (:,:)  = FillValueReal
            NewFatherGrid%YY_Cross (:,:)  = FillValueReal
            NewFatherGrid%ICross   (:,:)  = FillValueInt
            NewFatherGrid%JCross   (:,:)  = FillValueInt

            YY_Cross  => NewFatherGrid%YY_Cross
            XX_Cross  => NewFatherGrid%XX_Cross
            ICross    => NewFatherGrid%ICross
            JCross    => NewFatherGrid%JCross

        endif


        Xorig    = Me%Xorig      - ObjHorizontalGridFather%Xorig
        Yorig    = Me%Yorig      - ObjHorizontalGridFather%Yorig

        Rotation = Me%Grid_Angle - ObjHorizontalGridFather%Grid_Angle

        if (Rotation /= 0.) then
            write(*,*) 'The submodels do not work with referentials with different rotations'
            stop  'ConstructNewFatherGrid1D - HorizontalGrid - ERR22a'
        endif

        call RODAXY(0., 0., -Me%Grid_Angle, Xorig, Yorig)  


        if (Xorig < 0) stop 'ConstructNewFatherGrid1D - HorizontalGrid - ERR23'
        if (Yorig < 0) stop 'ConstructNewFatherGrid1D - HorizontalGrid - ERR24'

        !Compute points location ib the father grid
        if (NewFatherGrid%OkZ) then
            
            !Rotation of the points
do3 :       do j = JLBwork, JUBwork
do4 :       do i = ILBwork, IUBwork

                XX_Z(i, j) = Me%Compute%XX_Z(j)

                YY_Z(i, j) = Me%Compute%YY_Z(i)

                call RODAXY(Xorig, Yorig, Rotation, XX_Z(i, j), YY_Z(i, j))  

                call LocateCell (ObjHorizontalGridFather%Compute%XX_Z,               &
                                 ObjHorizontalGridFather%Compute%YY_Z,               &
                                 XX_Z(i, j), YY_Z(i, j),                             &
                                 ObjHorizontalGridFather%WorkSize%ILB,               &
                                 ObjHorizontalGridFather%WorkSize%IUB,               &
                                 ObjHorizontalGridFather%WorkSize%JLB,               &
                                 ObjHorizontalGridFather%WorkSize%JUB,               &
                                 IZ(i, j), JZ(i, j))
                
                !Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IZ(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IZ(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JZ(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JZ(i, j))

            end do do4
            end do do3

            nullify(XX_Z , YY_Z )
            nullify(IZ   , JZ   )

        endif 

        if (NewFatherGrid%OkU) then

            !Rotation of the points
do1 :       do j = JLBwork, JUBwork+1
do2 :       do i = ILBwork, IUBwork

                XX_U(i, j) = Me%Compute%XX_U(j)

                YY_U(i, j) = Me%Compute%YY_U(i)

                call RODAXY(Xorig, Yorig, Rotation, XX_U(i, j), YY_U(i, j))  

                call LocateCell (ObjHorizontalGridFather%Compute%XX_U,               &
                                 ObjHorizontalGridFather%Compute%YY_U,               &
                                 XX_U(i, j), YY_U(i, j),                             &
                                 ObjHorizontalGridFather%WorkSize%ILB,               &
                                 ObjHorizontalGridFather%WorkSize%IUB,               &
                                 ObjHorizontalGridFather%WorkSize%JLB,               &
                                 ObjHorizontalGridFather%WorkSize%JUB + 1,           &
                                 IU(i, j), JU(i, j))

                !Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IU(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IU(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JU(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JU(i, j))

            end do do2
            end do do1

            nullify(XX_U , YY_U ) 
            nullify(IU   , JU   )

        endif

        if (NewFatherGrid%OkV) then

            !Rotation of the points
do5:       do j = JLBwork, JUBwork
do6:       do i = ILBwork, IUBwork+1

                XX_V(i, j) = Me%Compute%XX_V(j)

                YY_V(i, j) = Me%Compute%YY_V(i)
            
                call RODAXY(Xorig, Yorig, Rotation, XX_V(i, j), YY_V(i, j))  

                call LocateCell (ObjHorizontalGridFather%Compute%XX_V,               &
                                 ObjHorizontalGridFather%Compute%YY_V,               &
                                 XX_V(i, j), YY_V(i, j),                             &
                                 ObjHorizontalGridFather%WorkSize%ILB,               &
                                 ObjHorizontalGridFather%WorkSize%IUB + 1,           &
                                 ObjHorizontalGridFather%WorkSize%JLB,               &
                                 ObjHorizontalGridFather%WorkSize%JUB,               &
                                 IV(i, j), JV(i, j))

                !MPI_Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IV(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IV(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JV(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JV(i, j))

            end do do6
            end do do5

            nullify(XX_V , YY_V )
            nullify(IV   , JV   )

        endif


        if (NewFatherGrid%OkCross) then

            !Rotation of the points
do7:       do j = JLBwork, JUBwork
do8:       do i = ILBwork, IUBwork

                XX_Cross(i, j) = Me%Compute%XX_Z(j)

                YY_Cross(i, j) = Me%Compute%YY_Z(i)
            
                call RODAXY(Xorig, Yorig, Rotation, XX_Cross(i, j), YY_Cross(i, j))  

                call LocateCell (ObjHorizontalGridFather%Compute%XX_Cross,           &
                                 ObjHorizontalGridFather%Compute%YY_Cross,           &
                                 XX_Cross(i, j), YY_Cross(i, j),                     &
                                 ObjHorizontalGridFather%WorkSize%ILB + 1,           &
                                 ObjHorizontalGridFather%WorkSize%IUB + 1,           &
                                 ObjHorizontalGridFather%WorkSize%JLB,               &
                                 ObjHorizontalGridFather%WorkSize%JUB,               &
                                 ICross(i, j), JCross(i, j))

                !MPI_Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, ICross(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, ICross(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JCross(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JCross(i, j))

            end do do8
            end do do7

            nullify(XX_Cross, YY_Cross)
            nullify(ICross  , JCross  )

        endif

        !Increase MPI Window Size by one (Fluxes at the upper boundary)
        NewFatherGrid%MPI_Window%IUB = NewFatherGrid%MPI_Window%IUB + 1
        NewFatherGrid%MPI_Window%JUB = NewFatherGrid%MPI_Window%JUB + 1
        
    end subroutine ConstructNewFatherGrid1D


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

    subroutine ConstructNewFatherGrid2D(ObjHorizontalGridFather, NewFatherGrid, GridID,    &
                                      OkZ, OkU, OkV, OkCross)

        !Arguments-------------------------------------------------------------
        type (T_HorizontalGrid), pointer          :: ObjHorizontalGridFather
        type (T_FatherGrid)    , pointer          :: NewFatherGrid
        integer,                intent (IN)       :: GridID
        logical,                intent (IN)       :: OkZ, OkU, OkV, OkCross

        !Local-----------------------------------------------------------------
        real                                      :: XPoint, YPoint  
        integer, dimension(:,:), pointer          :: IZ,  JZ,  IU,  JU,  IV,  JV,       &
                                                     ICross, JCross, DefinedPoint
        integer                                   :: ILB, IUB, JLB, JUB, i, j
        integer                                   :: ILBwork, IUBwork, JLBwork, JUBwork

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

        if (.not. Me%CornersXYInput) stop 'ConstructNewFatherGrid2D - ModuleHoriuzontalGrid - ERR01'

        ILB = Me%Size%ILB
        IUB = Me%Size%IUB

        JLB = Me%Size%JLB
        JUB = Me%Size%JUB

        ILBwork = Me%WorkSize%ILB
        IUBwork = Me%WorkSize%IUB

        JLBwork = Me%WorkSize%JLB
        JUBwork = Me%WorkSize%JUB

        allocate(NewFatherGrid)

        NewFatherGrid%GridID            = GridID
        NewFatherGrid%OkZ               = OkZ
        NewFatherGrid%OkU               = OkU
        NewFatherGrid%OkV               = OkV
        NewFatherGrid%OkCross           = OkCross
        NewFatherGrid%MPI_Window%ILB    = -null_int
        NewFatherGrid%MPI_Window%IUB    =  null_int
        NewFatherGrid%MPI_Window%JLB    = -null_int
        NewFatherGrid%MPI_Window%JUB    =  null_int

        NewFatherGrid%CornersXYInput    = .true.

        nullify (NewFatherGrid%Next )
        nullify (NewFatherGrid%Prev )

        write(*,*) 'The submodels in curvilinear grids must have identical distortion'
        write(*,*) 'of the coarser grid'
        write(*,*) 'ConstructNewFatherGrid2D - HorizontalGrid - WARN01'


        !Compute points location ib the father grid
        if (NewFatherGrid%OkZ) then

            !Allocates Variables
            allocate (NewFatherGrid%IZ   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JZ   (ILB:IUB, JLB:JUB))

            !Initialize Values

            NewFatherGrid%IZ   (:,:)  = FillValueInt
            NewFatherGrid%JZ   (:,:)  = FillValueInt

            IZ    => NewFatherGrid%IZ
            JZ    => NewFatherGrid%JZ

            NewFatherGrid%XX_Z => Me%Compute%XX2D_Z
            NewFatherGrid%YY_Z => Me%Compute%YY2D_Z

            DefinedPoint    => ObjHorizontalGridFather%DefineCellsMap

        endif

        if (NewFatherGrid%OkU) then

            allocate (NewFatherGrid%IU   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JU   (ILB:IUB, JLB:JUB))

            NewFatherGrid%IU   (:,:)  = FillValueInt
            NewFatherGrid%JU   (:,:)  = FillValueInt

            IU    => NewFatherGrid%IU
            JU    => NewFatherGrid%JU

            NewFatherGrid%XX_U => Me%Compute%XX2D_U
            NewFatherGrid%YY_U => Me%Compute%YY2D_U


            DefinedPoint    => ObjHorizontalGridFather%DefineFacesUMap
        
        endif

        if (NewFatherGrid%OkV) then

            allocate (NewFatherGrid%IV   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JV   (ILB:IUB, JLB:JUB))

            NewFatherGrid%IV   (:,:)  = FillValueInt
            NewFatherGrid%JV   (:,:)  = FillValueInt

            IV    => NewFatherGrid%IV
            JV    => NewFatherGrid%JV

            NewFatherGrid%XX_V => Me%Compute%XX2D_V
            NewFatherGrid%YY_V => Me%Compute%YY2D_V

            DefinedPoint    => ObjHorizontalGridFather%DefineFacesVMap

        endif


        if (NewFatherGrid%OkCross) then

            allocate (NewFatherGrid%ICross   (ILB:IUB, JLB:JUB))
            allocate (NewFatherGrid%JCross   (ILB:IUB, JLB:JUB))

            NewFatherGrid%ICross   (:,:)  = FillValueInt
            NewFatherGrid%JCross   (:,:)  = FillValueInt

            ICross    => NewFatherGrid%ICross
            JCross    => NewFatherGrid%JCross

            NewFatherGrid%XX_Cross => Me%XX_IE
            NewFatherGrid%YY_Cross => Me%YY_IE

            DefinedPoint    => ObjHorizontalGridFather%DefineCrossMap         

        endif


        !Compute points location ib the father grid
        if (NewFatherGrid%OkZ) then
            
            !Rotation of the points
do3 :       do j = JLBwork, JUBwork
do4 :       do i = ILBwork, IUBwork

                if (Me%CornersXYInput) then

                    XPoint = Me%Compute%XX2D_Z(i, j)
                    YPoint = Me%Compute%YY2D_Z(i, j)

                else

                    XPoint = Me%Compute%XX_Z(j)
                    YPoint = Me%Compute%YY_Z(i)

                endif


                call LocateCellPolygons(ObjHorizontalGridFather%Compute%XX2D_Z,         &
                                        ObjHorizontalGridFather%Compute%YY2D_Z,         &
                                        XPoint, YPoint, DefinedPoint,                   &
                                        ObjHorizontalGridFather%WorkSize%ILB,           &
                                        ObjHorizontalGridFather%WorkSize%IUB,           &
                                        ObjHorizontalGridFather%WorkSize%JLB,           &
                                        ObjHorizontalGridFather%WorkSize%JUB,           &
                                        IZ(i, j), JZ(i, j))

               
                !Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IZ(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IZ(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JZ(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JZ(i, j))

            end do do4
            end do do3

            nullify(IZ   , JZ   )

        endif 

        if (NewFatherGrid%OkU) then

            !Rotation of the points
do1 :       do j = JLBwork, JUBwork+1
do2 :       do i = ILBwork, IUBwork

                call LocateCellPolygons(ObjHorizontalGridFather%Compute%XX2D_U,         &
                                        ObjHorizontalGridFather%Compute%YY2D_U,         &
                                        Me%Compute%XX2D_U(i, j),                        &
                                        Me%Compute%YY2D_U(i, j),                        &
                                        DefinedPoint,                                   &
                                        ObjHorizontalGridFather%WorkSize%ILB,           &
                                        ObjHorizontalGridFather%WorkSize%IUB,           &
                                        ObjHorizontalGridFather%WorkSize%JLB,           &
                                        ObjHorizontalGridFather%WorkSize%JUB + 1,       &
                                        IU(i, j), JU(i, j))

                !Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IU(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IU(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JU(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JU(i, j))

            end do do2
            end do do1

            nullify(IU   , JU   )

        endif

        if (NewFatherGrid%OkV) then

            !Rotation of the points
do5:       do j = JLBwork, JUBwork
do6:       do i = ILBwork, IUBwork+1

                call LocateCellPolygons(ObjHorizontalGridFather%Compute%XX2D_V,         &
                                        ObjHorizontalGridFather%Compute%YY2D_V,         &
                                        Me%Compute%XX2D_V(i, j),                        &
                                        Me%Compute%YY2D_V(i, j),                        &
                                        DefinedPoint,                                   &
                                        ObjHorizontalGridFather%WorkSize%ILB,           &
                                        ObjHorizontalGridFather%WorkSize%IUB + 1,       &
                                        ObjHorizontalGridFather%WorkSize%JLB,           &
                                        ObjHorizontalGridFather%WorkSize%JUB,           &
                                        IV(i, j), JV(i, j))

                !MPI_Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IV(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IV(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JV(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JV(i, j))

            end do do6
            end do do5

            nullify(IV   , JV   )

        endif


        if (NewFatherGrid%OkCross) then

            !Rotation of the points
do7:       do j = JLBwork, JUBwork
do8:       do i = ILBwork, IUBwork

                call LocateCellPolygons(ObjHorizontalGridFather%XX_IE,                  &
                                        ObjHorizontalGridFather%YY_IE,                  &
                                        Me%XX_IE(i, j),                                 &
                                        Me%YY_IE(i, j),                                 &
                                        DefinedPoint,                                   &
                                        ObjHorizontalGridFather%WorkSize%ILB + 1,       &
                                        ObjHorizontalGridFather%WorkSize%IUB + 1,       &
                                        ObjHorizontalGridFather%WorkSize%JLB,           &
                                        ObjHorizontalGridFather%WorkSize%JUB,           &
                                        ICross(i, j), JCross(i, j))

                !MPI_Window
                NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, ICross(i, j))
                NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, ICross(i, j))
                NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JCross(i, j))
                NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JCross(i, j))

            end do do8
            end do do7

            nullify(ICross  , JCross  )

        endif

        !Increase MPI Window Size by one (Fluxes at the upper boundary)
        NewFatherGrid%MPI_Window%IUB = NewFatherGrid%MPI_Window%IUB + 1
        NewFatherGrid%MPI_Window%JUB = NewFatherGrid%MPI_Window%JUB + 1
        
    end subroutine ConstructNewFatherGrid2D
    !--------------------------------------------------------------------------
    ! This subroutine adds a new grid to the father grid List  

    subroutine Add_FatherGrid(NewFatherGrid)

        !Arguments-------------------------------------------------------------
        type(T_FatherGrid),         pointer     :: NewFatherGrid

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

        ! Add to the WaterFatherGrid List a new father grid
        if (.not.associated(Me%FirstFatherGrid)) then

            Me%FirstFatherGrid    => NewFatherGrid
            Me%LastFatherGrid     => NewFatherGrid
        else
            NewFatherGrid%Prev     => Me%LastFatherGrid
            Me%LastFatherGrid%Next => NewFatherGrid
            Me%LastFatherGrid      => NewFatherGrid

        end if 

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

    end subroutine Add_FatherGrid 

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

    subroutine Search_FatherGrid(NewFatherGrid, GridID)

        !Arguments-------------------------------------------------------------
        type(T_FatherGrid),         pointer     :: NewFatherGrid
        integer,                    intent (IN) :: GridID


        !Local-----------------------------------------------------------------

        NewFatherGrid => Me%FirstFatherGrid

do1 :   do while (associated(NewFatherGrid)) 

cd1 :       if (NewFatherGrid%GridID == GridID) then
                exit do1
            else
                NewFatherGrid => NewFatherGrid%Next                 
            end if cd1

        end do do1

        if (.not.Associated(NewFatherGrid))                                              &
            call SetError(FATAL_, INTERNAL_, 'Search_FatherGrid - HorizontalGrid - ERR01')

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

    end subroutine Search_FatherGrid

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

    subroutine ConstructGlobalVariables

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        integer                             :: STAT_CALL
        integer                             :: ObjEnterData = 0 
        integer, dimension(2)               :: AuxInt
        real,    dimension(2)               :: AuxReal
        real                                :: DX, DY
        integer                             :: flag
        integer                             :: ClientNumber
        logical                             :: BlockFound, ConstantSpacingX, ConstantSpacingY
        integer                             :: FirstLine, LastLine, line, i, j

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

        !Nullify T_Distances
        nullify (Me%XX )
        nullify (Me%YY )
        nullify (Me%DXX)
        nullify (Me%DYY)
        nullify (Me%DZX)
        nullify (Me%DZY)
        nullify (Me%DUX)
        nullify (Me%DUY)
        nullify (Me%DVX)
        nullify (Me%DVY)
        nullify (Me%XX_IE)
        nullify (Me%YY_IE)
        nullify (Me%Compute%XX_Z )
        nullify (Me%Compute%YY_Z )
        nullify (Me%Compute%XX_U )
        nullify (Me%Compute%YY_U )
        nullify (Me%Compute%XX_V )
        nullify (Me%Compute%YY_V )
        nullify (Me%Compute%XX_Cross)
        nullify (Me%Compute%YY_Cross)


        !Nullify Other
        nullify (Me%LatitudeZ    )
        nullify (Me%LongitudeZ   )
        nullify (Me%F            )
        nullify (Me%GridCellArea )
        nullify (Me%LatitudeConn )
        nullify (Me%LongitudeConn)

        

        !Opens Data File
        call ConstructEnterData(ObjEnterData, Me%FileName, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR01'

        !Reads ILB_IUB
        call GetData(AuxInt,ObjEnterData, flag,                                 &
                     keyword      = 'ILB_IUB',                                  &
                     ClientModule = 'HorizontalGrid',                     &
                     STAT         = STAT_CALL)           
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR02'
        if (flag      /= 2       ) stop 'ConstructGlobalVariables - HorizontalGrid - ERR03'

        Me%Size%ILB     = AuxInt(1)-1 
        Me%Size%IUB     = AuxInt(2)+1 

        Me%WorkSize%ILB = AuxInt(1)
        Me%WorkSize%IUB = AuxInt(2)


        !Reads JLB_JUB
        call GetData(AuxInt,ObjEnterData, flag,                                         &
                     keyword      = 'JLB_JUB',                                          &
                     ClientModule = 'HorizontalGrid',                             &
                     STAT         = STAT_CALL)           
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR04'
        if (flag      /= 2       ) stop 'ConstructGlobalVariables - HorizontalGrid - ERR05'

        Me%Size%JLB     = AuxInt(1)-1 
        Me%Size%JUB     = AuxInt(2)+1 

        Me%WorkSize%JLB = AuxInt(1)
        Me%WorkSize%JUB = AuxInt(2)


        !Reads COORD_TIP
        call GetData(Me%CoordType, ObjEnterData, flag,                                  &
                     keyword      = 'COORD_TIP',                                        &
                     ClientModule = 'HorizontalGrid',                                   &
                     STAT         = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR06'

        if (.NOT. ((Me%CoordType == GEOG_       ) .OR.                                  &
                   (Me%CoordType == UTM_        ) .OR.                                  &
                   (Me%CoordType == MIL_PORT_   ) .OR.                                  &
                   (Me%CoordType == SIMPLE_GEOG_) .OR.                                  &
                   (Me%CoordType == GRID_COORD_ ) .OR.                                  &
                   (Me%CoordType == CIRCULAR_   ) .OR.                                  &
                   (Me%CoordType == NLRD_       ))) then   
            call SetError (FATAL_, INTERNAL_, "Invalid Coordinates")
        endif    
 
        !Reads ZONE
        if (Me%CoordType == UTM_ ) then
            call GetData(Me%Zone, ObjEnterData, flag,                                   &
                         keyword      = 'ZONE',                                         &
                         ClientModule = 'HorizontalGrid',                               &
                         STAT         = STAT_CALL)           
            if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR07'
        else if (Me%CoordType .EQ. MIL_PORT_) then 
            Me%ZONE = PORTUGUESE_UTM_ZONE_
        else 
            Me%ZONE = null_int
        endif
        
        
        !Reads ORIGIN
        call GetData(AuxReal,ObjEnterData, flag,                                        &
                     keyword      = 'ORIGIN',                                           &
                     ClientModule = 'HorizontalGrid',                                   &
                     STAT         = STAT_CALL)           
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR08'
        if (flag      /= 2       ) stop 'ConstructGlobalVariables - HorizontalGrid - ERR09'

        Me%Xorig = AuxReal(1)
        Me%Yorig = AuxReal(2)

        
        !Reads Latitude
        call GetData(Me%Latitude, ObjEnterData, flag,                                   &
                     keyword      = 'LATITUDE',                                         &
                     ClientModule = 'HorizontalGrid',                                   &
                     STAT         = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR10'

        !Reads Longitude
        call GetData(Me%Longitude, ObjEnterData, flag,                                  &
                     keyword      = 'LONGITUDE',                                        &
                     ClientModule = 'HorizontalGrid',                                   &
                     STAT         = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR11'


        !Reads GridAngle
        call GetData(Me%Grid_Angle, ObjEnterData, flag,                                 &
                     keyword      = 'GRID_ANGLE',                                       &
                     default      = 0.0,                                                &
                     ClientModule = 'HorizontalGrid',                                   &
                     STAT         = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR12'


        !Allocates XX and YY
        allocate(Me%XX(Me%Size%JLB+1 : Me%Size%JUB+1))
        allocate(Me%YY(Me%Size%ILB+1 : Me%Size%IUB+1))

        allocate(Me%XX_IE(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
        allocate(Me%YY_IE(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
        
        allocate(Me%LatitudeConn(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
        allocate(Me%LongitudeConn(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))

        allocate(Me%DefineCellsMap(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))

        Me%XX    = null_real
        Me%YY    = null_real
        Me%XX_IE = null_real
        Me%YY_IE = null_real

        Me%LatitudeConn  = null_real
        Me%LongitudeConn = null_real
        

        Me%DefineCellsMap(:,:)  = 1

        !Reads XX_YE, YY_IE
        call RewindBuffer(ObjEnterData, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR13'

        call ExtractBlockFromBuffer(ObjEnterData, ClientNumber,                 &
                                    BeginCornersXY, EndCornersXY, BlockFound,   &
                                    FirstLine = FirstLine, LastLine = LastLine, &
                                    STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR14'

BF:     if (BlockFound) then                                                  

            Me%CornersXYInput = .true. 

            if (Me%CoordType == CIRCULAR_ .and. Me%CoordType == GEOG_) then   
                stop 'ConstructGlobalVariables - HorizontalGrid - ERR14a'
            endif    


            line = FirstLine

            do i = Me%WorkSize%JLB, Me%WorkSize%IUB+1 
            do j = Me%WorkSize%JLB, Me%WorkSize%JUB+1 

                line = line+1

                !Last line found before end?
                if (line == LastLine) then
                    write(*,*) 
                    write(*,*) 'Error reading Corners X Y'
                    stop 'ConstructGlobalVariables - HorizontalGrid - ERR15'
                end if

                !Reads XX_IE , YY_IE 
                call GetData(AuxReal, ObjEnterData, flag,                                &
                             Buffer_Line  = Line,                                        &
                             STAT         = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR16'

                Me%XX_IE(I, J) = AuxReal(1)
                Me%YY_IE(I, J) = AuxReal(2)

            end do
            end do


            if (Me%CoordType == CIRCULAR_ .or. Me%CornersXYInput) then
            
                Me%Distortion      = .true. 
                Me%RegularRotation = .false. 

            else

                Me%Distortion = .false. 

                if (abs(Me%Grid_Angle)>0.) then
                    Me%RegularRotation = .true.
                else
                    Me%RegularRotation = .false.
                endif

            endif


            allocate(Me%DefineFacesUMap(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
            allocate(Me%DefineFacesVMap(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
            allocate(Me%DefineCrossMap (Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))


            Me%DefineFacesUMap(:,:) = 1
            Me%DefineFacesVMap(:,:) = 1
            Me%DefineCrossMap (:,:) = 1

            do i = Me%WorkSize%JLB, Me%WorkSize%IUB
            do j = Me%WorkSize%JLB, Me%WorkSize%JUB + 1


                if (Me%XX_IE(i  , j  ) <= FillValueReal / 2.  .or.      &
                    Me%XX_IE(i+1, j  ) <= FillValueReal / 2. ) then
                    Me%DefineFacesUMap(i,j)  = 0
                endif
 

            end do
            end do

            do i = Me%WorkSize%JLB, Me%WorkSize%IUB + 1
            do j = Me%WorkSize%JLB, Me%WorkSize%JUB


                if (Me%XX_IE(i  , j+1) <= FillValueReal / 2.  .or.      &
                    Me%XX_IE(i  , j  ) <= FillValueReal / 2. ) then
                    Me%DefineFacesVMap(i,j)  = 0
                endif
 

            end do
            end do

            do i = Me%WorkSize%JLB, Me%WorkSize%IUB + 1
            do j = Me%WorkSize%JLB, Me%WorkSize%JUB + 1


                if (Me%XX_IE(i  , j  ) <= FillValueReal / 2. ) then
                    Me%DefineCrossMap(i, j)  = 0
                endif
 

            end do
            end do


            !If one corner has FillValue Real, the point is not to consider
            do i = Me%WorkSize%JLB, Me%WorkSize%IUB
            do j = Me%WorkSize%JLB, Me%WorkSize%JUB

                if (Me%XX_IE(i  , j  ) <= FillValueReal / 2. .or.       &
                    Me%XX_IE(i  , j+1) <= FillValueReal / 2. .or.       &
                    Me%XX_IE(i+1, j+1) <= FillValueReal / 2. .or.       &
                    Me%XX_IE(i+1, j  ) <= FillValueReal / 2.) then
                    Me%DefineCellsMap(i, j)  = 0
                endif
 

            end do
            end do



            call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL) 
            if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR17'

        else BF

           !Check if the spacing in Y is constant
            call GetData(ConstantSpacingY,                                                   &
                         ObjEnterData, flag,                                                 &
                         SearchType   = FromFile,                                            &
                         keyword      ='CONSTANT_SPACING_Y',                                 &
                         STAT         = STAT_CALL)        
            if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR47'

iconY:      if (ConstantSpacingY) Then
                !Get grid origin
                call GetData(DY,                                                             &
                             ObjEnterData, flag,                                             &
                             SearchType   = FromFile,                                        &
                             keyword      ='DY',                                             &
                             STAT         = STAT_CALL)        
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR57'
                
                Me%YY(1) = 0.
            
                do i = Me%WorkSize%ILB + 1, Me%WorkSize%IUB + 1 
                    Me%YY(i) = Me%YY(i - 1) + DY
                end do

            else iconY

                !Reads YY
                call RewindBuffer(ObjEnterData, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR23'

                call ExtractBlockFromBuffer(ObjEnterData, ClientNumber,                 &
                                            BeginYY, EndYY, BlockFound,                 &
                                            FirstLine = FirstLine, LastLine = LastLine, &
                                            STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR24'

                if (BlockFound) then                                                  

                    line = FirstLine

                    do i = Me%WorkSize%ILB, Me%WorkSize%IUB+1 
                        line = line+1

                        !Last line found before end?
                        if (line == LastLine) then
                            write(*,*) 
                            write(*,*) 'Error reading YY'
                            stop 'ConstructGlobalVariables - HorizontalGrid - ERR25'
                        end if

                        call GetData(Me%YY(i), ObjEnterData, flag,                           &
                                     Buffer_Line  = Line,                                    &
                                     STAT         = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR26'
                    end do
                end if

                call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR27'

            end if iconY

           !Check if the spacing in X is constant
            call GetData(ConstantSpacingX,                                                   &
                         ObjEnterData, flag,                                                 &
                         SearchType   = FromFile,                                            &
                         keyword      ='CONSTANT_SPACING_X',                                 &
                         STAT         = STAT_CALL)        
            if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR67'

iconX:      if (ConstantSpacingX) Then
                !Get grid spacing dx
                call GetData(DX,                                                             &
                             ObjEnterData, flag,                                             &
                             SearchType   = FromFile,                                        &
                             keyword      ='DX',                                             &
                             STAT         = STAT_CALL)        
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR77'
                
                Me%XX(1) = 0.
            
                do j = Me%WorkSize%JLB + 1, Me%WorkSize%JUB + 1 
                    Me%XX(j) = Me%XX(j - 1) + DX
                end do

            else iconX


                !Reads XX
                call RewindBuffer(ObjEnterData, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR18'

                call ExtractBlockFromBuffer(ObjEnterData, ClientNumber,                 &
                                            BeginXX, EndXX, BlockFound,                 &
                                            FirstLine = FirstLine, LastLine = LastLine, &
                                            STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR19'

                if (BlockFound) then                                                  

                    line = FirstLine

                    do j = Me%WorkSize%JLB, Me%WorkSize%JUB+1 
                        line = line+1

                        !Last line found before end?
                        if (line == LastLine) then
                            write(*,*) 
                            write(*,*) 'Error reading XX'
                            stop 'ConstructGlobalVariables - HorizontalGrid - ERR20'
                        end if

                        call GetData(Me%XX(J), ObjEnterData, flag,            &
                                     Buffer_Line  = Line,                                    &
                                     STAT         = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR21'
                    end do
                end if

                call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR22'


            endif iconX


        end if BF

        !Closes Data File
        call KillEnterData      (ObjEnterData, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR28'
        
     
        !Allocates variables common to the module
        call AllocateVariables

    end subroutine ConstructGlobalVariables

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

    subroutine AllocateVariables()

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        integer                             :: STATUS
        integer                             :: ILB, IUB, JLB, JUB

        ILB = Me%Size%ILB
        IUB = Me%Size%IUB

        JLB = Me%Size%JLB
        JUB = Me%Size%JUB

        !LatitudeZ
        allocate (Me%LatitudeZ(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR3'

        !LongitudeZ
        allocate (Me%LongitudeZ(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR4'

        !Allocate Distances
        allocate (Me%DXX(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR5'

        allocate (Me%DYY(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR6'

        allocate (Me%DZX(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR7'

        allocate (Me%DZY(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR8'

        allocate (Me%DUX(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR9'

        allocate (Me%DUY(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR10'

        allocate (Me%DVX(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR11'

        allocate (Me%DVY(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR12'

        !Coriolis
        allocate (Me%F(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR13'

        allocate (Me%GridCellArea(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR13a'


        !Compute points location

        if (.not. Me%CornersXYInput) then

            !XX_Z
            allocate (Me%Compute%XX_Z(JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR14'

            !YY_Z
            allocate (Me%Compute%YY_Z(ILB:IUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR15'

            !XX_U
            allocate (Me%Compute%XX_U(JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR16'

            !YY_V
            allocate (Me%Compute%YY_V(ILB:IUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR17'

            !XX_V
            allocate (Me%Compute%XX_V(JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR18'

            !YY_U
            allocate (Me%Compute%YY_U(ILB:IUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR19'

            !XX_Cross
            allocate (Me%Compute%XX_Cross(JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR20'

            !YY_Cross
            allocate (Me%Compute%YY_Cross(ILB:IUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR21'

            Me%Compute%XX_Z     (:) = FillValueReal
            Me%Compute%YY_Z     (:) = FillValueReal
            Me%Compute%XX_U     (:) = FillValueReal
            Me%Compute%YY_V     (:) = FillValueReal
            Me%Compute%XX_V     (:) = FillValueReal
            Me%Compute%YY_U     (:) = FillValueReal
            Me%Compute%XX_Cross (:) = FillValueReal
            Me%Compute%YY_Cross (:) = FillValueReal

        else

            !XX2D_Z
            allocate (Me%Compute%XX2D_Z(ILB:IUB, JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR22'

            !YY2D_Z
            allocate (Me%Compute%YY2D_Z(ILB:IUB, JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR23'

            !XX2D_U
            allocate (Me%Compute%XX2D_U(ILB:IUB, JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR24'

            !YY2D_V
            allocate (Me%Compute%YY2D_V(ILB:IUB, JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR25'

            !XX2D_V
            allocate (Me%Compute%XX2D_V(ILB:IUB, JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR26'

            !YY2D_U
            allocate (Me%Compute%YY2D_U(ILB:IUB, JLB:JUB), stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR27'

            Me%Compute%XX2D_Z     (:,:) = FillValueReal
            Me%Compute%YY2D_Z     (:,:) = FillValueReal
            Me%Compute%XX2D_U     (:,:) = FillValueReal
            Me%Compute%YY2D_V     (:,:) = FillValueReal
            Me%Compute%XX2D_V     (:,:) = FillValueReal
            Me%Compute%YY2D_U     (:,:) = FillValueReal

        endif

        !Initialize Values
        Me%LatitudeZ    (:,:)    = FillValueReal
        Me%LongitudeZ   (:,:)    = FillValueReal
        Me%DXX          (:,:)    = FillValueReal
        Me%DYY          (:,:)    = FillValueReal
        Me%DZX          (:,:)    = FillValueReal
        Me%DZY          (:,:)    = FillValueReal
        Me%DUX          (:,:)    = FillValueReal
        Me%DUY          (:,:)    = FillValueReal
        Me%DVX          (:,:)    = FillValueReal
        Me%DVY          (:,:)    = FillValueReal
        Me%F            (:,:)    = FillValueReal
        Me%GridCellArea (:,:)    = FillValueReal

        if (Me%CoordType == CIRCULAR_ .or. Me%CornersXYInput) then

            Me%Distortion      = .true. 

            allocate (Me%RotationX(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR28'

            Me%RotationX(:,:) = FillValueReal

            allocate (Me%RotationY(ILB:IUB, JLB:JUB), stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'AllocateVariables - HorizontalGrid - ERR29'

            Me%RotationY(:,:) = FillValueReal

        endif


    end subroutine AllocateVariables


    !--------------------------------------------------------------------------
    !This subroutine needs to be reviewed for the case that the type of the coordinates
    !is equal to one (Geographic coordinates) - Frank Abr 99.

    subroutine Mercator()

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        real, dimension(:,:), pointer       :: XX_IE, YY_IE
        real, dimension(:  ), pointer       :: XX, YY
        real, dimension(:,:), allocatable   :: XX_aux, YY_aux
        integer                             :: AuxCoordTip
        real(8)                             :: DB_LAT, DB_LONG
        real                                :: X_PONTO, Y_PONTO
        real(8), dimension(:), pointer      :: DLZONE
        integer, dimension(:), pointer      :: IZONE1
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j, STATUS

        real(8)                             :: Rad_Lat, radians, CosenLat, Radius, & 
                                               EarthRadius, DX, DY, Angle 

        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB

        XX_IE => Me%XX_IE
        YY_IE => Me%YY_IE

        XX    => Me%XX
        YY    => Me%YY

        !By default 
do5 :   do j = Me%WorkSize%JLB, Me%WorkSize%JUB
do6 :   do i = Me%WorkSize%ILB, Me%WorkSize%IUB
           Me%LatitudeZ (i, j) = Me%Latitude
           Me%LongitudeZ(i, j) = Me%Longitude
        end do do6
        end do do5


        !Projected coordenates
        !Calculates the coordinates in the lower left corner (IE)

cd1:    if (Me%CoordType == UTM_        .or.    &
            Me%CoordType == MIL_PORT_   .or.    &
            Me%CoordType == GRID_COORD_ .or.    & 
            Me%CoordType == NLRD_) then

Inp:        if (.not. Me%CornersXYInput) then

do1 :           do j = JLB, JUB + 1
do2 :           do i = ILB, IUB + 1
                    XX_IE(I, J) = XX(J)
                    YY_IE(I, J) = YY(I)
                end do do2
                end do do1


                !Rotation of the points
do3 :           do j = JLB, JUB + 1
do4 :           do i = ILB, IUB + 1
                    call RODAXY(Me%Xorig, Me%Yorig,            &
                                Me%Grid_Angle, XX_IE(i, j), YY_IE(i, j))  
                end do do4
                end do do3

            endif Inp

            !Computes the Geographic Coordinates in the center od the cells (to use later, 
            !when computing Coriolis and the elevation of the sun
cd2:        if ((Me%CoordType == UTM_) .OR. (Me%CoordType == MIL_PORT_)) then 
                
do7 :           do j = JLB, JUB
do8 :           do i = ILB, IUB

                    !X_PONTO = (XX_IE(i  , j+1) + XX_IE(i, j)) / 2.
                    !Y_PONTO = (YY_IE(i+1, j  ) + YY_IE(i, j)) / 2.

                    if (Me%DefineCellsMap(i, j) == 1) then

                        Y_PONTO  = (YY_IE(i  , j+1) + YY_IE(i, j) + YY_IE(i+1, j+1) + YY_IE(i+1, j)) / 4.
                        X_PONTO  = (XX_IE(i  , j+1) + XX_IE(i, j) + XX_IE(i+1, j+1) + XX_IE(i+1, j)) / 4.
                    else
                        Y_PONTO  = 0
                        X_PONTO  = 0
                    endif

                    !Passagem para o referencial geografico
                    AuxCoordTip = Me%CoordType
                    if (Me%CoordType == MIL_PORT_) AuxCoordTip = 5 !coord portuguesas para geograficas

                    call USCONVCO (AuxCoordTip, Me%Zone, DB_LAT, DB_LONG,        & 
                                   dble(X_PONTO), dble(Y_PONTO))

                    Me%LatitudeZ(i, j)  = real(DB_LAT)
                    Me%LongitudeZ(i, j) = real(DB_LONG)

                end do do8
                end do do7
            
            elseif(Me%CoordType == NLRD_)then cd2
                
                
                do j = JLB, JUB
                do i = ILB, IUB

                    if (Me%DefineCellsMap(i, j) == 1) then
                        Y_PONTO  = (YY_IE(i  , j+1) + YY_IE(i, j) + YY_IE(i+1, j+1) + YY_IE(i+1, j)) / 4.
                        X_PONTO  = (XX_IE(i  , j+1) + XX_IE(i, j) + XX_IE(i+1, j+1) + XX_IE(i+1, j)) / 4.
                    else
                        Y_PONTO  = 0.
                        X_PONTO  = 0.
                    endif

                    call rd_to_wgs84(dble(X_PONTO), dble(Y_PONTO), DB_LAT, DB_LONG)


                    Me%LatitudeZ(i, j)  = real(DB_LAT)
                    Me%LongitudeZ(i, j) = real(DB_LONG)

                end do 
                end do
                

            endif cd2

        !Geographic coordinates
        elseif (Me%CoordType == GEOG_ .or.                                &
                Me%CoordType == SIMPLE_GEOG_) then cd1

cd33 :      if (Me%Grid_Angle .NE. 0.0) then

                write(*,*) 'The Geographic coordinates can not have the bathymetry rotated'
                stop  'Subroutine Mercator; module HorizontalGrid. ERR35'

            endif cd33

            !Computes the Geographic Coordinates in the center od the cells (to use later, 
            !when computing Coriolis and the elevation of the sun

            if (.not. Me%CornersXYInput) then

do9  :          do j = JLB, JUB
do10 :          do i = ILB, IUB
                   Me%LatitudeZ (i, j) = Me%Yorig + (YY(i+1) + YY(i)) / 2.0
                   Me%LongitudeZ(i, j) = Me%Xorig + (XX(j+1) + XX(j)) / 2.0
                end do do10
                end do do9

            else

do19 :          do j = JLB, JUB
do20 :          do i = ILB, IUB
                   Me%LatitudeZ (i, j) = (YY_IE(i  , j+1) + YY_IE(i, j) + YY_IE(i+1, j+1) + YY_IE(i+1, j)) / 4.
                   Me%LongitudeZ(i, j) = (XX_IE(i  , j+1) + XX_IE(i, j) + XX_IE(i+1, j+1) + XX_IE(i+1, j)) / 4.
                end do do20
                end do do19

            endif

cd36:       if (Me%CoordType == GEOG_) then

                
Inp3:           if (.not. Me%CornersXYInput) then


                    !Convertes geodesic coordenates into UTM
                    nullify (DLZONE, IZONE1)
                    allocate(DLZONE(Me%Size%ILB:Me%Size%IUB))
                    allocate(IZONE1(Me%Size%JLB:Me%Size%JUB))

                    call ConvertCoordenates(XX, YY, Me%Xorig, Me%Yorig, IZONE1, DLZONE)

                    radians = 180.0 / PI
                    do i = ILB, IUB + 1
                        XX_IE(i, 1) = 0.
                        YY_IE(i, 1) = (Me%Yorig + YY(i)) * 6378000. / radians
                        do j = JLB+1, JUB+1
                            XX_IE(i,j) = abs((XX(j) - XX(1))* DLZONE(i)) / 6.
                            YY_IE(i,j) = YY_IE(I,1)
                        enddo
                    enddo

                    deallocate(DLZONE, stat = STATUS)
                        if (STATUS /= SUCCESS_) stop 'Mercator - HorizontalGrid - ERR13'
                    deallocate(IZONE1, stat = STATUS)
                        if (STATUS /= SUCCESS_) stop 'Mercator - HorizontalGrid - ERR14'
                    nullify(DLZONE,IZONE1)


                else
                    write(*,*)'Cannot define curvilinear grid in GEOGRAPHIC coordinates.'
                    write(*,*)'Please select SIMPLE GEOGRAPHIC coordinates.'
                    stop 'Mercator - HorizontalGrid - ERR14a'
                end if Inp3


            else if (Me%CoordType == SIMPLE_GEOG_) then cd36

                radians      = 180.0 / PI
                EarthRadius  = 6378000.

Inp2:           if (.not. Me%CornersXYInput) then

                    YY_IE(ILB, JLB : JUB + 1) = 0.

                    do  i = ILB + 1, IUB + 1

                        DY = abs(YY(i) - YY(i-1))  * EarthRadius / radians
                    
                        YY_IE(i, JLB) = YY_IE(i - 1, JLB) + DY

                        do j = JLB+1, JUB + 1
                            YY_IE(i,j) = YY_IE(i, JLB)
                        enddo

                    enddo

                    XX_IE(ILB : IUB + 1, JLB) = 0.

                    do i = ILB, IUB + 1
                        do j = JLB + 1, JUB + 1
                            Rad_Lat  = (Me%Yorig + YY(i))/ radians 
                            CosenLat = cos(Rad_Lat)
                            Radius   = CosenLat * EarthRadius
                            DX       = abs (XX(j) - XX(j-1)) * Radius / radians

                            XX_IE(i, j) = XX_IE(i, j - 1) + DX    

                        enddo

                    enddo


                else Inp2

                    allocate (XX_aux(ILB-1:IUB+1,JLB-1:JUB+1))
                    allocate (YY_aux(ILB-1:IUB+1,JLB-1:JUB+1))

                    Me%LatitudeConn  (:, :) = YY_IE(:, :)
                    Me%LongitudeConn (:, :) = XX_IE(:, :)


                    YY_aux(ILB, JLB : JUB + 1) = 0.

                    do  i = ILB + 1, IUB + 1

                        DY = abs(YY_IE(i, JLB) - YY_IE(i-1, JLB))  * EarthRadius / radians
                    
                        YY_aux(i, JLB) = YY_aux(i - 1, JLB) + DY

                    enddo


                    XX_aux(ILB : IUB + 1, JLB) = 0.

                    do i = ILB, IUB + 1
                        do j = JLB+1, JUB + 1
                            Rad_Lat  = (YY_IE(i, j))/ radians 
                            CosenLat = cos(Rad_Lat)
                            Radius   = CosenLat * EarthRadius
                            DX       = abs (XX_IE(i,j) - XX_IE(i, j-1)) * Radius / radians
                            DY       = abs (YY_IE(i,j) - YY_IE(i, j-1)) * EarthRadius / radians


                            XX_aux(i, j) = XX_aux(i, j - 1) + DX    
                            YY_aux(i, j) = YY_aux(i, j - 1) + DY

                        enddo
                    enddo

                    XX_IE(:, :) = XX_aux(:, :)
                    YY_IE(:, :) = YY_aux(:, :)

                    deallocate (XX_aux)
                    deallocate (YY_aux)

                endif Inp2

            endif cd36

        !Circular coordinates
        elseif (Me%CoordType == CIRCULAR_) then cd1
            
            radians      = PI / 180.0 

            
            do  i = ILB, IUB + 1

                Angle  = Me%Yorig * radians + YY(i)* radians
           
                do  j = JLB, JUB + 1

                    Radius      = Me%Xorig + XX(j)

                    XX_IE(i, j) = Radius * cos (Angle) 
                    YY_IE(i, j) = Radius * sin (Angle) 

                enddo
            enddo




        else

            stop 'Invalid type of coordinates - Mercator - HorizontalGrid - ERR02'

        endif cd1



        !In the case that the model uses Coordinates different then Grid Coordinates
        !or Circular coordinates write them to file
        if (Me%CoordType == GEOG_ .or. Me%CoordType == SIMPLE_GEOG_) then

            if (.not. Me%CornersXYInput) then

                do j = JLB, JUB + 1
                do i = ILB, IUB + 1
                   Me%LatitudeConn  (i, j) = Me%YORIG + Me%YY(i)
                   Me%LongitudeConn (i, j) = Me%XORIG + Me%XX(j)
                end do
                end do

            end if

        elseif (Me%CoordType == UTM_ .or. Me%CoordType == MIL_PORT_) then

            do j = JLB, JUB + 1
            do i = ILB, IUB + 1

                X_PONTO = Me%XX_IE(i, j)
                Y_PONTO = Me%YY_IE(i, j)

                if (X_PONTO < FillValueReal/2. .or. Y_PONTO < FillValueReal/2.) then
                    X_PONTO = 0.
                    Y_PONTO = 0.
                endif
            
                !Passagem para o referencial geografico
                AuxCoordTip = Me%CoordType
                !coord portuguesas para geograficas
                if (Me%CoordType == MIL_PORT_) AuxCoordTip = 5 

                call USCONVCO (AuxCoordTip, Me%Zone, DB_LAT, DB_LONG, &
                               dble(X_PONTO), dble(Y_PONTO))

                Me%LatitudeConn (i, j) = real(DB_LAT)
                Me%LongitudeConn(i, j) = real(DB_LONG)


            end do
            end do

        elseif(Me%CoordType == NLRD_)then

            do j = JLB, JUB + 1
            do i = ILB, IUB + 1


                X_PONTO = Me%XX_IE(i, j)
                Y_PONTO = Me%YY_IE(i, j)

                if (X_PONTO == FillValueReal .or. Y_PONTO == FillValueReal) then
                    X_PONTO = 0.
                    Y_PONTO = 0.
                endif

                call rd_to_wgs84(dble(X_PONTO), dble(Y_PONTO), DB_LAT, DB_LONG)

                Me%LatitudeConn (i, j) = real(DB_LAT)
                Me%LongitudeConn(i, j) = real(DB_LONG)

            enddo
            enddo

        endif

        nullify(XX, YY)
        nullify(XX_IE, YY_IE)


    end subroutine Mercator

    !***********************************************************************
    !                                                                      *
    !  Esta subroutina calcula 2 vectores e um escalar para a conversão    *
    !  de coordenadas geodésicas em UTM:                                   *
    !                                                                      *
    !        comprimento do fuso:                    DLZONE(I)             *
    !        nº do fuso:                             IZONE1(J)             *
    !        nº do fuso de origem:                   IZONE_ORIG            *
    !                                                                      *
    !  O meridiano de Greenwich está no fuso 31. Considera-se como origem  *
    !  o ínicio do 1º fuso da malha. Cada fuso tem 6 graus e há 60 fusos.  *
    !                                                                      *
    !                                                                      *
    !    AIRES: 30/6/1977                                                  *
    !    Pedro Montero 20/11/97   Se pasan las variables por argumentos    *
    !***********************************************************************
    
    subroutine ConvertCoordenates(XX, YY, XORIG, YORIG, IZONE1, DLZONE)

        !Arguments-------------------------------------------------------------
        real, dimension(:), pointer         :: XX, YY  
        real                                :: XORIG, YORIG
        integer, dimension(:), pointer      :: IZONE1
        real(8), dimension(:), pointer      :: DLZONE  

        !Local-----------------------------------------------------------------

        integer :: i, j
        integer :: ILB, IUB, JLB, JUB
        integer :: IZONE_ORIG

        real    :: radians, DB_LAT, DB_LONG, DB_LONG_ORIG 

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

        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB
        
        radians = 180.0 / PI

do1:    do i = ILB, IUB+1
            DB_LAT       = YORIG + YY(i)
            DB_LONG_ORIG = XORIG
            IZONE_ORIG   = int((DB_LONG_ORIG + 360.) / 6) + 31 - 60
            IZONE1(1)    = IZONE_ORIG
            DLZONE(i)    = 6378000. * COS(DB_LAT / radians) * 6. / radians

do2:        do j = JLB+1, JUB+1
                DB_LONG   = XORIG + XX(j)
                IZONE1(j) = int((DB_LONG + 360.) / 6) + 31 - 60
            enddo do2
        enddo do1

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

    end subroutine ConvertCoordenates

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

    subroutine GridPointsLocation1D

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------

        real, dimension(:  ), pointer       :: XX, YY
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j
        real, dimension(:  ), pointer       :: XX_Z, YY_Z, XX_U , YY_V , XX_V, YY_U, XX_Cross, YY_Cross


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

        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB

        XX        => Me%XX 
        YY        => Me%YY 
        XX_Z      => Me%Compute%XX_Z
        YY_Z      => Me%Compute%YY_Z
        XX_U      => Me%Compute%XX_U
        YY_V      => Me%Compute%YY_V
        XX_V      => Me%Compute%XX_V
        YY_U      => Me%Compute%YY_U
        XX_Cross  => Me%Compute%XX_Cross
        YY_Cross  => Me%Compute%YY_Cross


do1 :   do j = JLB, JUB

            XX_Z(J) = (XX(J) + XX(J+1)) / 2.

            XX_V(J) = XX_Z(J)

        end do do1


do2 :   do i = ILB, IUB

            YY_Z(I) = (YY(I) + YY(I+1)) / 2.
            
            YY_U(I) = YY_Z(I)

        end do do2

do3 :   do j = JLB, JUB + 1

            XX_U(J)     = XX(J)

            XX_Cross(J) = XX(J) 

        end do do3

do4 :   do i = ILB, IUB + 1

            YY_V    (I) = YY(I)

            YY_Cross(I) = YY(I) 

        end do do4


        nullify(XX   )
        nullify(YY   )
        nullify(XX_Z )
        nullify(YY_Z )
        nullify(XX_U )
        nullify(YY_V )
        nullify(XX_V )
        nullify(YY_U )

        nullify(XX_Cross)
        nullify(YY_Cross)


    end subroutine GridPointsLocation1D

    !--------------------------------------------------------------------------
    subroutine GridPointsLocation2D

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------

        real, dimension(:,:), pointer       :: XX, YY
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j
        real, dimension(:,:), pointer       :: XX_Z, YY_Z, XX_U , YY_V , XX_V, YY_U


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

        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB


        XX        => Me%XX_IE 
        YY        => Me%YY_IE 
        XX_Z      => Me%Compute%XX2D_Z
        YY_Z      => Me%Compute%YY2D_Z
        XX_U      => Me%Compute%XX2D_U
        YY_V      => Me%Compute%YY2D_V
        XX_V      => Me%Compute%XX2D_V
        YY_U      => Me%Compute%YY2D_U

do1 :   do j = JLB, JUB
do2 :   do i = ILB, IUB

            XX_Z(I,J) = (XX(I,J) + XX(I,J+1) + XX(I+1,J) + XX(I+1,J+1)) / 4.
            YY_Z(I,J) = (YY(I,J) + YY(I,J+1) + YY(I+1,J) + YY(I+1,J+1)) / 4.

        end do do2
        end do do1


do5:   do j = JLB, JUB + 1
do6:   do i = ILB, IUB

            XX_U(I,J) = (XX(I,J) + XX(I+1,J)) / 2.
            YY_U(I,J) = (YY(I,J) + YY(I+1,J)) / 2.

        end do do6
        end do do5

do7:   do j = JLB, JUB
do8:   do i = ILB, IUB + 1

            XX_V(I,J) = (XX(I,J) + XX(I,J+1)) / 2.
            YY_V(I,J) = (YY(I,J) + YY(I,J+1)) / 2.

        end do do8
        end do do7



        nullify(XX   )
        nullify(YY   )
        nullify(XX_Z )
        nullify(YY_Z )
        nullify(XX_U )
        nullify(YY_V )
        nullify(XX_V )
        nullify(YY_U )



    end subroutine GridPointsLocation2D

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

    subroutine ComputeDistances

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j
        real(8)                             :: AuxDble

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

        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB
         

cd1:    if (Me%Grid_Angle /= 0. .or. Me%CoordType == CIRCULAR_ .or. Me%CornersXYInput) then

            !DXX
            do j = JLB, JUB
            do i = ILB, IUB + 1
                AuxDble =                                                &
                    sqrt( (dble(Me%XX_IE(i, j+1))    -    &
                           dble(Me%XX_IE(i, j)))**2. +    &
                          (dble(Me%YY_IE(i, j+1))    -    &
                           dble(Me%YY_IE(i, j)))**2. )

                Me%DXX(i, j) = real (AuxDble)
            enddo
            enddo

            !DYY
            do j = JLB, JUB + 1
            do i = ILB, IUB
                AuxDble =                                                &
                    sqrt( (dble(Me%XX_IE(i+1, j))    -    &
                           dble(Me%XX_IE(i, j)))**2. +    &
                          (dble(Me%YY_IE(i+1, j))    -    &
                           dble(Me%YY_IE(i, j)))**2. )
                           
                Me%DYY(i, j) = real (AuxDble)
            enddo
            enddo

        else  cd1 !Grid_Angle = 0. No grid rotation. Or coordinate type not circular

            !DXX
            do j = JLB, JUB
            do i = ILB, IUB + 1
                Me%DXX(i, j) = Me%XX_IE(i  , j+1) -    &
                                              Me%XX_IE(i  , j  )
            enddo
            enddo

            !DYY
            do j = JLB, JUB + 1
            do i = ILB, IUB 
                Me%DYY(i, j) = Me%YY_IE(i+1, j  ) -    &
                                              Me%YY_IE(i  , j  )
            enddo
            enddo

        
        endif cd1

        !DUX, DVY
        do j = JLB, JUB
        do i = ILB, IUB
            Me%DUX(i, j) = (Me%DXX(i, j) +       &
                                           Me%DXX(i+1, j)) / 2.
            Me%DVY(i, j) = (Me%DYY(i, j) +       &
                                           Me%DYY(i, j+1)) / 2.
        enddo
        enddo
    
        !DVX
        do j = JLB, JUB - 1
        do i = ILB, IUB + 1
           Me%DVX(i, j) = (Me%DXX(i, j) +        &
                                          Me%DXX(i, j+1)) / 2.
        enddo
        enddo

        !DZX
        do j = JLB, JUB - 1
        do i = ILB, IUB
            Me%DZX(i, j) = (Me%DUX(i, j) +       &
                                           Me%DUX(i, j+1)) / 2.
        enddo
        enddo

        !DUY
        do j = JLB, JUB + 1
        do i = ILB, IUB - 1
            Me%DUY(i, j) = (Me%DYY(i, j) +       &
                                           Me%DYY(i+1, j)) / 2.
            
        enddo
        enddo

        !DZY
        do j = JLB, JUB
        do i = ILB, IUB - 1
            Me%DZY(i, j) = (Me%DVY(i, j) +       &
                                           Me%DVY(i+1, j)) / 2.
        enddo
        enddo

    end subroutine ComputeDistances


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

    subroutine ComputeRotation

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j
        real(8)                             :: x1, x2, y1, y2, dx, dy

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

        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB
         

cd1:    if (Me%CornersXYInput .or. Me%CoordType == CIRCULAR_) then


            do j = JLB, JUB
            do i = ILB, IUB

                !Rotation X

                x1 = (Me%XX_IE(i, j)   + Me%XX_IE(i+1, j)) / 2.
                x2 = (Me%XX_IE(i, j+1) + Me%XX_IE(i+1, j+1)) / 2.

                y1 = (Me%YY_IE(i, j)   + Me%YY_IE(i+1, j)) / 2.
                y2 = (Me%YY_IE(i, j+1) + Me%YY_IE(i+1, j+1)) / 2.

                dx = x2 - x1
                dy = y2 - y1

                Me%RotationX(i, j) = atan2(dy, dx)

                !Rotation Y

                x1 = (Me%XX_IE(i,   j) + Me%XX_IE(i,   j+1)) / 2.
                x2 = (Me%XX_IE(i+1, j) + Me%XX_IE(i+1, j+1)) / 2.

                y1 = (Me%YY_IE(i,   j) + Me%YY_IE(i  , j+1)) / 2.
                y2 = (Me%YY_IE(i+1, j) + Me%YY_IE(i+1, j+1)) / 2.

                dx = x2 - x1
                dy = y2 - y1

                Me%RotationY(i, j) = atan2(dy, dx)



            enddo
            enddo

         endif cd1

    end subroutine ComputeRotation

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

    subroutine ComputeGridCellArea

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j
       
        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB

        do j = JLB, JUB
        do i = ILB, IUB

            Me%GridCellArea(i, j) = Me%DUX(i, j) * Me%DVY(i, j)

        enddo
        enddo
        

    end subroutine ComputeGridCellArea

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

    subroutine ComputeCoriolis

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j
        real                                :: OMEGA

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

        !Worksize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB

        !OMEGA
        OMEGA = 2. * PI / (24.016 * 3600.)

        do j = JLB, JUB
        do i = ILB, IUB
            Me%F(i, j) = 2. * OMEGA * sin(Me%LatitudeZ(i, j)*PI / 180.)
        end do
        end do

    end subroutine ComputeCoriolis


    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    !SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR 

    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    subroutine GetHorizontalGridSize(HorizontalGridID, Size, WorkSize, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        type (T_Size2D), optional, intent(OUT)      :: Size
        type (T_Size2D), optional, intent(OUT)      :: WorkSize
        integer,         optional, intent(OUT)      :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            !Size
            if (present(Size       )) Size     = Me%Size
            if (present(WorkSize   )) WorkSize = Me%WorkSize

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

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

    end subroutine GetHorizontalGridSize

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

    subroutine GetHorizontalGrid(HorizontalGridID, XX_IE, YY_IE, XX_Z, YY_Z,             &
                                 XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross,             &
                                 DXX, DYY, DZX, DZY, DUX, DUY, DVX, DVY, XX, YY, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real, dimension(:, :), pointer, optional    :: XX_IE, YY_IE
        real, dimension(:   ), pointer, optional    :: XX_Z, YY_Z, XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross
        real, dimension(:, :), pointer, optional    :: DXX, DYY, DZX, DZY
        real, dimension(:, :), pointer, optional    :: DUX, DUY, DVX, DVY
        real, dimension(:   ), pointer, optional    :: XX, YY
        integer, optional,  intent(OUT)             :: STAT    

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            !XX_IE, YY_IE
            if (present(XX_IE)) then
                XX_IE => Me%XX_IE
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(YY_IE)) then
                YY_IE => Me%YY_IE
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            !XX_Z, YY_Z
            if (present(XX_Z)) then
                XX_Z => Me%Compute%XX_Z
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(YY_Z)) then
                YY_Z => Me%Compute%YY_Z
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif


            !XX_U, YY_U
            if (present(XX_U)) then
                XX_U => Me%Compute%XX_U
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(YY_U)) then
                YY_U => Me%Compute%YY_U
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif


            !XX_V, YY_V
            if (present(XX_V)) then
                XX_V => Me%Compute%XX_V
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(YY_V)) then
                YY_V => Me%Compute%YY_V
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            !XX_Cross, YY_Cross
            if (present(XX_Cross)) then
                XX_Cross => Me%Compute%XX_Cross
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(YY_Cross)) then
                YY_Cross => Me%Compute%YY_Cross
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif


            !DXX, DYY
            if (present(DXX)) then
                DXX => Me%DXX
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(DYY)) then
                DYY => Me%DYY
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            !DZX, DZY
            if (present(DZX)) then
                DZX => Me%DZX
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(DZY)) then
                DZY => Me%DZY
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            !DUX, DUY, DVX, DVY
            if (present(DUX)) then
                DUX => Me%DUX
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(DUY)) then
                DUY => Me%DUY
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(DVX)) then
                DVX => Me%DVX
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(DVY)) then
                DVY => Me%DVY
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            !XX, YY
            if (present(XX)) then
                XX => Me%XX
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif

            if (present(YY)) then
                YY => Me%YY
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            endif


            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) &
            STAT = STAT_

    end subroutine GetHorizontalGrid
    
    !--------------------------------------------------------------------------        

    subroutine GetGridCellArea (HorizontalGridID, GridCellArea, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real, dimension(:, :), pointer, optional    :: GridCellArea
        integer, optional,  intent(OUT)             :: STAT    

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            GridCellArea => Me%GridCellArea
            call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)


            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_

    end subroutine GetGridCellArea

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

    subroutine GetGridOrigin(HorizontalGridID, Xorig, Yorig, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real,              intent(OUT)              :: Xorig, Yorig
        integer, optional, intent(OUT)              :: STAT


        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            Xorig = Me%Xorig
            Yorig = Me%Yorig

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_

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

    end subroutine GetGridOrigin

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

    subroutine GetGridAngle(HorizontalGridID, Angle, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real,              intent(OUT)              :: Angle
        integer, optional, intent(OUT)              :: STAT


        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            Angle = Me%Grid_Angle

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_

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

    end subroutine GetGridAngle

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

    subroutine GetLatitudeLongitude(HorizontalGridID, Latitude, Longitude, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real,    optional, intent(OUT)              :: Latitude
        real,    optional, intent(OUT)              :: Longitude
        integer, optional, intent(OUT)              :: STAT


        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            if(present(Latitude)) Latitude  = Me%Latitude
            if(present(Longitude))Longitude = Me%Longitude

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_

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

    end subroutine GetLatitudeLongitude

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

    subroutine GetGridZone(HorizontalGridID, Zone, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        integer,           intent(OUT)              :: Zone
        integer, optional, intent(OUT)              :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            Zone = Me%Zone

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_

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

    end subroutine GetGridZone

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

    subroutine GetGridCoordType(HorizontalGridID, CoordType, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        integer,           intent(OUT)              :: CoordType
        integer, optional, intent(OUT)              :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            CoordType = Me%CoordType

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_

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

    end subroutine GetGridCoordType

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

    subroutine GetCoordTypeList(GEOG, UTM, MIL_PORT, SIMPLE_GEOG, GRID_COORD, CIRCULAR, NLRD, &
                                PORTUGUESE_UTM_ZONE)
        
        !Arguments-------------------------------------------------------------

        integer, optional, intent(OUT) :: GEOG                  !Coordenadas Geograficas       
        integer, optional, intent(OUT) :: SIMPLE_GEOG           !Coordenadas Geograficas Simplificadas     
        integer, optional, intent(OUT) :: UTM                   !Coordenadas (UTM)
        integer, optional, intent(OUT) :: MIL_PORT              !Coordenadas Militares Portuguesas
        integer, optional, intent(OUT) :: GRID_COORD            !Coordenadas da malha
        integer, optional, intent(OUT) :: CIRCULAR              !Coordenadas circulares (XX - raio, YY - anglo em graus)
        integer, optional, intent(OUT) :: NLRD                  !Coordenadas Netherlands RD
        integer, optional, intent(OUT) :: PORTUGUESE_UTM_ZONE   !

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

        if (present(GEOG                )) GEOG                 = GEOG_
        if (present(SIMPLE_GEOG         )) SIMPLE_GEOG          = SIMPLE_GEOG_
        if (present(UTM                 )) UTM                  = UTM_
        if (present(MIL_PORT            )) MIL_PORT             = MIL_PORT_
        if (present(GRID_COORD          )) GRID_COORD           = GRID_COORD_
        if (present(CIRCULAR            )) CIRCULAR             = CIRCULAR_
        if (present(NLRD                )) NLRD                 = NLRD_
        if (present(PORTUGUESE_UTM_ZONE )) PORTUGUESE_UTM_ZONE  = PORTUGUESE_UTM_ZONE_
                                                                                                 
        !----------------------------------------------------------------------

    end subroutine GetCoordTypeList

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


    subroutine GetGridMeanLatLong(HorizontalGridID, Latitude, Longitude, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real,              intent(OUT)              :: Latitude, Longitude
        integer, optional, intent(OUT)              :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
cd1 :   if ((ready_ .EQ. IDLE_ERR_     ) .OR. &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            Latitude  = Me%Latitude
            Longitude = Me%Longitude

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_

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

    end subroutine GetGridMeanLatLong

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

    subroutine GetCoriolisFrequency(HorizontalGridID, F, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real, dimension(:, :), pointer              :: F
        integer, optional,  intent(OUT)             :: STAT    

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

        if ((ready_ .EQ. IDLE_ERR_     ) .OR.                                            &
            (ready_ .EQ. READ_LOCK_ERR_)) then
            call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            F => Me%F
            STAT_ = SUCCESS_

        else 

            STAT_ = ready_

        end if

        if (present(STAT)) STAT = SUCCESS_

    end subroutine GetCoriolisFrequency

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

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

    subroutine GetDefineCellsMap(HorizontalGridID, DefineCellsMap, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        integer, dimension(:, :), pointer           :: DefineCellsMap
        integer, optional,  intent(OUT)             :: STAT    

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

        if ((ready_ .EQ. IDLE_ERR_     ) .OR.                                            &
            (ready_ .EQ. READ_LOCK_ERR_)) then
            call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            DefineCellsMap => Me%DefineCellsMap
            STAT_ = SUCCESS_

        else 

            STAT_ = ready_

        end if

        if (present(STAT)) STAT = SUCCESS_

    end subroutine GetDefineCellsMap

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



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

    subroutine GetCheckDistortion(HorizontalGridID, Distortion, CornersXYInput, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        logical          ,  intent(OUT)             :: Distortion
        logical, optional,  intent(OUT)             :: CornersXYInput
        integer, optional,  intent(OUT)             :: STAT    

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_

        !Begin-----------------------------------------------------------------
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

        if ((ready_ .EQ. IDLE_ERR_     ) .OR.                                            &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            if (present(CornersXYInput)) CornersXYInput = Me%CornersXYInput

            Distortion = Me%Distortion

            STAT_ = SUCCESS_

        else 

            STAT_ = ready_

        end if

        if (present(STAT)) STAT = SUCCESS_

    end subroutine GetCheckDistortion

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

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

    subroutine GetGridRotation(HorizontalGridID, RotationX, RotationY, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real, dimension(:, :), pointer              :: RotationX, RotationY
        integer, optional,  intent(OUT)             :: STAT    

        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_

        !Begin-----------------------------------------------------------------
        
        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

        if ((ready_ .EQ. IDLE_ERR_     ) .OR.                                            &
            (ready_ .EQ. READ_LOCK_ERR_)) then

            call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            RotationX  => Me%RotationX
             
            call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
            RotationY  => Me%RotationY

            STAT_ = SUCCESS_                

        else 

            STAT_ = ready_

        end if

        if (present(STAT)) STAT = SUCCESS_

    end subroutine GetGridRotation

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


    subroutine GetGridLatitudeLongitude(HorizontalGridID,                       &
                                        GridLatitude,                           &
                                        GridLongitude,                          &
                                        GridLatitudeConn,                       &
                                        GridLongitudeConn,                      &
                                        STAT)

        !Arguments---------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real, dimension(:,:), optional, pointer     :: GridLongitude, GridLatitude
        real, dimension(:,:), optional, pointer     :: GridLongitudeConn, GridLatitudeConn
        integer, optional,  intent(OUT)             :: STAT    

        !Local-------------------------------------------------------------------
        integer                                     :: STAT_             
        integer                                     :: ready_              

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

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
        if ((ready_ == IDLE_ERR_     ) .OR.                                     &
            (ready_ == READ_LOCK_ERR_)) then

cd1 :       if (present(GridLongitude)) then
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
                GridLongitude => Me%LongitudeZ
            end if cd1
            

cd2 :       if (present(GridLatitude)) then
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
                GridLatitude => Me%LatitudeZ
            end if cd2

cd3 :       if (present(GridLatitudeConn)) then
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
                GridLatitudeConn => Me%LatitudeConn
            end if cd3


cd4 :       if (present(GridLongitudeConn)) then
                call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
                GridLongitudeConn => Me%LongitudeConn
            end if cd4

            
            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if

        if (present(STAT))                                                      &
            STAT = STAT_

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

    end subroutine GetGridLatitudeLongitude

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

    subroutine GetComputeZUV(HorizontalGridID, ComputeZ, ComputeU, ComputeV,             &
                             ComputeCross,  ComputeZU, ComputeZV, STAT)

        !Arguments---------------------------------------------------------------
        integer                                     :: HorizontalGridID
        integer, optional,  intent(OUT)             :: ComputeZ, ComputeU, ComputeV
        integer, optional,  intent(OUT)             :: ComputeCross, ComputeZU, ComputeZV
        integer, optional,  intent(OUT)  :: STAT    

        !Local-------------------------------------------------------------------
        integer                                     :: STAT_             
        integer                                     :: ready_              
        
        !------------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
        if ((ready_ == IDLE_ERR_     ) .OR.                                              & 
            (ready_ == READ_LOCK_ERR_)) then

cd1 :       if (present(ComputeZ)) then
                ComputeZ = ComputeZ_
            end if cd1
            
cd2 :       if (present(ComputeU)) then
                ComputeU = ComputeU_
            end if cd2

cd3 :       if (present(ComputeV)) then
                ComputeV = ComputeV_
            end if cd3

cd4 :       if (present(ComputeCross)) then
                ComputeCross = ComputeCross_
            end if cd4

cd5 :       if (present(ComputeZU)) then
                ComputeZU = ComputeZU_
            end if cd5

cd6 :       if (present(ComputeZV)) then
                ComputeZV = ComputeZV_
            end if cd6
            
            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if

        if (present(STAT))                                                      &
            STAT = STAT_

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

    end subroutine GetComputeZUV

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

    subroutine GetFatherGridID(HorizontalGridID, LargeScaleModel, Assimila, SurfaceMM5, STAT)

        !Arguments---------------------------------------------------------------
        integer                                     :: HorizontalGridID
        integer, optional,  intent(OUT)             :: LargeScaleModel, Assimila, SurfaceMM5
        integer, optional,  intent(OUT)             :: STAT    
 
        !Local-------------------------------------------------------------------
        integer                                     :: STAT_             
        integer                                     :: ready_              
        
        !------------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
        if ((ready_ == IDLE_ERR_     ) .OR.                                              & 
            (ready_ == READ_LOCK_ERR_)) then

cd1 :       if (present(LargeScaleModel)) then
                LargeScaleModel = LargeScaleModel_
            end if cd1
            
cd2 :       if (present(Assimila)) then
                Assimila = Assimila_
            end if cd2

cd3 :       if (present(SurfaceMM5)) then
                SurfaceMM5 = SurfaceMM5_
            end if cd3

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if

        if (present(STAT))                                                      &
            STAT = STAT_

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

    end subroutine GetFatherGridID

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

    subroutine GetGridFileName (HorizontalGridID, FileName, STAT)

        !Arguments---------------------------------------------------------------
        integer                                     :: HorizontalGridID
        character(len=*)                            :: FileName
        integer, optional,  intent(OUT)             :: STAT    
 
        !Local-------------------------------------------------------------------
        integer                                     :: STAT_             
        integer                                     :: ready_              
        
        !------------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    
        
        if ((ready_ == IDLE_ERR_     ) .OR.                                              & 
            (ready_ == READ_LOCK_ERR_)) then
    
            FileName = Me%FileName
            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if

        if (present(STAT)) STAT = STAT_

    end subroutine GetGridFileName

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

    subroutine UngetHorizontalGrid1D(HorizontalGridID, Array, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real, pointer, dimension(:)                 :: Array
        integer, optional, intent(OUT)              :: STAT
   
        !Local-------------------------------------------------------------------
        integer                                     :: STAT_             
        integer                                     :: ready_              

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

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

        if ((ready_ == IDLE_ERR_     ) .OR.                                              & 
            (ready_ == READ_LOCK_ERR_)) then

            nullify(Array)
            call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UngetHorizontalGrid1D")

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if 


        if (present(STAT)) STAT = STAT_

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

    end subroutine UngetHorizontalGrid1D

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

    subroutine UngetHorizontalGrid2D(HorizontalGridID, Array, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        real, pointer, dimension(:,:)               :: Array
        integer, optional, intent(OUT)              :: STAT
   
        !Local-------------------------------------------------------------------
        integer                                     :: STAT_             
        integer                                     :: ready_              

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

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

        if ((ready_ == IDLE_ERR_     ) .OR.                                              & 
            (ready_ == READ_LOCK_ERR_)) then

            nullify(Array)

            call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UngetHorizontalGrid2D")

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if

        if (present(STAT)) STAT = STAT_

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

    end subroutine UngetHorizontalGrid2D

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

    subroutine UngetHorizontalGrid2DInt(HorizontalGridID, Array, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        integer, pointer, dimension(:,:)            :: Array
        integer, optional, intent(OUT)              :: STAT
   
        !Local-------------------------------------------------------------------
        integer                                     :: STAT_             
        integer                                     :: ready_              

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

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

        if ((ready_ == IDLE_ERR_     ) .OR.                                              & 
            (ready_ == READ_LOCK_ERR_)) then

            nullify(Array)

            call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UngetHorizontalGrid2DInt")

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if

        if (present(STAT)) STAT = STAT_

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

    end subroutine UngetHorizontalGrid2DInt


    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    !MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER  

    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

    subroutine InterpolRegularGrid2D(HorizontalGridSonID, HorizontalGridFatherID,        &
                                     Field2DFather, Field2DSon,                          &
                                     ComputeFather, Compute,                             &
                                     KUBFather, GridID, STAT)



        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridSonID
        integer                                     :: HorizontalGridFatherID
        real,    dimension(:,:  ), pointer          :: Field2DFather, Field2DSon
        integer, dimension(:,:,:), pointer, optional:: ComputeFather
        integer, optional, intent(IN)               :: Compute, KUBFather
        integer, optional, intent(IN)               :: GridID
        integer, optional, intent(OUT)              :: STAT

        !Local-----------------------------------------------------------------
        type(T_HorizontalGrid),    pointer          :: ObjHorizontalGridFather  
        type (T_FatherGrid),       pointer          :: FatherGrid
        real   , pointer, dimension(:,:)            :: XXSon, YYSon, DXFather, DYFather, XXFather2D, YYFather2D
        real   , pointer, dimension(:  )            :: XXFather, YYFather
        integer, pointer, dimension(:,:)            :: JSon, ISon


        real    :: YYUpper, YYLower, XXUpper, XXLower,                                   &
                   PropLowLeft,  PropUpLeft, PropLowRight, PropUpRight
        
        integer :: ONLowLeft, ONUpLeft, ONLowRight, ONUpRight
        integer :: ready_, GridID_
            
        integer :: JLBSon, JUBSon, ILBSon, IUBSon
        integer :: JLBFather, JUBFather, ILBFather, IUBFather
        integer :: Jlower, Jupper, Ilower, Iupper 
        integer :: STAT_            

        integer :: i, j
        integer :: Compute_, KUBFather_ 
        logical :: InterPolOK, MapFatherGrid

        real    :: XPosition, YPosition
        integer :: Dij, Dji

        !Begin------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridSonID, ready_)    

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then

             
            if (.NOT. associated(Field2DFather))                                         &
                call SetError(FATAL_, INTERNAL_,                                           &
                             "InterpolRegularGrid2D; HorizontalGrid. ERR01") 

            if (.NOT. associated(Field2DSon))                                            &
                call SetError(FATAL_, INTERNAL_,                                           &
                             "InterpolRegularGrid2D; HorizontalGrid. ERR02") 

            nullify(ObjHorizontalGridFather)
            call LocateObjFather (ObjHorizontalGridFather, HorizontalGridFatherID)

            JLBSon = Me%WorkSize%JLB
            JUBSon = Me%WorkSize%JUB
            ILBSon = Me%WorkSize%ILB
            IUBSon = Me%WorkSize%IUB

            if (present(GridID)) then 
                GridID_ = GridID
            else
                GridID_ = LargeScaleModel_
            endif

            if (present(ComputeFather)) then 
                MapFatherGrid = .true.
            else
                MapFatherGrid = .false.
            endif

            if (present(Compute)) then 
                Compute_ = Compute
            else
                Compute_ = ComputeZ_
            endif

            if (present(KUBFather)) then 
                KUBFather_ = KUBFather
            else
                KUBFather_ = 1
            endif

            call Search_FatherGrid (FatherGrid, GridID_)
            

            !Gets the bounds from the Father
            ILBFAther = ObjHorizontalGridFather%WorkSize%ILB 
            IUBFather = ObjHorizontalGridFather%WorkSize%IUB
            JLBFather = ObjHorizontalGridFather%WorkSize%JLB
            JUBFather = ObjHorizontalGridFather%WorkSize%JUB

            if      (Compute_ == ComputeU_)  then
   
                JUBFather = JUBFather + 1
                JUBSon    = JUBSon    + 1

                XXSon    => FatherGrid%XX_U
                YYSon    => FatherGrid%YY_U

                JSon     => FatherGrid%JU
                ISon     => FatherGrid%IU

                XXFather => ObjHorizontalGridFather%Compute%XX_U
                YYFather => ObjHorizontalGridFather%Compute%YY_U

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_U
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_U

                DXFather => ObjHorizontalGridFather%DXX
                DYFather => ObjHorizontalGridFather%DZY

                Dij = 0
                Dji = 1

            else if (Compute_ == ComputeV_)  then

                IUBFather = IUBFather + 1
                IUBSon    = IUBSon    + 1

                XXSon => FatherGrid%XX_V
                YYSon => FatherGrid%YY_V

                JSon     => FatherGrid%JV
                ISon     => FatherGrid%IV

                XXFather => ObjHorizontalGridFather%Compute%XX_V
                YYFather => ObjHorizontalGridFather%Compute%YY_V

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_V
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_V

                DXFather => ObjHorizontalGridFather%DZX
                DYFather => ObjHorizontalGridFather%DYY

                Dij = 1
                Dji = 0 

            else if (Compute_ == ComputeZ_)  then

                XXSon => FatherGrid%XX_Z
                YYSon => FatherGrid%YY_Z

                XXFather => ObjHorizontalGridFather%Compute%XX_Z
                YYFather => ObjHorizontalGridFather%Compute%YY_Z

                JSon     => FatherGrid%JZ
                ISon     => FatherGrid%IZ

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_Z
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_Z

                DXFather => ObjHorizontalGridFather%DVX
                DYFather => ObjHorizontalGridFather%DUY

                Dji = 1
                Dij = 1

            else if (Compute_ == ComputeCross_)  then

                IUBFather = IUBFather + 1
                JUBFather = JUBFather + 1

                IUBSon    = IUBSon

                XXSon => FatherGrid%XX_Z
                YYSon => FatherGrid%YY_Z

                XXFather => ObjHorizontalGridFather%Compute%XX_Cross
                YYFather => ObjHorizontalGridFather%Compute%YY_Cross

                JSon     => FatherGrid%JCross
                ISon     => FatherGrid%ICross

                XXFather2D => ObjHorizontalGridFather%XX_IE
                YYFather2D => ObjHorizontalGridFather%YY_IE

                DXFather => ObjHorizontalGridFather%DUX
                DYFather => ObjHorizontalGridFather%DVY

                Dji = 0
                Dij = 0

            endif

            if (.not.MapFatherGrid) then

                ONLowLeft   = 1
                ONUpLeft    = 1
                ONLowRight  = 1
                ONUpRight   = 1

            endif


doj:        do j = JLBSon, JUBSon
doi:        do i = ILBSon, IUBSon

                Jlower = JSon(i, j)
                Jupper = JSon(i, j) + 1

                Ilower = ISon(i, j)
                Iupper = ISon(i, j) + 1

                if (Jupper > JUBFather) Jupper = JUBFather

                if (Iupper > IUBFather) Iupper = IUBFather

                if (Jlower < JLBFather) Jlower = JLBFather

                if (Ilower < ILBFather) Ilower = ILBFather

                if (.not. ObjHorizontalGridFather%CornersXYInput) then

                    XXLower  = XXFather(Jlower)
                    XXUpper  = XXFather(Jupper)

                    YYLower  = YYFather(Ilower)
                    YYUpper  = YYFather(Iupper)

                    XPosition = XXSon(i, j)
                    YPosition = YYSon(i, j)

                else

                    call CellReferential(DXFather, DYFather, XXFather2D, YYFather2D,    &
                                         ObjHorizontalGridFather%RotationX,             &
                                         ObjHorizontalGridFather%RotationY,             &
                                         XXSon(i,j), YYSon(i,j),                        &
                                         Ilower, Iupper, Jlower, Jupper, dij, dji,      &
                                         XPosition, YPosition,                          &
                                         YYUpper, YYLower, XXUpper, XXLower)

                endif

                PropLowLeft = Field2DFather(Ilower, Jlower)
                PropUpLeft  = Field2DFather(Iupper, Jlower)
                PropLowRight= Field2DFather(Ilower, Jupper)
                PropUpRight = Field2DFather(Iupper, Jupper)           

                if (MapFatherGrid) then

                    ONLowLeft   = ComputeFather(Ilower, Jlower, KUBFather_)
                    ONUpLeft    = ComputeFather(Iupper, Jlower, KUBFather_)
                    ONLowRight  = ComputeFather(Ilower, Jupper, KUBFather_)
                    ONUpRight   = ComputeFather(Iupper, Jupper, KUBFather_)

                endif

                call InterpolPoint(XPosition, YPosition,                             &
                                   YYUpper, YYLower, XXUpper, XXLower,               &
                                   PropLowLeft, PropUpLeft, PropLowRight, PropUpRight,&
                                   ONLowLeft,   ONUpLeft,   ONLowRight,   ONUpRight, &
                                   Field2DSon(i, j), InterPolOK)

ifc:            if (.not. InterPolOK)   then
                    
                    Field2DSon(i, j) = FillValueReal

                endif ifc

            
            enddo doi
            enddo doj

            nullify(XXSon,    YYSon   )
            nullify(JSon,     ISon    )
            nullify(XXFather, YYFather)

           
            STAT_ = SUCCESS_
        else               
            STAT_ = ready_
        end if cd1

        if (present(STAT))                                                    &
            STAT = STAT_

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

    end subroutine InterpolRegularGrid2D

    !--------------------------------------------------------------------------
 
    subroutine InterpolRegularGrid3D(HorizontalGridSonID, HorizontalGridFatherID,         &
                                      Field3DFather, Field3DSon,                          &
                                      ComputeFather, Compute,                             &
                                      KLBFather, KUBFather, KUBSon,                       & 
                                      FluxType, GridID, STAT)



        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridSonID
        integer                                     :: HorizontalGridFatherID
        real(4), dimension(:,:,:), pointer          :: Field3DFather, Field3DSon
        integer, dimension(:,:,:), pointer          :: ComputeFather
        integer,           intent(IN)               :: Compute, KLBFather, KUBFather, KUBSon
        logical, optional, intent(IN)               :: FluxType
        integer, optional, intent(IN)               :: GridID
        integer, optional, intent(OUT)              :: STAT

        !Local-----------------------------------------------------------------
        type(T_HorizontalGrid),    pointer          :: ObjHorizontalGridFather  
        type (T_FatherGrid),       pointer          :: FatherGrid
        real,    dimension(:,:  ), pointer          :: DXFather, DYFather, XXFather2D, YYFather2D
        real   , pointer, dimension(:,:)            :: XXSon, YYSon, DYXFather
        real   , pointer, dimension(:  )            :: XXFather, YYFather
        integer, pointer, dimension(:,:)            :: JSon, ISon

        real    :: YYUpper, YYLower, XXUpper, XXLower,                                   &
                   PropLowLeft,  PropUpLeft, PropLowRight, PropUpRight, PropSon
        
        integer :: ONLowLeft, ONUpLeft, ONLowRight, ONUpRight
        integer :: ready_, GridID_
            
        integer :: JLBSon, JUBSon, ILBSon, IUBSon
        integer :: JLBFather, JUBFather, ILBFather, IUBFather
        integer :: Jlower, Jupper, Ilower, Iupper 
        integer :: STAT_            

        integer :: i, j, k
        logical :: InterPolOK, FluxType_

        real    :: XPosition, YPosition
        integer :: Dij, Dji


        !Begin------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridSonID, ready_)

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then

            if (.NOT. associated(Field3DFather))                                         &
                call SetError(FATAL_, INTERNAL_,                                           &
                             "InterpolRegularGrid3D; HorizontalGrid. ERR01") 

            if (.NOT. associated(Field3DSon))                                            &
                call SetError(FATAL_, INTERNAL_,                                           &
                             "InterpolRegularGrid3D; HorizontalGrid. ERR02") 

            nullify(ObjHorizontalGridFather)
            call LocateObjFather (ObjHorizontalGridFather, HorizontalGridFatherID)

            JLBSon = Me%WorkSize%JLB
            JUBSon = Me%WorkSize%JUB
            ILBSon = Me%WorkSize%ILB
            IUBSon = Me%WorkSize%IUB

            if (present(GridID)) then 

                GridID_ = GridID

            else

                GridID_ = LargeScaleModel_

            endif

            call Search_FatherGrid (FatherGrid, GridID_)
 

            !Gets the bounds from the Father
            ILBFAther = ObjHorizontalGridFather%WorkSize%ILB 
            IUBFather = ObjHorizontalGridFather%WorkSize%IUB
            JLBFather = ObjHorizontalGridFather%WorkSize%JLB
            JUBFather = ObjHorizontalGridFather%WorkSize%JUB


            if (present(FluxType)) then
                FluxType_ = FluxType
            else
                FluxType_ = .false.
            endif

            if      (Compute == ComputeU_)  then

                JUBFather = JUBFather + 1
                JUBSon    = JUBSon    + 1

                XXSon    => FatherGrid%XX_U
                YYSon    => FatherGrid%YY_U

                JSon     => FatherGrid%JU
                ISon     => FatherGrid%IU


                XXFather => ObjHorizontalGridFather%Compute%XX_U
                YYFather => ObjHorizontalGridFather%Compute%YY_U

                DYXFather=> ObjHorizontalGridFather%DYY

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_U
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_U

                DXFather => ObjHorizontalGridFather%DXX
                DYFather => ObjHorizontalGridFather%DZY

                Dij = 0
                Dji = 1

            else if (Compute == ComputeV_)  then

                IUBFather = IUBFather + 1
                IUBSon    = IUBSon    + 1

                XXSon => FatherGrid%XX_V
                YYSon => FatherGrid%YY_V

                JSon     => FatherGrid%JV
                ISon     => FatherGrid%IV

                XXFather => ObjHorizontalGridFather%Compute%XX_V
                YYFather => ObjHorizontalGridFather%Compute%YY_V

                DYXFather=> ObjHorizontalGridFather%DXX

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_V
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_V

                DXFather => ObjHorizontalGridFather%DZX
                DYFather => ObjHorizontalGridFather%DYY

                Dij = 1
                Dji = 0 


            else if (Compute == ComputeZ_)  then

                if (FluxType_) call SetError(FATAL_, INTERNAL_,                            &
                             "InterpolRegularGrid3D; HorizontalGrid. ERR04") 

                XXSon => FatherGrid%XX_Z
                YYSon => FatherGrid%YY_Z

                XXFather => ObjHorizontalGridFather%Compute%XX_Z
                YYFather => ObjHorizontalGridFather%Compute%YY_Z

                JSon     => FatherGrid%JZ
                ISon     => FatherGrid%IZ

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_Z
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_Z

                DXFather => ObjHorizontalGridFather%DVX
                DYFather => ObjHorizontalGridFather%DUY

                Dji = 1
                Dij = 1

            endif


            !INSERT CODE HERE
dok:        do k = KLBFather   , KUBFather
doj:        do j = JLBSon, JUBSon
doi:        do i = ILBSon, IUBSon
                
                InterPolOK    = .false.

                Jlower = JSon(i, j)
                Jupper = JSon(i, j) + 1

                Ilower = ISon(i, j)
                Iupper = ISon(i, j) + 1

                if (Jupper > JUBFather) Jupper = JUBFather

                if (Iupper > IUBFather) Iupper = IUBFather

                if (Jlower < JLBFather) Jlower = JLBFather

                if (Ilower < ILBFather) Ilower = ILBFather

                if (.not. ObjHorizontalGridFather%CornersXYInput) then

                    XXLower  = XXFather(Jlower)
                    XXUpper  = XXFather(Jupper)

                    YYLower  = YYFather(Ilower)
                    YYUpper  = YYFather(Iupper)

                    XPosition = XXSon(i,j)
                    YPosition = YYSon(i,j)


                else

                    call CellReferential(DXFather, DYFather, XXFather2D, YYFather2D,    &
                                         ObjHorizontalGridFather%RotationX,             &
                                         ObjHorizontalGridFather%RotationY,             &
                                         XXSon(i,j), YYSon(i,j),                        &
                                         Ilower, Iupper, Jlower, Jupper, dij, dji,      &
                                         XPosition, YPosition,                          &
                                         YYUpper, YYLower, XXUpper, XXLower)

                endif


                if (FluxType_) then 

                    PropLowLeft = Field3DFather(Ilower, Jlower, k) / DYXFather(Ilower, Jlower)
                    PropUpLeft  = Field3DFather(Iupper, Jlower, k) / DYXFather(Iupper, Jlower)
                    PropLowRight= Field3DFather(Ilower, Jupper, k) / DYXFather(Ilower, Jupper)
                    PropUpRight = Field3DFather(Iupper, Jupper, k) / DYXFather(Iupper, Jupper)

                else

                    PropLowLeft = Field3DFather(Ilower, Jlower, k)
                    PropUpLeft  = Field3DFather(Iupper, Jlower, k)
                    PropLowRight= Field3DFather(Ilower, Jupper, k)
                    PropUpRight = Field3DFather(Iupper, Jupper, k)           

                endif

                ONLowLeft   = ComputeFather(Ilower, Jlower, k)
                ONUpLeft    = ComputeFather(Iupper, Jlower, k)
                ONLowRight  = ComputeFather(Ilower, Jupper, k)
                ONUpRight   = ComputeFather(Iupper, Jupper, k)


                call InterpolPoint(XPosition, YPosition,                             &
                                   YYUpper, YYLower, XXUpper, XXLower,               &
                                   PropLowLeft, PropUpLeft, PropLowRight, PropUpRight,&
                                   ONLowLeft,   ONUpLeft,   ONLowRight,   ONUpRight, &
                                   PropSon, InterPolOK)

                if (.not. InterPolOK) then

                    PropSon      = FillValueReal

                endif

                if (KUBFather > 1) then

                    Field3DSon(i, j, k) = PropSon


                else

                    Field3DSon(i, j, KUBSon) = PropSon

                endif


            
            enddo doi
            enddo doj
            enddo dok


            nullify(XXSon,    YYSon   )
            nullify(JSon,     ISon    )
            nullify(XXFather, YYFather)
            nullify(DYXFather         )

           
            STAT_ = SUCCESS_
        else               
            STAT_ = ready_
        end if cd1


        if (present(STAT))                                                    &
            STAT = STAT_

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

    end subroutine InterpolRegularGrid3D


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

    subroutine InterpolRegularGrid3D8(HorizontalGridSonID, HorizontalGridFatherID,       &
                                      Field3DFather, Field3DSon,                         &
                                      ComputeFather, Compute,                            &
                                      KLBFather, KUBFather, KUBSon,                      &
                                      FluxType, GridID, STAT)



        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridSonID
        integer                                     :: HorizontalGridFatherID
        real(8), dimension(:,:,:), pointer          :: Field3DFather, Field3DSon
        integer, dimension(:,:,:), pointer          :: ComputeFather
        integer,           intent(IN)               :: Compute, KLBFather, KUBFather, KUBSon
        logical, optional, intent(IN)               :: FluxType
        integer, optional, intent(IN)               :: GridID
        integer, optional, intent(OUT)              :: STAT

        !Local-----------------------------------------------------------------
        type(T_HorizontalGrid),    pointer          :: ObjHorizontalGridFather  
        type (T_FatherGrid),       pointer          :: FatherGrid
        real,    dimension(:,:  ), pointer          :: DXFather, DYFather, XXFather2D, YYFather2D
        real   , pointer, dimension(:,:)            :: XXSon, YYSon, DYXFather
        real   , pointer, dimension(:  )            :: XXFather, YYFather
        integer, pointer, dimension(:,:)            :: JSon, ISon


        real    :: YYUpper, YYLower, XXUpper, XXLower,                                   &
                   PropLowLeft,  PropUpLeft, PropLowRight, PropUpRight, PropSon
        
        integer :: ONLowLeft, ONUpLeft, ONLowRight, ONUpRight
        integer :: ready_, GridID_ 
            
        integer :: JLBSon, JUBSon, ILBSon, IUBSon
        integer :: JLBFather, JUBFather, ILBFather, IUBFather
        integer :: Jlower, Jupper, Ilower, Iupper 
        integer :: STAT_            

        integer :: i, j, k
        logical :: InterPolOK, FluxType_

        real    :: XPosition, YPosition
        integer :: Dij, Dji
        !Begin------------------------------------------------------------------

        STAT_ = UNKNOWN_


        call Ready(HorizontalGridSonID, ready_)

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then

               
            if (.NOT. associated(Field3DFather))                                         &
                call SetError(FATAL_, INTERNAL_,                                           &
                             "InterpolRegularGrid3D8; HorizontalGrid. ERR01") 

            if (.NOT. associated(Field3DSon))                                            &
                call SetError(FATAL_, INTERNAL_,                                           &
                             "InterpolRegularGrid3D8; HorizontalGrid. ERR02") 

            nullify(ObjHorizontalGridFather)
            call LocateObjFather (ObjHorizontalGridFather, HorizontalGridFatherID)


            if (present(FluxType)) then

                FluxType_ = FluxType

            else

                FluxType_ = .false.

            endif

            JLBSon = Me%WorkSize%JLB
            JUBSon = Me%WorkSize%JUB
            ILBSon = Me%WorkSize%ILB
            IUBSon = Me%WorkSize%IUB

            if (present(GridID)) then 
                GridID_ = GridID
            else
                GridID_ = LargeScaleModel_
            endif


            call Search_FatherGrid (FatherGrid, GridID_)


            !Gets the bounds from the Father
            ILBFAther = ObjHorizontalGridFather%WorkSize%ILB 
            IUBFather = ObjHorizontalGridFather%WorkSize%IUB
            JLBFather = ObjHorizontalGridFather%WorkSize%JLB
            JUBFather = ObjHorizontalGridFather%WorkSize%JUB


            if      (Compute == ComputeU_)  then
   
                JUBFather = JUBFather + 1
                JUBSon    = JUBSon    + 1


                XXSon    => FatherGrid%XX_U
                YYSon    => FatherGrid%YY_U

                JSon     => FatherGrid%JU
                ISon     => FatherGrid%IU


                XXFather => ObjHorizontalGridFather%Compute%XX_U
                YYFather => ObjHorizontalGridFather%Compute%YY_U

                DYXFather=> ObjHorizontalGridFather%DYY

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_U
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_U

                DXFather => ObjHorizontalGridFather%DXX
                DYFather => ObjHorizontalGridFather%DZY

                Dij = 0
                Dji = 1

            else if (Compute == ComputeV_)  then

                IUBFather = IUBFather + 1
                IUBSon    = IUBSon    + 1

                XXSon => FatherGrid%XX_V
                YYSon => FatherGrid%YY_V

                JSon     => FatherGrid%JV
                ISon     => FatherGrid%IV

                XXFather => ObjHorizontalGridFather%Compute%XX_V
                YYFather => ObjHorizontalGridFather%Compute%YY_V

                DYXFather=> ObjHorizontalGridFather%DXX

            else if (Compute == ComputeZ_)  then

                if (FluxType_) call SetError(FATAL_, INTERNAL_,                            &
                             "InterpolRegularGrid3D8; HorizontalGrid. ERR04") 

                XXSon => FatherGrid%XX_Z
                YYSon => FatherGrid%YY_Z

                XXFather => ObjHorizontalGridFather%Compute%XX_Z
                YYFather => ObjHorizontalGridFather%Compute%YY_Z

                JSon     => FatherGrid%JZ
                ISon     => FatherGrid%IZ

                XXFather2D => ObjHorizontalGridFather%Compute%XX2D_Z
                YYFather2D => ObjHorizontalGridFather%Compute%YY2D_Z

                DXFather => ObjHorizontalGridFather%DVX
                DYFather => ObjHorizontalGridFather%DUY

                Dji = 1
                Dij = 1

            endif


dok:        do k = KLBFather, KUBFather
doj:        do j = JLBSon, JUBSon
doi:        do i = ILBSon, IUBSon
                
                InterPolOK    = .false.

                Jlower = JSon(i, j)
                Jupper = JSon(i, j) + 1

                Ilower = ISon(i, j)
                Iupper = ISon(i, j) + 1

                if (Jupper > JUBFather) Jupper = JUBFather
                if (Iupper > IUBFather) Iupper = IUBFather
                if (Jlower < JLBFather) Jlower = JLBFather
                if (Ilower < ILBFather) Ilower = ILBFather

                if (.not. ObjHorizontalGridFather%CornersXYInput) then
                

                    XXLower  = XXFather(Jlower)
                    XXUpper  = XXFather(Jupper)

                    YYLower  = YYFather(Ilower)
                    YYUpper  = YYFather(Iupper)

                    XPosition = XXSon(i,j)
                    YPosition = YYSon(i,j)

                else

                    call CellReferential(DXFather, DYFather, XXFather2D, YYFather2D,    &
                                         ObjHorizontalGridFather%RotationX,             &
                                         ObjHorizontalGridFather%RotationY,             &
                                         XXSon(i,j), YYSon(i,j),                        &
                                         Ilower, Iupper, Jlower, Jupper, dij, dji,      &
                                         XPosition, YPosition,                          &
                                         YYUpper, YYLower, XXUpper, XXLower)

                endif


                if (FluxType_) then 

                    PropLowLeft = Field3DFather(Ilower, Jlower, k) / DYXFather(Ilower, Jlower)
                    PropUpLeft  = Field3DFather(Iupper, Jlower, k) / DYXFather(Iupper, Jlower)
                    PropLowRight= Field3DFather(Ilower, Jupper, k) / DYXFather(Ilower, Jupper)
                    PropUpRight = Field3DFather(Iupper, Jupper, k) / DYXFather(Iupper, Jupper)

                else

                    PropLowLeft = Field3DFather(Ilower, Jlower, k)
                    PropUpLeft  = Field3DFather(Iupper, Jlower, k)
                    PropLowRight= Field3DFather(Ilower, Jupper, k)
                    PropUpRight = Field3DFather(Iupper, Jupper, k)           

                endif

                ONLowLeft   = ComputeFather(Ilower, Jlower, k)
                ONUpLeft    = ComputeFather(Iupper, Jlower, k)
                ONLowRight  = ComputeFather(Ilower, Jupper, k)
                ONUpRight   = ComputeFather(Iupper, Jupper, k)


                call InterpolPoint(XPosition, YPosition,                             &
                                   YYUpper, YYLower, XXUpper, XXLower,               &
                                   PropLowLeft, PropUpLeft, PropLowRight, PropUpRight,&
                                   ONLowLeft,   ONUpLeft,   ONLowRight,   ONUpRight, &
                                   PropSon, InterPolOK)


                if (.not. InterPolOK) then
                    PropSon     = FillValueReal
                endif


                if (KUBFather > 1) then
                    Field3DSon(i, j, k) =  dble(PropSon)
                else
                    Field3DSon(i, j, KUBSon) =  dble(PropSon)
                endif


            enddo doi
            enddo doj
            enddo dok


            nullify(XXSon,    YYSon   )
            nullify(JSon,     ISon    )
            nullify(XXFather, YYFather)
            nullify(DYXFather         )

            STAT_ = SUCCESS_
        else               
            STAT_ = ready_
        end if cd1


        if (present(STAT))                                                    &
            STAT = STAT_

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

    end subroutine InterpolRegularGrid3D8



    !--------------------------------------------------------------------------
    subroutine CellReferential(DXFather, DYFather, XXFather2D, YYFather2D,    &
                               RotationX, RotationY,                          &
                               XXSon, YYSon,                                  &
                               Ilower, Iupper, Jlower, Jupper, dij, dji,      &
                               XPosition, YPosition,                          &
                               YYUpper, YYLower, XXUpper, XXLower)
        
        !Arguments-------------------------------------------------------------
        real   , pointer, dimension(:,:)  :: DXFather, DYFather,              &
                                             XXFather2D, YYFather2D,          &
                                             RotationX, RotationY
        real                              :: YYUpper, YYLower, XXUpper, XXLower, &
                                             XXSon, YYSon, XPosition, YPosition
        integer                           :: Jlower, Jupper, Ilower, Iupper, dij, dji
        !Local-----------------------------------------------------------------
        real                              :: XlocalOrigin, YlocalOrigin, dx, dy, dteta


        !Begin------------------------------------------------------------------

            XXLower  = 0.
            XXUpper  = DXFather(Ilower + dji, Jlower)

            YYLower  = 0.
            YYUpper  = DYFather(Ilower, Jlower + dij)
            !West-East
            XlocalOrigin = (XXFather2D(Ilower, Jlower) + XXFather2D(Iupper, Jlower)) / 2.
            YlocalOrigin = (YYFather2D(Ilower, Jlower) + YYFather2D(Iupper, Jlower)) / 2.
            dx           = XXSon - XlocalOrigin
            dy           = YYSon - YlocalOrigin
            dteta        = atan2(dy, dx)
            dteta        = dteta - (RotationX(Ilower, Jlower) + RotationX(Ilower + dji, Jlower + dij)) / 2.
            XPosition    = sqrt(dx**2 + dy**2)
            XPosition    = XPosition * cos(dteta)
            XPosition    = abs(XPosition)

            !South-North
            XlocalOrigin = (XXFather2D(Ilower, Jlower) + XXFather2D(Ilower, Jupper)) / 2.
            YlocalOrigin = (YYFather2D(Ilower, Jlower) + YYFather2D(Ilower, Jupper)) / 2.
            dx           = XXSon - XlocalOrigin
            dy           = YYSon - YlocalOrigin
            dteta        = atan2(dy, dx)
            dteta        = dteta - (RotationY(Ilower, Jlower) + RotationY(Ilower + dji, Jlower + dij)) / 2.
            YPosition    = sqrt(dx**2 + dy**2)
            YPosition    = YPosition * cos(dteta)
            YPosition    = abs(YPosition)

    end subroutine CellReferential

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


    subroutine WriteHorizontalGrid (HorizontalGridID, ObjHDF5, OutputNumber, STAT)


        !Arguments-------------------------------------------------------------
        integer                                     :: HorizontalGridID
        integer                                     :: ObjHDF5
        integer, optional                           :: OutputNumber
        integer, optional                           :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: ILB, IUB
        integer                                     :: JLB, JUB
        integer                                     :: WorkILB, WorkIUB
        integer                                     :: WorkJLB, WorkJUB
        integer                                     :: STAT_, ready_, STAT_CALL

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then


            !Size
            ILB = Me%Size%ILB
            IUB = Me%Size%IUB

            JLB = Me%Size%JLB
            JUB = Me%Size%JUB


            !WorkSize
            WorkILB = Me%WorkSize%ILB
            WorkIUB = Me%WorkSize%IUB

            WorkJLB = Me%WorkSize%JLB
            WorkJUB = Me%WorkSize%JUB


            !Sets limits for next write operations
            call HDF5SetLimits   (ObjHDF5, WorkILB, WorkIUB+1, WorkJLB, WorkJUB+1,           &
                                  STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR02'


            if (present(OutputNumber)) then

                call HDF5WriteData   (ObjHDF5, "/Grid/ConnectionX", "ConnectionX", "m",          &
                                      Array2D = Me%XX_IE,                                        &
                                      OutputNumber = OutputNumber,                               &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR03'

                call HDF5WriteData   (ObjHDF5, "/Grid/ConnectionY", "ConnectionY", "m",          &
                                      Array2D = Me%YY_IE,                                        &
                                      OutputNumber = OutputNumber,                               &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR04'

                call HDF5WriteData   (ObjHDF5, "/Grid/Longitude", "Longitude", "º",              &
                                      Array2D = Me%LongitudeConn,                                &
                                      OutputNumber = OutputNumber,                               &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR05'

                call HDF5WriteData   (ObjHDF5, "/Grid/Latitude", "Latitude", "º",                &
                                      Array2D = Me%LatitudeConn,                                 &
                                      OutputNumber = OutputNumber,                               &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR06'

                if (Me%CornersXYInput) then

                    call HDF5WriteData   (ObjHDF5, "/Grid/Define_Cells", "Define Cells", "-",    &
                                          Array2D = Me%DefineCellsMap,                           &
                                          OutputNumber = OutputNumber,                           &
                                          STAT = STAT_CALL)
                    if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR09'


                endif

            else

                call HDF5WriteData   (ObjHDF5, "/Grid", "ConnectionX", "m",                      &
                                      Array2D = Me%XX_IE,                                        &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR03'

                call HDF5WriteData   (ObjHDF5, "/Grid", "ConnectionY", "m",                      &
                                      Array2D = Me%YY_IE,                                        &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR04'

                call HDF5WriteData   (ObjHDF5, "/Grid", "Longitude", "º",                        &
                                      Array2D = Me%LongitudeConn,                                &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR05'

                call HDF5WriteData   (ObjHDF5, "/Grid", "Latitude", "º",                         &
                                      Array2D = Me%LatitudeConn,                                 &
                                      STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR06'

                if (Me%CornersXYInput) then

                    call HDF5WriteData   (ObjHDF5, "/Grid", "Define Cells", "-",                 &
                                          Array2D = Me%DefineCellsMap,                           &
                                          STAT = STAT_CALL)
                    if (STAT_CALL /= SUCCESS_) stop 'WriteHorizontalGrid - HorizontalGrid - ERR09'


                endif

            endif

            STAT_ = SUCCESS_

        else               

            STAT_ = ready_

        end if cd1


        if (present(STAT))  STAT = STAT_


    end subroutine WriteHorizontalGrid 

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

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


    subroutine RotateVectorFieldToGrid2D(HorizontalGridID, VectorInX, VectorInY,        &
                                         VectorOutX, VectorOutY, WaterPoints2D,         &
                                         RotateX, RotateY, STAT)

        !Arguments-------------------------------------------------------------
        real,    dimension(:, :), pointer       :: VectorInX, VectorInY, VectorOutX, VectorOutY
        integer, dimension(:, :), pointer       :: WaterPoints2D
        integer                                 :: HorizontalGridID
        logical                                 :: RotateX, RotateY
        integer, optional                       :: STAT
        !Local-----------------------------------------------------------------
        real                                    :: XGrid, YGrid, AngleX, AngleY
        integer                                 :: ILB, IUB, JLB, JUB
        integer                                 :: i, j
        integer                                 :: STAT_, ready_
        !Begin--------------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then

            
        
if1:        if  (Me%Distortion .or. Me%RegularRotation) then

                !Bounds
                ILB = Me%WorkSize%ILB 
                IUB = Me%WorkSize%IUB 

                JLB = Me%WorkSize%JLB 
                JUB = Me%WorkSize%JUB 


                !Rotate
                do j = JLB, JUB
                do i = ILB, IUB
                        
if2:                if (WaterPoints2D(i, j) == WaterPoint) then

                        !VectorX = X_Grid * cos(AngleX) + Y_Grid * cos(AngleY)
                        !VectorY = X_Grid * sin(AngleX) + Y_Grid * sin(AngleY)

                        if      (Me%Distortion) then 
                
                            AngleX = Me%RotationX(i, j)
                            AngleY = Me%RotationY(i, j)

                        else if (Me%RegularRotation) then

                            AngleX = Me%Grid_Angle
                            AngleY = Me%Grid_Angle + Pi / 2.

                        endif

                        call FromCartesianToGrid (VectorInX(i, j), VectorInY(i, j),           &
                                                  AngleX, AngleY, Xgrid, Ygrid)

                        if (RotateX) then

                            VectorOutX(i, j) = Xgrid

                        endif

                        if (RotateY) then

                            VectorOutY(i, j) = Ygrid

                        endif


                    endif if2

                enddo
                enddo

            else  if1

                if (RotateX) then

                    VectorOutX(:, :) = VectorInX(:, :)

                endif

                if (RotateY) then

                    VectorOutY(:, :) = VectorInY(:, :)

                endif

            endif if1

            STAT_ = SUCCESS_

        else               

            STAT_ = ready_

        end if cd1


        if (present(STAT))  STAT = STAT_

    end subroutine RotateVectorFieldToGrid2D
    !--------------------------------------------------------------------------

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


    subroutine RotateVectorFieldToGrid3D(HorizontalGridID, VectorInX, VectorInY,        &
                                         VectorOutX, VectorOutY, WaterPoints3D,         &
                                         RotateX, RotateY, KLB, KUB, STAT)

        !Arguments-------------------------------------------------------------
        real,    dimension(:, :, :), pointer    :: VectorInX, VectorInY, VectorOutX, VectorOutY
        integer, dimension(:, :, :), pointer    :: WaterPoints3D
        logical                                 :: RotateX, RotateY
        integer                                 :: HorizontalGridID, KLB, KUB
        integer, optional                       :: STAT
        !Local-----------------------------------------------------------------
        real                                    :: XGrid, YGrid, AngleX, AngleY
        integer                                 :: ILB, IUB, JLB, JUB
        integer                                 :: i, j, k
        integer                                 :: STAT_, ready_
        !Begin--------------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then

        
        
if1:        if  (Me%Distortion .or. Me%RegularRotation) then

                !Bounds
                ILB = Me%WorkSize%ILB 
                IUB = Me%WorkSize%IUB 

                JLB = Me%WorkSize%JLB 
                JUB = Me%WorkSize%JUB 


                !Rotate
                do k = KLB, KUB
                do j = JLB, JUB
                do i = ILB, IUB
                        
if2:                if (WaterPoints3D(i, j, k) == WaterPoint) then

                        !VectorX = X_Grid * cos(AngleX) + Y_Grid * cos(AngleY)
                        !VectorY = X_Grid * sin(AngleX) + Y_Grid * sin(AngleY)

                        if      (Me%Distortion) then 
                
                            AngleX = Me%RotationX(i, j)
                            AngleY = Me%RotationY(i, j)

                        else if (Me%RegularRotation) then

                            AngleX = Me%Grid_Angle
                            AngleY = Me%Grid_Angle + Pi / 2.

                        endif

                        call FromGridToCartesian(VectorInX(i, j, k), VectorInY(i, j, k),           &
                                                  AngleX, AngleY, Xgrid, Ygrid)

                        if (RotateX) then

                            VectorOutX(i, j, k) = Xgrid

                        endif

                        if (RotateY) then

                            VectorOutY(i, j, k) = Ygrid

                        endif


                    endif if2

                enddo
                enddo
                enddo

            else  if1


                if (RotateX) then

                    VectorOutX(:, :, :) = VectorInX(:, :, :)

                endif

                if (RotateY) then

                    VectorOutY(:, :, :) = VectorInY(:, :, :)

                endif

            endif if1

            STAT_ = SUCCESS_

        else               

            STAT_ = ready_

        end if cd1


        if (present(STAT))  STAT = STAT_

    end subroutine RotateVectorFieldToGrid3D
    !--------------------------------------------------------------------------

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


    subroutine RotateVectorGridToField2D(HorizontalGridID, VectorInX, VectorInY,        &
                                         VectorOutX, VectorOutY, WaterPoints2D,         &
                                         RotateX, RotateY, STAT)

        !Arguments-------------------------------------------------------------
        real,    dimension(:, :), pointer       :: VectorInX, VectorInY, VectorOutX, VectorOutY
        integer, dimension(:, :), pointer       :: WaterPoints2D
        integer                                 :: HorizontalGridID
        logical                                 :: RotateX, RotateY
        integer, optional                       :: STAT
        !Local-----------------------------------------------------------------
        real                                    :: XGrid, YGrid, AngleX, AngleY
        integer                                 :: ILB, IUB, JLB, JUB
        integer                                 :: i, j
        integer                                 :: STAT_, ready_
        !Begin--------------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then

            
        
if1:        if  (Me%Distortion .or. Me%RegularRotation) then

                !Bounds
                ILB = Me%WorkSize%ILB 
                IUB = Me%WorkSize%IUB 

                JLB = Me%WorkSize%JLB 
                JUB = Me%WorkSize%JUB 


                !Rotate
                do j = JLB, JUB
                do i = ILB, IUB
                        
if2:                if (WaterPoints2D(i, j) == WaterPoint) then

                        !VectorX = X_Grid * cos(AngleX) + Y_Grid * cos(AngleY)
                        !VectorY = X_Grid * sin(AngleX) + Y_Grid * sin(AngleY)

                        if      (Me%Distortion) then 
                
                            AngleX = Me%RotationX(i, j)
                            AngleY = Me%RotationY(i, j)

                        else if (Me%RegularRotation) then

                            AngleX = Me%Grid_Angle
                            AngleY = Me%Grid_Angle + Pi / 2.

                        endif

                        call FromGridToCartesian (VectorInX(i, j), VectorInY(i, j),           &
                                                  AngleX, AngleY, Xgrid, Ygrid)

                        if (RotateX) then

                            VectorOutX(i, j) = Xgrid

                        endif

                        if (RotateY) then

                            VectorOutY(i, j) = Ygrid

                        endif


                    endif if2

                enddo
                enddo

            else  if1


                if (RotateX) then

                    VectorOutX(:, :) = VectorInX(:, :)

                endif

                if (RotateY) then

                    VectorOutY(:, :) = VectorInY(:, :)

                endif


            endif if1

            STAT_ = SUCCESS_

        else               

            STAT_ = ready_

        end if cd1


        if (present(STAT))  STAT = STAT_

    end subroutine RotateVectorGridToField2D
    !--------------------------------------------------------------------------

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


    subroutine RotateVectorGridToField3D(HorizontalGridID, VectorInX, VectorInY,        &
                                         VectorOutX, VectorOutY, WaterPoints3D,         &
                                         RotateX, RotateY, KLB, KUB, STAT)

        !Arguments-------------------------------------------------------------
        real,    dimension(:, :, :), pointer    :: VectorInX, VectorInY, VectorOutX, VectorOutY
        integer, dimension(:, :, :), pointer    :: WaterPoints3D
        logical                                 :: RotateX, RotateY
        integer                                 :: HorizontalGridID, KLB, KUB
        integer, optional                       :: STAT
        !Local-----------------------------------------------------------------
        real                                    :: XGrid, YGrid, AngleX, AngleY
        integer                                 :: ILB, IUB, JLB, JUB
        integer                                 :: i, j, k
        integer                                 :: STAT_, ready_
        !Begin--------------------------------------------------------------------------

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)

cd1 :   if (ready_ == IDLE_ERR_ .or. ready_ == READ_LOCK_ERR_) then

        
        
if1:        if  (Me%Distortion .or. Me%RegularRotation) then

                !Bounds
                ILB = Me%WorkSize%ILB 
                IUB = Me%WorkSize%IUB 

                JLB = Me%WorkSize%JLB 
                JUB = Me%WorkSize%JUB 


                !Rotate
                do k = KLB, KUB
                do j = JLB, JUB
                do i = ILB, IUB
                        
if2:                if (WaterPoints3D(i, j, k) == WaterPoint) then

                        !VectorX = X_Grid * cos(AngleX) + Y_Grid * cos(AngleY)
                        !VectorY = X_Grid * sin(AngleX) + Y_Grid * sin(AngleY)

                        if      (Me%Distortion) then 
                
                            AngleX = Me%RotationX(i, j)
                            AngleY = Me%RotationY(i, j)

                        else if (Me%RegularRotation) then

                            AngleX = Me%Grid_Angle
                            AngleY = Me%Grid_Angle + Pi / 2.

                        endif

                        call FromCartesianToGrid (VectorInX(i, j, k), VectorInY(i, j, k),           &
                                                  AngleX, AngleY, Xgrid, Ygrid)

                        if (RotateX) then

                            VectorOutX(i, j, k) = Xgrid

                        endif

                        if (RotateY) then

                            VectorOutY(i, j, k) = Ygrid

                        endif


                    endif if2

                enddo
                enddo
                enddo

            else  if1

                if (RotateX) then

                    VectorOutX(:, :, :) = VectorInX(:, :, :)

                endif

                if (RotateY) then

                    VectorOutY(:, :, :) = VectorInY(:, :, :)

                endif

            endif if1

            STAT_ = SUCCESS_

        else               

            STAT_ = ready_

        end if cd1


        if (present(STAT))  STAT = STAT_

    end subroutine RotateVectorGridToField3D
    !--------------------------------------------------------------------------

    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    !DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR  

    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    subroutine KillHorizontalGrid(HorizontalGridID, STAT)

        !Arguments-------------------------------------------------------------
        integer                             :: HorizontalGridID
        integer, optional                   :: STAT

        !Local-----------------------------------------------------------------
        integer                             :: STAT_, ready_, STATUS
        integer                             :: nUsers

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

        STAT_ = UNKNOWN_

        call Ready(HorizontalGridID, ready_)    

cd1 :   if (ready_ .NE. OFF_ERR_) then

            nUsers = DeassociateInstance(mHORIZONTALGRID_,  Me%InstanceID)

            if (nUsers == 0) then

                !XX
                deallocate (Me%XX, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR00'

                !YY
                deallocate (Me%YY, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR00'

                !XX_IE
                deallocate (Me%XX_IE, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR01'

                !YY_IE
                deallocate (Me%YY_IE, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR02'

                !LatitudeConn
                deallocate (Me%LatitudeConn, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR02a'

                !LongitudeConn
                deallocate (Me%LongitudeConn, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR02b'

                !LatitudeZ
                deallocate (Me%LatitudeZ, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR03'

                !LongitudeZ
                deallocate (Me%LongitudeZ, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR04'

                !Distances
                deallocate (Me%DXX, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR05'

                deallocate (Me%DYY, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR06'

                deallocate (Me%DZX, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR07'

                deallocate (Me%DZY, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR08'

                deallocate (Me%DUX, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR09'

                deallocate (Me%DUY, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR10'

                deallocate (Me%DVX, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR11'

                deallocate (Me%DVY, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR12'

                !Coriolis
                deallocate (Me%F, stat = STATUS)
                if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR13'

                !Compute points location

                if (.not. Me%CornersXYInput) then
                    !XX_Z
                    deallocate (Me%Compute%XX_Z, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR14'

                    !YY_Z
                    deallocate (Me%Compute%YY_Z, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR15'

                    !XX_U
                    deallocate (Me%Compute%XX_U, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR16'

                    !YY_V
                    deallocate (Me%Compute%YY_V, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR17'

        
                    !XX_V
                    deallocate (Me%Compute%XX_V, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR18'

                    !YY_U
                    deallocate (Me%Compute%YY_U, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR19'


                    !XX_Cross
                    deallocate (Me%Compute%XX_Cross, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR20'

                    !YY_Cross
                    deallocate (Me%Compute%YY_Cross, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR21'

                else

                    !XX2D_Z
                    deallocate (Me%Compute%XX2D_Z, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR140'

                    !YY2D_Z
                    deallocate (Me%Compute%YY2D_Z, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR150'

                    !XX2D_U
                    deallocate (Me%Compute%XX2D_U, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR160'

                    !YY2D_V
                    deallocate (Me%Compute%YY2D_V, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR170'

        
                    !XX2D_V
                    deallocate (Me%Compute%XX2D_V, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR180'

                    !YY2D_U
                    deallocate (Me%Compute%YY2D_U, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR190'

                endif

                !DefineCellsMap
                deallocate(Me%DefineCellsMap)

                if (Me%CoordType == CIRCULAR_ .or. Me%CornersXYInput) then

                    deallocate (Me%RotationX, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR22'

                    deallocate (Me%RotationY, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR23'

                endif


                if (Me%CornersXYInput) then

                    deallocate(Me%DefineFacesUMap, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR24'

                    deallocate(Me%DefineFacesVMap, stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR25'

                    deallocate(Me%DefineCrossMap,  stat = STATUS)
                    if (STATUS /= SUCCESS_) stop 'KillHorizontalGrid - HorizontalGrid - ERR26'

                endif


                call KillFatherGridList


                !ObjHorizontalGrid
                call DeallocateInstance
                
                HorizontalGridID = 0
                STAT_            = SUCCESS_

            end if

        else 

            STAT_ = ready_

        end if cd1


        if (present(STAT)) STAT = STAT_

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

    end subroutine KillHorizontalGrid

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

    subroutine DeallocateInstance

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        type (T_HorizontalGrid), pointer          :: AuxObjHorizontalGrid
        type (T_HorizontalGrid), pointer          :: PreviousObjHorizontalGrid

        !Updates pointers
        if (Me%InstanceID == FirstHorizontalGrid%InstanceID) then
            FirstHorizontalGrid => FirstHorizontalGrid%Next
        else
            PreviousObjHorizontalGrid => FirstHorizontalGrid
            AuxObjHorizontalGrid      => FirstHorizontalGrid%Next
            do while (AuxObjHorizontalGrid%InstanceID /= Me%InstanceID)
                PreviousObjHorizontalGrid => AuxObjHorizontalGrid
                AuxObjHorizontalGrid      => AuxObjHorizontalGrid%Next
            enddo

            !Now update linked list
            PreviousObjHorizontalGrid%Next => AuxObjHorizontalGrid%Next

        endif

        !Deallocates instance
        deallocate (Me)
        nullify    (Me) 
            
    end subroutine DeallocateInstance

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

    subroutine KillFatherGridList

        !Arguments-------------------------------------------------------------

        !Local-----------------------------------------------------------------
        type (T_FatherGrid)    , pointer    :: FatherGrid
        integer                             :: status

        ! Deallocates all the water properties 

        FatherGrid => Me%FirstFatherGrid

do1 :   do while(associated(FatherGrid))  
            
            if (FatherGrid%OkZ) then


                if (.not. FatherGrid%CornersXYInput) then

                    !XX_Z
                    deallocate(FatherGrid%XX_Z, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR10") 

                    !YY_Z
                    deallocate(FatherGrid%YY_Z, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR20") 

                endif

                nullify   (FatherGrid%XX_Z)
                nullify   (FatherGrid%YY_Z)

                !IZ
                deallocate(FatherGrid%IZ, STAT = status)
                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR30") 

                nullify   (FatherGrid%IZ)

                !JZ
                deallocate(FatherGrid%JZ, STAT = status)
                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR40") 

                nullify   (FatherGrid%JZ)

            endif

            if (FatherGrid%OkU) then

                if (.not. FatherGrid%CornersXYInput) then

                    !XX_U
                    deallocate(FatherGrid%XX_U, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR50") 

                    !YY_U
                    deallocate(FatherGrid%YY_U, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR60") 

                endif

                nullify   (FatherGrid%XX_U)
                nullify   (FatherGrid%YY_U)


                !IU
                deallocate(FatherGrid%IU, STAT = status)

                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR70") 

                nullify   (FatherGrid%IU)

                !JU
                deallocate(FatherGrid%JU, STAT = status)

                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR80") 

                nullify   (FatherGrid%JU)

            endif 

            if (FatherGrid%OkV) then

                if (.not. FatherGrid%CornersXYInput) then

                    !YY_V
                    deallocate(FatherGrid%YY_V, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR90") 

                    !XX_V
                    deallocate(FatherGrid%XX_V, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR100") 
                endif

                nullify   (FatherGrid%XX_V)
                nullify   (FatherGrid%YY_V)

                !IV
                deallocate(FatherGrid%IV, STAT = status)

                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR110") 

                nullify   (FatherGrid%IV)

                !JV
                deallocate(FatherGrid%JV, STAT = status)

                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR120") 

                nullify   (FatherGrid%JV)

            endif

            if (FatherGrid%OkCross) then

                if (.not. FatherGrid%CornersXYInput) then
                
                    !XX_Cross
                    deallocate(FatherGrid%XX_Cross, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR130") 

                    !YY_Cross
                    deallocate(FatherGrid%YY_Cross, STAT = status)

                    if (status /= SUCCESS_)                                             &
                        call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR140") 

                endif
                
                nullify   (FatherGrid%XX_Cross)
                nullify   (FatherGrid%YY_Cross)



                !ICross
                deallocate(FatherGrid%ICross, STAT = status)

                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR150") 

                nullify   (FatherGrid%ICross)

                !JCross
                deallocate(FatherGrid%JCross, STAT = status)

                if (status /= SUCCESS_)                                                 &
                    call SetError(FATAL_, INTERNAL_, "KillFatherGridList; HorizontalGrid. ERR160") 

                nullify   (FatherGrid%JCross)

            endif

            FatherGrid => FatherGrid%Next
                    
        end do do1
        
        nullify   (Me%FirstFatherGrid,Me%LastFatherGrid)

    end subroutine KillFatherGridList


    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    !MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEME

    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    subroutine Ready (ObjHorizontalGrid_ID, ready_) 

        !Arguments-------------------------------------------------------------
        integer                                     :: ObjHorizontalGrid_ID
        integer                                     :: ready_

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

        nullify (Me)

cd1:    if (ObjHorizontalGrid_ID > 0) then
            call LocateObjHorizontalGrid(ObjHorizontalGrid_ID)
            ready_ = VerifyReadLock (mHORIZONTALGRID_, Me%InstanceID)
        else
            ready_ = OFF_ERR_
        end if cd1
    
     
        !----------------------------------------------------------------------

    end subroutine Ready

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

    subroutine LocateObjHorizontalGrid (ObjHorizontalGridID)

        !Arguments-------------------------------------------------------------
        integer                                     :: ObjHorizontalGridID

        !Local-----------------------------------------------------------------

        Me => FirstHorizontalGrid
        do while (associated (Me))
            if (Me%InstanceID == ObjHorizontalGridID) exit
            Me => Me%Next
        enddo

        if (.not. associated(Me)) stop 'HorizontalGrid - LocateObjHorizontalGrid - ERR01'

    end subroutine LocateObjHorizontalGrid

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

    subroutine LocateObjFather (ObjHorizontalGrid, ObjHorizontalGridID)

        !Arguments-------------------------------------------------------------
        type (T_HorizontalGrid), pointer            :: ObjHorizontalGrid
        integer                                     :: ObjHorizontalGridID

        !Local-----------------------------------------------------------------

        ObjHorizontalGrid => FirstHorizontalGrid
        do while (associated (ObjHorizontalGrid))
            if (ObjHorizontalGrid%InstanceID == ObjHorizontalGridID) exit
            ObjHorizontalGrid => ObjHorizontalGrid%Next
        enddo

        if (.not. associated(ObjHorizontalGrid)) stop 'HorizontalGrid - LocateObjFather - ERR01'

    end subroutine LocateObjFather


    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    !AUXILARY AUXILARY AUXILARY AUXILARY AUXILARY AUXILARY AUXILARY AUXILARY AU

    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    subroutine wgs84_to_rd(lat_wgs84, lon_wgs84, x_rd, y_rd)
        
        !Source code converted from C code by Ejo Schrama <e.j.o.schrama@lr.tudelft.nl>
        !The source is provided as is and has only been tested within the range of
        !the RD coordinate system and the accuracy of the conversion is no better than 50 cm


        !Arguments-------------------------------------------------------------
        real(8), intent(in)                :: lat_wgs84, lon_wgs84
        real(8), intent(out)               :: x_rd, y_rd

        !Local-----------------------------------------------------------------
        real(8)                             :: lambes, phibes

        !Begin-----------------------------------------------------------------
        

        call wgs842bessel_(lat_wgs84, lon_wgs84, phibes, lambes)

        call bessel2rd_(phibes, lambes, x_rd, y_rd)

    end subroutine


    subroutine wgs842bessel_(phiwgs, lamwgs, phibes, lambes)

        !Source code converted from C code by Ejo Schrama <e.j.o.schrama@lr.tudelft.nl>
        !The source is provided as is and has only been tested within the range of
        !the RD coordinate system and the accuracy of the conversion is no better than 50 cm

        !Arguments-------------------------------------------------------------
        real(8), intent(in)                 :: phiwgs, lamwgs
        real(8), intent(out)                :: phibes, lambes

        !Local-----------------------------------------------------------------
        real(8)                             :: dlam, dphi, lamcor, phicor

        !Begin-----------------------------------------------------------------
        
        dphi = phiwgs - 52.
        dlam = lamwgs - 5.
        phicor = (-96.862 - dphi * 11.714 - dlam * .125) * 1e-5
        lamcor = (dphi * 0.329 - 37.902 - dlam * 14.667) * 1e-5
        phibes = phiwgs - phicor
        lambes = lamwgs - lamcor

    end subroutine wgs842bessel_


    subroutine rd_to_wgs84(x_rd, y_rd, lat_wgs84, lon_wgs84)
        
        !Source code converted from C code by Ejo Schrama <e.j.o.schrama@lr.tudelft.nl>
        !The source is provided as is and has only been tested within the range of
        !the RD coordinate system and the accuracy of the conversion is no better than 50 cm

        !Arguments-------------------------------------------------------------
        real(8), intent(in)                :: x_rd, y_rd
        real(8), intent(out)               :: lat_wgs84, lon_wgs84

        !Local-----------------------------------------------------------------
        real(8)                             :: phibes, lambes

        !Begin-----------------------------------------------------------------

        call rd2bessel_(x_rd, y_rd, phibes, lambes)

        call bessel2wgs84_(phibes, lambes, lat_wgs84, lon_wgs84)

    end subroutine



    subroutine rd2bessel_(x, y, phi, lambda)
        
        !Source code converted from C code by Ejo Schrama <e.j.o.schrama@lr.tudelft.nl>
        !The source is provided as is and has only been tested within the range of
        !the RD coordinate system and the accuracy of the conversion is no better than 50 cm


        !Arguments-------------------------------------------------------------
        real(8)                             :: x, y, phi, lambda

        !Local-----------------------------------------------------------------
        logical, save                       :: first = .true.
        real(8)                             :: d__1, d__2
        real(8)                             :: bigr, cpsi, spsi
        real(8)                             :: phiprime, a, b, e
        real(8)                             :: k, m, n, q, r, w
        real(8)                             :: b0, l0, x0, y0, ca, cb, dl, sa, sb
        real(8)                             :: dq, pi, lambda0, sdl, psi, phi0
        integer                             :: i

        !Begin-----------------------------------------------------------------
        
        if (first) then
            x0 = 1.55e5
            y0 = 4.63e5
            k = .9999079
            bigr = 6382644.571
            m = .003773953832
            n = 1.00047585668
            pi = atan(1.) * 4.
            lambda0 = pi * .029931327161111111
            phi0 = pi * .28975644753333335
            l0 = pi * .029931327161111111
            b0 = pi * .28956165138333334
            e = .08169683122
            a = 6377397.155
            first = .false.
        end if

        d__1 = x - x0
        d__2 = y - y0
        r = sqrt(d__1 * d__1 + d__2 * d__2)


        if(r .ne. 0.)then
            sa = (x - x0) / r
            ca = (y - y0) / r
        else
            sa = 0.
            ca = 0.
        endif


        psi = atan2(r, k * 2. * bigr) * 2.
        cpsi = cos(psi)
        spsi = sin(psi)
        sb = ca * cos(b0) * spsi + sin(b0) * cpsi
        d__1 = sb
        cb = sqrt(1. - d__1 * d__1)
        b = acos(cb)
        sdl = sa * spsi / cb
        dl = asin(sdl)
        lambda = dl / n + lambda0
        w = log(tan(b / 2. + pi / 4.))
        q = (w - m) / n
        phiprime = atan(exp(q)) * 2. - pi / 2.


        do i = 1, 4
            
            dq = e / 2. * log((e * sin(phiprime) + 1.) / (1. - e * sin(phiprime)))
            phi = atan(exp(q + dq)) * 2. - pi / 2.
            phiprime = phi

        enddo

        lambda = lambda / pi * 180.
        phi    = phi / pi * 180.


    end subroutine 


    subroutine bessel2wgs84_(phibes, lambes, phiwgs, lamwgs)

        !Source code converted from C code by Ejo Schrama <e.j.o.schrama@lr.tudelft.nl>
        !The source is provided as is and has only been tested within the range of
        !the RD coordinate system and the accuracy of the conversion is no better than 50 cm

        
        !Arguments-------------------------------------------------------------
        real(8)                             :: phibes, lambes, phiwgs, lamwgs

        !Local-----------------------------------------------------------------
        real(8)                             :: dlam, dphi, lamcor, phicor


        !Begin-----------------------------------------------------------------

        dphi = phibes - 52.
        dlam = lambes - 5.
        phicor = (-96.862 - dphi * 11.714 - dlam * .125) * 1e-5
        lamcor = (dphi * 0.329 - 37.902 - dlam * 14.667) * 1e-5
        phiwgs = phibes + phicor
        lamwgs = lambes + lamcor

    end subroutine bessel2wgs84_


    subroutine bessel2rd_(argphi, arglam, x, y)

        !Source code converted from C code by Ejo Schrama <e.j.o.schrama@lr.tudelft.nl>
        !The source is provided as is and has only been tested within the range of
        !the RD coordinate system and the accuracy of the conversion is no better than 50 cm

        
        !Arguments-------------------------------------------------------------
        real(8)                             :: argphi, arglam, x, y

        !Local-----------------------------------------------------------------
        logical, save                       :: first = .true.
        real(8)                             :: bigr, cpsi, cpsihalf, spsi
        real(8)                             :: spsihalf, tpsihalf, a, b, e
        real(8)                             :: s2psihalf, k, m, n, q, r, w
        real(8)                             :: b0, l0, x0, y0, ca, lambda, dl, sa
        real(8)                             :: dq, pi, qprime, lambda0, phi, phi0
        real(8)                             :: d__1, d__2

        !Begin-----------------------------------------------------------------

        if (first) then
            x0 = 1.55e5
            y0 = 4.63e5
            k = .9999079
            bigr = 6382644.571
            m = .003773953832
            n = 1.00047585668
            pi = atan(1.) * 4.
            lambda0 = pi * .029931327161111111
            phi0 = pi * .28975644753333335
            l0 = pi * .029931327161111111
            b0 = pi * .28956165138333334
            e = .08169683122
            a = 6377397.155
            first = .false.
        endif


        phi = argphi / 180. * pi
        lambda = arglam / 180. * pi
        qprime = log(tan(phi / 2. + pi / 4.))
        dq = e / 2. * log((e * sin(phi) + 1.) / (1. - e * sin(phi)))
        q = qprime - dq
        w = n * q + m
        b = atan(exp(w)) * 2. - pi / 2.
        dl = n * (lambda - lambda0)

        !/* Computing 2nd power */
        d__1 = sin((b - b0) / 2.)
        !/* Computing 2nd power */
        d__2 = sin(dl / 2.)
        s2psihalf = d__1 * d__1 + d__2 * d__2 * cos(b) * cos(b0)
        cpsihalf = sqrt(1. - s2psihalf)
        spsihalf = sqrt(s2psihalf)
        tpsihalf = spsihalf / cpsihalf
        spsi = spsihalf * 2. * cpsihalf
        cpsi = 1. - s2psihalf * 2.
        sa = sin(dl) * cos(b) / spsi
        ca = (sin(b) - sin(b0) * cpsi) / (cos(b0) * spsi)
        r = k * 2. * bigr * tpsihalf
        x = r * sa + x0
        y = r * ca + y0



    end subroutine bessel2rd_



!********************************************************************** 
!                                                                     * 
!     THIS PROGRAM IS SHAREWARE AND HAS BEEN DOWNLOADED FROM INTERNET * 
!                                                                     * 
!     THIS PROGRAM CONVERTS GPS TO UNIVERSIAL TRANSVERSE MERACTOR     * 
!     COORDINATES AND VICE VERSA FOR THE NAD27 AND NAD83 DATUM.       * 
!     THIS PROGRAM WAS WRITTEN BY E. CARLSON                          * 
!     SUBROUTINES TMGRID, TCONST, TMGEOD, TCONPC,                     * 
!     WERE WRITTEN BY T. VINCENTY, NGS, IN JULY 1984 .                * 
!     THE ORGINAL PROGRAM WAS WRITTEN IN SEPTEMBER OF 1988.           * 
!                                                                     * 
!     THIS PROGRAM WAS UPDATED ON FEBUARY 16, 1989.  THE UPDATE WAS   * 
!     HAVING THE OPTION OF SAVING AND *81* RECORD FILE.               * 
!                                                                     * 
!                                                                     * 
!     THIS PROGRAM WAS UPDATED ON APRIL 3, 1990.  THE FOLLOWING UPDATE* 
!     WERE MADE:                                                      * 
!      1. CHANGE FROM JUST A CHOICE OF NAD27 OF NAD83 REFERENCE       * 
!         ELLIPSOIDS TO; CLARKE 1866, GRS80/WGS84, INTERNATIONAL, AND * 
!         ALLOW FOR USER DEFINED OTHER.                               * 
!      2. ALLOW USE OF LATITUDES IN SOUTHERN HEMISPHERE AND LONGITUDES* 
!         UP TO 360 DEGREES WEST.                                     * 
!                                                                     * 
!********************************************************************** 
!                  DISCLAIMER                                         * 
!                                                                     * 
!   THIS PROGRAM AND SUPPORTING INFORMATION IS FURNISHED BY THE       * 
! GOVERNMENT OF THE UNITED STATES OF AMERICA, AND IS ACCEPTED AND     * 
! USED BY THE RECIPIENT WITH THE UNDERSTANDING THAT THE UNITED STATES * 
! GOVERNMENT MAKES NO WARRANTIES, EXPRESS OR IMPLIED, CONCERNING THE  * 
! ACCURACY, COMPLETENESS, RELIABILITY, OR SUITABILITY OF THIS         * 
! PROGRAM, OF ITS CONSTITUENT PARTS, OR OF ANY SUPPORTING DATA.       * 
!                                                                     * 
!   THE GOVERNMENT OF THE UNITED STATES OF AMERICA SHALL BE UNDER NO  * 
! LIABILITY WHATSOEVER RESULTING FROM ANY USE OF THIS PROGRAM.  THIS  * 
! PROGRAM SHOULD NOT BE RELIED UPON AS THE SOLE BASIS FOR SOLVING A   * 
! PROBLEM WHOSE INCORRECT SOLUTION COULD RESULT IN INJURY TO PERSON   * 
! OR PROPERTY.                                                        * 
!                                                                     * 
!   THIS PROGRAM IS PROPERTY OF THE GOVERNMENT OF THE UNITED STATES   * 
! OF AMERICA.  THEREFORE, THE RECIPIENT FURTHER AGREES NOT TO ASSERT  * 
! PROPRIETARY RIGHTS THEREIN AND NOT TO REPRESENT THIS PROGRAM TO     * 
! ANYONE AS BEING OTHER THAN A GOVERNMENT PROGRAM.                    * 
!                                                                     * 
!********************************************************************** 
! TRANSFORMADO EM SUBROTINA DE CHAMADA (FLAVIO)                         
!                                                                       
!     SIGNIFICADO DAS VARIAVEIS:                                        
!                                                                       
!     (in) icode: zona UTM para inum=2                                  
!     (in)valores de inum:                                              
!       ' 1  GEODETIC POSITIONS TO UTM COORDINATES '/,                  
!       ' 2  UTM COORDINATES TO GEODETIC POSITIONS'/,                   
!       ' 4  GEODETIC POSITIONS TO PORTUGUESE COORDINATES (HIDROMOD)'/, 
!       ' 5  PORTUGUESE COORDINATES TO GEODETIC POSITIONS (HIDROMOD)'/, 
!                                                                       
!     (in) or (out)  REAIS DE DUPLA PRECISAO                            
!            ext_lat LATITUDE EM GRAUS (POSITIVO PARA NORTE)            
!            ext_long LONGITUDE EM GRAUS (POSITIVO PARA ESTE)           
!            ext_x COORDENADA X EM METROS (POSITIVO PARA ESTE)          
!            ext_y COORDENADA Y EM METROS (POSITIVO PARA NORTE)         
!***********************************************************************
      SUBROUTINE USCONVCO (inum, icode, ext_lat, ext_long, ext_x, ext_y) 
!                                                                       
!      IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        real(8)                                     :: ext_lat, ext_long
        real(8)                                     :: ext_x, ext_y
        integer                                     :: inum, icode, LD, LM, LOD, LOM
        character(1)                                :: DATNUM, NORS, EORW 
        logical                                     :: MILP 
        COMMON / MIL / MILP 
        COMMON / CONST / RAD, ER, RF, ESQ, PI 
        COMMON / DATUM / DATNUM 
        COMMON / XY / NORTH, EAST 
                                                                        
        !Local-----------------------------------------------------------------
        real(8)                                     :: PI, RAD, SLAT, SLON
        real(8)                                     :: north, east, lat, lon
        real(8)                                     :: ER, RF, F, ESQ

      PI = 4.D0 * DATAN (1.D0) 
      RAD = 180.D0 / PI 
      MILP = .FALSE. 
! passagem para as variaveis internas                                   
      IF (INUM.EQ.1.OR.INUM.EQ.4) THEN 
         IF (ext_lat.gt.0.) then 
            NORS = 'N' 
         ELSE 
            NORS = 'S' 
         ENDIF 
         IF (dabs (ext_long) .lt.180) then 
            IF (ext_long.lt.0.) then 
               EORW = 'W' 
            ELSE 
               EORW = 'E' 
            ENDIF 
         ELSEIF (dabs (ext_long) .lt.360) then 
            IF (ext_long.lt.0.) then 
               ext_long = 360. - dabs (ext_long) 
               EORW = 'E' 
            ELSE 
               ext_long = 360. - dabs (ext_long) 
               EORW = 'W' 
            ENDIF 
         ELSE 
            WRITE ( * , * ) 'erro na longitude em usconvco' 
            STOP 
         ENDIF 
         ext_lat = DABS (ext_lat) / RAD 
         ext_long = DABS (ext_long) / RAD 
         CALL TODMS (ext_lat, LD, LM, SLAT) 
         CALL TODMS (ext_long, LOD, LOM, SLON) 
      ELSEIF (INUM.EQ.2.OR.INUM.EQ.5) then 
         north = ext_y 
         east = ext_x 
      ENDIF 
!***                                                                    
                                                     !escolha do elipsoi
      IF (INUM.NE.6) CALL DATUMM (ER, RF, F, ESQ, DATNUM) 
!                                                                       
!**   USE THE NUM TO DO THE CORRECT FUNCTION                            
                                                                        
      IF (INUM.EQ.1.OR.INUM.EQ.4) THEN 
                                   !Geodetic p. -> Portuguese Mil. Coord
         IF (INUM.EQ.4) MILP = .TRUE. 
         CALL GPUT83 (LD, LM, SLAT, NORS, LOD, LOM, SLON, EORW) 
      ELSEIF (INUM.EQ.2.OR.INUM.EQ.5) THEN 
                                   !Portuguese Mil. Coordinates -> Geode
         IF (INUM.EQ.5) MILP = .TRUE. 
         CALL UTGP83 (ICODE, lat, lon) 
      ENDIF 
      IF (INUM.EQ.1.OR.INUM.EQ.4) THEN 
!        if (eorw.eq.'W') then                                          
!       ext_x=-1.*east                                                    
!     else                                                                
!                                                                       
         ext_x = east 
!     endif                                                               
!        if (nors.eq.'S') then                                          
!         ext_y=-1.*north                                               
!     else                                                                
         ext_y = north 
!     endif                                                               
      ELSE 
         ext_lat = lat * rad 
         ext_long = lon * rad 
      ENDIF 
      RETURN 
      END SUBROUTINE USCONVCO    

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


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

    SUBROUTINE UTGP83 (ICODE, lat, lon) 

        !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        integer                                     :: ICODE
        real(8)                                     :: lat, lon 
        
        !Local-----------------------------------------------------------------
        real(8)                                     :: RAD, ER, RF, ESQ, PI 
        common / CONST / RAD, ER, RF, ESQ, PI 


        CALL IUTPC (ICODE, lat, lon) 

    END SUBROUTINE UTGP83      

    !--------------------------------------------------------------------------
      
                         
    !--------------------------------------------------------------------------
    
    SUBROUTINE GPUT83 (LD, LM, SLAT, NORS, LOD, LOM, SLON, EORW) 

        !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        real(8)                                     :: SLAT, SLON
        integer                                     :: LD, LM, LOD, LOM
        character(1)                                :: NORS, EORW

        !Local-----------------------------------------------------------------
        real(8)                                     :: RAD, ER, RF, ESQ, PI                                      
        COMMON / CONST / RAD, ER, RF, ESQ, PI 


        CALL IUTGP (LD, LM, SLAT, NORS, LOD, LOM, SLON, EORW) 

    END SUBROUTINE GPUT83                         

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


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

    SUBROUTINE IUTPC (ICODE, lat, lon) 

        !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        integer                                     :: ICODE
        real(8)                                     :: lat, lon 
        
        !Local-----------------------------------------------------------------
        logical                                     :: FILFLAG, FILPRT 
        character(1)                                :: DATNUM 
        character(30)                               :: NAME 
        character(80)                               :: CARDR 
        real(8)                                     :: north, east
        logical                                     :: MILP 

        COMMON / MIL / MILP 
        COMMON / XY / NORTH, EAST 
        COMMON / DATUM / DATNUM 


        FILFLAG = .FALSE. 
        FILPRT = .FALSE. 
        NAME = '                              ' 
        WRITE (CARDR, 41) NAME 
        41 FORMAT(T15,A30) 

        ! Only to give it a number!                            
        IF (MILP)                                                             & 
            ICODE = 29 

        CALL DRUTGP (ICODE, lat, lon) 

    END SUBROUTINE IUTPC                          

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


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

      SUBROUTINE IUTGP (LD, LM, SLAT, NORS, LOD, LOM, SLON, EORW) 

      !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        real(8)                                     :: SLAT, SLON
        integer                                     :: LD, LM, LOD, LOM, ISEC, JSEC
        character(1)                                :: EORW 
        character(1)                                :: NORS 
        character(11)                               :: CLAT 
        character(12)                               :: CLON 
        character(30)                               :: NAME 
        character(80)                               :: CARDR 
        logical FILFLAG 
        logical FIL81FL 
                                                                        
      FILFLAG = .FALSE. 
      FIL81FL = .FALSE. 
      NAME = '                              ' 
                                                                        
      ISEC = SLAT * 1.0D5 + 0.5D0 
      JSEC = SLON * 1.0D5 + 0.5D0 
                                                                        
      WRITE (CLAT, 67) LD, LM, ISEC 
   67 FORMAT(I6.2,I6.2,I14.7) 
                                                                        
      WRITE (CLON, 68) LOD, LOM, JSEC 
   68 FORMAT(I6.3,I6.2,I14.7) 
                                                                        
      WRITE (CARDR, 70) NAME, CLAT, NORS, CLON, EORW 
   70 FORMAT(T7,'*80*',T15,A30,T45,A11,A1,T57,A12,A1) 
                                                                        
      CALL DRGPUT (CARDR) 
                                                                        
      RETURN 
      END SUBROUTINE IUTGP                          
!***********************************************************************
      SUBROUTINE DRGPUT (CARDR) 
!*********************************************************************  
!                                                                       
!                                                                       
!      THIS IS THE DRIVER TO COMPUTE UTM NORTHINGS AND EASTINGS         
!      FOR EACH PRIMARY ZONE AND THE AJACENT ZONE IF THE LONGITUDE      
!      IS WITH 5 MINUTES OF THE ZONE BOUNDARIES                         
!                                                                       
!      THE OUTPUT IS FOR THE DATA SHEET PROGRAM                         
!                                                                       
!      VARIABLES                                                        
!      CARDR = A MODIFIED 80 RECORD CARD WITH A LENGTH OF 211 COLS      
!      ER = EQUATORIAL RADIUS OF THE ELLIPSOID (SEMI-MAJOR AXIS)        
!      RF = RECIPROCAL OF FLATTING OF THE ELLIPSOD                      
!      ESQ= E SQUARED                                                   
!      RAD = RADIAN CONVERSION FACTOR                                   
!      CM = CENTRAL MERIDIAN ( COMPUTED USING THE LONGITUDE)            
!      SF = SCALE FACTOR OF CENTRAL MERIDIAN ( ALWAYS .9996 FOR UTM)    
!      OR = SOUTHERNMOST PARALLEL OF LATITUDE ( ALWAYS ZERO FOR UTM)    
!      R, A, B, C, U, V, W = ELLIPSOID CONSTANTS USED FOR COMPUTING     
!                            MERIDIONAL DISTANCE FROM LATITUDE          
!      SO = MERIDIONAL DISTANCE (MULTIPLIED BY SCALE FACTOR )           
!           FROM THE EQUATOR TO THE SOUTHERNMOST PARALLEL OF LATITUDE   
!           ( ALWAYS ZERO FOR UTM)                                      
!                                                                       
                                                                        
      !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        character(80)                               :: CARDR 

        !Local-----------------------------------------------------------------
        character(1)                                :: EORW 
        character(1)                                :: NORS 
        character(4)                                :: ZONE 
        integer                                     :: LD, LM, LOD, LOM, FOUND, IZ, ICM
        real(8)                                     :: SLAT, SLON, FI, LAM, LCM, UCM 
        real(8)                                     :: KP 
        real(8)                                     :: TD, TM, ND, NM 
        real(8)                                     :: north, east
        real(8)                                     :: LOD1, R, RAD, CM, TOL, FN, FE, SF, OR
        real(8)                                     :: ER, RF, ESQ, EPS, A, B, C, U, V, W, SO
        real(8)                                     :: CONV, PI
        logical                                     :: MILP 
        COMMON / MIL / MILP 
        COMMON / CONST / RAD, ER, RF, ESQ, PI 
        COMMON / XY / NORTH, EAST 
                                                                        
      READ (CARDR, 50) LD, LM, SLAT, NORS, LOD, LOM, SLON, EORW 
   50 FORMAT (T45,I2,I2,F10.5,A1,I3,I2,F10.5,A1) 
                                                                        
                                                                        
!                                                                       
!      CONVERT THE LATITUDE AND LONGITUDE TO PI AND LAM                 
!                                                                       
      TD = DBLE (FLOAT (LD) ) 
      TM = DBLE (FLOAT (LM) ) 
      FI = (TD+ (TM + SLAT / 60.D0) / 60.D0) / RAD 
                                                                        
      ND = DBLE (FLOAT (LOD) ) 
      NM = DBLE (FLOAT (LOM) ) 
                                                                        
      IF ( (EORW.EQ.'E') .OR. (EORW.EQ.'e') ) THEN 
         LAM = (360.D0 - (ND+ (NM + SLON / 60.D0) / 60.D0) ) / RAD 
         LOD1 = (360.D0 - (ND+ (NM + SLON / 60.D0) / 60.D0) ) 
         LOD = DINT (LOD1) 
      ENDIF 
                                                                        
      IF ( (EORW.EQ.'W') .OR. (EORW.EQ.'w') ) THEN 
         LAM = (ND+ (NM + SLON / 60.D0) / 60.D0) / RAD 
         LOD = LOD 
      ENDIF 
                                                                        
                                                                        
                                                                        
!                                                                       
!      FIND THE ZONE FOR LONGITUDE LESS THAN 180 DEGREES                
!                                                                       
      IF (LOD.LT.180) THEN 
         IZ = LOD / 6 
         IZ = 30 - IZ 
                                                                        
         ICM = (183 - (6 * IZ) ) 
         CM = DBLE (FLOAT (ICM) ) / RAD 
         UCM = (ICM + 3) / RAD 
         LCM = (ICM - 3) / RAD 
      ENDIF 
!                                                                       
!       FIND THE ZONE FOR LONGITUDE GREATER THAN 180 DEGREES            
!                                                                       
      IF (LOD.GE.180) THEN 
         IZ = (LOD) / 6 
         IZ = 90 - IZ 
         ICM = (543 - (6 * IZ) ) 
         CM = DBLE (FLOAT (ICM) ) / RAD 
         UCM = (ICM + 3) / RAD 
         LCM = (ICM - 3) / RAD 
      ENDIF 
                                                                        
                                                                        
                                                                        
      TOL = (5.0D0 / 60.0D0) / RAD 
                                                                        
      FN = 0.D0 
                                                                        
      IF ( (NORS.EQ.'S') .OR. (NORS.EQ.'s') ) THEN 
         FN = 10000000.D0 
      ENDIF 
                                                                        
      IF ( (NORS.EQ.'N') .OR. (NORS.EQ.'n') ) THEN 
         FN = 0.D0 
      ENDIF 
                                                                        
                                                                        
      FE = 500000.0D0 
      SF = 0.9996D0 
      OR = 0.0D0 
                                                                        
      FOUND = 0 
                                                                        
      CALL TCONST (ER, RF, SF, OR, ESQ, EPS, R, A, B, C, U, V, W, SO) 
                                                                        
!                                                                       
!      COMPUTE THE NORTH AND EASTINGS                                   
!                                                                       
                                          ! Long. of origin of Portugues
      IF (MILP) CM = DBLE (8.13190611) / RAD 
                                          ! coordinates (8d 7' 54.862") 
!                                                                       
  200 CALL TMGRID (FI, LAM, NORTH, EAST, CONV, KP, ER, ESQ, EPS, CM, FE,&
      FN, SF, SO, R, A, B, C)                                           
!                                                                       
! --> Computes the transformations from UTM's to Portuguese Mil. coordin
!                                                                       
      IF (MILP) THEN 
         NORTH = NORTH / SF - 4092593.43D0 
         EAST = EAST / SF - 300200.D0 
      ENDIF 
!                                                                       
!      WRITE THE ZONE NUMBER                                            
!                                                                       
      IF (IZ.GT.9) THEN 
         WRITE (ZONE, 600) IZ 
  600 FORMAT   (1X,I2) 
      ELSE 
         WRITE (ZONE, 605) IZ 
  605 FORMAT   (1X,I6.2) 
      ENDIF 
                                                                        
!                                                                       
!      DO THE TEST TO SEE IF THE LONGITUDE IS WITHIN 5 MINUTES          
!      OF THE BOUNDARIES FOR THE ZONE AND IF SO COMPUTE THE             
!      NORHT AND EASTING FOR THE ADJACENT ZONE                          
!                                                                       
      IF (FOUND.NE.0) THEN 
         RETURN 
      ENDIF 
                                                                        
      IF (DABS (UCM - LAM) .LE.TOL) THEN 
         CM = DBLE (FLOAT (ICM + 6) ) / RAD 
         IZ = IZ - 1 
         IF (IZ.EQ.0) IZ = 60 
         FOUND = FOUND+1 
         GOTO 200 
      ENDIF 
                                                                        
      IF (DABS (LCM - LAM) .LE.TOL) THEN 
         CM = DBLE (FLOAT (ICM - 6) ) / RAD 
         IZ = IZ + 1 
         IF (IZ.EQ.61) IZ = 1 
         FOUND = FOUND+1 
         GOTO 200 
      ENDIF 
                                                                        
                                                                        
      RETURN 
      END SUBROUTINE DRGPUT                         
                                                                        
!******************************************************************     
      SUBROUTINE DRUTGP (ICODE, lat, lon) 
!                                                                       
!                                                                       
!                                                                       
!      THIS IS THE DRIVER TO COMPUTE LATITUDES AND LONGITUDES FROM      
!      THE UTMs FOR EACH ZONE                                           
!                                                                       
!                                                                       
!      VARIABLES                                                        
!      CARDR = A MODIFIED 80 RECORD CARD WITH A LENGTH OF 211 COLS      
!      ER = EQUATORIAL RADIUS OF THE ELLIPSOID (SEMI-MAJOR AXIS)        
!      RF = RECIPROCAL OF FLATTING OF THE ELLIPSOD                      
!      ESQ= E SQUARED                                                   
!      RAD = RADIAN CONVERSION FACTOR                                   
!      CM = CENTRAL MERIDIAN ( COMPUTED USEING THE LONGITUDE)           
!      SF = SCALE FACTOR OF CENTRAL MERIDIAN ( ALWAYS .9996 FOR UTM)    
!      OR = SOUTHERNMOST PARALLEL OF LATITUDE ( ALWAYS ZERO FOR UTM)    
!      R, A, B, C, U, V, W = ELLIPSOID CONSTANTS USED FOR COMPUTING     
!                            MERIDIONAL DISTANCE FROM LATITUDE          
!      SO = MERIDIONAL DISTANCE (MULTIPLIED BY SCALE FACTOR )           
!           FROM THE EQUATOR TO THE SOUTHERNMOST PARALLEL OF LATITUDE   
!           ( ALWAYS ZERO FOR UTM)                                      
!                                                                       
                                                                        
      !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        integer                                     :: ICODE
        real(8)                                     :: lat, lon 
        
        !Local-----------------------------------------------------------------
        integer                                     :: IZ, ICM
        character(1)                                :: NORS 
        character(4)                                :: ZONE 
        integer                                     :: FOUND 
        real(8)                                     :: LCM, UCM 
        real(8)                                     :: KP 
        real(8)                                     :: NORTH, EAST, R 
        character(1)                                :: EORW 
        logical                                     :: MILP
        real(8)                                     :: CM, RAD, TOL, FN, FE, SF, OR, CONV
        real(8)                                     :: EPS, SO, V0, V2, V4, V6, ER, ESQ, RF 
        real(8)                                     :: PI
        COMMON / MIL / MILP 
        COMMON / CONST / RAD, ER, RF, ESQ, PI 
        COMMON / XY / NORTH, EAST 


!                                                                       
!      FIND THE CENTRAL MERIDAIN IF THE ZONE NUMBER IS LESS THAN 30     
!                                                                       
      IF (ICODE.LT.30) THEN 
         IZ = ICODE 
         ICM = (183 - (6 * IZ) ) 
         CM = DBLE (FLOAT (ICM) ) / RAD 
         UCM = (ICM + 3) / RAD 
         LCM = (ICM - 3) / RAD 
      ENDIF 
!                                                                       
!       FIND THE CENTRAL MERIDAN IF THE ZONE NUMBER IS LARGER THAN 30   
!                                                                       
      IF (ICODE.GE.30) THEN 
         IZ = ICODE 
         ICM = (543 - (6 * IZ) ) 
         CM = DBLE (FLOAT (ICM) ) / RAD 
         UCM = (ICM + 3) / RAD 
         LCM = (ICM - 3) / RAD 
      ENDIF 
                                                                        
      TOL = (5.0D0 / 60.0D0) / RAD 
                                                                        
      IF (NORTH.GT.10000000.0) THEN 
         FN = 10000000.0D0 
         NORS = 'S' 
      ELSE 
         FN = 0.D0 
         NORS = 'N' 
      ENDIF 
                                                                        
      FE = 500000.0D0 
      SF = 0.9996D0 
      OR = 0.0D0 
                                                                        
      FOUND = 0 
                                                                        
      CALL TCONPC (SF, OR, EPS, R, SO, V0, V2, V4, V6, ER, ESQ, RF) 
                                                                        
!                                                                       
!      COMPUTE THE LATITUDES AND LONGITUDES                             
!                                                                       
                                          ! Long. of origin of Portugues
      IF (MILP) CM = DBLE (8.13190611) / RAD 
                                          ! coordinates (8d 7' 54.862") 
!                                                                       
!                                                                       
! --> Computes the transformations from UTM's to Portuguese Mil. coordin
!                                                                       
      IF (MILP) THEN 
         NORTH = (NORTH + 4092593.43D0) * SF 
         EAST = (EAST + 300200.D0) * SF 
      ENDIF 
  200 CALL TMGEOD (NORTH, EAST, LAT, LON, EPS, CM, FE, SF, SO, R, V0,   &
      V2, V4, V6, FN, ER, ESQ, CONV, KP)                                
                                                                        
!                                                                       
!      WRITE THE ZONE NUMBER                                            
!                                                                       
      IF (IZ.GT.9) THEN 
         WRITE (ZONE, 600) IZ 
  600 FORMAT   (1X,I2) 
      ELSE 
         WRITE (ZONE, 605) IZ 
  605 FORMAT   (1X,I6.2) 
      ENDIF 
                                                                        
      EORW = 'W' 
      IF (nors.eq.'s'.or.nors.eq.'S') lat = - 1.d0 * lat 
      IF (eorw.eq.'w'.or.eorw.eq.'W') lon = - 1.d0 * lon 
!                                                                       
!      DO THE TEST TO SEE IF THE LONGITUDE IS WITHIN 5 MINUTES          
!      OF THE BOUNDARIES FOR THE ZONE AND IF SO COMPUTE THE             
!      NORHT AND EASTING FOR THE ADJACENT ZONE                          
!                                                                       
      IF (FOUND.NE.0) THEN 
         RETURN 
      ENDIF 
                                                                        
      IF (DABS (UCM) .LE.TOL) THEN 
         CM = DBLE (FLOAT (ICM + 6) ) / RAD 
         IZ = IZ - 1 
         IF (IZ.EQ.0) IZ = 60 
         FOUND = FOUND+1 
         GOTO 200 
      ENDIF 
                                                                        
      IF (DABS (LCM) .LE.TOL) THEN 
         CM = DBLE (FLOAT (ICM - 6) ) / RAD 
         IZ = IZ + 1 
         IF (IZ.EQ.61) IZ = 1 
         FOUND = FOUND+1 
         GOTO 200 
      ENDIF 
                                                                        
                                                                        
      RETURN 
      END SUBROUTINE DRUTGP                         
                                                                        
!********************************************************************   
      SUBROUTINE TODMS (RAD, IDG, MIN, SEC) 
!********************************************************************   
!     RADIANS TO DEGREES,MINUTES AND SECONDS                            
!                                                                       
      INTEGER IDG, MIN
      DOUBLEPRECISION RAD, SEC, RHOSEC 
      DATA RHOSEC / 2.062648062471D05 / 
      SEC = RAD * RHOSEC 
      IDG = SEC / 3600.D0 
      SEC = SEC - DBLE (IDG * 3600) 
      MIN = SEC / 60.D0 
      SEC = SEC - DBLE (MIN * 60) 
      IF ( (60.D0 - DABS (SEC) ) .GT.5.D-6) GOTO 100 
      SEC = SEC - DSIGN (60.D0, SEC) 
      MIN = MIN + ISIGN (1, MIN) 
  100 IF (IABS (MIN) .LT.60) GOTO 101 
      MIN = MIN - ISIGN (60, MIN) 
      IDG = IDG + ISIGN (1, IDG) 
  101 MIN = IABS (MIN) 
      SEC = DABS (SEC) 
      IF (RAD.GE.0.D0) GOTO 102 
      IF (IDG.EQ.0) MIN = - MIN 
      IF (IDG.EQ.0.AND.MIN.EQ.0) SEC = - SEC 
  102 RETURN 
      END SUBROUTINE TODMS                          
!***********************************************************            
      SUBROUTINE DATUMM (ER, RF, F, ESQ, DATNUM) 
                                                                        
      CHARACTER(1) ANS, DATNUM 
      DOUBLEPRECISION ER 
      DOUBLEPRECISION RF 
      DOUBLEPRECISION F 
      DOUBLEPRECISION ESQ 
                                                                        
                                                                        
                                                                        
!   50 PRINT *, '                                          '            
!      PRINT *, '   WHICH ELLIPSOID DO YOU WANT ANSWER:    '            
!      PRINT *, '   1.  CLARKE 1866                       '             
!      PRINT *, '   2.  GRS80/WGS84                       '             
!      PRINT *, '   3.  INTERNATIONAL 1910                '             
!      PRINT *, '   4.  WGS72                             '             
!      PRINT *, '   5.  OTHER ELLIPSOID                   '             
!      PRINT *, '   TYPE NUMBER:  '                                     
!      READ(*,FMT='(A1)') ANS                                           
                !faz sempre para internacional                          
   50 ans = '3' 
!*                                                                      
!*     FIND THE RIGHT SEMI MAJOR AXIS AND FLATTING                      
!*                                                                      
      IF (ANS.EQ.'1') THEN 
!*      FOR THE NAD 27 DATUM                                            
!*                                                                      
         ER = 6378206.4D0 
         RF = 294.978698D0 
         F = 1.D0 / RF 
         ESQ = (F + F - F * F) 
                                                                        
      ELSEIF (ANS.EQ.'2') THEN 
!*     FOR THE NAD83 DATUM                                              
!*                                                                      
         ER = 6378137.D0 
         RF = 298.257222101D0 
         F = 1.D0 / RF 
         ESQ = (F + F - F * F) 
                                                                        
      ELSEIF (ANS.EQ.'3') THEN 
!*     FOR THE INT24 DATUM                                              
!*                                                                      
         ER = 6378388.D0 
         RF = 297.0D0 
         F = 1.D0 / RF 
         ESQ = (F + F - F * F) 
                                                                        
      ELSEIF (ANS.EQ.'4') THEN 
!*     FOR THE WGS72                                                    
!*                                                                      
         ER = 6378135.D0 
         RF = 298.26D0 
         F = 1.D0 / RF 
         ESQ = (F + F - F * F) 
                                                                        
      ELSEIF (ANS.EQ.'5') THEN 
!*     FOR THE OTHER DATUM                                              
!*                                                                      
   10 PRINT * , '                    ' 
      PRINT * , '  SEMIMAJOR AXIS (meters)  ' 
      PRINT * , '  TYPE VALUE NOW:  ' 
         READ ( * , FMT = '(F12.0)') ER 
                                                                        
         IF ( (ER.LE.6376400.D0) .OR. (ER.GT.6378500.D0) ) THEN 
      PRINT * , '                  ' 
            PRINT * , ' SEMIMAJOR AXIS IS OUT OF RANGE - DO YOU' 
            PRINT * , ' WANT TO TRY AGAIN ' 
            PRINT * , ' TYPE Y OR N ' 
            READ ( * , FMT = '(A1)') ANS 
            IF ( (ANS.EQ.'Y') .OR. (ANS.EQ.'y') ) THEN 
               GOTO 10 
            ELSE 
               GOTO 50 
            ENDIF 
         ENDIF 
                                                                        
   20 PRINT * , '                  ' 
      PRINT * , ' FLATTENING            ' 
      PRINT * , ' TYPE VALUE NOW:  ' 
         READ ( * , FMT = '(F11.0)') RF 
                                                                        
         IF ( (RF.LE.290.D0) .OR. (RF.GT.302.D0) ) THEN 
      PRINT * , '                  ' 
            PRINT * , ' FLATTENING IS OUT OF RANGE - DO YOU ' 
            PRINT * , ' WANT TO TRY AGAIN ' 
            PRINT * , ' TYPE Y OR N ' 
            READ ( * , FMT = '(A1)') ANS 
            IF ( (ANS.EQ.'Y') .OR. (ANS.EQ.'y') ) THEN 
               GOTO 20 
            ELSE 
               GOTO 50 
            ENDIF 
         ENDIF 
                                                                        
         F = 1.D0 / RF 
         ESQ = (F + F - F * F) 
                                                                        
      ELSE 
      PRINT * , ' YOU TYPED THE INCORRECT NUMBER   ' 
         PRINT * , ' SO LET''S TRY AGAIN ' 
      PRINT * , '                  ' 
         GOTO 50 
      ENDIF 
                                                                        
      DATNUM = ANS 
                                                                        
      RETURN 
      END SUBROUTINE DATUMM                         

      SUBROUTINE TMGRID (FI, LAM, NORTH, EAST, CONV, KP, ER, ESQ, EPS,  &
                         CM, FE, FN, SF, SO, R, A, B, C)                                   
                                                                        
        !IMPLICIT DOUBLEPRECISION (A - H, K - Z) 
        !Arguments-------------------------------------------------------------
        real(8)                                     :: FI, LAM, NORTH, EAST
        real(8)                                     :: CONV, KP, ER, ESQ, EPS
        real(8)                                     :: CM, FE, FN, SF, SO, R
        real(8)                                     :: A, B, C

        !Local-----------------------------------------------------------------
        REAL(8)                                     :: A2, A1, A4, A6, A3, A5, A7
        REAL(8)                                     :: C1, C3, C5, F2, F4
        real(8)                                     :: OM, S, SINFI, COSFI, TN, TS
        real(8)                                     :: ETS, L, LS, RN
!                                                                       
!*****  TRANSVERSE MERCATOR PROJECTION                                  
!       CONVERSION OF GEODETIC COORDINATES TO GRID COORDINATES          
!*****  Programmed by T. Vincenty, NGS, in July 1984.                   
!*****************  SYMBOLS AND DEFINITIONS *************************   
!   Latitude positive north, longitude positive west.  All angles are   
!     in radian measure.                                                
!   N, E are northing and easting coordinates respectively.             
!   LAT, LON are latitude and longitude respectively.                   
!   CONV is convergence.                                                
!   KP is point scale factor.                                           
!   ER is equatorial radius of the ellipsoid (= major semiaxis).        
!   ESQ is the square of first eccentricity of the ellipsoid.           
!   EPS is the square of second eccentricity of the ellipsoid.          
!   CM is the central meridian of the projection zone.                  
!   FE is false easting value at the central meridian.                  
!   FN is "false northing" at the southernmost latitude, usually zero.  
!   SF is scale factor at the central meridian.                         
!   SO is meridional distance (multiplied by the scale factor) from     
!     the equator to the southernmost parallel of latitude for the zone.
!   R is the radius of the rectifying sphere (used for computing        
!     meridional distance from latitude and vice versa).                
!   A, B, C, U, V, W are other precomputed constants for determination  
!     of meridional distance from latitude and vice versa.              
!                                                                       
!   The formula used in this subroutine gives geodetic accuracy within  
!   zones of 7 degrees in east-west extent.  Within State transverse    
!   Mercator projection zones, several minor terms of the equations     
!   may be omitted (see a separate NGS publication).  If programmed     
!   in full, the subroutine can be used for computations in surveys     
!   extending over two zones.                                           
!                                                                       
!*********************************************************************  
      OM = FI + A * SIN (2. * FI) + B * SIN (4. * FI) + C * SIN (6. *   &
      FI)                                                               
      S = R * OM * SF 
      SINFI = SIN (FI) 
      COSFI = COS (FI) 
      TN = SINFI / COSFI 
      TS = TN**2 
      ETS = EPS * COSFI**2 
      L = (LAM - CM) * COSFI 
      LS = L * L 
      RN = SF * ER / SQRT (1. - ESQ * SINFI**2) 
!                                                                       
      A2 = RN * TN / 2. 
      A4 = (5. - TS + ETS * (9. + 4. * ETS) ) / 12. 
      A6 = (61. + TS * (TS - 58.) + ETS * (270. - 330. * TS) ) / 360. 
      A1 = - RN 
      A3 = (1. - TS + ETS) / 6. 
      A5 = (5. + TS * (TS - 18.) + ETS * (14. - 58. * TS) ) / 120. 
      A7 = (61. - 479. * TS + 179. * TS**2 - TS**3) / 5040. 
      NORTH = S - SO + A2 * LS * (1. + LS * (A4 + A6 * LS) ) - FN 
      EAST = FE+A1 * L * (1. + LS * (A3 + LS * (A5 + A7 * LS) ) ) 
                                                                        
      IF (NORTH.LT.0.0) THEN 
         NORTH = NORTH * ( - 1.0D0) 
      ENDIF 
!                                                                       
!*** CONVERGENCE                                                        
      C1 = - TN 
      C3 = (1. + 3. * ETS + 2. * ETS**2) / 3. 
      C5 = (2. - TS) / 15. 
      CONV = C1 * L * (1. + LS * (C3 + C5 * LS) ) 
!                                                                       
!*** POINT SCALE FACTOR                                                 
      F2 = (1. + ETS) / 2. 
      F4 = (5. - 4. * TS + ETS * (9. - 24. * TS) ) / 12. 
      KP = SF * (1. + F2 * LS * (1. + F4 * LS) ) 
!                                                                       
      RETURN 
      END SUBROUTINE TMGRID                         
!********************************************************************   
      SUBROUTINE TCONST (ER, RF, SF, OR, ESQ, EPS, R, A, B, C, U, V, W, &
                         SO)                                                               
        !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        real(8)                                     :: ER, RF, SF, OR, ESQ, EPS
        real(8)                                     :: R, A, B, C, U, V, W, SO

        !Local-----------------------------------------------------------------
        real(8)                                     :: F, PR, EN, OMO
!                                                                       
!***** TRANSVERSE MERCATOR PROJECTION                                   
!      PRECOMPUTATION OF CONSTANTS                                      
!***** Programmed by T. Vincenty, NGS, in July 1984.                    
!******************** SYMBOLS AND DEFINITIONS  **********************   
!   ER is equatorial radius of the ellipsoid (= major semiaxis).        
!   RF is reciprocal of flattening of the ellipsoid.                    
!   SF is scale factor of the central meridian.                         
!   OR is southernmost parallel of latitude (in radians) for which      
!     the northing coordinate is zero at the central meridian.          
!   R, A, B, C, U, V, W are ellipsoid constants used for computing      
!     meridional distance from latitude and vice versa.                 
!   SO is meridional distance (multiplied by the scale factor) from     
!     the equator to the southernmost parallel of latitude.             
!******************************************************************     
!                                                                       
      F = 1. / RF 
      ESQ = (F + F - F**2) 
      EPS = ESQ / (1. - ESQ) 
      PR = (1. - F) * ER 
      EN = (ER - PR) / (ER + PR) 
      A = - 1.5D0 * EN + (9. / 16.) * EN**3 
      B = 0.9375D0 * EN**2 - (15. / 32.) * EN**4 
      C = - (35. / 48.) * EN**3 
      U = 1.5D0 * EN - (27. / 32.) * EN**3 
      V = 1.3125D0 * EN**2 - (55. / 32.) * EN**4 
      W = (151. / 96.) * EN**3 
      R = ER * (1. - EN) * (1. - EN**2) * (1. + 2.25D0 * EN**2 +        &
      (225. / 64.) * EN**4)                                             
      OMO = OR + A * SIN (2. * OR) + B * SIN (4. * OR) + C * SIN (6. *  &
      OR)                                                               
      SO = SF * R * OMO 
!                                                                       
      RETURN 
      END SUBROUTINE TCONST                         

      SUBROUTINE TCONPC (SF, OR, EPS, R, SO, V0, V2, V4, V6, ER, ESQ,   &
                         RF)                                                               

        !IMPLICIT DOUBLEPRECISION (A - H, O - Z) 
        !Arguments-------------------------------------------------------------
        real(8)                                     :: SF, OR, EPS, R, SO, RF
        real(8)                                     :: V0, V2, V4, V6, ER, ESQ
        
        !Local-----------------------------------------------------------------
        real(8)                                     :: F, PR, EN, EN2, EN3, EN4
        real(8)                                     :: C2, C4, C6, C8, U0, U2, U4, U6
        real(8)                                     :: COSOR, OMO
                                                                        
!**          TRANSVERSE MERCATOR PROJECTION               ***           
!** CONVERSION OF GRID COORDS TO GEODETIC COORDS                        
!** REVISED SUBROUTINE OF T. VINCENTY  FEB. 25, 1985                    
!************* SYMBOLS AND DEFINITIONS ***********************          
!** ER IS THE SEMI-MAJOR AXIS FOR GRS-80                                
!** SF IS THE SCALE FACTOR AT THE CM                                    
!** SO IS THE MERIDIONAL DISTANCE (TIMES THE SF) FROM THE               
!**       EQUATOR TO SOUTHERNMOST PARALLEL OF LAT. FOR THE ZONE         
!** R IS THE RADIUS OF THE RECTIFYING SPHERE                            
!** U0,U2,U4,U6,V0,V2,V4,V6 ARE PRECOMPUTED CONSTANTS FOR               
!**   DETERMINATION OF MERIDIONAL DIST. FROM LATITUDE                   
!** OR IS THE SOUTHERNMOST PARALLEL OF LATITUDE FOR WHICH THE           
!**       NORTHING COORD IS ZERO AT THE CM                              
!*************************************************************          
                                                                        
                                                                        
      F = 1.D0 / RF 
      EPS = ESQ / (1.D0 - ESQ) 
      PR = (1.D0 - F) * ER 
      EN = (ER - PR) / (ER + PR) 
      EN2 = EN * EN 
      EN3 = EN * EN * EN 
      EN4 = EN2 * EN2 
                                                                        
      C2 = - 3.D0 * EN / 2.D0 + 9.D0 * EN3 / 16.D0 
      C4 = 15.D0 * EN2 / 16.D0 - 15.D0 * EN4 / 32.D0 
      C6 = - 35.D0 * EN3 / 48.D0 
      C8 = 315.D0 * EN4 / 512.D0 
      U0 = 2.D0 * (C2 - 2.D0 * C4 + 3.D0 * C6 - 4.D0 * C8) 
      U2 = 8.D0 * (C4 - 4.D0 * C6 + 10.D0 * C8) 
      U4 = 32.D0 * (C6 - 6.D0 * C8) 
      U6 = 128.D0 * C8 
                                                                        
      C2 = 3.D0 * EN / 2.D0 - 27.D0 * EN3 / 32.D0 
      C4 = 21.D0 * EN2 / 16.D0 - 55.D0 * EN4 / 32.D0 
      C6 = 151.D0 * EN3 / 96.D0 
      C8 = 1097.D0 * EN4 / 512.D0 
      V0 = 2.D0 * (C2 - 2.D0 * C4 + 3.D0 * C6 - 4.D0 * C8) 
      V2 = 8.D0 * (C4 - 4.D0 * C6 + 10.D0 * C8) 
      V4 = 32.D0 * (C6 - 6.D0 * C8) 
      V6 = 128.D0 * C8 
                                                                        
      R = ER * (1.D0 - EN) * (1.D0 - EN * EN) * (1.D0 + 2.25D0 * EN *   &
      EN + (225.D0 / 64.D0) * EN4)                                      
      COSOR = DCOS (OR) 
      OMO = OR + DSIN (OR) * COSOR * (U0 + U2 * COSOR * COSOR + U4 *    &
      COSOR**4 + U6 * COSOR**6)                                         
      SO = SF * R * OMO 
                                                                        
      RETURN 
      END SUBROUTINE TCONPC                         

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

      SUBROUTINE TMGEOD (N, E, LAT, LON, EPS, CM, FE, SF, SO, R, V0, V2,&
                         V4, V6, FN, ER, ESQ, CONV, KP)                                    
        
        !IMPLICIT DOUBLEPRECISION (A - H, K - Z) 
        !Arguments-------------------------------------------------------------
        real(8)                                     :: N, E
        real(8)                                     :: LAT, LON, EPS, CM, FE
        real(8)                                     :: SF, SO, R, V0, V2, V4, V6
        real(8)                                     :: FN, ER, ESQ, CONV, KP

        !Local----------------------------------------------------------------
        real(8)                                     :: B1, B2, C1, TN, TS, ETS, RN, Q
        real(8)                                     :: OM, COSOM, FOOT, SINF, COSF 
        real(8)                                     :: QS, B4, B6, B3, B5, B7, L, FI, LAM
        real(8)                                     :: SINFI, COSFI, L1, LS, C3, C5, F2, F4
                                                                        
!**          TRANSVERSE MERCATOR PROJECTION               ***           
!** CONVERSION OF GRID COORDS TO GEODETIC COORDS                        
!** REVISED SUBROUTINE OF T. VINCENTY  FEB. 25, 1985                    
!************* SYMBOLS AND DEFINITIONS ***********************          
!** LATITUDE POSITIVE NORTH, LONGITUDE POSITIVE WEST.  ALL              
!**          ANGLES ARE IN RADIAN MEASURE.                              
!** LAT,LON ARE LAT. AND LONG. RESPECTIVELY                             
!** N,E ARE NORTHING AND EASTING COORDINATES RESPECTIVELY               
!** K IS POINT SCALE FACTOR                                             
!** ER IS THE SEMI-MAJOR AXIS OF THE ELLIPSOID                          
!** ESQ IS THE SQUARE OF THE 1ST ECCENTRICITY                           
!** E IS THE 1ST ECCENTRICITY                                           
!** CM IS THE CENTRAL MERIDIAN OF THE PROJECTION ZONE                   
!** FE IS THE FALSE EASTING VALUE AT THE CM                             
!** CONV IS CONVERGENCE                                                 
!** EPS IS THE SQUARE OF THE 2ND ECCENTRICITY                           
!** SF IS THE SCALE FACTOR AT THE CM                                    
!** SO IS THE MERIDIONAL DISTANCE (TIMES THE SF) FROM THE               
!**       EQUATOR TO SOUTHERNMOST PARALLEL OF LAT. FOR THE ZONE         
!** R IS THE RADIUS OF THE RECTIFYING SPHERE                            
!** U0,U2,U4,U6,V0,V2,V4,V6 ARE PRECOMPUTED CONSTANTS FOR               
!**   DETERMINATION OF MERIDIANAL DIST. FROM LATITUDE                   
!**                                                                     
!** THE FORMULA USED IN THIS SUBROUTINE GIVES GEODETIC ACCURACY         
!** WITHIN ZONES OF 7 DEGREES IN EAST-WEST EXTENT.  WITHIN STATE        
!** TRANSVERSE MERCATOR PROJECTION ZONES, SEVERAL MINOR TERMS OF        
!** THE EQUATIONS MAY BE OMITTED (SEE A SEPARATE NGS PUBLICATION).      
!** IF PROGRAMMED IN FULL, THE SUBROUTINE CAN BE USED FOR               
!** COMPUTATIONS IN SURVEYS EXTENDING OVER TWO ZONES.                   
!********************************************************************** 
      OM = (N - FN + SO) / (R * SF) 
      COSOM = DCOS (OM) 
      FOOT = OM + DSIN (OM) * COSOM * (V0 + V2 * COSOM * COSOM + V4 *   &
      COSOM**4 + V6 * COSOM**6)                                         
      SINF = DSIN (FOOT) 
      COSF = DCOS (FOOT) 
      TN = SINF / COSF 
      TS = TN * TN 
      ETS = EPS * COSF * COSF 
      RN = ER * SF / DSQRT (1.D0 - ESQ * SINF * SINF) 
      Q = (E-FE) / RN 
      QS = Q * Q 
      B2 = - TN * (1.D0 + ETS) / 2.D0 
      B4 = - (5.D0 + 3.D0 * TS + ETS * (1.D0 - 9.D0 * TS) - 4.D0 * ETS *&
      ETS) / 12.D0                                                      
      B6 = (61.D0 + 45.D0 * TS * (2.D0 + TS) + ETS * (46.D0 - 252.D0 *  &
      TS - 60.D0 * TS * TS) ) / 360.D0                                  
      B1 = 1.D0 
      B3 = - (1.D0 + TS + TS + ETS) / 6.D0 
      B5 = (5.D0 + TS * (28.D0 + 24.D0 * TS) + ETS * (6.D0 + 8.D0 * TS) &
      ) / 120.D0                                                        
      B7 = - (61.D0 + 662.D0 * TS + 1320.D0 * TS * TS + 720.D0 * TS**3) &
      / 5040.D0                                                         
      LAT = FOOT + B2 * QS * (1.D0 + QS * (B4 + B6 * QS) ) 
      L = B1 * Q * (1.D0 + QS * (B3 + QS * (B5 + B7 * QS) ) ) 
      LON = - L / COSF + CM 
                                                                        
!*********************************************************************  
!     COMPUTE CONVERENCE AND SCALE FACTOR                               
      FI = LAT 
      LAM = LON 
      SINFI = SIN (FI) 
      COSFI = COS (FI) 
      L1 = (LAM - CM) * COSFI 
      LS = L1 * L1 
                                                                        
      TN = SINFI / COSFI 
      TS = TN * TN 
!                                                                       
!*** CONVERGENCE                                                        
      C1 = - TN 
      C3 = (1. + 3. * ETS + 2. * ETS**2) / 3. 
      C5 = (2. - TS) / 15. 
      CONV = C1 * L1 * (1. + LS * (C3 + C5 * LS) ) 
!                                                                       
!*** POINT SCALE FACTOR                                                 
      F2 = (1. + ETS) / 2. 
      F4 = (5. - 4. * TS + ETS * (9. - 24. * TS) ) / 12. 
      KP = SF * (1. + F2 * LS * (1. + F4 * LS) ) 
                                                                        
      RETURN 
      END SUBROUTINE TMGEOD                         

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

!*************************************************************************
!                           SUBROTINA INTERF                             *
!*************************************************************************

    Subroutine Interf2D (FatherGrid, X, Y, PP, ILB, JLB, IUB, JUB, Abs1, Ord1, F)

        !Arguments-------------------------------------------------------------

        type (T_FatherGrid),     pointer  :: FatherGrid
        real, dimension (:,:),   pointer  :: F
        real, dimension (:  ),   pointer  :: Abs1, Ord1
        real                              :: X, Y, PP
        integer                           :: IUB, JUB, ILB, JLB

        !Local------------------------------------------------------------
        integer                           :: JX, IY, J1, J2, I1, I2
        real                              :: Fhc, RJ, PASX, RI, PASY, DF, FJ2, FJ1

        !Begin------------------------------------------------------------

        JX  = FatherGrid%JX
        IY  = FatherGrid%IY
        Fhc = FatherGrid%Fhc


      100 Continue
          IF (X.EQ.ABS1(JUB)) Then
             JX = JUB
             J1 = JUB-1
             J2 = JUB
             Goto 200
          Endif
          IF (X.LT.ABS1(JX)) Goto 110
          JX = JX+1
          IF (JX.GE.JUB) Goto 120
          Goto 100
      110 Continue
          JX = JX-1
          IF (JX.LT.JLB) Goto 330
          IF (X.LT.ABS1(JX)) Goto 110
          J1 = JX
          J2 = JX+1
          Goto 200
      120 Continue
          IF (JX.NE.JUB) Goto 330
          J1 = JX-1
          J2 = JX
      200 Continue
          IF (Y.EQ.ORD1(IUB)) Then
             IY = IUB
             I1 = IUB-1
             I2 = IUB
             Goto 300
          Endif
          IF (Y.LT.ORD1(IY)) Goto 210
          IY = IY+1
          IF (IY.GE.IUB) Goto 220
          Goto 200
      210 Continue
          IY = IY-1
          IF (IY.LT.ILB) Goto 330
          IF (Y.LT.ORD1(IY)) Goto 210
          I1 = IY
          I2 = IY+1
          Goto 300
      220 Continue
          IF (IY.NE.IUB) Goto 330
          I1 = IY-1
          I2 = IY
      300 Continue

          RJ   = X-ABS1(J1)
          PASX = ABS1(J2)-ABS1(J1)
          RI   = Y-ORD1(I1)
          PASY = ORD1(I2)-ORD1(I1)
          DF   = (F(I2,J2)-F(I1,J2))/PASY

          FJ2  = F(I1,J2)+DF*RI
          DF   = (F(I2,J1)-F(I1,J1))/PASY
          FJ1  = F(I1,J1)+DF*RI
          DF   = (FJ2-FJ1)/PASX
          PP   = FJ1+RJ*DF
          IF (JX.GT.JUB) JX=JUB
          IF (IY.GT.IUB) IY=IUB

          RETURN

      330 Continue
          PP = FHC
          IF (JX.LT.JLB) JX=JLB
          IF (JX.GT.JUB) JX=JUB
          IF (IY.LT.ILB) IY=ILB
          IF (IY.GT.IUB) IY=IUB

    end Subroutine Interf2D

!------------------------------------------------------------------------------
    
    Subroutine Interf3D (FatherGrid, X, Y, PP, ILB, JLB, IUB, JUB, K, Abs1, Ord1, F)

        !Arguments-------------------------------------------------------------

        type (T_FatherGrid),     pointer  :: FatherGrid
        real, dimension (:,:,:), pointer  :: F
        real, dimension (:  ),   pointer  :: Abs1, Ord1
        real                              :: X, Y, PP
        integer                           :: IUB, JUB, ILB, JLB, K

        !Local------------------------------------------------------------
        integer                           :: JX, IY, J1, J2, I1, I2
        real                              :: Fhc, RJ, PASX, RI, PASY, DF, FJ2, FJ1

        !Begin------------------------------------------------------------

        JX  = FatherGrid%JX
        IY  = FatherGrid%IY
        Fhc = FatherGrid%Fhc


      100 Continue
          IF (X.EQ.ABS1(JUB)) Then
             JX = JUB
             J1 = JUB-1
             J2 = JUB
             Goto 200
          Endif
          IF (X.LT.ABS1(JX)) Goto 110
          JX = JX+1
          IF (JX.GE.JUB) Goto 120
          Goto 100
      110 Continue
          JX = JX-1
          IF (JX.LT.JLB) Goto 330
          IF (X.LT.ABS1(JX)) Goto 110
          J1 = JX
          J2 = JX+1
          Goto 200
      120 Continue
          IF (JX.NE.JUB) Goto 330
          J1 = JX-1
          J2 = JX
      200 Continue
          IF (Y.EQ.ORD1(IUB)) Then
             IY = IUB
             I1 = IUB-1
             I2 = IUB
             Goto 300
          Endif
          IF (Y.LT.ORD1(IY)) Goto 210
          IY = IY+1
          IF (IY.GE.IUB) Goto 220
          Goto 200
      210 Continue
          IY = IY-1
          IF (IY.LT.ILB) Goto 330
          IF (Y.LT.ORD1(IY)) Goto 210
          I1 = IY
          I2 = IY+1
          Goto 300
      220 Continue
          IF (IY.NE.IUB) Goto 330
          I1 = IY-1
          I2 = IY
      300 Continue

          RJ   = X-ABS1(J1)
          PASX = ABS1(J2)-ABS1(J1)
          RI   = Y-ORD1(I1)
          PASY = ORD1(I2)-ORD1(I1)
          DF   = (F(I2,J2, K)-F(I1,J2, K))/PASY

          FJ2  = F(I1,J2,K)+DF*RI
          DF   = (F(I2,J1,K)-F(I1,J1,K))/PASY
          FJ1  = F(I1,J1,K)+DF*RI
          DF   = (FJ2-FJ1)/PASX
          PP   = FJ1+RJ*DF
          IF (JX.GT.JUB) JX=JUB
          IF (IY.GT.IUB) IY=IUB

          return

      330 Continue
          PP = FHC
          IF (JX.LT.JLB) JX=JLB
          IF (JX.GT.JUB) JX=JUB
          IF (IY.LT.ILB) IY=ILB
          IF (IY.GT.IUB) IY=IUB
    

        
    end Subroutine Interf3D
                          
!------------------------------------------------------------------------------

    Subroutine Interf3D8 (FatherGrid, X, Y, PP, ILB, JLB, IUB, JUB, K, Abs1, Ord1, F)

        !Arguments-------------------------------------------------------------

        type (T_FatherGrid),        pointer  :: FatherGrid
        real(8), dimension (:,:,:), pointer  :: F
        real, dimension (:  ),      pointer  :: Abs1, Ord1
        real                                 :: X, Y
        real(8)                              :: PP
        integer                              :: IUB, JUB, ILB, JLB, K

        !Local------------------------------------------------------------
        integer                              :: JX, IY, J1, J2, I1, I2
        real                                 :: Fhc, RJ, PASX, RI, PASY, DF, FJ2, FJ1

        !Begin------------------------------------------------------------

        JX  = FatherGrid%JX
        IY  = FatherGrid%IY
        Fhc = FatherGrid%Fhc


      100 Continue
          IF (X.EQ.ABS1(JUB)) Then
             JX = JUB
             J1 = JUB-1
             J2 = JUB
             Goto 200
          Endif
          IF (X.LT.ABS1(JX)) Goto 110
          JX = JX+1
          IF (JX.GE.JUB) Goto 120
          Goto 100
      110 Continue
          JX = JX-1
          IF (JX.LT.JLB) Goto 330
          IF (X.LT.ABS1(JX)) Goto 110
          J1 = JX
          J2 = JX+1
          Goto 200
      120 Continue
          IF (JX.NE.JUB) Goto 330
          J1 = JX-1
          J2 = JX
      200 Continue
          IF (Y.EQ.ORD1(IUB)) Then
             IY = IUB
             I1 = IUB-1
             I2 = IUB
             Goto 300
          Endif
          IF (Y.LT.ORD1(IY)) Goto 210
          IY = IY+1
          IF (IY.GE.IUB) Goto 220
          Goto 200
      210 Continue
          IY = IY-1
          IF (IY.LT.ILB) Goto 330
          IF (Y.LT.ORD1(IY)) Goto 210
          I1 = IY
          I2 = IY+1
          Goto 300
      220 Continue
          IF (IY.NE.IUB) Goto 330
          I1 = IY-1
          I2 = IY
      300 Continue

          RJ   = X-ABS1(J1)
          PASX = ABS1(J2)-ABS1(J1)
          RI   = Y-ORD1(I1)
          PASY = ORD1(I2)-ORD1(I1)
          DF   = (F(I2,J2, K)-F(I1,J2, K))/PASY

          FJ2  = F(I1,J2,K)+DF*RI
          DF   = (F(I2,J1,K)-F(I1,J1,K))/PASY
          FJ1  = F(I1,J1,K)+DF*RI
          DF   = (FJ2-FJ1)/PASX
          PP   = FJ1+RJ*DF
          IF (JX.GT.JUB) JX=JUB
          IF (IY.GT.IUB) IY=IUB

          return

      330 Continue
          PP = FHC
          IF (JX.LT.JLB) JX=JLB
          IF (JX.GT.JUB) JX=JUB
          IF (IY.LT.ILB) IY=ILB
          IF (IY.GT.IUB) IY=IUB
        
        

    end Subroutine Interf3D8
        

!------------------------------------------------------------------------------
        
    Subroutine LocateCell  (XX, YY, XPos, YPos,                                         &
                            ILB, IUB, JLB, JUB, ILower, Jlower)


        !Arguments---------------------------------------------------------
        real,   dimension(:)  , pointer      :: XX, YY
        real                  , intent (IN ) :: XPos, YPos
        integer               , intent (IN ) :: ILB, IUB, JLB, JUB
        integer               , intent (OUT) :: ILower, Jlower

        
        !Local-------------------------------------------------------------
        integer                              :: IMiddle, JMiddle, Dcd1, Dcd2, IUpper, JUpper 

        !Begin------------------------------------------------------------- 


        ILower = ILB
        IUpper = IUB

        Dcd1 = 2
        Dcd2 = 2

        if(Ypos < YY(ILB) .or. Ypos > YY(IUB))then

            ILower = null_int

        else

            do while (Dcd1 > 1)
                IMiddle = int((ILower + IUpper)/2)
                if (Ypos > YY(IMiddle)) then
                    ILower = IMiddle
                else 
                    IUpper = IMiddle
                endif
                Dcd1 = IUpper - ILower
            enddo

        end if


        Jlower = JLB
        JUpper = JUB

        if(Xpos < XX(JLB) .or. Xpos > XX(JUB))then

            Jlower = null_int

        else

            do while (Dcd2 > 1)
                JMiddle = int((Jlower + JUpper)/2)
                if (Xpos > XX(JMiddle)) then
                    Jlower = JMiddle
                else 
                    JUpper = JMiddle
                endif
                Dcd2 = JUpper - Jlower
            enddo

        end if

        if (.not. Dcd1>0 .and. Dcd2>0) then

            call SetError(FATAL_, INTERNAL_, 'LocateCell - HorizontalGrid - ERR01')

           
        endif

    end subroutine LocateCell


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

    real  function Bilinear (YYUpper, YYLower, YPos,                                     &
                            XXUpper, XXLower, XPos,                                      &  
                            PropLowLeft,  PropUpLeft, PropLowRight, PropUpRight,         &
                            ONLowLeft,    ONUpLeft,   ONLowRight,   ONUpRight)

        !Arguments---------------------------------------------------------
        real    :: YYUpper, YYLower, YPos
        real    :: XXUpper, XXLower, XPos
        real    :: PropLowLeft, PropUpLeft, PropLowRight, PropUpRight
        integer :: ONLowLeft,   ONUpLeft,   ONLowRight,   ONUpRight

        !Local-------------------------------------------------------------
        real    :: Prop1, Prop2, Dx3, Dx2, Dy2, Dy3, D1, D2, D3
        integer :: ON1, ON2

        !Begin------------------------------------------------------------- 

        Dy2   = YYUpper - YPos
        Dy3   = YPos    - YYLower

        Dx2   = XXUpper - XPos
        Dx3   = XPos    - XXLower

        D1    =  Dy2 * ONLowLeft  + Dy3 * ONUpLeft
        D2    =  Dy2 * ONLowRight + Dy3 * ONUpRight 

cd1:    if (D1 /= 0) then

            Prop1 = (Dy2 * ONLowLeft * PropLowLeft  + Dy3 * ONUpLeft * PropUpLeft) /  D1

        else if (ONLowLeft == 1) then cd1

            Prop1 = PropLowLeft

        else if (ONUpLeft  == 1) then cd1

            Prop1 = PropUpLeft

        else  cd1

            Prop1 = 0.

        endif cd1

cd2:    if (D2 /= 0) then

            Prop2 = (Dy2 * ONLowRight * PropLowRight  + Dy3 * ONUpRight * PropUpRight) /  D2

        else if (ONLowRight == 1) then cd2

            Prop2 = PropLowRight

        else if (ONUpRight  == 1) then cd2

            Prop2 = PropUpRight

        else  cd2

            Prop2 = 0.

        endif cd2


        ON1   = ONUpLeft   + ONLowLeft

        if (ON1 > 1) ON1 = 1

        ON2   = ONLowRight + ONUpRight

        if (ON2 > 1) ON2 = 1

        D3 = Dx2 * ON1 + Dx3 * ON2 


cd3:    if (D3 /= 0) then

            Bilinear  = (Dx2 * ON1 * Prop1 + Dx3 * ON2 * Prop2) / D3

        else if (ON1 == 1 ) then cd3

            Bilinear = Prop1

        else if (ON2  == 1) then cd3

            Bilinear = Prop2

        else  cd3

            Bilinear = 0.

        endif cd3

            

    end function Bilinear

!------------------------------------------------------------------------------
    
    subroutine InterpolPoint(XXSon, YYSon,                                              &
                             YYUpper, YYLower, XXUpper, XXLower,                        &
                             PropLowLeft,  PropUpLeft, PropLowRight, PropUpRight,       &
                             ONLowLeft,    ONUpLeft,   ONLowRight,   ONUpRight,         &
                             PropSon, OK)

        !Arguments ---------------------------------------------------------
        real   ,                 intent(IN)  :: XXSon, YYSon,                           &
                                                YYUpper, YYLower, XXUpper, XXLower,     &
                                                PropLowLeft,  PropUpLeft, PropLowRight, &
                                                PropUpRight
        integer,                 intent(IN ) :: ONLowLeft, ONUpLeft, ONLowRight,        &
                                                ONUpRight
        real   ,                 intent(OUT) :: PropSon
        logical                , intent(OUT) :: OK

        !Local -------------------------------------------------------------
        integer                            :: SumON

        !Begin -------------------------------------------------------------
        
        SumON = ONLowLeft + ONUpLeft + ONLowRight + ONUpRight

cd1:    if      (SumON > 0) then

            PropSon = Bilinear(YYUpper, YYLower, YYSon,                                 &
                               XXUpper, XXLower, XXSon,                                 &  
                               PropLowLeft,  PropUpLeft, PropLowRight, PropUpRight,     &
                               ONLowLeft,    ONUpLeft,   ONLowRight,   ONUpRight)

            OK = .True.

        else cd1
           
            OK = .False.            

        endif cd1


    end subroutine InterpolPoint



   !------------------------------------------------------------------------    
    
    subroutine LocateCellPolygons(XX2D_Z, YY2D_Z, XPoint, YPoint,                       &
                                  DefinedPoint, ILB, IUB, JLB, JUB, IZ, JZ)

        !Arguments ---------------------------------------------------------
        real,    dimension(:,:)  , pointer   :: XX2D_Z, YY2D_Z
        integer, dimension(:,:)  , pointer   :: DefinedPoint
        real   ,                 intent(IN ) :: XPoint, YPoint
        integer,                 intent(IN ) :: ILB, IUB, JLB, JUB
        integer,                 intent(OUT) :: IZ, JZ

        !Local -------------------------------------------------------------
        integer                              :: ICenter, JCenter, Iupper, ILower,       &
                                                Jupper, JLower
        logical                              :: CellFound
        type(T_PointF ),          pointer    :: Point
        type(T_Polygon),          pointer    :: Polygon
        integer                              :: i, j, pi
        integer                              :: I1, I2, I3, I4
        integer                              :: J1, J2, J3, J4
        !Begin -------------------------------------------------------------

        allocate(Polygon)
        allocate(Polygon%VerticesF(1: 2*(IUB-ILB+1+JUB-JLB+1)))

        allocate(Point)

        CellFound = .false.

        Iupper = IUB
        ILower = ILB

        Jupper = JUB
        JLower = JLB

        do while (.not. CellFound)

            ICenter = int(real((Iupper - ILower)/ 2.)) + ILower
            JCenter = int(real((Jupper - JLower)/ 2.)) + JLower

            !Construct 4 polygons NW, NE, SW, SE
            do pi = 1, 4
                
                !All polygons are defined anti-Clockwise 
                if      (pi==1) then

                    if (ILower == ICenter) cycle
                    if (JLower == JCenter) cycle
                    
                    I1 = ILower;  J1 = JLower
                    I2 = ILower;  J2 = JCenter
                    I3 = ICenter; J3 = JCenter
                    I4 = ICenter; J4 = JLower

                else if (pi==2) then

                    if (ILower == ICenter) cycle

                    I1 = ILower;  J1 = JCenter;
                    I2 = ILower;  J2 = Jupper;
                    I3 = ICenter; J3 = Jupper;
                    I4 = ICenter; J4 = JCenter;

                else if (pi==3) then

                    if (JLower == JCenter) cycle

                    I1 = ICenter; J1 = JLower;
                    I2 = ICenter; J2 = JCenter;
                    I3 = IUpper ; J3 = JCenter;
                    I4 = IUpper;  J4 = JLower;

                else if (pi==4) then


                    I1 = ICenter; J1 = JCenter;
                    I2 = ICenter; J2 = JUpper;
                    I3 = IUpper ; J3 = JUpper;
                    I4 = IUpper;  J4 = JCenter;

                endif
                
                Polygon%Count = 0.

                do j = J1, J2

                    if (DefinedPoint(I1,j) == 1) then

                        Polygon%Count = Polygon%Count + 1                   

                        Polygon%VerticesF(Polygon%Count)%X  = XX2D_Z(I1,j)
                        Polygon%VerticesF(Polygon%Count)%Y  = YY2D_Z(I1,j)

                    endif
                
                end do

                do i = I2, I3

                    if (DefinedPoint(i,J2) == 1) then

                        Polygon%Count = Polygon%Count + 1                   

                        Polygon%VerticesF(Polygon%Count)%X  = XX2D_Z(i,J2)
                        Polygon%VerticesF(Polygon%Count)%Y  = YY2D_Z(i,J2)

                    endif
                
                end do

                do j = J3, J4, -1

                    if (DefinedPoint(I3,j) == 1) then

                        Polygon%Count = Polygon%Count + 1                   

                        Polygon%VerticesF(Polygon%Count)%X  = XX2D_Z(I3,j)
                        Polygon%VerticesF(Polygon%Count)%Y  = YY2D_Z(I3,j)

                    endif
                
                end do

                do i = I4, I1, -1

                    if (DefinedPoint(i,J4) == 1) then

                        Polygon%Count = Polygon%Count + 1                   

                        Polygon%VerticesF(Polygon%Count)%X  = XX2D_Z(i,J4)
                        Polygon%VerticesF(Polygon%Count)%Y  = YY2D_Z(i,J4)

                    endif
                
                end do

                Polygon%Count = Polygon%Count + 1                   

                !Close polygon
                Polygon%VerticesF(Polygon%Count)%X  = Polygon%VerticesF(1)%X
                Polygon%VerticesF(Polygon%Count)%Y  = Polygon%VerticesF(1)%Y

                call SetLimits(Polygon)

                Point%X = XPoint
                Point%Y = YPoint

                if (IsPointInsidePolygon(Point, Polygon)) exit

            enddo

            Iupper  = I3
            ILower  = I1

            Jupper  = J2
            JLower  = J1

            !Test if the cell was founded
            if ((Iupper - ILower) == 1 .and. (Jupper - JLower) == 1) CellFound = .true.

        end do 

        IZ = ILower
        JZ = JLower

            
        deallocate(Polygon%VerticesF)        
        nullify   (Polygon%VerticesF)        

        deallocate(Polygon)
        nullify   (Polygon)

        deallocate(Point)
        nullify   (Point)


    end subroutine LocateCellPolygons

end module ModuleHorizontalGrid





