SimulationCore2/FluidFlow/Well_Pressure_Data_Transfer.f90

SUBROUTINE WellPressureDataTransfer
    
    !! This subroutine calculates pressure distribution in well
    !! We divide well hydraulic system to 4 zones. 1- pump to string zone; this zone is horizontal
    !! 2- string zone; this zone includes drill pipes and collars and is vertical or neraly vertical
    !! 3- annulus zone; this zone includes annulus space between string and openhole or casing
    !! 4- Openhole zone; this zone includes spaces under bit to bottom of the wellbore
    !! The last two zone includes annulus space and open hole, are influenced by kicks and have crucial role in 
    !! computing pressure distribution because reference point (back pressure in well head open condition 
    !! and bottomhole pressure in well head closed condition) are in one of them. because of this, pressure
    !! distribution calculates in two subroutines : PressureAnnAndOHDistribution and PressurePumptoBitDistribution
    !! pressure change through well consist of pressure change due to friction, change of elevation and acceleration
    !! frictional pressure drop calculated by subroutine Sub_Press_Drop_Calculator
    !! This subroutine includes hydrostatic pressure and later will include effect of corss-section area change
    !! which causes acceleration or decceleration of flow
    
    USE FricPressDropVars
    use MudSystemModule
    USE MudSystemVARIABLES
    USE GeoElements_FluidModule
    USE Fluid_Flow_Startup_Vars
    USE KickVariables
    USE CMudPropertiesVariables
    USE CBopStackVariables
    USE sROP_Variables
    USE CHOKEVARIABLES, PI2 => PI
    USE CChokeManifoldVariables
    USE CBitProblemsVariables
    USE CChokeProblemsVariables
    USE CStringConfigurationVariables
    USE CDrillWatchVariables
    USE CDataDisplayConsoleVariables    
   
    
    IMPLICIT NONE
    
    INTEGER :: i , j
    REAL :: InstantaneousTotalOpenChokeArea , VolumeDensityProduct
    
    ChokeIsClosing = .FALSE.
    Kchoke = 0.0
    KBit = 0.0
    
    
    !!!!!!!!!!!!!!!! bit is present
    BitTrue = StringConfiguration%StringConfigurations(1)%ComponentType  ==  Bit_ComponentType           ! Bit_ComponentType = 0
    BitCd = 0.98
    BitNozzDia = 32.0 * StringConfiguration%BitDefinition%BitNozzleSize       ! nozzle diameter in 1/32 in
    BitNozzleArea = 7.6699E-4 * BitNozzDia**2
    BitNozzleNum = StringConfiguration%BitDefinition%BitNozzleNo - BitJetsPlugged * BitProblems%PlugJetsCount
    BitTotNozzArea = BitNozzleNum * BitNozzleArea + 0.5 * BitJetsWashedOut * BitProblems%JetWashoutCount * BitNozzleArea
    IF (BitNozzleNum == 0) THEN
        BitTotallyPluged = .TRUE.
    ELSE
        BitTotallyPluged = .FALSE.
    END IF
    
    FloatValveIn = StringConfiguration%BitDefinition%FloatValve
    !FloatValveOpen = .TRUE.
    !IF (JetWashoutCount > 0 .OR. PlugJetsCount > 0) THEN
    !    WRITE (*,*) 'PlugJetsCount' , PlugJetsCount
    !    WRITE (*,*) 'JetWashoutCount', JetWashoutCount
    !END IF
    
    !WRITE (*,*) 'BitNozzleNum', BitNozzleNum
    !WRITE (*,*) 'BitTotNozzArea', BitTotNozzArea
    !WRITE (*,*) 'Float Valve??', BitDefinition%FloatValve
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    IF(ALLOCATED(GasPocketFlowEl))      DEALLOCATE(GasPocketFlowEl)
    
!!!!!!!!!!!!!!!!!!!!  Well Head Condition ( Open or Closed )
    ChokeKroneckerDelta = 0

    
    IF (MudSystem%WellToPitsOpen == .FALSE. .AND. MudSystem%WellToChokeManifoldOpen == .FALSE.) THEN
        
        WellHeadOpen = .FALSE.
        !WRITE (*,*) ' Well Head is closed '
    ELSE 
        WellHeadOpen = .TRUE.
    END IF
    
    IF (MudSystem%WellToPitsOpen == .FALSE. .AND. MudSystem%WellToChokeManifoldOpen)        ChokeKroneckerDelta = 1
    !WRITE (*,*) ' WelltoPitsOpen=' , WelltoPitsOpen, 'WellToChokeManifoldOpen= ' , WellToChokeManifoldOpen
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    
    NoHorizontalEl = MudSystem%NoHorizontalMudElements
    NoStringEl = MudSystem%NoStringMudElements
    NoAnnulusEl = MudSystem%NoCasingMudElements
    NoWellToChokeEl = 0
    IF (ChokeKroneckerDelta == 1)       NoWellToChokeEl = MudSystem%ChokeLine_Density%Length()
    NoOpenHoleEl = MudSystem%NoBottomHoleMudElements
    NumbEl = NoHorizontalEl + NoStringEl + NoAnnulusEl + NoWellToChokeEl + NoOpenHoleEl
    !WRITE (*,*) 'H, S, A, Ch, O', NoHorizontalEl , NoStringEl , NoAnnulusEl , NoWellToChokeEl , NoOpenHoleEl
    IF(ALLOCATED(FlowEl))       DEALLOCATE(FlowEl)
    ALLOCATE(FlowEl(NumbEl))
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  Transfering Data from module:MudSystemVARIABLES  to  module:FricPressDropVars to calculate Pressure distribution in flow path
    !!!!!!!!!!!!!!!!!!!!!!! Horizontal line
    !WRITE (*,*) 'H', NoHorizontalEl
    DO i = 1 , NoHorizontalEl
        FlowEl(i)%StartX = MudSystem%Xstart_MudElement%Array(i) 
        FlowEl(i)%EndX = MudSystem%Xend_MudElement%Array(i)
        FlowEl(i)%StartTVD = 0.d0
        FlowEl(i)%EndTVD = 0.d0
        FlowEl(i)%Id = 0.0
        FlowEl(i)%Od = MudSystem%PipeOD_MudElement%Array(i)
        FlowEl(i)%density = MudSystem%Density_MudElement%Array(i)
        
        FlowEl(i)%Dhyd = FlowEl(i)%Od
        FlowEl(i)%Length = ABS(FlowEl(i)%EndX - FlowEl(i)%StartX)
        FlowEl(i)%DepthDiff = 0.d0
        FlowEl(i)%Volume = PI / 4.0 * FlowEl(i)%Od**2 * REAL(FlowEl(i)%Length) / Convfttoinch**2
        IF (MudSystem%MudType_MudElement%Array(i) == 0) THEN           ! = 0 for mud, = 1 for gas kick, = 4 for air
            FlowEl(i)%MaterialType = 0
        ELSE IF (MudSystem%MudType_MudElement%Array(i) == 4) THEN
            FlowEl(i)%MaterialType = 4
        ELSE
            FlowEl(i)%MaterialType = 1
        END IF
        
        FlowEl(i)%dPdLgrav = 0.0
        FlowEl(i)%StaticPressDiff = 0.0

        !WRITE (*,*) 'H density, length, Type of ith element' , i, FlowEl(i)%density , FlowEl(i)%length, FlowEl(i)%MaterialType

        
    END DO
        !WRITE (*,*) 'inlet and outlet mud theta600' , FlowEl(1)%Theta600 , FlowEl(NumbEl)%Theta600
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    !!!!!!!!!!!!!!!!!!!!!!!!! String
    !WRITE (*,*) 'S', NoStringEl
    StringFirstEl = NoHorizontalEl + 1
    StringLastEl = NoHorizontalEl + NoStringEl
    DO i = NoHorizontalEl + 1 , NoHorizontalEl + NoStringEl
        !WRITE (*,*) 'ST , i' , i
        FlowEl(i)%StartX = MudSystem%Xstart_MudElement%Array(i)
        !WRITE (*,*) 'StartX', FlowEl(i)%StartX
        FlowEl(i)%EndX = MudSystem%Xend_MudElement%Array(i)
        !WRITE (*,*) 'EndX', FlowEl(i)%EndX
        FlowEl(i)%StartTVD = MudSystem%TVDstart_MudElement%Array(i)
        !WRITE (*,*) 'StartTVD', FlowEl(i)%StartTVD
        FlowEl(i)%EndTVD = MudSystem%TVDend_MudElement%Array(i)
        !WRITE (*,*) 'EndTVD', FlowEl(i)%EndTVD
        FlowEl(i)%Id = 0.0
        !WRITE (*,*) 'Id', FlowEl(i)%Id
        FlowEl(i)%Od = MudSystem%PipeOD_MudElement%Array(i)
        !WRITE (*,*) 'Od', FlowEl(i)%Od
        FlowEl(i)%Density = MudSystem%Density_MudElement%Array(i)
        !WRITE (*,*) 'Density', FlowEl(i)%Density
        
        
        FlowEl(i)%Dhyd = FlowEl(i)%Od
        FlowEl(i)%Length = ABS(FlowEl(i)%EndX - FlowEl(i)%StartX)
        !WRITE (*,*) 'Length', FlowEl(i)%Length
        FlowEl(i)%DepthDiff = ABS(FlowEl(i)%StartTVD - FlowEl(i)%EndTVD)
        !WRITE (*,*) 'DepthDiff', FlowEl(i)%DepthDiff
        FlowEl(i)%Area = PI / 4.0 * FlowEl(i)%Od**2 / Convfttoinch**2
        FlowEl(i)%Volume = FlowEl(i)%Area * REAL(FlowEl(i)%Length)
        IF (MudSystem%MudType_MudElement%Array(i) == 0) THEN           ! = 0 for mud, = 1 for gas kick, = 4 for air
            FlowEl(i)%MaterialType = 0
        ELSE IF (MudSystem%MudType_MudElement%Array(i) == 4) THEN
            FlowEl(i)%MaterialType = 4
        ELSE
            FlowEl(i)%MaterialType = 1
        END IF
        
        IF (FlowEl(i)%MaterialType == 1 .OR. FlowEl(i)%MaterialType == 4) THEN
            FlowEl(i)%dPdLgrav = 0.0
        ELSE
            FlowEl(i)%dPdLgrav = 0.052 * FlowEl(i)%Density
        END IF
        FlowEl(i)%StaticPressDiff = FlowEl(i)%dPdLgrav * REAL(FlowEl(i)%DepthDiff)
        
        !WRITE (*,*) 'S density, length, DeltaPStatic' , i, FlowEl(i)%density , FlowEl(i)%length, FlowEl(i)%StaticPressDiff, FlowEl(i)%MaterialType


        
    END DO
            !WRITE (*,*) 'inlet and outlet mud theta600' , FlowEl(1)%Theta600 , FlowEl(NumbEl)%Theta600
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!! Annulus
    !WRITE (*,*) 'A', NoAnnulusEl
    AnnulusFirstEl = NoHorizontalEl + NoStringEl + 1
    AnnulusLastEl = NoHorizontalEl + NoStringEl + NoAnnulusEl
    DO i = NoHorizontalEl + NoStringEl + 1 , NoHorizontalEl + NoStringEl + NoAnnulusEl
        FlowEl(i)%StartX = MudSystem%Xstart_MudElement%Array(i)
        FlowEl(i)%EndX = MudSystem%Xend_MudElement%Array(i)
        FlowEl(i)%StartTVD = MudSystem%TVDstart_MudElement%Array(i)
        FlowEl(i)%EndTVD = MudSystem%TVDend_MudElement%Array(i)
        FlowEl(i)%Id = MudSystem%PipeID_MudElement%Array(i)
        FlowEl(i)%Od = MudSystem%PipeOD_MudElement%Array(i)
        FlowEl(i)%Density = MudSystem%Density_MudElement%Array(i)        
        
        FlowEl(i)%Dhyd = FlowEl(i)%Od - FlowEl(i)%Id
        FlowEl(i)%Length = ABS(FlowEl(i)%EndX - FlowEl(i)%StartX)
        FlowEl(i)%DepthDiff = ABS(FlowEl(i)%StartTVD - FlowEl(i)%EndTVD)
        FlowEl(i)%Area = PI / 4.0 * (FlowEl(i)%Od**2 - FlowEl(i)%Id**2) / Convfttoinch**2
        FlowEl(i)%Volume = FlowEl(i)%Area * REAL(FlowEl(i)%Length)
        IF (MudSystem%MudType_MudElement%Array(i) == 0) THEN           ! = 0 for mud, = 1 for gas kick, = 4 for air
            FlowEl(i)%MaterialType = 0
        ELSE IF (MudSystem%MudType_MudElement%Array(i) == 4) THEN
            FlowEl(i)%MaterialType = 4
        ELSE
            FlowEl(i)%MaterialType = 1
        END IF
        
        IF (FlowEl(i)%MaterialType == 1 .OR. FlowEl(i)%MaterialType == 4) THEN
            FlowEl(i)%dPdLgrav = 0.0
        ELSE
            FlowEl(i)%dPdLgrav = 0.052 * FlowEl(i)%density
        END IF
        FlowEl(i)%StaticPressDiff = FlowEl(i)%dPdLgrav * REAL(FlowEl(i)%DepthDiff)

        !WRITE (*,*) 'A density, length, Type of ith element' , i, FlowEl(i)%density , FlowEl(i)%length, FlowEl(i)%MaterialType

        
    END DO
            !WRITE (*,*) 'inlet and outlet mud theta600' , FlowEl(1)%Theta600 , FlowEl(NumbEl)%Theta600
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    !!!!!!!!!!!!!!! Well to choke manifold path
    !ChokeFirstEl = AnnulusLastEl + 1
    !ChokeLastEl = OpenholeFirstEl - 1
    j = 1
    DO i = NoHorizontalEl + NoStringEl + NoAnnulusEl + 1 , NoHorizontalEl + NoStringEl + NoAnnulusEl + NoWelltoChokeEl
        FlowEl(i)%StartX = MudSystem%ChokeLine_Mud_Backhead_X%Array(j)
        FlowEl(i)%EndX = MudSystem%ChokeLine_Mud_Forehead_X%Array(j)
        FlowEl(i)%StartTVD = 0.d0
        FlowEl(i)%EndTVD = 0.d0
        FlowEl(i)%Id = 0.0
        FlowEl(i)%Od = BopStackSpecification%ChokeLineId
        FlowEl(i)%density = MudSystem%ChokeLine_Density%Array(j)
            
        FlowEl(i)%Dhyd = FlowEl(i)%Od
        FlowEl(i)%Length = ABS(FlowEl(i)%EndX - FlowEl(i)%StartX)
        FlowEl(i)%DepthDiff = 0.d0
        FlowEl(i)%Area = PI / 4.0 * FlowEl(i)%Od**2 / Convfttoinch**2
        FlowEl(i)%Volume = FlowEl(i)%Area * REAL(FlowEl(i)%Length)
        IF (MudSystem%ChokeLine_MudOrKick%Array(j) == 0 .OR. MudSystem%ChokeLine_MudOrKick%Array(j) == 4) THEN      ! = 0 for mud, = 1 for gas kick, = 4 for air
            FlowEl(i)%MaterialType = MudSystem%ChokeLine_MudOrKick%Array(j)
        ELSE
            FlowEl(i)%MaterialType = 1
        END IF
        
            
        FlowEl(i)%dPdLgrav = 0.0
        FlowEl(i)%StaticPressDiff = 0.0
        
        !WRITE (*,*) 'Ch density, length, Type' , i, FlowEl(i)%density , FlowEl(i)%length, FlowEl(i)%MaterialType

        j = j + 1
            
    END DO
            
                !WRITE (*,*) 'inlet and outlet mud theta600' , FlowEl(1)%Theta600 , FlowEl(NumbEl)%Theta600
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    !!!!!!!!!!!!!!! Open Hole
    OpenholeFirstEl = NoHorizontalEl + NoStringEl + NoAnnulusEl + NoWelltoChokeEl + 1
    j = 1
    DO i = OpenholeFirstEl , NumbEl ! = NoHorizontalEl + NoStringEl + NoCasingEl + NoOpenHoleEl
        FlowEl(i)%StartX = MudSystem%Xstart_OpMudElement%Array(j)
        FlowEl(i)%EndX = MudSystem%Xend_OpMudElement%Array(j)
        FlowEl(i)%StartTVD = MudSystem%TVDstart_OpMudElement%Array(j)
        FlowEl(i)%EndTVD = MudSystem%TVDend_OpMudElement%Array(j)
        FlowEl(i)%Id = 0.0
        FlowEl(i)%Od = MudSystem%PipeOD_OpMudElement%Array(j)
        FlowEl(i)%density = MudSystem%Density_OpMudElement%Array(j)
        FlowEl(i)%MaterialType = MudSystem%MudTypeOp_MudElement%Array(j)     ! = 0 for mud, = 1 for gas kick, = 4 for air
        
        FlowEl(i)%Dhyd = FlowEl(i)%Od
        FlowEl(i)%Length = ABS(FlowEl(i)%EndX - FlowEl(i)%StartX)
        FlowEl(i)%DepthDiff = ABS(FlowEl(i)%StartTVD - FlowEl(i)%EndTVD)
        FlowEl(i)%Area = PI / 4.0 * FlowEl(i)%Od**2 / Convfttoinch**2        
        FlowEl(i)%Volume = FlowEl(i)%Area * REAL(FlowEl(i)%Length)
        IF (MudSystem%MudTypeOp_MudElement%Array(j) == 0) THEN           ! = 0 for mud, = 1 for gas kick, = 4 for air
            FlowEl(i)%MaterialType = 0
        ELSE IF (MudSystem%MudTypeOp_MudElement%Array(j) == 4) THEN
            FlowEl(i)%MaterialType = 4
        ELSE
            FlowEl(i)%MaterialType = 1
        END IF
        
        j = j + 1
        IF (FlowEl(i)%MaterialType == 1) THEN
            FlowEl(i)%dPdLgrav = 0.0
        ELSE
            FlowEl(i)%dPdLgrav = 0.052 * FlowEl(i)%density
        END IF
        FlowEl(i)%StaticPressDiff = FlowEl(i)%dPdLgrav * REAL(FlowEl(i)%DepthDiff)
    
        !WRITE (*,*) 'O density, length, Type of ith element' , i, FlowEl(i)%density , FlowEl(i)%length, FlowEl(i)%MaterialType

    END DO
            !WRITE (*,*) 'inlet and outlet mud theta600' , FlowEl(1)%Theta600 , FlowEl(NumbEl)%Theta600
    
    !DO i = NoHorizontalEl + 1, NoHorizontalEl + NoStringEl
    !    WRITE (*,*) 'element No, Start , Length, DeltaPSt', i, FlowEl(i)%StartX, FlowEl(i)%Length, FlowEl(i)%StaticPressDiff
    !END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
!!!!!!!!!!!!!!!!!!  Mud Element Viscosity Calculation
    
    DO i = 1 , NumbEl
        IF (FlowEl(i)%MaterialType /= 0)    CYCLE
        IF (FlowEl(i)%density < 8.33) THEN
            WRITE (*,*) 'H, S, A, Ch, O', NoHorizontalEl , StringFirstEl , AnnulusFirstEl , NoWellToChokeEl , OpenHoleFirstEl
            
            WRITE (*,*) 'Element No, Volume(gal), Density(ppg), Type'
            DO j = 1 , NumbEl
                WRITE (*,*) j, FlowEl(j)%Volume * ConvFt3toUSGal , FlowEl(j)%Density, FlowEl(j)%MaterialType
            END DO
        
            CALL ErrorStop('Density below water density (8.33)', i)
            CYCLE
        END IF
        
        FlowEl(i)%Theta600 = Theta600Refrence + 11.0 * (MAX(FlowEl(i)%density , 8.33) - DensityRefrence)
        FlowEl(i)%Theta300 = Theta300Refrence + 6.0 * (MAX(FlowEl(i)%density , 8.33) - DensityRefrence)
        !WRITE (*,*) 'i, Theta600, Theta600', i, FlowEl(i)%Theta600, FlowEl(i)%Theta300
    END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    
    FlowEl(:)%FlowRate = 0.0
    FlowEl(:)%StartPress = 0.0
    FlowEl(:)%EndPress = 0.0
    
    !WRITE (*,*) 'NoGasPocket=' , SIZE(GasPocketWeight%Array)
    
    
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    MudSystem%OpenArea32 = ((ChokeManifold%LeftManualChoke * (1.0 - ManChoke1Plug * REAL(ChokeProblems%ManualChoke1PluggedPercent)) / 100.0) + 0.5 * ManChoke1Washout) * Choke%ChokeAreaFullyOpen
    MudSystem%OpenArea33 = CHOOKE(1)%AreaChokeFinal * Convfttoinch**2              !  (1.d0 - CHOOKE(1)%PercentClose) * ChokeAreaFullyOpen
    MudSystem%OpenArea34 = CHOOKE(2)%AreaChokeFinal * Convfttoinch**2              !  (1.d0 - CHOOKE(2)%PercentClose) * ChokeAreaFullyOpen
    MudSystem%OpenArea35 = ((ChokeManifold%RightManualChoke * (1.0 - ManChoke2Plug * REAL(ChokeProblems%ManualChoke2PluggedPercent)) / 100.0) + 0.5 * ManChoke2Washout) * Choke%ChokeAreaFullyOpen
    !write(*,*) 'OpenArea32=' , OpenArea32, active32, ManChoke1Plug, ManualChoke1PluggedPercent, ManChoke1Washout
    !write(*,*) 'OpenArea33=' , OpenArea33, active33
    !write(*,*) 'OpenArea34=' , OpenArea34, active34, !HydraulicChoke2PluggedPercent
    !write(*,*) 'OpenArea35=' , OpenArea35, active35, ManChoke2Plug, ManualChoke2PluggedPercent, ManChoke2Washout
    InstantaneousTotalOpenChokeArea = MudSystem%OpenArea32 * active32 + MudSystem%OpenArea33 * active33 + MudSystem%OpenArea34 * active34 + MudSystem%OpenArea35 * active35 + ChokeBypassArea * active29
    OldTotalOpenChokeArea = TotalOpenChokeArea
    !WRITE (*,*) 'Instantaneous / Old TotalOpenChokeArea', InstantaneousTotalOpenChokeArea, OldTotalOpenChokeArea
    IF (OldTotalOpenChokeArea <= 0.01 * Choke%ChokeAreaFullyOpen)     OldTotalOpenChokeArea = 0.01 * Choke%ChokeAreaFullyOpen
    IF (InstantaneousTotalOpenChokeArea <= 0.01 * Choke%ChokeAreaFullyOpen) THEN
        MudSystem%WellToChokeManifoldOpen = .FALSE.
        OldTotalOpenChokeArea = 0.01 * Choke%ChokeAreaFullyOpen
        !WRITE (*,*) ' Choke Controler Here 2'
        
    ELSE IF (InstantaneousTotalOpenChokeArea > 0.5 * Choke%ChokeAreaFullyOpen .OR. MudSystem%WellToPitsOpen) THEN
                    ! mud flows through well to bell nipple, or choke is rather open
        TotalOpenChokeArea = InstantaneousTotalOpenChokeArea
        !WRITE (*,*) ' Choke Controler Here 1'
    ELSE
        IF (InstantaneousTotalOpenChokeArea > OldTotalOpenChokeArea) THEN
            TotalOpenChokeArea = MIN(1.1 * OldTotalOpenChokeArea , InstantaneousTotalOpenChokeArea)
            WRITE (*,*) ' Choke is opening' !, TotalOpenChokeArea, OldTotalOpenChokeArea, InstantaneousTotalOpenChokeArea
        ELSE IF (InstantaneousTotalOpenChokeArea < OldTotalOpenChokeArea) THEN
            TotalOpenChokeArea = MAX(0.9 * OldTotalOpenChokeArea , InstantaneousTotalOpenChokeArea)
            WRITE (*,*) ' Choke is closing' !, TotalOpenChokeArea, OldTotalOpenChokeArea, InstantaneousTotalOpenChokeArea
            !ExitMass = 0.1 * ExitMass       ! for decreasing mass exit due to expansion
        ELSE
            TotalOpenChokeArea = InstantaneousTotalOpenChokeArea
        END IF
    END IF
    !WRITE (*,*) 'TotalOpenChokeArea=', TotalOpenChokeArea 
    
        
    
    
            
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!            
    
!!!!!!!!!!!!!!!!!!!!  Well Head Condition ( Open or Closed )
    ChokeKroneckerDelta = 0
    !WRITE (*,*) ' WelltoPitsOpen = ', WelltoPitsOpen
    !WRITE (*,*) ' WellToChokeManifoldOpen = ', WellToChokeManifoldOpen
    IF (MudSystem%WellToPitsOpen == .FALSE. .AND. MudSystem%WellToChokeManifoldOpen == .FALSE.) THEN
        !WRITE (*,*) ' WelltoPitsOpen = ', WelltoPitsOpen
        !WRITE (*,*) ' WellToChokeManifoldOpen = ', WellToChokeManifoldOpen
        
        WellHeadOpen = .FALSE.
        !WRITE (*,*) ' Well Head is closed '
    ELSE
        WellHeadOpen = .TRUE.
    END IF
    
    IF (MudSystem%WellToPitsOpen == .FALSE. .AND. MudSystem%WellToChokeManifoldOpen)        ChokeKroneckerDelta = 1
    !WRITE (*,*) ' WelltoPitsOpen=' , WelltoPitsOpen, 'WellToChokeManifoldOpen= ' , WellToChokeManifoldOpen
    !WRITE (*,*) 'WellToChokeOpen , WellToChokeWasOpen', WellToChokeManifoldOpen, WellToChokeManifoldWasOpen
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    IF (NoGasPocket > 0) THEN
        !IF (NoWelltoChokeEl > 0 .AND. (FlowEl(OpenholeFirstEl - 1)%MaterialType == 1 .OR. KickWasExitingThroughChoke) .AND. WellHeadOpen) THEN
        IF (MudSystem%Choke_Kick_Saved_Volume_Final > 0.d0) THEN
            !WRITE (*,*) 'Choke_Kick_Saved_Volume_Final=', Choke_Kick_Saved_Volume_Final, (GasPocketModifiedVol%Array(NoGasPocket) * Convft3ToUSgal)
            ExitMass = REAL(MudSystem%Choke_Kick_Saved_Volume_Final / (GasPocketModifiedVol%Array(NoGasPocket) * Convft3ToUSgal)) * GasPocketWeight%Array(NoGasPocket)
            !IF (NoGasPocket > 1)        ExitMass = ExitMass + (SUM(GasPocketDeltaVol%Array(1 : NoGasPocket - 1)) / GasPocketNewVol%Array(NoGasPocket)) * GasPocketWeight%Array(NoGasPocket) 
            !WRITE (*,*) 'Total Exit mass=' , ExitMass !, DeltaVolumePipe, MudSystem%StringFlowRate * dt / ConvMinToSec
            IF (FlowEl(OpenholeFirstEl - 1)%MaterialType /= 1) THEN
                CALL RemoveGasPocket(NoGasPocket)
                WRITE (*,*) 'Choke_Kick_Saved_Volume_Final=', MudSystem%Choke_Kick_Saved_Volume_Final
            END IF
            
        END IF
        
        
        !NoGasPocket = SIZE(GasPocketWeight%Array)
        IF (NoGasPocket > 0)         CALL GasPocketFlowElementTransformer
    END IF
    
    IF (NoWelltoChokeEl > 0) THEN
        OldChokeDensity = ChokeDensity
        !ChokeDensity = SUM(FlowEl(AnnulusLastEl + 1 : OpenholeFirstEl - 1)%Density * FlowEl(AnnulusLastEl + 1 : OpenholeFirstEl - 1)%Volume) / SUM(FlowEl(AnnulusLastEl + 1 : OpenholeFirstEl - 1)%Volume)
        !ChokeDensity = MAX(ChokeDensity , ChokeMinDensity)
        VolumeDensityProduct = 0.0
        Do i = AnnulusLastEl + 1 , OpenholeFirstEl - 1
            VolumeDensityProduct = VolumeDensityProduct + ((MAX(FlowEl(i)%Density , ChokeMinDensity)) * FlowEl(i)%Volume)
        End Do
        ChokeDensity = VolumeDensityProduct / SUM(FlowEl(AnnulusLastEl + 1 : OpenholeFirstEl - 1)%Volume)
        
        !IF ((ChokeDensity - OldChokeDensity) > (MaxChokeDensityChange * dt / Convmintosec)) THEN
        !    ChokeDensity = OldChokeDensity + (MaxChokeDensityChange * dt / Convmintosec)
        !ELSE IF ((OldChokeDensity - ChokeDensity) > (MaxChokeDensityChange * dt / Convmintosec)) THEN
        !    ChokeDensity = OldChokeDensity - (MaxChokeDensityChange * dt / Convmintosec)
        !END IF
                    
        !WRITE (*,*) 'ChokeDensity=', ChokeDensity
        Kchoke = (ChokeDensity / ((2.0 * 89158.0) * (0.26 * 0.61 * TotalOpenChokeArea)**2)) * 4.d0    ! *4.d0: seyyed gofte
    END IF
    
    IF (BitTotallyPluged == .FALSE.) THEN
        KBit = FlowEl(StringLastEl)%density / 12042.0 / BitCd**2 / BitTotNozzArea**2
    END IF
    
    
    
END SUBROUTINE WellPressureDataTransfer
    
    
SUBROUTINE GasPocketFlowElementTransformer

    !! This Subroutine makes relationship between fluid flow elements and gas pocket elements
    !! The main variable of this subroutine is GasPocketFlowEl
    !! For each gas pocket, GasPocketFlowEl has a row, which columns are numbers of flow elements which are contains that gas pocket
    !! For example, if we have two gas pockets and they are in (7,8) and (12,3,4) flow elements respectively,
    !! GasPocketFlowEl(1,:) = [12 , 3 , 4] and GasPocketFlowEl(2,:) = [7 , 8 , 0]


    USE FricPressDropVars
    USE MudSystemVARIABLES
    USE GeoElements_FluidModule
    USE Fluid_Flow_Startup_Vars
    USE KickVariables
    USE CMudPropertiesVariables
    USE CBopStackVariables
    USE CDownHoleVariables
    USE CError
    
    IMPLICIT NONE
    
    INTEGER :: i , j , k
    REAL :: PressureCorrection
    PressureCorrection = 1.0
    
    ALLOCATE(GasPocketFlowEl(NoGasPocket , 1))
    GasPocketFlowEl(:,:) = 0
    !WRITE (*,*) 'NoGasPocket=' , NoGasPocket
    
    i = 1
    j = 1
    DO k = OpenholeFirstEl , NumbEl
        IF (FlowEl(k)%MaterialType == 1) THEN
            IF (j > 1) THEN
                IF (GasPocketFlowEl(i , j - 1) /= k - 1) THEN
                    i = i + 1
                    j = 1
                END IF
            END IF
            
            IF (SIZE(GasPocketFlowEl , dim = 2) < j) THEN
                ALLOCATE(tempGasPocketFlowEl(SIZE(GasPocketFlowEl , dim = 1) , SIZE(GasPocketFlowEl , dim = 2)))
                tempGasPocketFlowEl = GasPocketFlowEl
                DEALLOCATE(GasPocketFlowEl)
                ALLOCATE(GasPocketFlowEl(SIZE(tempGasPocketFlowEl , dim = 1) , SIZE(tempGasPocketFlowEl , dim = 2) + 1))
                GasPocketFlowEl(: , 1 : SIZE(tempGasPocketFlowEl , dim = 2)) = tempGasPocketFlowEl
                GasPocketFlowEl(: , SIZE(GasPocketFlowEl , dim = 2)) = 0
                DEALLOCATE(tempGasPocketFlowEl)
            END IF
            
            GasPocketFlowEl(i , j) = k
            j = j + 1
        END IF
    END DO
    !WRITE (*,*) 'GasPocketFlowEl=' , GasPocketFlowEl
    
    DO k = AnnulusFirstEl , OpenholeFirstEl - 1
        IF (FlowEl(k)%MaterialType == 1) THEN
            IF (j > 1) THEN
                IF (k == AnnulusFirstEl  .AND.  GasPocketFlowEl(i , j - 1) == NumbEl) THEN
                    !WRITE (*,*) 'Kick is Around Bit'
                ELSE IF (GasPocketFlowEl(i , j - 1) /= k - 1) THEN
                    i = i + 1
                    j = 1
                    !WRITE (*,*) 'i, j, k', i, j, k, FlowEl(k)%MaterialType
                END IF
            END IF
            
            IF (j > SIZE(GasPocketFlowEl , dim = 2)) THEN
                ALLOCATE(tempGasPocketFlowEl(SIZE(GasPocketFlowEl , dim = 1) , SIZE(GasPocketFlowEl , dim = 2)))
                tempGasPocketFlowEl = GasPocketFlowEl
                DEALLOCATE(GasPocketFlowEl)
                ALLOCATE(GasPocketFlowEl(SIZE(tempGasPocketFlowEl , dim = 1) , SIZE(tempGasPocketFlowEl , dim = 2) + 1))
                GasPocketFlowEl(: , 1 : SIZE(tempGasPocketFlowEl , dim = 2)) = tempGasPocketFlowEl
                GasPocketFlowEl(: , SIZE(GasPocketFlowEl , dim = 2)) = 0
                DEALLOCATE(tempGasPocketFlowEl)
            END IF
            IF (i > NoGasPocket) THEN
                WRITE (*,*) 'GasPocketFlowEl', GasPocketFlowEl
                WRITE (*,*) 'i , j , k', i, j, k
                WRITE (*,*) 'H, S, A, Ch, O', NoHorizontalEl , StringFirstEl , AnnulusFirstEl , NoWellToChokeEl , OpenHoleFirstEl
            
                WRITE (*,*) 'Element No, Volume(gal), Density(ppg),         Type'
                DO j = 1 , NumbEl
                    WRITE (*,*) j, FlowEl(j)%Volume * ConvFt3toUSGal , FlowEl(j)%Density, FlowEl(j)%MaterialType
                END DO
                !WRITE (*,*) 'Op Last: Volume(gal), Type', NumbEl, FlowEl(NumbEl)%Volume * ConvFt3ToUSgal , FlowEl(NumbEl)%MaterialType
                !WRITE (*,*) 'Ann First: Volume(gal), Type', AnnulusFirstEl, FlowEl(AnnulusFirstEl)%Volume * ConvFt3ToUSgal , FlowEl(AnnulusFirstEl)%MaterialType
                !IF (ChokeKroneckerDelta == 1) THEN
                !    WRITE (*,*) 'Ann Last: Volume(gal), Type', AnnulusLastEl, FlowEl(AnnulusLastEl)%Volume * ConvFt3ToUSgal , FlowEl(AnnulusLastEl)%MaterialType
                !    WRITE (*,*) 'Ch First: Volume(gal), Type', AnnulusLastEl + 1 , FlowEl(AnnulusLastEl + 1)%Volume * ConvFt3ToUSgal , FlowEl(AnnulusLastEl + 1)%MaterialType
                !END IF
                
                CALL ErrorStop('Error in calculating GasPocketFlowEl')
            END IF

            GasPocketFlowEl(i , j) = k
            !WRITE (*,*) 'GasPocketFlowEl=' , GasPocketFlowEl
            j = j + 1
            
        END IF
        !IF (i > NoGasPocket)        CALL ErrorStop('Error in GasPocketFlowEl' , i)
    END DO
    
    !i = SIZE(GasPocketFlowEl , dim = 1)
    !j = SIZE(GasPocketFlowEl , dim = 2)
    !WRITE(*,*) 'GasPocketFlowEl Data' , GasPocketFlowEl
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
    GasPocketModifiedVol%Array(:) = 0.d0
    
    DO i = 1 , NoGasPocket
        !WRITE(*,*) 'GasPocketFlowEl Data' , i , GasPocketFlowEl(: , i)
        DO j = 1 , SIZE(GasPocketFlowEl , dim = 2)
            IF (GasPocketFlowEl(i , j) > 0)     GasPocketModifiedVol%Array(i) = GasPocketModifiedVol%Array(i) + FlowEl(GasPocketFlowEl(i , j))%Volume
        END DO
    END DO
    i = AnnulusLastEl
    j = OpenholeFirstEl - 1
    !WRITE (*,*) ' b) GasPocketModifiedVol (ft3)=', GasPocketModifiedVol(1)
    !WRITE (*,*) ' b) GasPocketModifiedVol (gal)=' , GasPocketModifiedVol(1) * 7.48
    !IF (NoWelltoChokeEl > 0 .AND. ANY(FlowEl(i : j)%MaterialType == 1) .AND. FlowEl(j)%MaterialType == 0) THEN
    !        BehindKickMudDensity = FlowEl(GasPocketFlowEl(NoGasPocket , 1) - 1)%Density
    !        FrontKickMudDensity = FlowEl(MAXVAL(GasPocketFlowEl(NoGasPocket , :)) + 1)%Density
    !        KickVolBeforeChoke = REAL(GasPocketModifiedVol%Array(NoGasPocket))
    !        WRITE (*,*) 'BehindKickMudDensity', BehindKickMudDensity
    !        WRITE (*,*) 'FrontKickMudDensity' , FrontKickMudDensity
    !        WRITE (*,*) 'KickVolBeforeChoke', KickVolBeforeChoke
    !END IF
    
    
    GasPocketOldVol%Array(:) = GasPocketModifiedVol%Array(:)
    
    GasPocketOldPress%Array(:) = GasPocketNewPress%Array(:)
    !WRITE (*,*) ' GasPocketweight%Array(1)', GasPocketweight%Array(1)
    IF (GasPocketWeight%Array(1) < 0.0)     CALL ErrorStop(' Error in Calculating Mass ' , KickmdotBCoef)
    
    
    GasPocketDensity%Array(:) = (GasPocketweight%Array(:) / GasPocketModifiedVol%Array(:)) / convft3toUSgal     ! [lbm/ft^3 to ppg]
    DO i = 1 , SIZE(GasPocketFlowEl , dim = 1)
        DO j = 1 , SIZE(GasPocketFlowEl , dim = 2)
            IF (GasPocketFlowEl(i , j) > 0)     FlowEl(GasPocketFlowEl(i , j))%Density = GasPocketDensity%Array(i)
        END DO
    END DO
    
    
    !WRITE (*,*) 'Kick density (ppg)=' , GasPocketDensity(1)
    DownHole%InfluxRate = MAX(((KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1))) / GasPocketDensity%Array(1) * ConvMinToSec) , 0.0)
    !WRITE (*,*) ' InfluxRate (gpm) =', InfluxRate
    
    i = OpenholeFirstEl - 1
    !WRITE (*,*) 'MaterialType, NoWelltoChokeEl', FlowEl(i)%MaterialType , NoWelltoChokeEl
    !IF (NoWelltoChokeEl > 0 .AND. FlowEl(OpenholeFirstEl - 1)%MaterialType == 1 .AND. WellHeadOpen) THEN     ! kick is last element in choke line and does not exit.
    IF (MudSystem%Choke_Kick_Saved_Volume_Final > 0.d0 .AND. WellHeadOpen) THEN
        WRITE (*,*) 'Kick is exiting through the choke'
        !ExitMass = GasPocketDensity%Array(NoGasPocket) * (SUM(GasPocketFlowInduced%Array(:)) + (DeltaVolumePipe * ConvMinToSec / dt) + MudSystem%StringFlowRate) / ConvMinToSec * dt
        GasPocketWeight%Array(NoGasPocket) = GasPocketweight%Array(NoGasPocket) - ExitMass
        IF (GasPocketWeight%Array(NoGasPocket) > 0.0) THEN
            GasPocketOldPress%Array(NoGasPocket) = GasPocketWeight%Array(NoGasPocket) * GasType(KickGasType)%GasConstant * &
                                                    GasPocketCompressibility%Array(NoGasPocket) * GasPocketNewTemp%Array(NoGasPocket) / GasPocketModifiedVol%Array(NoGasPocket)
            GasPocketDensity%Array(NoGasPocket) = (GasPocketweight%Array(NoGasPocket) / GasPocketModifiedVol%Array(NoGasPocket)) / convft3toUSgal
            GasPocketNewVol%Array(:) = GasPocketOldVol%Array(:)
            !FlowEl(i)%Density = GasPocketDensity%Array(NoGasPocket)    ! for choke pressure drop calculation
            !WRITE (*,*) 'kick weight, pressure & density', GasPocketWeight%Array(NoGasPocket), GasPocketOldPress%Array(NoGasPocket), FlowEl(i)%Density
        ELSE    ! gas pocket is escaped from the well completely
            WRITE (*,*) ' Last Pocket Removed'
            WRITE (*,*) ' GasPocketFlowEl', GasPocketFlowEl
            CALL RemoveGasPocket(NoGasPocket)
            
            !WRITE (*,*) ' GasPocketFlowEl', GasPocketFlowEl
            
        END IF
        
    ELSE IF (NoGasPocket == 1 .OR. WellHeadOpen) THEN        ! kick is not last element of choke line or well head is closed
        !WRITE (*,*) ' Gas Kick Pressure (psi)  (before) = ' , GasPocketNewPress(1)
        !DO j = 1 , NoGasPocket
            GasPocketNewPress%Array(:) = GasPocketNewPress%Array(:) * REAL(GasPocketNewVol%Array(:) / GasPocketOldVol%Array(:))
            GasPocketOldPress%Array(:) = GasPocketNewPress%Array(:)
            GasPocketNewVol%Array(:) = GasPocketOldVol%Array(:)
        !WRITE (*,*) ' Gas Kick Pressure (psi) (after) = ' , GasPocketOldPress(1), REAL(GasPocketNewVol(1) / GasPocketOldVol(1))
        !END DO
    
    ELSE    ! Many gas pocket is in the well and well head is closed
        
103     FORMAT (2I, 4X, (F9.3), 5X, (F9.3))
        
        !WRITE (*,*) 'Before 0'
        !DO i = 1 , NoGasPocket
        !    WRITE (*,*)  i, REAL(GasPocketNewPress%Array(i)), REAL(GasPocketNewVol%Array(i) * Convft3toUSgal), REAL(GasPocketOldVol%Array(i) * Convft3toUSgal)
        !END DO
        GasPocketNewPress%Array(NoGasPocket) = GasPocketOldPress%Array(NoGasPocket)
        GasPocketNewVol%Array(NoGasPocket) = GasPocketCompressibility%Array(NoGasPocket) * GasType(KickGasType)%GasConstant * &
                            GasPocketNewTemp%Array(NoGasPocket) * GasPocketWeight%Array(NoGasPocket) / GasPocketNewPress%Array(NoGasPocket)
        
        
        !GasPocketNewPress%Array(:) = GasPocketNewPress%Array(:) * REAL(GasPocketNewVol%Array(:) / GasPocketOldVol%Array(:))
        !GasPocketOldPress%Array(:) = GasPocketNewPress%Array(:)
        
        !WRITE (*,*) 'Before 1'
        !DO i = 1 , NoGasPocket
        !    WRITE (*,*)  i, REAL(GasPocketNewPress%Array(i)), REAL((GasPocketNewVol%Array(i) * Convft3toUSgal))
        !END DO
        
                !WRITE (*,*) 'Mid'
        DO WHILE (ABS(PressureCorrection * 10.0) > KickConvergenceTolerance)
            DO j = NoGasPocket - 1 , 1 , -1
                CALL KickFunctionsCalculator(GasPocketNewPress%Array(j) , j , 5)
                GasPocketNewVol%Array(j) = GasPocketCompressibility%Array(j) * GasType(KickGasType)%GasConstant * &
                            GasPocketNewTemp%Array(j) * GasPocketWeight%Array(j) / GasPocketNewPress%Array(j)
                !WRITE (*,*) j, REAL(GasPocketNewPress%Array(j)), REAL(GasPocketNewVol%Array(j) * Convft3toUSgal)
            END DO
            PressureCorrection = (SUM(GasPocketNewVol%Array(:)) - SUM(GasPocketOldVol%Array(:))) / SUM(GasPocketNewVol%Array(:) / GasPocketNewPress%Array(:))
            GasPocketNewPress%Array(NoGasPocket) = GasPocketNewPress%Array(NoGasPocket) + PressureCorrection
            GasPocketNewVol%Array(NoGasPocket) = GasPocketCompressibility%Array(NoGasPocket) * GasType(KickGasType)%GasConstant * &
                            GasPocketNewTemp%Array(NoGasPocket) * GasPocketWeight%Array(NoGasPocket) / GasPocketNewPress%Array(NoGasPocket)
            j = NoGasPocket
            !WRITE (*,*) j, REAL(GasPocketNewPress%Array(j)), REAL(GasPocketNewVol%Array(j) * Convft3toUSgal)
        END DO
        
        !WRITE (*,*) 'Pressure Correction = ', PressureCorrection
        WRITE (*,*) 'After'
        DO i = 1 , NoGasPocket
            WRITE (*,*) i, REAL(GasPocketNewPress%Array(i)), REAL((GasPocketNewVol%Array(i) * Convft3toUSgal))
        END DO
        
        GasPocketOldPress%Array(:) = GasPocketNewPress%Array(:)
        
    END IF
    
    
    IF (GasPocketFlowEl(NoGasPocket , 1) == 0) THEN
        CALL RemoveGasPocket(1)
        
        !WRITE (*,*) 'GasPocketFlowEl', GasPocketFlowEl
        WRITE (*,*) 'First Pocket Removed'
        !WRITE (*,*) 'GasPocketFlowEl', GasPocketFlowEl
    END IF
    
    
    
    DownHole%KickVolume = INT(SUM(GasPocketOldVol%Array(:)) * convft3toUSgal / 42. * 10.0) / 10.0
    
    
    !WRITE (*,*) ' Gas Kick Pressure (psi) = ' , GasPocketOldPress(1) , INT((GasPocketNewVol(1) / GasPocketOldVol(1)) * 1000.d0) / 1000.d0
    !WRITE (*,*) ' b) GasPocketNewVol (ft^3) = ' , GasPocketOldVol(1)
        
        
    
    END SUBROUTINE
    
    
SUBROUTINE RemoveGasPocket(ilocal)
    
    USE KickVARIABLES
    USE MudSystemVARIABLES
    
    IMPLICIT NONE
    
    INTEGER :: ilocal

    DEALLOCATE(KickJacobian , OldKickJacobian , KickVandPFunction , KickUnknownVector , KickCorrectionVector)
    
        
    WRITE (*,*) 'Gas Pocket Will be Removed', ilocal
    call GasPocketOldPress%Remove(ilocal)
    call GasPocketNewPress%Remove(ilocal)
    call GasPocketOldTemp%Remove(ilocal)
    call GasPocketNewTemp%Remove(ilocal)
    call GasPocketOldVol%Remove(ilocal)
    call GasPocketNewVol%Remove(ilocal)
    call GasPocketDeltaVol%Remove(ilocal)
    call GasPocketFlowInduced%Remove(ilocal)
    call GasPocketModifiedVol%Remove(ilocal)
    call GasPocketWeight%Remove(ilocal)
    call GasPocketDensity%Remove(ilocal)
    call GasPocketCompressibility%Remove(ilocal)
    
    NoGasPocket = NoGasPocket - 1
     
    IF (NoGasPocket > 0) THEN
        ALLOCATE(KickJacobian(2 * NoGasPocket , 2 * NoGasPocket) , OldKickJacobian(2 * NoGasPocket , 2 * NoGasPocket) , KickVandPFunction(2 * NoGasPocket))
        ALLOCATE(KickUnknownVector(2 * NoGasPocket) , KickCorrectionVector(2 * NoGasPocket))
        
        IF (ALLOCATED(GasPocketFlowEl)) THEN
            ALLOCATE(tempGasPocketFlowEl(NoGasPocket , SIZE(GasPocketFlowEl , dim = 2)))
            tempGasPocketFlowEl = GasPocketFlowEl(1 : NoGasPocket , :)
            DEALLOCATE(GasPocketFlowEl)
            ALLOCATE(GasPocketFlowEl(SIZE(tempGasPocketFlowEl , dim = 1) , SIZE(tempGasPocketFlowEl , dim = 2)))
            GasPocketFlowEl = tempGasPocketFlowEl
            DEALLOCATE(tempGasPocketFlowEl)
        END IF
        
    ELSE        ! NoGasPocket = 0
        
        
        CALL GasPocketOldPress%Empty
        CALL GasPocketNewPress%Empty
        CALL GasPocketOldTemp%Empty
        CALL GasPocketNewTemp%Empty
        CALL GasPocketOldVol%Empty
        CALL GasPocketNewVol%Empty
        CALL GasPocketDeltaVol%Empty
        CALL GasPocketFlowInduced%Empty
        CALL GasPocketModifiedVol%Empty
        CALL GasPocketWeight%Empty
        CALL GasPocketDensity%Empty
        
        DownHole%InfluxRate = 0.0
        DownHole%KickVolume = 0.0
                
    END IF         
    
    !WRITE (*,*) 'Gas Pocket Removed', NoGasPocket + 1, ALLOCATED(GasPocketWeight%Array)
    
END SUBROUTINE