!------------------------------------------------------------------------------
!        IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE         : Mohid Model
! PROJECT       : Mohid Base 2
! MODULE        : Geometry
! URL           : http://www.mohid.com
! AFFILIATION   : IST/MARETEC, Marine Modelling Group
! DATE          : May 2003
! REVISION      : Frank Braunschweig - v4.0
! DESCRIPTION   : Module to calculate the vertical Geometry 
!
!------------------------------------------------------------------------------
! IMPERMEABILITY                : 0/1               -           !Consider impermeable cell faces
! IMPER_COEF_U                  : real             [1]          !
! IMPER_COEFX_U                 : real             [0]          !
! IMPER_COEF_V                  : real             [1]          !
! IMPER_COEFX_V                 : real             [0]          !

!<begindomain>
!   ID                          : int               -           !Domain ID
!   TYPE                        : char              -           !Type of vertical coordinate of the domain
!                                                               !Multiple options: FIXSPACING, SIGMA,
!                                                               !LAGRANGIAN, CARTESIAN, HARMONIC, FIXSEDIMENT
!   LAYERS                      : int               -           !Number of layers
!   EQUIDISTANT                 : real             [0]          !Equidistant layers spacing in meters
!   LAYERTHICKNESS              : real vector       -           !If not equidistant specifies layers thickness
!                                                               !starting from bottom layer (e.g. 50. 20. 10. 5.)
!   TOLERANCEDEPTH              : real            [0.05]        !Just for SIGMA,ISOPYCNIC, LAGRANGIAN coordinates
!   TOTALTHICKNESS              : real              -           !Total domain thickness 
!                                                               !(Just for FIXSPACING, FIXSEDIMENT, SOIL_TOPLAYER)
!   EMPTY_TOP_LAYERS            : int              [0]          !Number of empty layers counting from top
!   DOMAINDEPTH                 : real

Module ModuleGeometry

    use ModuleGlobalData
    use ModuleTime                 
    use ModuleEnterData          
    use ModuleGridData,         only: GetGridData, UngetGridData, GetMaximumValue,      &
                                      GetGridDataFileName, WriteGridData,               &
                                      GetGridDataEvolution     
    use ModuleHorizontalMap,    only: GetWaterPoints2D, GetExteriorBoundaryFaces,       &
                                      UnGetHorizontalMap
    use ModuleHorizontalGrid,   only: GetHorizontalGridSize, GetHorizontalGrid,         &
                                      GetGridCoordType, GetGridAngle, GetGridZone,      &
                                      GetLatitudeLongitude, GetGridOrigin,              &
                                      GetCheckDistortion, UnGetHorizontalGrid
    use ModuleHDF5

    implicit none
    private

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


    !Constructor
    public  :: ConstructGeometry
    private ::      AllocateInstance
    private ::      ConstructGlobalVariables
    private ::          GetDomainsFromFile
    private ::              Add_Domain
    private ::          ComputeLayers
    private ::          AllocateVariables
    private ::          VerifyBathymetry
    private ::      ConstructKFloor


    !Modifier
    public  :: ComputeInitialGeometry       !To call when the initial SurfaceElevation is known
    public  :: ComputeVerticalGeometry      !To call in the "working cycle"
    private ::      ComputeSZZ
    private ::          ComputeFixSpacing
    private ::          ComputeFixSediment
    private ::          ComputeCartesian
    private ::          ComputeHarmonic
    private ::          ComputeSigma
    private ::          ComputeLagrangian
    private ::      ComputeZCellCenter
    private ::      ComputeDistances
    private ::      ComputeAreas
    private ::      ComputeVolumes
    private ::      StoreVolumeZOld


    !Input / Output
    public  :: ReadGeometry
    public  :: WriteGeometry

    !Selector
    public  :: GetGeometryDistances                 !SZZ, DZZ, DWZ, ZCellCenter, DUZ, DVZ
    public  :: GetGeometryAreas
    public  :: GetGeometryVolumes
    public  :: GetGeometryKFloor
    public  :: GetGeometryKTop
    public  :: GetGeometryWaterColumn               !HT
    public  :: GetGeometryInputFile
    public  :: GetGeometrySize
    public  :: GetGeometryMinWaterColumn
    public  :: UnGetGeometry


    !Destructor
    public  :: KillGeometry
    private ::      DeallocateInstance
    private ::      DeallocateVariables


    !Management
    private ::      Ready
    private ::          LocateObjGeometry

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

    private :: UnGetGeometry2Dreal
    private :: UnGetGeometry3Dreal4
    private :: UnGetGeometry3Dreal8
    private :: UnGetGeometry2Dinteger
    interface  UnGetGeometry
        module procedure UnGetGeometry2Dreal
        module procedure UnGetGeometry3Dreal4
        module procedure UnGetGeometry3Dreal8
        module procedure UnGetGeometry2Dinteger
    end interface  UnGetGeometry

    private :: ReadGeometryHDF
    private :: ReadGeometryBin
    interface  ReadGeometry
        module procedure ReadGeometryBin
        module procedure ReadGeometryHDF
    end interface

    private :: WriteGeometryHDF
    private :: WriteGeometryBin
    interface  WriteGeometry
        module procedure WriteGeometryBin
        module procedure WriteGeometryHDF
    end interface 

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

    !STAT
    real   , parameter :: FillValueDouble = -9.9e15

    !Domain types
    integer, parameter :: FixSpacing            = 1
    integer, parameter :: Sigma                 = 2
    integer, parameter :: Isopycnic             = 3
    integer, parameter :: Lagrangian            = 4
    integer, parameter :: Cartesian             = 5
    integer, parameter :: Harmonic              = 6
    integer, parameter :: FixSediment           = 7
    integer, parameter :: SigmaTop              = 8
    integer, parameter :: CartesianTop          = 9


    !Other
    integer, parameter :: INITIALGEOMETRY       = 1
    integer, parameter :: TRANSIENTGEOMETRY     = 2

    !Faces Option
    integer, parameter :: MinTickness           = 3
    integer, parameter :: AverageTickness       = 2


    !Type----------------------------------------------------------------------

    type T_Domain
        integer                                 :: ID                        = FillValueInt
        integer                                 :: DomainType                = FillValueInt
        integer                                 :: NumberOfLayers            = FillValueInt
        integer                                 :: EmptyTopLayers            = FillValueInt
        integer                                 :: UpperLayer, LowerLayer    = FillValueInt
        real                                    :: DomainDepth               = FillValueReal
        real                                    :: TotalThickness            = FillValueReal 
        real, dimension(:), pointer             :: LayerThickness
        real                                    :: ToleranceDepth            = FillValueReal 
        real                                    :: MinInitialLayerThickness  = FillValueReal
        real                                    :: MinEvolveLayerThickness   = FillValueReal
        real                                    :: MinEsp                    = FillValueReal
        real                                    :: GridMovementDump          = FillValueReal
        real                                    :: DisplacementLimit         = FillValueReal 
        integer                                 :: InitializationMethod      = FillValueInt
        real                                    :: Equidistant               = FillValueReal        
        type (T_Domain), pointer                :: Next, Prev
    end type T_Domain

    type T_Distances
        real, dimension(:, :, :), pointer       :: SZZ 
        real, dimension(:, :, :), pointer       :: DZZ 
        real, dimension(:, :, :), pointer       :: DWZ, DUZ, DVZ, DZI, DZE 
        real, dimension(:, :, :), pointer       :: InitialSZZ 
        real, dimension(:, :, :), pointer       :: ZCellCenter  !Distance from the refernce level, 
                                                                ! center of cells, positive upwards
    end type T_Distances

    type T_Areas
        real, dimension(:, :, :), pointer       :: AreaU, AreaV
        logical                                 :: Impermeability = .false.
        real, dimension(      :), pointer       :: Coef_U, CoefX_U, Coef_V, CoefX_V  
    end type T_Areas

    type T_Volumes
        real(8), dimension(:, :, :), pointer    :: VolumeZ, VolumeU, VolumeV, VolumeW, VolumeZOld
        logical                                 :: FirstVolW = .true.
    end type T_Volumes

    type T_KFloor
        integer, dimension(:, :), pointer       :: Z, U, V, Domain
    end type T_KFloor

    type T_KTop
        integer, dimension(:, :), pointer       :: Z
    end type T_KTop


    type T_WaterColumn
        real, dimension(:, :), pointer          :: U, V
        real, dimension(:, :), pointer          :: Z                    !Former HT
        real                                    :: Zmin = FillValueReal !Former Hmin
    end type T_WaterColumn

    type T_External
        logical                                 :: ContinuesCompute = .false.
        logical                                 :: NonHydrostatic   = .false.
        real                                    :: BathymTopoFactor = 1.0
        real, dimension(:,:,:), pointer         :: DecayTime
    end type T_External

    type T_Geometry
        integer                                 :: InstanceID
        integer                                 :: FacesOption
        type (T_External)                       :: ExternalVar
        type (T_Distances)                      :: Distances
        type (T_Areas)                          :: Areas
        type (T_Volumes)                        :: Volumes
        type (T_WaterColumn)                    :: WaterColumn      
        type (T_Domain), pointer                :: FirstDomain
        type (T_Domain), pointer                :: LastDomain
        type (T_KFloor)                         :: KFloor
        type (T_KTop)                           :: KTop

        type (T_Time)                           :: ActualTime  
        type (T_Size3D)                         :: Size
        type (T_Size3D)                         :: WorkSize
       
        character(len=Pathlength)               :: InputFile

        !Instance of other modules
        integer                                 :: ObjTopography        = 0
        integer                                 :: ObjHorizontalGrid    = 0
        integer                                 :: ObjHorizontalMap     = 0

        type (T_Geometry), pointer              :: Next

    end type T_Geometry

    !Global Module Variables
    type (T_Geometry), pointer                  :: FirstGeometry
    type (T_Geometry), pointer                  :: Me
   

    contains

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

    !CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR 

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

    subroutine ConstructGeometry(GeometryID, GridDataID, HorizontalGridID,              &
                                 HorizontalMapID, ActualTime,                           &
                                 NewDomain, SurfaceElevation, BathymTopoFactor, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        integer                                     :: GridDataID
        integer                                     :: HorizontalGridID
        integer                                     :: HorizontalMapID
        type (T_Time),                 optional     :: ActualTime
        character (len=*), intent(IN), optional     :: NewDomain
        real, dimension(:,:), pointer, optional     :: SurfaceElevation
        real,    intent(in),           optional     :: BathymTopoFactor
        integer,  intent(OUT),         optional     :: STAT
        
        !Local-----------------------------------------------------------------
        integer                                     :: STAT_, ready_

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

        STAT_ = UNKNOWN_

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

        call Ready(GeometryID, ready_)    

        if (ready_ .EQ. OFF_ERR_) then

            !Allocates Instance
            call AllocateInstance

            !Associates External Instances
            Me%ObjTopography     = AssociateInstance (mGRIDDATA_,       GridDataID      )
            Me%ObjHorizontalMap  = AssociateInstance (mHORIZONTALMAP_,  HorizontalMapID )
            Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)

            !Actualize the time
            if (present(ActualTime)) then
                Me%ActualTime = ActualTime
            endif

            !Factor For MOHID Water / MOHID Land
            if (present(BathymTopoFactor)) then
                Me%ExternalVar%BathymTopoFactor = BathymTopoFactor
            else
                Me%ExternalVar%BathymTopoFactor = 1.0
            endif

            !Construct the variable common to all module
            if (present(NewDomain)) then
                if (present(SurfaceElevation)) then
                    call ConstructGlobalVariables(NewDomain = NewDomain, SurfaceElevation = SurfaceElevation)
                else
                    call ConstructGlobalVariables(NewDomain = NewDomain)
                endif
            else
                if (present(SurfaceElevation)) then
                    call ConstructGlobalVariables (SurfaceElevation = SurfaceElevation)
                else
                    call ConstructGlobalVariables
                endif
            endif

            !Constructs the Matrixes which contains the Ks - KFloor, etc...
            if (present(SurfaceElevation)) then
                call ConstructKFloor (SurfaceElevation)
            else
                call ConstructKFloor
            endif    

            !Returns ID
            GeometryID    = Me%InstanceID

            STAT_ = SUCCESS_

        else 
            
            stop 'Geometry - ConstructGeometry - ERR99' 

        end if 


        if (present(STAT)) STAT = STAT_


    end subroutine ConstructGeometry
 
    !--------------------------------------------------------------------------

    subroutine AllocateInstance

        !Arguments-------------------------------------------------------------
    
        !Local-----------------------------------------------------------------
        type (T_Geometry), pointer                  :: NewObjGeometry
        type (T_Geometry), pointer                  :: PreviousObjGeometry


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

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

        Me%InstanceID = RegisterNewInstance (mGEOMETRY_)

    end subroutine AllocateInstance

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

    subroutine ConstructGlobalVariables(NewDomain, SurfaceElevation)

        !Arguments-------------------------------------------------------------
        character (len=*), intent(IN), optional                 :: NewDomain
        real, dimension(:,:), pointer, optional                 :: SurfaceElevation
        
        !External--------------------------------------------------------------
        integer                                                 :: STAT_CALL
        character(len = StringLength)                           :: Message
        type (T_Size2D)                                         :: Size2D, WorkSize2D

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

        !Nullify ObjGeometryVariables
        nullify (Me%FirstDomain)

        !Nullify T_Volumes
        nullify (Me%Volumes%VolumeZ)
        nullify (Me%Volumes%VolumeU)
        nullify (Me%Volumes%VolumeV)
        nullify (Me%Volumes%VolumeW)
        nullify (Me%Volumes%VolumeZOld)

        !Nullify T_Areas
        nullify (Me%Areas%AreaU)
        nullify (Me%Areas%AreaV)

        !Nullify T_Distances
        nullify (Me%Distances%SZZ        )
        nullify (Me%Distances%DZZ        )
        nullify (Me%Distances%DWZ        )
        nullify (Me%Distances%DUZ        )
        nullify (Me%Distances%DZI        )
        nullify (Me%Distances%DZE        )
        nullify (Me%Distances%DVZ        )
        nullify (Me%Distances%InitialSZZ )
        nullify (Me%Distances%ZCellCenter)

        !Nullify T_KFloor
        nullify (Me%KFloor%Z)
        nullify (Me%KFloor%U)
        nullify (Me%KFloor%V)
        nullify (Me%KFloor%Domain)
       
        !Gets horizontal size from the Bathymetry
        call GetHorizontalGridSize(Me%ObjHorizontalGrid,                                 &
                                   Size     = Size2D,                                    &
                                   WorkSize = WorkSize2D,                                &
                                   STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - Geometry - ERR02'

        Me%Size%ILB = Size2D%ILB
        Me%Size%IUB = Size2D%IUB
        Me%Size%JLB = Size2D%JLB
        Me%Size%JUB = Size2D%JUB

        Me%WorkSize%ILB = WorkSize2D%ILB
        Me%WorkSize%IUB = WorkSize2D%IUB
        Me%WorkSize%JLB = WorkSize2D%JLB
        Me%WorkSize%JUB = WorkSize2D%JUB

        if (present(NewDomain)) then
            Me%InputFile = NewDomain
        else
            !Reads the file name of the domain properties
            Message   ='File of the domain properties.'
            call ReadFileName('DOMAIN', Me%InputFile, Message = Message, STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - Geometry - ERR03'
        endif 

        !Reads the Domain Properties
        call GetDomainsFromFile

        !Calculates / Verifies Layers / Computes the bounds in K
        call ComputeLayers

        call ConstructImpermeability

        !Checks if Bathymetry is consistent with the tolerance depth - Fromer REBAIXA
        !and if changes bathymetry if cartasian domain type exists
        if (present(SurfaceElevation)) then
            call VerifyBathymetry (SurfaceElevation)
        else
            call VerifyBathymetry
        endif

        !Allocates variables
        call AllocateVariables

    end subroutine ConstructGlobalVariables

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

    subroutine ConstructImpermeability

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

        !Local-----------------------------------------------------------------
        integer                                     :: ObjEnterData = 0   
        integer                                     :: STATUS
        integer                                     :: iflag
        integer                                     :: WorkKLB, WorkKUB

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

        call ConstructEnterData(ObjEnterData, Me%InputFile, STAT = STATUS)           
        if (STATUS /= SUCCESS_) stop "ConstructImpermeability - Geometry - ERR01"

            !Searches for the FacesOption
        call GetData(Me%Areas%ImPermeability,                                           &
                     ObjEnterData, iflag,                                               &
                     SearchType     = FromFile,                                         &
                     keyword        = 'FACES_IMPERMEABILITY',                           &
                     Default        = .false.,                                          &  
                     ClientModule   = 'ModuleGeometry',                                 &
                     STAT           = STATUS)

        if (Me%Areas%ImPermeability) then

            allocate (Me%Areas%Coef_U (Me%Size%KLB:Me%Size%KUB))
            allocate (Me%Areas%CoefX_U(Me%Size%KLB:Me%Size%KUB))

            allocate (Me%Areas%Coef_V (Me%Size%KLB:Me%Size%KUB))
            allocate (Me%Areas%CoefX_V(Me%Size%KLB:Me%Size%KUB))

            WorkKLB = Me%WorkSize%KLB
            WorkKUB = Me%WorkSize%KUB

            call GetData(Me%Areas%Coef_U(WorkKLB:WorkKUB),                              &
                         ObjEnterData, iflag,                                           &
                         SearchType     = FromFile,                                     &
                         keyword        = 'IMPER_COEF_U',                               &
                         Default        = 1.,                                           &  
                         ClientModule   ='ModuleGeometry',                              &
                         STAT           = STATUS)

            call GetData(Me%Areas%CoefX_U(WorkKLB:WorkKUB),                             &
                         ObjEnterData, iflag,                                           &
                         SearchType     = FromFile,                                     &
                         keyword        = 'IMPER_COEFX_U',                              &
                         Default        = 0.,                                           &  
                         ClientModule   ='ModuleGeometry',                              &
                         STAT           = STATUS)

            call GetData(Me%Areas%Coef_V(WorkKLB:WorkKUB),                              &
                         ObjEnterData, iflag,                                           &
                         SearchType     = FromFile,                                     &
                         keyword        = 'IMPER_COEF_V',                               &
                         Default        = 1.,                                           &  
                         ClientModule   ='ModuleGeometry',                              &
                         STAT           = STATUS)

            call GetData(Me%Areas%CoefX_V(WorkKLB:WorkKUB),                             &
                         ObjEnterData, iflag,                                           &
                         SearchType     = FromFile,                                     &
                         keyword        = 'IMPER_COEFX_V',                              &
                         Default        = 0.,                                           &  
                         ClientModule   ='ModuleGeometry',                              &
                         STAT           = STATUS)

        endif

        call KillEnterData(ObjEnterData, STAT = STATUS)

    end subroutine ConstructImpermeability

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

    subroutine AllocateVariables

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

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

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

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

        KLB = Me%Size%KLB
        KUB = Me%Size%KUB

        !Allocates T_Volumes
        allocate (Me%Volumes%VolumeZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR01'
        Me%Volumes%VolumeZ = FillValueDouble

        allocate (Me%Volumes%VolumeU(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR02'
        Me%Volumes%VolumeU = FillValueDouble

        allocate (Me%Volumes%VolumeV(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR03'
        Me%Volumes%VolumeV = FillValueDouble

        allocate (Me%Volumes%VolumeZOld(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR04'
        Me%Volumes%VolumeZOld = FillValueDouble

        !Allocate T_Areas
        allocate (Me%Areas%AreaU(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR05'
        Me%Areas%AreaU = FillValueReal

        allocate (Me%Areas%AreaV(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR06'
        Me%Areas%AreaV = FillValueReal

        !Allocate T_Distances
        allocate (Me%Distances%SZZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR07'
        Me%Distances%SZZ = FillValueReal

        allocate (Me%Distances%DZZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR08'
        Me%Distances%DZZ = FillValueReal

        allocate (Me%Distances%DWZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR09'
        Me%Distances%DWZ = FillValueReal

        allocate (Me%Distances%DUZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR09a'
        Me%Distances%DUZ = FillValueReal

        allocate (Me%Distances%DVZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR09b'
        Me%Distances%DVZ = FillValueReal

        allocate (Me%Distances%InitialSZZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR09c'
        Me%Distances%SZZ = FillValueReal

        allocate (Me%Distances%ZCellCenter(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR09d'
        Me%Distances%ZCellCenter = FillValueReal

        allocate (Me%Distances%DZI(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR09e'
        Me%Distances%DZI = FillValueReal

        allocate (Me%Distances%DZE(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR09f'
        Me%Distances%DZE = FillValueReal


        !Allocate T_KFloor
        allocate (Me%KFloor%Z(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR10'
        Me%KFloor%Z = FillValueInt

        allocate (Me%KFloor%U(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR11'
        Me%KFloor%U = FillValueInt

        allocate (Me%KFloor%V(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR12'
        Me%KFloor%V = FillValueInt

        allocate (Me%KFloor%Domain(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR13'
        Me%KFloor%Domain = FillValueInt

        allocate (Me%WaterColumn%Z(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR14'
        Me%WaterColumn%Z(:,:) = FillValueReal

        allocate (Me%WaterColumn%U(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR15'
        Me%WaterColumn%U(:,:) = FillValueReal

        allocate (Me%WaterColumn%V(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR16'
        Me%WaterColumn%V(:,:) = FillValueReal
        
        allocate (Me%KTop%Z(ILB:IUB, JLB:JUB), stat = STATUS)
        if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR17'
        Me%KTop%Z(:,:) = FillValueInt

    end subroutine AllocateVariables

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

    subroutine GetDomainsFromFile

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

        !Local-----------------------------------------------------------------
        integer                                     :: STATUS
        integer                                     :: FromBlock, FromFile
        logical                                     :: BlockFound
        logical                                     :: BlockLayersFound
        integer                                     :: iflag, NLayer
        integer                                     :: LBo, UBo
        integer                                     :: ClientNumber

        character(len = StringLength), parameter    :: block_begin = '<begindomain>'
        character(len = StringLength), parameter    :: block_end   = '<enddomain>'

        character(len = StringLength), parameter    :: beginlayers = '<<beginlayers>>'
        character(len = StringLength), parameter    :: endlayers   = '<<endlayers>>'


        character(len = StringLength)               :: DomainType
        real, dimension(:), allocatable             :: AuxVector
        type (T_Domain), pointer                    :: NewDomain
        integer                                     :: ObjEnterData = 0
        integer                                     :: i, ActualID, ID, LastLine, FirstLine, Line

        
        !Get Enter data parameter
        call GetExtractType(FromBlock = FromBlock, FromFile = FromFile)

        call ConstructEnterData(ObjEnterData, Me%InputFile, STAT = STATUS)           
        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR10"


        !Searches for the MinWaterColumn
        call GetData(Me%WaterColumn%Zmin, ObjEnterData, iflag,                          &
                     SearchType     = FromFile,                                         &
                     keyword        = 'MINIMUMDEPTH',                                   &
                     Default        =  0.1,                                             &  
                     ClientModule   = 'ModuleGeometry',                                 &
                     STAT           = STATUS)
        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR20"

        !Searches for the FacesOption
        call GetData(Me%FacesOption, ObjEnterData, iflag,                               &
                     SearchType     = FromFile,                                         &
                     keyword        = 'FACES_OPTION',                                   &
                     Default        = AverageTickness,                                  & 
                     ClientModule   ='ModuleGeometry',                                  &
                     STAT = STATUS)
        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR30"

        if (Me%FacesOption /= AverageTickness .and. Me%FacesOption /= MinTickness) then
            
            write(*,*) "The option FACES_OPTION : ", int(Me%FacesOption), " is not valid"
            write(*,*) "The only valid options are FACES_OPTION : ", AverageTickness, " advised option"
            write(*,*) "The only valid options are FACES_OPTION : ", MinTickness, " advanced user option"
            stop "GetDomainsFromFile - Geometry - ERR40"

        endif

        !The MinTickness option gives bad results espeially in the large scales.
        ! At the estuary scale bad results have also been identified. 
        !It is necessary in the future to understand the reason of this bad results.

        !Rewinds Buffer
        call RewindBuffer(ObjEnterData, STAT = STATUS)
        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR50"


        !Read the defined domains
        NLayer   = 0
        ActualID = 1
        BlockFound = .true.
        do while (BlockFound)

            call ExtractBlockFromBuffer(ObjEnterData, ClientNumber,                     &
                                        block_begin, block_end, BlockFound,             &
                                        STAT = STATUS)
            if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR60"
            
Block:      if (BlockFound) then

                !Searches for the ID of the domain
                call GetData(ID, ObjEnterData, iflag,                                   &
                             keyword        = 'ID',                                     &
                             SearchType     = FromBlock,                                &
                             ClientModule   = 'ModuleGeometry',                         &
                             STAT           = STATUS)
                if (STATUS /= SUCCESS_ .or. iflag == 0)                                 &
                    stop "GetDomainsFromFile - Geometry - ERR70"

CorretID:       if (ID == ActualID) then

                    allocate (NewDomain)
                    nullify  (NewDomain%Next)
                    nullify  (NewDomain%Prev)
                    NewDomain%ID = ActualID

                    !Searches for the type of the domain        
                    call GetData(DomainType, ObjEnterData, iflag,                       &
                                 keyword       = 'TYPE',                                &
                                 SearchType    = FromBlock,                             &
                                 ClientModule  = 'ModuleGeometry',                      &
                                 STAT          = STATUS)
                    if (STATUS /= SUCCESS_ .or. iflag == 0)                             &
                        stop "GetDomainsFromFile - Geometry - ERR80"

                    select case (trim(adjustl(DomainType)))

                    case ("FIXSPACING", "Fixspacing", "fixspacing")
                        NewDomain%DomainType = FixSpacing   
                    case ("SIGMA", "Sigma", "sigma")
                        NewDomain%DomainType = Sigma        
                    case ("LAGRANGIAN", "Lagrangian", "lagrangian")
                        NewDomain%DomainType = Lagrangian   
                    case ("CARTESIAN", "Cartesian", "cartesian")
                        NewDomain%DomainType = Cartesian
                    case ("HARMONIC", "Harmonic", "harmonic")
                        NewDomain%DomainType = Harmonic
                    case ("FIXSEDIMENT", "Fixsediment", "fixsediment")
                        NewDomain%DomainType = FixSediment
                    case ("SIGMATOP", "Sigmatop", "sigmatop")
                        NewDomain%DomainType = SigmaTop
                    case ("CARTESIANTOP", "Cartesiantop", "cartesiantop")
                        NewDomain%DomainType = CartesianTop
                    case default
                        stop "GetDomainsFromFile - Geometry - ERR90"
                    end select

                    !Searches for the number of Layers
                    call GetData(NewDomain%NumberOfLayers, ObjEnterData, iflag,         &
                                 keyword        = 'LAYERS',                             &
                                 SearchType     = FromBlock,                            &
                                 ClientModule   = 'ModuleGeometry',                     &
                                 STAT           = STATUS)
                    if (STATUS /= SUCCESS_ .or. iflag == 0)                             &
                        stop "GetDomainsFromFile - Geometry - ERR100"

                    NLayer = NLayer + NewDomain%NumberOfLayers

                    !Allocates LayerThickness
                    nullify  (NewDomain%LayerThickness)
                    LBo = NLayer - NewDomain%NumberOfLayers + 1                     !This way is usefull
                    UBo = NLayer                                                    !in the ComputeSZZ routines
                    allocate (NewDomain%LayerThickness(LBo:UBo), STAT = STATUS)
                    if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR110"
                    
                    ! Allows the definition of equidistant layers, layer thickness = constant
                    call GetData(NewDomain%Equidistant,                                 &
                                     ObjEnterData, iflag,                               &
                                     keyword      = 'EQUIDISTANT',                      &
                                     SearchType   = FromBlock,                          &
                                     ClientModule = 'ModuleGeometry',                   &
                                     DEFAULT      = 0. ,                                &
                                     STAT         = STATUS)
                    if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR120"

cd0:                if(iflag == 0) then

                        call ExtractBlockFromBlock(ObjEnterData, ClientNumber,          &
                                                   beginlayers, endlayers,              &
                                                   BlockLayersFound,                    &
                                                   FirstLine = FirstLine,               &
                                                   LastLine  = LastLine,                &
                                                   STAT      = STATUS)
                            
cd1 :                   if (STATUS .EQ. SUCCESS_      ) then    
cd2 :                       if (BlockLayersFound) then                        

                                if ((UBo - LBo + 1)/= (LastLine - FirstLine - 1))       &
                                    stop "GetDomainsFromFile - Geometry - ERR130"

                                Line = FirstLine + 1

                                do  i = UBo, LBo, -1

                                    call GetData(NewDomain%LayerThickness(i), ObjEnterData, &
                                                 iflag, Buffer_Line  = Line, STAT = STATUS)

                                    if (STATUS /= SUCCESS_ .or. iflag /= 1)             &
                                        stop "GetDomainsFromFile - Geometry - ERR140"

                                    Line = Line + 1

                                enddo


                            else

                                allocate (AuxVector(1:UBo-LBo+1), STAT = STATUS)

                                !Searches for the Layer Thickness
                                call GetData(AuxVector, ObjEnterData, iflag,            &
                                            SearchType   = FromBlock,                   &
                                            keyword      = 'LAYERTHICKNESS',            &
                                            ClientModule = 'ModuleGeometry',            &
                                            STAT = STATUS)
                                if (STATUS /= SUCCESS_ .or. iflag == 0)                 &
                                    stop "GetDomainsFromFile - Geometry - ERR150"

                                do i = LBo, UBo
                                    NewDomain%LayerThickness(i) = AuxVector(i-LBo+1)
                                enddo

                                deallocate (AuxVector)

                            endif cd2

                        else if (STATUS .EQ. BLOCK_END_ERR_) then cd1

                            stop "GetDomainsFromFile - Geometry - ERR160"

                        endif cd1

                    else  cd0
                     
                        do i = LBo, UBo
                            NewDomain%LayerThickness(i) = NewDomain%Equidistant
                        end do

                    endif cd0

                    !Searches for the Tolerance Depth (TYPE == SIGMA, ISOPYCNIC or Lagrangian)
                    if ((NewDomain%DomainType == Sigma)     .or.                        &
                        (NewDomain%DomainType == IsoPycnic) .or.                        &
                        (NewDomain%DomainType == Lagrangian)) then
                        
                        call GetData(NewDomain%ToleranceDepth, ObjEnterData, iflag,     &
                                     SearchType   = FromBlock,                          &  
                                     keyword      = 'TOLERANCEDEPTH',                   &
                                     ClientModule = 'ModuleGeometry',                   &
                                     Default      = 0.05,                               &
                                     STAT         = STATUS)
                        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR170"
                    endif

                    !Searches for Total Thickness / Bottom Surface
                    if (NewDomain%DomainType == FixSpacing  .or.                        &
                        NewDomain%DomainType == FixSediment .or.                        &
                        NewDomain%DomainType == SigmaTop) then

                        call GetData(NewDomain%TotalThickness, ObjEnterData, iflag,     &
                                     SearchType     = FromBlock,                        &
                                     keyword        = 'TOTALTHICKNESS',                 &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_ .or. iflag == 0)                         &
                            stop "GetDomainsFromFile - Geometry - ERR180"

                        call GetData(NewDomain%EmptyTopLayers, ObjEnterData, iflag,     &
                                     SearchType     = FromBlock,                        &
                                     keyword        = 'EMPTY_TOP_LAYERS',               &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_)                                         &
                            stop "GetDomainsFromFile - Geometry - ERR190"

                    else

                        if (NewDomain%DomainType /= CartesianTop) then
                            call GetData(NewDomain%DomainDepth, ObjEnterData, iflag,    &
                                         SearchType     = FromBlock,                    &
                                         keyword        = 'DOMAINDEPTH',                &
                                         ClientModule   = 'ModuleGeometry',             &
                                         STAT           = STATUS)
                            if (STATUS /= SUCCESS_ .or. iflag == 0)                     &
                                stop "GetDomainsFromFile - Geometry - ERR200"
                        endif

                    endif

                    !Searches for the coeficient which indicates how much a Lagrangian layer 
                    !can deform (MinEvolveLayerThickness)
                    if (NewDomain%DomainType == Lagrangian) then
                        
                        call GetData(NewDomain%MinEvolveLayerThickness,                 &
                                     ObjEnterData, iflag,                               &
                                     SearchType     = FromBlock,                        &
                                     keyword        = 'MINEVOLVELAYERTHICKNESS',        &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     Default        = 0.0,                              &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_ .or. iflag == 0)                         &
                            stop "GetDomainsFromFile - Geometry - ERR210"

                        call GetData(NewDomain%GridMovementDump,                        &
                                     ObjEnterData, iflag,                               &
                                     SearchType     = FromBlock,                        &
                                     keyword        = 'GRIDMOVEMENTDUMP',               &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     Default        = 0.0,                              &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_)                                         &
                            stop "GetDomainsFromFile - Geometry - ERR220"


                        call GetData(NewDomain%DisplacementLimit,                       &
                                     ObjEnterData, iflag,                               &
                                     keyword        = 'DISPLACEMENT_LIMIT',             &
                                     Default        = 1000.,                            &
                                     SearchType     = FromBlock,                        &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_)                                         &
                            stop "GetDomainsFromFile - Geometry - ERR230"

                       
                        call GetData(DomainType,                                        &
                                     ObjEnterData, iflag,                               &
                                     SearchType     = FromBlock,                        &
                                     keyword        = 'INITIALIZATION_METHOD ',         &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     Default        = "SIGMA",                          &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_)                                         &
                            stop "GetDomainsFromFile - Geometry - ERR240"


                        if     (DomainType .EQ. "SIGMA"     ) then
                            NewDomain%InitializationMethod = Sigma
                        elseif (DomainType .EQ. "CARTESIAN" ) then
                            NewDomain%InitializationMethod = Cartesian
                        else
                            write (*,*) "Initialization Method invalid"
                            stop "GetDomainsFromFile - Geometry - ERR250"
                        endif

                    endif


                   !Searches for the coeficient which indicates the minimal thickness in % 
                   !of the bottom cells if MinInitialLayerThickness = 1 - classic cartesian coordinates
                   !if < 1 - cartesian with shave cells
                   !In cartesian coordinates the reference layer thickness is compute for the 
                   !maximum depth
                    if (NewDomain%DomainType == Cartesian     .or.                      &
                        NewDomain%DomainType == Harmonic      .or.                      &
                        (NewDomain%DomainType == Lagrangian .and.                       &
                         NewDomain%InitializationMethod == Cartesian)) then

                        call GetData(NewDomain%MinInitialLayerThickness,                &
                                     ObjEnterData, iflag,                               &
                                     SearchType     = FromBlock,                        &
                                     keyword        = 'MININITIALLAYERTHICKNESS',       &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     Default        = 0.05,                             &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR270"
                    endif

                    if (NewDomain%DomainType == CartesianTop) then
                        NewDomain%MinInitialLayerThickness = 0.0
                    endif


                    !Seraches for the minimum thickness of colapsing cells of the Harmonic domain
                    if (NewDomain%DomainType == Harmonic) then

                        call GetData(NewDomain%MinEsp,                                  &
                                     ObjEnterData, iflag,                               &
                                     SearchType     = FromBlock,                        &
                                     keyword        = 'MIN_TOP_THICKNESS',              &
                                     ClientModule   = 'ModuleGeometry',                 &
                                     Default        = Me%WaterColumn%ZMin,              &
                                     STAT           = STATUS)
                        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR280"

                    endif

                    


                    !Inserts new domain into the domain list
                    !call AddDomainToGeometry(NewDomain)

                    call Add_Domain(NewDomain)

                    !Next ID to search for
                    ActualID = ActualID + 1

                    !Rewinds Buffer
                    call Block_Unlock(ObjEnterData, ClientNumber)
                    call RewindBuffer(ObjEnterData, STAT = STATUS)
                    if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR290"

                endif CorretID

            endif Block


        enddo 

        call Block_Unlock(ObjEnterData, ClientNumber)


        call KillEnterData(ObjEnterData, STAT = STATUS)
        if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR300"

    end subroutine GetDomainsFromFile

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

    subroutine ComputeLayers

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

        !Local-----------------------------------------------------------------
        type (T_Domain), pointer        :: CurrentDomain
        real                            :: RelativSum, Error, DomainDif, Sum,           &
                                           BottomDepth, TopDepth
        integer                         :: iLayer, LayersBelow
        !Begin-----------------------------------------------------------------

        !Computes limits for K
        Me%WorkSize%KUB = 0
        CurrentDomain => Me%FirstDomain
        do while (associated(CurrentDomain))
            Me%WorkSize%KUB = Me%WorkSize%KUB + CurrentDomain%NumberOfLayers
            CurrentDomain => CurrentDomain%Next
        enddo

        Me%WorkSize%KLB = 1
        Me%Size%KLB     = 0
        Me%Size%KUB     = Me%WorkSize%KUB + 1

        LayersBelow = 0        
        CurrentDomain => Me%FirstDomain
        do while (associated(CurrentDomain))

            !Computes UpperLayer and LowerLayer
            CurrentDomain%LowerLayer = LayersBelow + 1
            CurrentDomain%UpperLayer = LayersBelow + CurrentDomain%NumberOfLayers
            LayersBelow = LayersBelow + CurrentDomain%NumberOfLayers

            !Verifies and corrects relativ layer thickness
            !This is just done in the case of Sigma, Fixspacing, FixSediment and Lagrangian with Sigma
            !initialization
            if (CurrentDomain%DomainType  == Sigma              .or.                     &
                CurrentDomain%DomainType  == Fixspacing         .or.                     &
                CurrentDomain%DomainType  == FixSediment        .or.                     &
                CurrentDomain%DomainType  == SigmaTop           .or.                     &
                CurrentDomain%DomainType  == IsoPycnic          .or.                     &
                (CurrentDomain%DomainType == Lagrangian         .and.                    &
                 CurrentDomain%InitializationMethod == Sigma)) then

                RelativSum = 0.
                do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
                    RelativSum = RelativSum + CurrentDomain%LayerThickness(iLayer)
                enddo
                Error = abs(1.-RelativSum)
                if ( Error > 0.01) then

                    write(*,*)'Layers incorrectly defined - Domain number ',CurrentDomain%ID
                    stop 'ComputeLayers - Geometry - ERR10.'
!PCL
                else if (Error<=0.01 .and. Error > 1.e-6) then
                    do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
                        CurrentDomain%LayerThickness(iLayer) =                               &
                                                      CurrentDomain%LayerThickness(iLayer) - &
                                                      (1. - RelativSum) /                    &
                                                      real(CurrentDomain%NumberOfLayers)
                    enddo
                endif

            else if (CurrentDomain%DomainType  == Cartesian) then 

                Sum = 0.
                do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
                    Sum = Sum + CurrentDomain%LayerThickness(iLayer)
                enddo

                if (CurrentDomain%ID == Me%LastDomain%ID) then
                    TopDepth = 0.
                else
                    TopDepth = CurrentDomain%DomainDepth
                endif

                
                if (CurrentDomain%ID > 1) then
                    BottomDepth = CurrentDomain%Prev%DomainDepth
                else
                    call GetMaximumValue(Me%ObjTopography, BottomDepth)
                endif
                
                DomainDif = BottomDepth - TopDepth   

                !Error = abs(DomainDif - Sum)

                if (DomainDif > Sum) then
                    write(*,*)'Layers incorrectly defined - Domain number ',CurrentDomain%ID
                    stop 'ComputeLayers - Geometry - ERR20.'

                endif

            endif

            !Verifies if there is a fixspacing domain, other then the lower domain
            if (CurrentDomain%ID /= 1 .and. CurrentDomain%DomainType == FixSpacing) then
                write(*,*)'Just the first domain can be of the type FixSpacing'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR30.'
            endif

            !Verifies if there is a fixsediment domain, other then the lower domain
            if (CurrentDomain%ID /= 1 .and. CurrentDomain%DomainType == FixSediment) then
                write(*,*)'Just the first domain can be of the type FixSediment and it cant coexist with other domains'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR40.'
            endif

            !Verifies if a cartesian domain is above a fixspacing domain
            if (CurrentDomain%ID > 1 .and. CurrentDomain%DomainType == Cartesian .and.   &
                                 (Me%FirstDomain%DomainType == FixSpacing .or.           &
                                  Me%FirstDomain%DomainType == FixSediment)) then
                write(*,*)'A Cartesian domain type cant overlay a Fixspacing domain'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR50.'
            endif

            !Verifies if a harmonic domain is above a fixspacing domain
            if (CurrentDomain%ID > 1 .and. CurrentDomain%DomainType == Harmonic .and.    &
                                 (Me%FirstDomain%DomainType == FixSpacing .or.           &
                                  Me%FirstDomain%DomainType == FixSediment)) then
                write(*,*)'A Harmonic domain type cant overlay a Fixspacing or FixSediment domain'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR60.'
            endif

            !Verifies if a Lagrangian domain is above a fixspacing domain
            if (CurrentDomain%ID > 1 .and. CurrentDomain%DomainType == Lagrangian .and.  &
                                 (Me%FirstDomain%DomainType == FixSpacing .or.           &
                                  Me%FirstDomain%DomainType == FixSediment)) then
                write(*,*)'A Largragian domain type cant overlay a Fixspacing or FixSediment domain'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR70.'
            endif

            !Verifies if SigmaTop domain is topmost domain
            if (CurrentDomain%DomainType == SigmaTop .and.                               &
                CurrentDomain%ID         /= Me%LastDomain%ID) then
                write(*,*)'A SigmaTop domain type can only be the topmost domain'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR10.'
            endif

            !Verifies if there is a sigma domain below a SigmaTop domain (just one layer?)
            if (CurrentDomain%DomainType      == SigmaTop .and. CurrentDomain%ID == 1) then
                write(*,*)'A SigmaTop domain needs a Sigma domain below'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR10.'
            endif

            !Verifies if there is a sigma domain below a SigmaTop domain (sigma below)
            if (CurrentDomain%DomainType      == SigmaTop .and.                          &
                CurrentDomain%Prev%DomainType /= SIGMA) then
                write(*,*)'A SigmaTop domain needs a Sigma domain below'
                write(*,*)'Verify domain number :',CurrentDomain%ID
                stop 'ComputeLayers - Geometry - ERR11.'
            endif

            !Verifies if CartesianTop domain is topmost domain
            if (CurrentDomain%DomainType == CartesianTop) then
                if (CurrentDomain%ID /= Me%LastDomain%ID .or. CurrentDomain%ID /= 1) then
                    write(*,*)'A CartesianTop domain type can only be the topmost domain'
                    write(*,*)'A CartesianTop domain type can only exist alone'
                    write(*,*)'Verify domain number :',CurrentDomain%ID
                    stop 'ComputeLayers - Geometry - ERR11a.'
                endif
            endif



            !Verifies if the Minimal Layer thickness of a lagrangian domain is smalles then
            !50% of the initial layer thickness
            if (CurrentDomain%DomainType == Lagrangian) then
                    Error = 1.e10
                    do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
                        if (Error > CurrentDomain%LayerThickness(iLayer)) then
                            Error = CurrentDomain%LayerThickness(iLayer)
                        endif
                    enddo
                    if (CurrentDomain%MinEvolveLayerThickness > 0.50) then
                        write(*,*)'The MinimalThickness of the layers should not be greater then 50% of'
                        write(*,*)'of the initial layer thickness.'
                        write(*,*)'Verify domain number :',CurrentDomain%ID
                        write(*,*)'Layer                :',iLayer
                        stop 'ComputeLayers - Geometry - ERR80.'
                    else if (CurrentDomain%DisplacementLimit > 0.499 * Error) then
                        write(*,*)'The DisplacementLimit of the layers is higher than half the minimum LayerThickness'
                        stop 'ComputeLayers - Geometry - ERR90.'
                    endif
            endif

            !Verifies if the inital Layer thickness of the cartesian coordinates is not small then 
            !then 1% 
            if  (CurrentDomain%DomainType           == Cartesian .or.                 &
                 CurrentDomain%DomainType           == Harmonic  .or.                 &
                (CurrentDomain%DomainType           == Lagrangian .and.               &
                 CurrentDomain%InitializationMethod == Cartesian)) then

                    if (CurrentDomain%MinInitialLayerThickness < 0.01) then
                        write(*,*)'The MinimalThickness of the layers should not be smaller then 1% of'
                        write(*,*)'of the reference layer thickness'
                        write(*,*)'Verify domain number :',CurrentDomain%ID
                        write(*,*)'Layer                :',iLayer
                        stop 'ComputeLayers - Geometry - ERR100.'
                    endif
            endif


        CurrentDomain => CurrentDomain%Next
        enddo


    end subroutine ComputeLayers

    !--------------------------------------------------------------------------
    !In the case that the domain is Sigma, Isopycnic or Lagrangian the Bathymetry is tested in
    !order to avoid layers with zero height.
    !In the case of a cartesian domain, the Bathymetry is rearranged, in order to avoid layers to thin
    subroutine VerifyBathymetry (SurfaceElevation)

        !Arguments-------------------------------------------------------------
        real, dimension(:, :), pointer, optional    :: SurfaceElevation

        !Local-----------------------------------------------------------------
        real, dimension(:, :), pointer              :: Bathymetry, NewBathymetry, XX_IE, YY_IE
        real, dimension(:), pointer                 :: XX, YY
        real                                        :: GRID_ANGLE, Latitude, Longitude, XORIG, YORIG
        integer                                     :: ICOORD_TIP, Zone
                                                    
        integer                                     :: ILB, IUB, JLB, JUB
        integer                                     :: i, j, iLayer, STATUS
        real                                        :: DomainDepth, Tolerance
        real                                        :: TopDepth, BottomDepth, DomainThickness
        real                                        :: LayerTopDepth, LayerMinBottomDepth
        real                                        :: MinimalThickness, AllmostZero_ 
        real                                        :: ThicknessError 
        character(len=StringLength)                 :: BathymetryFile
        character(len=StringLength)                 :: Comment1, Comment2
        logical                                     :: WriteNewBathymetry = .false., Distortion
        type (T_Domain), pointer                    :: CurrentDomain
        type (T_Size2D)                             :: Size2D
        integer                                     :: LengthWithoutExt

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

        ILB = Me%WorkSize%ILB;  Size2D%ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB;  Size2D%IUB = Me%WorkSize%IUB
        JLB = Me%WorkSize%JLB;  Size2D%JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB;  Size2D%JUB = Me%WorkSize%JUB

        call GetGridData(Me%ObjTopography, Bathymetry, STAT = STATUS)
        if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR01'
        
        nullify (NewBathymetry)
        allocate(NewBathymetry(Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB))
        
        !Copy bathymetry array
        NewBathymetry(:,:) = Bathymetry(:,:)


        CurrentDomain => Me%FirstDomain
do1:    do while (associated(CurrentDomain))
            
cd1:        if ((CurrentDomain%DomainType == Sigma)      .or.                            &
                (CurrentDomain%DomainType == IsoPycnic)  .or.                            &
                (CurrentDomain%DomainType == Lagrangian .and.                            &
                 CurrentDomain%InitializationMethod == Sigma)) then

                if (CurrentDomain%ID == Me%LastDomain%ID) exit do1 !Dont correct the upper domain
                DomainDepth = CurrentDomain%DomainDepth
                Tolerance   = CurrentDomain%ToleranceDepth
                do j = JLB, JUB
                do i = ILB, IUB
                    if ((Bathymetry(i, j) >  DomainDepth) .and. &
                        (Bathymetry(i, j) <  DomainDepth + Tolerance)) then
                            write(*,*) NewBathymetry(i, j), DomainDepth, Tolerance
                            NewBathymetry(i, j) = DomainDepth
                            WriteNewBathymetry = .true.
                    endif
                enddo
                enddo

            !Correction in the case of cartesian domain
            else if (CurrentDomain%DomainType           == Cartesian .or.                &
                     CurrentDomain%DomainType           == Harmonic  .or.                &
                    (CurrentDomain%DomainType           == Lagrangian .and.              &
                     CurrentDomain%InitializationMethod == Cartesian)) then

                !Gets the upper and the lower depth of the domain and calculates its
                !thickness
                if (CurrentDomain%ID == Me%LastDomain%ID) then
                    TopDepth = 0.
                else
                    TopDepth = CurrentDomain%DomainDepth
                endif
                if (CurrentDomain%ID > 1) then
                    BottomDepth = CurrentDomain%Prev%DomainDepth
                else
                    call GetMaximumValue(Me%ObjTopography, BottomDepth)
                endif
                DomainThickness = BottomDepth - TopDepth

                !Verifies if the depth of the Bathymetry is between the upper and the lower
                !part of the layer, if so corrects the Bathymetry
                LayerTopDepth   = TopDepth
doLayer:        do iLayer = CurrentDomain%UpperLayer, CurrentDomain%LowerLayer, -1
                    
                    !MinimalThickness is now defined as a percentage of the initial thickness of each layer  
                    MinimalThickness = CurrentDomain%MinInitialLayerThickness

                    !The LayerThickness is now given in meters, so dont multiply by DomainTickness
                    !Frank - Jan 2001
                    !LayerMinBottomDepth = LayerTopDepth + CurrentDomain%LayerThickness(iLayer) * &
                    !                      DomainThickness * MinimalThickness
                    LayerMinBottomDepth = LayerTopDepth + CurrentDomain%LayerThickness(iLayer) * &
                                          MinimalThickness
                    
doj:                do j = JLB, JUB
doi:                do i = ILB, IUB
                        AllmostZero_ = AllmostZeroFraction * CurrentDomain%LayerThickness(iLayer)
                        ! Avoid round-off errors - Hernâni Sept 2004
                        ! if ((Bathymetry(i, j) > LayerTopDepth) .and. &
                        !     (Bathymetry(i, j) < LayerMinBottomDepth)) then
                        if ((Bathymetry(i, j) - LayerTopDepth > AllmostZero_) .and. &
                            (Bathymetry(i, j) - LayerMinBottomDepth < -AllmostZero_)) then
                            NewBathymetry(i, j) = LayerTopDepth
                            write(*,*)'i              = ', i
                            write(*,*)'j              = ', j
                            write(*,*)'Bathymetry     = ', Bathymetry(i, j)
                            write(*,*)'New Bathymetry = ',NewBathymetry(i, j)
                            WriteNewBathymetry = .true.
                                 
                        endif

                    enddo doi
                    enddo doj
                    
                    !
                    !The Layer thickness is now given in meters, so dont multiply by DomainThickness
                    !Frank Jan 2001
                    LayerTopDepth = LayerTopDepth + CurrentDomain%LayerThickness(iLayer)
                enddo doLayer

            elseif (CurrentDomain%DomainType == CartesianTop) then

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

                    TopDepth        = SurfaceElevation(i, j)
                    BottomDepth     = Bathymetry(i, j)
                    DomainThickness = TopDepth - BottomDepth

                    iLayer          = CurrentDomain%UpperLayer        
                    LayerTopDepth   = TopDepth
                    do while (iLayer >= CurrentDomain%LowerLayer)
                        LayerTopDepth = LayerTopDepth - CurrentDomain%LayerThickness(iLayer)
                        AllmostZero_ = AllmostZeroFraction * CurrentDomain%LayerThickness(iLayer)
                        if (abs(BottomDepth - LayerTopDepth) < 0.01 .and. abs(BottomDepth - LayerTopDepth) > AllmostZero_) then
                            NewBathymetry(i, j)       = LayerTopDepth
                            write(*,*)'i              = ', i
                            write(*,*)'j              = ', j
                            write(*,*)'Bathymetry     = ',Bathymetry(i, j)
                            write(*,*)'New Bathymetry = ',NewBathymetry(i, j)
                            WriteNewBathymetry = .true.
                            exit
                        endif
                        iLayer        = iLayer - 1
                    enddo


                enddo
                enddo
                

            elseif(CurrentDomain%DomainType == FixSediment)then

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

                    ThicknessError = abs(Bathymetry(i, j) - CurrentDomain%TotalThickness)
                    
                    if (ThicknessError > 1e-4 .and. Bathymetry(i, j) .lt. -99.) then
                        write(*,*) 'Inconsistent sediment domain file in cells(i,j)'
                        write(*,*) i, j
                        write(*,*) 'Domain total thickness = ', CurrentDomain%TotalThickness
                        write(*,*) 'Grid data file value   = ', NewBathymetry(i, j)
                        NewBathymetry(i, j) = CurrentDomain%TotalThickness
                        WriteNewBathymetry = .true.
                    endif
                enddo
                enddo

            endif cd1

            CurrentDomain => CurrentDomain%Next

        enddo do1

        if (WriteNewBathymetry) then
            
            !Gets XX and YY
            call GetHorizontalGrid(Me%ObjHorizontalGrid, XX = XX, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR02'

            call GetHorizontalGrid(Me%ObjHorizontalGrid, YY = YY, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR03'

            !Gets XX_IE and YY_IE
            call GetHorizontalGrid(Me%ObjHorizontalGrid, XX_IE = XX_IE, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR02a'

            call GetHorizontalGrid(Me%ObjHorizontalGrid, YY_IE = YY_IE, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR03a'

            !Gets the type of Coordinates
            call GetGridCoordType(Me%ObjHorizontalGrid, ICOORD_TIP, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR04'


            !Gets the Grid angle
            call GetGridAngle    (Me%ObjHorizontalGrid, GRID_ANGLE, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR05'

            !Gets the Latitude / Longitude
            call GetLatitudeLongitude(Me%ObjHorizontalGrid, Latitude, Longitude, STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR06'

            !Gets Zone of the Bathymetry
            call GetGridZone         (Me%ObjHorizontalGrid, Zone, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR07'

            !Gets Origin of the Bathymetry
            call GetGridOrigin      (Me%ObjHorizontalGrid, Xorig, Yorig, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR08'

            !Writes new Bathymetry
            call GetGridDataFileName(Me%ObjTopography, FileName = BathymetryFile, STAT= STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR09'

            BathymetryFile  = adjustl(BathymetryFile)
            Comment1        = "Automatic Generated Grid Data File"
            Comment2        = "Based On "//trim(BathymetryFile)
            LengthWithoutExt= len_trim(BathymetryFile) - 4
            BathymetryFile  = BathymetryFile(1:LengthWithoutExt)//"_"//".new"

            call GetCheckDistortion (Me%ObjHorizontalGrid, Distortion, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR09a'

            if (.not. Distortion) then

                call WriteGridData  (FileName       = BathymetryFile,       &
                                     XX             = XX,                   &
                                     YY             = YY,                   &
                                     COMENT1        = Comment1,             &
                                     COMENT2        = Comment2,             &
                                     WorkSize       = Size2D,               & 
                                     CoordType      = ICOORD_TIP,           &
                                     Xorig          = Xorig,                &
                                     Yorig          = Yorig,                &
                                     Zone           = Zone,                 &
                                     GRID_ANGLE     = GRID_ANGLE,           &
                                     Latitude       = Latitude,             &
                                     Longitude      = Longitude,            &
                                     GridData2D_Real= NewBathymetry,        &
                                     FillValue      = -99.,                 &
                                     Overwrite      = ON,                   &
                                     STAT           = STATUS)
                if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR10'

            else

                call WriteGridData  (FileName       = BathymetryFile,       &
                                     ConnectionX    = XX_IE,                &
                                     ConnectionY    = YY_IE,                &
                                     COMENT1        = "***",                &
                                     COMENT2        = "***",                &
                                     WorkSize       = Size2D,               & 
                                     CoordType      = ICOORD_TIP,           &
                                     Xorig          = Xorig,                &
                                     Yorig          = Yorig,                &
                                     Zone           = Zone,                 &
                                     GRID_ANGLE     = GRID_ANGLE,           &
                                     Latitude       = Latitude,             &
                                     Longitude      = Longitude,            &
                                     GridData2D_Real= NewBathymetry,        &
                                     FillValue      = -99.,                 &
                                     Overwrite      = ON,                   &
                                     STAT           = STATUS)
                if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR10a'
                
            endif

            deallocate(NewBathymetry)
            nullify   (NewBathymetry)

            !Displays Message to inform
            write(*,*)'A new Bathymetry has been created, which consists with the geometry'
            write(*,*)'Modify the file Nomfich.dat and Re-run the model'
            write(*,*)'New Bathymetry file : ', trim(BathymetryFile)
            stop
        endif

        !Disposes pointer to the Bathymetry
        call UngetGridData(Me%ObjTopography, Bathymetry, STAT = STATUS)
            if (STATUS /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR11'

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

    end subroutine VerifyBathymetry


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

    subroutine ConstructKFloor (SurfaceElevation)

        !Arguments-------------------------------------------------------------
        real, dimension(:, :), pointer, optional    :: SurfaceElevation

        !Local-----------------------------------------------------------------
        real, dimension(:, :), pointer              :: Bathymetry
        integer, dimension(:, :), pointer           :: WaterPoints2D
        integer                                     :: ILB, IUB, JLB, JUB
        integer                                     :: iLayer, i, j, LayersBelow
        real(8)                                     :: BottomDepth, TopDepth, LayerTopDepth
        real(8)                                     :: LayerBottomDepthMin, LayerBottomDepthMax 
        real(8)                                     :: MinimalThickness, DomainThickness
        real(8)                                     :: LayerThicknessMin, LayerThicknessMax
        logical                                     :: FoundKFloor
        type (T_Domain), pointer                    :: CurrentDomain
        integer, dimension(:,:), pointer            :: ExteriorFacesU, ExteriorFacesV
        integer                                     :: STAT_CALL
        logical                                     :: NeedToStop 
        real(8)                                     :: MaxDomainThickness
        real(8)                                     :: TotalLayerThickness, AuxDepth
        real                                        :: Aux4, AllmostZero_
                                                    
        !WorkSize
        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

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

        !Gets a pointer to the Bathymetry
        call GetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructKFloor - Geometry - ERR01'

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructKFloor - Geometry - ERR02'

        !Exterior Boundaries
        call GetExteriorBoundaryFaces(Me%ObjHorizontalMap,                  &
                                      BoundaryPointsFaceU = ExteriorFacesU, &
                                      BoundaryPointsFaceV = ExteriorFacesV, &
                                      STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR03'


        !Checks which is the domain close to the bottom
        if (Me%FirstDomain%DomainType == FixSpacing  .or.                   &
            Me%FirstDomain%DomainType == FixSediment .or.                   &
            Me%FirstDomain%DomainType == CartesianTop) then

            do j = JLB, JUB
            do i = ILB, IUB
    
                if (WaterPoints2D(i, j) == WaterPoint) then
                    Me%KFloor%Domain(i, j) = Me%FirstDomain%ID
                endif

            enddo
            enddo

        else
            
            do j = JLB, JUB
            do i = ILB, IUB
    
                if (WaterPoints2D(i, j) == WaterPoint) then

                    Me%KFloor%Domain(i, j) = Me%LastDomain%ID
                    CurrentDomain => Me%FirstDomain
                    do while (associated(CurrentDomain))
                        if (Bathymetry(i, j) >  CurrentDomain%DomainDepth) then
                            Me%KFloor%Domain(i, j) = CurrentDomain%ID
                            exit
                        endif
                        CurrentDomain => CurrentDomain%Next
                    enddo

                endif

            enddo
            enddo

        endif

         
        !Computes for KFloor%Z
        NeedToStop         = .false.
        MaxDomainThickness =  null_real
doj:    do j = JLB, JUB
doi:    do i = ILB, IUB
                
            if (WaterPoints2D(i, j) == WaterPoint) then
                
                LayersBelow = 0
                CurrentDomain => Me%FirstDomain

                do while (CurrentDomain%ID /= Me%KFloor%Domain(i, j))
                    LayersBelow = LayersBelow + CurrentDomain%NumberOfLayers
                    CurrentDomain => CurrentDomain%Next
                enddo

                
                !Set Current domain so that it points to the domain close to the bottom.
                CurrentDomain => Me%FirstDomain
                do while (CurrentDomain%ID /= Me%KFloor%Domain(i, j))
                    CurrentDomain => CurrentDomain%Next
                enddo

                !In the case of a sigma, lagrangian (init sigma) or fixspacing domain
                !KFloorZ is always equal to LayersBelow + 1
                if ((CurrentDomain%DomainType /= Cartesian   )     .and.                 &
                    (CurrentDomain%DomainType /= Harmonic    )     .and.                 &
                    (CurrentDomain%DomainType /= CartesianTop)     .and.                 &
                    (.not. (CurrentDomain%DomainType == Lagrangian .and.                 &
                            CurrentDomain%InitializationMethod == Cartesian))) then
                    Me%KFloor%Z(i, j) = 1 + LayersBelow
                else
                    !Domain Top Depth
                    if (CurrentDomain%ID == Me%LastDomain%ID) then
                        if (CurrentDomain%DomainType /= CartesianTop) then
                            TopDepth = 0.
                        else
                            TopDepth = dble(-1.0 * SurfaceElevation(i, j))
                        endif
                    else
                        TopDepth = dble(CurrentDomain%DomainDepth)
                    endif

                    !Domain Bottom Depth
                    if (CurrentDomain%ID > 1) then
                        BottomDepth = dble(CurrentDomain%Prev%DomainDepth)
                    else
                        if (CurrentDomain%DomainType /= CartesianTop) then 
                            call GetMaximumValue(Me%ObjTopography, Aux4)
                            BottomDepth = dble(Aux4)
                        else

                            !To minimize roundoff errors when BathymTopoFactor = 1
                            if (Me%ExternalVar%BathymTopoFactor/=1.0) then
                                BottomDepth = dble(Me%ExternalVar%BathymTopoFactor) * dble(Bathymetry(i, j))
                            else
                                BottomDepth = dble(Bathymetry(i, j))
                            endif
                        endif
                    endif
                    
                    !Domain Thickness
                    DomainThickness = BottomDepth - TopDepth

                    !Verifies if the depth of the Bathymetry is between the upper and the lower
                    !part of the layer
                    FoundKFloor = .false.
                    LayerTopDepth = TopDepth
                    iLayer = CurrentDomain%UpperLayer
                    do while (iLayer >= CurrentDomain%LowerLayer .and. .not. FoundKFloor)

                        !MinimalThickness is now defined as a percentage of the initial thickness of each layer  
                        MinimalThickness = dble(CurrentDomain%MinInitialLayerThickness)

                        
                        !In the case of Cartesian, Harmonic or Lagrangian(Cartesian) domain
                        !the layer thickness is given in meters
                        if (CurrentDomain%DomainType            == Cartesian    .or.     &
                            CurrentDomain%DomainType            == Harmonic     .or.     &
                            CurrentDomain%DomainType            == CartesianTop .or.     &
                            (CurrentDomain%DomainType           == Lagrangian .and.      &
                             CurrentDomain%InitializationMethod == Cartesian)) then

                            LayerThicknessMax   = dble(CurrentDomain%LayerThickness(iLayer))
                            LayerThicknessMin   = dble(CurrentDomain%LayerThickness(iLayer)) * &
                                                  MinimalThickness
                        else
                            LayerThicknessMax   = dble(CurrentDomain%LayerThickness(iLayer)) * DomainThickness
                            LayerThicknessMin   = dble(CurrentDomain%LayerThickness(iLayer)) * &
                                                  DomainThickness * MinimalThickness
                        endif

                        AllmostZero_ = AllmostZeroFraction * CurrentDomain%LayerThickness(iLayer)

                        LayerBottomDepthMin = LayerTopDepth + LayerThicknessMin - AllmostZero_
                        LayerBottomDepthMax = LayerTopDepth + LayerThicknessMax + AllmostZero_

                        !To minimize roundoff errors when BathymTopoFactor = 1
                        if (Me%ExternalVar%BathymTopoFactor/=1.0) then

                            AuxDepth = dble(Me%ExternalVar%BathymTopoFactor) * dble(Bathymetry(i, j))
                        
                        else

                            AuxDepth = dble(Bathymetry(i, j))

                        endif

                        if (LayerBottomDepthMin <= AuxDepth .and. LayerBottomDepthMax >= AuxDepth) then

                            Me%KFloor%Z(i, j) = iLayer
                            FoundKFloor = .true.

                        else

                            LayerTopDepth = LayerTopDepth +  LayerThicknessMax

                        endif
                        iLayer = iLayer - 1
                    enddo

                    !Points with Bathymetry > TopDepth
                    if (.not. FoundKFloor) then
                        if (CurrentDomain%DomainType /= CartesianTop) then
                            Me%KFloor%Z(i, j) = CurrentDomain%UpperLayer
                        else
                            MaxDomainThickness  = max(MaxDomainThickness, DomainThickness)
                            TotalLayerThickness = LayerTopDepth - TopDepth
                            NeedToStop = .true.
                        endif
                    endif
                
                endif

            endif

        enddo doi
        enddo doj

        if (NeedToStop) then
            write(*,*)'Domain Thickness Exceeds Layer Definition'
            write(*,*)'Maximum DomainThickness : ',MaxDomainThickness
            write(*,*)'Total Layer Thickness   : ',TotalLayerThickness
            stop 'ConstructKFloor - ModuleGeometry - ERR99'
        endif

        !Computes KFloor%U
        do j = JLB + 1, JUB
        do i = ILB, IUB

            !Water in (i,j) and in (i,j-1)? or current face is an exterior face
            if ( ((Me%KFloor%Z(i, j    ) > FillValueInt) .and.           &
                  (Me%KFloor%Z(i, j - 1) > FillValueInt)) .or.           &
                 (ExteriorFacesU(i,j) == 1) ) then
                    Me%KFloor%U(i, j) = max(Me%KFloor%Z(i, j), Me%KFloor%Z(i, j - 1))
            endif
        
        enddo
        enddo


        !Computes KFloor%V
        do j = JLB, JUB
        do i = ILB + 1, IUB

            !Water in (i,j) and in (i-1,j)? or current face is an exterior face
            if ( ((Me%KFloor%Z(i,     j) > FillValueInt) .and.           &
                  (Me%KFloor%Z(i - 1, j) > FillValueInt)) .or.           &
                 (ExteriorFacesV(i,j) == 1)) then
                    Me%KFloor%V(i, j) = max(Me%KFloor%Z(i, j), Me%KFloor%Z(i - 1, j))
                                                 

            endif
        
        enddo
        enddo

        call UnGetHorizontalMap(Me%ObjHorizontalMap, ExteriorFacesU, STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR04'
        
        call UnGetHorizontalMap(Me%ObjHorizontalMap, ExteriorFacesV, STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR05'


        !Disposes pointer to the Bathymetry
        call UngetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR06'


        !Disposes pointer to WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR07'

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

    end subroutine ConstructKFloor

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

    subroutine Add_Domain(NewDomain)

        !Arguments-------------------------------------------------------------
        type(T_Domain),           pointer     :: NewDomain

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

        ! Add to the WaterDomain List a new property
        if (.not.associated(Me%FirstDomain)) then
            Me%FirstDomain       => NewDomain
            Me%LastDomain        => NewDomain
        else
            NewDomain%Prev       => Me%LastDomain
            Me%LastDomain%Next   => NewDomain
            Me%LastDomain        => NewDomain
        end if 

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

    end subroutine Add_Domain 


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

    !MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MO 

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

    subroutine ComputeInitialGeometry(GeometryID, WaterPoints3D, SurfaceElevation,       &
                                      ActualTime, ContinuesCompute, NonHydrostatic,      &
                                      SZZ, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, dimension(:, :, :), pointer        :: WaterPoints3D
        real, dimension(:, :), pointer              :: SurfaceElevation
        type(T_Time),        optional               :: ActualTime
        logical, intent(in), optional               :: ContinuesCompute
        logical, intent(in), optional               :: NonHydrostatic
        real, dimension(:, :, :), pointer, optional :: SZZ
        integer, intent(out), optional              :: STAT

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

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

        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

            !Actualize the time
            if (present(ActualTime)) then
                Me%ActualTime = ActualTime
            endif

            if (present(NonHydrostatic)) then
                Me%ExternalVar%NonHydrostatic = NonHydrostatic
            else
                Me%ExternalVar%NonHydrostatic = .false.
            endif

            if (present(ContinuesCompute)) then
                Me%ExternalVar%ContinuesCompute = ContinuesCompute
            else
                Me%ExternalVar%ContinuesCompute = .false.
            endif

cd2 :       if (Me%ExternalVar%ContinuesCompute) then 

                !Computes the Distances
                call ComputeDistances   (WaterPoints3D)

                !Computes the Areas
                call ComputeAreas      

                !Computes the Volumes
                call ComputeVolumes     

                !It is necessary for the Soil model
                call ComputeZCellCenter 

            else cd2

                if (present(SZZ)) then    
                    Me%Distances%SZZ (:,:,:) = SZZ(:,:,:)
                else
                   !Constructs SZZ with the initial surface elevation
                    call ComputeSZZ         (SurfaceElevation, INITIALGEOMETRY)
                endif

                !Computes the Distances
                call ComputeDistances   (WaterPoints3D)

                !Computes the Volumes
                call ComputeVolumes     

                !Computes the Areas
                call ComputeAreas       

                !Stores VolumeZOLD
                !In this case VolumeZOld will be equal to VolumeZ
                call StoreVolumeZOld    

                !It is necessary for the Soil model
                call ComputeZCellCenter 

            end if cd2

            !Computes the WaterColumn
            call ComputeWaterColumn (SurfaceElevation)

            STAT_ = SUCCESS_

        else cd1              

            STAT_ = ready_

        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine ComputeInitialGeometry

    !--------------------------------------------------------------------------
    ! This subroutine computes the VerticalGeometry
    ! It should be called during the working cycle

    subroutine ComputeVerticalGeometry(GeometryID, WaterPoints3D, SurfaceElevation,      &
                                       ActualTime, VerticalVelocity, DT_Waterlevel,      &
                                       SZZ, DecayTime, KTop, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, dimension(:, :, :), pointer        :: WaterPoints3D
        real, dimension(:, :), pointer, optional    :: SurfaceElevation
        type (T_Time),            optional          :: ActualTime
        real, dimension(:, :, :), pointer, optional :: VerticalVelocity             !Gives the vertical variation
        real, intent(in), optional                  :: DT_Waterlevel                !for the lagragean coordinate
        real, dimension(:, :, :), pointer, optional :: SZZ, DecayTime
        integer, dimension(:, :), pointer, optional :: KTop
        integer, intent(out), optional              :: STAT

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

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


        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

            !Actualize the time
            if (present(ActualTime)) then
                Me%ActualTime = ActualTime
            endif


            if (present(SZZ)) then    
                Me%Distances%SZZ (:,:,:) = SZZ(:,:,:)
            else
               !Computes SZZ
                call ComputeSZZ(SurfaceElevation, TRANSIENTGEOMETRY, VerticalVelocity, DT_Waterlevel)
            endif  

            if(present(KTop))then
                Me%KTop%Z(:,:) = KTop(:,:)
            endif

            if(present(DecayTime)) then
                Me%ExternalVar%DecayTime => DecayTime
            endif

            !Stores VolumeZOld
            call StoreVolumeZOld

            !Computes Distances
            call ComputeDistances(WaterPoints3D)

            !Computes Volumes
            call ComputeVolumes

            !Computes Areas
            call ComputeAreas

            !It is necessary for the Soil model
            call ComputeZCellCenter

            !Computes the WaterColumn
            if (Me%FirstDomain%DomainType /= FixSediment) then
                call ComputeWaterColumn(SurfaceElevation)
            endif

            nullify(Me%Externalvar%DecayTime)


            STAT_ = SUCCESS_

        else cd1              

            STAT_ = ready_

        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine ComputeVerticalGeometry

    !--------------------------------------------------------------------------
    
    subroutine ComputeKTop (SedimentDomain)

        !Arguments-------------------------------------------------------------
        type(T_Domain), pointer             :: SedimentDomain

        !External--------------------------------------------------------------
        integer                             :: STAT_CALL

        !Local-----------------------------------------------------------------
        integer, dimension(:,:), pointer    :: WaterPoints2D
        integer                             :: ILB, IUB, JLB, JUB
        integer                             :: i, j

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

        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB
        JLB = Me%WorkSize%JLB
        JUB = Me%WorkSize%JUB

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeKTop - ModuleGeometry - ERR01'

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

            if (WaterPoints2D(i, j) == WaterPoint) then
                Me%KTop%Z(i, j) = Me%WorkSize%KUB - SedimentDomain%EmptyTopLayers
            endif

        enddo
        enddo

        !Disposes pointer to WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeKTop - ModuleGeometry - ERR02'

    end subroutine ComputeKTop

    !--------------------------------------------------------------------------
    !Computes the WaterColumn (Bathymetry + SurfaceElevation)

    subroutine ComputeWaterColumn(SurfaceElevation)

        !Arguments-------------------------------------------------------------
        real, dimension(:, :), pointer          :: SurfaceElevation

        !Local-----------------------------------------------------------------
        real, dimension(:,:,:), pointer         :: DUZ, DVZ
        real, dimension(:, : ), pointer         :: Bathymetry,   WaterColumnZ,           &
                                                   WaterColumnU, WaterColumnV

        integer, dimension(:, :), pointer       :: WaterPoints2D, KFloorU, KFloorV
        integer                                 :: ILB, IUB, JLB, JUB, KUB
        integer                                 :: i, j, k, STAT_CALL, kbottom

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

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

        KUB = Me%WorkSize%KUB


        WaterColumnZ => Me%WaterColumn%Z
        WaterColumnU => Me%WaterColumn%U
        WaterColumnV => Me%WaterColumn%V

        KFloorU      => Me%KFloor%U
        KFloorV      => Me%KFloor%V

        DUZ          => Me%Distances%DUZ
        DVZ          => Me%Distances%DVZ



        !Gets Bathymetry
        call GetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeWaterColumn - Geometry - ERR01'

        !Gets WaterPoints2D
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                                       &
            stop 'ComputeWaterColumn - Geometry - ERR02'


        !Computes WaterColumn
        do j = JLB, JUB
        do i = ILB, IUB
            if (WaterPoints2D(i, j) ==  WaterPoint)                                      &
                WaterColumnZ(i, j) = Bathymetry(i, j) + SurfaceElevation(i, j)
        enddo
        enddo

        !Computes WaterColumnU
        do j = JLB+1, JUB
        do i = ILB  , IUB

            if (WaterPoints2D(i, j - 1) == WaterPoint .and.                              &
                WaterPoints2D(i, j    ) == WaterPoint) then

                WaterColumnU (i, j) = 0.

                kbottom             = KFloorU(i, j)

                do k = kbottom, KUB
                
                    WaterColumnU (i, j) =  WaterColumnU (i, j) + DUZ(i, j, k)

                enddo

            endif

        enddo
        enddo


        !Computes WaterColumnV
        do j = JLB  , JUB
        do i = ILB+1, IUB

            if (WaterPoints2D(i - 1, j) == WaterPoint .and.                              &
                WaterPoints2D(i    , j) == WaterPoint) then

                WaterColumnV (i, j) = 0.

                kbottom             = KFloorV(i, j)

                do k = kbottom, KUB
                
                    WaterColumnV (i, j) =  WaterColumnV (i, j) + DVZ(i, j, k)

                enddo

            endif

        enddo
        enddo


        nullify(WaterColumnZ, WaterColumnU, WaterColumnV)
        nullify(KFloorU, KFloorV)
        nullify(DUZ, DVZ)

        !UnGets WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                                       &
            stop 'ComputeWaterColumn - Geometry - ERR03'


        !UnGets Bathymetry
        call UngetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                                       &
            stop 'ComputeWaterColumn - Geometry - ERR04'

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

     end subroutine ComputeWaterColumn

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


    !--------------------------------------------------------------------------
    !Computes the Distances (DWZ + DZZ + DUZ + DVZ + DZI + DZE)

    subroutine ComputeDistances(WaterPoints3D)

        !Arguments-------------------------------------------------------------
        integer, dimension(:, :, :), pointer    :: WaterPoints3D

        !Local-----------------------------------------------------------------
        real   , dimension(:, :, :), pointer    :: DWZ, DZZ, DUZ, DVZ, SZZ, DZI, DZE 
        real   , dimension(:, :   ), pointer    :: DUX, DVY
        integer                                 :: FacesOption
        integer                                 :: ILB, IUB, JLB, JUB, KLB, KUB
        integer                                 :: i, j, k, STAT_CALL
        !Begin-----------------------------------------------------------------

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

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

        KLB = Me%WorkSize%KLB
        KUB = Me%WorkSize%KUB

        DZZ => Me%Distances%DZZ
        DWZ => Me%Distances%DWZ
        DUZ => Me%Distances%DUZ
        DVZ => Me%Distances%DVZ
        SZZ => Me%Distances%SZZ
        DZI => Me%Distances%DZI
        DZE => Me%Distances%DZE

        FacesOption = Me%FacesOption

        !Computes DWZ
        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB

            if (WaterPoints3D(i, j, k) == WaterPoint) then

                DWZ(i, j, k) = SZZ(i, j, k-1) - SZZ(i, j, k  )

            endif
        enddo
        enddo
        enddo

cd1:    if (FacesOption == MinTickness) then

            !Computes DUZ
            do k = KLB  , KUB
            do j = JLB+1, JUB
            do i = ILB  , IUB

                if (WaterPoints3D(i, j - 1, k) == WaterPoint .and.                            &
                    WaterPoints3D(i, j    , k) == WaterPoint) then
                    DUZ(i, j, k) =  min(DWZ(i, j-1, k), DWZ(i, j, k)) 
                endif

            enddo
            enddo
            enddo

            !Computes DVZ
            do k = KLB  , KUB
            do j = JLB  , JUB
            do i = ILB+1, IUB

                if (WaterPoints3D(i - 1, j, k) == WaterPoint .and.                            &
                    WaterPoints3D(i    , j, k) == WaterPoint) then

                    DVZ(i, j, k) =  min(DWZ(i-1, j, k), DWZ(i, j, k)) 

                endif

            enddo
            enddo
            enddo


        else if (FacesOption == AverageTickness) then cd1

            !Gets DZX, DZY
            call GetHorizontalGrid(Me%ObjHorizontalGrid, DUX = DUX, DVY = DVY, & 
                                   STAT = STAT_CALL)

            if (STAT_CALL /= SUCCESS_)                                                  &
                stop 'ComputeAreas - Geometry - ERR01'


            !Computes DUZ
            do k = KLB  , KUB
            do j = JLB+1, JUB
            do i = ILB  , IUB

                if (WaterPoints3D(i, j - 1, k) == WaterPoint .and.                            &
                    WaterPoints3D(i, j    , k) == WaterPoint) then

                    DUZ(i, j, k) =  (DWZ(i, j, k)     * DUX(i, j - 1) +                  &
                                     DWZ(i, j - 1, k) * DUX(i, j))    /                  &
                                    (DUX(i, j - 1)    + DUX(i, j)) 
                
                endif

            enddo
            enddo
            enddo

            !Computes DVZ
            do k = KLB  , KUB
            do j = JLB  , JUB
            do i = ILB+1, IUB

                if (WaterPoints3D(i - 1, j, k) == WaterPoint .and.                            &
                    WaterPoints3D(i    , j, k) == WaterPoint) then

                    DVZ(i, j, k) =  (DWZ(i, j, k)     * DVY(i - 1, j) +                  &
                                     DWZ(i - 1, j, k) * DVY(i, j))    /                  &
                                    (DVY(i - 1, j)    + DVY(i, j))

                endif

            enddo
            enddo
            enddo

            !Nullifies auxilary pointers
            call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DUX, STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_)                                                   &
                stop 'ComputeDistances - Geometry - ERR02'


            call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DVY, STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_)                                                   &
                stop 'ComputeDistances - Geometry - ERR03'

        endif cd1

        !Computes DZZ
        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB

            if (WaterPoints3D(i, j, k) == WaterPoint) then
                DZZ(i, j, k) =  (DWZ(i, j, k+1) + DWZ(i, j, k)) / 2.0 
            endif

        enddo
        enddo
        enddo


        !Computes DZE
        do k = KLB  , KUB
        do j = JLB+1, JUB
        do i = ILB  , IUB

            if (WaterPoints3D(i, j - 1, k) == WaterPoint .and.                            &
                WaterPoints3D(i, j    , k) == WaterPoint) then

                DZE(i, j, k) =  (DUZ(i, j, k+1) + DUZ(i, j, k)) / 2.
            
            endif

        enddo
        enddo
        enddo

        !Computes DZI
        do k = KLB  , KUB
        do j = JLB  , JUB
        do i = ILB+1, IUB

            if (WaterPoints3D(i - 1, j, k) == WaterPoint .and.                            &
                WaterPoints3D(i    , j, k) == WaterPoint) then

                DZI(i, j, k) =  (DVZ(i, j, k+1) + DVZ(i, j, k)) / 2.
            
            endif

        enddo
        enddo
        enddo


        nullify(DWZ, DZZ, DUZ, DVZ, DZI, DZE, SZZ)

    end subroutine ComputeDistances

    !--------------------------------------------------------------------------
    !Computes the Volumes (VolumeZ, VolumeU, VolumeV)

    subroutine ComputeVolumes

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

        !Local-----------------------------------------------------------------
        real, dimension(:, :), pointer          :: DUX, DVY, DZX, DZY
        integer                                 :: ILB, IUB, JLB, JUB, KLB, KUB
        integer                                 :: i, j, k, STAT_CALL

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

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

        KLB = Me%WorkSize%KLB
        KUB = Me%WorkSize%KUB

        !Gets DUX, DVY
        call GetHorizontalGrid(Me%ObjHorizontalGrid, DUX = DUX, DVY = DVY,      &
                               DZX = DZX, DZY = DZY, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeVolumes - Geometry - ERR01'


        !Computes VolumeZ
        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB
            !Calculo dos volumes das celulas dos Z's como a area projectada * a altura media
            ! de cada prisma 
            Me%Volumes%VolumeZ(i, j ,k) = dble(Me%Distances%DWZ(i, j, k)) * &
                                          dble(DUX(i, j)) * dble(DVY(i, j))
        enddo
        enddo
        enddo

        !Computes VolumeU
        do k = KLB  , KUB
        do j = JLB+1, JUB
        do i = ILB  , IUB
            !Calculo dos volumes das celulas dos U's como a area projectada * a altura media
            Me%Volumes%VolumeU(i, j, k) = (Me%Volumes%VolumeZ(i, j-1, k)  + &
                                           Me%Volumes%VolumeZ(i, j, k  )) / &
                                           2.0 
        enddo
        enddo
        enddo

        !Computes VolumeV
        do k = KLB  , KUB
        do j = JLB  , JUB
        do i = ILB+1, IUB
            !Calculo dos volumes das celulas dos V's 
            Me%Volumes%VolumeV(i, j, k) = (Me%Volumes%VolumeZ(i-1, j, k)  + &
                                           Me%Volumes%VolumeZ(i, j, k  )) / &
                                           2.0
        enddo
        enddo
        enddo

        if (Me%ExternalVar%NonHydrostatic) then

            if (Me%Volumes%FirstVolW) then

                allocate (Me%Volumes%VolumeW(Me%Size%ILB:Me%Size%IUB, &
                                             Me%Size%JLB:Me%Size%JUB, &
                                             Me%Size%KLB:Me%Size%KUB),&
                                             stat = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'ComputeVolumes - Geometry - ERR42'
                Me%Volumes%VolumeW = FillValueDouble

                Me%Volumes%FirstVolW = .false.
            endif

            !Computes VolumeW
            do k = KLB + 1 , KUB
            do j = JLB     , JUB
            do i = ILB     , IUB
                !Calculo dos volumes das celulas dos W's 
                Me%Volumes%VolumeW(i, j, k) = (Me%Volumes%VolumeZ(i, j, k-1)  + &
                                               Me%Volumes%VolumeZ(i, j, k  )) / &
                                               2.0
            enddo
            enddo
            enddo

        endif



        !Nullifies auxilary pointers
        call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DUX, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeVolumes - Geometry - ERR02'


        call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DVY, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeVolumes - Geometry - ERR03'


        !Nullifies auxilary pointers
        call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DZX, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeVolumes - Geometry - ERR04'


        call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DZY, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeVolumes - Geometry - ERR05'

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

    end subroutine ComputeVolumes

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

    !--------------------------------------------------------------------------
    !Computes the Areas(AreaU, AreaV)

    subroutine ComputeAreas

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

        !Local-----------------------------------------------------------------
        real, dimension(:, :, :), pointer       :: DUZ, DVZ
        real, dimension(:, :),    pointer       :: DXX, DYY
        real                                    :: WaterLevel
        integer                                 :: ILB, IUB, JLB, JUB, KLB, KUB
        integer                                 :: i, j, k, STAT_CALL

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

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

        KLB = Me%WorkSize%KLB
        KUB = Me%WorkSize%KUB

        DUZ => Me%Distances%DUZ
        DVZ => Me%Distances%DVZ

        !Gets DZX, DZY
        call GetHorizontalGrid(Me%ObjHorizontalGrid, DXX = DXX, DYY = DYY, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)stop 'ComputeAreas - Geometry - ERR01'
            

        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB
            !Calculo das areas nas faces nos pontos de calculo U's (AREA_U(I,J,K))
!PCL
            Me%Areas%AreaU(i, j, k) = DUZ(i, j, k) * DYY (i, j)  

!PCL
            Me%Areas%AreaV(i, j, k) = DVZ(i, j, k) * DXX (i, j)

            if (Me%Areas%Impermeability) then

                WaterLevel = max(-Me%Distances%SZZ(i, j, KUB), -Me%Distances%SZZ(i, j+1, KUB))

                Me%Areas%AreaU(i, j, k) = Me%Areas%AreaU(i, j, k) *             &
                                         (Me%Areas%Coef_U(k) + WaterLevel * Me%Areas%CoefX_U(k)) 
                                                
                WaterLevel = max(-Me%Distances%SZZ(i, j, KUB), -Me%Distances%SZZ(i+1, j, KUB))

                Me%Areas%AreaV(i, j, k) = Me%Areas%AreaV(i, j, k) *             &
                                         (Me%Areas%Coef_V(k) + WaterLevel * Me%Areas%CoefX_V(k))
                                                   


            endif

        enddo
        enddo
        enddo

        !Nullifies auxilary pointers
        call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DXX, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)stop 'ComputeAreas - Geometry - ERR02'
            


        call UnGetHorizontalGrid(Me%ObjHorizontalGrid, DYY, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeAreas - Geometry - ERR03'
            
        
        !nullify local pointer
        nullify(DUZ, DVZ)

        !----------------------------------------------------------------------
    
    end subroutine ComputeAreas

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


    !--------------------------------------------------------------------------
    !Stores the VolumeZOld (VolumeZOLD = VolumeZ)

    subroutine StoreVolumeZOld

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

        !Local-----------------------------------------------------------------
        integer                                 :: ILB, IUB, JLB, JUB, KLB, KUB
        integer                                 :: i, j, k

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

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

        KLB = Me%WorkSize%KLB
        KUB = Me%WorkSize%KUB

        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB
            Me%Volumes%VolumeZOld(i, j ,k) = Me%Volumes%VolumeZ(i, j ,k) 
        enddo
        enddo
        enddo

    end subroutine StoreVolumeZOld


    !--------------------------------------------------------------------------
    !For every domain calls the respective computation rotine
    subroutine ComputeSZZ (SurfaceElevation, ComputionType, VerticalVelocity, DT_Waterlevel)

        !Parameter-------------------------------------------------------------
        real, dimension(:, :), pointer              :: SurfaceElevation
        integer                                     :: ComputionType
        real, dimension(:, :, :), optional, pointer :: VerticalVelocity
        real, intent(in), optional                  :: DT_Waterlevel

        !Esternal--------------------------------------------------------------

        integer :: STAT_CALL

        !Local-----------------------------------------------------------------
        real,    dimension(:, :), pointer           :: Bathymetry
        integer, dimension(:, :), pointer           :: WaterPoints2D
        integer                                     :: i, j, STAT, KLandFace
        logical                                     :: BathymEvolution
        type (T_Domain), pointer                    :: CurrentDomain

        !Gets Bathymetry
        call GetGridData(Me%ObjTopography, Bathymetry, STAT)
        if (STAT /= SUCCESS_) stop 'ComputeSZZ - Geometry - ERR01'

        call GetGridDataEvolution(Me%ObjTopography, BathymEvolution, STAT)
        if (STAT /= SUCCESS_) stop 'ComputeSZZ - Geometry - ERR01'

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeSZZ - Geometry - ERR02'



        if (ComputionType == INITIALGEOMETRY .or. BathymEvolution) then

            do j = Me%WorkSize%JLB, Me%WorkSize%JUB
            do i = Me%WorkSize%ILB, Me%WorkSize%IUB

                if (WaterPoints2D(i, j) == WaterPoint) then

                    KLandFace = Me%KFloor%Z(i, j) - 1

                    !To minimize roundoff errors when BathymTopoFactor = 1
                    if (Me%ExternalVar%BathymTopoFactor/=1.0) then

                        Me%Distances%SZZ(i, j, KLandFace) = Me%ExternalVar%BathymTopoFactor * Bathymetry(i, j)

                    else

                        Me%Distances%SZZ(i, j, KLandFace) = Bathymetry(i, j)

                    endif

                endif

            enddo
            enddo
        endif

        !When this rotine is called from the constructor, all domains are initialized, so the
        !variable ComputionType should be equal to INITIALGEOMETRY. In the case that this 
        !rotine is called from the modifier, the variable ComputionType should be equal to 
        !TRANSIENTGEOMETRY 

        CurrentDomain => Me%FirstDomain
        do while (associated(CurrentDomain))

            select case (CurrentDomain%DomainType)

                case (FixSpacing)
                    if (ComputionType == INITIALGEOMETRY)                                &
                        call ComputeFixSpacing  (CurrentDomain)

                case (FixSediment)
                    if (ComputionType == INITIALGEOMETRY) then
                        call ComputeKTop        (CurrentDomain)
                        call ComputeFixSpacing  (CurrentDomain)
                    else
                        !call ComputeFixSediment (CurrentDomain)
                    endif

                case (SigmaTop)

                    if (ComputionType == INITIALGEOMETRY) then
                        call ComputeSigma (SurfaceElevation, CurrentDomain)
                    else
                        !Do Nothing
                    endif
                    
                case (Sigma)
                    if ((ComputionType == INITIALGEOMETRY) .or.                          &
                        (CurrentDomain%ID == Me%LastDomain%ID))                 &
                        call ComputeSigma(SurfaceElevation, CurrentDomain)  
                                          

                case (Isopycnic)
!                    call ConstructIsopycnic(ObjGeometry, CurrentDomain)

                case (Lagrangian)
                    if (ComputionType == INITIALGEOMETRY) then
                        if (CurrentDomain%InitializationMethod == Sigma) then
                            call ComputeSigma(SurfaceElevation, CurrentDomain)
                        else
                            call ComputeCartesian(SurfaceElevation, CurrentDomain, ComputionType)
                        endif
                    else
                        call ComputeLagrangian(SurfaceElevation, VerticalVelocity,       &
                                               DT_Waterlevel, CurrentDomain)
                    endif

                case (Cartesian)
                    if ((ComputionType     == INITIALGEOMETRY) .or.             &
                        (CurrentDomain%ID  == Me%LastDomain%ID))                & 
                        call ComputeCartesian(SurfaceElevation, CurrentDomain, ComputionType)

                case (Harmonic)
                    call ComputeHarmonic(SurfaceElevation, CurrentDomain)

                case (CartesianTop)
                    if (ComputionType == INITIALGEOMETRY) then
                        call ComputeCartesian(SurfaceElevation, CurrentDomain, ComputionType)
                    else
                        !Do Nothing
                    endif

                case default           
                    write(*,*)'UNKNOWN_ Domain type - review data file'
                    stop 'ComputeSZZ - Geometry - ERR01.'

            end select

            CurrentDomain => CurrentDomain%Next

        enddo

        !Stores the initial SZZ
        if (ComputionType == INITIALGEOMETRY) then
            Me%Distances%InitialSZZ(:, :, :) = Me%Distances%SZZ(:, :, :)
        endif

        !Disposes pointer to the Bathymetry
        call UngetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeSZZ - Geometry - ERR03'

       !UnGets WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeSZZ - Geometry - ERR04'


    end subroutine ComputeSZZ

    
    !--------------------------------------------------------------------------
    !Computes SZZ for a FixSpacing Domain

    subroutine ComputeFixSpacing(Domain)

        !Parameter-------------------------------------------------------------
        type (T_Domain), pointer                :: Domain

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

        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

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

        do k = Domain%LowerLayer, Domain%UpperLayer

            LayerThickness = Domain%TotalThickness * Domain%LayerThickness(k)

            do j = JLB, JUB
            do i = ILB, IUB
                Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k-1) - LayerThickness
                                                     
            enddo
            enddo

        enddo

    end subroutine ComputeFixSpacing

    !--------------------------------------------------------------------------
    !Computes SZZ for a FixSediment Domain

    subroutine ComputeFixSediment (Domain, VerticalVelocity, DT)


        !Parameter-------------------------------------------------------------
        type (T_Domain), pointer                :: Domain
        real, dimension(:, :, :), pointer       :: VerticalVelocity
        real                                    :: DT

        !Local-----------------------------------------------------------------
        integer                                 :: i, j, k, ILB, IUB, JLB, JUB
        integer                                 :: STAT_CALL
        integer, dimension(:,:), pointer        :: WaterPoints2D
        real                                    :: Displacement
        real                                    :: IntegratedDisplacement

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeFixSediment - Geometry - ERR01'


        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

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

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

            if (WaterPoints2D(i, j) == WaterPoint) then

                IntegratedDisplacement = 0.
                do k = Domain%LowerLayer, Domain%UpperLayer 
                    Displacement            = -1. * VerticalVelocity(i, j, k+1) * DT
                    IntegratedDisplacement  = IntegratedDisplacement + Displacement
                    Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k) + &
                                                         IntegratedDisplacement
                enddo

            endif

        enddo
        enddo
        

        !UnGets WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeFixSediment - Geometry - ERR02'
        
    end subroutine ComputeFixSediment

    
    !--------------------------------------------------------------------------
    !Computes SZZ for a Cartesian Domain

    subroutine ComputeCartesian(SurfaceElevation, Domain, ComputionType)

        !Parameter-------------------------------------------------------------
        real, dimension(:, :), pointer          :: SurfaceElevation
        type (T_Domain), pointer                :: Domain
        integer                                 :: ComputionType

        !Local-----------------------------------------------------------------
        integer, dimension(:,:), pointer        :: WaterPoints2D
        integer                                 :: i, j, k, ILB, IUB, JLB, JUB
        integer                                 :: STAT_CALL, kbottom
        real                                    :: TopDepth
        real                                    :: LayerThickness 

        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

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

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeCartesian - Geometry - ERR01'


cd1:    if (ComputionType == INITIALGEOMETRY) then

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

                if (WaterPoints2D(i, j) == WaterPoint .and.                              &
                    Domain%UpperLayer >= Me%KFloor%Z(i, j)) then

                    !Domain Top Depth
                    if (Domain%ID == Me%LastDomain%ID) then
                        if (Domain%DomainType /= CartesianTop) then
                            TopDepth = 0.
                        else
                            TopDepth = -1. * SurfaceElevation(i, j)
                        endif
                    else
                        TopDepth = Domain%DomainDepth
                    endif

                    !Inits SZZ
                    Me%Distances%SZZ(i, j, Domain%UpperLayer) = TopDepth

                    kbottom = max(Domain%LowerLayer, Me%KFloor%Z(i, j))

                    do k = Domain%UpperLayer - 1, kbottom, -1

                        !The Layerthickness is now given in meters
                        !Frank - Jan 2001
                        LayerThickness = Domain%LayerThickness(k+1)
                        Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k+1) + LayerThickness

                   enddo

                endif

            enddo
            enddo

        endif  cd1

        if (Domain%UpperLayer == Me%WorkSize%KUB) then 

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

                if (WaterPoints2D(i, j) == WaterPoint)  then
                                
                      !In the cartesian domain the layer close to the surface have variable thickness
                      Me%Distances%SZZ(i, j, Me%WorkSize%KUB) = -1. * SurfaceElevation(i, j)

                endif

            enddo
            enddo

        endif


       !UnGets WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeCartesian - Geometry - ERR02'



    end subroutine ComputeCartesian
    

    !--------------------------------------------------------------------------
    !Computes SZZ for a Harmonic Domain

    subroutine ComputeHarmonic(SurfaceElevation, Domain)

        !Parameter-------------------------------------------------------------
        real, dimension(:, :), pointer          :: SurfaceElevation
        type (T_Domain), pointer                :: Domain

        !Local-----------------------------------------------------------------
        integer, dimension(:,:), pointer        :: WaterPoints2D
        integer                                 :: i, j, k, ILB, IUB, JLB, JUB, KUB
        integer                                 :: STAT_CALL, kbottom
        real                                    :: MinEsp, MeanLayerThickness, LayerThickness

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

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

        MinEsp = Domain%MinEsp

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeHarmonic - Geometry - ERR01'

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

            if (WaterPoints2D(i, j) == WaterPoint) then

                !Inits SZZ
                Me%Distances%SZZ(i, j, Domain%UpperLayer) = 0.0

                !Lowest layer in the domain
                kbottom = max(Domain%LowerLayer, Me%KFloor%Z(i, j))

                !Upper layer in the domain
                KUB     = Domain%UpperLayer

                !Calculates first SZZ approach
                do k = Domain%UpperLayer - 1, kbottom, -1

                    LayerThickness = Domain%LayerThickness(k+1)
                    Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k+1) + LayerThickness

                enddo

                !Sets SZZ(KUB) igual to waterlevel
                Me%Distances%SZZ(i, j, KUB)   = -1. * SurfaceElevation(i, j)

                !Calculates harmonic down
                k              = KUB - 1
                LayerThickness = Me%Distances%SZZ(i, j, k) - Me%Distances%SZZ(i, j, k+1)
                do while (LayerThickness < MinEsp .and. k >= kbottom)
                    Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k+1) + MinEsp
                    k              = k - 1
                    LayerThickness = Me%Distances%SZZ(i, j, k) - Me%Distances%SZZ(i, j, k+1)
                enddo

                !If the thickness of the cell close to the bottom is smaller then MinEsp
                !cells have to be redistributed. Two Options:
                !1. n Layers * MinEsp <= WaterColumn => Redistribute all cells equally
                !2. Increase Bottom Cell until every cell thickness > MinEsp
                if (Me%Distances%SZZ(i, j, kbottom-1) - Me%Distances%SZZ(i, j, kbottom) < MinEsp) then

                    if (Me%Distances%SZZ(i, j, kbottom-1)  - Me%Distances%SZZ(i, j, KUB) <= (KUB - kbottom + 1) * MinEsp) then
                        MeanLayerThickness = (Me%Distances%SZZ(i, j, kbottom-1)  -           &
                                              Me%Distances%SZZ(i, j, KUB)      ) /           &
                                              float(KUB - kbottom + 1)
                        do k = KUB - 1, kbottom, -1
                            Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k+1) + MeanLayerThickness
                        enddo
                    else
                        Me%Distances%SZZ(i, j, kbottom) = Me%Distances%SZZ(i, j, kbottom-1) - MinEsp
                        k = kbottom
                        do while (Me%Distances%SZZ(i, j, k) - Me%Distances%SZZ(i, j, k + 1) <  MinEsp)
                            Me%Distances%SZZ(i, j, k + 1) = Me%Distances%SZZ(i, j, k) - MinEsp 
                            k = k + 1
                        enddo
                    endif
                endif
            endif
        enddo
        enddo

        !UnGets WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeHarmonic - Geometry - ERR02'


    end subroutine ComputeHarmonic


    !--------------------------------------------------------------------------
    !Computes SZZ for a Sigma Domain

    subroutine ComputeSigma(SurfaceElevation, Domain)

        !Parameter-------------------------------------------------------------
        real, dimension(:, :), pointer          :: SurfaceElevation
        type (T_Domain), pointer                :: Domain

        !Local-----------------------------------------------------------------
        integer, dimension(:, :), pointer       :: WaterPoints2D
        integer                                 :: i, j, k, ILB, IUB, JLB, JUB
        integer                                 :: LowerLayer, UpperLayer
        integer                                 :: LowerLayerTwoDomain, UpperLayerTwoDomain
        integer                                 :: STAT_CALL
        real                                    :: TopDepth, BottomDepth, DomainThickness
        real                                    :: MinorThickness, LayerThickness
        real                                    :: MeanLayerThicknessBelow, TwoDomainThickness

        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

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

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeSigma - Geometry - ERR01'


        !Lower and Upper Layer of the Domain
        LowerLayer = Domain%LowerLayer
        UpperLayer = Domain%UpperLayer

doj:    do j = JLB, JUB
doi:    do i = ILB, IUB

            
cd4:        if (WaterPoints2D(i, j) == WaterPoint .and. &
                LowerLayer >= Me%KFloor%Z(i, j) ) then

                !Upper Limit (m)
                if (Domain%ID < Me%LastDomain%ID) then
                    if (associated(Domain%Next)) then
                        if (Domain%Next%DomainType == SigmaTop) then
                            TopDepth = -1.* SurfaceElevation(i, j) + Domain%Next%TotalThickness
                        else
                            TopDepth = Domain%DomainDepth
                        endif
                    else
                        TopDepth = Domain%DomainDepth
                    endif
                else
                    TopDepth = -1.* SurfaceElevation(i, j)
                endif

                !Lower Limit
                BottomDepth = Me%Distances%SZZ(i, j, LowerLayer-1)
        
                !DomainThickness 
                DomainThickness = BottomDepth - TopDepth

                !Frank
                ! The test below is now done in the hydrodynamic module. It can't be done here,
                ! due to the Volumebalance of the cell. 
                ! Paulo Chambel Leitão
                !OLD WAY - is correct because in this way the volume is never zero. 
                !if the volume is zero all the operations that devided the volume will have a overflow
!                if (DomainThickness <= Me%WaterColumn%ZMin) then
!                    write(*,*)'Low domain thickness', i, j
!                    DomainThickness = Me%WaterColumn%ZMin
!                endif

                ! Paulo Chambel Leitão
                !NEW WAY - is wrong            
                ! NOTE: If the Waterlevel in a certain cell is lower then the WaterColumn%ZMin, water
                ! can just enter to the cell. So it will be no "danger" (I hope) to calculate the
                ! volumes of the cell, since the waterlevel is positive. 
                ! If the WaterLevel is negative, the volumes and areas are set to zero, instead of 
                ! setting them to WaterColumn%ZMin (like until now), so the erro is smaller. 
                ! Frank May 99
                !if (DomainThickness < 0.) DomainThickness = 0.

                MinorThickness = 1.e5
                !Normal computation
                do k = LowerLayer, UpperLayer
                    LayerThickness = DomainThickness * Domain%LayerThickness(k)
                    Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k-1) - LayerThickness
                    if (LayerThickness < MinorThickness) MinorThickness = LayerThickness
                enddo

                !Redefinition of the current domain and the lower domain in the case that the lower
                !one is of the type FixSpacing and, at the same time (the DomainThickness is less 
                !or equal to the WaterColumn%ZMin) or (the Minor layer thichness is less or equal to
                !the mean layer thickness of the fixspacing domain). 
                !In the case this happens, the thickness of the two layers is added and then
                !evenly distributed
                !
                !The if's are separated, to avoid that the array indexes run out of bound
                !

cd1:            if (Domain%ID > 1) then
cd2:                if (Domain%Prev%DomainType == FixSpacing) then
                        MeanLayerThicknessBelow = Domain%Prev%TotalThickness /           &
                                                  float(Domain%Prev%NumberOfLayers)
                        
cd3:                    if ( (DomainThickness <= Me%WaterColumn%ZMin) .or.      &
                             (MinorThickness  <= MeanLayerThicknessBelow) ) then

                            !The total thickness of the two domains
                            LowerLayerTwoDomain = Domain%Prev%LowerLayer
                            UpperLayerTwoDomain = Domain%UpperLayer
                            TwoDomainThickness  = Me%Distances%SZZ(i, j, LowerLayerTwoDomain-1) &
                                                - Me%Distances%SZZ(i, j, UpperLayerTwoDomain) 

                            !Redefinition of the FixSpacing Domain
                            do k = Domain%Prev%LowerLayer, Domain%Prev%UpperLayer
                                LayerThickness = TwoDomainThickness * &
                                                 Domain%Prev%LayerThickness(k) / 2.
                                Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k-1) - LayerThickness
                            enddo

                            !Redefinition of the Sigma Domain
                            do k = Domain%LowerLayer, Domain%UpperLayer
                                LayerThickness = TwoDomainThickness * &
                                                 Domain%LayerThickness(k) / 2.
                                Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k-1) - LayerThickness
                            enddo

                        endif cd3
                    endif cd2   
                endif cd1

            endif cd4

        enddo doi
        enddo doj

       !UnGets WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                             &
            stop 'ComputeSigma - Geometry - ERR02'

    end subroutine ComputeSigma


    !--------------------------------------------------------------------------
    !Computes SZZ for a Lagrangian Domain
    subroutine ComputeLagrangian(SurfaceElevation, VerticalVelocity, DT_Waterlevel, Domain)

        !Parameter-------------------------------------------------------------
        real, dimension(:, :), pointer          :: SurfaceElevation
        real, dimension(:, :, :), pointer       :: VerticalVelocity
        real, intent(in)                        :: DT_Waterlevel
        type (T_Domain), pointer                :: Domain

        !Local-----------------------------------------------------------------
        integer                                 :: i, j, k, ILB, IUB, JLB, JUB
        integer                                 :: LowerLayer, UpperLayer
        integer, dimension(:, :), pointer       :: WaterPoints2D
        real                                    :: TopDepth, BottomDepth, DomainThickness
        real                                    :: MinimalThickness, GridMovementDump
        real                                    :: MaximumDisplacementUp, MaximumSzzUp, MaximumUndumpSzzUp
        real                                    :: MaximumDisplacementDown, MaximumSzzDown, MaximumUndumpSzzDown
        real                                    :: FreeSZZ, NewSZZ, MaximumDisplacement, FreeGridVelocity
        real                                    :: DisplacementLimit
        integer                                 :: STAT_CALL

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

        !Shorten variables name 

        ILB = Me%WorkSize%ILB
        IUB = Me%WorkSize%IUB

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

        !Dumping
        GridMovementDump = Domain%GridMovementDump

        !Upper Layer of the Domain
        UpperLayer       = Domain%UpperLayer

        !Gets a pointer to 2D WaterPoints
        call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'ComputeLagrangian - Geometry - ERR10'

doj:    do j = JLB, JUB
doi:    do i = ILB, IUB

cd0 :       if (WaterPoints2D(i, j) == WaterPoint) then

                !Upper Limit (m)
                if (Domain%ID < Me%LastDomain%ID) then
                    TopDepth = Domain%DomainDepth
                else
                    TopDepth = -1.* SurfaceElevation(i, j)
                endif

                !Lower Layer
                LowerLayer = max(Domain%LowerLayer, Me%KFloor%Z(i, j))

                !Lower Limit
                BottomDepth = Me%Distances%SZZ(i, j, LowerLayer-1)
            
                !DomainThickness 
                DomainThickness = BottomDepth - TopDepth
            
                !See ComputeSigma for more information about the next statement
                !This correction is now done in the hydrodynamic module - Frank
!                if (DomainThickness <= Me%WaterColumn%ZMin) &
!                    DomainThickness = Me%WaterColumn%ZMin

                !Resizes the layers, according to the vertical velocity
                !This process is done from the BOTTOM to the TOP

                !MinimalThickness is now defined as a percentage of the initial thickness of each layer  
                MinimalThickness   = Domain%MinEvolveLayerThickness


                !DisplacementLimit is the maximum displacement that the model allow cell 
                !faces to move vertically in meters
                DisplacementLimit  = Domain%DisplacementLimit

                do k = LowerLayer, UpperLayer-1

                    MaximumDisplacementUp   = (Me%Distances%InitialSZZ(i, j, k) -   &
                                               Me%Distances%InitialSZZ(i, j, k+1))* &
                                               (1.-MinimalThickness)/2.



                    MaximumDisplacementDown = (Me%Distances%InitialSZZ(i, j, k-1) - &
                                               Me%Distances%InitialSZZ(i, j, k)) *  &
                                              (1.-MinimalThickness)/2.

                    MaximumDisplacement     = min(MaximumDisplacementUp,                     &
                                                  MaximumDisplacementDown,                   &
                                                  DisplacementLimit)


                    MaximumSzzUp            = Me%Distances%InitialSZZ(i, j, k) -    &
                                              MaximumDisplacement

                    MaximumUndumpSzzUp      = Me%Distances%InitialSZZ(i, j, k) -    &
                                              MaximumDisplacement   *                        &
                                              (1.0 - GridMovementDump)

                    MaximumSzzDown          = Me%Distances%InitialSZZ(i, j, k) +    &
                                              MaximumDisplacement

                    MaximumUndumpSzzDown    = Me%Distances%InitialSZZ(i, j, k) +    &
                                              MaximumDisplacement     *                      &
                                              (1.0 - GridMovementDump)

                    !For the Layer k the Vertical Velocity to consider is the velocity in 
                    !k+1(Upper Face) once the vertical velocity has a different index than 
                    !the SZZ
                    FreeGridVelocity = VerticalVelocity(i, j, k+1)

                    
                    FreeSZZ = Me%Distances%SZZ(i, j, k) -                           &
                              FreeGridVelocity * DT_Waterlevel

                    if (FreeSZZ < MaximumUndumpSzzUp) then !Dumping to the maximum level

                        NewSZZ = MaximumUndumpSzzUp - (MaximumUndumpSzzUp - MaximumSzzUp) *  &
                                 (1. - 1./(MaximumUndumpSzzUp - FreeSZZ + 1))

                    elseif (FreeSZZ > MaximumUndumpSzzDown) then !Dumping to the minimum level

                        NewSZZ = MaximumUndumpSzzDown + (MaximumSzzDown -                    &
                                 MaximumUndumpSzzDown) * (1. - 1./(FreeSZZ -                 &
                                 MaximumUndumpSzzDown + 1.))
                    else

                        NewSZZ = FreeSZZ

                    endif

                    Me%Distances%SZZ(i, j, k) = NewSZZ

                    !Nudging
                    if (associated(Me%ExternalVar%DecayTime)) then
                        Me%Distances%SZZ(i, j, k) =  Me%Distances%SZZ(i, j, k) +        &
                                                    (Me%Distances%InitialSZZ(i, j, k) - &
                                                     Me%Distances%SZZ(i, j, k))       * &
                                                     DT_Waterlevel / Me%ExternalVar%DecayTime(i, j, k)
                                                    
                    endif

                enddo
                !Actualizes SZZ at the surface if the domain exists (DomainThickness>=0) 
                !in this I,J cell 
                if (DomainThickness>=0)                                                 &
                    Me%Distances%SZZ(i, j, UpperLayer) = TopDepth

            endif cd0

        enddo doi
        enddo doj

        !UnGets WaterPoints2D
        call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_)                                                      &
            stop 'ComputeLagrangian - Geometry - ERR50'


    end subroutine ComputeLagrangian


    !--------------------------------------------------------------------------
    !Reads SZZ and VolumeZOld

    subroutine ReadGeometryBin(GeometryID, Unit, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, intent(in)                         :: Unit
        integer, intent(out), optional              :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: ready_   
        integer                                     :: STAT_
        integer                                     :: i, j, k, ILB, IUB, JLB, JUB, KLB, KUB

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


        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

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

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

            KLB = Me%Size%KLB
            KUB = Me%Size%KUB

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                read(Unit) Me%Distances%SZZ(i, j, k)
            enddo
            enddo
            enddo

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                read(Unit) Me%Distances%InitialSZZ(i, j, k)
            enddo
            enddo
            enddo

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                read(Unit) Me%Distances%DWZ(i, j, k)
            enddo
            enddo
            enddo

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                read(Unit) Me%Volumes%VolumeZOld(i, j, k)
            enddo
            enddo
            enddo

            STAT_ = SUCCESS_

        else 

            STAT_ = ready_

        end if cd1


        if (present(STAT)) &
            STAT = STAT_

    end subroutine ReadGeometryBin


    !--------------------------------------------------------------------------
    !Reads SZZ and VolumeZOld

    subroutine WriteGeometryBin(GeometryID, Unit, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, intent(in)                         :: Unit
        integer, intent(out), optional              :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: ready_   
        integer                                     :: STAT_
        integer                                     :: i, j, k, ILB, IUB, JLB, JUB, KLB, KUB

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


        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

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

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

            KLB = Me%Size%KLB
            KUB = Me%Size%KUB

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                write(Unit) Me%Distances%SZZ(i, j, k)
            enddo
            enddo
            enddo

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                write(Unit) Me%Distances%InitialSZZ(i, j, k)
            enddo
            enddo
            enddo

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                write(Unit) Me%Distances%DWZ(i, j, k)
            enddo
            enddo
            enddo

            do k = KLB, KUB
            do j = JLB, JUB
            do i = ILB, IUB
                write(Unit)Me%Volumes%VolumeZOld(i, j, k)
            enddo
            enddo
            enddo

            STAT_ = SUCCESS_

        else 

            STAT_ = ready_

        end if cd1


        if (present(STAT)) &
            STAT = STAT_

    end subroutine WriteGeometryBin
    
    !----------------------------------------------------------------------------


        !--------------------------------------------------------------------------
    !Reads SZZ and VolumeZOld

    subroutine WriteGeometryHDF(GeometryID, ObjHDF5, WriteHDF, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, intent(in)                         :: ObjHDF5
        logical                                     :: WriteHDF
        integer, intent(out), optional              :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: ready_   
        integer                                     :: STAT_
        integer                                     :: ILB, IUB, JLB, JUB, KLB, KUB
        integer                                     :: STAT_CALL

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

        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then
            
            !WorkSize
            ILB = Me%WorkSize%ILB
            IUB = Me%WorkSize%IUB

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

            KLB = Me%WorkSize%KLB
            KUB = Me%WorkSize%KUB
            
            if (WriteHDF) then
            
                !Sets limits for next write operations
                call HDF5SetLimits(ObjHDF5, ILB, IUB, JLB, JUB, KLB-1, KUB, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR01'

                call HDF5WriteData(ObjHDF5,"/Geometry", "VerticalZ", "m",               &
                                   Array3D = Me%Distances%SZZ,                          &
                                   STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR02'

                call HDF5WriteData(ObjHDF5,"/Geometry", "InitialSZZ", "m",              &
                                   Array3D = Me%Distances%InitialSZZ,                   &
                                   STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR03'

                !Sets limits for next write operations
                call HDF5SetLimits(ObjHDF5, ILB, IUB, JLB, JUB, KLB, KUB, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR01'


                call HDF5WriteData(ObjHDF5,"/Geometry", "DWZ", "m",                     &
                                   Array3D = Me%Distances%DWZ,                          &
                                   STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR04'

                call HDF5WriteData(ObjHDF5,"/Geometry", "VolumeZ", "m3",                &
                                   Array3D = Me%Volumes%VolumeZOld,                     &
                                   STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR05'
                
                call HDF5WriteData(ObjHDF5,"/Geometry", "KTop", "-",                    &
                                   Array2D = Me%KTop%Z,                                 &
                                   STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR06'

                call HDF5FlushMemory (ObjHDF5, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'WriteGeometryHDF - Geometry - ERR07'

            end if


            STAT_ = SUCCESS_

        else 

            STAT_ = ready_

        end if cd1


        if (present(STAT)) &
            STAT = STAT_

    end subroutine WriteGeometryHDF
    
    !----------------------------------------------------------------------------

    subroutine ReadGeometryHDF(GeometryID, HDF5FileName, STAT)

        !Parameter-------------------------------------------------------------
        integer                                 :: GeometryID
        character(len=*)                        :: HDF5FileName
        integer, intent(out), optional          :: STAT

        !Local-----------------------------------------------------------------
        integer                                 :: ready_   
        integer                                 :: STAT_
        integer                                 :: ILB, IUB, JLB, JUB, KLB, KUB
        integer                                 :: STAT_CALL
        integer                                 :: ObjHDF5 = 0
        integer                                 :: HDF5_READ

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


        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

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

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

            KLB = Me%WorkSize%KLB
            KUB = Me%WorkSize%KUB
            
            !Gets File Access Code
            call GetHDF5FileAccess  (HDF5_READ = HDF5_READ)

            !Opens HDF File
            call ConstructHDF5      (ObjHDF5,                                            &
                                     trim(HDF5FileName),                                 &
                                     HDF5_READ, STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR01'
            
            !Sets limits for next read operations
            call HDF5SetLimits(ObjHDF5, ILB, IUB, JLB, JUB, KLB-1, KUB, STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR02'

            call HDF5ReadData(ObjHDF5,"/Geometry", "VerticalZ",                         &
                               Array3D = Me%Distances%SZZ,                              &
                               STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR03'

            call HDF5ReadData(ObjHDF5,"/Geometry", "InitialSZZ",                        &
                               Array3D = Me%Distances%InitialSZZ,                       &
                               STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR04'

            !Sets limits for next read operations
            call HDF5SetLimits(ObjHDF5, ILB, IUB, JLB, JUB, KLB, KUB, STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR05'

            call HDF5ReadData(ObjHDF5,"/Geometry", "DWZ",                               &
                               Array3D = Me%Distances%DWZ,                              &
                               STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR06'

            call HDF5ReadData(ObjHDF5,"/Geometry", "VolumeZ",                           &
                               Array3D = Me%Volumes%VolumeZOld,                         &
                               STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR07'


            call HDF5ReadData(ObjHDF5,"/Geometry", "KTop",                              &
                               Array2D = Me%KTop%Z,                                     &
                               STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR08'


            call KillHDF5 (ObjHDF5, STAT = STAT_CALL)
            if (STAT_CALL /= SUCCESS_) stop 'ReadGeometryHDF - Geometry - ERR09'

            STAT_ = SUCCESS_

        else 

            STAT_ = ready_

        end if cd1


        if (present(STAT)) &
            STAT = STAT_

    end subroutine ReadGeometryHDF
    
    !----------------------------------------------------------------------------



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

    subroutine ComputeZCellCenter

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

        !Local-------------------------------------------------------------------
        real, dimension(:, :, :), pointer       :: SZZ 
        real, dimension(:, :, :), pointer       :: ZCellCenter 
        integer :: I, J, K

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

        SZZ         => Me%Distances%SZZ
        ZCellCenter => Me%Distances%ZCellCenter

do1 :   do K = Me%WorkSize%KLB, Me%WorkSize%KUB
do2 :   do J = Me%WorkSize%JLB, Me%WorkSize%JUB
do3 :   do I = Me%WorkSize%ILB, Me%WorkSize%IUB

            ZCellCenter(I,J,K) = -1.0 * (SZZ(I,J,K-1) + SZZ(I,J,K)) / 2.0

        end do do3
        end do do2
        end do do1

        nullify(SZZ        )
        nullify(ZCellCenter)

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

    end subroutine ComputeZCellCenter

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

    !SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR 

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

    subroutine GetGeometrySize(GeometryID , Size, WorkSize, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        type (T_Size3D), optional                   :: Size
        type (T_Size3D), optional                   :: WorkSize
        integer, optional, intent(OUT)              :: STAT

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

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

        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

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

            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 GetGeometrySize
    
    !--------------------------------------------------------------------------

    subroutine GetGeometryMinWaterColumn(GeometryID, MinWaterColumn, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        real                                        :: MinWaterColumn
        integer, intent(out), optional              :: STAT

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

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

        STAT_ = UNKNOWN_

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

             MinWaterColumn = Me%WaterColumn%ZMin

             STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryMinWaterColumn

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

    subroutine GetGeometryDistances(GeometryID, SZZ, DZZ, DWZ, DUZ, DVZ, DZI, DZE,        &
                                    ZCellCenter, ActualTime, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        real, dimension(:,:,:), pointer, optional   :: SZZ, DZZ, DWZ, DUZ, DVZ 
        real, dimension(:,:,:), pointer, optional   :: DZI, DZE, ZCellCenter
        type(T_Time),                    optional   :: ActualTime
        integer, intent(OUT),            optional   :: STAT

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

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

        STAT_ = UNKNOWN_

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

            !Verifies if the distances are up to date
            if (present(ActualTime)) then
                if (ActualTime .ne. Me%ActualTime) &
                    stop 'GetGeometryDistances - Geometry - ERR01'
            endif

            !SZZ
            if (present(SZZ)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                SZZ => Me%Distances%SZZ
            endif

            !DZZ
            if (present(DZZ)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                DZZ => Me%Distances%DZZ
            endif

            !DWZ
            if (present(DWZ)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                DWZ => Me%Distances%DWZ
            endif


            !DUZ
            if (present(DUZ)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                DUZ => Me%Distances%DUZ
            endif

            !DVZ
            if (present(DVZ)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                DVZ => Me%Distances%DVZ
            endif

            !DZI
            if (present(DZI)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                DZI => Me%Distances%DZI
            endif

            !DZE
            if (present(DZE)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                DZE => Me%Distances%DZE
            endif


            !ZCellCenter
            if (present(ZCellCenter)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                ZCellCenter => Me%Distances%ZCellCenter
            endif

            STAT_ = SUCCESS_

        else cd1

            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryDistances

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

    subroutine GetGeometryAreas(GeometryID, AreaU, AreaV, ActualTime, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        real, dimension(:, :, :), pointer, optional :: AreaU, AreaV
        type (T_Time), optional                     :: ActualTime
        integer, intent(out), optional              :: STAT

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

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


        STAT_ = UNKNOWN_

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

            !Verifies if the areas are up to date
            if (present(ActualTime)) then
                if (ActualTime .ne. Me%ActualTime) &
                    stop 'GetGeometryAreas - Geometry - ERR01'
            endif

            !AreaU
            if (present(AreaU)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                AreaU => Me%Areas%AreaU
            endif

            !AreaV
            if (present(AreaV)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                AreaV => Me%Areas%AreaV
            endif

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryAreas

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

    subroutine GetGeometryVolumes(GeometryID, VolumeZ, VolumeU, VolumeV, &
                                  VolumeW, VolumeZOld, ActualTime, STAT)

        !Parameter-------------------------------------------------------------
        integer                                         :: GeometryID
        real(8), dimension(:, :, :), pointer, optional  :: VolumeZ, VolumeU, VolumeV, VolumeW, VolumeZOld
        type (T_Time), optional                         :: ActualTime
        integer, intent(out), optional                  :: STAT

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

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

        STAT_ = UNKNOWN_

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

            !Verifies if the volumes are up to date
            if (present(ActualTime)) then
                if (ActualTime .ne. Me%ActualTime) &
                    stop 'GetGeometryVolumes - Geometry - ERR01'
            endif

            !VolumeZ
            if (present(VolumeZ)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                VolumeZ => Me%Volumes%VolumeZ
            endif

            !VolumeU
            if (present(VolumeU)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                VolumeU => Me%Volumes%VolumeU
            endif

            !VolumeV
            if (present(VolumeV)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                VolumeV => Me%Volumes%VolumeV
            endif

            !VolumeW
            if (present(VolumeW)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                VolumeW => Me%Volumes%VolumeW
            endif

            !VolumeZ
            if (present(VolumeZOld)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                VolumeZOld => Me%Volumes%VolumeZOld
            endif

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryVolumes


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

    subroutine GetGeometryKTop(GeometryID, KTopZ, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, dimension(:, :), pointer, optional :: KTopZ
        integer, intent(out), optional              :: STAT

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

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

        STAT_ = UNKNOWN_

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

            !KTopZ
            if (present(KTopZ)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                KTopZ => Me%KTop%Z
            endif

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryKTop
    
    !--------------------------------------------------------------------------

    subroutine GetGeometryKFloor(GeometryID, Z, U, V, Domain, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, dimension(:, :), pointer, optional :: Z, U, V, Domain
        integer, intent(out), optional              :: STAT

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

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

        STAT_ = UNKNOWN_

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

            !Z
            if (present(Z)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                Z => Me%KFloor%Z
            endif

            !U
            if (present(U)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                U => Me%KFloor%U
            endif

            !V
            if (present(V)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                V => Me%KFloor%V
            endif

            !Domain
            if (present(Domain)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                Domain => Me%KFloor%Domain
            endif

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryKFloor

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

    subroutine GetGeometryWaterColumn(GeometryID, WaterColumn , WaterColumnU,           &
                                                   WaterColumnV, ActualTime, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        real, dimension(:, :), pointer, optional    :: WaterColumn
        real, dimension(:, :), pointer, optional    :: WaterColumnU, WaterColumnV
        type (T_Time), optional                     :: ActualTime
        integer, intent(out), optional              :: STAT

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

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

        STAT_ = UNKNOWN_

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

            !Verifies if the areas are up to date
            if (present(ActualTime)) then
                if (ActualTime .ne. Me%ActualTime) &
                    stop 'GetGeometryWaterColumn - Geometry - ERR01'
            endif

            if (present(WaterColumn)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                WaterColumn => Me%WaterColumn%Z
            endif

            if (present(WaterColumnU)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                WaterColumnU => Me%WaterColumn%U
            endif

            if (present(WaterColumnV)) then
                call Read_Lock(mGEOMETRY_, Me%InstanceID)
                WaterColumnV => Me%WaterColumn%V
            endif

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryWaterColumn

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

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

    subroutine GetGeometryInputFile(GeometryID, InputFile, STAT)

        !Parameter-------------------------------------------------------------
        integer                                     :: GeometryID
        character(len=*)                            :: InputFile
        integer, intent(out), optional              :: STAT

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

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

        STAT_ = UNKNOWN_

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

            InputFile = Me%InputFile

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

    end subroutine GetGeometryInputFile

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

    subroutine UnGetGeometry2Dinteger(GeometryID, Array, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, pointer, dimension(:,:)            :: Array
        integer, optional, intent(OUT)              :: STAT

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

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

        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. READ_LOCK_ERR_) then

            nullify(Array)
            call Read_UnLock(mGEOMETRY_, Me%InstanceID, "UnGetGeometry2Dinteger")

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT))                                                    &
            STAT = STAT_

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

    end subroutine UnGetGeometry2Dinteger

    !--------------------------------------------------------------------------
    
    subroutine UnGetGeometry2Dreal(GeometryID, Array, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        real, pointer, dimension(:,:)               :: Array
        integer, optional, intent(OUT)              :: STAT



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

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

        STAT_ = UNKNOWN_


        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. READ_LOCK_ERR_) then

            nullify(Array)
            call Read_UnLock(mGEOMETRY_, Me%InstanceID, "UnGetGeometry2Dreal")

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

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

    end subroutine UnGetGeometry2Dreal

    !--------------------------------------------------------------------------
    
    subroutine UnGetGeometry3Dreal4(GeometryID, Array, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        real(4), pointer, dimension(:,:,:)          :: Array
        integer, optional, intent(OUT)              :: STAT

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

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

        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. READ_LOCK_ERR_) then

            nullify(Array)
            call Read_Unlock(mGEOMETRY_, Me%InstanceID, "UnGetGeometry3Dreal4")

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

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

    end subroutine UnGetGeometry3Dreal4

    !--------------------------------------------------------------------------
    
    subroutine UnGetGeometry3Dreal8(GeometryID, Array, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        real(8), pointer, dimension(:,:,:)          :: Array
        integer, optional, intent(OUT)              :: STAT

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

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

        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .EQ. READ_LOCK_ERR_) then

            nullify(Array)
            call Read_UnLock(mGEOMETRY_, Me%InstanceID, "UnGetGeometry3Dreal8")

            STAT_ = SUCCESS_
        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_

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

    end subroutine UnGetGeometry3Dreal8


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

    !DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCT 

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

    subroutine KillGeometry(GeometryID, STAT)

        !Arguments-------------------------------------------------------------
        integer                                     :: GeometryID
        integer, intent(out), optional              :: STAT

        !Local-----------------------------------------------------------------
        integer                                     :: ready_   
        integer                                     :: STAT_
        type (T_Domain), pointer                    :: CurrentDomain, DomainToDelete
        integer                                     :: nUsers

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

        STAT_ = UNKNOWN_

        call Ready(GeometryID, ready_)    

cd1 :   if (ready_ .NE. OFF_ERR_) then

            nUsers = DeassociateInstance(mGEOMETRY_,  Me%InstanceID)

            if (nUsers == 0) then

                !Deassociates External Instances
                nUsers = DeassociateInstance (mGRIDDATA_,       Me%ObjTopography)
                if (nUsers == 0) stop 'KillGeometry - Geometry - ERR10'

                nUsers = DeassociateInstance (mHORIZONTALMAP_,  Me%ObjHorizontalMap)
                if (nUsers == 0) stop 'KillGeometry - Geometry - ERR20'

                nUsers = DeassociateInstance (mHORIZONTALGRID_, Me%ObjHorizontalGrid)
                if (nUsers == 0) stop 'KillGeometry - Geometry - ERR30'

                !Deallocates local variables
                call DeallocateVariables

                !Delete all domains
                CurrentDomain => Me%FirstDomain
                do while (associated(CurrentDomain))
                    DomainToDelete => CurrentDomain
                    CurrentDomain  => CurrentDomain%Next
                    call DeleteDomain(DomainToDelete)
                enddo

                if (Me%Areas%Impermeability) then

                    deallocate (Me%Areas%Coef_U )
                    deallocate (Me%Areas%CoefX_U)

                    deallocate (Me%Areas%Coef_V )
                    deallocate (Me%Areas%CoefX_V)

                endif

                call DeallocateInstance

                GeometryID = 0
                STAT_      = SUCCESS_

            end if
        else 
            STAT_ = ready_
        end if cd1


        if (present(STAT)) STAT = STAT_
              
        !----------------------------------------------------------------------

    end subroutine KillGeometry

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

    subroutine DeallocateInstance 

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

        !Local-----------------------------------------------------------------
        type (T_Geometry), pointer          :: AuxObjGeometry
        type (T_Geometry), pointer          :: PreviousObjGeometry

        !Updates pointers
        if (Me%InstanceID == FirstGeometry%InstanceID) then
            FirstGeometry => FirstGeometry%Next
        else
            PreviousObjGeometry => FirstGeometry
            AuxObjGeometry      => FirstGeometry%Next
            do while (AuxObjGeometry%InstanceID /= Me%InstanceID)
                PreviousObjGeometry => AuxObjGeometry
                AuxObjGeometry      => AuxObjGeometry%Next
            enddo

            !Now update linked list
            PreviousObjGeometry%Next => AuxObjGeometry%Next

        endif
            
        !Deallocates instance
        deallocate (Me)
        nullify    (Me) 

    end subroutine DeallocateInstance

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

    subroutine DeallocateVariables

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

        !Local-----------------------------------------------------------------
        integer                         :: STATUS

        !Deallocates T_Volumes
        if (associated(Me%Volumes%VolumeZ)) then
            deallocate (Me%Volumes%VolumeZ, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR01'
        endif

        if (associated(Me%Volumes%VolumeU)) then
            deallocate (Me%Volumes%VolumeU, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR02'
        endif


        if (associated(Me%Volumes%VolumeV)) then
            deallocate (Me%Volumes%VolumeV, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR03'
        endif

        if (associated(Me%Volumes%VolumeW)) then
            deallocate (Me%Volumes%VolumeW, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR03a'
        endif


        if (associated(Me%Volumes%VolumeZOld)) then
            deallocate (Me%Volumes%VolumeZOld, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR04'
        endif

        !Allocate T_Areas
        if (associated(Me%Areas%AreaU)) then
            deallocate (Me%Areas%AreaU, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR05'
        endif

        if (associated(Me%Areas%AreaV)) then
            deallocate (Me%Areas%AreaV, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR06'
        endif

        !Allocate T_Distances
        if (associated(Me%Distances%SZZ)) then
            deallocate (Me%Distances%SZZ, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR07'
        endif

        if (associated(Me%Distances%DZZ)) then
            deallocate (Me%Distances%DZZ, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR08'
        endif

        if (associated(Me%Distances%DWZ)) then
            deallocate (Me%Distances%DWZ, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR09'
        endif

        if (associated(Me%Distances%DUZ)) then
            deallocate (Me%Distances%DUZ, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR09a'
        endif

        if (associated(Me%Distances%DVZ)) then
            deallocate (Me%Distances%DVZ, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR09b'
        endif

        if (associated(Me%Distances%DZI)) then
            deallocate (Me%Distances%DZI, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR09c'
        endif

        if (associated(Me%Distances%DZE)) then
            deallocate (Me%Distances%DZE, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR09d'
        endif


        if (associated(Me%Distances%InitialSZZ)) then
            deallocate (Me%Distances%InitialSZZ, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR09e'
        endif


        if (associated(Me%Distances%ZCellCenter)) then
            deallocate (Me%Distances%ZCellCenter, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR09f'
        endif


        !deallocate T_KFloor
        if (associated(Me%KFloor%Z)) then
            deallocate (Me%KFloor%Z, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR10'
        endif

        if (associated(Me%KFloor%U)) then
            deallocate (Me%KFloor%U, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR11'
        endif

        if (associated(Me%KFloor%V)) then
            deallocate (Me%KFloor%V, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR12'
        endif

        if (associated(Me%KFloor%Domain)) then
            deallocate (Me%KFloor%Domain, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR13'
        endif

        !deallocate T_KTop
        if (associated(Me%KTop%Z)) then
            deallocate (Me%KTop%Z, stat = STATUS)
            if (STATUS /= SUCCESS_) stop 'DeallocateVariables - Geometry - ERR20'
        endif
        

    end subroutine DeallocateVariables

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

    subroutine DeleteDomain (DomainToDelete)

        !Arguments-------------------------------------------------------------
        type (T_Domain), pointer                    :: DomainToDelete

        !Local-----------------------------------------------------------------
        type (T_Domain), pointer                    :: CurrentDomain
        type (T_Domain), pointer                    :: Previous
        type (T_Domain), pointer                    :: Next


        CurrentDomain => Me%FirstDomain
        do while (associated(CurrentDomain))
            if (CurrentDomain%ID == DomainToDelete%ID) then

                Previous       => CurrentDomain%Prev
                Next           => CurrentDomain%Next

                !Updates foward pointer
                if (associated(CurrentDomain%Prev)) then
                    Previous%Next  => Next
                endif

                !Updates backward pointer
                if (associated(CurrentDomain%Next)) then
                    Next%Prev      => Previous
                endif

                !Updates first domain
                if (DomainToDelete%ID == Me%FirstDomain%ID) then
                    Me%FirstDomain => Next
                endif

                !Updates last domain
                if (DomainToDelete%ID == Me%LastDomain%ID) then
                    Me%LastDomain => Previous
                endif

                !Deletes domain
                if (associated(DomainToDelete%LayerThickness)) then
                    deallocate(DomainToDelete%LayerThickness)
                    nullify   (DomainToDelete%LayerThickness)
                endif

                nullify   (DomainToDelete%Next)
                nullify   (DomainToDelete%Prev)
                deallocate(DomainToDelete)
                nullify   (DomainToDelete)

                return

            endif

            CurrentDomain => CurrentDomain%Next
        enddo

    end subroutine DeleteDomain



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

    !MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEME

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


    subroutine Ready (ObjGeometry_ID, ready_) 

        !Arguments-------------------------------------------------------------
        integer                                     :: ObjGeometry_ID
        integer                                     :: ready_

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

        nullify (Me)

cd1:    if (ObjGeometry_ID > 0) then
            call LocateObjGeometry(ObjGeometry_ID)
            ready_ = VerifyReadLock (mGEOMETRY_, Me%InstanceID)
        else
            ready_ = OFF_ERR_
        end if cd1

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

    end subroutine Ready

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

    subroutine LocateObjGeometry (ObjGeometryID)

        !Arguments-------------------------------------------------------------
        integer                                     :: ObjGeometryID

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

        Me => FirstGeometry
        do while (associated (Me))
            if (Me%InstanceID == ObjGeometryID) exit
            Me => Me%Next
        enddo

        if (.not. associated(Me)) stop 'Geometry - LocateObjGeometry - ERR01'

    end subroutine LocateObjGeometry

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

end module ModuleGeometry



