coneq.f

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

               DO 225 I = 1,NBLM1
                  ODT(I+1,1) = COE3(I)+COE4(I)+COE5(I)
  225          CONTINUE
            ENDIF
C
C     X-DISPERSION ACCUMULATE OVER PHASES
C
            DO 230 I = 1,NBLM1
               COEX(I+1) = COEX(I+1)+COE3(I)+COE4(I)+COE5(I)
  230       CONTINUE
C
C     Y-DIRECTION DERIVATIVES
C
            IF (NY.GT.1) THEN 
C
C     FLUX AND CONCENTRATION GRADIENT FOR THE CALCULATION
C     OF Y-X CROSS DISPERSION 
C
               IF (ICOORD.EQ.4) THEN
C
C  ---FLEXI GRID 
C     DCDX UPWIND IN Y DIR'N
C
                  DO 240 I = 2,NBLM2
                     LWKSP1(I)=.TRUE.
                     COE1(I) = 0.5*(VX(LY(I,L),L)+VX(LY(I,L)-1,L))     
                     CG1(I) = 2.*(C(LY(I,L)+1,KC,L)-C(LY(I,L)-1,KC,L))
     *                     /(DDX(LY(I,L))*(ONE+ALX1(LY(I,L)))+
     *                     DDX(LY(I,L)-1)*(ONE+ALX1(LY(I,L)-1)))
                     IF (S(LY(I,L)+1,L).LE.ONEM8.AND.S(LY(I,L)-1,L).GT. 
     *                  ONEM8) THEN 
                        CG1(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-1,KC,L))
     *                     /(DDX(LY(I,L)-1)*(ONE+ALX1(LY(I,L)-1)))
                        LWKSP1(I)=.FALSE.
                     ENDIF
  240             CONTINUE
                  DO 245 I = 2,NBLM2
                     IF (LWKSP1(I).AND.S(LY(I,L)-1,L).LE.ONEM8
     *                  .AND.S(LY(I,L)+1,L).GT.ONEM8) THEN 
                        CG1(I) = 2.*(C(LY(I,L)+1,KC,L)-C(LY(I,L),KC,L))
     *                     /(DDX(LY(I,L))*(ONE+ALX1(LY(I,L))))
                     ENDIF
  245             CONTINUE
C
C     B.C.   DCDX = 0 | X = L.H.S
C
                  DO 250 I = 1,NBLM2,NX
                     COE1(I) = 0.5*VX(LY(I,L),L)
                     CG1(I) = 2.*(C(LY(I,L)+1,KC,L)-C(LY(I,L),KC,L))
     *                    /(DDX(LY(I,L))*(ONE+ALX1(LY(I,L))))
                     IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)+1,L).LE.
     *               ONEM8) CG1(I)=0.0
  250             CONTINUE
C
C     B.C.   DCDX = 0 | X = R.H.S
C
                  DO 260 I = NX,NBLM2,NX
                     COE1(I) = 0.5*VX(LY(I,L)-1,L)
                     CG1(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-1,KC,L))
     *                    /(DDX(LY(I,L)-1)*(ONE+ALX1(LY(I,L)-1)))
                     IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)-1,L).LE.
     *               ONEM8) CG1(I)=0.0  
  260             CONTINUE
              ELSE
                  DO 243 I = 2,NBLM2
                     LWKSP1(I)=.TRUE.
                     COE1(I) = 0.5*(VX(LY(I,L),L)+VX(LY(I,L)-1,L))
                     CG1(I) = (C(LY(I,L)+1,KC,L)-C(LY(I,L)-1,KC,L))/
     *                (.5*(DDX(LY(I,L)+1)+DDX(LY(I,L)-1))+DDX(LY(I,L)))
                     IF (S(LY(I,L)+1,L).LE.ONEM8.AND.S(LY(I,L)-1,L).GT. 
     *                  ONEM8) THEN 
                        CG1(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-1,KC,L)
     *                    )/(DDX(LY(I,L))+DDX(LY(I,L)-1))             
                        LWKSP1(I)=.FALSE.
                     ENDIF
  243             CONTINUE
                  DO 248 I = 2,NBLM2
                     IF (LWKSP1(I).AND.S(LY(I,L)-1,L).LE.ONEM8
     *                  .AND.S(LY(I,L)+1,L).GT.ONEM8) THEN 
                        CG1(I) = 2.*(C(LY(I,L)+1,KC,L)-C(LY(I,L),KC,L))
     *                    /(DDX(LY(I,L))+DDX(LY(I,L)+1))
                     ENDIF
  248             CONTINUE
                  DO 255 I = 1,NBLM2,NX
                     COE1(I) = 0.5*VX(LY(I,L),L)
                     CG1(I) = 2.*(C(LY(I,L)+1,KC,L)-C(LY(I,L),KC,L))/
     *                       (DDX(LY(I,L))+DDX(LY(I,L)+1))
                     IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)+1,L).LE.
     *               ONEM8) CG1(I)=0.0
  255             CONTINUE
                  DO 265 I = NX,NBLM2,NX
                     COE1(I) = 0.5*VX(LY(I,L)-1,L)
                     CG1(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-1,KC,L))/
     *                      (DDX(LY(I,L))+DDX(LY(I,L)-1))
                     IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)-1,L).LE.
     *               ONEM8) CG1(I)=0.0  
  265             CONTINUE
               ENDIF
C
C     FLUX AND CONCENTRATION GRADIENT FOR THE CALCULATION
C     OF Y-Z CROSS DISPERSION 
C
               IF (NZ.GT.1) THEN
                  IF (ICOORD.EQ.4) THEN
C
C  ---FLEXI GRID 
C
                     DO 270 I = 1,NXNY
                        COE2(I) = 0.5*VZ(LY(I,L),L)
                        CG2(I) = 2.*(C(LY(I,L)+NXNY,KC,L)-C(LY(I,L)
     *                     ,KC,L))/(DDZ(LY(I,L))*(ONE+ALZ1(LY(I,L))))
                        IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)+NXNY,L)
     *                    .LE.ONEM8) CG2(I)=0.0
  270                CONTINUE
                     DO 280 I = (NZ-1)*NXNY+1,NBLM2
                        COE2(I) = 0.5*VZ(LY(I,L)-NXNY,L)
                        CG2(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-NXNY,KC,
     *                 L))/(DDZ(LY(I,L)-NXNY)*(ONE+ALZ1(LY(I,L)-NXNY)))
                        IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)-NXNY,L)
     *                     .LE.ONEM8) CG2(I)=0.0  
  280                CONTINUE
                     IF (NZ.GT.2) THEN
                        DO 290 I = NXNY+1,NBLM3
                           LWKSP1(I)=.TRUE.
                           COE2(I) = 0.5*(VZ(LY(I,L),L)+
     *                               VZ(LY(I,L)-NXNY,L)) 
                           CG2(I)=2.*
     *                      (C(LY(I,L)+NXNY,KC,L)-C(LY(I,L)-NXNY,KC,L))
     *                      /(DDZ(LY(I,L))*(ONE+ALZ1(LY(I,L)))+
     *                      DDZ(LY(I,L)-NXNY)*(ONE+ALZ1(LY(I,L)-NXNY)))
                           IF (S(LY(I,L)+NXNY,L).LE.ONEM8.AND.S(LY(I,L)
     *                        -NXNY,L).GT.ONEM8) THEN 
                              CG2(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-
     *               NXNY,KC,L))/(DDZ(LY(I,L)-NXNY)*(ONE+ALZ1(LY(I,L)
     *               -NXNY)))  
                              LWKSP1(I)=.FALSE.
                           ENDIF
  290                   CONTINUE 
                        DO 295 I = NXNY+1,NBLM3
                           IF (LWKSP1(I).AND.S(LY(I,L)-NXNY,L).LE.ONEM8
     *                        .AND.S(LY(I,L)+NXNY,L).GT.ONEM8) THEN 
                              CG2(I) = 2.*(C(LY(I,L)+NXNY,KC,L)-
     *                           C(LY(I,L),KC,L))/(DDZ(LY(I,L))*(ONE+
     *                           ALZ1(LY(I,L))))
                           ENDIF
  295                   CONTINUE
                     ENDIF
                  ELSE
C
C     FLUX AND CONCENTRATION GRADIENT FOR THE CALCULATION
C     OF Y-Z CROSS DISPERSION 
C
                     DO 275 I = 1,NXNY
                        COE2(I) = 0.5*VZ(LY(I,L),L)
                        CG2(I) = 2.*(C(LY(I,L)+NXNY,KC,L)-C(LY(I,L)
     *                      ,KC,L))/(DDZ(LY(I,L))+DDZ(LY(I,L)+NXNY))
                        IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)+NXNY,L)
     *                      .LE.ONEM8) CG2(I)=0.0
  275                CONTINUE
                     DO 285 I = (NZ-1)*NXNY+1,NBLM2
                        COE2(I) = 0.5*VZ(LY(I,L)-NXNY,L)
                        CG2(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-NXNY,
     *                     KC,L))/(DDZ(LY(I,L))+DDZ(LY(I,L)-NXNY))
                        IF (S(LY(I,L),L).LE.ONEM8.OR.S(LY(I,L)-NXNY,L)
     *                     .LE.ONEM8) CG2(I)=0.0  
  285                CONTINUE
                     IF (NZ.GT.2) THEN
                        DO 293 I = NXNY+1,NBLM3
                           LWKSP1(I)=.TRUE.
                           COE2(I) = 0.5*(VZ(LY(I,L),L)+VZ(LY(I,L)-
     *                               NXNY,L))
                           CG2(I) = (C(LY(I,L)+NXNY,KC,L)-C(LY(I,L)-
     *                        NXNY,KC,L))/(.5*(DDZ(LY(I,L)+NXNY)+
     *                        DDZ(LY(I,L)-NXNY))+DDZ(LY(I,L)))
                           IF (S(LY(I,L)+NXNY,L).LE.ONEM8.AND.S(LY(I,L)
     *                        -NXNY,L).GT.ONEM8) THEN 
                              CG2(I) = 2.*(C(LY(I,L),KC,L)-C(LY(I,L)-
     *                         NXNY,KC,L))/(DDZ(LY(I,L))+
     *                         DDZ(LY(I,L)-NXNY)) 
                               LWKSP1(I)=.FALSE.
                           ENDIF
  293                   CONTINUE 
                        DO 297 I = NXNY+1,NBLM3
                           IF (LWKSP1(I).AND.S(LY(I,L)-NXNY,L).LE.ONEM8
     *                        .AND.S(LY(I,L)+NXNY,L).GT.ONEM8) THEN 
                              CG2(I) = 2.*(C(LY(I,L)+NXNY,KC,L)-C(LY(I
     *                      ,L),KC,L))/(DDZ(LY(I,L))+DDZ(LY(I,L)+NXNY))
                           ENDIF
  297                   CONTINUE
                     ENDIF
                  ENDIF
               ELSE 
                  DO 300 I = 1,NBLM2
                     COE2(I) = 0.0
  300             CONTINUE
               ENDIF
C
C     CALCULATE DISPERSION DIAGONAL TERMS
C 
               DO 310 I = 1,NBLM2
                  COE3(I) = VY(I,L)**2
                  IF (ABS(VY(I,L)).LE.ONEM10)THEN
                    COE1(I)=0.0
                    COE2(I)=0.0
                  ENDIF
                  COE4(I) = COE1(I)**2
                  COE5(I) = COE2(I)**2
                  COE6(I) = SQRT(COE3(I)+COE4(I)+COE5(I))
                  COE6(I) = CVMGT(ONE199,COE6(I),COE6(I).LE. ONEM10)
                  COE6(I) = ONE/COE6(I)                                
  310          CONTINUE
               DO 320 I = 1,NBLM2
                  IF(KC.GE.MTWM1+1.AND.KC.LE.MTWM1+NTW)THEN  
                     COE7(I) = SF(LY(I,L),L)*PORC(LY(I,L))
                  ELSE
                     COE7(I) = S(LY(I,L),L)*PORC(LY(I,L))
                  ENDIF
  320          CONTINUE
C
C     ADD THE MOLECULAR DIFFUSION
C
               DO 330 I = 1,NBLM2
                  COE7(I) = COE7(I)*D(KC,L)+(ALPHAL(L)*COE3(I)+
     *                      ALPHAT(L)*(COE4(I)+COE5(I)))*COE6(I)
                  COE1(I) = ALPHAX*(COE1(I)*VY(I,L))*COE6(I)
                  COE2(I) = ALPHAX*(COE2(I)*VY(I,L))*COE6(I)
  330          CONTINUE
C
C     CONC. GRADIENT IN Y-DIRECTION
C
            DO 335 I = 1,NBLM2
               CG3(I) = 2.*(C(I+NX,KC,L)-C(I,KC,L))/(DDY(I+NX)+DDY(I))
               IF (S(I+NX,L).LE.ONEM8.OR.S(I,L).LE.ONEM8) CG3(I)=0.0
  335       CONTINUE          
C
C     CALCULATE NUMERICAL DISPERSION COEFFICIENT TERM
C
            DO 337 I = NX+1,NBLM2
               COE6(I) = 0.0
               IF (ITC.GE.1) COE6(I) = -.5*DT*VY(I,L)*VY(I,L)/PORC(I)
  337       CONTINUE
            DO 338 I = 1,NX
               COE6(I) = 0.0
               IF (ITC.GE.1) COE6(I) = -.5*DT*VY(I,L)*VY(I,L)/PORC(I)
  338       CONTINUE
C  
C     CALCULATE NUMERICAL DISPERSION CONTROL AND
C     DIAGONAL DISPERSION TERMS IN Y-DIRECTION
C
            DO 340 I = 1,NBLM2
               COE3(I) =(COE6(I)-COE7(I))*CG3(I) 
  340       CONTINUE
C
C     CALCULATE Y-X CROSS DISPERSION TERM
C
            IF (NY.GT.1) THEN 
               DO 350 I = 1,NBLM2 
                  COE4(I) = -COE1(I)*CG1(I)
                  IF (ITC.GE.1) COE4(I) = COE4(I)-.5*DT*VX(I,L)*VY(I,L)
     *               /PORC(I)*CG1(I)
  350          CONTINUE
            ELSE
               DO 360 I = 1,NBLM2
                  COE4(I) = 0.0
  360          CONTINUE
            ENDIF
C  
C     CALCLATE Y-Z CROSS DISPERSION TERM
C
            DO 370 I = 1,NBLM2
               COE5(I) = 0.0
  370       CONTINUE
            IF (NZ.GT.1) THEN 
               DO 380 I = 1,NBLM2 
                  COE5(I) = -COE2(I)*CG2(I)
                  IF (ITC.GE.1) COE5(I)=COE5(I)-.5*DT*VY(I,L)*VZ(I,L)
     *             /PORC(I)*CG2(I)
  380          CONTINUE
            ENDIF
C                             
C     Y-DISPERSION FOR OIL COMPONENT IN WATER PHASE OR SINGLE M.E
C     ME (TYPE II(-))
            IF (IMASS.EQ.2.AND.(KC.EQ.2.OR.(KC.GT.MMOM1.AND.
     *              KC.LE.MMOM1+NO)).AND.L.EQ.1) THEN
               DO 385 I = 1,NBLM2
                  ODT(I+NX,2) = COE3(I)+COE4(I)+COE5(I)
  385          CONTINUE
            ENDIF
C                             
C     Y-DISPERSION ACCUMULATE OVER PHASES
C
            DO 400 I = 1,NBLM2
               COEYU(I+NX) = COEYU(I+NX)+COE3(I)+COE4(I)+COE5(I)
  400       CONTINUE
         ENDIF
C
C     Z-DIRECTION DERIVATIVES
C
            IF (NZ.GT.1) THEN 
C
C     FLUX AND CONCENTRATION GRADIENT FOR THE CALCULATION
C     OF Z-X CROSS DISPERSION 
C
               IF (ICOORD.EQ.2) THEN
                  DO 405 I = 2,NBL-1
                     LWKSP1(I) = .TRUE.
                     COE1(I) = 0.5*(VX(LZ(I,L),L)+VX(LZ(I,L)-1,L))
                     CG1(I) = 2.*RR(I)*(C(LZ(I,L)+1,KC,L)-
     *                 C(LZ(I,L)-1,KC,L))/(RRSQ(I+1)-RRSQ(I-1))
                     IF (S(LZ(I,L)+1,L).LE.ONEM8.AND.S(LZ(I,L)-1,L).GT.
     *                 ONEM8) THEN 
                       CG1(I) = 2.*RR(I)*(C(LZ(I,L),KC,L)-
     *                   C(LZ(I,L)-1,KC,L))/(RRSQ(I)-RRSQ(I-1))             
                       LWKSP1(I) = .FALSE.
                     END IF
  405             CONTINUE
                  DO 410 I = 2,NBL-1
                     IF (LWKSP1(I).AND.S(LZ(I,L)-1,L).LE.ONEM8
     *                  .AND.S(LZ(I,L)+1,L).GT.ONEM8) THEN 
                       CG1(I) = 2.*RR(I)*(C(LZ(I,L)+1,KC,L)-
     *                        C(LZ(I,L),KC,L))/(RRSQ(I+1)-RRSQ(I))
                     ENDIF
  410             CONTINUE
                  DO 415 I = 1,NBL,NX
                     COE1(I) = 0.5*VX(LZ(I,L),L)
                       CG1(I) = 2.*RR(I)*(C(LZ(I,L)+1,KC,L)-
     *                        C(LZ(I,L),KC,L))/(RRSQ(I+1)-RRSQ(I))
                        IF (S(LZ(I,L),L).LE.ONEM8.OR.S(LZ(I,L)+1,L).LE.
     *                  ONEM8) CG1(I) = 0.0
  415             CONTINUE
                  DO 417 I = NX,NBL,NX
                     COE1(I) = 0.5*VX(LZ(I,L)-1,L)
                       CG1(I) = 2.*RR(I)*(C(LZ(I,L),KC,L)-
     *                   C(LZ(I,L)-1,KC,L))/(RRSQ(I)-RRSQ(I-1))             
                       IF (S(LZ(I,L),L).LE.ONEM8.OR.S(LZ(I,L)-1,L).LE.
     *                  ONEM8) CG1(I)=0.0  
  417             CONTINUE
               ELSE
                  IF (ICOORD.EQ.4) THEN
C
C  ---FLEXI GRID 
C
                     DO 420 I = 2,NBLM3
                        LWKSP1(I) = .TRUE.