!------------------------------------------------------------------------------
!        IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE         : Mohid Model
! PROJECT       : Mohid Base 2
! MODULE        : Interface 
! URL           : http://www.mohid.com
! AFFILIATION   : IST/MARETEC, Marine Modelling Group
! DATE          : May 2003
! REVISION      : Luis Fernandes - v4.0
! DESCRIPTION   : Module to serve as link to zero-dimensional sinks and sources modules
!
!------------------------------------------------------------------------------

Module ModuleInterface

    use ModuleGlobalData
    use ModuleTime
    use ModuleWaterQuality
    use ModuleSedimentQuality
    use ModuleCEQUALW2
    use ModuleLife
    use ModuleBenthos
    use ModuleEnterData, only: ReadFileName

    implicit none

    private

    !Subroutines & Functions---------------------------------------------------

    !Constructor
    public  :: ConstructInterface
    private ::      AllocateInstance
    private ::      ReadInterfaceFilesName
    private ::      StartSinksSourcesModel
    private ::      AllocateVariables                                                
    private ::      Check_Options
    private ::          FindProperty

    !Modifier 
    public  :: Modify_Interface
    private ::      FillMassTempSalinity
    private ::      InputData
    private ::      PropertyIndexNumber
    private ::      UnfoldMatrix

    !Selector
    public  :: GetRateFlux

    !Destructor 
    public  :: KillInterface
    private ::      DeallocateInstance

    !Management
    private ::      Ready
    private ::          LocateObjInterface

    !Interfaces
    private :: ConstructInterface1D
    private :: ConstructInterface2D
    private :: ConstructInterface3D
    interface  ConstructInterface
        module procedure ConstructInterface1D
        module procedure ConstructInterface2D
        module procedure ConstructInterface3D
    end interface  ConstructInterface


    private :: Modify_Interface1D
    private :: Modify_Interface2D
    private :: Modify_Interface3D
    interface  Modify_Interface
        module procedure Modify_Interface1D
        module procedure Modify_Interface2D
        module procedure Modify_Interface3D
    end interface  Modify_Interface

    
    private :: FillMassTempSalinity1D
    private :: FillMassTempSalinity2D
    private :: FillMassTempSalinity3D
    interface  FillMassTempSalinity
        module procedure FillMassTempSalinity1D
        module procedure FillMassTempSalinity2D
        module procedure FillMassTempSalinity3D
    end interface  FillMassTempSalinity

    private :: InputData1D
    private :: InputData2D
    private :: InputData3D
    interface  InputData
        module procedure InputData1D
        module procedure InputData2D
        module procedure InputData3D
    end interface  InputData

    private :: GetRateFlux2D
    private :: GetRateFlux3D
    interface  GetRateFlux
        module procedure GetRateFlux2D
        module procedure GetRateFlux3D
    end interface  GetRateFlux
    
    private :: UnfoldMatrix1D_I
    private :: UnfoldMatrix1D_R
    private :: UnfoldMatrix2D_I
    private :: UnfoldMatrix2D_R
    private :: UnfoldMatrix3D_I
    private :: UnfoldMatrix3D_R
    interface  UnfoldMatrix
        module procedure UnfoldMatrix1D_I
        module procedure UnfoldMatrix1D_R
        module procedure UnfoldMatrix2D_I
        module procedure UnfoldMatrix2D_R
        module procedure UnfoldMatrix3D_I
        module procedure UnfoldMatrix3D_R
    end interface  UnfoldMatrix


    !Type----------------------------------------------------------------------
    type       T_External
        type(T_Time     )                       :: Now
        integer, dimension(:      ), pointer    :: RiverPoints1D
        integer, dimension(:, :   ), pointer    :: WaterPoints2D
        integer, dimension(:, :, :), pointer    :: WaterPoints3D
        integer, dimension(:      ), pointer    :: OpenPoints1D
        integer, dimension(:, :   ), pointer    :: OpenPoints2D
        integer, dimension(:, :, :), pointer    :: OpenPoints3D
        real,    dimension(:      ), pointer    :: DWZ1D
        real,    dimension(:, :, :), pointer    :: DWZ
        logical                                 :: Vertical1D = .false.
    end type T_External

    type      T_Interface
        private
        integer                                 :: InstanceID
        character(StringLength)                 :: FileName
        character(StringLength)                 :: SinksSourcesModel
        type(T_Size3D)                          :: Size3D
        type(T_Size2D)                          :: Size2D
        type(T_Size1D)                          :: Size1D
        type(T_Size1D)                          :: Array
        type(T_Size1D)                          :: Prop
        type(T_External)                        :: ExternalVar
        real,    pointer, dimension(:,:  )      :: Mass
        real,    pointer, dimension(:,:  )      :: ConcentrationIncrement
        real,    pointer, dimension(:    )      :: Salinity
        real,    pointer, dimension(:    )      :: Alkalinity
        real,    pointer, dimension(:    )      :: Temperature
        real,    pointer, dimension(:    )      :: Oxygen
        real,    pointer, dimension(:    )      :: ShortWaveRadiation
        real,    pointer, dimension(:    )      :: LightExtCoefField
        real,    pointer, dimension(:    )      :: Thickness
        real,    pointer, dimension(:    )      :: FishFood
        real,    pointer, dimension(:    )      :: WaterPercentage
        real,    pointer, dimension(:    )      :: DissolvedToParticulate
        integer, pointer, dimension(:    )      :: OpenPoints
        logical, pointer, dimension(:    )      :: AddedProperties

        !Instance of ModuleTime
        integer                                 :: ObjTime              = 0

        !Instance of ModuleWaterQuality         
        integer                                 :: ObjWaterQuality      = 0

        !Instance of ModuleSedimentQuality      
        integer                                 :: ObjSedimentQuality   = 0

        !Instance of ModuleCEQUALW2
        integer                                 :: ObjCEQUALW2          = 0

        !Instance of ModuleBenthos
        integer                                 :: ObjBenthos           = 0

        !Instance of ModuleLife
        integer                                 :: ObjLife              = 0

        !Collection of instances                
        type(T_Interface),          pointer     :: Next        

    end type T_Interface

        
    !Global Module Variables
    type (T_Interface), pointer                 :: FirstObjInterface
    type (T_Interface), pointer                 :: Me

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

    contains

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

    !CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CO

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

    subroutine ConstructInterface3D(InterfaceID,                           &
                                    TimeID,                                &
                                    SinksSourcesModel,                     &
                                    DT,PropertiesList,                     &
                                    WaterPoints3D,                         &
                                    Size3D,                                &
                                    Vertical1D,                            &
                                    STAT)

        !Arguments-------------------------------------------------------------
        integer                                                 :: InterfaceID
        integer                                                 :: TimeID
        character(len=StringLength)                             :: SinksSourcesModel
        integer, dimension(:    ), pointer                      :: PropertiesList
        real,   intent (OUT)                                    :: DT
        integer,                     dimension(:,:,:), pointer  :: WaterPoints3D
        type(T_Size3D)                                          :: Size3D
        logical,intent (IN),  optional                          :: Vertical1D
        integer,intent (OUT), optional                          :: STAT     

        !External--------------------------------------------------------------
        integer                                                 :: ready_

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

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

        STAT_ = UNKNOWN_

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

        call Ready(InterfaceID, ready_)    

cd0 :   if (ready_ .EQ. OFF_ERR_) then
            
            !Allocates a new Instance
            call AllocateInstance

            !Associates External Instances
            Me%ObjTime                      =  AssociateInstance (mTIME_, TimeID)
            Me%SinksSourcesModel            =  SinksSourcesModel
            Me%ExternalVar%WaterPoints3D    => WaterPoints3D
            Me%Size3D                       =  Size3D 
            Me%Array%ILB                    = 1
            Me%Array%IUB                    = sum(Me%ExternalVar%WaterPoints3D)

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

            !Path to data file
            call ReadInterfaceFilesName

            !Start sinks and sources model
            call StartSinksSourcesModel(DT)

            !Verify compute options
            call Check_Options(PropertiesList)

            !Allocate variables global to whole module
            call AllocateVariables

            STAT_ = SUCCESS_
            
            !Returns ID
            InterfaceID = Me%InstanceID

            STAT_ = SUCCESS_

        else cd0
            
            stop 'ModuleInterface - ConstructInterface3D - ERR01' 

        end if cd0

        if (present(STAT)) STAT = STAT_

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

    end subroutine ConstructInterface3D

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

    subroutine ConstructInterface2D(InterfaceID,                           &
                                    TimeID,                                &
                                    SinksSourcesModel,                     &
                                    DT,PropertiesList,                     &
                                    WaterPoints2D,                         &
                                    Size2D,                                &
                                    STAT)

        !Arguments-------------------------------------------------------------
        integer                                                 :: InterfaceID
        integer                                                 :: TimeID
        character(len=StringLength)                             :: SinksSourcesModel
        integer, dimension(:    ), pointer                      :: PropertiesList
        real,   intent (OUT)                                    :: DT
        integer,                     dimension(:,:  ), pointer  :: WaterPoints2D
        type(T_Size2D)                                          :: Size2D
        integer,intent (OUT), optional                          :: STAT     

        !External--------------------------------------------------------------
        integer                                                 :: ready_

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

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

        STAT_ = UNKNOWN_

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


        call Ready(InterfaceID, ready_)    

cd0 :   if (ready_ .EQ. OFF_ERR_) then
            
            !Allocates a new Instance
            call AllocateInstance

            !Associates External Instances
            Me%ObjTime                      =  AssociateInstance (mTIME_, TimeID)
            Me%SinksSourcesModel            =  SinksSourcesModel
            Me%ExternalVar%WaterPoints2D    => WaterPoints2D
            Me%Size2D                       =  Size2D
            Me%Array%ILB                    = 1
            Me%Array%IUB                    = sum(Me%ExternalVar%WaterPoints2D)

            !Path to data file
            call ReadInterfaceFilesName

            !Start sinks and sources model
            call StartSinksSourcesModel(DT)

            !Verify compute options
            call Check_Options(PropertiesList)

            !Allocate variables global to whole module
            call AllocateVariables

            STAT_ = SUCCESS_
            
            !Returns ID
            InterfaceID = Me%InstanceID

            STAT_ = SUCCESS_

        else cd0
            
            stop 'ModuleInterface - ConstructInterface2D - ERR01' 

        end if cd0

        if (present(STAT)) STAT = STAT_

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

    end subroutine ConstructInterface2D

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

    subroutine ConstructInterface1D(InterfaceID,                           &
                                    TimeID,                                &
                                    SinksSourcesModel,                     &
                                    DT,PropertiesList,                     &
                                    RiverPoints1D,                         &
                                    Size1D,                                &
                                    STAT)

        !Arguments-------------------------------------------------------------
        integer                                                 :: InterfaceID
        integer                                                 :: TimeID
        character(len=StringLength)                             :: SinksSourcesModel
        integer, dimension(:    ), pointer                      :: PropertiesList
        real,   intent (OUT)                                    :: DT
        integer, dimension(:), pointer                          :: RiverPoints1D
        type(T_Size1D)                                          :: Size1D
        integer,intent (OUT), optional                          :: STAT     

        !External--------------------------------------------------------------
        integer                                                 :: ready_

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

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

        STAT_ = UNKNOWN_

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

        call Ready(InterfaceID, ready_)    

cd0 :   if (ready_ .EQ. OFF_ERR_) then
            
            !Allocates a new Instance
            call AllocateInstance

            !Associates External Instances
            Me%ObjTime                      =  AssociateInstance (mTIME_, TimeID)
            Me%SinksSourcesModel            =  SinksSourcesModel
            Me%ExternalVar%RiverPoints1D    => RiverPoints1D
            Me%Size1D                       =  Size1D
            Me%Array%ILB                    = 1
            Me%Array%IUB                    = sum(Me%ExternalVar%RiverPoints1D)

            !Path to data file
            call ReadInterfaceFilesName

            !Start sinks and sources model
            call StartSinksSourcesModel(DT)

            !Verify compute options
            call Check_Options(PropertiesList)

            !Allocate variables global to whole module
            call AllocateVariables

            STAT_ = SUCCESS_
            
            !Returns ID
            InterfaceID = Me%InstanceID

            STAT_ = SUCCESS_

        else cd0
            
            stop 'ModuleInterface - ConstructInterface2D - ERR01' 

        end if cd0

        if (present(STAT)) STAT = STAT_

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

    end subroutine ConstructInterface1D

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

    subroutine AllocateInstance

        !Local-----------------------------------------------------------------
        type (T_Interface), pointer           :: NewInterface
        type (T_Interface), pointer           :: PreviousInterface

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

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

        Me%InstanceID = RegisterNewInstance (mINTERFACE_)

    end subroutine AllocateInstance


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

    
    logical function  FindProperty(PropertiesList, PropertyID)

        !Arguments-------------------------------------------------------------
        integer, pointer, dimension(:)              :: PropertiesList
        integer, intent (IN)                        :: PropertyID

        !Local-----------------------------------------------------------------
        integer                                     :: i
        
        !Begin-----------------------------------------------------------------
        
        FindProperty = .false.

        do i = 1, size(PropertiesList)  
            if (PropertyID == PropertiesList(i))then
                FindProperty =.true.
                exit
            end if
        enddo

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

    end function FindProperty


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

    
    subroutine AllocateVariables

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

        !Local-----------------------------------------------------------------
        integer                       :: ArrayLB, ArrayUB
        integer                       :: PropLB, PropUB

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

        PropLB  = Me%Prop%ILB
        PropUB  = Me%Prop%IUB

        ArrayLB = Me%Array%ILB
        ArrayUB = Me%Array%IUB

        
        !Allocate global variables--------------------------------------------- 
        allocate(Me%Mass(PropLB:PropUB, ArrayLB:ArrayUB), STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR01'

        allocate(Me%ConcentrationIncrement(PropLB:PropUB, ArrayLB:ArrayUB), &
                                                     STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR02'

        allocate(Me%Temperature(ArrayLB:ArrayUB), STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR03'

        allocate(Me%OpenPoints(ArrayLB:ArrayUB), STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR04'
        
        allocate(Me%AddedProperties(PropLB:PropUB), STAT = STAT_CALL)
        if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR05'
        
        
        select case (Me%SinksSourcesModel)

            case (WaterQualityModel)

                allocate(Me%Salinity (ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR07'
        
                allocate(Me%LightExtCoefField(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR08'
        
                allocate(Me%ShortWaveRadiation(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR09'

                allocate(Me%Thickness(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR10'

                allocate(Me%FishFood(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR11'
                       
                Me%Salinity           = FillValueReal
                Me%FishFood           = FillValueReal
                Me%LightExtCoefField  = FillValueReal
                Me%ShortWaveRadiation = FillValueReal
                Me%Thickness          = FillValueReal


            case (SedimentQualityModel)
                
                allocate (Me%WaterPercentage(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR12'
                
                Me%WaterPercentage = FillValueReal

                allocate (Me%DissolvedToParticulate(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR13'
                
                Me%DissolvedToParticulate = FillValueReal
         

            case (CEQUALW2Model)

                allocate(Me%Salinity (ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR14'
        
                allocate(Me%LightExtCoefField(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR15'
        
                allocate(Me%ShortWaveRadiation(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR16'

                allocate(Me%Thickness(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR17'

                allocate(Me%Alkalinity(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR18'
                       
                Me%Salinity           = FillValueReal
                Me%Alkalinity         = FillValueReal
                Me%LightExtCoefField  = FillValueReal
                Me%ShortWaveRadiation = FillValueReal
                Me%Thickness          = FillValueReal



            case (LifeModel)

                allocate(Me%Salinity (ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR19'
        
                allocate(Me%LightExtCoefField(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR20'
        
                allocate(Me%ShortWaveRadiation(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR21'

                allocate(Me%Thickness(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR22'

                allocate(Me%Alkalinity(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR23'
                       
                Me%Salinity           = FillValueReal
                Me%Alkalinity         = FillValueReal
                Me%LightExtCoefField  = FillValueReal
                Me%ShortWaveRadiation = FillValueReal
                Me%Thickness          = FillValueReal



            case (BenthicCEQUALW2Model)
            
            
            case (BenthosModel)
            
                allocate(Me%Oxygen(ArrayLB:ArrayUB), STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'AllocateVariables - ModuleInterface - ERR24'

            
            case default
                write(*,*) 
                write(*,*) 'Defined sinks and sources model was not recognised.'
                stop 'AllocateVariables - ModuleInterface - ERR25'
        end select

        Me%Mass                       = FillValueReal
        Me%ConcentrationIncrement     = FillValueReal
        Me%Temperature                = FillValueReal
        Me%OpenPoints                 = FillValueInt
        Me%AddedProperties            = .false.

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

    end subroutine AllocateVariables   
    

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


    subroutine ReadInterfaceFilesName 

        !External--------------------------------------------------------------
        integer                                     :: STAT_CALL
        
        !Local-----------------------------------------------------------------
        character(Len = StringLength)               :: Message

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

        select case (Me%SinksSourcesModel)

            case (WaterQualityModel)

                Message  =trim('Water Quality Data File')

                call ReadFileName('WQDATA',Me%FileName, Message = Message, STAT = STAT_CALL)

cd1 :           if(STAT_CALL .EQ. KEYWORD_NOT_FOUND_ERR_) then

                    call ReadFileName('DISPQUAL',                    &
                                      Me%FileName,                   &
                                      Message = Message,             &
                                      STAT    = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_)stop 'ReadInterfaceFilesName - ModuleInterface - ERR01' 

                else if (STAT_CALL .EQ. SUCCESS_) then cd1
                    continue
                else
                    stop 'ReadInterfaceFilesName - ModuleInterface - ERR02' 
                end if cd1


            case (SedimentQualityModel)

                Message  = trim('SedimentQuality Data File')
                call ReadFileName('SQ_DATA',Me%FileName, Message = Message, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'ReadInterfaceFilesName - Module Interface - ERR03' 
            

            case(CEQUALW2Model)
                
                Message  = trim('CEQUALW2 Data File')
                call ReadFileName('CEQUALW2_DATA',Me%FileName, Message = Message, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'ReadInterfaceFilesName - Module Interface - ERR04' 
            

            case(LifeModel)
                
                Message  = trim('Life Data File')
                call ReadFileName('LIFE_DATA',Me%FileName, Message = Message, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'ReadInterfaceFilesName - Module Interface - ERR05' 


            case(BenthicCEQUALW2Model)

                Message  = trim('Benthic CEQUALW2 Data File')
                call ReadFileName('BENTHIC_CEQUALW2_DATA',Me%FileName, Message = Message, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'ReadInterfaceFilesName - Module Interface - ERR06' 

            case(BenthosModel)

                Message  = trim('Benthos Data File')
                call ReadFileName('BENTHOS_DATA',Me%FileName, Message = Message, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'ReadInterfaceFilesName - Module Interface - ERR07' 

            case default
                write(*,*) 
                write(*,*) 'Defined sinks and sources model was not recognised.'
                stop 'ReadInterfaceFilesName - ModuleInterface - ERR08' 
        end select


        Me%FileName = trim(Me%FileName)

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

    end subroutine ReadInterfaceFilesName 

    !--------------------------------------------------------------------------
    
    subroutine StartSinksSourcesModel (DT)

        !Arguments-------------------------------------------------------------
        real, intent(OUT)                       :: DT
        
        !External--------------------------------------------------------------
        integer                                 :: STAT_CALL
        integer                                 :: PropLB, PropUB

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

        select case (Me%SinksSourcesModel)

            case (WaterQualityModel)

                !Construct WaterQuality Model
                call StartWaterQuality(Me%ObjWaterQuality,                     &
                                       Me%FileName,                            &
                                       STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR01'
 
                !Construct mass fluxes between properties
                call Construct_WQRateFlux(Me%ObjWaterQuality,                  &
                                          Me%Array%ILB,                        &
                                          Me%Array%IUB,                        &
                                          STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR02'
               
                !Number of properties involved
                call GetWaterQualitySize(Me%ObjWaterQuality,                   &
                                         PropLB = PropLB,                      &
                                         PropUB = PropUB,                      &
                                         STAT   = STAT_CALL)                            
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR03'

                !Number of properties involved
                Me%Prop%ILB     = PropLB
                Me%Prop%IUB     = PropUB

                !Get DT from WaterQuality model to exit as argument to WaterProperties
                call GetDTWQM(Me%ObjWaterQuality, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR04'


            case (SedimentQualityModel)

                !Construct SedimentQuality Model
                call StartSedimentQuality(Me%ObjSedimentQuality,               &
                                          Me%FileName,                         &
                                          STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR05'
                
                !Construct mass fluxes between properties
                call Construct_SQRateFlux(Me%ObjSedimentQuality,               &
                                          Me%Array%ILB,                        &
                                          Me%Array%IUB,                        &
                                         STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR06'
                
                !Number of properties involved
                call GetSedimentQualitySize(Me%ObjSedimentQuality,             &
                                            PropLB = PropLB,                   &
                                            PropUB = PropUB,                   &
                                            STAT   = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR07'

                
                !Number of properties involved
                Me%Prop%ILB     = PropLB
                Me%Prop%IUB     = PropUB

                !Get DT from SedimentQuality model to exit as argument to SedimentProperties
                call GetDTSedimentQuality(Me%ObjSedimentQuality, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR08'


            case(LifeModel)

                !Construct Life Model
                call ConstructLife(Me%ObjLife,                                  &
                                   FileName = Me%FileName,                      &
                                   STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR09'
 
                !Number of properties involved
                call GetLifeSize(Me%ObjLife,                                    &
                                     PropLB = PropLB,                           &
                                     PropUB = PropUB,                           &
                                     STAT   = STAT_CALL)                            
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR010'

                !Number of properties involved
                Me%Prop%ILB = PropLB
                Me%Prop%IUB = PropUB

                !Get DT from Life model to exit as argument to WaterProperties
                call GetDTLife(Me%ObjLife, DT = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR011'


            case(CEQUALW2Model)

                !Construct CEQUALW2 Model
                call StartCEQUALW2(Me%ObjCEQUALW2,                              &
                                   Me%Array%ILB,                                &
                                   Me%Array%IUB,                                &
                                   FileName = Me%FileName,                      &
                                   Model = CEQUALW2Model,                       & 
                                   STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR09'
 
                !Number of properties involved
                call GetCEQUALW2Size(Me%ObjCEQUALW2,                            &
                                     PropLB = PropLB,                           &
                                     PropUB = PropUB,                           &
                                     STAT   = STAT_CALL)                            
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR010'

                !Number of properties involved
                Me%Prop%ILB = PropLB
                Me%Prop%IUB = PropUB

                !Get DT from CEQUALW2 model to exit as argument to WaterProperties
                call GetDTCEQUALW2(Me%ObjCEQUALW2, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR011'


            case(BenthicCEQUALW2Model)

                !Construct CEQUALW2 Model
                call StartCEQUALW2(Me%ObjCEQUALW2,                              &
                                   Me%Array%ILB,                                &
                                   Me%Array%IUB,                                &
                                   FileName = Me%FileName,                      &
                                   Model = BenthicCEQUALW2Model,                & 
                                   STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR012'
 
                !Number of properties involved
                call GetCEQUALW2Size(Me%ObjCEQUALW2,                            &
                                     PropLB = PropLB,                           &
                                     PropUB = PropUB,                           &
                                     STAT   = STAT_CALL)                            
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR012'

                !Number of properties involved
                Me%Prop%ILB = PropLB
                Me%Prop%IUB = PropUB

                !Get DT from CEQUALW2 model to exit as argument to WaterProperties
                call GetDTCEQUALW2(Me%ObjCEQUALW2, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR014'
            
            case(BenthosModel)


                !Construct Benthos Model
                call StartBenthos (Me%ObjBenthos,                               &
                                   FileName = Me%FileName,                      &
                                   STAT = STAT_CALL) 
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR015'
 
                !Number of properties involved
                call GetBenthosSize (Me%ObjBenthos,                             &
                                     PropLB = PropLB,                           &
                                     PropUB = PropUB,                           &
                                     STAT   = STAT_CALL)                            
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR016'

                !Number of properties involved
                Me%Prop%ILB = PropLB
                Me%Prop%IUB = PropUB

                !Get DT from CEQUALW2 model to exit as argument to WaterProperties
                call GetDTBenthos(Me%ObjBenthos, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'StartSinksSourcesModel - ModuleInterface - ERR017'



            case default
                write(*,*) 
                write(*,*) 'Defined sinks and sources model was not recognised.'
                stop 'StartSinksSourcesModel - ModuleInterface - ERR00' 
        end select

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


    end subroutine StartSinksSourcesModel    
    
    
    !--------------------------------------------------------------------------
   
   
    subroutine Check_Options(PropertiesList)

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

        !External--------------------------------------------------------------
        type(T_Time)                                         :: EndTime
        logical                                              :: Zoo, Phyto
        logical                                              :: Diatoms 
        logical                                              :: Nitrogen, Phosphorus
        logical                                              :: Silica  
        logical                                              :: Oxygen, BOD
        logical                                              :: Carbon
        logical                                              :: Bacteria, Ciliate 
        logical                                              :: Larvae
        real                                                 :: DT
        integer                                              :: STAT_CALL

        !Local-----------------------------------------------------------------
        real                                                 :: ErrorAux, auxFactor 
        real                                                 :: RunPeriod, Dtlag
        integer, dimension(:), pointer                       :: CEQUALW2List
        integer, dimension(:), pointer                       :: LifeList, BenthosList
        integer                                              :: i,PropLB, PropUB
        !----------------------------------------------------------------------
        
        call GetComputeCurrentTime(Me%ObjTime,                  &
                                   Me%ExternalVar%Now,          &
                                   STAT = STAT_CALL)                    
        if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR00' 

        !Get end time
        call GetComputeTimeLimits(Me%ObjTime, EndTime = EndTime, STAT = STAT_CALL)
        if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR01' 
            
        select case (Me%SinksSourcesModel)

            case (WaterQualityModel)
                !Get water quality model time step
                call GetDTWQM(Me%ObjWaterQuality, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'Check_Options - ModuleInterface - ERR02' 
        
                !Run period in seconds
                RunPeriod = EndTime - Me%ExternalVar%Now

                !The run period must be a multiple of the WQ DT
                auxFactor = RunPeriod / DT

                ErrorAux  = auxFactor - int(auxFactor)
                
                if (ErrorAux /= 0) then
                    Dtlag = int(ErrorAux * DT)
                    write(*,*) 
                    write(*,*) 'DTSECONDS is not multiple of the run period.'
                    write(*,*) 'Water Quality wont be computed in the last', Dtlag, ' seconds.'
                    write(*,*) 'Check_Options - ModuleInterface - WRN01.'
                endif 

                !Sinks and sources compute options 
                call GetWQOptions(Me%ObjWaterQuality,   Phyto            = Phyto,           &
                                                        Nitrogen         = Nitrogen,        &
                                                        Diatoms          = Diatoms,         &
                                                        Silica           = Silica,          & 
                                                        Phosphorus       = Phosphorus,      &
                                                        Zoo              = Zoo,             &
                                                        BOD              = BOD,             &
                                                        Oxygen           = Oxygen,          &
                                                        Bacteria         = Bacteria,        &
                                                        Ciliate          = Ciliate,         &
                                                        Larvae           = Larvae,          &
                                                        STAT             = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR03' 


cd1 :           if (Phyto) then 
                    if (.not.FindProperty(PropertiesList, Phytoplankton_))          &
                        stop 'WQM needs property phytoplankton - Check_Options'
                end if cd1

cd2 :           if (Zoo) then
                    if (.not.FindProperty(PropertiesList, Zooplankton_))            &
                        stop 'WQM needs property zooplankton - Check_Options'
                end if cd2

cd3 :           if (Larvae) then
                    if (.not.FindProperty(PropertiesList, Larvae_))                 &
                        stop 'WQM needs property Larvae - Check_Options'
                end if cd3

cd5 :           if (Phosphorus) then

                    if (.not.FindProperty(PropertiesList, POP_))                    &
                        stop 'WQM needs property particulate organic phosphorus - Check_Options'

                    if (.not.FindProperty(PropertiesList, DOPRefractory_))          &
                        stop 'WQM needs property dissolved refractory organic phosphorus - Check_Options'

                    if (.not. FindProperty(PropertiesList, DOPNon_Refractory_))     &
                        stop 'WQM needs property dissolved non-refractory organic phosphorus - Check_Options'

                    if (.not.FindProperty(PropertiesList, Inorganic_Phosphorus_))   &
                        stop 'WQM needs property inorganic phosphorus - Check_Options'
                end if cd5

cd6 :           if (Nitrogen) then
                    if (.not.FindProperty(PropertiesList, PON_))                    &
                        stop 'WQM needs property particulate organic nitrogen - Check_Options'

                    if (.not.FindProperty(PropertiesList, DONRefractory_))          &
                        stop 'WQM needs property dissolved refractory organic nitrogen - Check_Options'

                    if (.not. FindProperty(PropertiesList, DONNon_Refractory_))     &
                        stop 'WQM needs property dissolved non-refractory organic nitrogen - Check_Options'

                    if (.not.FindProperty(PropertiesList, Ammonia_))                &
                        stop 'WQM needs property ammonia - Check_Options'

                    if (.not.FindProperty(PropertiesList, Nitrate_))                &
                        stop 'WQM needs property nitrate - Check_Options'

                    if (.not.FindProperty(PropertiesList, Nitrite_))                &
                        stop 'WQM needs property nitrite - Check_Options'
                end if cd6

cd7 :           if (BOD) then
                    if (.not.FindProperty(PropertiesList, BOD_))                    &
                        stop 'WQM needs property biochemical oxygen demand - Check_Options'
                end if cd7


cd8 :           if (Bacteria) then
                    if (.not.FindProperty(PropertiesList, Bacteria_))               &
                        stop 'WQM needs property bacteria - Check_Options'
                end if cd8


cd9 :           if (Ciliate) then
                    if (.not.FindProperty(PropertiesList, Ciliate_))                &
                        stop 'WQM needs property ciliate - Check_Options'
                end if cd9


cd12 :          if (Diatoms) then
                    if (.not.FindProperty(PropertiesList, Diatoms_))                &
                        stop 'WQM needs property diatoms - Check_Options'
                end if cd12

cd13 :          if (Silica) then
                    if (.not.FindProperty(PropertiesList, DSilica_))                &
                        stop 'WQM needs property diatoms - Check_Options'

                    if (.not.FindProperty(PropertiesList, BioSilica_))               &
                        stop 'WQM needs property diatoms - Check_Options'
                end if cd13

                !If one wants to run just the age property oxygen is not needed - Frank 09-2003
                if (Phosphorus .or. BOD .or. Ciliate .or. Nitrogen .or. Bacteria) then
                    if (.not.FindProperty(PropertiesList, Oxygen_))                     &
                        stop 'WQM needs property oxygen - Check_Options'
                endif

            
            case (SedimentQualityModel)

                !Get SedimentQuality model time step
                call GetDTSedimentQuality(Me%ObjSedimentQuality, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR04' 
        
                !Run period in seconds
                RunPeriod = EndTime - Me%ExternalVar%Now

                !The run period must be a multiple of the SedimentQuality DT
                auxFactor = RunPeriod / DT

                ErrorAux  = auxFactor - int(auxFactor)
                
                if (ErrorAux /= 0) then
                    Dtlag = int(ErrorAux * DT)
                    write(*,*) 
                    write(*,*) 'DTSECONDS is not multiple of the run period.'
                    write(*,*) 'SedimentQuality wont be computed in the last', Dtlag, ' seconds.'
                    write(*,*) 'Check_Options - ModuleInterface - WRN02.'
                endif
                
              !Sinks and sources compute options 
                call GetSQOptions(Me%ObjSedimentQuality,  Nitrogen    = Nitrogen  ,   &
                                                                    Carbon      = Carbon    ,   &
                                                                    Oxygen      = Oxygen    ,   &
                                                                    STAT        = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR05' 

cd10 :           if (Nitrogen) then
                   
                    if (.not.FindProperty(PropertiesList, Ammonia_)              )   &
                        stop 'SQM needs property ammonia - Check_Options'

                    if (.not.FindProperty(PropertiesList, Nitrate_)              )   &
                        stop 'SQM needs property nitrate - Check_Options'

                    !if (.not.FindProperty(PropertiesList, AdsorbedAmmonia_)      )   &
                     !   STOP 'WQM needs property Adsorbed Ammonia - Check_Options'
                    
                    if (.not.FindProperty(PropertiesList, RefreactaryOrganicN_)  )   &
                        stop 'SQM needs property Refreactary N - Check_Options'

                    if (.not.FindProperty(PropertiesList, PON_           )       )   &
                        stop 'SQM needs property Labile N - Check_Options'

                    if (.not.FindProperty(PropertiesList, Ngas_)                 )   &
                        stop 'SQM needs property N gas - Check_Options'

                    if (.not.FindProperty(PropertiesList, HeterotrophicN_)       )   &
                        stop 'SQM needs property Heterotrophic N - Check_Options'

                    if (.not.FindProperty(PropertiesList, AnaerobicN_)           )   &
                        stop 'SQM needs property Anaerobic N - Check_Options'

                    if (.not.FindProperty(PropertiesList, AutotrophicN_)         )   &
                        STOP 'SQM needs property Autotrophic N - Check_Options'

                end if cd10

cd11 :           if (Carbon) then
                   
                    if (.not.FindProperty(PropertiesList, LabileOrganicC_)       )   &
                        stop 'SQM needs property Labile Organic C - Check_Options'

                    if (.not.FindProperty(PropertiesList, RefreactaryOrganicC_)  )   &
                        stop 'SQM needs property Refreactary Organic C - Check_Options'

                    if (.not.FindProperty(PropertiesList, HeterotrophicC_)       )   &
                        stop 'SQM needs property Heterotrophic C - Check_Options'

                    if (.not.FindProperty(PropertiesList, AnaerobicC_)           )   &
                        stop 'SQM needs property Anaerobic C - Check_Options'
                    
                    if (.not.FindProperty(PropertiesList, AutotrophicC_)         )   &
                        stop 'SQM needs property Autotrophic C - Check_Options'


                end if cd11

                if (.not.FindProperty(PropertiesList, Oxygen_))                     &
                    stop 'SQM needs property oxygen - Check_Options'         
            
             
            case(CEQUALW2Model, BenthicCEQUALW2Model)
                
                !Get CEQUALW2 model time step
                call GetDTCEQUALW2(Me%ObjCEQUALW2, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR05' 
        
                !Run period in seconds
                RunPeriod = EndTime - Me%ExternalVar%Now

                !The run period must be a multiple of the SedimentQuality DT
                auxFactor = RunPeriod / DT

                ErrorAux  = auxFactor - int(auxFactor)
                
                if (ErrorAux /= 0) then
                    Dtlag = int(ErrorAux * DT)
                    write(*,*) 
                    write(*,*) 'DTSECONDS is not multiple of the run period.'
                    write(*,*) 'CEQUALW2 wont be computed in the last', Dtlag, ' seconds.'
                    write(*,*) 'Check_Options - ModuleInterface - WRN03.'
                endif
                
                call GetCEQUALW2PropertyList (Me%ObjCEQUALW2, CEQUALW2List, STAT=STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR06'
                
                !Number of properties involved
                PropLB = Me%Prop%ILB
                PropUB = Me%Prop%IUB

                do i = PropLB, PropUB
                    if (.not.FindProperty(PropertiesList, CEQUALW2List(i))) then
                        write(*,*) 'Property ',GetPropertyName(CEQUALW2List(i)),' not found'
                        stop 'Properties lists inconsistent  - Check_Options- ModuleInterface- ERR07'    
                    end if
                end do

                call UngetCEQUALW2 (Me%ObjCEQUALW2, CEQUALW2List, STAT=STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR08'


            case(LifeModel)

                !Get Life model time step
                call GetDTLife(Me%ObjLife, DT = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR09' 
        
                !Run period in seconds
                RunPeriod = EndTime - Me%ExternalVar%Now

                !The run period must be a multiple of the SedimentQuality DT
                auxFactor = RunPeriod / DT

                ErrorAux  = auxFactor - int(auxFactor)
                
                if (ErrorAux /= 0) then
                    Dtlag = int(ErrorAux * DT)
                    write(*,*) 
                    write(*,*) 'DTSECONDS is not multiple of the run period.'
                    write(*,*) 'Life wont be computed in the last', Dtlag, ' seconds.'
                    write(*,*) 'Check_Options - ModuleInterface - WRN03.'
                endif
                
                call GetLifePropertyList (Me%ObjLife, LifeList, STAT=STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR10'
                
                !Number of properties involved
                PropLB = Me%Prop%ILB
                PropUB = Me%Prop%IUB

                do i = PropLB, PropUB
                    if (.not.FindProperty(PropertiesList, LifeList(i))) then
                        write(*,*) 'Property ',GetPropertyName(LifeList(i)),' not found'
                        stop 'Properties lists inconsistent  - Check_Options- ModuleInterface- ERR11'    
                    end if
                end do

                call UngetLife (Me%ObjLife, LifeList, STAT=STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR12'


            case(BenthosModel)

                !Get Benthos model time step
                call GetDTBenthos(Me%ObjBenthos, DTSecond = DT, STAT = STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR13' 
        
                !Run period in seconds
                RunPeriod = EndTime - Me%ExternalVar%Now

                !The run period must be a multiple of the SedimentQuality DT
                auxFactor = RunPeriod / DT

                ErrorAux  = auxFactor - int(auxFactor)
                
                if (ErrorAux /= 0) then
                    Dtlag = int(ErrorAux * DT)
                    write(*,*) 
                    write(*,*) 'DTSECONDS is not multiple of the run period.'
                    write(*,*) 'Life wont be computed in the last', Dtlag, ' seconds.'
                    write(*,*) 'Check_Options - ModuleInterface - WRN04.'
                endif
                
                call GetBenthosPropertyList (Me%ObjBenthos, BenthosList, STAT=STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR14'
                
                !Number of properties involved
                PropLB = Me%Prop%ILB
                PropUB = Me%Prop%IUB

                do i = PropLB, PropUB
                    if (.not.FindProperty(PropertiesList, BenthosList(i))) then
                        write(*,*) 'Property ',GetPropertyName(BenthosList(i)),' not found'
                        stop 'Properties lists inconsistent  - Check_Options- ModuleInterface- ERR15'    
                    end if
                end do

                call UngetBenthos (Me%ObjBenthos, BenthosList, STAT=STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR12'

            case default
                write(*,*) 
                write(*,*) 'Defined sinks and sources model was not recognised.'
                if (STAT_CALL /= SUCCESS_) stop 'Check_Options - ModuleInterface - ERR13' 
        end select


        call null_time   (Me%ExternalVar%Now)

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

    end subroutine Check_Options

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

    !SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR SELECTOR 

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

    !--------------------------------------------------------------------------
   
    subroutine GetRateFlux3D(InterfaceID,                            & 
                             FirstProp,                              &
                             SecondProp,                             &
                             RateIndex,                              &
                             RateFlux3D,                             &
                             WaterPoints3D,                          &
                             STAT)

        !Arguments--------------------------------------------------------------
        integer                                         :: InterfaceID
        real,               dimension(:,:,:),   pointer :: RateFlux3D
        integer,            dimension(:,:,:),   pointer :: WaterPoints3D
        integer, optional,  intent(IN)                  :: FirstProp
        integer, optional,  intent(IN)                  :: SecondProp
        integer, optional,  intent(IN)                  :: RateIndex
        integer, optional,  intent(OUT)                 :: STAT

        !External--------------------------------------------------------------
        integer                                         :: ready_, STAT_CALL          

        !Local-----------------------------------------------------------------
        integer                                         :: STAT_
        integer                                         :: nFirstProp, nSecondProp          
        integer                                         :: ILB, IUB, JLB, JUB, KLB, KUB
        integer                                         :: Index
        integer                                         :: i, j, k
        real,    dimension(:), pointer                  :: RateFlux

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

        STAT_ = UNKNOWN_

        call Ready(InterfaceID, ready_)   
         
        ILB     = Me%Size3D%ILB     
        IUB     = Me%Size3D%IUB     
        JLB     = Me%Size3D%JLB    
        JUB     = Me%Size3D%JUB    
        KLB     = Me%Size3D%KLB   
        KUB     = Me%Size3D%KUB   

        nullify (RateFlux)

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

            if (present(FirstProp)) then

                !Number indexed to each property 
                nFirstProp  = PropertyIndexNumber(FirstProp )
                nSecondProp = PropertyIndexNumber(SecondProp)

                !WaterPoints3D
                Me%ExternalVar%WaterPoints3D => WaterPoints3D
                
                select case (Me%SinksSourcesModel)

                    case (WaterQualityModel)

                        call GetWQPropRateFlux( Me%ObjWaterQuality,                          &
                                                nFirstProp, nSecondProp,                     &
                                                RateFlux, STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'GetRateFlux3D - ModuleInterface - ERR01'

                    case (SedimentQualityModel)
                    
                        call GetPropRateFlux(Me%ObjSedimentQuality,                      &
                                             nFirstProp, nSecondProp,                    &
                                             RateFlux, STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'GetRateFlux3D - ModuleInterface - ERR02' 

                    case(LifeModel)

                end select

            elseif (present(RateIndex)) then


                call GetCEQUALW2RateFlux( Me%ObjCeQualW2,                                  &
                                          RateIndex, RateFlux,                             &
                                          STAT_CALL)
                if (STAT_CALL /= SUCCESS_) stop 'GetRateFlux3D - ModuleInterface - ERR03'


            endif

            
            !Number indexed to 3D cell in the vector 
            Index = 0 

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

                if (Me%ExternalVar%WaterPoints3D(i, j, k) == 1) then

                    Index = Index + 1
                    RateFlux3D(i, j, k) = RateFlux (Index)

                end if
            end do
            end do
            end do


            select case (Me%SinksSourcesModel)

                case (WaterQualityModel)

                    call UnGetWQPropRateFlux(Me%ObjWaterQuality, RateFlux, STAT_CALL)
                    if (STAT_CALL /= SUCCESS_) stop 'GetRateFlux3D - ModuleInterface - ERR03' 
                
                case (SedimentQualityModel)

                    call UngetPropRateFlux(Me%ObjSedimentQuality, RateFlux, STAT_CALL)
                    if (STAT_CALL /= SUCCESS_) stop 'GetRateFlux3D - ModuleInterface - ERR04' 

                case(CEQUALW2Model, BenthicCEQUALW2Model)
               
                   call UnGetCEQUALW2RateFlux(Me%ObjCeQualW2, RateFlux, STAT_CALL)
                   if (STAT_CALL /= SUCCESS_) stop 'GetRateFlux3D - ModuleInterface - ERR06'
             
            end select

            STAT_ = SUCCESS_
        else
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_


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

    end subroutine GetRateFlux3D

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


    subroutine GetRateFlux2D(InterfaceID,                            & 
                             FirstProp,                              &
                             SecondProp,                             &
                             RateFlux2D,                             &
                             WaterPoints2D,                          &
                             STAT)

        !Arguments--------------------------------------------------------------
        integer                                         :: InterfaceID
        real,               dimension(:,:  ),   pointer :: RateFlux2D
        integer,            dimension(:,:  ),   pointer :: WaterPoints2D
        integer,            intent(IN)                  :: FirstProp
        integer,            intent(IN)                  :: SecondProp
        integer, optional,  intent(OUT)                 :: STAT

        !External--------------------------------------------------------------
        integer                                         :: ready_         

        !Local-----------------------------------------------------------------
        integer                                         :: STAT_
        integer                                         :: nFirstProp, nSecondProp          
        integer                                         :: ILB, IUB, JLB, JUB
        integer                                         :: Index
        integer                                         :: i, j
        real,    dimension(:), pointer                  :: RateFlux

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

        STAT_ = UNKNOWN_

        call Ready(InterfaceID, ready_)   
         
        ILB     = Me%Size2D%ILB     
        IUB     = Me%Size2D%IUB     
        JLB     = Me%Size2D%JLB    
        JUB     = Me%Size2D%JUB    

        nullify (RateFlux  )
        nullify (RateFlux2D)

        !Number indexed to each property 
        nFirstProp  = PropertyIndexNumber(FirstProp )
        nSecondProp = PropertyIndexNumber(SecondProp)


cd1 :   if ((ready_ .EQ. IDLE_ERR_) .OR. (ready_ .EQ. READ_LOCK_ERR_)) then
            
            call Read_Lock(mINTERFACE_, Me%InstanceID)
            
            !WaterPoints3D
            Me%ExternalVar%WaterPoints2D => WaterPoints2D
                
            !Number indexed to 3D cell in the vector 
           Index = 0 

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

                if (Me%ExternalVar%WaterPoints2D(i, j) == 1) then

                    Index = Index + 1
                    RateFlux2D(i, j) = RateFlux (Index)

                end if
            end do
            end do


            STAT_ = SUCCESS_
        else

            STAT_ = ready_

        end if cd1

        if (present(STAT)) STAT = STAT_

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

    end subroutine GetRateFlux2D

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

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

    !MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER  

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

    subroutine Modify_Interface3D(InterfaceID, PropertyID, Concentration,         &
                                  WaterPoints3D, OpenPoints3D, DWZ,               &
                                  ShortWaveRadiation, LightExtCoefField,          &
                                  WaterPercentage, DissolvedToParticulate3D,      &
                                  DTProp, STAT)

        !Arguments-------------------------------------------------------------
        integer                                         :: InterfaceID
        integer, intent(IN)                             :: PropertyID
        real,              dimension(:,:,:), pointer    :: Concentration
        real,    optional, dimension(:,:,:), pointer    :: ShortWaveRadiation
        real,    optional, dimension(:,:,:), pointer    :: LightExtCoefField
        real   , optional, dimension(:,:,:), pointer    :: WaterPercentage
        real,    optional, dimension(:,:,:), pointer    :: DissolvedToParticulate3D
        integer,           dimension(:,:,:), pointer    :: WaterPoints3D
        integer, optional, dimension(:,:,:), pointer    :: OpenPoints3D
        real,    optional, dimension(:,:,:), pointer    :: DWZ
        real,    optional,  intent(IN)                  :: DTProp
        integer, optional,  intent(OUT)                 :: STAT

        !External--------------------------------------------------------------
        real                                            :: DT
        integer                                         :: STAT_CALL
        integer                                         :: ready_
        logical                                         :: ReadyToCompute

        !Local-----------------------------------------------------------------
        integer                                         :: STAT_
        integer                                         :: nProperty
        integer                                         :: Index 
        integer                                         :: i, j, k
        integer                                         :: prop, JulDay
        integer                                         :: ILB, IUB, JLB, JUB, KLB, KUB     
        integer                                         :: PropLB, PropUB, ArrayLB, ArrayUB 
        real                                            :: DTProp_
        logical                                         :: Increment

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

        STAT_ = UNKNOWN_

        call Ready(InterfaceID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

            ILB     = Me%Size3D%ILB
            IUB     = Me%Size3D%IUB
            JLB     = Me%Size3D%JLB
            JUB     = Me%Size3D%JUB
            KLB     = Me%Size3D%KLB
            KUB     = Me%Size3D%KUB

            PropLB  = Me%Prop%ILB
            PropUB  = Me%Prop%IUB

            ArrayLB = Me%Array%ILB
            ArrayUB = Me%Array%IUB

            if(present(OpenPoints3D ))Me%ExternalVar%OpenPoints3D     => OpenPoints3D
            if(present(DWZ          ))Me%ExternalVar%DWZ              => DWZ
            
            Me%ExternalVar%WaterPoints3D => WaterPoints3D
            
            ReadyToCompute = .false.

            Increment = OFF

cd7 :       if (present(DTProp))then
                DTProp_   = DTProp
                Increment = ON  !If Increment is ON  then the concentrations is will be changed
            else
                DTProp_   = FillValueReal
                Increment = OFF !If Increment is OFF then the sinks and sources module will be computed
            end if cd7


cd5 :       if (.not. Increment) then 

                call FillMassTempSalinity(PropertyID, Concentration, ReadyToCompute)

cd4 :           if (ReadyToCompute) then

                    call UnfoldMatrix(Me%ExternalVar%OpenPoints3D, Me%OpenPoints, Vertical1D = Me%ExternalVar%Vertical1D)

                    !Stores the concentration before changing them
                    Me%ConcentrationIncrement = Me%Mass

                    select case (Me%SinksSourcesModel)

                        case(WaterQualityModel)

                            call UnfoldMatrix(ShortWaveRadiation,Me%ShortWaveRadiation)
                            call UnfoldMatrix(LightExtCoefField, Me%LightExtCoefField)
                            call UnfoldMatrix(Me%ExternalVar%DWZ,Me%Thickness)


                            call WaterQuality(Me%ObjWaterQuality,                   &
                                              Me%Salinity,                          &
                                              Me%Temperature,                       &
                                              Me%ShortWaveRadiation,                &
                                              Me%LightExtCoefField,                 &
                                              Me%Thickness,                         &
                                              Me%Mass,                              &
                                              Me%Array%ILB,                         &
                                              Me%Array%IUB,                         &
                                              Me%OpenPoints,                        &
                                              FishFood = Me%FishFood,               &
                                              STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR05'
                        
                        case(SedimentQualityModel)

                            call UnfoldMatrix(WaterPercentage, Me%WaterPercentage)
                            call UnfoldMatrix(DissolvedToParticulate3D, Me%DissolvedToParticulate)

                            call SedimentQuality(Me%ObjSedimentQuality,             &
                                                 Me%Temperature,                    &
                                                 Me%Mass,                           &
                                                 Me%WaterPercentage,                &
                                                 Me%DissolvedToParticulate,         &
                                                 Me%Array%ILB,                      &
                                                 Me%Array%IUB,                      &
                                                 STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR06'

                        case(LifeModel)

                            call UnfoldMatrix(ShortWaveRadiation,Me%ShortWaveRadiation)
                            call UnfoldMatrix(LightExtCoefField, Me%LightExtCoefField)
                            call UnfoldMatrix(Me%ExternalVar%DWZ,Me%Thickness)
                            
                            call GetComputeCurrentTime(Me%ObjTime,                  &
                                                       Me%ExternalVar%Now,          &
                                                       STAT = STAT_CALL)                    
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERRXX' 


                            call JulianDay(Me%ExternalVar%Now, JulDay)

                            call ModifyLife(Me%ObjLife,                     &
                                      Me%Salinity,                          &
                                      Me%Temperature,                       &
                                      Me%ShortWaveRadiation,                &
                                      Me%LightExtCoefField,                 &
                                      Me%Thickness,                         &
                                      Me%Mass,                              &
                                      Me%OpenPoints,                        &
                                      Me%Array,                             &
                                      JulDay,                               &
                                    STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR07'



                        case(CEQUALW2Model)

                            call UnfoldMatrix(ShortWaveRadiation,Me%ShortWaveRadiation)
                            call UnfoldMatrix(LightExtCoefField, Me%LightExtCoefField)
                            call UnfoldMatrix(Me%ExternalVar%DWZ,Me%Thickness)
                            

                            call CEQUALW2(Me%ObjCEQUALW2,                       &
                                          Me%Salinity,                          &
                                          Me%Temperature,                       &
                                          Me%Alkalinity,                        &
                                          Me%ShortWaveRadiation,                &
                                          Me%LightExtCoefField,                 &
                                          Me%Thickness,                         &
                                          Me%Mass,                              &
                                          Me%OpenPoints,                        &
                                          STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR08'
                        

                    end select

do7 :               do prop  = PropLB,  PropUB
do6 :               do index = ArrayLB, ArrayUB
                        Me%ConcentrationIncrement(prop, index) = Me%Mass(prop, index) - &
                                                                 Me%ConcentrationIncrement(prop, index)
                    end do do6
                    end do do7
                
                end if cd4

            elseif (Increment) then cd5

                Index = 0

                nProperty = PropertyIndexNumber(PropertyID)

                select case (Me%SinksSourcesModel)

                    case(WaterQualityModel)
                        
                        call GetDTWQM(Me%ObjWaterQuality, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR09'

do1 :                   do k = KLB, KUB
do2 :                   do j = JLB, JUB
do3 :                   do i = ILB, IUB
cd2 :                       if (Me%ExternalVar%WaterPoints3D(i, j, k) == 1) then
                                Index = Index + 1
                                !Concentrations are only actualized in OpenPoints because of instability
                                !in waterpoints that are not openpoints
                                if (Me%ExternalVar%OpenPoints3D(i, j, k) == 1) then
                                    Concentration(i, j, k) = Concentration( i, j, k)      + &
                                    Me%ConcentrationIncrement(nProperty, Index) * DTProp / DT 
                                end if
                            end if cd2
                        end do do3
                        end do do2
                        end do do1

                    case(SedimentQualityModel)
                        
                        call GetDTSedimentQuality(Me%ObjSedimentQuality, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR10'                

do8 :                   do k = KLB, KUB
do9 :                   do j = JLB, JUB
do10:                   do i = ILB, IUB
cd8 :                       if (Me%ExternalVar%WaterPoints3D(i, j, k) == 1) then
                                Index = Index + 1
                                Concentration(i, j, k) = Concentration( i, j, k)                     + &
                                                         Me%ConcentrationIncrement(nProperty, Index) * &
                                                         DTProp / DT 
                            end if cd8
                        end do do10
                        end do do9
                        end do do8



                    case(CEQUALW2Model)

                        call GetDTCEQUALW2(Me%ObjCEQUALW2, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR11'

do11 :                   do k = KLB, KUB
do12 :                   do j = JLB, JUB
do13 :                   do i = ILB, IUB
cd9 :                       if (Me%ExternalVar%WaterPoints3D(i, j, k) == 1) then
                                Index = Index + 1
                                !Concentrations are only actualized in OpenPoints because of instability
                                !in waterpoints that are not openpoints
                                if (Me%ExternalVar%OpenPoints3D(i, j, k) == 1) then
                                    Concentration(i, j, k) = Concentration( i, j, k)      + &
                                    Me%ConcentrationIncrement(nProperty, Index) * DTProp / DT 
                                end if
                            end if cd9
                        end do do13
                        end do do12
                        end do do11


                   
                    case(LifeModel)

                         call GetDTLife(Me%ObjLife, DT = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface3D - ModuleInterface - ERR12'

do15 :                   do k = KLB, KUB
do16 :                   do j = JLB, JUB
do17 :                   do i = ILB, IUB
cd14 :                       if (Me%ExternalVar%WaterPoints3D(i, j, k) == 1) then
                                Index = Index + 1
                                !Concentrations are only actualized in OpenPoints because of instability
                                !in waterpoints that are not openpoints
                                if (Me%ExternalVar%OpenPoints3D(i, j, k) == 1) then
                                    Concentration(i, j, k) = Concentration( i, j, k)      + &
                                    Me%ConcentrationIncrement(nProperty, Index) * DTProp / DT 
                                end if
                            end if cd14
                        end do do17
                        end do do16
                        end do do15



                end select

            end if cd5

            STAT_ = SUCCESS_
        else               
            STAT_ = ready_
        end if cd1


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

    end subroutine Modify_Interface3D

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

    subroutine Modify_Interface2D(InterfaceID, PropertyID, Concentration,               &
                                  WaterPoints2D, OpenPoints2D, Oxygen2D,                &
                                  DTProp, STAT)

        !Arguments-------------------------------------------------------------
        integer                                         :: InterfaceID
        integer,            intent(IN)                  :: PropertyID
        real, optional,     intent(IN)                  :: DTProp
        real,              dimension(:,:  ), pointer    :: Concentration         
        integer,           dimension(:,:  ), pointer    :: WaterPoints2D
        integer, optional, dimension(:,:  ), pointer    :: OpenPoints2D
        real   , optional, dimension(:,:  ), pointer    :: Oxygen2D

        integer, optional,  intent(OUT)                 :: STAT

        !External--------------------------------------------------------------
        integer                                         :: ready_, STAT_CALL
        logical                                         :: ReadyToCompute

        !Local-----------------------------------------------------------------
        integer                                         :: STAT_
        integer                                         :: nProperty
        integer                                         :: Index 
        integer                                         :: prop
        integer                                         :: ILB, IUB, JLB, JUB, i, j  
        integer                                         :: PropLB, PropUB, ArrayLB, ArrayUB 
        real                                            :: DTProp_, DT
        logical                                         :: Increment

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

        STAT_ = UNKNOWN_

        call Ready(InterfaceID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

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

            PropLB  = Me%Prop%ILB
            PropUB  = Me%Prop%IUB

            ArrayLB = Me%Array%ILB
            ArrayUB = Me%Array%IUB

            Me%ExternalVar%WaterPoints2D        => WaterPoints2D
            
            if(present(OpenPoints2D))then
                Me%ExternalVar%OpenPoints2D     => OpenPoints2D
            end if

            if(present(Oxygen2D))then
                call UnfoldMatrix(Oxygen2D, Me%Oxygen)
            end if

            Increment = OFF
            ReadyToCompute = .false.

cd7 :       if (present(DTProp))then
                DTProp_   = DTProp
                Increment = ON  !If Increment is ON  then the concentrations is will be changed
            else
                DTProp_   = FillValueReal
                Increment = OFF !If Increment is OFF then the sinks and sources module will be computed
            end if cd7


cd5 :       if (.not. Increment) then 

                call FillMassTempSalinity(PropertyID, Concentration, ReadyToCompute)

cd4 :           if (ReadyToCompute) then

                    call UnfoldMatrix(Me%ExternalVar%OpenPoints2D, Me%OpenPoints)

                    !Stores the concentration before changing them
                    Me%ConcentrationIncrement = Me%Mass

                    select case (Me%SinksSourcesModel)

                        case(BenthicCEQUALW2Model)

                            call CEQUALW2Benthic(Me%ObjCEQUALW2,                       &
                                                 Me%Temperature,                       &
                                                 Me%Mass,                              &
                                                 Me%OpenPoints,                        &
                                                 STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface2D - ModuleInterface - ERR01'
                        
                        case(BenthosModel)

                            call ModifyBenthos  (Me%ObjBenthos,                        &
                                                 Me%Temperature,                       &
                                                 Me%Oxygen,                            &
                                                 Me%OpenPoints,                        &
                                                 Me%Mass,                              &
                                                 Me%Array%ILB,                         &
                                                 Me%Array%IUB,                         &
                                                 STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface2D - ModuleInterface - ERR02'

                    end select

do7 :               do prop  = PropLB,  PropUB
do6 :               do index = ArrayLB, ArrayUB
                        Me%ConcentrationIncrement(prop, index) = Me%Mass(prop, index) - &
                                                                 Me%ConcentrationIncrement(prop, index)
                    end do do6
                    end do do7
                
                end if cd4

            elseif (Increment) then cd5

                Index = 0

                nProperty = PropertyIndexNumber(PropertyID)

                select case (Me%SinksSourcesModel)

                    case(BenthicCEQUALW2Model)

                        call GetDTCEQUALW2(Me%ObjCEQUALW2, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface2D - ModuleInterface - ERR03'


                    case(BenthosModel)

                        call GetDTBenthos(Me%ObjBenthos, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface2D - ModuleInterface - ERR04'

                end select

                do j = JLB, JUB
                do i = ILB, IUB
                    if (Me%ExternalVar%WaterPoints2D(i, j) == 1) then
                        Index = Index + 1
                        
                        if (Me%ExternalVar%OpenPoints2D(i, j) == 1) then
                            Concentration(i, j) = Concentration(i, j)      + &
                            Me%ConcentrationIncrement(nProperty, Index) * DTProp / DT 
                        end if

                    end if
                end do
                end do


            end if cd5

            STAT_ = SUCCESS_
        else               
            STAT_ = ready_
        end if cd1


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

    end subroutine Modify_Interface2D
    
    !--------------------------------------------------------------------------

    subroutine Modify_Interface1D(InterfaceID, PropertyID, Concentration,         &
                                  RiverPoints1D, OpenPoints1D, DWZ,               &
                                  ShortWaveRadiation, LightExtCoefField,          &
                                  WaterPercentage, DissolvedToParticulate1D,      &
                                  DTProp, STAT)

        !Arguments-------------------------------------------------------------
        integer                                         :: InterfaceID
        integer, intent(IN)                             :: PropertyID
        real,              dimension(:), pointer        :: Concentration
        real,    optional, dimension(:), pointer        :: ShortWaveRadiation
        real,    optional, dimension(:), pointer        :: LightExtCoefField
        real   , optional, dimension(:), pointer        :: WaterPercentage
        real,    optional, dimension(:), pointer        :: DissolvedToParticulate1D
        integer,           dimension(:), pointer        :: RiverPoints1D
        integer, optional, dimension(:), pointer        :: OpenPoints1D
        real,    optional, dimension(:), pointer        :: DWZ
        real,    optional,  intent(IN)                  :: DTProp
        integer, optional,  intent(OUT)                 :: STAT

        !External--------------------------------------------------------------
        real                                            :: DT
        integer                                         :: STAT_CALL
        integer                                         :: ready_
        logical                                         :: ReadyToCompute

        !Local-----------------------------------------------------------------
        integer                                         :: STAT_
        integer                                         :: nProperty
        integer                                         :: Index 
        integer                                         :: i
        integer                                         :: prop, JulDay
        integer                                         :: ILB, IUB
        integer                                         :: PropLB, PropUB, ArrayLB, ArrayUB 
        real                                            :: DTProp_
        logical                                         :: Increment

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

        STAT_ = UNKNOWN_

        call Ready(InterfaceID, ready_)    

cd1 :   if (ready_ .EQ. IDLE_ERR_) then

            ILB     = Me%Size1D%ILB
            IUB     = Me%Size1D%IUB

            PropLB  = Me%Prop%ILB
            PropUB  = Me%Prop%IUB

            ArrayLB = Me%Array%ILB
            ArrayUB = Me%Array%IUB

            if(present(OpenPoints1D ))Me%ExternalVar%OpenPoints1D     => OpenPoints1D
            if(present(DWZ          ))Me%ExternalVar%DWZ1D            => DWZ
            
            Me%ExternalVar%RiverPoints1D => RiverPoints1D

            Increment = OFF
            ReadyToCompute = .false.

cd7 :       if (present(DTProp))then
                DTProp_   = DTProp
                Increment = ON  !If Increment is ON  then the concentrations is will be changed
            else
                DTProp_   = FillValueReal
                Increment = OFF !If Increment is OFF then the sinks and sources module will be computed
            end if cd7


cd5 :       if (.not. Increment) then 

                call FillMassTempSalinity(PropertyID, Concentration, ReadyToCompute)

cd4 :           if (ReadyToCompute) then

                    call UnfoldMatrix(Me%ExternalVar%OpenPoints1D, Me%OpenPoints)

                    !Stores the concentration before changing them
                    Me%ConcentrationIncrement = Me%Mass

                    select case (Me%SinksSourcesModel)

                        case(WaterQualityModel)

                            call UnfoldMatrix(ShortWaveRadiation,   Me%ShortWaveRadiation)
                            call UnfoldMatrix(LightExtCoefField,    Me%LightExtCoefField)
                            call UnfoldMatrix(Me%ExternalVar%DWZ1D, Me%Thickness)


                            call WaterQuality(Me%ObjWaterQuality,                   &
                                              Me%Salinity,                          &
                                              Me%Temperature,                       &
                                              Me%ShortWaveRadiation,                &
                                              Me%LightExtCoefField,                 &
                                              Me%Thickness,                         &
                                              Me%Mass,                              &
                                              Me%Array%ILB,                         &
                                              Me%Array%IUB,                         &
                                              Me%OpenPoints,                        &
                                              FishFood = Me%FishFood,               &
                                              STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR05'
                        
                        case(SedimentQualityModel)

                            call UnfoldMatrix(WaterPercentage, Me%WaterPercentage)
                            call UnfoldMatrix(DissolvedToParticulate1D, Me%DissolvedToParticulate)

                            call SedimentQuality(Me%ObjSedimentQuality,             &
                                                 Me%Temperature,                    &
                                                 Me%Mass,                           &
                                                 Me%WaterPercentage,                &
                                                 Me%DissolvedToParticulate,         &
                                                 Me%Array%ILB,                      &
                                                 Me%Array%IUB,                      &
                                                 STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR06'

                        case(LifeModel)

                            call UnfoldMatrix(ShortWaveRadiation,   Me%ShortWaveRadiation)
                            call UnfoldMatrix(LightExtCoefField,    Me%LightExtCoefField)
                            call UnfoldMatrix(Me%ExternalVar%DWZ1D, Me%Thickness)
                            
                            call GetComputeCurrentTime(Me%ObjTime,                  &
                                                       Me%ExternalVar%Now,          &
                                                       STAT = STAT_CALL)                    
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERRXX' 


                            call JulianDay(Me%ExternalVar%Now, JulDay)

                            call ModifyLife(Me%ObjLife,                     &
                                      Me%Salinity,                          &
                                      Me%Temperature,                       &
                                      Me%ShortWaveRadiation,                &
                                      Me%LightExtCoefField,                 &
                                      Me%Thickness,                         &
                                      Me%Mass,                              &
                                      Me%OpenPoints,                        &
                                      Me%Array,                             &
                                      JulDay,                               &
                                    STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR07'



                        case(CEQUALW2Model)

                            call UnfoldMatrix(ShortWaveRadiation,    Me%ShortWaveRadiation)
                            call UnfoldMatrix(LightExtCoefField,     Me%LightExtCoefField)
                            call UnfoldMatrix(Me%ExternalVar%DWZ1D,  Me%Thickness)
                            

                            call CEQUALW2(Me%ObjCEQUALW2,                       &
                                          Me%Salinity,                          &
                                          Me%Temperature,                       &
                                          Me%Alkalinity,                        &
                                          Me%ShortWaveRadiation,                &
                                          Me%LightExtCoefField,                 &
                                          Me%Thickness,                         &
                                          Me%Mass,                              &
                                          Me%OpenPoints,                        &
                                          STAT = STAT_CALL)
                            if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR08'
                        

                    end select

do7 :               do prop  = PropLB,  PropUB
do6 :               do index = ArrayLB, ArrayUB
                        Me%ConcentrationIncrement(prop, index) = Me%Mass(prop, index) - &
                                                                 Me%ConcentrationIncrement(prop, index)
                    end do do6
                    end do do7
                
                end if cd4

            elseif (Increment) then cd5

                Index = 0

                nProperty = PropertyIndexNumber(PropertyID)

                select case (Me%SinksSourcesModel)

                    case(WaterQualityModel)
                        
                        call GetDTWQM(Me%ObjWaterQuality, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR09'

do3 :                   do i = ILB, IUB
cd2 :                       if (Me%ExternalVar%RiverPoints1D(i) == 1) then
                                Index = Index + 1
                                !Concentrations are only actualized in OpenPoints because of instability
                                !in waterpoints that are not openpoints
                                if (Me%ExternalVar%OpenPoints1D(i) == 1) then
                                    Concentration(i) = Concentration( i)      + &
                                    Me%ConcentrationIncrement(nProperty, Index) * DTProp / DT 
                                end if
                            end if cd2
                        end do do3

                    case(SedimentQualityModel)
                        
                        call GetDTSedimentQuality(Me%ObjSedimentQuality, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR10'                

do10:                   do i = ILB, IUB
cd8 :                       if (Me%ExternalVar%RiverPoints1D(i) == 1) then
                                Index = Index + 1
                                Concentration(i) = Concentration( i)                     + &
                                                         Me%ConcentrationIncrement(nProperty, Index) * &
                                                         DTProp / DT 
                            end if cd8
                        end do do10



                    case(CEQUALW2Model)

                        call GetDTCEQUALW2(Me%ObjCEQUALW2, DTSecond = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR11'

do13 :                   do i = ILB, IUB
cd9 :                       if (Me%ExternalVar%RiverPoints1D(i) == 1) then
                                Index = Index + 1
                                !Concentrations are only actualized in OpenPoints because of instability
                                !in waterpoints that are not openpoints
                                if (Me%ExternalVar%OpenPoints1D(i) == 1) then
                                    Concentration(i) = Concentration( i)      + &
                                    Me%ConcentrationIncrement(nProperty, Index) * DTProp / DT 
                                end if
                            end if cd9
                        end do do13


                   
                    case(LifeModel)

                         call GetDTLife(Me%ObjLife, DT = DT, STAT = STAT_CALL)
                        if (STAT_CALL /= SUCCESS_) stop 'Modify_Interface1D - ModuleInterface - ERR12'

do17 :                   do i = ILB, IUB
cd14 :                       if (Me%ExternalVar%RiverPoints1D(i) == 1) then
                                Index = Index + 1
                                !Concentrations are only actualized in OpenPoints because of instability
                                !in waterpoints that are not openpoints
                                if (Me%ExternalVar%OpenPoints1D(i) == 1) then
                                    Concentration(i) = Concentration(i)      + &
                                    Me%ConcentrationIncrement(nProperty, Index) * DTProp / DT 
                                end if
                            end if cd14
                        end do do17



                end select

            end if cd5

            STAT_ = SUCCESS_
        else               
            STAT_ = ready_
        end if cd1


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

    end subroutine Modify_Interface1D

    !--------------------------------------------------------------------------
    
    subroutine FillMassTempSalinity3D(PropertyID, Concentration, Ready) 

        !Arguments-------------------------------------------------------------
        integer,                intent(IN )     :: PropertyID
        real, dimension(:,:,:), pointer         :: Concentration
        logical,                intent(OUT)     :: Ready

        !External--------------------------------------------------------------
        integer                                 :: numZoo, numPhyto, numDiatoms 
        integer                                 :: numSiBio, numSiDiss 
        integer                                 :: numAmmonia, numNitrate, numNitrite
        integer                                 :: numDONRefractory
        integer                                 :: numDONNonRefractory
        integer                                 :: numPartOrganicNitrogen
        integer                                 :: numOxygen, numBOD
        integer                                 :: numDOPRefractory, numDOPNonRefractory
        integer                                 :: numPartOrganicPhosphorus, numInorganicPhosphorus 
        integer                                 :: numBacteria, numCiliate
        integer                                 :: numLarvae
        integer                                 :: numHetrotrophicN, numHetrotrophicC
        integer                                 :: numAutotrotrophicN, numAutotrotrophicC
        integer                                 :: numAnaerobicN, numAnaerobicC        
        integer                                 :: numLabil_OM_N, numLabil_OM_C
        integer                                 :: numRefractOM_N, numRefractOM_C
        integer                                 :: numNgas 
        integer                                 :: STAT_CALL
        logical                                 :: lZoo, lPhyto, lDiatoms 
        logical                                 :: lNitrogen, lPhosphorus, lSilica 
        logical                                 :: lOxygen, lSalinity, lBOD
        logical                                 :: lBacteria, lCiliate
        logical                                 :: lLarvae
        logical                                 :: lCarbon

        !Local-----------------------------------------------------------------
        integer                                 :: i, PropLB, PropUB, IndexNumber 
        logical                                 :: TemperatureAdded = .false.
        logical                                 :: SalinityAdded = .false.
        logical                                 :: FishFoodAdded = .false.
        logical                                 :: AlkalinityAdded = .false.
        !----------------------------------------------------------------------

        PropLB  = Me%Prop%ILB
        PropUB  = Me%Prop%IUB

        select case (Me%SinksSourcesModel)

            case (WaterQualityModel)

                call GetWQPropIndex(Me%ObjWaterQuality,                                         &
                                     Zoo                             = numZoo,                  &
                                     Larvae                          = numLarvae,               &
                                     Phyto                           = numPhyto,                &
                                     Diatoms                         = numDiatoms,              &
                                     Ammonia                         = numAmmonia,              &
                                     Nitrate                         = numNitrate,              &
                                     Nitrite                         = numNitrite,              &
                                     BiogenicSilica                  = numSiBio,                & 
                                     DissolvedSilica                 = numSiDiss,               & 
                                     DissOrganicNitrogenRefractory   = numDONRefractory,        &
                                     DONNonRefractory                = numDONNonRefractory,     &
                                     PartOrganicNitrogen             = numPartOrganicNitrogen,  &
                                     Oxygen                          = numOxygen,               &
                                     BOD                             = numBOD,                  &
                                     DissOrganicPhosphorusRefractory = numDOPRefractory,        &
                                     DOPNonRefractory                = numDOPNonRefractory,     &
                                     PartOrganicPhosphorus           = numPartOrganicPhosphorus,&
                                     InorganicPhosphorus             = numInorganicPhosphorus,  &
                                     Bacteria                        = numBacteria,             &
                                     Ciliate                         = numCiliate,              &
                                     STAT                            = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity3D - ModuleInterface - ERR01' 

                call GetWQOptions(Me%ObjWaterQuality,                                           &
                                   Zoo            = lZoo,                                       &
                                   Larvae         = lLarvae,                                    &
                                   Phyto          = lPhyto,                                     &
                                   Diatoms        = lDiatoms,                                   & 
                                   Silica         = lSilica,                                    & 
                                   Nitrogen       = lNitrogen,                                  &
                                   Phosphorus     = lPhosphorus,                                &
                                   Oxygen         = lOxygen,                                    &
                                   Salinity       = lSalinity,                                  &
                                   BOD            = lBOD,                                       & 
                                   Bacteria       = lBacteria,                                  &
                                   Ciliate        = lCiliate,                                   &
                                   STAT           = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity3D - ModuleInterface - ERR02'


cd16 :          if (lPhyto) then 
cd17 :          if (PropertyID == Phytoplankton_) then
                    call InputData(Concentration,numPhyto)
                    Me%AddedProperties(numPhyto)    = .TRUE.
                end if cd17
                end if cd16

cd18 :          if (lZoo) then
cd19 :          if (PropertyID == Zooplankton_) then
                    call InputData(Concentration,numZoo)
                    Me%AddedProperties(numZoo)      = .TRUE.
                end if cd19
                end if cd18

cd118 :         if (lLarvae) then
cd119 :         if (PropertyID == Larvae_) then
                    call InputData(Concentration,numLarvae)
                    Me%AddedProperties(numLarvae)   = .TRUE.
                end if cd119
                end if cd118

cd20 :          if (lPhosphorus) then

                    if (PropertyID == POP_) then
                        call InputData(Concentration, numPartOrganicPhosphorus)
                        Me%AddedProperties(numPartOrganicPhosphorus)    = .TRUE.
                    end if

                    if (PropertyID== DOPRefractory_) then
                        call InputData(Concentration, numDOPRefractory)
                        Me%AddedProperties(numDOPRefractory)            = .TRUE.
                    end if

                    if (PropertyID== DOPNon_Refractory_) then
                        call InputData(Concentration, numDOPNonRefractory)
                        Me%AddedProperties(numDOPNonRefractory)         = .TRUE.
                    end if

cd22 :              if (PropertyID== Inorganic_Phosphorus_) then
                        call InputData(Concentration, numInorganicPhosphorus)
                        Me%AddedProperties(numInorganicPhosphorus)      = .TRUE.
                    end if cd22
                end if cd20


cd9 :           if (lNitrogen) then
cd10 :              if (PropertyID== PON_) then
                        call InputData(Concentration, numPartOrganicNitrogen)
                        Me%AddedProperties(numPartOrganicNitrogen)      = .TRUE.
                    end if cd10

cd11 :              if (PropertyID== DONRefractory_) then
                        call InputData(Concentration, numDONRefractory)
                        Me%AddedProperties(numDONRefractory)            = .TRUE.
                    end if cd11

cd12 :              if (PropertyID== DONNon_Refractory_) then
                        call InputData(Concentration,numDONNonRefractory)
                        Me%AddedProperties(numDONNonRefractory)         = .TRUE.
                    end if cd12

cd13 :              if (PropertyID== Ammonia_) then
                        call InputData(Concentration,numAmmonia)
                        Me%AddedProperties(numAmmonia)                  = .TRUE.
                    end if cd13

cd14 :              if (PropertyID== Nitrate_) then
                        call InputData(Concentration,numNitrate)
                        Me%AddedProperties(numNitrate)                  = .TRUE.
                    end if cd14

cd15 :              if (PropertyID== Nitrite_) then
                        call InputData(Concentration,numNitrite)
                        Me%AddedProperties(numNitrite)                  = .TRUE.
                    end if cd15
                end if cd9


cd7 :           if (lBOD) then
cd8 :           if (PropertyID== BOD_) then
                    call InputData(Concentration,numBOD)
                    Me%AddedProperties(numBOD)      = .TRUE.
                end if cd8
                end if cd7



cd23 :          if (lBacteria) then
cd24 :          if (PropertyID== Bacteria_) then
                    call InputData(Concentration,numBacteria)
                    Me%AddedProperties(numBacteria) = .TRUE.
                end if cd24
                end if cd23



cd25 :          if (lCiliate) then
cd26 :          if (PropertyID== Ciliate_) then
                    call InputData(Concentration,numCiliate)
                    Me%AddedProperties(numCiliate)  = .TRUE.
                end if cd26
                end if cd25


                !The oxygen is always compute. If it's not defined in the waterproperties 
                ! module then the saturation value is compute
cd6 :           if (PropertyID== Oxygen_) then
                    call InputData(Concentration,numOxygen)
                    Me%AddedProperties(numOxygen)   = .TRUE.
                end if cd6

                if (lSalinity) then
                
cd4 :               if (PropertyID== Salinity_) then
                        call UnfoldMatrix(Concentration, Me%Salinity)
                        SalinityAdded       =.TRUE.
                    end if cd4
                end if

cd3 :           if (PropertyID== Temperature_) then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded    =.TRUE.
                end if cd3

cd99:           if (PropertyID== FishFood_) then
                    call UnfoldMatrix(Concentration, Me%FishFood)
                    FishFoodAdded       =.TRUE.
                end if cd99

cd116 :         if (lDiatoms) then 
cd117 :         if (PropertyID == Diatoms_) then
                    call InputData(Concentration,numDiatoms)
                    Me%AddedProperties(numDiatoms)    = .TRUE.
                end if cd117
                end if cd116


cd150 :          if (lSilica) then

                    if (PropertyID == DSilica_) then
                        call InputData(Concentration, numSiDiss)
                        Me%AddedProperties(numSiDiss)    = .TRUE.
                    end if

                    if (PropertyID== BioSilica_) then
                        call InputData(Concentration, numSiBio)
                        Me%AddedProperties(numSiBio)            = .TRUE.
                    end if

                end if cd150




cd2 :           if ((SalinityAdded .or. (.not. lSalinity)) .AND. TemperatureAdded) then
                    Ready = .TRUE.

do1 :               do i = PropLB, PropUB
cd1 :               if (.NOT. Me%AddedProperties(i)) then
                        Ready = .FALSE.
                        exit do1
                    end if cd1
                    end do do1

                    if (Ready) Me%AddedProperties = .FALSE.
                end if cd2        
        
        case (SedimentQualityModel)

            call GetPropIndex(  Me%ObjSedimentQuality,                                      & 
                                HetrotrophicN                   = numHetrotrophicN,         &
                                HetrotrophicC                   = numHetrotrophicC,         &
                                AutotrotrophicN                 = numAutotrotrophicN,       &
                                AutotrotrophicC                 = numAutotrotrophicC,       &
                                AnaerobicN                      = numAnaerobicN,            &
                                AnaerobicC                      = numAnaerobicC,            &
                                Labil_OM_C                      = numLabil_OM_C,            &
                                Labil_OM_N                      = numLabil_OM_N,            &
                                RefractOM_C                     = numRefractOM_C,           &
                                RefractOM_N                     = numRefractOM_N,           &
                                Ammonia                         = numAmmonia,               &
                                Nitrate                         = numNitrate,               &
                                Ngas                            = numNgas,                  &
                                Oxygen                          = numOxygen,                &
                                STAT                            = STAT_CALL )          
             if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity3D - ModuleInterface - ERR03' 


             call GetSQOptions( Me%ObjSedimentQuality,                                      &
                                Nitrogen                        = lNitrogen,                &
                                Carbon                          = lCarbon,                  &
                                Oxygen                          = lOxygen,                  &
                                STAT                            = STAT_CALL  )

             if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity - ModuleInterface - ERR04' 



cd27 :          if (lNitrogen) then

cd28 :              if (PropertyID== PON_                  ) then
                        call InputData(Concentration, numLabil_OM_N)
                        Me%AddedProperties(numLabil_OM_N)         = .TRUE.
                    end if cd28

cd29 :              if (PropertyID== RefreactaryOrganicN_  ) then
                        call InputData(Concentration, numRefractOM_N)
                        Me%AddedProperties(numRefractOM_N)        = .TRUE.
                    end if cd29

cd30 :              if (PropertyID== Ammonia_              ) then
                        call InputData(Concentration,numAmmonia)
                        Me%AddedProperties(numAmmonia)            = .TRUE.
                    end if cd30

cd31 :              if (PropertyID== Nitrate_              ) then
                        call InputData(Concentration,numNitrate)
                        Me%AddedProperties(numNitrate)            = .TRUE.
                    end if cd31

cd32 :              if (PropertyID== Ngas_                 ) then
                        call InputData(Concentration,numNgas)
                        Me%AddedProperties(numNgas)               = .TRUE.
                    end if cd32

cd33 :              if (PropertyID== HeterotrophicN_       ) then
                        call InputData(Concentration,numHetrotrophicN)
                        Me%AddedProperties(numHetrotrophicN)      = .TRUE.
                    end if cd33

cd34 :              if (PropertyID== AutotrophicN_         ) then
                        call InputData(Concentration,numAutotrotrophicN)
                        Me%AddedProperties(numAutotrotrophicN)    = .TRUE.
                    end if cd34

cd35 :              if (PropertyID== AnaerobicN_           ) then
                        call InputData(Concentration,numAnaerobicN)
                        Me%AddedProperties(numAnaerobicN)         = .TRUE.
                    end if cd35

                
                end if cd27


cd36 :          if (lCarbon) then

cd37 :              if (PropertyID== LabileOrganicC_       ) then
                        call InputData(Concentration, numLabil_OM_C)
                        Me%AddedProperties(numLabil_OM_C)         = .TRUE.
                    end if cd37

cd38 :              if (PropertyID== RefreactaryOrganicC_  ) then
                        call InputData(Concentration, numRefractOM_C)
                        Me%AddedProperties(numRefractOM_C)        = .TRUE.
                    end if cd38

cd39 :              if (PropertyID== HeterotrophicC_       ) then
                        call InputData(Concentration,numHetrotrophicC)
                        Me%AddedProperties(numHetrotrophicC)      = .TRUE.
                    end if cd39

cd40 :              if (PropertyID== AnaerobicC_           ) then
                        call InputData(Concentration,numAnaerobicC)
                        Me%AddedProperties(numAnaerobicC)         = .TRUE.
                    end if cd40

cd41 :              if (PropertyID== AutotrophicC_         ) then
                        call InputData(Concentration,numAutotrotrophicC)
                        Me%AddedProperties(numAutotrotrophicC)    = .TRUE.
                    end if cd41

                end if cd36

                !The oxygen is always compute. If it's not defined in the waterproperties 
                ! module then the saturation value is compute
cd42 :          if (PropertyID== Oxygen_)then
                    call InputData(Concentration,numOxygen)
                    Me%AddedProperties(numOxygen)                 = .TRUE.
                end if cd42


cd43 :          if (PropertyID== Temperature_)then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded =.TRUE.
                end if cd43


cd44 :          if (TemperatureAdded) then
                    Ready = .TRUE.

do2  :              do i = PropLB, PropUB
cd45 :                  if (.NOT. Me%AddedProperties(i)) then
                            Ready = .FALSE.
                            exit do2
                    end if cd45
                    end do do2
                    if (Ready) Me%AddedProperties = .FALSE.
                
                end if cd44

            case(LifeModel)
                call GetLifePropIndex(Me%ObjLife,PropertyID,IndexNumber,STAT=STAT_CALL)
                if (STAT_CALL==SUCCESS_) then
                    call InputData(Concentration, IndexNumber)
                    Me%AddedProperties(IndexNumber) =.TRUE.
                else if (STAT_CALL.NE.SUCCESS_ .AND. STAT_CALL.NE.NOT_FOUND_ERR_) then
                    stop 'FillMassTempSalinity3D - ModuleInterface - ERR05' 
                end if
                     
                               
                if (PropertyID== Salinity_) then
                    call UnfoldMatrix(Concentration, Me%Salinity)
                    SalinityAdded       =.TRUE.
                end if 

                if (PropertyID== Temperature_) then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded    =.TRUE.
                end if 

       
                if (SalinityAdded .AND. TemperatureAdded) then
                    Ready = .TRUE.

                    do i = PropLB, PropUB
                    if (.NOT. Me%AddedProperties(i)) then
                            Ready = .FALSE.
                            exit 
                    end if 
                    end do 
                    if (Ready) Me%AddedProperties = .FALSE.
                end if



            case(CEQUALW2Model) 
                
                call GetCEQUALW2PropIndex(Me%ObjCEQUALW2,PropertyID,IndexNumber,STAT=STAT_CALL)
                if (STAT_CALL==SUCCESS_) then
                    call InputData(Concentration, IndexNumber)
                    Me%AddedProperties(IndexNumber) =.TRUE.
                else if (STAT_CALL.NE.SUCCESS_ .AND. STAT_CALL.NE.NOT_FOUND_ERR_) then
                    stop 'FillMassTempSalinity3D - ModuleInterface - ERR05' 
                end if
                     
                               
                if (PropertyID== Salinity_) then
                    call UnfoldMatrix(Concentration, Me%Salinity)
                    SalinityAdded       =.TRUE.
                end if 

                if (PropertyID== Temperature_) then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded    =.TRUE.
                end if 

                if (PropertyID== Alkalinity_) then
                    call UnfoldMatrix(Concentration, Me%Alkalinity)
                    AlkalinityAdded       =.TRUE.
                end if 

                if (SalinityAdded .AND. TemperatureAdded .AND. AlkalinityAdded) then
                    Ready = .TRUE.

                    do i = PropLB, PropUB
                    if (.NOT. Me%AddedProperties(i)) then
                            Ready = .FALSE.
                            exit 
                    end if 
                    end do 
                    if (Ready) Me%AddedProperties = .FALSE.
                end if

        end select

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

    end Subroutine FillMassTempSalinity3D

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

    subroutine FillMassTempSalinity2D(PropertyID, Concentration, Ready) 

        !Arguments-------------------------------------------------------------
        integer,                intent(IN )     :: PropertyID
        real, dimension(:,:),   pointer         :: Concentration
        logical,                intent(OUT)     :: Ready

        !External--------------------------------------------------------------

        !Local-----------------------------------------------------------------
        integer                                 :: PropLB, PropUB, IndexNumber, STAT_CALL
        logical                                 :: TemperatureAdded = .false.
        integer                                 :: i
        !----------------------------------------------------------------------

        PropLB  = Me%Prop%ILB
        PropUB  = Me%Prop%IUB

        select case (Me%SinksSourcesModel)

            case(BenthicCeQualW2Model)

                call GetCEQUALW2PropIndex(Me%ObjCEQUALW2,PropertyID,IndexNumber,STAT=STAT_CALL)

            case(BenthosModel)

                call GetBenthosPropIndex(Me%ObjBenthos, PropertyID, IndexNumber, STAT = STAT_CALL)

        end select
        
        if (STAT_CALL==SUCCESS_) then
            call InputData(Concentration, IndexNumber)
            Me%AddedProperties(IndexNumber) =.TRUE.
        else if (STAT_CALL.NE.SUCCESS_ .AND. STAT_CALL.NE.NOT_FOUND_ERR_) then
            stop 'FillMassTempSalinity3D - ModuleInterface - ERR05' 
        end if

        if (PropertyID== Temperature_) then
            call UnfoldMatrix(Concentration, Me%Temperature)
            TemperatureAdded    =.TRUE.
        end if  

        if (TemperatureAdded) then
            Ready = .TRUE.

            do i = PropLB, PropUB
                if (.NOT. Me%AddedProperties(i)) then
                    Ready = .FALSE.
                    exit 
                end if 
            end do 
            if (Ready) Me%AddedProperties = .FALSE.

        end if


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

    end Subroutine FillMassTempSalinity2D

    !--------------------------------------------------------------------------
    
    subroutine FillMassTempSalinity1D(PropertyID, Concentration, Ready) 

        !Arguments-------------------------------------------------------------
        integer,            intent(IN )         :: PropertyID
        real, dimension(:), pointer             :: Concentration
        logical,            intent(OUT)         :: Ready

        !External--------------------------------------------------------------
        integer                                 :: numZoo, numPhyto, numDiatoms 
        integer                                 :: numSiBio, numSiDiss 
        integer                                 :: numAmmonia, numNitrate, numNitrite
        integer                                 :: numDONRefractory
        integer                                 :: numDONNonRefractory
        integer                                 :: numPartOrganicNitrogen
        integer                                 :: numOxygen, numBOD
        integer                                 :: numDOPRefractory, numDOPNonRefractory
        integer                                 :: numPartOrganicPhosphorus, numInorganicPhosphorus 
        integer                                 :: numBacteria, numCiliate
        integer                                 :: numLarvae
        integer                                 :: numHetrotrophicN, numHetrotrophicC
        integer                                 :: numAutotrotrophicN, numAutotrotrophicC
        integer                                 :: numAnaerobicN, numAnaerobicC        
        integer                                 :: numLabil_OM_N, numLabil_OM_C
        integer                                 :: numRefractOM_N, numRefractOM_C
        integer                                 :: numNgas 
        integer                                 :: STAT_CALL
        logical                                 :: lZoo, lPhyto, lDiatoms 
        logical                                 :: lNitrogen, lPhosphorus, lSilica 
        logical                                 :: lOxygen, lSalinity, lBOD
        logical                                 :: lBacteria, lCiliate
        logical                                 :: lLarvae
        logical                                 :: lCarbon

        !Local-----------------------------------------------------------------
        integer                                 :: i, PropLB, PropUB, IndexNumber 
        logical                                 :: TemperatureAdded
        logical                                 :: SalinityAdded
        logical                                 :: FishFoodAdded
        logical                                 :: AlkalinityAdded
        !----------------------------------------------------------------------

        PropLB  = Me%Prop%ILB
        PropUB  = Me%Prop%IUB

        select case (Me%SinksSourcesModel)

            case (WaterQualityModel)

                call GetWQPropIndex(Me%ObjWaterQuality,                                         &
                                     Zoo                             = numZoo,                  &
                                     Larvae                          = numLarvae,               &
                                     Phyto                           = numPhyto,                &
                                     Diatoms                         = numDiatoms,              & 
                                     Ammonia                         = numAmmonia,              &
                                     Nitrate                         = numNitrate,              &
                                     Nitrite                         = numNitrite,              &
                                     BiogenicSilica                  = numSiBio,                & 
                                     DissolvedSilica                 = numSiDiss,               & 
                                     DissOrganicNitrogenRefractory   = numDONRefractory,        &
                                     DONNonRefractory                = numDONNonRefractory,     &
                                     PartOrganicNitrogen             = numPartOrganicNitrogen,  &
                                     Oxygen                          = numOxygen,               &
                                     BOD                             = numBOD,                  &
                                     DissOrganicPhosphorusRefractory = numDOPRefractory,        &
                                     DOPNonRefractory                = numDOPNonRefractory,     &
                                     PartOrganicPhosphorus           = numPartOrganicPhosphorus,&
                                     InorganicPhosphorus             = numInorganicPhosphorus,  &
                                     Bacteria                        = numBacteria,             &
                                     Ciliate                         = numCiliate,              &
                                     STAT                            = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity3D - ModuleInterface - ERR01' 

                call GetWQOptions(Me%ObjWaterQuality,                                           &
                                   Zoo            = lZoo,                                       &
                                   Larvae         = lLarvae,                                    &
                                   Phyto          = lPhyto,                                     &
                                   Diatoms        = lDiatoms,                                   & 
                                   Silica         = lSilica,                                    & 
                                   Nitrogen       = lNitrogen,                                  &
                                   Phosphorus     = lPhosphorus,                                &
                                   Oxygen         = lOxygen,                                    &
                                   Salinity       = lSalinity,                                  &
                                   BOD            = lBOD,                                       & 
                                   Bacteria       = lBacteria,                                  &
                                   Ciliate        = lCiliate,                                   &
                                   STAT           = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity3D - ModuleInterface - ERR02'


cd16 :          if (lPhyto) then 
cd17 :          if (PropertyID == Phytoplankton_) then
                    call InputData(Concentration,numPhyto)
                    Me%AddedProperties(numPhyto)    = .TRUE.
                end if cd17
                end if cd16

cd18 :          if (lZoo) then
cd19 :          if (PropertyID == Zooplankton_) then
                    call InputData(Concentration,numZoo)
                    Me%AddedProperties(numZoo)      = .TRUE.
                end if cd19
                end if cd18

cd118 :         if (lLarvae) then
cd119 :         if (PropertyID == Larvae_) then
                    call InputData(Concentration,numLarvae)
                    Me%AddedProperties(numLarvae)   = .TRUE.
                end if cd119
                end if cd118

cd20 :          if (lPhosphorus) then

                    if (PropertyID == POP_) then
                        call InputData(Concentration, numPartOrganicPhosphorus)
                        Me%AddedProperties(numPartOrganicPhosphorus)    = .TRUE.
                    end if

                    if (PropertyID== DOPRefractory_) then
                        call InputData(Concentration, numDOPRefractory)
                        Me%AddedProperties(numDOPRefractory)            = .TRUE.
                    end if

                    if (PropertyID== DOPNon_Refractory_) then
                        call InputData(Concentration, numDOPNonRefractory)
                        Me%AddedProperties(numDOPNonRefractory)         = .TRUE.
                    end if

cd22 :              if (PropertyID== Inorganic_Phosphorus_) then
                        call InputData(Concentration, numInorganicPhosphorus)
                        Me%AddedProperties(numInorganicPhosphorus)      = .TRUE.
                    end if cd22
                end if cd20


cd9 :           if (lNitrogen) then
cd10 :              if (PropertyID== PON_) then
                        call InputData(Concentration, numPartOrganicNitrogen)
                        Me%AddedProperties(numPartOrganicNitrogen)      = .TRUE.
                    end if cd10

cd11 :              if (PropertyID== DONRefractory_) then
                        call InputData(Concentration, numDONRefractory)
                        Me%AddedProperties(numDONRefractory)            = .TRUE.
                    end if cd11

cd12 :              if (PropertyID== DONNon_Refractory_) then
                        call InputData(Concentration,numDONNonRefractory)
                        Me%AddedProperties(numDONNonRefractory)         = .TRUE.
                    end if cd12

cd13 :              if (PropertyID== Ammonia_) then
                        call InputData(Concentration,numAmmonia)
                        Me%AddedProperties(numAmmonia)                  = .TRUE.
                    end if cd13

cd14 :              if (PropertyID== Nitrate_) then
                        call InputData(Concentration,numNitrate)
                        Me%AddedProperties(numNitrate)                  = .TRUE.
                    end if cd14

cd15 :              if (PropertyID== Nitrite_) then
                        call InputData(Concentration,numNitrite)
                        Me%AddedProperties(numNitrite)                  = .TRUE.
                    end if cd15
                end if cd9


cd7 :           if (lBOD) then
cd8 :           if (PropertyID== BOD_) then
                    call InputData(Concentration,numBOD)
                    Me%AddedProperties(numBOD)      = .TRUE.
                end if cd8
                end if cd7



cd23 :          if (lBacteria) then
cd24 :          if (PropertyID== Bacteria_) then
                    call InputData(Concentration,numBacteria)
                    Me%AddedProperties(numBacteria) = .TRUE.
                end if cd24
                end if cd23



cd25 :          if (lCiliate) then
cd26 :          if (PropertyID== Ciliate_) then
                    call InputData(Concentration,numCiliate)
                    Me%AddedProperties(numCiliate)  = .TRUE.
                end if cd26
                end if cd25


                !The oxygen is always compute. If it's not defined in the waterproperties 
                ! module then the saturation value is compute
cd6 :           if (PropertyID== Oxygen_) then
                    call InputData(Concentration,numOxygen)
                    Me%AddedProperties(numOxygen)   = .TRUE.
                end if cd6

                if (lSalinity) then
                
cd4 :               if (PropertyID== Salinity_) then
                        call UnfoldMatrix(Concentration, Me%Salinity)
                        SalinityAdded       =.TRUE.
                    end if cd4
                end if

cd3 :           if (PropertyID== Temperature_) then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded    =.TRUE.
                end if cd3

cd99:           if (PropertyID== FishFood_) then
                    call UnfoldMatrix(Concentration, Me%FishFood)
                    FishFoodAdded       =.TRUE.
                end if cd99


cd116 :         if (lDiatoms) then 
cd117 :         if (PropertyID == Diatoms_) then
                    call InputData(Concentration,numDiatoms)
                    Me%AddedProperties(numDiatoms)    = .TRUE.
                end if cd117
                end if cd116


cd150 :          if (lSilica) then

                    if (PropertyID == DSilica_) then
                        call InputData(Concentration, numSiDiss)
                        Me%AddedProperties(numSiBio)    = .TRUE.
                    end if

                    if (PropertyID== BioSilica_) then
                        call InputData(Concentration, numSiBio)
                        Me%AddedProperties(numSiBio)            = .TRUE.
                    end if

                end if cd150


cd2 :           if ((SalinityAdded .or. (.not. lSalinity)) .AND. TemperatureAdded) then
                    Ready = .TRUE.

do1 :               do i = PropLB, PropUB
cd1 :               if (.NOT. Me%AddedProperties(i)) then
                        Ready = .FALSE.
                        exit do1
                    end if cd1
                    end do do1

                    if (Ready) Me%AddedProperties = .FALSE.
                end if cd2        
        
        case (SedimentQualityModel)

            call GetPropIndex(  Me%ObjSedimentQuality,                                      & 
                                HetrotrophicN                   = numHetrotrophicN,         &
                                HetrotrophicC                   = numHetrotrophicC,         &
                                AutotrotrophicN                 = numAutotrotrophicN,       &
                                AutotrotrophicC                 = numAutotrotrophicC,       &
                                AnaerobicN                      = numAnaerobicN,            &
                                AnaerobicC                      = numAnaerobicC,            &
                                Labil_OM_C                      = numLabil_OM_C,            &
                                Labil_OM_N                      = numLabil_OM_N,            &
                                RefractOM_C                     = numRefractOM_C,           &
                                RefractOM_N                     = numRefractOM_N,           &
                                Ammonia                         = numAmmonia,               &
                                Nitrate                         = numNitrate,               &
                                Ngas                            = numNgas,                  &
                                Oxygen                          = numOxygen,                &
                                STAT                            = STAT_CALL )          
             if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity3D - ModuleInterface - ERR03' 


             call GetSQOptions( Me%ObjSedimentQuality,                                      &
                                Nitrogen                        = lNitrogen,                &
                                Carbon                          = lCarbon,                  &
                                Oxygen                          = lOxygen,                  &
                                STAT                            = STAT_CALL  )

             if (STAT_CALL .NE. SUCCESS_)stop 'FillMassTempSalinity - ModuleInterface - ERR04' 



cd27 :          if (lNitrogen) then

cd28 :              if (PropertyID== PON_                  ) then
                        call InputData(Concentration, numLabil_OM_N)
                        Me%AddedProperties(numLabil_OM_N)         = .TRUE.
                    end if cd28

cd29 :              if (PropertyID== RefreactaryOrganicN_  ) then
                        call InputData(Concentration, numRefractOM_N)
                        Me%AddedProperties(numRefractOM_N)        = .TRUE.
                    end if cd29

cd30 :              if (PropertyID== Ammonia_              ) then
                        call InputData(Concentration,numAmmonia)
                        Me%AddedProperties(numAmmonia)            = .TRUE.
                    end if cd30

cd31 :              if (PropertyID== Nitrate_              ) then
                        call InputData(Concentration,numNitrate)
                        Me%AddedProperties(numNitrate)            = .TRUE.
                    end if cd31

cd32 :              if (PropertyID== Ngas_                 ) then
                        call InputData(Concentration,numNgas)
                        Me%AddedProperties(numNgas)               = .TRUE.
                    end if cd32

cd33 :              if (PropertyID== HeterotrophicN_       ) then
                        call InputData(Concentration,numHetrotrophicN)
                        Me%AddedProperties(numHetrotrophicN)      = .TRUE.
                    end if cd33

cd34 :              if (PropertyID== AutotrophicN_         ) then
                        call InputData(Concentration,numAutotrotrophicN)
                        Me%AddedProperties(numAutotrotrophicN)    = .TRUE.
                    end if cd34

cd35 :              if (PropertyID== AnaerobicN_           ) then
                        call InputData(Concentration,numAnaerobicN)
                        Me%AddedProperties(numAnaerobicN)         = .TRUE.
                    end if cd35

                
                end if cd27


cd36 :          if (lCarbon) then

cd37 :              if (PropertyID== LabileOrganicC_       ) then
                        call InputData(Concentration, numLabil_OM_C)
                        Me%AddedProperties(numLabil_OM_C)         = .TRUE.
                    end if cd37

cd38 :              if (PropertyID== RefreactaryOrganicC_  ) then
                        call InputData(Concentration, numRefractOM_C)
                        Me%AddedProperties(numRefractOM_C)        = .TRUE.
                    end if cd38

cd39 :              if (PropertyID== HeterotrophicC_       ) then
                        call InputData(Concentration,numHetrotrophicC)
                        Me%AddedProperties(numHetrotrophicC)      = .TRUE.
                    end if cd39

cd40 :              if (PropertyID== AnaerobicC_           ) then
                        call InputData(Concentration,numAnaerobicC)
                        Me%AddedProperties(numAnaerobicC)         = .TRUE.
                    end if cd40

cd41 :              if (PropertyID== AutotrophicC_         ) then
                        call InputData(Concentration,numAutotrotrophicC)
                        Me%AddedProperties(numAutotrotrophicC)    = .TRUE.
                    end if cd41

                end if cd36

                !The oxygen is always compute. If it's not defined in the waterproperties 
                ! module then the saturation value is compute
cd42 :          if (PropertyID== Oxygen_)then
                    call InputData(Concentration,numOxygen)
                    Me%AddedProperties(numOxygen)                 = .TRUE.
                end if cd42


cd43 :          if (PropertyID== Temperature_)then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded =.TRUE.
                end if cd43


cd44 :          if (TemperatureAdded) then
                    Ready = .TRUE.

do2  :              do i = PropLB, PropUB
cd45 :                  if (.NOT. Me%AddedProperties(i)) then
                            Ready = .FALSE.
                            exit do2
                    end if cd45
                    end do do2
                    if (Ready) Me%AddedProperties = .FALSE.
                
                end if cd44

            case(LifeModel)
                call GetLifePropIndex(Me%ObjLife,PropertyID,IndexNumber,STAT=STAT_CALL)
                if (STAT_CALL==SUCCESS_) then
                    call InputData(Concentration, IndexNumber)
                    Me%AddedProperties(IndexNumber) =.TRUE.
                else if (STAT_CALL.NE.SUCCESS_ .AND. STAT_CALL.NE.NOT_FOUND_ERR_) then
                    stop 'FillMassTempSalinity3D - ModuleInterface - ERR05' 
                end if
                     
                               
                if (PropertyID== Salinity_) then
                    call UnfoldMatrix(Concentration, Me%Salinity)
                    SalinityAdded       =.TRUE.
                end if 

                if (PropertyID== Temperature_) then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded    =.TRUE.
                end if 

       
                if (SalinityAdded .AND. TemperatureAdded) then
                    Ready = .TRUE.

                    do i = PropLB, PropUB
                    if (.NOT. Me%AddedProperties(i)) then
                            Ready = .FALSE.
                            exit 
                    end if 
                    end do 
                    if (Ready) Me%AddedProperties = .FALSE.
                end if



            case(CEQUALW2Model) 
                
                call GetCEQUALW2PropIndex(Me%ObjCEQUALW2,PropertyID,IndexNumber,STAT=STAT_CALL)
                if (STAT_CALL==SUCCESS_) then
                    call InputData(Concentration, IndexNumber)
                    Me%AddedProperties(IndexNumber) =.TRUE.
                else if (STAT_CALL.NE.SUCCESS_ .AND. STAT_CALL.NE.NOT_FOUND_ERR_) then
                    stop 'FillMassTempSalinity3D - ModuleInterface - ERR05' 
                end if
                     
                               
                if (PropertyID== Salinity_) then
                    call UnfoldMatrix(Concentration, Me%Salinity)
                    SalinityAdded       =.TRUE.
                end if 

                if (PropertyID== Temperature_) then
                    call UnfoldMatrix(Concentration, Me%Temperature)
                    TemperatureAdded    =.TRUE.
                end if 

                if (PropertyID== Alkalinity_) then
                    call UnfoldMatrix(Concentration, Me%Alkalinity)
                    AlkalinityAdded       =.TRUE.
                end if 

                if (SalinityAdded .AND. TemperatureAdded .AND. AlkalinityAdded) then
                    Ready = .TRUE.

                    do i = PropLB, PropUB
                    if (.NOT. Me%AddedProperties(i)) then
                            Ready = .FALSE.
                            exit 
                    end if 
                    end do 
                    if (Ready) Me%AddedProperties = .FALSE.
                end if

        end select

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

    end subroutine FillMassTempSalinity1D
    
    !----------------------------------------------------------------------

   
    subroutine InputData3D (Concentration, nProperty)

        !Arguments-------------------------------------------------------------
        real, dimension(:,:,:), pointer         :: Concentration
        integer, intent (in)                    :: nProperty 

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

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

        ILB = Me%Size3D%ILB     
        IUB = Me%Size3D%IUB     

        JLB = Me%Size3D%JLB    
        JUB = Me%Size3D%JUB    

        KLB = Me%Size3D%KLB   
        KUB = Me%Size3D%KUB   
        
        !Number indexed to 3D cell in the vector 
        Index = 0

        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB
            if (Me%ExternalVar%WaterPoints3D(i,j,k)==1) then
                Index = Index + 1
                Me%Mass(nProperty,Index) = Concentration(i,j,k)
            endif    
        enddo
        enddo
        enddo

        if ((Index)> Me%Array%IUB) stop 'InputData3D - ModuleInterface - ERR01'

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

    end subroutine InputData3D

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

    subroutine InputData2D (Concentration, nProperty)

        !Arguments-------------------------------------------------------------
        real, dimension(:,:), pointer           :: Concentration
        integer, intent (in)                    :: nProperty 

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

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

        ILB = Me%Size2D%ILB     
        IUB = Me%Size2D%IUB     

        JLB = Me%Size2D%JLB    
        JUB = Me%Size2D%JUB    

        
        !Number indexed to 3D cell in the vector 
        Index = 0

        do j = JLB, JUB
        do i = ILB, IUB
            if (Me%ExternalVar%WaterPoints2D(i,j)==1) then
                Index = Index + 1
                Me%Mass(nProperty,Index) = Concentration(i,j)
            endif    
        enddo
        enddo

        if ((Index)> Me%Array%IUB) stop 'InputData2D - ModuleInterface - ERR01'

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

    end subroutine InputData2D

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

    subroutine InputData1D (Concentration, nProperty)

        !Arguments-------------------------------------------------------------
        real, dimension(:), pointer             :: Concentration
        integer, intent (in)                    :: nProperty 

        !Local-----------------------------------------------------------------
        integer                                 :: Index
        integer                                 :: i
        integer                                 :: ILB, IUB      

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

        ILB = Me%Size1D%ILB     
        IUB = Me%Size1D%IUB     

        !Number indexed to 1D cell in the vector 
        Index = 0

        do i = ILB, IUB
            if (Me%ExternalVar%RiverPoints1D(i) == WaterPoint) then
                Index = Index + 1
                Me%Mass(nProperty,Index) = Concentration(i)
            endif    
        enddo

        if ((Index)> Me%Array%IUB) stop 'InputData1D - ModuleInterface - ERR01'

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

    end subroutine InputData1D

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

    subroutine UnfoldMatrix3D_R (Matrix3D, Vector)

        !Arguments-------------------------------------------------------------
        real, dimension(:,:,:), pointer     :: Matrix3D
        real, dimension(:    ), pointer     :: Vector

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

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

        ILB = Me%Size3D%ILB     
        IUB = Me%Size3D%IUB     

        JLB = Me%Size3D%JLB    
        JUB = Me%Size3D%JUB    

        KLB = Me%Size3D%KLB   
        KUB = Me%Size3D%KUB   
        
        !Number indexed to 3D cell in the vector 
        Index = 0

        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB
            if (Me%ExternalVar%WaterPoints3D(i,j,k)==1) then
                Index         = Index + 1
                Vector(Index) = Matrix3D(i,j,k)
            endif    
        enddo
        enddo
        enddo

        if ((Index) > Me%Array%IUB) stop 'UnfoldMatrix3D_R - ModuleInterface - ERR01'
            
        !----------------------------------------------------------------------

    end subroutine UnfoldMatrix3D_R

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

    subroutine UnfoldMatrix3D_I (Matrix3D, Vector, Vertical1D)

        !Arguments-------------------------------------------------------------
        integer, dimension(:,:,:), pointer      :: Matrix3D
        integer, dimension(:    ), pointer      :: Vector
        logical, optional, intent(IN)           :: Vertical1D

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

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

        ILB = Me%Size3D%ILB     
        IUB = Me%Size3D%IUB     

        JLB = Me%Size3D%JLB    
        JUB = Me%Size3D%JUB    

        KLB = Me%Size3D%KLB   
        KUB = Me%Size3D%KUB   

        if (present(Vertical1D)) then
            Vertical1D_Aux = Vertical1D
        else
            Vertical1D_Aux = .false.
        endif
        
        !Number indexed to 3D cell in the vector 
        Index = 0

        do k = KLB, KUB
        do j = JLB, JUB
        do i = ILB, IUB
            if (Me%ExternalVar%WaterPoints3D(i,j,k)==1) then
                Index         = Index + 1
                Vector(Index) = Matrix3D(i,j,k)
                if (Vertical1D_Aux .and. .not.(i==2 .and. j==2)) Vector(Index) = 0
            endif    
        enddo
        enddo
        enddo

        if ((Index) > Me%Array%IUB) stop 'UnfoldMatrix3D_I - ModuleInterface - ERR01'
            

    end subroutine UnfoldMatrix3D_I

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

    subroutine UnfoldMatrix2D_R (Matrix2D, Vector)

        !Arguments-------------------------------------------------------------
        real, dimension(:,:), pointer       :: Matrix2D
        real, dimension(:  ), pointer       :: Vector

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

        ILB = Me%Size2D%ILB     
        IUB = Me%Size2D%IUB     

        JLB = Me%Size2D%JLB    
        JUB = Me%Size2D%JUB    

          
        
        !Number indexed to 3D cell in the vector 
        Index = 0

        do j = JLB, JUB
        do i = ILB, IUB
            if (Me%ExternalVar%WaterPoints2D(i,j)==1) then

                    Index         = Index + 1
                    
             if (Me%ExternalVar%OpenPoints2D(i,j)==1) then
                
                Vector(Index) = Matrix2D(i,j)
             
             endif
                    
                    
            endif    
        enddo
        enddo

        if ((Index) > Me%Array%IUB) stop 'UnfoldMatrix2D_R - ModuleInterface - ERR01'
            

    end subroutine UnfoldMatrix2D_R
        
    !----------------------------------------------------------------------

    subroutine UnfoldMatrix2D_I (Matrix2D, Vector)

        !Arguments-------------------------------------------------------------
        integer, dimension(:,:), pointer        :: Matrix2D
        integer, dimension(:  ), pointer        :: Vector

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

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

        ILB = Me%Size2D%ILB     
        IUB = Me%Size2D%IUB     

        JLB = Me%Size2D%JLB    
        JUB = Me%Size2D%JUB    
        
        !Number indexed to 3D cell in the vector 
        Index = 0

        do j = JLB, JUB
        do i = ILB, IUB
            if (Me%ExternalVar%WaterPoints2D(i,j)==1) then
             Index         = Index + 1
             if (Me%ExternalVar%OpenPoints2D(i,j)==1) then
                
                Vector(Index) = Matrix2D(i,j)
             
             endif
            endif    
        enddo
        enddo

        if ((Index) > Me%Array%IUB) stop 'UnfoldMatrix2D_I - ModuleInterface - ERR01'
            

    end subroutine UnfoldMatrix2D_I

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

    subroutine UnfoldMatrix1D_R (Matrix1D, Vector)

        !Arguments-------------------------------------------------------------
        real, dimension(:  ), pointer       :: Matrix1D
        real, dimension(:  ), pointer       :: Vector

        !Local-----------------------------------------------------------------
        integer                             :: Index
        integer                             :: i
        integer                             :: ILB, IUB      
           
        !----------------------------------------------------------------------

        ILB = Me%Size1D%ILB     
        IUB = Me%Size1D%IUB     

        !Number indexed to 1D cell in the vector 
        Index = 0

        do i = ILB, IUB
            if (Me%ExternalVar%RiverPoints1D(i)==1) then

                    Index         = Index + 1
                    
             if (Me%ExternalVar%OpenPoints1D(i)==1) then
                
                Vector(Index) = Matrix1D(i)
             
             endif
                    
                    
            endif    
        enddo

        if ((Index) > Me%Array%IUB) stop 'UnfoldMatrix1D_R - ModuleInterface - ERR01'
            

    end subroutine UnfoldMatrix1D_R
        
    !----------------------------------------------------------------------

    subroutine UnfoldMatrix1D_I (Matrix1D, Vector)

        !Arguments-------------------------------------------------------------
        integer, dimension(:  ), pointer        :: Matrix1D
        integer, dimension(:  ), pointer        :: Vector

        !Local-----------------------------------------------------------------
        integer                                 :: Index
        integer                                 :: i
        integer                                 :: ILB, IUB      
        !----------------------------------------------------------------------

        ILB = Me%Size1D%ILB     
        IUB = Me%Size1D%IUB     

        !Number indexed to 3D cell in the vector 
        Index = 0

        do i = ILB, IUB
            if (Me%ExternalVar%RiverPoints1D(i)==1) then
             Index         = Index + 1
             if (Me%ExternalVar%OpenPoints1D(i)==1) then
                
                Vector(Index) = Matrix1D(i)
             
             endif
            endif    
        enddo

        if ((Index) > Me%Array%IUB) stop 'UnfoldMatrix1D_I - ModuleInterface - ERR01'
            

    end subroutine UnfoldMatrix1D_I

    !--------------------------------------------------------------------------
        
    integer function PropertyIndexNumber(PropertyID)

        !Arguments-------------------------------------------------------------
        integer                         :: PropertyID

        !External--------------------------------------------------------------
        integer                         :: STAT_CALL
        integer                         :: numZoo, numPhyto, numDiatoms 
        integer                         :: numLarvae
        integer                         :: numAmmonia, numNitrate, numNitrite
        integer                         :: numSiBio, numSiDiss 
        integer                         :: numDONRefractory, numDONNonRefractory
        integer                         :: numPartOrganicNitrogen
        integer                         :: numOxygen, numBOD
        integer                         :: numDOPRefractory, numDOPNonRefractory
        integer                         :: numPartOrganicPhosphorus, numInorganicPhosphorus 
        integer                         :: numBacteria, numCiliate

        integer                         :: numHetrotrophicN, numHetrotrophicC
        integer                         :: numAutotrotrophicN, numAutotrotrophicC
        integer                         :: numAnaerobicN, numAnaerobicC        
        integer                         :: numLabil_OM_N, numLabil_OM_C
        integer                         :: numRefractOM_N, numRefractOM_C
        integer                         :: numNgas
        
        !Local-----------------------------------------------------------------
        integer                         :: nProperty
        logical                         :: CheckName

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

        select case (Me%SinksSourcesModel)

            case (WaterQualityModel)
        
                call GetWQPropIndex(Me%ObjWaterQuality,                                  &
                             Zoo                              = numZoo,                  &
                             Larvae                           = numLarvae,               &
                             Phyto                            = numPhyto,                &
                             Diatoms                          = numDiatoms,              & 
                             Ammonia                          = numAmmonia,              &
                             Nitrate                          = numNitrate,              &
                             Nitrite                          = numNitrite,              &
                             BiogenicSilica                   = numSiBio,                & 
                             DissolvedSilica                  = numSiDiss,               & 
                             DissOrganicNitrogenRefractory    = numDONRefractory,        &
                             DONNonRefractory                 = numDONNonRefractory,     &
                             PartOrganicNitrogen              = numPartOrganicNitrogen,  &
                             Oxygen                           = numOxygen,               &
                             BOD                              = numBOD,                  &
                             DissOrganicPhosphorusRefractory  = numDOPRefractory,        &
                             DOPNonRefractory                 = numDOPNonRefractory,     &
                             PartOrganicPhosphorus            = numPartOrganicPhosphorus,&
                             InorganicPhosphorus              = numInorganicPhosphorus,  &
                             Bacteria                         = numBacteria,             &
                             Ciliate                          = numCiliate,              &
                             STAT                             = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_) stop 'PropertyIndexNumber - ModuleInterface - ERR01'
                   
   
cd1 :           if      (PropertyID== Phytoplankton_       ) then
                    nProperty = numPhyto

                else if (PropertyID== Zooplankton_         ) then cd1
                    nProperty = numZoo
 
                else if (PropertyID== Larvae_              ) then cd1
                    nProperty = numLarvae

                else if (PropertyID== POP_                 ) then cd1
                    nProperty = numPartOrganicPhosphorus

                else if (PropertyID== DOPRefractory_       ) then cd1
                    nProperty = numDOPRefractory

                else if (PropertyID== DOPNon_Refractory_   ) then cd1
                    nProperty = numDOPNonRefractory

                else if (PropertyID== Inorganic_Phosphorus_) then cd1
                    nProperty = numInorganicPhosphorus

                else if (PropertyID== PON_                 ) then cd1
                    nProperty = numPartOrganicNitrogen

                else if (PropertyID== DONRefractory_       ) then cd1
                    nProperty = numDONRefractory

                else if (PropertyID== DONNon_Refractory_   ) then cd1
                    nProperty = numDONNonRefractory 

                else if (PropertyID== Ammonia_             ) then cd1
                    nProperty = numAmmonia

                else if (PropertyID== Nitrate_             ) then cd1
                    nProperty = numNitrate

                else if (PropertyID== Nitrite_             ) then cd1
                    nProperty = numNitrite

                else if (PropertyID== Bacteria_            ) then cd1
                    nProperty = numBacteria

                else if (PropertyID== Ciliate_             ) then cd1
                    nProperty = numCiliate

                else if (PropertyID== BOD_                 ) then cd1
                    nProperty = numBOD

                else if (PropertyID== Oxygen_              ) then cd1
                    nProperty = numOxygen
                
                else if (PropertyID== Diatoms_            ) then cd1
                    nProperty = numDiatoms 

                else if (PropertyID== DSilica_            ) then cd1
                    nProperty = numSiDiss 

                else if (PropertyID== BioSilica_           ) then cd1
                    nProperty = numSiBio

                else if (PropertyID== GrossProd_           ) then cd1 
                    CheckName = CheckPropertyName('grossprod', number = nProperty)
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR03'

                else if (PropertyID== NutrientLim_         ) then cd1   
                    CheckName = CheckPropertyName('nutrientlim', number = nProperty)     
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR04'

                else if (PropertyID== NLim_                ) then cd1   
                    CheckName = CheckPropertyName('nitrogenlim', number = nProperty)     
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR05'

                else if (PropertyID== PLim_                ) then cd1   
                    CheckName = CheckPropertyName('phosphoruslim', number = nProperty)     
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR06'

                else if (PropertyID== LightLim_            ) then cd1   
                    CheckName = CheckPropertyName('lightlim', number = nProperty) 
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR07'

                else if (PropertyID== TemperatureLim_      ) then cd1 
                    CheckName = CheckPropertyName('temperaturelim', number = nProperty)    
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR08'

                else if (PropertyID== DiaGrossProd_           ) then cd1 
                    CheckName = CheckPropertyName('diagrossprod', number = nProperty)
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR09'

                else if (PropertyID== DiaNutrientLim_         ) then cd1   
                    CheckName = CheckPropertyName('dianutrientlim', number = nProperty)     
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR10'

                else if (PropertyID== DiaNLim_         ) then cd1   
                    CheckName = CheckPropertyName('dianitrogenlim', number = nProperty)     
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR11'

                else if (PropertyID== DiaPLim_         ) then cd1   
                    CheckName = CheckPropertyName('diaphosphoruslim', number = nProperty)     
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR12'

                else if (PropertyID== DiaSiLim_         ) then cd1   
                    CheckName = CheckPropertyName('diasilicalim', number = nProperty)     
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR13'

                else if (PropertyID== DiaLightLim_            ) then cd1   
                    CheckName = CheckPropertyName('dialightlim', number = nProperty) 
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR14'

                else if (PropertyID== DiaTemperatureLim_      ) then cd1 
                    CheckName = CheckPropertyName('diatemperaturelim', number = nProperty)    
                    if (.NOT.CheckName) stop 'PropertyIndexNumber - ModuleInterface - ERR15'



                else  cd1
                    write(*,*) 
                    write(*,*) 'Inconsistency between Interface and WaterQuality.'
                    stop       'PropertyIndexNumber - ModuleInterface - ERR16'
                end if cd1


            case (SedimentQualityModel)
                
                call GetPropIndex(  Me%ObjSedimentQuality,                                  & 
                                    HetrotrophicN                   = numHetrotrophicN,     &
                                    HetrotrophicC                   = numHetrotrophicC,     &
                                    AutotrotrophicN                 = numAutotrotrophicN,   &
                                    AutotrotrophicC                 = numAutotrotrophicC,   &
                                    AnaerobicN                      = numAnaerobicN,        &
                                    AnaerobicC                      = numAnaerobicC,        &
                                    Labil_OM_C                      = numLabil_OM_C,        &
                                    Labil_OM_N                      = numLabil_OM_N,        &
                                    RefractOM_C                     = numRefractOM_C,       &
                                    RefractOM_N                     = numRefractOM_N,       &
                                    Ammonia                         = numAmmonia,           &
                                    Nitrate                         = numNitrate,           &
                                    Ngas                            = numNgas,              &
                                    Oxygen                          = numOxygen,            &
                                    STAT                            = STAT_CALL )          
                if (STAT_CALL .NE. SUCCESS_) stop 'PropertyIndexNumber - ModuleInterface - ERR17'


                if      (PropertyID== PON_                 ) then            
                    nProperty = numLabil_OM_N
            
                else if (PropertyID== RefreactaryOrganicN_ ) then
                    nProperty = numRefractOM_N

                else if (PropertyID== Ammonia_             ) then
                    nProperty = numAmmonia

                else if (PropertyID== Nitrate_             ) then
                    nProperty = numNitrate

                else if (PropertyID== Ngas_                ) then
                    nProperty = numNgas

                else if (PropertyID== HeterotrophicN_      ) then
                    nProperty = numHetrotrophicN

                else if (PropertyID== AutotrophicN_        ) then
                    nProperty = numAutotrotrophicN

                else if (PropertyID== AnaerobicN_          ) then
                    nProperty = numAnaerobicN

                else if (PropertyID== LabileOrganicC_      ) then
                    nProperty = numLabil_OM_C

                else if (PropertyID== RefreactaryOrganicC_ ) then
                    nProperty = numRefractOM_C

                else if (PropertyID== HeterotrophicC_      ) then
                    nProperty = numHetrotrophicC

                else if (PropertyID== AnaerobicC_          ) then
                    nProperty = numAnaerobicC

                else if (PropertyID== AutotrophicC_        ) then
                    nProperty = numAutotrotrophicC

                else if (PropertyID== Oxygen_              ) then
                    nProperty = numOxygen

                else
                    write(*,*) 
                    write(*,*) 'Inconsistency between Interface and SedimentQuality.'
                    stop       'PropertyIndexNumber - ModuleInterface - ERR18'

                end if

            case(LifeModel)
                call GetLifePropIndex(Me%ObjLife,PropertyID,nProperty,STAT_CALL)
                
                if (STAT_CALL.NE.SUCCESS_) then 
                    
                    write(*,*) 
                    write(*,*) 'Inconsistency between Interface and Life.'
                    stop       'PropertyIndexNumber - ModuleInterface - ERR19'

                end if  


            case(CEQUALW2Model, BenthicCEQUALW2Model)

                call GetCEQUALW2PropIndex(Me%ObjCEQUALW2,PropertyID,nProperty,STAT_CALL)
                
                if (STAT_CALL.NE.SUCCESS_) then 
                    
                    write(*,*) 
                    write(*,*) 'Inconsistency between Interface and CEQUALW2.'
                    stop       'PropertyIndexNumber - ModuleInterface - ERR20'

                end if  

            case(BenthosModel)
                
                call GetBenthosPropIndex(Me%ObjBenthos, PropertyID, nProperty, STAT_CALL)
                
                if (STAT_CALL.NE.SUCCESS_) then 
                    
                    write(*,*) 
                    write(*,*) 'Inconsistency between Interface and Benthos Model.'
                    stop       'PropertyIndexNumber - ModuleInterface - ERR21'

                end if  

                
        end select

        PropertyIndexNumber = nProperty

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

    end Function PropertyIndexNumber

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


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

    !DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR DESTRUCTOR  

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



    Subroutine KillInterface(InterfaceID, STAT)

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

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

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

        STAT_ = UNKNOWN_

        call Ready(InterfaceID, ready_) 

cd1 :   if (ready_ .NE. OFF_ERR_) then

            nUsers = DeassociateInstance(mINTERFACE_, Me%InstanceID)

            if (nUsers == 0) then

                nUsers = DeassociateInstance (mTIME_,Me%ObjTime)
                if (nUsers == 0) stop 'KillInterface - ModuleInterface - ERR01'

                select case(Me%SinksSourcesModel)
                    
                    case(WaterQualityModel)

                        call KillWaterQuality(Me%ObjWaterQuality, STAT = STAT_CALL)
                        if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR02'
                    
                    case(SedimentQualityModel)

                        call KillSedimentQuality(Me%ObjSedimentQuality, STAT = STAT_CALL)
                        if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR03'

                    case(LifeModel)

                        call KillLife(Me%ObjLife, STAT = STAT_CALL)
                        if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR04.1'

                    case(CEQUALW2Model, BenthicCEQUALW2Model)

                        call KillCEQUALW2(Me%ObjCEQUALW2, STAT = STAT_CALL)
                        if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR04.2'

                    case(BenthosModel)
                        
                        call KillBenthos(Me%ObjBenthos, STAT = STAT_CALL)
                        if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR04.3'

                end select

                if(associated(Me%Salinity))then
                    deallocate(Me%Salinity, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR05'
                end if
                
                if(associated(Me%Temperature))then
                    deallocate(Me%Temperature, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR06'
                end if

                if(associated(Me%Oxygen))then
                    deallocate(Me%Oxygen, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR06A'
                end if                
                
                if(associated(Me%ShortWaveRadiation))then
                    deallocate(Me%ShortWaveRadiation, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR07'
                end if
                
                
                if(associated(Me%LightExtCoefField))then
                    deallocate(Me%LightExtCoefField, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR08'
                end if
                
                if(associated(Me%Thickness))then
                    deallocate(Me%Thickness, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR09'
                end if
                
                deallocate(Me%ConcentrationIncrement, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR10'

                deallocate(Me%OpenPoints, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR11'
                
                
                if(associated(Me%FishFood))then
                    deallocate(Me%FishFood, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR12'
                end if

                if(associated(Me%Alkalinity))then
                    deallocate(Me%Alkalinity, STAT = STAT_CALL)
                    if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR13'
                end if
                
                
                deallocate(Me%Mass, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR14'


                deallocate(Me%AddedProperties, STAT = STAT_CALL)
                if (STAT_CALL .NE. SUCCESS_) stop 'KillInterface - ModuleInterface - ERR15'

                !Deallocates Instance
                call DeallocateInstance

                InterfaceID = 0

                STAT_ = SUCCESS_

            end if

        else 
            STAT_ = ready_
        end if cd1

        if (present(STAT)) STAT = STAT_


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

    end subroutine KillInterface


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

    subroutine DeallocateInstance

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

        !Local-----------------------------------------------------------------
        type (T_Interface), pointer           :: AuxInterface
        type (T_Interface), pointer           :: PreviousInterface

        !Updates pointers
        if (Me%InstanceID == FirstObjInterface%InstanceID) then
            FirstObjInterface => FirstObjInterface%Next
        else
            PreviousInterface => FirstObjInterface
            AuxInterface      => FirstObjInterface%Next
            do while (AuxInterface%InstanceID /= Me%InstanceID)
                PreviousInterface => AuxInterface
                AuxInterface      => AuxInterface%Next
            enddo

            !Now update linked list
            PreviousInterface%Next => AuxInterface%Next

        endif

        !Deallocates instance
        deallocate (Me)
        nullify    (Me)

            
    end subroutine DeallocateInstance





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

    !MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEMENT MANAGEME

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


    subroutine Ready (InterfaceID, ready_) 

        !Arguments-------------------------------------------------------------
        integer                            :: InterfaceID
        integer                            :: ready_

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

        nullify (Me)

cd1:    if (InterfaceID > 0) then
            call LocateObjInterface(InterfaceID)
            ready_ = VerifyReadLock (mINTERFACE_, Me%InstanceID)
        else
            ready_ = OFF_ERR_
        end if cd1

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

    end subroutine Ready

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

    subroutine LocateObjInterface (InterfaceID)

        !Arguments-------------------------------------------------------------
        integer                            :: InterfaceID

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

        Me => FirstObjInterface
        do while (associated (Me))
            if (Me%InstanceID == InterfaceID) exit
            Me => Me%Next
        enddo

        if (.not. associated(Me))                                        &
            stop 'ModuleInterface - LocateObjInterface - ERR01'

    end subroutine LocateObjInterface

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

end Module ModuleInterface

