coneq.f

油田化学驱模拟的经典 fortran源代码

C
      SUBROUTINE CONEQ
      USE MODULE1,ONLY :
     &    AD21,AD22,C2ADSS,IADSO
     &  , IREV
     &  , NALK1,NALK2,NALK3,NALK4,NALK5,NALK6
     &  , D,ALPHAL,ALPHAT,IMES
     &  , QL,IBOUND,IBL,IBR,IBT,IZONE
     &  , CD,SD,SF
     &  , CM,CTF,FFL,FFH,SLP
     &  , ICAPFL
     &  , PR,POLD,PORC,CORFP,CORFF
     &  , ZERO,ONE,ONEM,ONEM4,ONEM5,ONEM6,ONEM7,ONEM8,ONEM9
     &  , ONEM10,ONEM12,ONEP12,ONEM50,ONEP50,ONEM5M,PONEM,ONE199,PRCSN
     &  , PIE,F1P8
     &  , LX,LY,LZ
     &  , VBM,CSION,TSURF,PSURF
     &  , IREACT,ICAP,ICWM,IMODE,IBIO
     &  , COE6=>F6S,COE7=>F7S,COEX=>F8S,COEYU=>CA71
     &  , COEZT=>CA72,DDX=>CA73,DDY=>CA81,DDZ=>CA82,DC=>CA83
     &  , ALX1,BLX1,ALY1,BLY1
     &  , ALZ1,BLZ1,RPERMX                               
     &  , RPERMY,RPERMZ
     &  , AG1,AG2,CRG,AGK,BGK  
     &  , A15D,B15D,C15ADS,C14ADS,QW
     &  , CRNAK,HNAK,C16ADS,CNAADS,A14D,B14=>B14D,C160,ICREX
     &  , AK1,AK2,SCR,X4,X14,X16,WM4,GKIN,GRATE
     &  , AK1T,AK2T,X13,KGOPT
     &  , P4RW,E4W,S4RC,P4RC,E4C,T44,S4RW
     &  , SGI,IGAS,NPHAS
     &  , NGP1,NGP2,NGP3,NGP4,NG,NGP5,EPHI1
     &  ,  LWKSP2,LWKSP3,LWKSP4
     &  , CDC
     &  , TQLOS,RINO,RINU,TEM
     &  , TVHC,TCONO,TCONU,IHLOS,IANAL
     &  , CRTC,CVSPR,CVSPL,DENS,CUMHI,CUMHP,TEMPI,IENG
     &  , BVI,TEMINJ,DENO,DENU,CVSPO,CVSPU
     &  , ENTHE,TEMREF
     &  , CTOT,C,CSE,S
     &  , CE  
     &  , BRK,CRK,VIS1,VIS2,AP1,AP2,AP3,GAMMAC,EPHI4,EPHI3
     &  , GAMHF,SSLOPE,POWN,CSE1,ALPHAV,BETAP,TSTAND,VISW,IPOLYM
     &  , VIS,RPERM,PERMX,PERMY
     &  , PERMZ,QI,QB,Q,QT
     &  , CUMQI,CUMQP,PWF
     &  , EL,DX,DY,DZ,R,RP
     &  , RPSQ,IJKPOS,IPRESS
     &  , DT,CURANT,NXM1,NX,NY,NZ,NXNY,NBL,NBLW,N
     &  , PORMY,P,CPC,EPC
     &  , SSTAR,IOW 
     &  , TRSX,TRSY,TRSZ,TX
     &  , TY,TZ,CONVX,CONVY
     &  , CONVZ
     &  , VX=>VELX,VY=>VELY,VZ=>VELZ
     &  , AD31,AD32,B3D,C3ADSS,A3DS,FADS,EFFK
     &  , IADS1,IADK
     &  , AD41,AD42,B4D,C4ADSS,AD4
     &  , POR,RKF 
     &  , PSTAND,COMPR,COMPC,CTERM
     &  , C6JO,C6ADSS,C6HATS,QV,XKC,XKS,EQW
     &  , QVV
     &  , DDYY
     &  , LWKSP1
     &  , MTWM1,MTAM1,MGLM1
     &  , LTW,LTA,LGC,LACD,LGL
     &  , MGCM1,NGC
     &  , MMOM1,NO,LMO
     &  , AKOCK,C2ADSK,DCOK
     &  , LCAL2
     &  , DCO,WSOL,CNEM2,IMASS,ISOL,ICOR
     &  , FOREC,FOREC1,RESPV,RESATK,BTO,VB
     &  , CUMI,CUMP,OIP,OP,T,TINJ,WHPV
     &  , PRF,ICNT,IINJ,INEC,IRST
     &  , DCMAX,IDISPC,ICF,ICOORD,ITC,IUNIT
     &  , DNX,DNZ
     &  , TITLE,RELERR,SWI
     &  , TK,TTRP,TREC,TBT
     &  , RDC,TRD,RET,CINJT,TMW,TDEN
     &  , TAK,RCOR,ITRU,NT,ITREAC,NRT
     &  , TKS,TKT,TAKT,C5INI,NTW,NTA
     &  , DELC,DCMIN,ICOUNT
     &  , COE1=>WKSP1,COE2=>WKSP2,COE3=>WKSP3
     &  , COE4=>WKSP4,COE5=>WKSP5
     &  , RR=>AW,RRSQ=>AE,RRP=>AN,RRPSQ=>AS
     &  , REDT=>AT,GELT=>AB,TEMPK=>AC,DUM8=>BV
     &  , TWS1
     &  , CONVO=>TWS2,ODT=>TWS3
     &  , TWS4,TWS5
     &  , CG1=>DUM1,CG2=>DUM2,CG3=>DUM3,HION=>DUM4
     &  , DUM5,DUM6,DUM7
     &  , NWELL,IRO,IDW,IFLAG,IDIR
     &  , IWC,JWC,KWC,NWBC
     &  , RW,SWELL,KPRF,IPRF
     &  , IW,JW,IFIRST,ILAST
     &  , NH,NNA,NCA,NMG,NCARB,NACD,NALU,NSILI,NOXY      
     &  , NELET,NFLD,NSLD,NSORB,NACAT,NIAQ,NEX
     &  , NKEX,NFLASH,NSLEL,NSORBX,NIND,NIND1
     &  , ICHRGE,NSURF1,NSURF2,IRSPS 
          USE MODULE1,ONLY :
     &    AR,BR,DR
     &  , ER,BB,EXSLD
     &  , CHARGE,SCHARG,CHACAT
     &  , EXCAI,EXCA,EXEX
     &  , EXACAT,EQK,EXK
     &  , SPK,EXKL,SPKL,ACATK
     &  , CSLDI,CSORBI,CELAQI
     &  , CAC2I,CI,ACIDIS,CSLDT,CSORBT
     &  , CAQI,C50,C60,CAQSP,CACATT
     &  , CACAT,C1I,C2I
     &  , INCMEL
     &  , ELCRG
     &  , NALK8,NHFD,NCRFD
     &  , VFRACM,DMX,DMY,DMZ
     &  , XLSUB,YLSUB,ZLSUB
     &  , VSUB,PORCM
     &  , CSUB
     &  , DXL,DYL,DZL,VNUM
     &  , CT,THM,TVM
C
C     ------------------------------------------------------------
C     PURPOSE : THIS ROUTINE SOLVES THE CONVECTION-DIFFUSION 
C               (CONCENTRATION) EQUATION EXPLICITLY AND CALCULATES 
C               THE RELATIVE ERROR FOR EACH COMPONENT.
C
C     CALL    : GEL
C     ------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     ------------------------------------------------------------------
      DIMENSION TEMP(N),TAKNEW(NBLW,2)
      DIMENSION DC2(NBL)
      EXTERNAL ISAMAX
   
C
      EPS = ONEM8
      EPS1 = ONEM10
      EPS2 = ONEM12
      WMCR = 51.99
      WMOH = 17.0
      THREE = 3.0
C
C     CUTOFF FOR DNAPL CONC. IN WATER BY EPA
C
      EPAC = 3.08*ONEM9
      IJK = 0 
      NBLM1 = NBL-1
      NBLM2 = NBL-NX
      NBLM3 = NBL-NXNY
C
C
C     1) NORMALIZE ALL CONCENTRATIONS FROM PER UNIT FLUID VOLUME
C        BASIS TO PER UNIT PORE VOLUME BASIS IF THERE IS SURFACTANT
C        (COMPONENT 3) ADSORPTION OR ORGANIC ADSORPTION.
C        IN SUBROUTINE ADSORB, WE CHANGE BACK TO UNIT VOL. BASIS
      IF (ICNT.EQ.1) GO TO 1132
      DO 1100 I = 1,NBL
         COEX(I) = 1.
1100  CONTINUE
C
      IF (IADSO.GT.0) THEN
         DO 1105  I = 1,NBL
            COEX(I) = COEX(I)-C2ADSS(I)
1105     CONTINUE
      ENDIF
      IF (ICF(3).NE.0.AND.ABS(AD31).GE.PRCSN) THEN
          DO 1110 I = 1,NBL
            COEX(I) = COEX(I)-C3ADSS(I)
1110      CONTINUE
      ENDIF
C
      DO 1115  KC = 1,N
         IF (ICF(KC).EQ.0) GO TO 1115
         DO 1120 I = 1,NBL
            CTOT(I,KC) = CTOT(I,KC)*COEX(I) 
1120     CONTINUE                      
1115  CONTINUE
C
C     2) ADD THE ADSORBED AMOUNT TO THE FLUID PART OF POLYMER
C        AND SURFACTACNT BEFORE SOLVING THE CONCENTRATION EQ.
C
      IF (ICF(3).NE.0.AND.ABS(AD31).GE.PRCSN) THEN
         DO 3 I = 1,NBL
            CTOT(I,3) = CTOT(I,3)+C3ADSS(I)
  3      CONTINUE
      ENDIF
C
      IF (ICF(4).EQ.1) THEN
         DO 4 I = 1,NBL
            CTOT(I,4) = CTOT(I,4)+C4ADSS(I)
  4      CONTINUE
      ENDIF
C
      IF (IADSO.GT.0) THEN
         DO 14 I = 1,NBL
            CTOT(I,2) = CTOT(I,2)+C2ADSS(I)
14       CONTINUE
      ENDIF
      IF (IADSO.GT.0.AND.LMO) THEN
         DO 1130 K = 1,NO
            KC = K + MMOM1
            DO 1125 I = 1,NBL
               CTOT(I,KC) = CTOT(I,KC)+C2ADSK(I,K)
 1125       CONTINUE
 1130    CONTINUE
      ENDIF
 1132 CONTINUE
C
C    INITIATE TOTAL OIL FOR MULTIPLE OIL OPTION
C
      LCAL2=.FALSE.
      IF (LMO) THEN
         CUMP(2) = 0.
         CUMI(2) = 0.
         OP(2) = 0.
         DO 1135 I = 1,NBL
            CTOT(I,2) = 0.
            DC2(I) = 0.
 1135    CONTINUE
      ENDIF
C
C     3) FOR THE GEL OPTION, CALCULATE THE SUM OF FLUID AND ADSORBED 
C        CHROMIUM, GEL, AND HYDROGEN CONC. BEFORE SOLVING THE 
C        CONSERVATION EQUATIONS
C
      IF (IREACT.EQ.1.OR.IREACT.EQ.4) THEN
         DO 5 I = 1,NBL
            IF (IREACT.EQ.4) THEN
               CTOT(I,NGP3) = CAQSP(I,NCRFD)*(1000.*WMCR)*CTOT(I,1)
            ENDIF
            CTOT(I,NGP3) = CTOT(I,NGP3)+C14ADS(I)
            CTOT(I,NGP4) = CTOT(I,NGP4)+C15ADS(I)
  5      CONTINUE
C
C     TEMP. DEPENDENT REACTION RATES FOR GEL OPTION
C                 
         DO 801 I = 1,NBL
            GELT(I) = AK2
            REDT(I) = AK1
 801     CONTINUE
C
         IF (IENG.EQ.1 ) THEN
            IF (ABS(AK2T).GT.ZERO) THEN
               DO 802 I = 1,NBL
                  TEMPK(I) = (TEM(I)+459.67)/F1P8
                  GELT(I) = AK2*EXP(AK2T*(ONE/TEMPK(I)-ONE/TEMREF))
 802           CONTINUE
            ENDIF
            IF (ABS(AK1T).GT.ZERO) THEN
               DO 803 I = 1,NBL
                  TEMPK(I) = (TEM(I)+459.67)/F1P8
                  REDT(I) = AK1*EXP(AK1T*(ONE/TEMPK(I)-ONE/TEMREF))
 803           CONTINUE
            ENDIF
         ENDIF
      ENDIF
C
      IF (IREACT.EQ.1.AND.ICREX.EQ.1) THEN
         DO 6 I = 1,NBL
            CTOT(I,NGP5) = CTOT(I,NGP5)+C16ADS(I)
  6      CONTINUE                 
      ENDIF
C
C     ASSIGN HYDROGEN ION CONC. (MEQ/ML WATER) TO HION VARIABLE
C
      IF (IREACT.EQ.1) THEN
         DO 7 I = 1,NBL
            HION(I) = C(I,NGP5,1)
7        CONTINUE
      ELSEIF (IREACT.EQ.4) THEN 
         DO 8 I = 1,NBL
            HION(I) =  CAQSP(I,NHFD)
            C(I,NGP3,1) = CAQSP(I,NCRFD)*(1000.*WMCR)
8        CONTINUE
      ENDIF
C
      IF (ICAP.EQ.0) THEN
         DO 9 L = 1,NPHAS
         DO 9 I = 1,NBL
            SF(I,L) = S(I,L)
    9    CONTINUE
      ENDIF
C
C     REACTING TRACERS
C
      IF (ITREAC.EQ.1) THEN
         DO 2000 IT = 1,NRT
            IF (IENG.EQ.1.AND.ABS(TAKT(IT)).GT.ZERO) THEN
               DO 2010 I = 1,NBL
                  TEMPK(I) = (TEM(I)+459.67)/F1P8
                  TAKNEW(I,IT) = TAK(IT)*EXP(TAKT(IT)*(ONE/TEMPK(I)
     *                            -ONE/TEMREF))
 2010         CONTINUE
            ELSE
              DO 2020 I = 1,NBL
                 TAKNEW(I,IT) = TAK(IT)
 2020         CONTINUE
            ENDIF
 2000    CONTINUE
      ENDIF
C
      IF (ICOORD.NE.3) THEN
C
C  ---FLEXI GRID 
C
         IF (ICOORD.EQ.4) THEN
            DO 11 J = 1,NY
               IK = 0
            DO 11 K = 1,NZ
            DO 11 I = 1,NX
               IJK = I + (K-1)*NXNY + (J-1)*NX
               IK = IK+1
               DDX(IJK) = DNX(IK)
               DDY(IJK) = DY(J)
               DDZ(IJK) = DNZ(IK)
  11       CONTINUE
         ELSE 
           DO 10 K = 1,NZ
           DO 10 J = 1,NY
           DO 10 I = 1,NX
              IJK = IJK+1
              DDX(IJK) = DX(I)
              DDY(IJK) = DY(J)
              DDZ(IJK) = DZ(K)
   10      CONTINUE
         ENDIF
      ELSE 
         IJK = 0 
         DO 12 K = 1,NZ
         DO 12 I = 1,NX
            IJK = IJK+1
            DDX(IJK) = DX(I)
            DDY(IJK) = DDYY(I)
            DDZ(IJK) = DZ(K)
   12    CONTINUE
      ENDIF
C
C     CALCULATE THE DISTANCES FOR THE RADIAL GEOMETRY
C
      IF (ICOORD.EQ.2) THEN
         IK = 0
         DO 15 K = 1,NZ
            DO 15 I = 1,NX
               IK = IK+1
               RR(IK) = R(I)