SimulationCore2/FluidFlow/kick/Gas_Kick_Calculator.f90

SUBROUTINE GasKickCalculator
    
    USE KickVariables
    Use TD_WellGeometry
    Use CReservoirVariables
    Use CFormationVariables
    USE Fluid_Flow_Startup_Vars
    USE PressureDisplayVARIABLES
    USE FricPressDropVars
    USE MudSystemVARIABLES
    USE CMudPropertiesVariables
    USE CError
    USE , INTRINSIC :: IEEE_ARITHMETIC
    
    
    !! Note: a subject that may be confusing is that when we use gauge pressure and when're using absolute pressure?!
    !! all stated pressure are gauge pressure, so I do like this.
    !! only when we want to use a state equation (like ideal gas equation), we should use absolute equation and so I do this.
    !! Thus gas pocket pressure are all absolute pressure.
    
    IMPLICIT NONE
    
    INTEGER :: i , j , k , l
    
    SolvingEquationError = .FALSE.
    KickVandPFunction(:) = 0.d0
    KickJacobian(: , :) = 0.d0
    
!====================================================
!            Gas Properties (Methane Gas)
!====================================================
    !GasPocketNewTemp%Array(1) = 600
    BottomHoleTemperature = 600
    KickFluxAvgPressure = (BottomHolePress + FormPressure) / 2 + StandardPress
    KickFluxAvgTemperature = (FormTemperature + BottomHoleTemperature) / 2
    KickFluxAvgCompressibility = 0.98d0
        
    K_Aa = (5.8742362 * 10.**(-3) * KickFluxAvgTemperature**1.2288) / (511.1728532 + KickFluxAvgTemperature)
    K_Bb = 5.5565586 + (1000.01 / KickFluxAvgTemperature)
    K_Cc = 2.47862 - 0.12294 * K_Bb
    GasKickSIDensity = KickFluxAvgPressure / (KickFluxAvgCompressibility * &
                                KickFluxAvgTemperature * GasType(KickGasType)%GasConstant) * Convpcftogpcm3
    
    GasKickViscosity = K_Aa * EXP(K_Bb * GasKickSIDensity**K_Cc)
    !!!!!!!!!!!!!!!!!!!!!!!!!
    
    !!!!!!!!!!!!!! Calculating compressibility for bottom hole condition
    !K_BHTpr = BottomHoleTemperature / KickTc
    !K_BHPpr = (BottomHolePress + StandardPress) / KickPc
    
    !K_A_Bottomhole = 3.53 * K_BHPpr
    !K_B_Bottomhole = 10.0**(0.9813 * K_BHTpr)
    !K_C_Bottomhole = 0.274 * (K_BHPpr**2)
    !K_D_Bottomhole = 10.0**(0.8157 * K_BHTpr)
    
    BottomHoleCompressibility = 0.98d0 !1. - (K_A_Bottomhole / K_B_Bottomhole) + (K_C_Bottomhole / K_D_Bottomhole)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    GasKickBg = 0.00504 * KickFluxAvgCompressibility * KickFluxAvgTemperature / KickFluxAvgPressure  ![bbl/SCF]


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !WRITE (*,*) 'Gas Kick Top'
!===> Kick Flow Rate Calculation
    IF (FormPressure > BottomHolePress) THEN
        KickmdotACoef = 10.0**(-8) * 1.15741d0 * 7.080 * FormPermeability * REAL(KickFormLength) * GasType(KickGasType)%StDensity / &
                            (GasKickViscosity * GasKickBg * LOG(10000.0))
        !IF (WellHeadOpen .AND. NoGasPocket == 1)       KickmdotACoef = (1.d0 + 2.d0) * KickmdotACoef
    ELSE
        KickmdotACoef = 0.0
    END IF
    
    i = StringLastEl
    j = OpenholeFirstEl - 1
    k = GasPocketFlowEl(1 , 1)
    KickmdotBCoef = FormPressure + StandardPress    !! - Sum(static and friction pressure loss) of flow elements below gas pocket, see below
    IF (FormPressure > BottomHolePress) THEN
    
        !WRITE (*,*) 'k , i, j' , k , i, j
        IF (k >= OpenholeFirstEl) THEN      ! Bottom of active kick is in openhole
            KickmdotBCoef = KickmdotBCoef - (SUM(FlowEl(OpenholeFirstEl : k)%StaticPressDiff))      !+ SUM(FlowEl(j + 1 : GasPocketFlowEl(1 , 1) - 1)%FricPressLoss 
            !WRITE (*,*) '1 SUM(FlowEl(j + 1 : k)%FricPressLoss', k, SUM(FlowEl(j + 1 : k)%FricPressLoss)
        ELSE IF (k < OpenholeFirstEl) THEN       ! bottom of 1st gas pocket (active kick) is in annulus ond/or choke line only
            KickmdotBCoef = KickmdotBCoef - SUM(FlowEl(OpenholeFirstEl : NumbEl)%StaticPressDiff) &
                                - (SUM(FlowEl(AnnulusFirstEl : k)%StaticPressDiff) + SUM(FlowEl(AnnulusFirstEl : k)%FricPressLoss))
            !WRITE (*,*) '2 SUM(FlowEl(j + 1 : NumbEl)%FricPressLoss', k, SUM(FlowEl(j + 1 : NumbEl)%FricPressLoss)
        END IF
        !WRITE (*,*) ' KickmdotBCoef=', KickmdotBCoef
    END IF
    !WRITE (*,*) 'Kick A, B', KickmdotACoef, KickmdotBCoef

    DO l = 1 , NoGasPocket
        KickUnknownVector(2 * l - 1) = GasPocketNewVol%Array(l)
        KickUnknownVector(2 * l) = GasPocketNewPress%Array(l)
    END DO
    
    IF (WellHeadOpen) THEN
    !!!!!!!!!! Calculation of functions of pocket Pressure and gas Volumes
        !IF (GasPocketElementNo(1) > 0) WRITE (*,*) ' GasPocketElementNo(1) ' , GasPocketElementNo(1)
        !WRITE (*,*) ' Kick Unknown Vector' , KickUnknownVector!(1) , KickUnknownVector(2)
        IF (KickIteration == 1) THEN    ! updating initial guess based on previous time step data
            DO l = 1 , NoGasPocket
                KickUnknownVector(2 * l - 1) = KickUnknownVector(2 * l - 1) + GasPocketDeltaVol%Array(l)
            END DO
        END IF
        
        KickVandPFunction(1) = KickUnknownVector(1) - GasPocketCompressibility%Array(1) * GasType(KickGasType)%GasConstant * &      ! VandP(1) = V(1)
                                GasPocketNewTemp%Array(1) * (GasPocketWeight%Array(1) + KickmdotACoef * MAX(((KickmdotBCoef - KickUnknownVector(2)) * dt) , 0.0)) / KickUnknownVector(2)
        !WRITE (*,*) 'KickVandPFunction(1)',KickVandPFunction(1)
        l = 2 * NoGasPocket
        CALL KickFunctionsCalculator(KickVandPFunction(l) , NoGasPocket , 2)    ! VandP(last) = P(last)
        !WRITE (*,*) 'KickVandPFunction(l)', l, KickVandPFunction(l)
        DO l = 2 , NoGasPocket      ! VandP(Odd) = V(l, l > 1)
            KickVandPFunction(2 * l - 1) = KickUnknownVector(2 * l - 1) - GasPocketCompressibility%Array(l) * GasType(KickGasType)%GasConstant * &
                                GasPocketNewTemp%Array(l) * GasPocketWeight%Array(l) / KickUnknownVector(2 * l)
            !WRITE(*,*) 'KickVandPFunction(V)', l, KickVandPFunction(2 * l - 1)
        END DO
        
        DO l = NoGasPocket - 1 , 1 , -1
            CALL KickFunctionsCalculator(KickVandPFunction(2 * l) , l , 1)    ! VandP(Even) = P(l, l < NoGasPocket)
            !WRITE(*,*) 'KickVandPFunction(P)', l , KickVandPFunction(2 * l)
        END DO
        !!!!!!!!!! END - Calculation of functions of pocket Pressure and gas Volumes
        
        !!!!!!!!!! Calculation of Jacobian
        DO k = 1 , 2 * NoGasPocket      ! Main Diagonal
            KickJacobian(k , k) = 1.d0
        END DO
        
        KickJacobian(1,2) = (GasPocketCompressibility%Array(1) * GasType(KickGasType)%GasConstant * GasPocketNewTemp%Array(1) &
                                * (GasPocketWeight%Array(1) + KickmdotACoef * KickmdotBCoef * dt) / KickUnknownVector(2)**2)      ! Row 1 Finished
        IF (KickJacobian(1,2) == 0.d0) THEN
            CALL Error('KickJacobian(1,2) = 0.0')
            KickJacobian(1,2) = OldKickJacobian(1,2)
        END IF
        
        !WRITE(*,*) 'KickJacobian(1,2)', KickJacobian(1,2)
        
        l = NoGasPocket
        CALL KickFunctionsCalculator(KickJacobian(2 * l , 2 * l - 1) , NoGasPocket , 4)     ! for last Row
        IF (KickJacobian(2 * l , 2 * l - 1) == 0.d0) THEN
            CALL Error ('KickJacobian(Last,Odd) = 0.0')
            KickJacobian(2 * l , 2 * l - 1) = OldKickJacobian(2 * l , 2 * l - 1)
        END IF
        
        DO k = NoGasPocket - 1 , 1 , -1
            KickJacobian(2 * l , 2 * k - 1) = KickJacobian(2 * l , 2 * l - 1)
        END DO      ! Last Row Finished
        !WRITE(*,*) 'KickJacobian(2,1)', KickJacobian(2,1)
        
        
        DO k = 2 , NoGasPocket
            KickJacobian(2 * k - 1 , 2 * k) = GasPocketCompressibility%Array(k) * GasType(KickGasType)%GasConstant * GasPocketNewTemp%Array(k) &
                                * GasPocketWeight%Array(k) / KickUnknownVector(2 * k)**2
            
        END DO      ! Odd Rows (V equations) Finished
        
        DO k = 1 , NoGasPocket - 1
            KickJacobian(2 * k , 2 * k + 2) = -1.d0
        END DO      ! Even Rows (P equations) effect of upper pocket
        
        DO k = 2 , 2 * (NoGasPocket - 1) , 2
            DO l = 1 , k - 1 , 2
                CALL KickFunctionsCalculator(KickJacobian(k , l) , k / 2 , 3)
                IF (KickJacobian(k , l) == 0.d0) THEN
                    WRITE (*,*) 'Jacobian Array = 0.0', k , l
                    CALL Error ('KickJacobian(k , l) = 0.0')
                    KickJacobian(k , l) = OldKickJacobian(k , l)
                END IF
            END DO
        END DO
        
        IF (ANY(IEEE_IS_NaN(KickJacobian)))     CALL ErrorStop ('NaN in calculating Kick Jacobian, Call your Service Provider')
        
    !!!!!!!!!! Solving linear equation in order to finding correction vector for correcting pocket pressure and gas induced flowrates
        KickVandPFunction = -1.d0 * KickVandPFunction
        !WRITE (*,*) 'Max Remainder', MAXVAL(ABS(KickVandPFunction)) , MAXLOC(ABS(KickVandPFunction))
        !WRITE (*,*) 'SIZE(A , dim = 1), SIZE(A , dim = 2), SIZE(b)', SIZE(KickJacobian , dim = 1), SIZE(KickJacobian , dim = 2), SIZE(KickVandPFunction)
        CALL SOLVE_LINEAR_EQUATIONS(KickJacobian , KickCorrectionVector , KickVandPFunction , SolvingEquationError, SIZE(KickCorrectionVector))
        IF (SolvingEquationError)       CALL ErrorStop( ' Error in solving kick equation ' )
        
        KickUnknownVector = KickUnknownVector + KickCorrectionUnderRelaxation * KickCorrectionVector
        
        DO l = 1 , NoGasPocket
            GasPocketNewVol%Array(l) = KickUnknownVector(2 * l - 1)
            IF (IEEE_IS_NaN(GasPocketNewVol%Array(l)))      CALL ErrorStop('Volume of this pocket is Not a Number:', l)
            IF (GasPocketNewVol%Array(l) <= 0.d0)           CALL Error('Volume of this pocket is Negative or Zero:', l)
            GasPocketNewPress%Array(l) = KickUnknownVector(2 * l)
            IF (IEEE_IS_NaN(GasPocketNewPress%Array(l)))        CALL ErrorStop('Pressure of this Pocket is Not a Number:', l)
            IF (GasPocketNewPress%Array(l) <= 0.d0)             CALL ErrorStop('Pressure of this Pocket is Negative or Zero:', l)
        END DO
        
        
        !WRITE(*,*) 'GasPocketDeltaVol (gal)' , GasPocketDeltaVol(1) * 7.48
        !WRITE (*,*) 'GasPocketFlowInduced (gpm), GasPocketNewVol' , GasPocketFlowInduced(1), GasPocketNewVol(1)
        !IF (Kchoke > 0.0)       
        !WRITE(*,*) ' New Vol (ft3), New Press (psi)', GasPocketNewVol(1), GasPocketNewPress(1)
        
    ELSE       ! well haed is closed, so build up process or migration occurs
        
        !WRITE (*,*) 'GasPocketOldPress (before)' , GasPocketOldPress(1)
        GasPocketNewPress%Array(1) = GasPocketOldPress%Array(1) * &
                    (REAL((GasPocketWeight%Array(1) + KickmdotACoef * KickmdotBCoef * dt) / (GasPocketWeight%Array(1) + KickmdotACoef * GasPocketOldPress%Array(1) * dt)))
        !WRITE (*,*) 'GasPocketNewPress (after)' , GasPocketNewPress(1), ((GasPocketWeight(1) + KickmdotACoef * KickmdotBCoef * dt) / (GasPocketWeight(1) + KickmdotACoef * GasPocketOldPress(1) * dt))
                
        
        !WRITE (*,*) ' Well head is closed, GasPocketNewPress =' , GasPocketNewPress(1)
        !WRITE (*,*) 'Old Press, Weight, A, B' , GasPocketOldPress(1), GasPocketWeight(1), KickmdotACoef, KickmdotBCoef
        !WRITE (*,*) 'Numerator and denumerator Gas kick' , KickmdotACoef * KickmdotBCoef * dt , KickmdotACoef * GasPocketOldPress(1) * dt
        
        !WRITE (*,*) ' Gas Kick Volume (ft^3) = ' , GasPocketNewVol(1)
        
    END IF
        
    !DO l = 1 , NoGasPocket
        GasPocketDeltaVol%Array(:) = GasPocketNewVol%Array(:) - GasPocketOldVol%Array(:)
        GasPocketFlowInduced%Array(:) = (GasPocketDeltaVol%Array(:)) / dt * 448.8       ! gpm
    !END DO
    
    GasKickPumpFlowRate = 0.0
    IF (NOT(KickOffBottom) .AND. WellHeadOpen)     GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
    
        !WRITE (*,*) ' No Iteration, KickCorrectionVector =' , KickIteration , KickCorrectionVector(1) , KickCorrectionVector(2)
        !WRITE (*,*) ' Kick Jacobian ', REAL(KickJacobian)
        !WRITE (*,*) ' KickVandPFunction = ' , REAL(KickVandPFunction)
        !WRITE (*,*) ' Kick Unknown Vector = ' , REAL(-KickUnknownVector)

            

    !WRITE (*,*) 'Gas Kick Bottom'


    END SUBROUTINE
    
    
    SUBROUTINE KickFunctionsCalculator(ExitValue , GasPocketNo , CalcMode)
    
    
    USE KickVARIABLES
    USE FricPressDropVars
    USE Fluid_Flow_Startup_Vars
    USE CError
    USE , INTRINSIC :: IEEE_Arithmetic
    
    
    IMPLICIT NONE
    
    REAL(8) :: ExitValue
    INTEGER :: GasPocketNo , CalcMode
    INTEGER :: x
    INTEGER :: y
    INTEGER :: z , i , j
    
    x = GasPocketFlowEl(GasPocketNo , 1)
    IF (GasPocketNo < NoGasPocket)      y = GasPocketFlowEl(GasPocketNo + 1 , 1)
    i = StringLastEl
    j = OpenholeFirstEl - 1
    
    ! Case 1: gas pocket is completely in OP and STARTX of upper gas pocket is also
    ! Case 2: gas pocket is completely in OP and STARTX of upper gas pocket is above Bit
    ! Case 3: gas pocket is AROUNDBIT and so upper gas pocket is in ANN (or Choke line)
    ! Case 4: gas pocket is completely in ANN and upper gas pocket is also
    
    ! CalcMode 1: KickVandPFunction between 2 pocket
    ! CalcMode 2: KickVandPFunction for top gas pocket
    ! CalcMode 3: KickJacobian between 2 Pocket
    ! CalcMode 4: KickJacobian for top (last) gas pocket
    ! CalcMode 5: Calculating pressure of bottom pocket based on upper pocket
    
    
    IF (CalcMode == 1) THEN     ! calculating pressure difference between two pocket, include static pressure difference and frictional 
                                ! pressure difference, use in calculating 'KickVandPFunction'
        ExitValue = KickUnknownVector(2 * GasPocketNo) - KickUnknownVector(2 * GasPocketNo + 2)
        IF (x >= OpenholeFirstEl .AND. y < OpenholeFirstEl) THEN            ! Case 2 , Case 3
            ExitValue = ExitValue - SUM(FlowEl(x : NumbEl)%StaticPressDiff) - SUM(FlowEl(x : NumbEl)%FricPressLoss)   &
                                  - SUM(FlowEl(AnnulusFirstEl : y)%StaticPressDiff) - SUM(FlowEl(AnnulusFirstEl : y)%FricPressLoss)
        ELSE     ! Case 1 , Case 4
            ExitValue = ExitValue - SUM(FlowEl(x : y)%StaticPressDiff) - SUM(FlowEl(x : y)%FricPressLoss)
        END IF
        
        
    ELSE IF (CalcMode == 2) THEN
        ExitValue = KickUnknownVector(2 * GasPocketNo) - StandardPress - Kchoke * FlowEl(OpenholeFirstEl - 1)%Flowrate**2
        IF (x >= OpenholeFirstEl) THEN     ! Gas Pocket is in Openhole
            ExitValue = ExitValue - SUM(FlowEl(x : NumbEl)%StaticPressDiff) - SUM(FlowEl(x : NumbEl)%FricPressLoss)   &
                                    - SUM(FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%StaticPressDiff) - SUM(FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FricPressLoss)
        ELSE                ! Gas Pocket is in Annulus
            ExitValue = ExitValue - SUM(FlowEl(x : OpenholeFirstEl - 1)%StaticPressDiff) - SUM(FlowEl(x : OpenholeFirstEl - 1)%FricPressLoss)
        END IF
        
        
    ELSE IF (CalcMode == 3) THEN    ! calculating derivative of pressure difference between two pocket, relative to change in flowrate
                                    ! use in calculating 'KickJacobian'
        IF (x >= OpenholeFirstEl .AND. y < OpenholeFirstEl) THEN            ! Top kick STARTX is in Annulus
            DO z = x , NumbEl       ! open hole elements
                CALL PartialDerivativeFricToFlowRate(z)
                IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
                    WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
                    WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
                    CALL ErrorStop('NaN in calculating partial derivative')
                END IF
            END DO
            DO z = AnnulusFirstEl , y       ! Annulus elements
                CALL PartialDerivativeFricToFlowRate(z)
                IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
                    WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
                    WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
                    CALL ErrorStop('NaN in calculating partial derivative')
                END IF
            END DO
            ExitValue = ExitValue - (SUM(FlowEl(x : NumbEl)%FricToQPartialDiff) + SUM(FlowEl(AnnulusFirstEl : y)%FricToQPartialDiff)) * 448.8 / dt
        ELSE        ! both pockets are one side of bit
            DO z = x , y
                CALL PartialDerivativeFricToFlowRate(z)
                IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
                    WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
                    WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
                    CALL ErrorStop('NaN in calculating partial derivative')
                END IF
            END DO
            ExitValue = ExitValue - SUM(FlowEl(x : y)%FricToQPartialDiff) * 448.8 / dt
        END IF
        
        
    ELSE IF (CalcMode == 4) THEN    ! partial derivative of frictional pressure drop relative to flowrate for top gas pocket
        ExitValue = - 2.d0 * Kchoke * FlowEl(OpenholeFirstEl - 1)%Flowrate * 448.8 / dt
        IF (x >= OpenholeFirstEl) THEN            ! kick STARTX is in openhole
            DO z = x , NumbEl       ! open hole elements
                CALL PartialDerivativeFricToFlowRate(z)
                IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
                    WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
                    WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
                    CALL ErrorStop('NaN in calculating partial derivative')
                END IF
            END DO
            DO z = AnnulusFirstEl , OpenholeFirstEl - 1       ! Annulus elements
                CALL PartialDerivativeFricToFlowRate(z)
                IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
                    WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
                    WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
                    CALL ErrorStop('NaN in calculating partial derivative')
                END IF
            END DO
            ExitValue = ExitValue - (SUM(FlowEl(x : NumbEl)%FricToQPartialDiff) + SUM(FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FricToQPartialDiff)) * 448.8 / dt
        ELSE
            DO z = x , OpenholeFirstEl - 1       ! Annulus elements
                CALL PartialDerivativeFricToFlowRate(z)
                IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
                    WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
                    WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
                    CALL ErrorStop('NaN in calculating partial derivative')
                END IF
            END DO
            ExitValue = ExitValue - SUM(FlowEl(x : OpenholeFirstEl - 1)%FricToQPartialDiff) * 448.8 / dt
        END IF
    
    
    ELSE IF (CalcMode == 5) THEN
        IF (x >= OpenholeFirstEl .AND. y < OpenholeFirstEl) THEN     ! Gas Pocket is in Openhole and upper pocket is in annulus
            !WRITE (*,*) 'x , y 1' , x, y
            ExitValue = GasPocketNewPress%Array(GasPocketNo + 1) + SUM(FlowEl(x : NumbEl)%StaticPressDiff) + SUM(FlowEl(AnnulusFirstEl : y)%StaticPressDiff)
        ELSE                ! Both gas pockets are in Annulus or openhole
            !WRITE (*,*) 'x , y 2' , x, y
            ExitValue = GasPocketNewPress%Array(GasPocketNo + 1) + SUM(FlowEl(x : y)%StaticPressDiff)
        END IF
    END IF
        
  
    
    END SUBROUTINE
    
    
    
SUBROUTINE NewGasKick
    
    USE KickVariables
    Use TD_WellGeometry
    Use CReservoirVariables
    Use CFormationVariables
    USE Fluid_Flow_Startup_Vars
    USE PressureDisplayVARIABLES
    USE FricPressDropVars
    USE MudSystemVARIABLES
    USE CMudPropertiesVariables
    USE CError
    USE , INTRINSIC :: IEEE_ARITHMETIC
    
    
    !! Note: a subject that may be confusing is that when we use gauge pressure and when using absolute pressure?!
    !! all stated pressure are gauge pressure, so I do like this.
    !! only when we want to use a state equation (like ideal gas equation), we should use absolute equation and so I do this.
    !! Thus gas pocket pressure are all absolute pressure.
    
    IMPLICIT NONE
    
    INTEGER :: i , j , k , l
    
    IF (NOT(ALLOCATED(GasPocketWeight%Array))) THEN     ! 1st kick
            WRITE (*,*) ' New Influx 1'
            
            NoGasPocket = 1
            MudSystemDotNewInfluxNumber = MudSystemDotNewInfluxNumber + 1
            
            MudSystemDotNewInfluxElementCreated = 0
            KickOffBottom = .FALSE.
            

            CALL GasPocketOldPress%AddToFirst((BottomHolePress + StandardPress) * 1.d0)
            CALL GasPocketNewPress%AddToFirst((BottomHolePress + StandardPress) * 1.d0)
            CALL GasPocketOldTemp%AddToFirst(600.0)
            CALL GasPocketNewTemp%AddToFirst(600.0)
            CALL GasPocketOldVol%AddToFirst(0.d0)
            CALL GasPocketNewVol%AddToFirst(0.d0)
            CALL GasPocketdeltaVol%AddToFirst(0.0)
            CALL GasPocketFlowInduced%AddToFirst(0.0)
            CALL GasPocketModifiedVol%AddToFirst(0.0)
            CALL GasPocketWeight%AddToFirst(0.0)
            CALL GasPocketDensity%AddToFirst(2.0)
            CALL GasPocketCompressibility%AddToFirst(0.98)
            
            ALLOCATE(KickJacobian(2 , 2) , OldKickJacobian(2 , 2) , KickVandPFunction(2) , KickUnknownVector(2) , KickCorrectionVector(2))
            
            BottomHoleTemperature = 600
            KickFluxAvgPressure = (BottomHolePress + FormPressure) / 2 + StandardPress
            KickFluxAvgTemperature = (FormTemperature + BottomHoleTemperature) / 2
            KickFluxAvgCompressibility = 0.98
        
            !K_Aa = (5.8742362 * 10.**(-3) * KickFluxAvgTemperature**1.2288) / (511.1728532 + KickFluxAvgTemperature)
            !K_Bb = 5.5565586 + (1000.01 / KickFluxAvgTemperature)
            !K_Cc = 2.47862 - 0.12294 * K_Bb
            GasKickSIDensity = KickFluxAvgPressure / (KickFluxAvgCompressibility * &
                                        KickFluxAvgTemperature * GasType(KickGasType)%GasConstant) * Convpcftogpcm3
            GasKickDensity = GasKickSIDensity * 8.3523
    
            !GasKickViscosity = K_Aa * EXP(K_Bb * GasKickSIDensity**K_Cc)
            !
            !K_BHTpr = BottomHoleTemperature / KickTc
            !K_BHPpr = (BottomHolePress + StandardPress) / KickPc
            !
            !K_A_Bottomhole = 3.53 * K_BHPpr
            !K_B_Bottomhole = 10.0**(0.9813 * K_BHTpr)
            !K_C_Bottomhole = 0.274 * (K_BHPpr**2)
            !K_D_Bottomhole = 10.0**(0.8157 * K_BHTpr)
            !
            !BottomHoleCompressibility = 0.98 !1. - (K_A_Bottomhole / K_B_Bottomhole) + (K_C_Bottomhole / K_D_Bottomhole)
            !
            !GasKickBg = 0.00504 * KickFluxAvgCompressibility * KickFluxAvgTemperature / KickFluxAvgPressure  ![bbl/SCF]

            
            !KickmdotACoef = 10.**(-8) * 1.15741d0 * 7.08d0 * FormPermeability * REAL(KickFormLength) * GasType(KickGasType)%StDensity / &
                                !(GasKickViscosity * GasKickBg * LOG(10000.d0))
            !IF (WellHeadOpen)       KickmdotACoef = (1.0 + 2.0) * KickmdotACoef

            !KickmdotBCoef = FormPressure + StandardPress    !! - Sum(static and friction pressure loss) of flow elements below gas pocket, see below

            !GasPocketWeight%Array(1) = GasKickDensity * 0.05     !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            GasPocketWeight%Array(1) = GasKickDensity * MinKickVol   !1.0:seyyed gofte     !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            
            GasPocketNewVol%Array(1) = GasPocketCompressibility%Array(1) * GasType(KickGasType)%GasConstant * &
                                GasPocketNewTemp%Array(1) * GasPocketWeight%Array(1) / GasPocketNewPress%Array(1)
            
            GasPocketDeltaVol%Array(1) = 0.05    !GasPocketNewVol%Array(1)
            GasPocketFlowInduced%Array(1) = (GasPocketDeltaVol%Array(1)) / dt * 448.8       ! gpm
            GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
            
            WRITE (*,*) ' FormPressure , BottomHolePress' , FormPressure , BottomHolePress, GasKickDensity
            WRITE (*,*) ' No  Press(psia)  Vol(gal)  Weight(lbm)  Flow Induced(gpm)'
            DO i = 1 , NoGasPocket
                WRITE (*,102) i , GasPocketNewPress%Array(i), GasPocketNewVol%Array(i) * Convft3toUSgal, GasPocketWeight%Array(i), GasPocketFlowInduced%Array(i)
            END DO
            
    !ELSE IF (NoGasPocket < MaxGasPocket .AND. KickOffBottom .AND. (GasPocketNewVol%Array(1) > MinAllowableKickVol .OR. KickWasExitingThroughChoke)) THEN
    ELSE IF (NoGasPocket < MaxGasPocket .AND. KickOffBottom .AND. (GasPocketNewVol%Array(1) > MinAllowableKickVol .OR. ANY(GasPocketFlowEl(1 , :) == OpenholeFirstEl - 1))) THEN
            WRITE (*,*) ' New Influx', NoGasPocket + 1
            
102     FORMAT (I2, 4X, (F8.1), 3X, (F8.3), 2X, (F8.3), 8X, (F8.3))
            
            
            NoGasPocket = NoGasPocket + 1
            MudSystemDotNewInfluxNumber = MudSystemDotNewInfluxNumber + 1
            
            MudSystemDotNewInfluxElementCreated = 0
            KickOffBottom = .FALSE.

            CALL GasPocketOldPress%AddToFirst((BottomHolePress + StandardPress) * 1.d0)
            CALL GasPocketNewPress%AddToFirst((BottomHolePress + StandardPress) * 1.d0)
            CALL GasPocketOldTemp%AddToFirst(600.0)
            CALL GasPocketNewTemp%AddToFirst(600.0)
            CALL GasPocketOldVol%AddToFirst(0.d0)
            CALL GasPocketNewVol%AddToFirst(0.d0)
            CALL GasPocketdeltaVol%AddToFirst(0.0)
            CALL GasPocketFlowInduced%AddToFirst(0.0)
            CALL GasPocketModifiedVol%AddToFirst(0.0)
            CALL GasPocketWeight%AddToFirst(0.0)
            CALL GasPocketDensity%AddToFirst(2.0)
            CALL GasPocketCompressibility%AddToFirst(0.98)
            
            DEALLOCATE(KickJacobian , OldKickJacobian , KickVandPFunction , KickUnknownVector , KickCorrectionVector)
            
            ALLOCATE(KickJacobian(2 * NoGasPocket , 2 * NoGasPocket) , OldKickJacobian(2 * NoGasPocket , 2 * NoGasPocket))
            ALLOCATE(KickUnknownVector(2 * NoGasPocket) , KickCorrectionVector(2 * NoGasPocket) , KickVandPFunction(2 * NoGasPocket))
    
            
            BottomHoleTemperature = 600
            KickFluxAvgPressure = (BottomHolePress + FormPressure) / 2 + StandardPress
            KickFluxAvgTemperature = (FormTemperature + BottomHoleTemperature) / 2
            KickFluxAvgCompressibility = 0.98
        
            !K_Aa = (5.8742362 * 10.**(-3) * KickFluxAvgTemperature**1.2288) / (511.1728532 + KickFluxAvgTemperature)
            !K_Bb = 5.5565586 + (1000.01 / KickFluxAvgTemperature)
            !K_Cc = 2.47862 - 0.12294 * K_Bb
            GasKickSIDensity = KickFluxAvgPressure / (KickFluxAvgCompressibility * &
                                        KickFluxAvgTemperature * GasType(KickGasType)%GasConstant) * Convpcftogpcm3
            GasKickDensity = GasKickSIDensity * 8.3523

            !GasKickViscosity = K_Aa * EXP(K_Bb * GasKickSIDensity**K_Cc)
            !
            !K_BHTpr = BottomHoleTemperature / KickTc
            !K_BHPpr = (BottomHolePress + StandardPress) / KickPc
            !
            !K_A_Bottomhole = 3.53 * K_BHPpr
            !K_B_Bottomhole = 10.0**(0.9813 * K_BHTpr)
            !K_C_Bottomhole = 0.274 * (K_BHPpr**2)
            !K_D_Bottomhole = 10.0**(0.8157 * K_BHTpr)
            !
            !BottomHoleCompressibility = 0.98 !1. - (K_A_Bottomhole / K_B_Bottomhole) + (K_C_Bottomhole / K_D_Bottomhole)
            !
            !GasKickBg = 0.00504 * KickFluxAvgCompressibility * KickFluxAvgTemperature / KickFluxAvgPressure  ![bbl/SCF]
            !
            !
            !KickmdotACoef = 10.**(-8) * 1.15741d0 * 7.08d0 * FormPermeability * REAL(KickFormLength) * GasType(KickGasType)%StDensity / &
            !                    (GasKickViscosity * GasKickBg * LOG(10000.d0))
            !IF (WellHeadOpen)       KickmdotACoef = (1.0 + 2.0) * KickmdotACoef
            !
            !KickmdotBCoef = FormPressure + StandardPress    !! - Sum(static and friction pressure loss) of flow elements below gas pocket, see below

            !GasPocketWeight%Array(1) = GasKickDensity * 0.05        !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            GasPocketWeight%Array(1) = GasKickDensity * MinKickVol   !1.0:seyyed gofte        !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            
            GasPocketNewVol%Array(1) = GasPocketCompressibility%Array(1) * GasType(KickGasType)%GasConstant * &
                                GasPocketNewTemp%Array(1) * GasPocketWeight%Array(1) / GasPocketNewPress%Array(1)
            
            GasPocketDeltaVol%Array(1) = 0.05    !GasPocketNewVol%Array(1)
            GasPocketFlowInduced%Array(1) = (GasPocketDeltaVol%Array(1)) / dt * 448.8       ! gpm
            GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
            
            WRITE (*,*) ' FormPressure , BottomHolePress' , FormPressure , BottomHolePress, GasKickDensity
            WRITE (*,*) ' No  Press(psia)  Vol(gal)  Weight(lbm)  Flow Induced(gpm)'
            DO i = 1 , NoGasPocket
                WRITE (*,102) i , GasPocketNewPress%Array(i), GasPocketNewVol%Array(i) * Convft3toUSgal, GasPocketWeight%Array(i), GasPocketFlowInduced%Array(i)
            END DO
    
    ELSE    ! no new kick, so mass of 1st kick should increase        
        GasPocketWeight%Array(1) = GasPocketweight%Array(1) + KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
        GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
        
        IF (NoGasPocket > 1 .OR. SecondaryKickWeight > 0.0) THEN
            SecondaryKickWeight = SecondaryKickWeight + KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            SecondaryKickVol = SecondaryKickWeight / GasReservoirDensity / 42.0     ! 42 USGal = 1bbl
        END IF
        
    END IF
    
    
END SUBROUTINE