SimulationCore2/FluidFlow/kick/Gas_Kick_Calculator.f90

SUBROUTINE GasKickCalculator
    use SimulationVariables
    use KickVARIABLESModule
    use SimulationVariables !@
    Use CReservoirVariables
    Use CFormationVariables
    USE Fluid_Flow_Startup_Vars
    use PressureDisplayVARIABLESModule
    USE FricPressDropVarsModule
    USE MudSystemVARIABLES
use SimulationVariables !@@@
    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
    
    KickVARIABLES%SolvingEquationError = .FALSE.
    KickVARIABLES%KickVandPFunction(:) = 0.d0
    KickVARIABLES%KickJacobian(: , :) = 0.d0
    
!====================================================
!            Gas Properties (Methane Gas)
!====================================================
    !GasPocketNewTemp%Array(1) = 600
    KickVARIABLES%BottomHoleTemperature = 600
    KickVARIABLES%KickFluxAvgPressure = (KickVARIABLES%BottomHolePress + KickVARIABLES%FormPressure) / 2 + StandardPress
    KickVARIABLES%KickFluxAvgTemperature = (KickVARIABLES%FormTemperature + KickVARIABLES%BottomHoleTemperature) / 2
    KickVARIABLES%KickFluxAvgCompressibility = 0.98d0
        
    KickVARIABLES%K_Aa = (5.8742362 * 10.**(-3) * KickVARIABLES%KickFluxAvgTemperature**1.2288) / (511.1728532 + KickVARIABLES%KickFluxAvgTemperature)
    KickVARIABLES%K_Bb = 5.5565586 + (1000.01 / KickVARIABLES%KickFluxAvgTemperature)
    KickVARIABLES%K_Cc = 2.47862 - 0.12294 * KickVARIABLES%K_Bb
    KickVARIABLES%GasKickSIDensity = KickVARIABLES%KickFluxAvgPressure / (KickVARIABLES%KickFluxAvgCompressibility * &
                                KickVARIABLES%KickFluxAvgTemperature * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant) * Convpcftogpcm3
    
    KickVARIABLES%GasKickViscosity = KickVARIABLES%K_Aa * EXP(KickVARIABLES%K_Bb * KickVARIABLES%GasKickSIDensity**KickVARIABLES%K_Cc)
    !!!!!!!!!!!!!!!!!!!!!!!!!
    
    !!!!!!!!!!!!!! Calculating compressibility for bottom hole condition

    
    KickVARIABLES%BottomHoleCompressibility = 0.98d0 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    KickVARIABLES%GasKickBg = 0.00504 * KickVARIABLES%KickFluxAvgCompressibility * KickVARIABLES%KickFluxAvgTemperature / KickVARIABLES%KickFluxAvgPressure  ![bbl/SCF]


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

    DO l = 1 , KickVARIABLES%NoGasPocket
        KickVARIABLES%KickUnknownVector(2 * l - 1) = GasPocketNewVol%Array(l)
        KickVARIABLES%KickUnknownVector(2 * l) = GasPocketNewPress%Array(l)
    END DO
    
    IF (KickVARIABLES%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 , KickVARIABLES%NoGasPocket
                KickVARIABLES%KickUnknownVector(2 * l - 1) = KickVARIABLES%KickUnknownVector(2 * l - 1) + GasPocketDeltaVol%Array(l)
            END DO
        END IF
        
        KickVARIABLES%KickVandPFunction(1) = KickVARIABLES%KickUnknownVector(1) - GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * &      ! VandP(1) = V(1)
                                GasPocketNewTemp%Array(1) * (GasPocketWeight%Array(1) + KickVARIABLES%KickmdotACoef * MAX(((KickVARIABLES%KickmdotBCoef - KickVARIABLES%KickUnknownVector(2)) * dt) , 0.0)) / KickVARIABLES%KickUnknownVector(2)
        !WRITE (*,*) 'KickVandPFunction(1)',KickVandPFunction(1)
        l = 2 * KickVARIABLES%NoGasPocket
        CALL KickFunctionsCalculator(KickVARIABLES%KickVandPFunction(l) , KickVARIABLES%NoGasPocket , 2)    ! VandP(last) = P(last)
        !WRITE (*,*) 'KickVandPFunction(l)', l, KickVandPFunction(l)
        DO l = 2 , KickVARIABLES%NoGasPocket      ! VandP(Odd) = V(l, l > 1)
            KickVARIABLES%KickVandPFunction(2 * l - 1) = KickVARIABLES%KickUnknownVector(2 * l - 1) - GasPocketCompressibility%Array(l) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * &
                                GasPocketNewTemp%Array(l) * GasPocketWeight%Array(l) / KickVARIABLES%KickUnknownVector(2 * l)
            !WRITE(*,*) 'KickVandPFunction(V)', l, KickVandPFunction(2 * l - 1)
        END DO
        
        DO l = KickVARIABLES%NoGasPocket - 1 , 1 , -1
            CALL KickFunctionsCalculator(KickVARIABLES%KickVandPFunction(2 * l) , l , 1)    ! VandP(Even) = P(l, l < KickVARIABLES%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 * KickVARIABLES%NoGasPocket      ! Main Diagonal
            KickVARIABLES%KickJacobian(k , k) = 1.d0
        END DO
        
        KickVARIABLES%KickJacobian(1,2) = (GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * GasPocketNewTemp%Array(1) &
                                * (GasPocketWeight%Array(1) + KickVARIABLES%KickmdotACoef * KickVARIABLES%KickmdotBCoef * dt) / KickVARIABLES%KickUnknownVector(2)**2)      ! Row 1 Finished
        IF (KickVARIABLES%KickJacobian(1,2) == 0.d0) THEN
            CALL Error('KickJacobian(1,2) = 0.0')
            KickVARIABLES%KickJacobian(1,2) = KickVARIABLES%OldKickJacobian(1,2)
        END IF
        
        !WRITE(*,*) 'KickJacobian(1,2)', KickJacobian(1,2)
        
        l = KickVARIABLES%NoGasPocket
        CALL KickFunctionsCalculator(KickVARIABLES%KickJacobian(2 * l , 2 * l - 1) , KickVARIABLES%NoGasPocket , 4)     ! for last Row
        IF (KickVARIABLES%KickJacobian(2 * l , 2 * l - 1) == 0.d0) THEN
            CALL Error ('KickJacobian(Last,Odd) = 0.0')
            KickVARIABLES%KickJacobian(2 * l , 2 * l - 1) = KickVARIABLES%OldKickJacobian(2 * l , 2 * l - 1)
        END IF
        
        DO k = KickVARIABLES%NoGasPocket - 1 , 1 , -1
            KickVARIABLES%KickJacobian(2 * l , 2 * k - 1) = KickVARIABLES%KickJacobian(2 * l , 2 * l - 1)
        END DO      ! Last Row Finished
        !WRITE(*,*) 'KickJacobian(2,1)', KickJacobian(2,1)
        
        
        DO k = 2 , KickVARIABLES%NoGasPocket
            KickVARIABLES%KickJacobian(2 * k - 1 , 2 * k) = GasPocketCompressibility%Array(k) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * GasPocketNewTemp%Array(k) &
                                * GasPocketWeight%Array(k) / KickVARIABLES%KickUnknownVector(2 * k)**2
            
        END DO      ! Odd Rows (V equations) Finished
        
        DO k = 1 , KickVARIABLES%NoGasPocket - 1
            KickVARIABLES%KickJacobian(2 * k , 2 * k + 2) = -1.d0
        END DO      ! Even Rows (P equations) effect of upper pocket
        
        DO k = 2 , 2 * (KickVARIABLES%NoGasPocket - 1) , 2
            DO l = 1 , k - 1 , 2
                CALL KickFunctionsCalculator(KickVARIABLES%KickJacobian(k , l) , k / 2 , 3)
                IF (KickVARIABLES%KickJacobian(k , l) == 0.d0) THEN
                    WRITE (*,*) 'Jacobian Array = 0.0', k , l
                    CALL Error ('KickJacobian(k , l) = 0.0')
                    KickVARIABLES%KickJacobian(k , l) = KickVARIABLES%OldKickJacobian(k , l)
                END IF
            END DO
        END DO
        
        IF (ANY(IEEE_IS_NaN(KickVARIABLES%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
        KickVARIABLES%KickVandPFunction = -1.d0 * KickVARIABLES%KickVandPFunction
        CALL SOLVE_LINEAR_EQUATIONS(KickVARIABLES%KickJacobian , KickVARIABLES%KickCorrectionVector , KickVARIABLES%KickVandPFunction , KickVARIABLES%SolvingEquationError, SIZE(KickVARIABLES%KickCorrectionVector))
        IF (KickVARIABLES%SolvingEquationError)       CALL ErrorStop( ' Error in solving kick equation ' )
        
        KickVARIABLES%KickUnknownVector = KickVARIABLES%KickUnknownVector + KickVARIABLES%KickCorrectionUnderRelaxation * KickVARIABLES%KickCorrectionVector
        
        DO l = 1 , KickVARIABLES%NoGasPocket
            GasPocketNewVol%Array(l) = KickVARIABLES%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) = KickVARIABLES%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
        
        
        
    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) + KickVARIABLES%KickmdotACoef * KickVARIABLES%KickmdotBCoef * dt) / (GasPocketWeight%Array(1) + KickVARIABLES%KickmdotACoef * GasPocketOldPress%Array(1) * dt)))
                
        
        
    END IF
        
    !DO l = 1 , KickVARIABLES%NoGasPocket
        GasPocketDeltaVol%Array(:) = GasPocketNewVol%Array(:) - GasPocketOldVol%Array(:)
        GasPocketFlowInduced%Array(:) = (GasPocketDeltaVol%Array(:)) / dt * 448.8       ! gpm
    !END DO
    
    KickVARIABLES%GasKickPumpFlowRate = 0.0
    IF (NOT(KickVARIABLES%KickOffBottom) .AND. KickVARIABLES%WellHeadOpen)     KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
    


    END SUBROUTINE
    
    
    SUBROUTINE KickFunctionsCalculator(ExitValue , GasPocketNo , CalcMode)
    
    use SimulationVariables
    use KickVARIABLESModule
    USE FricPressDropVarsModule
    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 = KickVARIABLES%GasPocketFlowEl(GasPocketNo , 1)
    IF (GasPocketNo < KickVARIABLES%NoGasPocket)      y = KickVARIABLES%GasPocketFlowEl(GasPocketNo + 1 , 1)
    i = data%State%FricPressDrop%StringLastEl
    j = data%State%FricPressDrop%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 = KickVARIABLES%KickUnknownVector(2 * GasPocketNo) - KickVARIABLES%KickUnknownVector(2 * GasPocketNo + 2)
        IF (x >= data%State%FricPressDrop%OpenholeFirstEl .AND. y < data%State%FricPressDrop%OpenholeFirstEl) THEN            ! Case 2 , Case 3
            ExitValue = ExitValue - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%StaticPressDiff) - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%FricPressLoss)   &
                                  - SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : y)%StaticPressDiff) - SUM(FlowEl(data%State%FricPressDrop%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 = KickVARIABLES%KickUnknownVector(2 * GasPocketNo) - StandardPress - data%State%FricPressDrop%Kchoke * FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%Flowrate**2
        IF (x >= data%State%FricPressDrop%OpenholeFirstEl) THEN     ! Gas Pocket is in Openhole
            ExitValue = ExitValue - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%StaticPressDiff) - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%FricPressLoss)   &
                                    - SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%StaticPressDiff) - SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FricPressLoss)
        ELSE                ! Gas Pocket is in Annulus
            ExitValue = ExitValue - SUM(FlowEl(x : data%State%FricPressDrop%OpenholeFirstEl - 1)%StaticPressDiff) - SUM(FlowEl(x : data%State%FricPressDrop%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 >= data%State%FricPressDrop%OpenholeFirstEl .AND. y < data%State%FricPressDrop%OpenholeFirstEl) THEN            ! Top kick STARTX is in Annulus
            DO z = x , data%State%FricPressDrop%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 = data%State%FricPressDrop%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 : data%State%FricPressDrop%NumbEl)%FricToQPartialDiff) + SUM(FlowEl(data%State%FricPressDrop%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 * data%State%FricPressDrop%Kchoke * FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%Flowrate * 448.8 / dt
        IF (x >= data%State%FricPressDrop%OpenholeFirstEl) THEN            ! kick STARTX is in openhole
            DO z = x , data%State%FricPressDrop%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 = data%State%FricPressDrop%AnnulusFirstEl , data%State%FricPressDrop%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 : data%State%FricPressDrop%NumbEl)%FricToQPartialDiff) + SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FricToQPartialDiff)) * 448.8 / dt
        ELSE
            DO z = x , data%State%FricPressDrop%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 : data%State%FricPressDrop%OpenholeFirstEl - 1)%FricToQPartialDiff) * 448.8 / dt
        END IF
    
    
    ELSE IF (CalcMode == 5) THEN
        IF (x >= data%State%FricPressDrop%OpenholeFirstEl .AND. y < data%State%FricPressDrop%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 : data%State%FricPressDrop%NumbEl)%StaticPressDiff) + SUM(FlowEl(data%State%FricPressDrop%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 SimulationVariables
    use KickVARIABLESModule
    use SimulationVariables !@
    Use CReservoirVariables
    Use CFormationVariables
    USE Fluid_Flow_Startup_Vars
    use PressureDisplayVARIABLESModule
    USE FricPressDropVarsModule
    USE MudSystemVARIABLES
use SimulationVariables !@@@
    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'
            
            KickVARIABLES%NoGasPocket = 1
            data%State%MudSystem%NewInfluxNumber = data%State%MudSystem%NewInfluxNumber + 1
            
            data%State%MudSystem%NewInfluxElementCreated = 0
            KickVARIABLES%KickOffBottom = .FALSE.
            

            CALL GasPocketOldPress%AddToFirst((KickVARIABLES%BottomHolePress + StandardPress) * 1.d0)
            CALL GasPocketNewPress%AddToFirst((KickVARIABLES%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(KickVARIABLES%KickJacobian(2 , 2) , KickVARIABLES%OldKickJacobian(2 , 2) , KickVARIABLES%KickVandPFunction(2) , KickVARIABLES%KickUnknownVector(2) , KickVARIABLES%KickCorrectionVector(2))
            
            KickVARIABLES%BottomHoleTemperature = 600
            KickVARIABLES%KickFluxAvgPressure = (KickVARIABLES%BottomHolePress + KickVARIABLES%FormPressure) / 2 + StandardPress
            KickVARIABLES%KickFluxAvgTemperature = (KickVARIABLES%FormTemperature + KickVARIABLES%BottomHoleTemperature) / 2
            KickVARIABLES%KickFluxAvgCompressibility = 0.98
        

            KickVARIABLES%GasKickSIDensity = KickVARIABLES%KickFluxAvgPressure / (KickVARIABLES%KickFluxAvgCompressibility * &
                                        KickVARIABLES%KickFluxAvgTemperature * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant) * Convpcftogpcm3
            KickVARIABLES%GasKickDensity = KickVARIABLES%GasKickSIDensity * 8.3523
            GasPocketWeight%Array(1) = KickVARIABLES%GasKickDensity * KickVARIABLES%MinKickVol   !1.0:seyyed gofte     !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            
            GasPocketNewVol%Array(1) = GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%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
            KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
            
            WRITE (*,*) ' FormPressure , BottomHolePress' , KickVARIABLES%FormPressure , KickVARIABLES%BottomHolePress, KickVARIABLES%GasKickDensity
            WRITE (*,*) ' No  Press(psia)  Vol(gal)  Weight(lbm)  Flow Induced(gpm)'
            DO i = 1 , KickVARIABLES%NoGasPocket
                WRITE (*,102) i , GasPocketNewPress%Array(i), GasPocketNewVol%Array(i) * Convft3toUSgal, GasPocketWeight%Array(i), GasPocketFlowInduced%Array(i)
            END DO
            
    !ELSE IF (NoGasPocket < MaxGasPocket .AND. KickVARIABLES%KickOffBottom .AND. (GasPocketNewVol%Array(1) > KickVARIABLES%MinAllowableKickVol .OR. KickWasExitingThroughChoke)) THEN
    ELSE IF (KickVARIABLES%NoGasPocket < KickVARIABLES%MaxGasPocket .AND. KickVARIABLES%KickOffBottom .AND. (GasPocketNewVol%Array(1) > KickVARIABLES%MinAllowableKickVol .OR. ANY(KickVARIABLES%GasPocketFlowEl(1 , :) == data%State%FricPressDrop%OpenholeFirstEl - 1))) THEN
            WRITE (*,*) ' New Influx', KickVARIABLES%NoGasPocket + 1
            
102     FORMAT (I2, 4X, (F8.1), 3X, (F8.3), 2X, (F8.3), 8X, (F8.3))
            
            
            KickVARIABLES%NoGasPocket = KickVARIABLES%NoGasPocket + 1
            data%State%MudSystem%NewInfluxNumber = data%State%MudSystem%NewInfluxNumber + 1
            
            data%State%MudSystem%NewInfluxElementCreated = 0
            KickVARIABLES%KickOffBottom = .FALSE.

            CALL GasPocketOldPress%AddToFirst((KickVARIABLES%BottomHolePress + StandardPress) * 1.d0)
            CALL GasPocketNewPress%AddToFirst((KickVARIABLES%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(KickVARIABLES%KickJacobian , KickVARIABLES%OldKickJacobian , KickVARIABLES%KickVandPFunction , KickVARIABLES%KickUnknownVector , KickVARIABLES%KickCorrectionVector)
            
            ALLOCATE(KickVARIABLES%KickJacobian(2 * KickVARIABLES%NoGasPocket , 2 * KickVARIABLES%NoGasPocket) , KickVARIABLES%OldKickJacobian(2 * KickVARIABLES%NoGasPocket , 2 * KickVARIABLES%NoGasPocket))
            ALLOCATE(KickVARIABLES%KickUnknownVector(2 * KickVARIABLES%NoGasPocket) , KickVARIABLES%KickCorrectionVector(2 * KickVARIABLES%NoGasPocket) , KickVARIABLES%KickVandPFunction(2 * KickVARIABLES%NoGasPocket))
    
            
            KickVARIABLES%BottomHoleTemperature = 600
            KickVARIABLES%KickFluxAvgPressure = (KickVARIABLES%BottomHolePress + KickVARIABLES%FormPressure) / 2 + StandardPress
            KickVARIABLES%KickFluxAvgTemperature = (KickVARIABLES%FormTemperature + KickVARIABLES%BottomHoleTemperature) / 2
            KickVARIABLES%KickFluxAvgCompressibility = 0.98
        

            KickVARIABLES%GasKickSIDensity = KickVARIABLES%KickFluxAvgPressure / (KickVARIABLES%KickFluxAvgCompressibility * &
                                        KickVARIABLES%KickFluxAvgTemperature * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant) * Convpcftogpcm3
            KickVARIABLES%GasKickDensity = KickVARIABLES%GasKickSIDensity * 8.3523

            !GasPocketWeight%Array(1) = KickVARIABLES%GasKickDensity * 0.05        !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            GasPocketWeight%Array(1) = KickVARIABLES%GasKickDensity * KickVARIABLES%MinKickVol   !1.0:seyyed gofte        !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            
            GasPocketNewVol%Array(1) = GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%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
            KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
            
            WRITE (*,*) ' FormPressure , BottomHolePress' , KickVARIABLES%FormPressure , KickVARIABLES%BottomHolePress, KickVARIABLES%GasKickDensity
            WRITE (*,*) ' No  Press(psia)  Vol(gal)  Weight(lbm)  Flow Induced(gpm)'
            DO i = 1 , KickVARIABLES%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) + KickVARIABLES%KickmdotACoef * (KickVARIABLES%KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
        KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
        
        IF (KickVARIABLES%NoGasPocket > 1 .OR. KickVARIABLES%SecondaryKickWeight > 0.0) THEN
            KickVARIABLES%SecondaryKickWeight = KickVARIABLES%SecondaryKickWeight + KickVARIABLES%KickmdotACoef * (KickVARIABLES%KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
            KickVARIABLES%SecondaryKickVol = KickVARIABLES%SecondaryKickWeight / KickVARIABLES%GasReservoirDensity / 42.0     ! 42 USGal = 1bbl
        END IF
        
    END IF
    
    
END SUBROUTINE