      SUBROUTINE PARLOD                                                   B  10
      DOUBLE PRECISION DMIN1,DEXP,DLOG,DABS                               B  20
      DOUBLE PRECISION TMRX,VPRM,HI,HR,HC,HK,WT,REC,RECH,TIM,AOPT,TITLE   B  30
      DOUBLE PRECISION XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,     B  40
     1                 HMIN,PYR                                           B  41
      DOUBLE PRECISION FCTR,TIMX,TINIT,PIES,YNS,XNS,RAT,HMX,HMY           B  50
      DOUBLE PRECISION TMSUM,ANTIM,TDEL                                   B  56
      DOUBLE PRECISION DXINV,DYINV,ARINV,PORINV                           B  57
      DOUBLE PRECISION TINT,ALPHA1,ANITP                                  B  60
      DOUBLE PRECISION TOTLQ,TOTLQI,TPIN,TPOUT                            B  62
      INTEGER OVERRD                                                      B  65
      COMMON /PRMJ/ NTIM,NPMP,NPNT,NITP,N,NX,NY,NP,NREC,INT,NNX,NNY,NUMO  B  70
     1BS,NMOV,IMOV,NPMAX,ITMAX,NZCRIT,IPRNT,NPTPND,NPNTMV,NPNTVL,NPNTD,N  B  80
     2PNCHV,NPDELC,ICHK                                                   B  91
      COMMON /PRMC/ NODEID(40,40),NPCELL(20,20),NPOLD(20,20),LIMBO(500), SB  95
     1IXOBS(5),IYOBS(5)                                                   B  96
      COMMON /HEDA/ THCK(40,40),PERM(40,40),TMWL(5,50),TMOBS(50),ANFCTR  SB 120
      COMMON /HEDB/ TMRX(40,40,2),VPRM(40,40),HI(40,40),HR(40,40),HC(40, SB 130
     140),HK(40,40),WT(40,40),REC(40,40),RECH(40,40),TIM(100),AOPT(20),  SB 140
     2TITLE(10),XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,HMIN,PYR    B 150
      COMMON /HEDC/ MX,MY,MMX,MMY,NMX,NMY,MCHK                           SB 155
      COMMON /HEDD/ TINIT,TIMX                                            B 157
      COMMON /CHMA/ PART(3,6400),CONC(20,20),TMCN(5,50),VX(40,40),VY(40, SB 161
     140),CONINT(20,20),CNRECH(20,20),POROS,SUMTCH,BETA,TIMV,STORM,STORM SB 170
     2I,CMSIN,CMSOUT,FLMIN,FLMOT,SUMIO,CELDIS,DLTRAT,CSTORM               B 180
      COMMON /CHMR/ RF,DK,RHOB,THALF,DECAY,ADSORB,SORBI,DMASS1,CSTM2      B 182R
      COMMON /BALM/ TOTLQ,TOTLQI,TPIN,TPOUT                               B 186
      COMMON /XINV/ DXINV,DYINV,ARINV,PORINV                              B 200
      COMMON /CHMC/ SUMC(20,20),VXBDY(40,40),VYBDY(40,40),PERM2(40,40)   SB 211
C     ***************************************************************     B 220
      IF (INT.GT.1) GO TO 10                                              B 230
      WRITE (6,750)                                                       B 240
      READ (5,720) TITLE                                                  B 250
      WRITE (6,730) TITLE                                                 B 260
C     ***************************************************************     B 270
C     ---INITIALIZE TEST AND CONTROL VARIABLES---                         B 280
      STORMI=0.0                                                          B 290
      SORBI=0                                                             B 295R
      RHOB=0.0                                                            B 296R
      DK=0.0                                                              B 297R
      TEST=0.0                                                            B 300
      TOTLQ=0.0                                                           B 310
      TOTLQI=0.0                                                          B 315
      TPIN=0.0                                                            B 317
      TPOUT=0.0                                                           B 318
      SUMT=0.0                                                            B 320
      SUMTCH=0.0                                                          B 330
      INT=0                                                               B 340
      IPRNT=0                                                             B 350
      MX=1                                                               SB 353
      MY=1                                                               SB 354
      MCHK=0                                                             SB 355
      NCA2=0                                                             SB 356
      NCA=0                                                               B 360
      N=0                                                                 B 370
      IMOV=0                                                              B 380
      NMOV=0                                                              B 390
      ICHK=0                                                              B 395
      DMASS1=0.0                                                          B 396R
      RF=1.0                                                              B 397R
      DECAY=0.0                                                           B 398R
      THALF=0.0                                                           B 399R
C     ***************************************************************     B 400
C     ---LOAD CONTROL PARAMETERS---                                       B 410
      READ (5,740) NTIM,NPMP,NX,NY,NPMAX,NPNT,NITP,NUMOBS,ITMAX,NREC,NPT  B 420
     1PND,NCODES,NPNTMV,NPNTVL,NPNTD,NPDELC,NPNCHV,NREACT                 B 431R
      NMX=NX                                                             SB 432
      NMY=NY                                                             SB 433
      MMX=NX-1                                                           SB 434
      MMY=NY-1                                                           SB 435
C        ---READ UPPER LEFT AND LOWER RIGHT NODAL COORDS. OF             SB 436
C           TRANSPORT SUBGRID, IN FREE FORMAT, IF NX.LT.0---             SB 437
      IF (NX.GT.0) GO TO 5                                               SB 438A
      NX=-NX                                                             SB 438B
      MCHK=1                                                             SB 438C
      READ (5,*) MX,MY,MMX,MMY                                           SB 438D
      NMX=MMX-MX+1                                                       SB 439A
      NMY=MMY-MY+1                                                       SB 439B
    5 CONTINUE                                                           SB 439H
      READ (5,800) PINT,TOL,POROS,BETA,S,TIMX,TINIT,XDEL,YDEL,DLTRAT,CEL  B 440
     1DIS,ANFCTR                                                          B 450
C        ---READ REACTION TERMS IN FREE FORMAT---                         B 454R
      IF (NREACT.EQ.1) READ (5,*) DK,RHOB,THALF                           B 455R
      NNX=NX-1                                                            B 470
      NNY=NY-1                                                            B 480
      NP=NPMAX                                                            B 490
      DXINV=1.0/XDEL                                                      B 500
      DYINV=1.0/YDEL                                                      B 510
      ARINV=DXINV*DYINV                                                   B 520
      PORINV=1.0/POROS                                                    B 530
      RF=1.0+(DK*RHOB/POROS)                                              B 532R
      IF (THALF.GT.0.0) DECAY=ALOG(2.0)/THALF                             B 534R
C     ---PRINT CONTROL PARAMETERS---                                      B 540
      WRITE (6,755)                                                      SB 545
      WRITE (6,760)                                                       B 550
      WRITE (6,770) NX,NY,XDEL,YDEL                                       B 560
      IF (MCHK.GT.0) WRITE (6,775) NMX,NMY,MX,MY,MMX,MMY                 SB 565
      WRITE (6,780) NTIM,NPMP,PINT,TIMX,TINIT                             B 570
      WRITE (6,790) S,POROS,BETA,DLTRAT,ANFCTR                            B 580
      WRITE (6,870) NITP,TOL,ITMAX,CELDIS,NPMAX,NPTPND                    B 590
      IF (NPTPND.NE.4.AND.NPTPND.NE.5.AND.NPTPND.NE.8.AND.NPTPND.NE.9.    B 601
     1.AND.NPTPND.NE.16.AND.NPTPND.NE.1) WRITE(6,880)                     B 612
      IF (NPTPND.EQ.1) WRITE (6,882)                                      B 614
      IF((NX.EQ.3.OR.NY.EQ.3).AND.NPTPND.NE.1) WRITE(6,883)               B 614A
      IF (NITP.LE.0) WRITE (6,885)                                        B 615
      WRITE (6,888)                                                      SB 617
      WRITE (6,890) NPNT,NPNTMV,NPNTVL,NPNTD,NUMOBS,NREC,NCODES,NPNCHV,N  B 620
     1PDELC                                                               B 630
      WRITE (6,895) DK,RHOB,RF,THALF,DECAY                                B 635R
      GO TO 20                                                            B 650
C     ***************************************************************     B 660
C     ---READ DATA TO REVISE TIME STEPS AND STRESSES FOR SUBSEQUENT       B 670
C        PUMPING PERIODS---                                               B 680
   10 READ (5,1060) ICHK                                                  B 690
      IF (ICHK.LE.0) WRITE (6,1110) INT                                   B 695
      IF (ICHK.LE.0) GO TO 20                                             B 701
      READ (5,1070) NTIM,NPNT,NITP,ITMAX,NREC,NPNTMV,NPNTVL,NPNTD,NPDELC  B 710
     1,NPNCHV,PINT,TIMX,TINIT                                             B 720
      WRITE (6,1080) INT                                                  B 730
      WRITE (6,1090) NTIM,NPNT,NITP,ITMAX,NREC,NPNTMV,NPNTVL,NPNTD,NPDEL  B 740
     1C,NPNCHV,PINT,TIMX,TINIT                                            B 750
C     ***************************************************************     B 760
C     ---LIST TIME INCREMENTS---                                          B 770
   20 DO 30 J=1,100                                                       B 780
      TIM(J)=0.0                                                          B 790
   30 CONTINUE                                                            B 800
      PYR=PINT*86400.0*365.25                                             B 805
      TIM(1)=TINIT                                                        B 810
      IF (NPNTMV.EQ.0) NPNTMV=999                                         B 815
      IF (S.EQ.0.0) GO TO 50                                              B 820
      DO 40 K=2,NTIM                                                      B 830
   40 TIM(K)=TIMX*TIM(K-1)                                                B 840
      WRITE (6,470)                                                       B 850
      WRITE (6,490) (TIM(K),K=1,NTIM)                                     B 861
      IF (TINIT.GT.PYR) WRITE (6,475)                                     B 865
      GO TO 60                                                            B 870
   50 ANTIM=NTIM                                                          B 882
      DO 55 K=1,NTIM                                                      B 884
   55 TIM(K)=PYR/ANTIM                                                    B 886
      WRITE (6,480) TIM(1)                                                B 890
C     ***************************************************************     B 900
C     ---INITIALIZE MATRICES---                                           B 910
   60 IF (INT.GT.1) GO TO 100                                             B 920
      DO 70 IY=1,NY                                                       B 930
      DO 70 IX=1,NX                                                       B 940
      VPRM(IX,IY)=0.0                                                     B 950
      PERM(IX,IY)=0.0                                                     B 960
      PERM2(IX,IY)=0.0                                                    B 965
      THCK(IX,IY)=0.0                                                     B 970
      RECH(IX,IY)=0.0                                                     B 980
      REC(IX,IY)=0.0                                                      B1000
      NODEID(IX,IY)=0                                                     B1010
      TMRX(IX,IY,1)=0.0                                                   B1020
      TMRX(IX,IY,2)=0.0                                                   B1030
      HI(IX,IY)=0.0                                                       B1040
      HR(IX,IY)=0.0                                                       B1050
      HC(IX,IY)=0.0                                                       B1060
      HK(IX,IY)=0.0                                                       B1070
      WT(IX,IY)=0.0                                                       B1080
      VX(IX,IY)=0.0                                                       B1090
      VY(IX,IY)=0.0                                                       B1100
      VXBDY(IX,IY)=0.0                                                    B1110
      VYBDY(IX,IY)=0.0                                                    B1120
   70 CONTINUE                                                            B1160
      DO 75 IY=1,NMY                                                     SB1161
      DO 75 IX=1,NMX                                                     SB1162
      CNRECH(IX,IY)=0.0                                                  SB1163
      CONC(IX,IY)=0.0                                                    SB1164
      CONINT(IX,IY)=0.0                                                  SB1165
      SUMC(IX,IY)=0.0                                                    SB1166
      NPCELL(IX,IY)=0                                                    SB1167
   75 CONTINUE                                                           SB1168
C     ****************************************************************    B1170
C     ---READ OBSERVATION WELL LOCATIONS---                               B1180
      IF (NUMOBS.LE.0) GO TO 100                                          B1190
      WRITE (6,755)                                                       B1195
      WRITE (6,900)                                                       B1200
      DO 80 J=1,NUMOBS                                                    B1210
      READ (5,700) IX,IY                                                  B1220
      WRITE (6,810) J,IX,IY                                               B1230
      IXOBS(J)=IX                                                         B1240
   80 IYOBS(J)=IY                                                         B1250
      DO 90 I=1,NUMOBS                                                    B1260
      DO 90 J=1,50                                                        B1270
      TMWL(I,J)=0.0                                                       B1280
   90 TMCN(I,J)=0.0                                                       B1290
C     ****************************************************************    B1300
C     ---READ PUMPAGE DATA -- (X-Y COORDINATES AND RATE IN CFS)---        B1310
C        ---SIGNS : WITHDRAWAL = POS.;  INJECTION = NEG.---               B1320
C        ---IF INJ. WELL, ALSO READ CONCENTRATION OF INJECTED WATER---    B1330
  100 IF (NREC.LE.0) GO TO 120                                            B1340
      IF (INT.GT.1.AND.ICHK.LE.0) RETURN                                  B1345
      WRITE (6,755)                                                       B1347
      WRITE (6,910)                                                       B1350
      DO 110 I=1,NREC                                                     B1360
      READ (5,710) IX,IY,FCTR,CNREC                                       B1370
      JX=IX-MX+1                                                         SB1372
      JY=IY-MY+1                                                         SB1373
      IF (JX.LT.1.OR.JY.LT.1.OR.JX.GT.NMX.OR.JY.GT.NMY) GO TO 105        SB1374
      IF (FCTR.LT.0.0) CNRECH(JX,JY)=CNREC                               SB1382
  105 REC(IX,IY)=FCTR                                                    SB1392
  110 WRITE (6,820) IX,IY,REC(IX,IY),CNREC                               SB1402
C     ****************************************************************    B1410
  120 IF (INT.GT.1) RETURN                                                B1420
      AREA=XDEL*YDEL                                                      B1430
      WRITE (6,755)                                                       B1435
      WRITE (6,690) AREA                                                  B1440
      WRITE (6,600)                                                       B1450
      WRITE (6,610) XDEL                                                  B1460
      WRITE (6,610) YDEL                                                  B1470
C     ****************************************************************    B1480
C     ---READ TRANSMISSIVITY IN FT**2/SEC INTO VPRM ARRAY---              B1490
C        ---FCTR = TRANSMISSIVITY MULTIPLIER  --->  FT**2/SEC---          B1500
      WRITE (6,530)                                                       B1510
      WRITE (6,755)                                                       B1515
      READ (5,550) INPUT,FCTR                                             B1520
      DO 160 IY=1,NY                                                      B1530
      IF (INPUT.EQ.1) READ (5,560) (VPRM(IX,IY),IX=1,NX)                  B1540
      DO 150 IX=1,NX                                                      B1550
      IF (INPUT.NE.1) GO TO 130                                           B1560
      VPRM(IX,IY)=VPRM(IX,IY)*FCTR                                        B1570
      GO TO 140                                                           B1580
  130 VPRM(IX,IY)=FCTR                                                    B1590
  140 IF (IX.EQ.1.OR.IX.EQ.NX) VPRM(IX,IY)=0.0                            B1600
      IF (IY.EQ.1.OR.IY.EQ.NY) VPRM(IX,IY)=0.0                            B1610
  150 CONTINUE                                                            B1620
  160 WRITE (6,840) (VPRM(IX,IY),IX=1,NX)                                 B1631
C     ****************************************************************    B1640
C     ---SET UP COEFFICIENT MATRIX --- BLOCK-CENTERED GRID---             B1650
C     ---AVERAGE TRANSMISSIVITY --- HARMONIC MEAN---                      B1660
      IF (ANFCTR.NE.0.0) GO TO 170                                        B1670
      WRITE (6,1050)                                                      B1680
      ANFCTR=1.0                                                          B1690
  170 PIES=3.1415927*3.1415927/2.0                                        B1700
      YNS=NY*NY                                                           B1710
      XNS=NX*NX                                                           B1720
      HMIN=2.0                                                            B1730
      DO 180 IY=2,NNY                                                     B1740
      DO 180 IX=2,NNX                                                     B1750
      IF (VPRM(IX,IY).EQ.0.0) GO TO 180                                   B1760
      TMRX(IX,IY,1)=2.0*VPRM(IX,IY)*VPRM(IX+1,IY)/(VPRM(IX,IY)*XDEL+VPRM  B1770
     1(IX+1,IY)*XDEL)                                                     B1780
      TMRX(IX,IY,2)=2.0*VPRM(IX,IY)*VPRM(IX,IY+1)/(VPRM(IX,IY)*YDEL+VPRM  B1790
     1(IX,IY+1)*YDEL)                                                     B1800
C        ---ADJUST COEFFICIENT FOR ANISOTROPY---                          B1810
      TMRX(IX,IY,2)=TMRX(IX,IY,2)*ANFCTR                                  B1820
C        ---COMPUTE MINIMUM ITERATION PARAMETER (HMIN)---                 B1830
      IF (TMRX(IX,IY,1).EQ.0.0) GO TO 180                                 B1840
      IF (TMRX(IX,IY,2).EQ.0.0) GO TO 180                                 B1850
      RAT=TMRX(IX,IY,1)*YDEL/(TMRX(IX,IY,2)*XDEL)                         B1860
      HMX=PIES/(XNS*(1.0+RAT))                                            B1870
      HMY=PIES/(YNS*(1.0+(1.0/RAT)))                                      B1880
      IF (HMX.LT.HMIN) HMIN=HMX                                           B1890
      IF (HMY.LT.HMIN) HMIN=HMY                                           B1900
  180 CONTINUE                                                            B1910
      IF (NX.GT.3.AND.NY.GT.3) GO TO 185                                  B1912
      HMX=PIES/XNS                                                        B1913
      HMY=PIES/YNS                                                        B1914
      HMIN=DMIN1(HMIN,HMX,HMY)                                            B1915
  185 CONTINUE                                                            B1916
C     ****************************************************************    B1920
C     ---READ AQUIFER THICKNESS---                                        B1930
      WRITE (6,510)                                                       B1940
      WRITE (6,755)                                                       B1945
      READ (5,550) INPUT,FCTR                                             B1950
      DO 210 IY=1,NY                                                      B1960
      IF (INPUT.EQ.1) READ (5,540) (THCK(IX,IY),IX=1,NX)                  B1970
      DO 200 IX=1,NX                                                      B1980
      IF (INPUT.NE.1) GO TO 190                                           B1990
      THCK(IX,IY)=THCK(IX,IY)*FCTR                                        B2000
      GO TO 200                                                           B2010
  190 IF (VPRM(IX,IY).NE.0.0) THCK(IX,IY)=FCTR                            B2020
  200 CONTINUE                                                            B2030
  210 WRITE (6,500) (THCK(IX,IY),IX=1,NX)                                 B2040
C     ****************************************************************    B2050
C     ---READ DIFFUSE RECHARGE AND DISCHARGE---                           B2060
      WRITE (6,830)                                                       B2070
      WRITE (6,755)                                                       B2075
      READ (5,550) INPUT,FCTR                                             B2080
      DO 240 IY=1,NY                                                      B2090
      IF (INPUT.EQ.1) READ (5,560) (RECH(IX,IY),IX=1,NX)                  B2100
      DO 230 IX=1,NX                                                      B2110
      IF (INPUT.NE.1) GO TO 220                                           B2120
      RECH(IX,IY)=RECH(IX,IY)*FCTR                                        B2130
      GO TO 230                                                           B2140
  220 IF (THCK(IX,IY).NE.0.0) RECH(IX,IY)=FCTR                            B2150
  230 CONTINUE                                                            B2160
  240 WRITE (6,840) (RECH(IX,IY),IX=1,NX)                                 B2170
C     ****************************************************************    B2180
C     ---COMPUTE PERMEABILITY FROM TRANSMISSIVITY---                      B2190
C     ---COUNT NO. OF CELLS IN AQUIFER---                                 B2200
C     ---SET NZCRIT = 2% OF THE NO. OF CELLS IN THE AQUIFER---            B2210
      DO 250 IX=1,NX                                                      B2220
      DO 250 IY=1,NY                                                      B2230
      IF (THCK(IX,IY).EQ.0.0) GO TO 250                                   B2240
      PERM(IX,IY)=VPRM(IX,IY)/THCK(IX,IY)                                 B2250
      NCA=NCA+1                                                           B2260
      IF (MCHK.EQ.0) GO TO 250                                           SB2262
      IF (IX.LT.MX.OR.IX.GT.MMX) GO TO 250                               SB2263
      IF (IY.LT.MY.OR.IY.GT.MMY) GO TO 250                               SB2264
      NCA2=NCA2+1                                                        SB2265
  250 VPRM(IX,IY)=0.0                                                     B2270
C                                                                         B2280
      AAQ=NCA*AREA                                                        B2290
      NZCRIT=(NCA+25)/50                                                  B2300
      IF (MCHK.GT.0) NZCRIT=(NCA2+25)/50                                 SB2301
      IF (NZCRIT.EQ.0) NZCRIT=1                                          SB2305
      WRITE (6,620)                                                       B2310
      WRITE (6,755)                                                       B2315
      DO 260 IY=1,NY                                                      B2320
  260 WRITE (6,840) (PERM(IX,IY),IX=1,NX)                                 B2331
      IF (MCHK.EQ.0) GO TO 265                                           SB2332
      AAQ2=NCA2*AREA                                                     SB2333
      WRITE (6,633)                                                      SB2334
      WRITE (6,635) NCA,AAQ                                              SB2335
      WRITE (6,634)                                                      SB2336
      WRITE (6,630) NCA2,AAQ2,NZCRIT                                     SB2337
      GO TO 267                                                          SB2338
  265 CONTINUE                                                           SB2339
      WRITE (6,630) NCA,AAQ,NZCRIT                                        B2340
  267 CONTINUE                                                           SB2345
C     ****************************************************************    B2350
C     ---READ NODE IDENTIFICATION CARDS---                                B2360
C        ---SET VERT. PERM., SOURCE CONC., AND DIFFUSE RECHARGE---        B2370
C           ---SPECIFY CODES TO FIT YOUR NEEDS---                         B2380
      WRITE (6,570)                                                       B2390
      READ (5,550) INPUT,FCTR                                             B2400
      DO 280 IY=1,NY                                                      B2410
      IF (INPUT.EQ.1) READ (5,640) (NODEID(IX,IY),IX=1,NX)                B2420
      DO 270 IX=1,NX                                                      B2430
  270 IF (INPUT.NE.1.AND.THCK(IX,IY).NE.0.0) NODEID(IX,IY)=FCTR           B2440
  280 WRITE (6,580) (NODEID(IX,IY),IX=1,NX)                               B2450
      WRITE (6,920) NCODES                                                B2460
      IF (NCODES.LE.0) GO TO 310                                          B2470
      WRITE (6,930)                                                       B2480
      DO 300 IJ=1,NCODES                                                  B2490
      READ (5,850) ICODE,FCTR1,FCTR2,FCTR3,OVERRD                         B2500
      DO 290 IX=1,NX                                                      B2510
      DO 290 IY=1,NY                                                      B2520
      IF (NODEID(IX,IY).NE.ICODE) GO TO 290                               B2530
      VPRM(IX,IY)=FCTR1                                                   B2540
      JX=IX-MX+1                                                         SB2544
      JY=IY-MY+1                                                         SB2545
      IF (JX.LT.1.OR.JY.LT.1.OR.JX.GT.NMX.OR.JY.GT.NMY) GO TO 285        SB2546
      CNRECH(JX,JY)=FCTR2                                                SB2551
  285 CONTINUE                                                           SB2555
      IF (OVERRD.NE.0) RECH(IX,IY)=FCTR3                                  B2560
  290 CONTINUE                                                            B2570
      WRITE (6,860) ICODE,FCTR1,FCTR2                                     B2580
  300 IF (OVERRD.NE.0) WRITE (6,1100) FCTR3                               B2590
  310 WRITE (6,590)                                                       B2600
      WRITE (6,755)                                                       B2615
      DO 320 IY=1,NY                                                      B2610
  320 WRITE (6,840) (VPRM(IX,IY),IX=1,NX)                                 B2621
C     ***************************************************************     B2630
C     ---READ WATER-TABLE ELEVATION---                                    B2640
      WRITE (6,670)                                                       B2650
      WRITE (6,755)                                                       B2655
      READ (5,550) INPUT,FCTR                                             B2660
      DO 350 IY=1,NY                                                      B2670
      IF (INPUT.EQ.1) READ (5,660) (WT(IX,IY),IX=1,NX)                    B2680
      DO 340 IX=1,NX                                                      B2690
      IF (INPUT.NE.1) GO TO 330                                           B2700
      WT(IX,IY)=WT(IX,IY)*FCTR                                            B2710
      GO TO 340                                                           B2720
  330 IF (THCK(IX,IY).NE.0.0) WT(IX,IY)=FCTR                              B2730
  340 CONTINUE                                                            B2740
  350 WRITE (6,680) (WT(IX,IY),IX=1,NX)                                   B2750
C     ****************************************************************    B2760
C     ---SET INITIAL HEADS---                                             B2770
      DO 360 IX=1,NX                                                      B2780
      DO 360 IY=1,NY                                                      B2790
      HI(IX,IY)=WT(IX,IY)                                                 B2800
      HC(IX,IY)=HI(IX,IY)                                                 B2810
      HR(IX,IY)=HI(IX,IY)                                                 B2820
  360 HK(IX,IY)=HI(IX,IY)                                                 B2830
C                                                                         B2840
      CALL OUTPT                                                          B2850
C     ***************************************************************     B2860
C     ---COMPUTE ITERATION PARAMETERS---                                  B2870
      DO 370 ID=1,20                                                      B2880
      AOPT(ID)=0.0                                                        B2890
  370 CONTINUE                                                            B2900
      ANITP=NITP-1                                                        B2910
      ALPHA1=DEXP(DLOG(1.0/HMIN)/ANITP)                                   B2920
      AOPT(1)=HMIN                                                        B2930
      DO 380 IP=2,NITP                                                    B2940
  380 AOPT(IP)=AOPT(IP-1)*ALPHA1                                          B2950
C                                                                         B2960
      WRITE (6,450)                                                       B2970
      WRITE (6,460) (AOPT(IP),IP=1,NITP)                                  B2981
C     ****************************************************************    B2990
C     ---READ INITIAL CONCENTRATIONS AND COMPUTE INITIAL MASS STORED---   B3000
C        ---IF USING SMALLER SUBGRID FOR TRANSPORT, READ INITIAL         SB3005
C           CONCENTRATION ARRAY FOR SUBGRID NODES ONLY---                SB3006
      READ (5,550) INPUT,FCTR                                             B3010
      CFCTR2=DK*RHOB                                                      B3015R
      DO 420 IY=1,NMY                                                    SB3021
      JY=IY+MY-1                                                         SB3025
      IF (INPUT.EQ.1) READ (5,660) (CONC(IX,IY),IX=1,NMX)                SB3031
      DO 410 IX=1,NMX                                                    SB3041
      JX=IX+MX-1                                                         SB3045
      IF (INPUT.NE.1) GO TO 390                                           B3050
      CONC(IX,IY)=CONC(IX,IY)*FCTR                                        B3060
      GO TO 400                                                           B3070
  390 IF (THCK(JX,JY).NE.0.0) CONC(IX,IY)=FCTR                           SB3081
  400 CONINT(IX,IY)=CONC(IX,IY)                                           B3090
      CFCTR=CONINT(IX,IY)*THCK(JX,JY)*AREA                               SB3093R
      SORBI=SORBI+CFCTR*CFCTR2                                            B3094R
  410 STORMI=STORMI+CFCTR*POROS                                           B3096R
  420 CONTINUE                                                            B3110
C     ****************************************************************    B3120
C     ---CHECK DATA SETS FOR INTERNAL CONSISTENCY---                      B3130
      DO 440 IX=1,NX                                                      B3140
      DO 440 IY=1,NY                                                      B3150
      IF (THCK(IX,IY).GT.0.0) THEN                                        B3161
       IF (IX.EQ.1.OR.IY.EQ.1.OR.IX.EQ.NX.OR.IY.EQ.NY)                    B3162
     1     WRITE (6,935) IX,IY                                            B3163
       GO TO 430                                                          B3164
      END IF                                                              B3165
      IF (TMRX(IX,IY,1).GT.0.0) WRITE (6,940) IX,IY                       B3170
      IF (TMRX(IX,IY,2).GT.0.0) WRITE (6,950) IX,IY                       B3180
      IF (IX.NE.1) THEN                                                   B3182
       IF (TMRX(IX-1,IY,1).GT.0.0) WRITE (6,940) IX,IY                    B3183
      END IF                                                              B3184
      IF (IY.NE.1) THEN                                                   B3185
       IF (TMRX(IX,IY-1,2).GT.0.0) WRITE (6,950) IX,IY                    B3186
      END IF                                                              B3187
      IF (NODEID(IX,IY).GT.0) WRITE (6,960) IX,IY                         B3190
      IF (WT(IX,IY).NE.0.0) WRITE (6,970) IX,IY                           B3200
      IF (RECH(IX,IY).NE.0.0) WRITE (6,980) IX,IY                         B3210
      IF (REC(IX,IY).NE.0.0) WRITE (6,990) IX,IY                          B3220
  430 IF (PERM(IX,IY).GT.0.0) GO TO 440                                   B3230
      IF (NODEID(IX,IY).GT.0.0) WRITE (6,1000) IX,IY                      B3240
      IF (WT(IX,IY).NE.0.0) WRITE (6,1010) IX,IY                          B3250
      IF (RECH(IX,IY).NE.0.0) WRITE (6,1020) IX,IY                        B3260
      IF (REC(IX,IY).NE.0.0) WRITE (6,1030) IX,IY                         B3270
      IF (THCK(IX,IY).GT.0.0) WRITE (6,1040) IX,IY                        B3280
  440 CONTINUE                                                            B3290
C     COMPUTE HARMONIC MEAN PERMEABILITY                                  B3291
      DO 445 IX=2,NNX                                                     B3292
      DO 445 IY=2,NNY                                                     B3293
      IF (THCK(IX,IY).EQ.0.0) GO TO 445                                   B3294
      PNODE=PERM(IX,IY)                                                   B3295
      PERM(IX,IY)=2.0*PNODE*PERM(IX+1,IY)/(PNODE+PERM(IX+1,IY))           B3296
      PERM2(IX,IY)=2.0*PNODE*PERM(IX,IY+1)/(PNODE+PERM(IX,IY+1))          B3297
  445 CONTINUE                                                            B3298
C     ****************************************************************    B3300
      RETURN                                                              B3310
C     ****************************************************************    B3320
C                                                                         B3330
C                                                                         B3340
C                                                                         B3350
  450 FORMAT (1H1,20HITERATION PARAMETERS)                                B3360
  460 FORMAT (3H   ,1G20.6)                                               B3370
  470 FORMAT (1H1,27HTIME INTERVALS (IN SECONDS))                         B3380
  475 FORMAT (1H0,5X,65H*** WARNING ***   INITIAL TIME STEP IS LONGER TH  B3384
     1AN PUMPING PERIOD/25X,34H***ADJUST EITHER TINIT OR PINT.***/)       B3385
  480 FORMAT (1H1,15X,17HSTEADY-STATE FLOW//5X,57HTIME INTERVAL (IN SEC)  B3390
     1 FOR SOLUTE-TRANSPORT SIMULATION = ,G12.5)                          B3400
  490 FORMAT (3H   ,10G12.5)                                              B3410
  500 FORMAT (3H   ,25F5.1)                                              SB3421
  510 FORMAT (1H1,22HAQUIFER THICKNESS (FT))                              B3430
  530 FORMAT (1H1,30HTRANSMISSIVITY MAP (FT*FT/SEC))                      B3450
  540 FORMAT (20G3.0)                                                     B3460
  550 FORMAT (I1,G10.0)                                                   B3470
  560 FORMAT (20G4.1)                                                     B3480
  570 FORMAT (1H1,23HNODE IDENTIFICATION MAP//)                           B3490
  580 FORMAT (1H ,40I3)                                                  SB3501
  590 FORMAT (1H1,45HVERTICAL PERMEABILITY/THICKNESS (FT/(FT*SEC)))       B3510
  600 FORMAT (1H0,10X,12HX-Y SPACING:)                                    B3520
  610 FORMAT (1H ,12X,10G12.5)                                            B3530
  620 FORMAT (1H1,24HPERMEABILTY MAP (FT/SEC))                            B3540
  630 FORMAT (1H0,05X,05X,44HNO. OF FINITE-DIFFERENCE CELLS IN AQUIFER = SB3551
     1 ,I4//10X,28HAREA OF AQUIFER IN MODEL  = ,G12.5,10H   SQ. FT.////1  B3560
     20X,47HNZCRIT   (MAX. NO. OF CELLS THAT CAN BE VOID OF/20X,56HPARTI  B3570
     3CLES;  IF EXCEEDED, PARTICLES ARE REGENERATED)   = ,I4/)            B3580
  633 FORMAT (1H0,////2X,26HFLOW MODEL (PRIMARY GRID):/)                 SB3582
  634 FORMAT (1H0,///2X,18HTRANSPORT SUBGRID:/)                          SB3583
  635 FORMAT (1H0,05X,05X,44HNO. OF FINITE-DIFFERENCE CELLS IN AQUIFER = SB3584
     1 ,I4//10X,28HAREA OF AQUIFER IN MODEL  = ,G12.5,10H   SQ. FT./)    SB3585
  640 FORMAT (20I1)                                                       B3590
  660 FORMAT (20G4.0)                                                     B3610
  670 FORMAT (1H1,11HWATER TABLE)                                         B3620
  680 FORMAT (1H ,20F5.0)                                                 B3630
  690 FORMAT (1H0,10X,19HAREA OF ONE CELL = ,G12.4)                       B3640
  700 FORMAT (2I2)                                                        B3650
  710 FORMAT (2I2,2G8.2)                                                  B3660
  720 FORMAT (10A8)                                                       B3670
  730 FORMAT (1H0,10A8)                                                   B3680
  740 FORMAT (18I4)                                                       B3691R
  750 FORMAT (1H1,77HU.S.G.S. METHOD-OF-CHARACTERISTICS MODEL FOR SOLUTE  B3700
     1 TRANSPORT IN GROUND WATER)                                         B3710
  755 FORMAT (1H )                                                       SB3715
  760 FORMAT (1H0,21X,21HI N P U T     D A T A)                           B3720
  770 FORMAT (1H0,23X,16HGRID DESCRIPTORS//13X,30HNX    (NUMBER OF COLUM  B3730
     1NS)  =  ,I4/13X,28HNY    (NUMBER OF ROWS)     =,I6/13X,29HXDEL  (X  B3740
     2-DISTANCE IN FEET) = ,F7.1/13X,29HYDEL  (Y-DISTANCE IN FEET) = ,F7  B3750
     3.1)                                                                 B3760
  775 FORMAT (1H0,18X,31HSECONDARY SUBGRID FOR TRANSPORT//16X,30HNMX   ( SB3763
     1NUMBER OF COLUMNS)  =  ,I4/16X,30HNMY   (NUMBER OF ROWS)     =  ,I SB3764
     24//16X,38HCROSS-REF. TO PRIMARY GRID      IX  IY/46X,8H --- ---/18 SB3765
     3X,28HFIRST NODE (UPPER LEFT) AT: ,2I4/18X,28HLAST NODE (LOWER RIGH SB3766
     4T) AT: ,2I4)                                                       SB3767
  780 FORMAT (1H0,23X,16HTIME  PARAMETERS//13X,40HNTIM   (MAX. NO. OF TI  B3770
     1ME STEPS)       = ,I6/13X,40HNPMP   (NO. OF PUMPING PERIODS)        B3780
     2 = ,I6/13X,39HPINT   (PUMPING PERIOD IN YEARS)      =,F11.3/13X,39  B3791
     3HTIMX   (TIME INCREMENT MULTIPLIER)    =,F10.2/13X,39HTINIT  (INIT  B3800
     4IAL TIME STEP IN SEC.)    =,F8.0)                                   B3810
  790 FORMAT (1H0,14X,34HHYDROLOGIC AND CHEMICAL PARAMETERS//13X,1HS,7X,  B3820
     129H(STORAGE COEFFICIENT)       =,5X,F9.6/13X,28HPOROS   (EFFECTIVE  B3830
     2 POROSITY),8X,3H=  ,F9.3/13X,39HBETA    (LONGITUDINAL DISPERSIVITY SB3841
     3) =  ,F7.1/13X,31HDLTRAT  (RATIO OF TRANSVERSE TO/21X,30HLONGITUDI SB3851
     4NAL DISPERSIVITY)  = ,F9.2/13X,39HANFCTR  (RATIO OF T-YY TO T-XX)   B3860
     5    =  ,F12.6)                                                      B3870
  800 FORMAT (12G5.0)                                                     B3880
  810 FORMAT (1H ,16X,I2,5X,I2,4X,I2)                                     B3890
  820 FORMAT (1H ,7X,2I4,3X,F9.4,3X,F8.2)                                 B3895
  830 FORMAT (1H1,39HDIFFUSE RECHARGE AND DISCHARGE (FT/SEC))             B3910
  840 FORMAT (1H ,1P10E10.2)                                              B3920
  850 FORMAT (I2,3G10.2,I2)                                               B3930
  860 FORMAT (1H0,7X,I2,7X,E10.3,5X,F9.2)                                 B3940
  870 FORMAT (1H0,21X,20HEXECUTION PARAMETERS//13X,39HNITP   (NO. OF ITE  B3950
     1RATION PARAMETERS) = ,I4/13X,39HTOL    (CONVERGENCE CRITERIA - ADI  B3960
     2P) = ,E9.2/13X,39HITMAX  (MAX.NO.OF ITERATIONS - ADIP) = ,I4/13X,3  B3971
     34HCELDIS (MAX.CELL DISTANCE PER MOVE/24X,28HOF PARTICLES - M.O.C.)  B3980
     4    = ,F8.3/13X,30HNPMAX  (MAX. NO. OF PARTICLES),7X,2H= ,I4/12X,3  B3990
     52H NPTPND (NO. PARTICLES PER NODE),6X,3H=  ,I4)                     B4000
  880 FORMAT (1H0,5X,47H*** WARNING *** NPTPND MUST =  1,4,5,8,9, OR 16)  B4011A
  882 FORMAT (1H0,5X,58H*** CAUTION *** USE OF NPTPND=1 MAY CAUSE LOSS O  B4012A
     1F ACCURACY)                                                         B4012B
  883 FORMAT (1H0,5X,56H*** ONE-DIMENSIONAL *** WILL USE ONLY 1 ROW OF P  B4013A
     1ARTICLES/13X,35HUSE 2 PARTICLES FOR NPTPND = 4 OR 5/13X,35HUSE 3 P  B4013B
     2ARTICLES FOR NPTPND = 8 OR 9/13X,31HUSE 4 PARTICLES FOR NPTPND = 1  B4013C
     36)                                                                  B4013D
  885 FORMAT (1H0,5X,38H*** WARNING ***  NITP MUST BE POSITIVE)           B4015
  888 FORMAT (1H1)                                                       SB4017
  890 FORMAT (1H0,23X,15HPROGRAM OPTIONS//13X,30HNPNT   (TIME STEP INTER  B4020
     1VAL FOR/21X,18HCOMPLETE PRINTOUT),7X,3H=  ,I4/13X,31HNPNTMV (MOVE   B4030
     2INTERVAL FOR CHEM./21X,28HCONCENTRATION PRINTOUT)  =  ,I4/13X,29HN  B4040
     3PNTVL (PRINT OPTION-VELOCITY/21X,24H0=NO; 1=FIRST TIME STEP;/21X,1  B4050
     47H2=ALL TIME STEPS),8X,3H=  ,I4/13X,31HNPNTD  (PRINT OPTION-DISP.C  B4060
     5OEF./21X,24H0=NO; 1=FIRST TIME STEP;/21X,17H2=ALL TIME STEPS),8X,3  B4070
     6H=  ,I4/13X,32HNUMOBS (NO. OF OBSERVATION WELLS/21X,28HFOR HYDROGR  B4080
     7APH PRINTOUT) =  ,I4/13X,35HNREC   (NO. OF PUMPING WELLS)    = ,I5  B4090
     8/13X,24HNCODES (FOR NODE IDENT.),9X,2H= ,I5/13X,25HNPNCHV (PUNCH V  B4100
     9ELOCITIES),8X,2H= ,I5/13X,36HNPDELC (PRINT OPT.-CONC. CHANGE) =  ,  B4110
     $I4)                                                                 B4120
  895 FORMAT (1H0,23X,14HREACTION TERMS//13X,37HDK      (DISTRIBUTION CO  B4121R
     1EFFICIENT) = ,E12.5/13X,37HRHOB    (BULK DENSITY OF SOLIDS)   = ,E  B4122R
     212.5/13X,37HRF      (RETARDATION FACTOR)       = ,E12.5/13X,37HTHA  B4123R
     3LF   (HALF LIFE OF DECAY,IN SEC)= ,E12.5/13X,37HDECAY   (DECAY CON  B4124R
     4STANT=LN 2/THALF)= ,E12.5)                                          B4125R
  900 FORMAT (1H0,10X,29HLOCATION OF OBSERVATION WELLS//17X,3HNO.,5X,1HX  B4130
     1,5X,1HY/)                                                           B4140
  910 FORMAT (1H0,10X,28HLOCATION  OF  PUMPING  WELLS//11X,28HX   Y   RA  B4150
     1TE(IN CFS)   CONC./)                                                B4160
  920 FORMAT (1H0,5X,37HNO. OF NODE IDENT. CODES SPECIFIED = ,I2)         B4170
  930 FORMAT (1H0,10X,41HTHE FOLLOWING ASSIGNMENTS HAVE BEEN MADE:/5X,51  B4180
     1HCODE NO.     LEAKANCE     SOURCE CONC.     RECHARGE)               B4190
  935 FORMAT (1H ,5X,54H*** WARNING ***   THCK.NE.0.0 ON BOUNDARY AT NOD  B4195
     1E IX =,I4,6H, IY =,I4)                                              B4196
  940 FORMAT (1H ,5X,61H*** WARNING ***   THCK.EQ.0.0 AND TMRX(X).GT.0.0  B4200
     1 AT NODE IX =,I4,6H, IY =,I4)                                       B4210
  950 FORMAT (1H ,5X,61H*** WARNING ***   THCK.EQ.0.0 AND TMRX(Y).GT.0.0  B4220
     1 AT NODE IX =,I4,6H, IY =,I4)                                       B4230
  960 FORMAT (1H ,5X,61H*** WARNING ***    THCK.EQ.0.0 AND NODEID.GT.0.0  B4240
     1 AT NODE IX =,I4,6H, IY =,I4)                                       B4250
  970 FORMAT (1H ,5X,56H*** WARNING ***   THCK.EQ.0.0 AND WT.NE.0.0 AT N  B4260
     1ODE IX =,I4,6H, IY =,I4)                                            B4270
  980 FORMAT (1H ,5X,58H*** WARNING ***   THCK.EQ.0.0 AND RECH.NE.0.0 AT  B4280
     1 NODE IX =,I4,6H, IY =,I4)                                          B4290
  990 FORMAT (1H ,5X,58H*** WARNING ***    THCK.EQ.0.0 AND REC.NE.0.0 AT  B4300
     1 NODE IX =,I4,6H, IY =,I4)                                          B4310
 1000 FORMAT (1H ,5X,61H*** WARNING ***    PERM.EQ.0.0 AND NODEID.GT.0.0  B4320
     1 AT NODE IX =,I4,6H, IY =,I4)                                       B4330
 1010 FORMAT (1H ,5X,56H*** WARNING ***   PERM.EQ.0.0 AND WT.NE.0.0 AT N  B4340
     1ODE IX =,I4,6H, IY =,I4)                                            B4350
 1020 FORMAT (1H ,5X,58H*** WARNING ***   PERM.EQ.0.0 AND RECH.NE.0.0 AT  B4360
     1 NODE IX =,I4,6H, IY =,I4)                                          B4370
 1030 FORMAT (1H ,5X,58H*** WARNING ***    PERM.EQ.0.0 AND REC.NE.0.0 AT  B4380
     1 NODE IX =,I4,6H, IY =,I4)                                          B4390
 1040 FORMAT (1H ,5X,58H*** WARNING ***   PERM.EQ.0.0 AND THCK.GT.0.0 AT  B4400
     1 NODE IX =,I4,6H, IY =,I4)                                          B4410
 1050 FORMAT (1H0,5X,45H*** WARNING ***   ANFCTR WAS SPECIFIED AS 0.0/23  B4420
     1X,34HDEFAULT ACTION: RESET ANFCTR = 1.0)                            B4430
 1060 FORMAT (I1)                                                         B4440
 1070 FORMAT (10I4,3G5.0)                                                 B4450
 1080 FORMAT (1H1,5X,25HSTART PUMPING PERIOD NO. ,I2//2X,75HTHE FOLLOWIN  B4460
     1G TIME STEP, PUMPAGE, AND PRINT PARAMETERS HAVE BEEN REDEFINED:/)   B4470
 1090 FORMAT (1H0,14X,9HNTIM   = ,I4/15X,9HNPNT   = ,I4/15X,9HNITP   = ,  B4480
     1I4/15X,9HITMAX  = ,I4/15X,9HNREC   = ,I4/15X,9HNPNTMV = ,I4/15X,9H  B4490
     2NPNTVL = ,I4/15X,9HNPNTD  = ,I4/15X,9HNPDELC = ,I4/15X,9HNPNCHV =   B4500
     3,I4/15X,9HPINT   = ,F10.3/15X,9HTIMX   = ,F10.3/15X,9HTINIT  = ,F1  B4510
     40.3/)                                                               B4520
 1100 FORMAT (1H ,46X,E10.3)                                              B4530
 1110 FORMAT (1H1,5X,25HSTART PUMPING PERIOD NO. ,I2//2X,23HNO PARAMETER  B4532
     1S REDEFINED/)                                                       B4533
      END                                                                 B4540-
      SUBROUTINE ITERAT                                                   C  10
      DOUBLE PRECISION DMIN1,DEXP,DLOG,DABS                               C  20
      DOUBLE PRECISION TMRX,VPRM,HI,HR,HC,HK,WT,REC,RECH,TIM,AOPT,TITLE   C  30
      DOUBLE PRECISION XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,     C  40
     1                 HMIN,PYR                                           C  41
      DOUBLE PRECISION DXINV,DYINV,ARINV,PORINV                           C  45
      DOUBLE PRECISION B,G,W,A,C,E,F,DR,DC,TBAR,TMK,COEF,BLH,BRK,CHK,QL,  C  50
     1                 BRH                                                C  51
      DOUBLE PRECISION TOTLQ,TOTLQI,TPIN,TPOUT,PQIN,PQOUT,DELQ            C  52
      COMMON /PRMJ/ NTIM,NPMP,NPNT,NITP,N,NX,NY,NP,NREC,INT,NNX,NNY,NUMO  C  60
     1BS,NMOV,IMOV,NPMAX,ITMAX,NZCRIT,IPRNT,NPTPND,NPNTMV,NPNTVL,NPNTD,N  C  70
     2PNCHV,NPDELC,ICHK                                                   C  81
      COMMON /HEDA/ THCK(40,40),PERM(40,40),TMWL(5,50),TMOBS(50),ANFCTR  SC 110
      COMMON /HEDB/ TMRX(40,40,2),VPRM(40,40),HI(40,40),HR(40,40),HC(40, SC 120
     140),HK(40,40),WT(40,40),REC(40,40),RECH(40,40),TIM(100),AOPT(20),  SC 130
     2TITLE(10),XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,HMIN,PYR    C 140
      COMMON /BALM/ TOTLQ,TOTLQI,TPIN,TPOUT                               C 146
      COMMON /XINV/ DXINV,DYINV,ARINV,PORINV                              C 160
      DIMENSION W(40), B(40), G(40)                                      SC 171
C     ***************************************************************     C 180
      KOUNT=0                                                             C 190
      PQIN=0.0                                                            C 192
      PQOUT=0.0                                                           C 193
C     ---COMPUTE ROW AND COLUMN---                                        C 200
C        ---CALL NEW ITERATION PARAMETER---                               C 210
   10 REMN=MOD(KOUNT,NITP)                                                C 220
      IF (REMN.EQ.0) NTH=0                                                C 230
      NTH=NTH+1                                                           C 240
      PARAM=AOPT(NTH)                                                     C 250
C     ***************************************************************     C 260
C     ---ROW COMPUTATIONS---                                              C 270
      TEST=0.0                                                            C 280
      RHO=S/TIM(N)                                                        C 290
      BRK=-RHO                                                            C 300
      DO 50 IY=1,NY                                                       C 310
      DO 20 M=1,NX                                                        C 320
      W(M)=0.0                                                            C 330
      B(M)=0.0                                                            C 340
      G(M)=0.0                                                            C 350
   20 CONTINUE                                                            C 360
      DO 30 IX=1,NX                                                       C 370
      IF (THCK(IX,IY).EQ.0.0) GO TO 30                                    C 380
      COEF=VPRM(IX,IY)                                                    C 390
      QL=-COEF*WT(IX,IY)                                                  C 400
      A=TMRX(IX-1,IY,1)*DXINV                                             C 410
      C=TMRX(IX,IY,1)*DXINV                                               C 420
      E=TMRX(IX,IY-1,2)*DYINV                                             C 430
      F=TMRX(IX,IY,2)*DYINV                                               C 440
      TBAR=A+C+E+F                                                        C 450
      TMK=TBAR*PARAM                                                      C 460
      BLH=-A-C-RHO-COEF-TMK                                               C 470
      BRH=E+F-TMK                                                         C 480
      DR=BRH*HC(IX,IY)+BRK*HK(IX,IY)-E*HC(IX,IY-1)-F*HC(IX,IY+1)+QL+RECH  C 490
     1(IX,IY)+REC(IX,IY)*ARINV                                            C 500
      W(IX)=BLH-A*B(IX-1)                                                 C 510
      B(IX)=C/W(IX)                                                       C 520
      G(IX)=(DR-A*G(IX-1))/W(IX)                                          C 530
   30 CONTINUE                                                            C 540
C                                                                         C 550
C        ---BACK SUBSTITUTION---                                          C 560
      DO 40 J=2,NX                                                        C 570
      IJ=J-1                                                              C 580
      IS=NX-IJ                                                            C 590
   40 HR(IS,IY)=G(IS)-B(IS)*HR(IS+1,IY)                                   C 600
   50 CONTINUE                                                            C 610
C     ***************************************************************     C 620
C     ---COLUMN COMPUTATIONS---                                           C 630
      DO 90 IX=1,NX                                                       C 640
      DO 60 M=1,NY                                                        C 650
      W(M)=0.0                                                            C 660
      B(M)=0.0                                                            C 670
   60 G(M)=0.0                                                            C 680
      DO 70 IY=1,NY                                                       C 690
      IF (THCK(IX,IY).EQ.0.0) GO TO 70                                    C 700
      COEF=VPRM(IX,IY)                                                    C 710
      QL=-COEF*WT(IX,IY)                                                  C 720
      A=TMRX(IX,IY-1,2)*DYINV                                             C 730
      C=TMRX(IX,IY,2)*DYINV                                               C 740
      E=TMRX(IX-1,IY,1)*DXINV                                             C 750
      F=TMRX(IX,IY,1)*DXINV                                               C 760
      TBAR=A+C+E+F                                                        C 770
      TMK=TBAR*PARAM                                                      C 780
      BLH=-A-C-RHO-COEF-TMK                                               C 790
      BRH=E+F-TMK                                                         C 800
      DC=BRH*HR(IX,IY)+BRK*HK(IX,IY)-E*HR(IX-1,IY)-F*HR(IX+1,IY)+QL+RECH  C 810
     1(IX,IY)+REC(IX,IY)*ARINV                                            C 820
      W(IY)=BLH-A*B(IY-1)                                                 C 830
      B(IY)=C/W(IY)                                                       C 840
      G(IY)=(DC-A*G(IY-1))/W(IY)                                          C 850
   70 CONTINUE                                                            C 860
C                                                                         C 870
C        ---BACK SUBSTITUTION---                                          C 880
      DO 80 J=2,NY                                                        C 890
      IJ=J-1                                                              C 900
      IB=NY-IJ                                                            C 910
      HC(IX,IB)=G(IB)-B(IB)*HC(IX,IB+1)                                   C 920
      IF (THCK(IX,IB).EQ.0.0) GO TO 80                                    C 930
      CHK=DABS(HC(IX,IB)-HR(IX,IB))                                       C 940
      IF (CHK.GT.TOL) TEST=1.0                                            C 950
   80 CONTINUE                                                            C 960
   90 CONTINUE                                                            C 970
C     ***************************************************************     C 980
      KOUNT=KOUNT+1                                                       C 990
      IF (TEST.EQ.0.0) GO TO 120                                          C1000
      IF (KOUNT.GE.ITMAX) GO TO 100                                       C1010
      GO TO 10                                                            C1020
C     ****************************************************************    C1030
C     ---TERMINATE PROGRAM -- ITMAX EXCEEDED---                           C1040
  100 WRITE (6,160)                                                       C1050
      DO 110 IX=1,NX                                                      C1060
      DO 110 IY=1,NY                                                      C1070
  110 HK(IX,IY)=HC(IX,IY)                                                 C1080
      CALL OUTPT                                                          C1090
      STOP                                                                C1100
C     ****************************************************************    C1110
C     ---SET NEW HEAD (HK)---                                             C1120
  120 DO 130 IY=1,NY                                                      C1130
      DO 130 IX=1,NX                                                      C1140
      IF (THCK(IX,IY).EQ.0.0) GO TO 130                                   C1150
      HR(IX,IY)=HK(IX,IY)                                                 C1160
      HK(IX,IY)=HC(IX,IY)                                                 C1170
C                                                                         C1180
C     ---CUMULATE PUMPAGE AND RECHARGE FOR MASS BALANCE---                C1181
      IF (REC(IX,IY).GT.0.0) GO TO 32                                     C1182
      PQIN=PQIN+REC(IX,IY)                                                C1183
      GO TO 34                                                            C1184
   32 PQOUT=PQOUT+REC(IX,IY)                                              C1185
   34 IF (RECH(IX,IY).GT.0.0) GO TO 36                                    C1186
      PQIN=PQIN+RECH(IX,IY)*AREA                                          C1187
      GO TO 38                                                            C1188
   36 PQOUT=PQOUT+RECH(IX,IY)*AREA                                        C1189
C     ---COMPUTE LEAKAGE FOR MASS BALANCE---                              C1190
   38 IF (VPRM(IX,IY).EQ.0.0) GO TO 130                                   C1201
      DELQ=VPRM(IX,IY)*AREA*(WT(IX,IY)-HK(IX,IY))                         C1210
      IF (DELQ.GT.0.0) GO TO 125                                          C1215
      TOTLQ=TOTLQ+DELQ*TIM(N)                                             C1220
      GO TO 130                                                           C1222
  125 TOTLQI=TOTLQI+DELQ*TIM(N)                                           C1224
  130 CONTINUE                                                            C1230
      TPIN=PQIN*TIM(N)+TPIN                                               C1232
      TPOUT=PQOUT*TIM(N)+TPOUT                                            C1233
C                                                                         C1240
      WRITE (6,140) N                                                     C1250
      WRITE (6,150) KOUNT                                                 C1260
C     ****************************************************************    C1270
      RETURN                                                              C1280
C     ****************************************************************    C1290
C                                                                         C1300
C                                                                         C1310
C                                                                         C1320
  140 FORMAT (1H0//3X,4HN = ,1I4)                                         C1330
  150 FORMAT (1H ,2X,23HNUMBER OF ITERATIONS = ,1I4)                      C1340
  160 FORMAT (1H0,5X,64H***   EXECUTION TERMINATED -- MAX. NO. ITERATION  C1350
     1S EXCEEDED   ***/26X,21HFINAL OUTPUT FOLLOWS:)                      C1360
      END                                                                 C1370-
      SUBROUTINE GENPT                                                    D  10
      DOUBLE PRECISION TMRX,VPRM,HI,HR,HC,HK,WT,REC,RECH,TIM,AOPT,TITLE   D  20
      DOUBLE PRECISION XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,     D  30
     1                 HMIN,PYR                                           D  31
      INTEGER PTID                                                        D  36
      COMMON /PRMJ/ NTIM,NPMP,NPNT,NITP,N,NX,NY,NP,NREC,INT,NNX,NNY,NUMO  D  40
     1BS,NMOV,IMOV,NPMAX,ITMAX,NZCRIT,IPRNT,NPTPND,NPNTMV,NPNTVL,NPNTD,N  D  50
     2PNCHV,NPDELC,ICHK                                                   D  61
      COMMON /PRMC/ NODEID(40,40),NPCELL(20,20),NPOLD(20,20),LIMBO(500), SD  65
     1IXOBS(5),IYOBS(5)                                                   D  66
      COMMON /HEDA/ THCK(40,40),PERM(40,40),TMWL(5,50),TMOBS(50),ANFCTR  SD  90
      COMMON /HEDB/ TMRX(40,40,2),VPRM(40,40),HI(40,40),HR(40,40),HC(40, SD 100
     140),HK(40,40),WT(40,40),REC(40,40),RECH(40,40),TIM(100),AOPT(20),  SD 110
     2TITLE(10),XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,HMIN,PYR    D 120
      COMMON /HEDC/ MX,MY,MMX,MMY,NMX,NMY,MCHK                           SD 125
      COMMON /CHMA/ PART(3,6400),CONC(20,20),TMCN(5,50),VX(40,40),VY(40, SD 131
     140),CONINT(20,20),CNRECH(20,20),POROS,SUMTCH,BETA,TIMV,STORM,STORM SD 140
     2I,CMSIN,CMSOUT,FLMIN,FLMOT,SUMIO,CELDIS,DLTRAT,CSTORM               D 150
      COMMON /CHMP/ PTID(6400)                                            D 156
      DIMENSION RPT(16),RNT(16),RP(16),RN(16),IPT(16)                     D 161
C     ***************************************************************     D 170
      IONED=0                                                             D 172
      IF(NX.EQ.3.OR.NY.EQ.3) IONED=1                                      D 173
      F1=0.30                                                             D 180
      F2=1.0/3.0                                                          D 190
      IF (NPTPND.EQ.4) F1=0.25                                            D 200
      IF (NPTPND.EQ.5.AND.IONED.EQ.1) F1=0.25                             D 202
      IF (NPTPND.EQ.9) F1=1.0/3.0                                         D 210
      IF (NPTPND.EQ.8.AND.IONED.NE.1) F2=0.25                             D 221
      IF (NPTPND.EQ.8.AND.IONED.EQ.1) F1=1.0/3.0                          D 222
      IF (NPTPND.EQ.16) F1=0.25                                           D 223
      IF (NPTPND.EQ.16) F2=0.125                                          D 224
      NCHK=NPTPND                                                         D 230
      IF (NPTPND.EQ.5.OR.NPTPND.EQ.9) NCHK=NPTPND-1                       D 240
      IF (TEST.GT.98.) GO TO 10                                           D 250
C     ***************************************************************     D 260
C     ---INITIALIZE VALUES---                                             D 270
      STORM=0.0                                                           D 280
      CMSIN=0.0                                                           D 290
      CMSOUT=0.0                                                          D 300
      FLMIN=0.0                                                           D 310
      FLMOT=0.0                                                           D 320
      SUMIO=0.0                                                           D 330
C     ***************************************************************     D 340
   10 DO 20 IN=1,NPMAX                                                    D 345
      PTID(IN)=0                                                          D 355
      DO 20 ID=1,3                                                        D 365
   20 PART(ID,IN)=0.0                                                     D 370
      DO 30 IA=1,16                                                       D 381
      RP(IA)=0.0                                                          D 390
      RN(IA)=0.0                                                          D 400
      RPT(IA)=0.0                                                         D 401
      RNT(IA)=0.0                                                         D 402
   30 IPT(IA)=0                                                           D 410
C     ---SET UP LIMBO ARRAY---                                            D 420
      DO 40 IN=1,500                                                      D 430
   40 LIMBO(IN)=0                                                         D 441
      IND=1                                                               D 450
C     ***************************************************************     D 490
C     ---INSERT PARTICLES---                                              D 500
C        ---TRACK PARTICLE LOCATIONS IN COORDINATES OF PRIMARY GRID---   SD 505
      DO 410 IX=1,NMX                                                    SD 512
      JX=IX+MX-1                                                         SD 515
      DO 410 IY=1,NMY                                                    SD 522
      JY=IY+MY-1                                                         SD 525
      IF (THCK(JX,JY).EQ.0.0) GO TO 410                                  SD 531
      KR=0                                                                D 540
      KR2=0                                                               D 541
      TEST2=0.0                                                           D 550
      METH=1                                                              D 560
      NPCELL(IX,IY)=0                                                     D 570
      NPOLD(IX,IY)=NPTPND                                                 D 575
      C1=CONC(IX,IY)                                                      D 580
      IF (NPTPND.EQ.1) GO TO 150                                          D 585
      IF (C1.LE.1.0E-05) TEST2=1.0                                        D 590
      IF (IX.EQ.1.OR.IX.EQ.NMX.OR.IY.EQ.1.OR.IY.EQ.NMY) TEST2=1.0        SD 595
      IF (VPRM(JX,JY).GT.0.09) TEST2=1.0                                 SD 601
      IF (REC(JX,JY).NE.0.0) TEST2=1.0                                   SD 611
      IF (THCK(JX+1,JY+1).EQ.0.0.OR.THCK(JX+1,JY-1).EQ.0.0.OR.THCK(JX-1, SD 621
     1JY+1).EQ.0.0.OR.THCK(JX-1,JY-1).EQ.0.0) TEST2=1.0                  SD 631
      IF ((THCK(JX,JY+1).EQ.0.0.OR.THCK(JX,JY-1).EQ.0.0.OR.THCK(JX+1,JY) SD 641
     1.EQ.0.0.OR.THCK(JX-1,JY).EQ.0.0).AND.NPTPND.GT.5) TEST2=1.0        SD 651
      CNODE=C1*(1.0-F1)                                                   D 660
      IF (TEST.LT.98.0.OR.TEST2.GT.0.0) GO TO 70                          D 670
      SUMC=CONC(IX+1,IY)+CONC(IX-1,IY)+CONC(IX,IY+1)+CONC(IX,IY-1)        D 680
      IF (NCHK.EQ.4) GO TO 60                                             D 690
      SUMC=SUMC+CONC(IX+1,IY+1)+CONC(IX+1,IY-1)+CONC(IX-1,IY+1)+CONC(IX-  D 700
     11,IY-1)                                                             D 710
   60 AVC=SUMC/NCHK                                                       D 720
      IF (AVC.GT.C1) METH=2                                               D 730
C                                                                         D 740
C        ---PUT 4 PARTICLES ON CELL DIAGONALS---                          D 750
   70 DO 140 IT=1,2                                                       D 760
      EVET=(-1.0)**IT                                                     D 770
      DO 140 IS=1,2                                                       D 780
      IF(IONED.EQ.1.AND.IT.EQ.1.AND.IS.EQ.2) GO TO 140                    D 782
      IF(IONED.EQ.1.AND.IT.EQ.2.AND.IS.EQ.1) GO TO 140                    D 783
      EVES=(-1.0)**IS                                                     D 790
      KR=KR+1                                                             D 793
      IF (NPTPND.EQ.16) GO TO 72                                          D 795
      PART(1,IND)=JX+F1*EVET                                             SD 801
      PART(2,IND)=JY+F1*EVES                                             SD 811
      PART(2,IND)=-PART(2,IND)                                            D 820
      PART(3,IND)=C1                                                      D 830
      IPT(KR)=IND                                                         D 833
      PTID(IND)=KR                                                        D 834
      IF(IONED.EQ.1.AND.IS.EQ.2) PTID(IND)=4                              D 834A
      GO TO 76                                                            D 835
   72 IF (TEST.LT.98.0.OR.TEST2.GT.0.0) GO TO 135                         D 841
   76 IF (TEST.LT.98.0.OR.TEST2.GT.0.0) GO TO 139                         D 845
      IXD=IX+EVET                                                         D 850
      IYD=IY+EVES                                                         D 860
      IF (METH.EQ.2) GO TO 80                                             D 890
      PARTC=CNODE+CONC(IXD,IYD)*F1                                        D 901
      GO TO 90                                                            D 910
   80 PARTC=2.0*C1*CONC(IXD,IYD)/(C1+CONC(IXD,IYD))                       D 921
   90 IF (C1-CONC(IXD,IYD)) 100,110,120                                   D 930
  100 RPT(KR)=CONC(IXD,IYD)-PARTC                                         D 941
      RNT(KR)=C1-PARTC                                                    D 951
      GO TO 130                                                           D 960
  110 RPT(KR)=0.0                                                         D 971
      RNT(KR)=0.0                                                         D 981
      GO TO 130                                                           D 990
  120 RPT(KR)=C1-PARTC                                                    D1001
      RNT(KR)=CONC(IXD,IYD)-PARTC                                         D1011
  130 IF (NPTPND.EQ.16) GO TO 135                                         D1012
      PART(3,IND)=PARTC                                                   D1013
      RP(KR)=RPT(KR)                                                      D1014
      RN(KR)=RNT(KR)                                                      D1015
      GO TO 139                                                           D1016
  135 DO 138 ITT = 1,2                                                    D1017
      EVET2=(-1.0)**ITT                                                   D1018
      DO 138 ISS = 1,2                                                    D1019
      IF(IONED.EQ.1.AND.ITT.EQ.1.AND.ISS.EQ.2) GO TO 138                  D1019A
      IF(IONED.EQ.1.AND.ITT.EQ.2.AND.ISS.EQ.1) GO TO 138                  D1019B
      EVES2=(-1.0)**ISS                                                   D1021
      PART(1,IND)=(JX+F1*EVET)+F2*EVET2                                  SD1022A
      PART(2,IND)=(JY+F1*EVES)+F2*EVES2                                  SD1023A
      PART(2,IND)=-PART(2,IND)                                            D1024
      KR2=KR2+1                                                           D1025
      IF (TEST.LT.98.0.OR.TEST2.GT.0.0) GO TO 136                         D1026
      PART(3,IND) = PARTC                                                 D1027
      RP(KR2) = RPT(KR)                                                   D1028
      RN(KR2) = RNT(KR)                                                   D1029
      IPT(KR2) = IND                                                      D1031
      GO TO 137                                                           D1032
  136 PART(3,IND) = C1                                                    D1033
  137 PTID(IND)=KR2                                                       D1034
      IF(IONED.EQ.1) THEN                                                 D1034A
       IF(IT.EQ.1.AND.ISS.EQ.2) PTID(IND)=4                               D1034B
       IF(IT.EQ.2) THEN                                                   D1034C
        IF(ISS.EQ.1) PTID(IND)=13                                         D1034D
        IF(ISS.EQ.2) PTID(IND)=16                                         D1034E
       END IF                                                             D1034F
      END IF                                                              D1034G
      IND=IND+1                                                           D1035
  138 CONTINUE                                                            D1036
      GO TO 140                                                           D1037
  139 IND=IND+1                                                           D1038
  140 CONTINUE                                                            D1039
      IF (NPTPND.EQ.16) GO TO 290                                         D1045
      IF ((NPTPND.EQ.5.AND.IONED.NE.1).OR.NPTPND.EQ.9) GO TO 150          D1051
      IF (NPTPND.EQ.8.AND.IONED.EQ.1) GO TO 150                           D1052
      GO TO 160                                                           D1060
C        ---PUT ONE PARTICLE AT CENTER OF CELL---                         D1070
  150 PART(1,IND)=JX                                                     SD1076
      PART(2,IND)=-JY                                                    SD1091
      PART(3,IND)=C1                                                      D1100
      PTID(IND)=5                                                         D1105
      IND=IND+1                                                           D1110
      IF (NPTPND.EQ.1) GO TO 410                                          D1115
C        ---PLACE NORTH, SOUTH, EAST, AND WEST PARTICLES---               D1120
  160 IF (NPTPND.LT.8) GO TO 290                                          D1130
      IF(IONED.EQ.1) GO TO 290                                            D1132
      CNODE=C1*(1.0-F2)                                                   D1140
      DO 280 IT=1,2                                                       D1150
      EVET=(-1.0)**IT                                                     D1160
      PART(1,IND)=JX+F2*EVET                                             SD1171
      PART(2,IND)=-JY                                                    SD1181
      PART(3,IND)=C1                                                      D1190
      IF (EVET.LT.0) PTID(IND)=6                                          D1192
      IF (EVET.GT.0) PTID(IND)=8                                          D1194
      IF (TEST.LT.98.0.OR.TEST2.GT.0.0) GO TO 220                         D1200
      IXD=IX+EVET                                                         D1210
      KR=KR+1                                                             D1220
      IPT(KR)=IND                                                         D1230
      IF (METH.EQ.2) GO TO 170                                            D1240
      PART(3,IND)=CNODE+CONC(IXD,IY)*F2                                   D1250
      GO TO 180                                                           D1260
  170 PART(3,IND)=2.0*C1*CONC(IXD,IY)/(C1+CONC(IXD,IY))                   D1270
  180 IF (C1-CONC(IXD,IY)) 190,200,210                                    D1280
  190 RP(KR)=CONC(IXD,IY)-PART(3,IND)                                     D1290
      RN(KR)=C1-PART(3,IND)                                               D1300
      GO TO 220                                                           D1310
  200 RP(KR)=0.0                                                          D1320
      RN(KR)=0.0                                                          D1330
      GO TO 220                                                           D1340
  210 RP(KR)=C1-PART(3,IND)                                               D1350
      RN(KR)=CONC(IXD,IY)-PART(3,IND)                                     D1360
  220 IND=IND+1                                                           D1370
      PART(1,IND)=JX                                                     SD1381
      PART(2,IND)=JY+F2*EVET                                             SD1391
      PART(2,IND)=-PART(2,IND)                                            D1400
      PART(3,IND)=C1                                                      D1410
      IF (EVET.LT.0) PTID(IND)=7                                          D1412
      IF (EVET.GT.0) PTID(IND)=9                                          D1414
      IF (TEST.LT.98.0.OR.TEST2.GT.0.0) GO TO 280                         D1420
      IYD=IY+EVET                                                         D1430
      KR=KR+1                                                             D1440
      IPT(KR)=IND                                                         D1450
      IF (METH.EQ.2) GO TO 230                                            D1460
      PART(3,IND)=CNODE+CONC(IX,IYD)*F2                                   D1470
      GO TO 240                                                           D1480
  230 PART(3,IND)=2.0*C1*CONC(IX,IYD)/(C1+CONC(IX,IYD))                   D1490
  240 IF (C1-CONC(IX,IYD)) 250,260,270                                    D1500
  250 RP(KR)=CONC(IX,IYD)-PART(3,IND)                                     D1510
      RN(KR)=C1-PART(3,IND)                                               D1520
      GO TO 280                                                           D1530
  260 RP(KR)=0.0                                                          D1540
      RN(KR)=0.0                                                          D1550
      GO TO 280                                                           D1560
  270 RP(KR)=C1-PART(3,IND)                                               D1570
      RN(KR)=CONC(IX,IYD)-PART(3,IND)                                     D1580
  280 IND=IND+1                                                           D1590
C                                                                         D1600
  290 IF (TEST.LT.98.0.OR.TEST2.GT.0.0) GO TO 410                         D1610
      SUMPT=0.0                                                           D1620
C        ---COMPUTE CONC. GRADIENT WITHIN CELL---                         D1630
      DO 300 KPT=1,NCHK                                                   D1640
      IK=IPT(KPT)                                                         D1650
  300 SUMPT=PART(3,IK)+SUMPT                                              D1660
      CBAR=SUMPT/NCHK                                                     D1670
C        ---CHECK MASS BALANCE WITHIN CELL AND ADJUST PT. CONCS.---       D1680
      SUMPT=0.0                                                           D1690
      IF (CBAR-C1) 310,410,330                                            D1700
  310 CRCT=1.0-(CBAR/C1)                                                  D1710
      IF (METH.EQ.1) CRCT=CBAR/C1                                         D1720
      DO 320 KPT=1,NCHK                                                   D1730
      IK=IPT(KPT)                                                         D1740
      PART(3,IK)=PART(3,IK)+RP(KPT)*CRCT                                  D1750
  320 SUMPT=SUMPT+PART(3,IK)                                              D1760
      CBARN=SUMPT/NCHK                                                    D1770
      GO TO 350                                                           D1780
  330 CRCT=1.0-(C1/CBAR)                                                  D1790
      IF (METH.EQ.1) CRCT=C1/CBAR                                         D1800
      DO 340 KPT=1,NCHK                                                   D1810
      IK=IPT(KPT)                                                         D1820
      PART(3,IK)=PART(3,IK)+RN(KPT)*CRCT                                  D1830
  340 SUMPT=SUMPT+PART(3,IK)                                              D1840
      CBARN=SUMPT/NCHK                                                    D1850
  350 IF (CBARN.EQ.C1) GO TO 410                                          D1860
C        ---CORRECT FOR OVERCOMPENSATION---                               D1870
      CRCT=C1/CBARN                                                       D1880
      DO 380 KPT=1,NCHK                                                   D1890
      IK=IPT(KPT)                                                         D1900
      PART(3,IK)=PART(3,IK)*CRCT                                          D1910
C        ---CHECK CONSTRAINTS---                                          D1920
      IF (PART(3,IK)-C1) 360,380,370                                      D1930
  360 CLIM=C1-RP(KPT)+RN(KPT)                                             D1940
      IF (PART(3,IK).LT.CLIM) GO TO 390                                   D1950
      GO TO 380                                                           D1960
  370 CLIM=C1+RP(KPT)-RN(KPT)                                             D1970
      IF (PART(3,IK).GT.CLIM) GO TO 390                                   D1980
  380 CONTINUE                                                            D1990
      GO TO 410                                                           D2000
  390 TEST2=1.0                                                           D2010
      DO 400 KPT=1,NCHK                                                   D2020
      IK=IPT(KPT)                                                         D2030
  400 PART(3,IK)=C1                                                       D2040
  410 CONTINUE                                                            D2050
      NP=IND-1                                                            D2061
      IF (INT.EQ.0) CALL CHMOT                                            D2070
C     ****************************************************************    D2080
      RETURN                                                              D2090
C     ****************************************************************    D2100
      END                                                                 D2110-
