SimulationCore2/TorqueDrag/TD_DrillingSubs/MeshGeneration_FluidModule.f90

subroutine MeshGeneration_FluidModule

   Use TD_DrillStemComponents
   use SimulationVariables !@
   use SimulationVariables !@
   Use TD_StringConnectionData
   Use GeoElements_FluidModule
   use CPumpsVariables
    use CPumps
   Use CStringConfigurationVariables
   Use CBopStackVariables
   !use ConfigurationVariables !@ , only: Drawworks

   implicit none

   Integer    :: ii , jj , semijj , kk , k , m , s
   Integer    :: ElementsCount , StringConfigCount
   REAL(8)    :: mm , nn , dl , StartAngle , EndAngle
   REAL(8)    :: A(30)    !?????????10
   REAL(8)    :: TD_ElementLength
  
   
   
   
   A    = 0.d0
   A(1) = data%Configuration%BopStack%AboveAnnularHeight   ! WellHead[ft]
   jj   = 1
   

   
   
   
   
   
!====================================================
!      Mesh Generation of DrillStem Components
!====================================================
   
   !data%State%TD_String%StringConfigurationCount = StringConfigurationCount    !???????????
   StringConfigCount = data%State%TD_String%StringConfigurationCount
   if (allocated(data%State%F_String)) Deallocate(data%State%F_String)
   if (data%State%TD_DrillStem(1)%ComponentType==0) then
       Allocate (data%State%F_String(StringConfigCount-1))
   else
       Allocate (data%State%F_String(StringConfigCount))
   end if
   
   
   if (data%State%TD_DrillStem(1)%ComponentType==0) then
       ElementsCount = 1
       Do ii=1,(StringConfigCount-1)
           data%State%F_String(ii)%ID           = data%State%TD_DrillStem(ii+1)%Id*12.d0    ![inch]
           data%State%F_String(ii)%OD           = data%State%TD_DrillStem(ii+1)%Od*12.d0    ![inch]
           data%State%F_String(ii)%FirstElement = ElementsCount+1
           data%State%F_String(ii)%LastElement  = data%State%F_String(ii)%FirstElement+(data%State%TD_DrillStem(ii+1)%Numbs-1)
           data%State%F_String(ii)%ElType       = data%State%TD_DrillStem(ii+1)%ComponentType
           ElementsCount             = data%State%F_String(ii)%LastElement
       End Do
       StringConfigCount = StringConfigCount-1
   else
       ElementsCount = 0
       Do ii=1,StringConfigCount
           data%State%F_String(ii)%ID           = data%State%TD_DrillStem(ii)%Id*12.d0      ![inch]
           data%State%F_String(ii)%OD           = data%State%TD_DrillStem(ii)%Od*12.d0      ![inch]
           data%State%F_String(ii)%FirstElement = ElementsCount+1
           data%State%F_String(ii)%LastElement  = data%State%F_String(ii)%FirstElement+(data%State%TD_DrillStem(ii)%Numbs-1)
           data%State%F_String(ii)%ElType       = data%State%TD_DrillStem(ii)%ComponentType
           ElementsCount             = data%State%F_String(ii)%LastElement
       End Do
   end if
   
   
   Do ii=1,StringConfigCount
       data%State%F_String(ii)%TopDepth  = data%State%TD_DrillStems(data%State%F_String(ii)%LastElement)%TopDepthIni    ![ft]
       data%State%F_String(ii)%DownDepth = data%State%TD_DrillStems(data%State%F_String(ii)%FirstElement)%DownDepthIni  ![ft]
       
       if (data%State%F_String(ii)%DownDepth>A(1)) then
           jj    = jj+1
           A(jj) = data%State%F_String(ii)%DownDepth
       end if
   End Do
   
   
   
   !=========> Removed Volume Calculation in DrillStem
   TD_ElementLength = data%State%F_String(StringConfigCount)%DownDepth-data%State%F_String(StringConfigCount)%TopDepth
   if ( StringConfigCount==data%State%TD_Vol%PreCount ) then
       if ( data%State%TD_Vol%PreElementLength>TD_ElementLength ) then
           data%State%TD_Vol%RemoveVolume = (data%State%TD_Vol%PreElementLength-TD_ElementLength)*((pi*((data%State%F_String(StringConfigCount)%ID/12.d0)**2))/4.d0)  ![ft^3]
       else
           data%State%TD_Vol%RemoveVolume = 0.d0
       end if
   else if ( StringConfigCount<data%State%TD_Vol%PreCount ) then
       data%State%TD_Vol%RemoveVolume = data%State%TD_Vol%PreElementVolume
   else
       data%State%TD_Vol%RemoveVolume = 0.d0
   end if
   data%State%TD_Vol%PreCount         = StringConfigCount
   data%State%TD_Vol%PreElementVolume = TD_ElementLength*(((pi*((data%State%F_String(StringConfigCount)%ID/12.d0)**2))/4.d0))   ![ft^3]
   data%State%TD_Vol%PreElementLength = TD_ElementLength
   
   
   !print*, 'data%State%TD_Vol%RemoveVolume=' , data%State%TD_Vol%RemoveVolume
   
   
   
   
   
!====================================================
!      Mesh Generation of Annulus Components
!====================================================
   
   if (data%State%TD_Casing(1)%Length>0.) then
       jj    = jj+1
       A(jj) = data%State%TD_Casing(1)%DownDepth
       data%State%OD_Annulus(4)%StartMD = data%State%TD_Casing(1)%DownDepth      ![ft]
       data%State%OD_Annulus(4)%EndMD   = data%State%TD_Casing(1)%TopDepth       ![ft]
       data%State%OD_Annulus(4)%ODValue = data%State%TD_Casing(1)%Id*12.d0       ![inch]
   !else
   !    data%State%OD_Annulus(4)%StartMD = 0.                      !????????????????????
   !    data%State%OD_Annulus(4)%EndMD   = data%State%OD_Annulus(4)%StartMD   !????????????????????
   !    data%State%OD_Annulus(4)%ODValue = data%State%TD_Casing(1)%Id         !????????????????????
   end if
   
   if (data%State%TD_Liner(1)%Length>0.) then
       jj    = jj+1
       A(jj) = data%State%TD_Liner(1)%DownDepth
       data%State%OD_Annulus(3)%StartMD = data%State%TD_Liner(1)%DownDepth       ![ft]
       data%State%OD_Annulus(3)%EndMD   = data%State%TD_Liner(1)%TopDepth        ![ft]
       data%State%OD_Annulus(3)%ODValue = data%State%TD_Liner(1)%Id*12.d0        ![inch]
   else
       data%State%OD_Annulus(3)%StartMD = data%State%TD_Casing(1)%DownDepth   !????????????????????
       data%State%OD_Annulus(3)%EndMD   = data%State%OD_Annulus(3)%StartMD    !????????????????????
       data%State%OD_Annulus(3)%ODValue = data%State%OD_Annulus(4)%ODValue    !????????????????????
   end if
   
   if (data%State%TD_OpenHole(1)%Length>0.) then
       !jj    = jj+1
       !A(jj) = data%State%TD_OpenHole(1)%DownDepth
       data%State%OD_Annulus(2)%StartMD = data%State%TD_OpenHole(1)%DownDepth      ![ft]
       data%State%OD_Annulus(2)%EndMD   = data%State%TD_OpenHole(1)%TopDepth       ![ft]
       data%State%OD_Annulus(2)%ODValue = data%State%TD_OpenHole(1)%Id*12.d0       ![inch]
       do ii = 1,data%State%TD_WellGeneral%WellIntervalsCount-1
           if ( data%State%TD_WellGeo(ii)%DownDepth>data%State%TD_OpenHole(1)%TopDepth ) then
               jj    = jj+1
               A(jj) = data%State%TD_WellGeo(ii)%DownDepth
           end if
       end do
   else
       data%State%OD_Annulus(2)%StartMD = data%State%TD_Liner(1)%DownDepth   !??????????????
       data%State%OD_Annulus(2)%EndMD   = data%State%OD_Annulus(2)%StartMD   !??????????????
       data%State%OD_Annulus(2)%ODValue = data%State%OD_Annulus(3)%ODValue   !??????????????
   end if
   
   if (data%State%TD_ROPHole(1)%Length>0.) then
       jj    = jj+1
       A(jj) = data%State%TD_ROPHole(1)%DownDepth
       data%State%OD_Annulus(1)%StartMD = data%State%TD_ROPHole(1)%DownDepth      ![ft]
       data%State%OD_Annulus(1)%EndMD   = data%State%TD_ROPHole(1)%TopDepth       ![ft]
       data%State%OD_Annulus(1)%ODValue = data%State%TD_ROPHole(1)%Id*12.d0       ![inch]
   else
       data%State%OD_Annulus(1)%StartMD = data%State%TD_OpenHole(1)%DownDepth  !??????????????
       data%State%OD_Annulus(1)%EndMD   = data%State%OD_Annulus(1)%StartMD     !??????????????
       data%State%OD_Annulus(1)%ODValue = data%State%OD_Annulus(2)%ODValue     !??????????????
   end if
   
   !print*, 'A(jj)2=' , A
   
   
   
   
   
   
   
!====================================================
!           Mesh Generation of Well
!====================================================
   !Do ii=2,data%State%TD_WellGeneral%WellIntervalsCount-1
   !    jj = jj+1
   !    A(jj) = data%State%TD_WellGeo(ii)%TopDepth
   !End Do
   
   
   !print*, 'A(jj)3=' , A
   
   !print*, 'A=' , A
   !print*, 'jj=' , jj
   
   semijj = jj
   !print*, 'semijj=' , semijj
   Do ii=1,jj
       do kk = 1,jj
           if ( A(ii)==A(KK) .and. ii/=kk .and. A(KK)/=0. ) then !A(KK)/=0. ?????????????????
               !print*, 'indoA=' , semijj , ii , kk , A(KK) , A(ii)
               A(kk)  = 0.d0
               semijj = semijj-1
               !print*, 'indoA2=' , semijj , ii , kk , A(KK) , A(ii)
           end if
       end do
   End Do
   !print*, 'semijj2=' , semijj
   jj = semijj
   
   !print*, 'jj2=' , jj
   !print*, 'AA=' , A
   
   
   
   
   
   if (allocated(data%State%TDGeo%MD)) Deallocate(data%State%TDGeo%MD)
   Allocate (data%State%TDGeo%MD(jj))
   data%State%TDGeo%MD = 0.0d0
   Do ii=1,jj
        data%State%TDGeo%MD(ii)=A(1)
       do kk = 2,jj
           if (A(KK)>data%State%TDGeo%MD(ii)) then
               data%State%TDGeo%MD(ii) = A(kk)
               k=kk
           end if
       end do
       !!data%State%TDGeo%MD(ii) = MAXVAL(A)
       !!kk     = MAXLOC(A)
       !!A(kk)  = 0.0d0
        A(k)    = 0.0d0
   End Do

   
   
   
   
   
   
   
!====================================================
!  Geometrical Intervals Configuration (for fluid module)
!====================================================
   
   if (allocated(data%State%F_Interval)) Deallocate(data%State%F_Interval)
   Allocate (data%State%F_Interval(jj-1+StringConfigCount+1))
   !data%State%F_Interval            = 0.0d0
   data%State%F_Counts%IntervalsTotalCounts = jj-1+StringConfigCount+1
   
!=========> String Intervals Configuration
   ! pump output
   ii = 1
   data%State%F_Interval(ii)%Number      = ii
   data%State%F_Interval(ii)%GeoType     = 0
   data%State%F_Interval(ii)%EndDepth    = data%State%F_String(StringConfigCount)%TopDepth
   data%State%F_Interval(ii)%StartDepth  = data%State%F_String(StringConfigCount)%TopDepth-265.d0
   data%State%F_Interval(ii)%ID          = 0.0d0
   data%State%F_Interval(ii)%OD          = dmax1(data%Configuration%Pumps%MudPump1Output,data%Configuration%Pumps%MudPump2Output)   !???????????
   data%State%F_Interval(ii)%HydDiameter = data%State%F_Interval(ii)%OD    ![inch]
   data%State%F_Interval(ii)%Volume      = ((pi*((data%State%F_Interval(ii)%OD/12.d0)**2))/4.d0) &
                                *(data%State%F_Interval(ii)%EndDepth-data%State%F_Interval(ii)%StartDepth)*7.48051948d0     ![gal-us]
   
   
   Do ii=2,StringConfigCount+1
       data%State%F_Interval(ii)%Number      = ii
       data%State%F_Interval(ii)%GeoType     = 0
       data%State%F_Interval(ii)%EndDepth    = data%State%F_String(StringConfigCount+1-ii+1)%DownDepth
       data%State%F_Interval(ii)%StartDepth  = data%State%F_String(StringConfigCount+1-ii+1)%TopDepth
       data%State%F_Interval(ii)%ID          = 0.0d0
       data%State%F_Interval(ii)%OD          = data%State%F_String(StringConfigCount+1-ii+1)%ID
       data%State%F_Interval(ii)%HydDiameter = data%State%F_Interval(ii)%OD
       data%State%F_Interval(ii)%Volume      = ((pi*((data%State%F_Interval(ii)%OD/12.d0)**2))/4.d0) &
                                    *(data%State%F_Interval(ii)%EndDepth-data%State%F_Interval(ii)%StartDepth)*7.48051948d0  ![gal-us]
   End Do
   
   data%State%F_Counts%StringIntervalCounts = StringConfigCount+1
   
   !if (data%State%F_String(1)%ElType==0) then   !bit
   !    !data%State%F_Interval(StringConfigCount)%Number     = ii
   !    data%State%F_Interval(StringConfigCount)%EndDepth    = data%State%F_String(1)%DownDepth
   !    data%State%F_Interval(StringConfigCount)%StartDepth  = data%State%F_String(1)%TopDepth
   !    data%State%F_Interval(StringConfigCount)%ID          = 0.0d0
   !    data%State%F_Interval(StringConfigCount)%OD          = data%State%F_String(1)%OD
   !    data%State%F_Interval(StringConfigCount)%HydDiameter = data%State%F_Interval(StringConfigCount)%OD
   !    data%State%F_Interval(StringConfigCount)%Volume      = ((pi*((data%State%F_Interval(StringConfigCount)%OD/12.)**2))/4.)*(data%State%F_Interval(StringConfigCount)%EndDepth-data%State%F_Interval(StringConfigCount)%StartDepth)*7.48051948
   !end if

!=========> Annulus Intervals Configuration
   data%State%F_Counts%BottomHoleIntervalCounts = 0
   data%State%F_Counts%AnnulusIntervalCounts    = 0
   Do ii=(StringConfigCount+2),(jj-1+StringConfigCount+1)
       data%State%F_Interval(ii)%Number     = ii
       data%State%F_Interval(ii)%StartDepth = data%State%TDGeo%MD(ii-StringConfigCount-1)
       data%State%F_Interval(ii)%EndDepth   = data%State%TDGeo%MD(ii+1-StringConfigCount-1)
       if (data%State%F_Interval(ii)%EndDepth>=data%State%F_Interval(StringConfigCount+1)%EndDepth) then
           data%State%F_Interval(ii)%GeoType     = 1
           data%State%F_Counts%BottomHoleIntervalCounts = data%State%F_Counts%BottomHoleIntervalCounts+1
       else
           data%State%F_Interval(ii)%GeoType  = 2
           data%State%F_Counts%AnnulusIntervalCounts = data%State%F_Counts%AnnulusIntervalCounts+1
       end if
   End Do
   
   
   
   
   
   
   
   
   
   
   
   
!====================================================
!      Diameter Calculation of Annulus Intervals
!====================================================
   
!=========> OD
   m = 1
   Do ii = (StringConfigCount+2),data%State%F_Counts%IntervalsTotalCounts
       do kk = m,4
           if (data%State%F_Interval(ii)%StartDepth<=data%State%OD_Annulus(kk)%StartMD .and. data%State%F_Interval(ii)%EndDepth>=data%State%OD_Annulus(kk)%EndMD) then
               data%State%F_Interval(ii)%OD = data%State%OD_Annulus(kk)%ODValue
               m = kk
               exit
           end if
       end do
   End Do
   
   
   
!=========> ID
   m = 1
   Do ii = (StringConfigCount+2),data%State%F_Counts%IntervalsTotalCounts
       do kk = m,StringConfigCount
           if (data%State%F_Interval(ii)%StartDepth<=data%State%F_String(kk)%DownDepth .and. data%State%F_Interval(ii)%EndDepth>=data%State%F_String(kk)%TopDepth) then
               data%State%F_Interval(ii)%ID = data%State%F_String(kk)%OD
               m = kk
               exit
           else
               data%State%F_Interval(ii)%ID = 0.0d0
           end if
       end do
       data%State%F_Interval(ii)%HydDiameter = data%State%F_Interval(ii)%OD-data%State%F_Interval(ii)%ID
       data%State%F_Interval(ii)%Volume      = ((pi*(((data%State%F_Interval(ii)%OD/12.d0)**2)-((data%State%F_Interval(ii)%ID/12.d0)**2)))/4.d0) &
                                    *abs(data%State%F_Interval(ii)%EndDepth-data%State%F_Interval(ii)%StartDepth)*7.48051948d0     ![gal-us]
   End Do
   
   
   !print*, 'StringConfigCount=' , StringConfigCount
   !print*, 'data%State%F_Counts%IntervalsTotalCounts=' , data%State%F_Counts%IntervalsTotalCounts
   !print*, 'data%State%F_Counts%StringIntervalCounts=' , data%State%F_Counts%StringIntervalCounts
   !print*, 'data%State%F_Counts%BottomHoleIntervalCounts=' , data%State%F_Counts%BottomHoleIntervalCounts
   !print*, 'data%State%F_Counts%AnnulusIntervalCounts=' , data%State%F_Counts%AnnulusIntervalCounts
   !Do ii=1,(data%State%F_Counts%IntervalsTotalCounts)
   !    print*, 'Number/' , 'StartDepth/' , 'EndDepth/' 
   !    print*, 'ID/' , 'OD/' , 'GeoType/' 
   !    print*, 'HydDiameter/' , 'Volume/'
   !    print*, data%State%F_Interval(ii)%Number , data%State%F_Interval(ii)%StartDepth , data%State%F_Interval(ii)%EndDepth 
   !    print*, data%State%F_Interval(ii)%ID , data%State%F_Interval(ii)%OD , data%State%F_Interval(ii)%GeoType 
   !    print*, data%State%F_Interval(ii)%HydDiameter , data%State%F_Interval(ii)%Volume
   !End Do
   
   
   
   
   
   
   
!====================================================
!      TVD Calculation of Geo Intervals
!====================================================
   
   !=========> Out of Well Intervals
   data%State%F_Interval(1)%StartTVD  = 0.0d0   !??????????????
   data%State%F_Interval(1)%EndTVD    = -data%State%TD_StConn%ConnectionHeight
   data%State%F_Interval(1)%StartAngle = 0.d0   !???????????????????????
   data%State%F_Interval(1)%EndAngle   = 0.d0   !???????????????????????
   data%State%F_Counts%OutOfWellIntervalCounts = 1
   Do ii= 2 , StringConfigCount+1
       if ( data%State%F_Interval(ii)%StartDepth<=0.d0 ) then
           data%State%F_Interval(ii)%StartTVD = data%State%F_Interval(ii)%StartDepth
           data%State%F_Interval(ii)%StartAngle = 0.0d0
       end if
       if ( data%State%F_Interval(ii)%EndDepth<=0.d0 ) then
           data%State%F_Interval(ii)%EndTVD   = data%State%F_Interval(ii)%EndDepth
           data%State%F_Interval(ii)%EndAngle   = 0.0d0
           data%State%F_Counts%OutOfWellIntervalCounts = data%State%F_Counts%OutOfWellIntervalCounts+1
       end if
   End Do
   
   
   !=========>
   if (allocated(data%State%TDGeo%TVD)) Deallocate(data%State%TDGeo%TVD) 
   Allocate (data%State%TDGeo%TVD(jj))
   if (allocated(data%State%TDGeo%Angle)) Deallocate(data%State%TDGeo%Angle) 
   Allocate (data%State%TDGeo%Angle(jj))
   data%State%TDGeo%TVD(1)  = data%State%TD_WellGeneral%WellTotalVerticalLength !?????
   !print*, 'data%State%TD_WellGeneral%WellTotalVerticalLength=' ,data%State%TD_WellGeneral%WellTotalVerticalLength
   data%State%TDGeo%TVD(jj) = data%Configuration%BopStack%AboveAnnularHeight   ! WellHead[ft]      !0.d0     ????????????????
   
   k        = 1
   mm       = 0.d0
   nn       = 0.d0 
   EndAngle = data%State%TD_WellGeo(1)%StartAngle    !???????????????
   data%State%TDGeo%Angle(jj) = EndAngle
   Do ii = jj-1,1,-1   !???1or2
       do kk = k, data%State%TD_WellGeneral%WellIntervalsCount
           StartAngle = EndAngle
           if ( data%State%TDGeo%MD(ii)>data%State%TD_WellGeo(kk)%TopDepth ) then
               if ( data%State%TDGeo%MD(ii)>data%State%TD_WellGeo(kk)%DownDepth ) then
                   dl = data%State%TD_WellGeo(kk)%DownDepth-nn                       ![ft]
                   nn = data%State%TD_WellGeo(kk)%DownDepth
                   if ( data%State%TD_WellGeo(kk)%HoleType==0 ) then
                       EndAngle = StartAngle
                       !print*, 'StartAngle1=' ,ii , kk,  StartAngle
                       !print*, 'EndAngle1=' ,ii , kk,  EndAngle
                       data%State%TDGeo%TVD(ii)  = mm+(dl*cos(data%State%TD_WellGeo(kk)%StartAngle))  ![ft]
                       data%State%TDGeo%Angle(ii) = EndAngle
                       !print*, 'data%State%TDGeo%TVD(ii)1=' ,ii , kk,  data%State%TDGeo%TVD(ii)
                   else if ( data%State%TD_WellGeo(kk)%HoleType==1 ) then
                       EndAngle = StartAngle+(dl/data%State%TD_WellGeo(kk)%RCurvature) !?????????????????
                       !print*, 'StartAngle2=' ,ii , kk,  StartAngle
                       !print*, 'EndAngle2=' ,ii , kk,  EndAngle
                       data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(EndAngle)-abs(StartAngle))*cos(abs(StartAngle)))-(data%State%TD_WellGeo(kk)%RCurvature*(1.-cos(abs(EndAngle)-abs(StartAngle)))*sin(abs(StartAngle)))
                       !data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(EndAngle)-abs(StartAngle)))
                       data%State%TDGeo%Angle(ii) = EndAngle
                       !print*, 'data%State%TDGeo%TVD(ii)2=' , ii , kk, data%State%TDGeo%TVD(ii)
                   else if ( data%State%TD_WellGeo(kk)%HoleType==2 ) then
                       EndAngle = StartAngle-(dl/data%State%TD_WellGeo(kk)%RCurvature) !?????????????????
                       !print*, 'StartAngle22=' ,ii , kk,  StartAngle
                       !print*, 'EndAngle22=' ,ii , kk,  EndAngle
                       data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(abs(EndAngle)-abs(StartAngle)))*cos(abs(StartAngle)))+(data%State%TD_WellGeo(kk)%RCurvature*(1.-cos(abs(abs(EndAngle)-abs(StartAngle))))*sin(abs(StartAngle)))
                       !data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(EndAngle)-abs(StartAngle)))
                       data%State%TDGeo%Angle(ii) = EndAngle
                       !print*, 'data%State%TDGeo%TVD(ii)22=' , ii , kk, data%State%TDGeo%TVD(ii)
                   end if
                   mm = data%State%TDGeo%TVD(ii)
               else
                   dl = data%State%TDGeo%MD(ii)-nn
                   nn = data%State%TDGeo%MD(ii)
                   if ( data%State%TD_WellGeo(kk)%HoleType==0 ) then
                       EndAngle = StartAngle
                       !print*, 'StartAngle3=' ,ii , kk,  StartAngle
                       !print*, 'EndAngle3=' ,ii , kk,  EndAngle
                       data%State%TDGeo%TVD(ii)  = mm+(dl*cos(data%State%TD_WellGeo(kk)%StartAngle))
                       data%State%TDGeo%Angle(ii) = EndAngle
                       !print*, 'data%State%TDGeo%TVD(ii)3=' ,ii , kk,  data%State%TDGeo%TVD(ii)
                   else if ( data%State%TD_WellGeo(kk)%HoleType==1 ) then
                       EndAngle = StartAngle+(dl/data%State%TD_WellGeo(kk)%RCurvature) !???????????????????
                       !print*, 'StartAngle4=' ,ii , kk,  StartAngle
                       !print*, 'EndAngle4=' ,ii , kk,  EndAngle
                       data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(EndAngle)-abs(StartAngle))*cos(abs(StartAngle)))-(data%State%TD_WellGeo(kk)%RCurvature*(1.-cos(abs(EndAngle)-abs(StartAngle)))*sin(abs(StartAngle)))
                       data%State%TDGeo%Angle(ii) = EndAngle
                       !data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(EndAngle)-abs(StartAngle)))
                       !print*, 'mm=' ,mm
                       !print*, 'data%State%TDGeo%TVD(ii)4=' ,ii , kk,  data%State%TDGeo%TVD(ii)
                   else if ( data%State%TD_WellGeo(kk)%HoleType==2 ) then
                       EndAngle = StartAngle-(dl/data%State%TD_WellGeo(kk)%RCurvature) !???????????????????
                       !print*, 'StartAngle44=' ,ii , kk,  StartAngle
                       !print*, 'EndAngle44=' ,ii , kk,  EndAngle
                       data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(abs(EndAngle)-abs(StartAngle)))*cos(abs(StartAngle)))+(data%State%TD_WellGeo(kk)%RCurvature*(1.-cos(abs(abs(EndAngle)-abs(StartAngle))))*sin(abs(StartAngle)))
                       data%State%TDGeo%Angle(ii) = EndAngle
                       !data%State%TDGeo%TVD(ii)  = mm+(data%State%TD_WellGeo(kk)%RCurvature*sin(abs(EndAngle)-abs(StartAngle)))
                       !print*, 'mm=' ,mm , (data%State%TD_WellGeo(kk)%RCurvature*sin(abs(abs(EndAngle)-abs(StartAngle)))*cos(abs(StartAngle)))+(data%State%TD_WellGeo(kk)%RCurvature*(1.-cos(abs(abs(EndAngle)-abs(StartAngle))))*sin(abs(StartAngle)))
                       !print*, 'data%State%TDGeo%TVD(ii)44=' ,ii , kk,  data%State%TDGeo%TVD(ii)
                   end if
                   mm = data%State%TDGeo%TVD(ii)
                   k  = kk
                   exit
               end if
           end if
       end do
   End Do
        
    
   
   !if (data%State%F_Counts%OutOfWellIntervalCounts==1) then
   !    s = 2
   !else
   !    s = data%State%F_Counts%OutOfWellIntervalCounts
   !end if
   !
   
    Do ii = (data%State%F_Counts%OutOfWellIntervalCounts+1),data%State%F_Counts%IntervalsTotalCounts
        do kk = jj,1,-1
            if ( data%State%F_Interval(ii)%StartDepth==data%State%TDGeo%MD(kk) ) then
                data%State%F_Interval(ii)%StartTVD = data%State%TDGeo%TVD(kk)      ![ft]
                data%State%F_Interval(ii)%StartAngle = data%State%TDGeo%Angle(kk)  ![rad]
            end if
            if ( data%State%F_Interval(ii)%EndDepth==data%State%TDGeo%MD(kk) ) then
                data%State%F_Interval(ii)%EndTVD = data%State%TDGeo%TVD(kk)        ![ft]
                data%State%F_Interval(ii)%EndAngle = data%State%TDGeo%Angle(kk)    ![rad]
            end if
        end do
    End Do
    
    
   
   !Do ii=1,(data%State%F_Counts%IntervalsTotalCounts)
   !    print*, 'data%State%F_Interval(ii)%StartTVD=' , ii , data%State%F_Interval(ii)%StartTVD
   !    print*, 'data%State%F_Interval(ii)%EndTVD=' , ii ,  data%State%F_Interval(ii)%EndTVD
   !    print*, 'data%State%F_Interval(ii)%StartAngle=' , ii , data%State%F_Interval(ii)%StartAngle
   !    print*, 'data%State%F_Interval(ii)%EndAngle=' , ii ,  data%State%F_Interval(ii)%EndAngle
   !end do
   !    
   ! 
   !Do ii=1,data%State%TD_WellGeneral%WellIntervalsCount
   !    print*, 'data%State%TD_WellGeo(kk)%TopDepth=' , ii ,   data%State%TD_WellGeo(ii)%TopDepth
   !    print*, 'data%State%TD_WellGeo(kk)%DownDepth=' , ii ,   data%State%TD_WellGeo(ii)%DownDepth
   !    print*, 'data%State%TD_WellGeo(kk)%HoleType=' , ii ,   data%State%TD_WellGeo(ii)%HoleType
   !    print*, 'data%State%TD_WellGeo(kk)%RCurvature=' , ii ,   data%State%TD_WellGeo(ii)%RCurvature
   !    print*, 'data%State%TD_WellGeo(kk)%EndAngle=' , ii ,   data%State%TD_WellGeo(ii)%EndAngle
   !    print*, 'data%State%TD_WellGeo(kk)%StartAngle=' , ii ,   data%State%TD_WellGeo(ii)%StartAngle
   !end do 
   
   
   
   
   
end subroutine