      SUBROUTINE MOVE                                                     F  10
      DOUBLE PRECISION TMRX,VPRM,HI,HR,HC,HK,WT,REC,RECH,TIM,AOPT,TITLE   F  20
      DOUBLE PRECISION XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,     F  30
     1                 HMIN,PYR                                           F  31
      DOUBLE PRECISION DXINV,DYINV,ARINV,PORINV                           F  32
      DOUBLE PRECISION DCYFCT,DCYT,DCYT2                                  F  33
      INTEGER PTID                                                        F  36
      COMMON /PRMJ/ NTIM,NPMP,NPNT,NITP,N,NX,NY,NP,NREC,INT,NNX,NNY,NUMO  F  40
     1BS,NMOV,IMOV,NPMAX,ITMAX,NZCRIT,IPRNT,NPTPND,NPNTMV,NPNTVL,NPNTD,N  F  50
     2PNCHV,NPDELC,ICHK                                                   F  61
      COMMON /PRMC/ NODEID(40,40),NPCELL(20,20),NPOLD(20,20),LIMBO(500), SF  65
     1IXOBS(5),IYOBS(5)                                                   F  66
      COMMON /HEDA/ THCK(40,40),PERM(40,40),TMWL(5,50),TMOBS(50),ANFCTR  SF  90
      COMMON /HEDB/ TMRX(40,40,2),VPRM(40,40),HI(40,40),HR(40,40),HC(40, SF 100
     140),HK(40,40),WT(40,40),REC(40,40),RECH(40,40),TIM(100),AOPT(20),  SF 110
     2TITLE(10),XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,HMIN,PYR    F 120
      COMMON /HEDC/ MX,MY,MMX,MMY,NMX,NMY,MCHK                           SF 125
      COMMON /XINV/ DXINV,DYINV,ARINV,PORINV                              F 130
      COMMON /CHMA/ PART(3,6400),CONC(20,20),TMCN(5,50),VX(40,40),VY(40, SF 141
     140),CONINT(20,20),CNRECH(20,20),POROS,SUMTCH,BETA,TIMV,STORM,STORM SF 150
     2I,CMSIN,CMSOUT,FLMIN,FLMOT,SUMIO,CELDIS,DLTRAT,CSTORM               F 160
      COMMON /CHMC/ SUMC(20,20),VXBDY(40,40),VYBDY(40,40),PERM2(40,40)   SF 171
      COMMON /CHMR/ RF,DK,RHOB,THALF,DECAY,ADSORB,SORBI,DMASS1,CSTM2      F 172R
      COMMON /CHMP/ PTID(6400)                                            F 176
C     ***************************************************************     F 190
      WRITE (6,680) NMOV                                                  F 200
      IF (THALF.GT.0.0.AND.THALF.LT.TIMV) WRITE (6,685)                   F 205R
      SUMTCH=SUMT-TIM(N)                                                  F 210
      IONED=0                                                             F 211A
      IF(NX.EQ.3.OR.NY.EQ.3) IONED=1                                      F 211B
      F1=0.30                                                             F 212
      F2=1.0/3.0                                                          F 214
      IF (NPTPND.EQ.4) F1=0.25                                            F 216
      IF (NPTPND.EQ.5.AND.IONED.EQ.1) F1=0.25                             F 217
      IF (NPTPND.EQ.9) F1=F2                                              F 218
      IF (NPTPND.EQ.8.AND.IONED.NE.1) F2=0.25                             F 222A
      IF (NPTPND.EQ.8.AND.IONED.EQ.1) F1=F2                               F 222B
      IF (NPTPND.EQ.16) F1=0.25                                           F 223
      IF (NPTPND.EQ.16) F2=0.125                                          F 224
      CONST1=TIMV*DXINV/RF                                                F 251R
      CONST2=TIMV*DYINV/RF                                                F 261R
      DCYFCT=TIMV*DECAY                                                   F 265R
      DCYT=1.D0                                                           F 268A
      DCYT2=1.D0                                                          F 268B
      IF(DECAY.NE.0.0) THEN                                               F 268C
       DCYT=DEXP(-DCYFCT)                                                 F 268D
       DCYT2=DEXP(-DCYFCT*0.5D0)                                          F 268E
      END IF                                                              F 268F
C     ---MOVE PARTICLES 'NMOV' TIMES---                                   F 270
      DO 650 IMOV=1,NMOV                                                  F 280
   10 NPTM=NP                                                             F 290
C        ---MOVE EACH PARTICLE---                                         F 300
      DO 590 IN=1,NP                                                      F 310
      IF (PART(1,IN).EQ.0.0) GO TO 590                                    F 320
      KFLAG=0                                                             F 330
      NLOC=0                                                             SF 335
C     ***************************************************************     F 340
C           ---COMPUTE OLD LOCATION---                                    F 350
   20 XOLD=PART(1,IN)                                                     F 400
      IX=XOLD+0.5                                                         F 410
      JX=IX-MX+1                                                         SF 415
      IFLAG=1                                                             F 420
      IF (PART(2,IN).GE.0.0) GO TO 30                                     F 430
      IFLAG=-1                                                            F 440
      PART(2,IN)=-PART(2,IN)                                              F 450
   30 YOLD=PART(2,IN)                                                     F 460
      IY=YOLD+0.5                                                         F 470
      JY=IY-MY+1                                                         SF 475
      IF (THCK(IX,IY).EQ.0.0) GO TO 590                                   F 482
C     ***************************************************************     F 490
C           ---COMPUTE NEW LOCATION AND LOCATE CLOSEST NODE---            F 500
C           ---LOCATE NORTHWEST CORNER---                                 F 510
      IVX=XOLD                                                            F 520
      IVY=YOLD                                                            F 530
      IXE=IVX+1                                                           F 540
      IYS=IVY+1                                                           F 550
C     ***************************************************************     F 560
C        ---LOCATE QUADRANT, VEL. AT 4 CORNERS, CHECK FOR BOUNDARIES---   F 570
      CELDX=XOLD-IX                                                       F 580
      CELDY=YOLD-IY                                                       F 590
      ICD=9                                                               F 595
      IF (CELDX.EQ.0.0.AND.CELDY.EQ.0.0) GO TO 280                        F 600
      IF (CELDX.GE.0.0.OR.CELDY.GE.0.0) GO TO 70                          F 610
C              ---PT. IN NW QUADRANT---                                   F 620
      VXNW=VXBDY(IVX,IVY)                                                 F 630
      VXNE=VX(IXE,IVY)                                                    F 640
      VXSW=VXBDY(IVX,IYS)                                                 F 650
      VXSE=VX(IXE,IYS)                                                    F 660
      VYNW=VYBDY(IVX,IVY)                                                 F 670
      VYNE=VYBDY(IXE,IVY)                                                 F 680
      VYSW=VY(IVX,IYS)                                                    F 690
      VYSE=VY(IXE,IYS)                                                    F 700
      ICD=1                                                               F 705
      IF (THCK(IVX,IVY).EQ.0.0) GO TO 50                                  F 710
      IF (REC(IXE,IVY).EQ.0.0.AND.VPRM(IXE,IVY).LT.0.09) GO TO 40         F 720
      IF (THCK(IXE+1,IVY).NE.0.0) VXNE=VXNW                               F 731
   40 IF (REC(IVX,IYS).EQ.0.0.AND.VPRM(IVX,IYS).LT.0.09) GO TO 50         F 740
      IF (THCK(IVX,IYS+1).NE.0.0) VYSW=VYNW                               F 751
   50 IF (REC(IXE,IYS).EQ.0.0.AND.VPRM(IXE,IYS).LT.0.09) GO TO 270        F 760
      IF (THCK(IVX,IYS).EQ.0.0) GO TO 60                                  F 770
      IF (THCK(IXE+1,IYS).GT.0.0) VXSE=VXSW                               F 782
   60 IF (THCK(IXE,IVY).EQ.0.0) GO TO 270                                 F 790
      IF (THCK(IXE,IYS+1).GT.0.0) VYSE=VYNE                               F 802
      GO TO 270                                                           F 810
C                                                                         F 820
   70 IF (CELDX.LE.0.0.OR.CELDY.GE.0.0) GO TO 130                         F 830
C              ---PT. IN NE QUADRANT---                                   F 840
   80 VXNW=VX(IVX,IVY)                                                    F 850
      VXNE=VXBDY(IVX,IVY)                                                 F 860
      VXSW=VX(IVX,IYS)                                                    F 870
      VXSE=VXBDY(IVX,IYS)                                                 F 880
      VYNW=VYBDY(IVX,IVY)                                                 F 890
      VYNE=VYBDY(IXE,IVY)                                                 F 900
      VYSW=VY(IVX,IYS)                                                    F 910
      VYSE=VY(IXE,IYS)                                                    F 920
      ICD=2                                                               F 925
      IF (CELDX.EQ.0.0) GO TO 120                                         F 930
      IF (THCK(IXE,IVY).EQ.0.0) GO TO 100                                 F 940
      IF (REC(IVX,IVY).EQ.0.0.AND.VPRM(IVX,IVY).LT.0.09) GO TO 90         F 950
      IF (THCK(IVX-1,IVY).NE.0.0) VXNW=VXNE                               F 961
   90 IF (REC(IXE,IYS).EQ.0.0.AND.VPRM(IXE,IYS).LT.0.09) GO TO 100        F 970
      IF (THCK(IXE,IYS+1).NE.0.0) VYSE=VYNE                               F 981
  100 IF (REC(IVX,IYS).EQ.0.0.AND.VPRM(IVX,IYS).LT.0.09) GO TO 270        F 990
      IF (THCK(IXE,IYS).EQ.0.0) GO TO 110                                 F1000
      IF (THCK(IVX-1,IYS).GT.0.0) VXSW=VXSE                               F1012
  110 IF (THCK(IVX,IVY).EQ.0.0) GO TO 270                                 F1020
      IF (THCK(IVX,IYS+1).GT.0.0) VYSW=VYNW                               F1032
      GO TO 270                                                           F1040
  120 IF (REC(IVX,IYS).EQ.0.0.AND.VPRM(IVX,IYS).LE.0.09) GO TO 270        F1050
      IF (THCK(IVX,IVY).EQ.0.0) GO TO 270                                 F1060
      IF (THCK(IVX,IYS+1).GT.0.0) VYSW=VYNW                               F1072
      GO TO 270                                                           F1080
C                                                                         F1090
  130 IF (CELDY.LE.0.0.OR.CELDX.GE.0.0) GO TO 190                         F1100
C              ---PT. IN SW QUADRANT---                                   F1110
  140 VXNW=VXBDY(IVX,IVY)                                                 F1120
      VXNE=VX(IXE,IVY)                                                    F1130
      VXSW=VXBDY(IVX,IYS)                                                 F1140
      VXSE=VX(IXE,IYS)                                                    F1150
      VYNW=VY(IVX,IVY)                                                    F1160
      VYNE=VY(IXE,IVY)                                                    F1170
      VYSW=VYBDY(IVX,IVY)                                                 F1180
      VYSE=VYBDY(IXE,IVY)                                                 F1190
      ICD=3                                                               F1195
      IF (CELDY.EQ.0.0) GO TO 180                                         F1200
      IF (THCK(IVX,IYS).EQ.0.0) GO TO 160                                 F1210
      IF (REC(IVX,IVY).EQ.0.0.AND.VPRM(IVX,IVY).LT.0.09) GO TO 150        F1220
      IF (THCK(IVX,IVY-1).NE.0.0) VYNW=VYSW                               F1231
  150 IF (REC(IXE,IYS).EQ.0.0.AND.VPRM(IXE,IYS).LT.0.09) GO TO 160        F1240
      IF (THCK(IXE+1,IYS).NE.0.0) VXSE=VXSW                               F1251
  160 IF (REC(IXE,IVY).EQ.0.0.AND.VPRM(IXE,IVY).LT.0.09) GO TO 270        F1260
      IF (THCK(IVX,IVY).EQ.0.0) GO TO 170                                 F1270
      IF (THCK(IXE+1,IVY).GT.0.0) VXNE=VXNW                               F1282
  170 IF (THCK(IXE,IYS).EQ.0.0) GO TO 270                                 F1290
      IF (THCK(IXE,IVY-1).GT.0.0) VYNE=VYSE                               F1302
      GO TO 270                                                           F1310
  180 IF (REC(IXE,IVY).EQ.0.0.AND.VPRM(IXE,IVY).LE.0.09) GO TO 270        F1320
      IF (THCK(IVX,IVY).EQ.0.0) GO TO 270                                 F1330
      IF (THCK(IXE+1,IVY).GT.0.0) VXNE=VXNW                               F1342
      GO TO 270                                                           F1350
C                                                                         F1360
  190 IF (CELDY.LE.0.0.OR.CELDX.LE.0.0) GO TO 260                         F1370
C              ---PT. IN SE QUADRANT---                                   F1380
  200 VXNW=VX(IVX,IVY)                                                    F1390
      VXNE=VXBDY(IVX,IVY)                                                 F1400
      VXSW=VX(IVX,IYS)                                                    F1410
      VXSE=VXBDY(IVX,IYS)                                                 F1420
      VYNW=VY(IVX,IVY)                                                    F1430
      VYNE=VY(IXE,IVY)                                                    F1440
      VYSW=VYBDY(IVX,IVY)                                                 F1450
      VYSE=VYBDY(IXE,IVY)                                                 F1460
      ICD=4                                                               F1465
      IF (CELDY.EQ.0.0) GO TO 240                                         F1470
      IF (CELDX.EQ.0.0) GO TO 250                                         F1480
      IF (THCK(IXE,IYS).EQ.0.0) GO TO 220                                 F1490
      IF (REC(IXE,IVY).EQ.0.0.AND.VPRM(IXE,IVY).LT.0.09) GO TO 210        F1500
      IF (THCK(IXE,IVY-1).NE.0.0) VYNE=VYSE                               F1511
  210 IF (REC(IVX,IYS).EQ.0.0.AND.VPRM(IVX,IYS).LT.0.09) GO TO 220        F1520
      IF (THCK(IVX-1,IYS).NE.0.0) VXSW=VXSE                               F1531
  220 IF (REC(IVX,IVY).EQ.0.0.AND.VPRM(IVX,IVY).LT.0.09) GO TO 270        F1540
      IF (THCK(IXE,IVY).EQ.0.0) GO TO 230                                 F1550
      IF (THCK(IVX-1,IVY).GT.0.0) VXNW=VXNE                               F1562
  230 IF (THCK(IVX,IYS).EQ.0.0) GO TO 270                                 F1570
      IF (THCK(IVX,IVY-1).GT.0.0) VYNW=VYSW                               F1582
      GO TO 270                                                           F1590
  240 IF (REC(IVX,IVY).EQ.0.0.AND.VPRM(IVX,IVY).LE.0.09) GO TO 270        F1600
      IF (THCK(IXE,IVY).EQ.0.0) GO TO 270                                 F1610
      IF (THCK(IVX-1,IVY).GT.0.0) VXNW=VXNE                               F1622
      GO TO 270                                                           F1630
  250 IF (REC(IVX,IVY).EQ.0.0.AND.VPRM(IVX,IVY).LE.0.09) GO TO 270        F1640
      IF (THCK(IVX,IYS).EQ.0.0) GO TO 270                                 F1650
      IF (THCK(IVX,IVY-1).GT.0.0) VYNW=VYSW                               F1662
      GO TO 270                                                           F1670
C                                                                         F1680
  260 IF (CELDX.EQ.0.0.AND.CELDY.LT.0.0) GO TO 80                         F1690
      IF (CELDX.LT.0.0.AND.CELDY.EQ.0.0) GO TO 140                        F1700
      IF (CELDX.GT.0.0.AND.CELDY.EQ.0.0) GO TO 200                        F1710
      IF (CELDX.EQ.0.0.AND.CELDY.GT.0.0) GO TO 200                        F1720
      WRITE (6,690) IN,IX,IY                                              F1730
  270 CONTINUE                                                            F1740
C              --- CHECK FOR ADJACENT NO-FLOW BOUNDARIES---               F1741A
      GO TO (1270,1275,1280,1285,1290) ICD                                F1741B
      GO TO 1290                                                          F1741C
 1270 IF (THCK(IXE,IVY).EQ.0.0) GO TO 1272                                F1742A
      IF (THCK(IVX,IYS).EQ.0.0) GO TO 1273                                F1742B
      IF (THCK(IVX,IVY).EQ.0.0) GO TO 1274                                F1742C
      GO TO 1290                                                          F1742D
 1272 VXNE=VXSE                                                           F1742E
      IF (THCK(IVX,IYS).GT.0.0) GO TO 1274                                F1742F
 1273 VYSW=VYSE                                                           F1742G
 1274 VXNW=VXSW                                                           F1742H
      VYNW=VYNE                                                           F1742I
      GO TO 1290                                                          F1742J
 1275 IF (THCK(IVX,IVY).EQ.0.0) GO TO 1277                                F1744A
      IF (THCK(IXE,IYS).EQ.0.0) GO TO 1278                                F1744B
      IF (THCK(IXE,IVY).EQ.0.0) GO TO 1279                                F1744C
      GO TO 1290                                                          F1744D
 1277 VXNW=VXSW                                                           F1744E
      IF (THCK(IXE,IYS).GT.0.0) GO TO 1279                                F1744F
 1278 VYSE=VYSW                                                           F1744G
 1279 VXNE=VXSE                                                           F1744H
      VYNE=VYNW                                                           F1744I
      GO TO 1290                                                          F1744J
 1280 IF (THCK(IXE,IYS).EQ.0.0) GO TO 1282                                F1746A
      IF (THCK(IVX,IVY).EQ.0.0) GO TO 1283                                F1746B
      IF (THCK(IVX,IYS).EQ.0.0) GO TO 1284                                F1746C
      GO TO 1290                                                          F1746D
 1282 VXSE=VXNE                                                           F1746E
      IF (THCK(IVX,IVY).GT.0.0) GO TO 1284                                F1746F
 1283 VYNW=VYNE                                                           F1746G
 1284 VXSW=VXNW                                                           F1746H
      VYSW=VYSE                                                           F1746I
      GO TO 1290                                                          F1746J
 1285 IF (THCK(IVX,IYS).EQ.0.0) GO TO 1287                                F1748A
      IF (THCK(IXE,IVY).EQ.0.0) GO TO 1288                                F1748B
      IF (THCK(IXE,IYS).EQ.0.0) GO TO 1289                                F1748C
      GO TO 1290                                                          F1748D
 1287 VXSW=VXNW                                                           F1748E
      IF (THCK(IXE,IVY).GT.0.0) GO TO 1289                                F1748F
 1288 VYNE=VYNW                                                           F1748G
 1289 VYSE=VYSW                                                           F1748H
      VXSE=VXNE                                                           F1748I
 1290 CONTINUE                                                            F1749A
C     ****************************************************************    F1750
C        ---BILINEAR INTERPOLATION---                                     F1760
      CELXD=XOLD-IVX                                                      F1770
      CELDXH=AMOD(CELXD,0.5)                                              F1780
      CELDX=CELDXH*2.0                                                    F1790
      CELDY=YOLD-IVY                                                      F1800
C     ****************************************************************    F1810
C           ---X VELOCITY---                                              F1820
      VXN=VXNW*(1.0-CELDX)+VXNE*CELDX                                     F1830
      VXS=VXSW*(1.0-CELDX)+VXSE*CELDX                                     F1850
      XVEL=VXN*(1.0-CELDY)+VXS*CELDY                                      F1870
C           ---Y VELOCITY---                                              F1900
      CELDYH=AMOD(CELDY,0.5)                                              F1910
      CELYD=CELDYH*2.0                                                    F1921
      VYW=VYNW*(1.0-CELYD)+VYSW*CELYD                                     F1931
      VYE=VYNE*(1.0-CELYD)+VYSE*CELYD                                     F1951
      YVEL=VYW*(1.0-CELXD)+VYE*CELXD                                      F1970
C                                                                         F2000
      GO TO 290                                                           F2010
  280 XVEL=VX(IX,IY)                                                      F2020
      YVEL=VY(IX,IY)                                                      F2030
  290 DISTX=XVEL*CONST1                                                   F2040
      DISTY=YVEL*CONST2                                                   F2050
C     ****************************************************************    F2060
C        ---BOUNDARY CONDITIONS---                                        F2070
      TEMPX=XOLD+DISTX                                                    F2080
      TEMPY=YOLD+DISTY                                                    F2090
      INX=TEMPX+0.5                                                       F2100
      INY=TEMPY+0.5                                                       F2110
      IF (THCK(INX,INY).GT.0.0) GO TO 330                                 F2120
C     ****************************************************************    F2130
C           ---X BOUNDARY---                                              F2140
      IF (THCK(INX,IY).EQ.0.0) GO TO 300                                  F2150
      PART(1,IN)=TEMPX                                                    F2160
      GO TO 310                                                           F2170
  300 BEYON=TEMPX-IX                                                      F2180
      IF (BEYON.LT.0.0) BEYON=BEYON+0.5                                   F2190
      IF (BEYON.GT.0.0) BEYON=BEYON-0.5                                   F2200
      PART(1,IN)=TEMPX-2.0*BEYON                                          F2210
      INX=PART(1,IN)+0.5                                                  F2220
      TEMPX=PART(1,IN)                                                    F2230
C     ****************************************************************    F2240
C           ---Y BOUNDARY---                                              F2250
  310 IF (THCK(INX,INY).EQ.0.0) GO TO 320                                 F2260
      PART(2,IN)=TEMPY                                                    F2270
      GO TO 340                                                           F2280
C     ****************************************************************    F2290
  320 BEYON=TEMPY-IY                                                      F2300
      IF (BEYON.LT.0.0) BEYON=BEYON+0.5                                   F2310
      IF (BEYON.GT.0.0) BEYON=BEYON-0.5                                   F2320
      PART(2,IN)=TEMPY-2.0*BEYON                                          F2330
      INY=PART(2,IN)+0.5                                                  F2340
      TEMPY=PART(2,IN)                                                    F2350
      GO TO 340                                                           F2360
  330 PART(1,IN)=TEMPX                                                    F2370
      PART(2,IN)=TEMPY                                                    F2380
  340 CONTINUE                                                            F2390
      JNX=INX-MX+1                                                       SF2395
      JNY=INY-MY+1                                                       SF2396
      IF (MCHK.EQ.0) GO TO 342                                           SF2397
      IF (JNX.LT.1.OR.JNX.GT.NMX.OR.JNY.LT.1.OR.JNY.GT.NMY) NLOC=1       SF2398
C     ***************************************************************     F2400
      IF (NLOC.GT.0) GO TO 345                                           SF2405
  342 CONTINUE                                                           SF2407
C        ---SUM CONCENTRATIONS AND COUNT PARTICLES---                     F2410
C           ---DECAY PARTICLES---                                         F2414R
      PART(3,IN)=PART(3,IN)*DCYT                                          F2416R
      SUMC(JNX,JNY)=SUMC(JNX,JNY)+PART(3,IN)                             SF2421
      NPCELL(JNX,JNY)=NPCELL(JNX,JNY)+1                                  SF2431
C     ****************************************************************    F2440
  345 CONTINUE                                                           SF2445
C        ---CHECK FOR CHANGE IN CELL LOCATION---                          F2450
      IF (IX.EQ.INX.AND.IY.EQ.INY) GO TO 580                              F2460
C           ---CHECK FOR CONST.-HEAD BDY. OR SOURCE AT OLD LOCATION---    F2470
      IBD=0                                                              SF2478
C     IF (REC(IX,IY).LT.-0.1) GO TO 350                                    FXXX
C     IF (REC(IX,IY).GT.0.1) GO TO 360                                    FXXX
      IF (REC(IX,IY).LT.0.0) GO TO 350                                    F2480
      IF (REC(IX,IY).GT.0.0) GO TO 360                                    F2490
      IF (VPRM(IX,IY).LT.0.09) GO TO 348                                 SF2501
      IF (WT(IX,IY).GT.HK(IX,IY)) GO TO 350                               F2510
      IF (WT(IX,IY).LT.HK(IX,IY)) GO TO 360                               F2520
  348 IF (MCHK.EQ.0) GO TO 540                                           SF2522
      IF (JX.EQ.1.AND.VXBDY(IX-1,IY).GT.0.0.AND.JNX.GT.JX) IBD=1         SF2523
      IF (JX.EQ.NMX.AND.VXBDY(IX,IY).LT.0.0.AND.JNX.LT.JX) IBD=1         SF2524
      IF (JY.EQ.1.AND.VYBDY(IX,IY-1).GT.0.0.AND.JNY.GT.JY) IBD=1         SF2525
      IF (JY.EQ.NMY.AND.VYBDY(IX,IY).LT.0.0.AND.JNY.LT.JY) IBD=1         SF2526
      IF (IBD.EQ.1) GO TO 350                                            SF2527
      GO TO 540                                                           F2530
C     ****************************************************************    F2540
C           ---CREATE NEW PARTICLES AT BOUNDARIES---                      F2550
  350 IF (IFLAG.GT.0) GO TO 550                                           F2560
      KFLAG=1                                                             F2570
  360 CONTINUE                                                            F2581
C     ****************************************************************    F2630
C           ---GENERATE NEW TEMPORARY PARTICLE---                         F2641
      IF (NPTM.EQ.NPMAX) GO TO 600                                        F2650
      NPTM=NPTM+1                                                         F2660
      IP=NPTM                                                             F2670
C                                                                         F2700
  390 IF (KFLAG.EQ.0) GO TO 398                                           F2705
      IF (IBD.EQ.1) GO TO 525                                            SF2708
      ITEM=PTID(IN)                                                       F2845
      GO TO 399                                                           F2855
  398 SUMC(JX,JY)=SUMC(JX,JY)+CONC(JX,JY)*DCYT2                          SF2868R
      NPCELL(JX,JY)=NPCELL(JX,JY)+1                                      SF2876
      IF (NPOLD(JX,JY).GT.0) NPOLD(JX,JY)=NPOLD(JX,JY)-1                 SF2886
      IF (IFLAG.GT.0) GO TO 441                                           F2895
      IF (KFLAG.EQ.0) GO TO 441                                           F2899
  399 IF (NPTPND.EQ.16) GO TO 400                                         F2901
      GO TO (401,411,421,431,441,451,461,471,481),ITEM                    F2906
      GO TO 441                                                           F2915
  400 GO TO (482,483,484,485,486,487,488,489,490,491,492,493,494,495,     F2921
     1       496,497),ITEM                                                F2922
      GO TO 441                                                           F2924
  401 PART(1,IP)=IX-F1                                                    F2925
      PART(2,IP)=IY-F1                                                    F2935
      PTID(IP)=1                                                          F2945
      GO TO 530                                                           F2955
  411 PART(1,IP)=IX-F1                                                    F2965
      PART(2,IP)=IY+F1                                                    F2975
      PTID(IP)=2                                                          F2985
      GO TO 530                                                           F2995
  421 PART(1,IP)=IX+F1                                                    F3005
      PART(2,IP)=IY-F1                                                    F3015
      PTID(IP)=3                                                          F3025
      GO TO 530                                                           F3035
  431 PART(1,IP)=IX+F1                                                    F3045
      PART(2,IP)=IY+F1                                                    F3055
      PTID(IP)=4                                                          F3065
      GO TO 530                                                           F3075
  441 PART(1,IP)=IX                                                       F3085
      PART(2,IP)=IY                                                       F3095
      PTID(IP)=5                                                          F3105
      GO TO 530                                                           F3115
  451 PART(1,IP)=IX-F2                                                    F3125
      PART(2,IP)=IY                                                       F3135
      PTID(IP)=6                                                          F3145
      GO TO 530                                                           F3155
  461 PART(1,IP)=IX                                                       F3165
      PART(2,IP)=IY-F2                                                    F3175
      PTID(IP)=7                                                          F3185
      GO TO 530                                                           F3195
  471 PART(1,IP)=IX+F2                                                    F3205
      PART(2,IP)=IY                                                       F3215
      PTID(IP)=8                                                          F3225
      GO TO 530                                                           F3235
  481 PART(1,IP)=IX                                                       F3245
      PART(2,IP)=IY+F2                                                    F3255
      PTID(IP)=9                                                          F3265
      GO TO 530                                                           F3266
  482 PART(1,IP)=IX-F1-F2                                                 F3267
      PART(2,IP)=IY-F1-F2                                                 F3268
      PTID(IP)=1                                                          F3269
      GO TO 530                                                           F3271
  483 PART(1,IP)=IX-F1-F2                                                 F3272
      PART(2,IP)=IY-F1+F2                                                 F3273
      PTID(IP)=2                                                          F3274
      GO TO 530                                                           F3275
  484 PART(1,IP)=IX-F1+F2                                                 F3276
      PART(2,IP)=IY-F1-F2                                                 F3277
      PTID(IP)=3                                                          F3278
      GO TO 530                                                           F3279
  485 PART(1,IP)=IX-F1+F2                                                 F3281
      PART(2,IP)=IY-F1+F2                                                 F3282
      PTID(IP)=4                                                          F3283
      GO TO 530                                                           F3284
  486 PART(1,IP)=IX-F1-F2                                                 F3285
      PART(2,IP)=IY+F1-F2                                                 F3286
      PTID(IP)=5                                                          F3287
      GO TO 530                                                           F3288
  487 PART(1,IP)=IX-F1-F2                                                 F3289
      PART(2,IP)=IY+F1+F2                                                 F3291
      PTID(IP)=6                                                          F3292
      GO TO 530                                                           F3293
  488 PART(1,IP)=IX-F1+F2                                                 F3294
      PART(2,IP)=IY+F1-F2                                                 F3295
      PTID(IP)=7                                                          F3296
      GO TO 530                                                           F3297
  489 PART(1,IP)=IX-F1+F2                                                 F3298
      PART(2,IP)=IY+F1+F2                                                 F3299
      PTID(IP)=8                                                          F3301
      GO TO 530                                                           F3302
  490 PART(1,IP)=IX+F1-F2                                                 F3303
      PART(2,IP)=IY-F1-F2                                                 F3304
      PTID(IP)=9                                                          F3305
      GO TO 530                                                           F3306
  491 PART(1,IP)=IX+F1-F2                                                 F3307
      PART(2,IP)=IY-F1+F2                                                 F3308
      PTID(IP)=10                                                         F3309
      GO TO 530                                                           F3311
  492 PART(1,IP)=IX+F1+F2                                                 F3312
      PART(2,IP)=IY-F1-F2                                                 F3313
      PTID(IP)=11                                                         F3314
      GO TO 530                                                           F3315
  493 PART(1,IP)=IX+F1+F2                                                 F3316
      PART(2,IP)=IY-F1+F2                                                 F3317
      PTID(IP)=12                                                         F3318
      GO TO 530                                                           F3319
  494 PART(1,IP)=IX+F1-F2                                                 F3321
      PART(2,IP)=IY+F1-F2                                                 F3322
      PTID(IP)=13                                                         F3323
      GO TO 530                                                           F3324
  495 PART(1,IP)=IX+F1-F2                                                 F3325
      PART(2,IP)=IY+F1+F2                                                 F3326
      PTID(IP)=14                                                         F3327
      GO TO 530                                                           F3328
  496 PART(1,IP)=IX+F1+F2                                                 F3329
      PART(2,IP)=IY+F1-F2                                                 F3331
      PTID(IP)=15                                                         F3332
      GO TO 530                                                           F3333
  497 PART(1,IP)=IX+F1+F2                                                 F3334
      PART(2,IP)=IY+F1+F2                                                 F3335
      PTID(IP)=16                                                         F3336
      GO TO 530                                                           F3337
  525 DLX=INX-IX                                                         SF3425
      PART(1,IP)=TEMPX-DLX                                               SF3435
      DLY=INY-IY                                                         SF3455
      PART(2,IP)=TEMPY-DLY                                               SF3465
      PART(2,IP)=-PART(2,IP)                                             SF3475
      PART(3,IP)=CONC(JX,JY)*DCYT                                        SF3478
      SUMC(JX,JY)=SUMC(JX,JY)+PART(3,IP)                                 SF3479
      NPCELL(JX,JY)=NPCELL(JX,JY)+1                                      SF3481
      GO TO 540                                                          SF3483
C                                                                         F3510
  530 PART(2,IP)=-PART(2,IP)                                              F3520
      PART(3,IP)=CONC(JX,JY)                                             SF3531
C     ****************************************************************    F3550
C           ---CHECK FOR DISCHARGE BOUNDARY AT NEW LOCATION---            F3560
  540 IFLAG=1                                                             F3571
  550 IF (NLOC.GT.0) GO TO 565                                           SF3575
      IF (VPRM(INX,INY).GT.0.09.AND.WT(INX,INY).LT.HK(INX,INY)) GO TO 56 SF3581
     10                                                                   F3590
C     IF (REC(INX,INY).GT.0.1) GO TO 560                                  FXXX
      IF (REC(INX,INY).GT.0.0) GO TO 560                                  F3600
      GO TO 590                                                           F3610
C     ****************************************************************    F3620
C           ---PUT PT. IN LIMBO IF PT. DENSITY NOT INCREASED---           F3625
  560 IF (NPOLD(JNX,JNY).LE.0) GO TO 590                                 SF3636
  565 CONTINUE                                                           SF3638
      PART(1,IN)=0.0                                                      F3645
      PART(2,IN)=0.0                                                      F3650
      PART(3,IN)=0.0                                                      F3660
      IF (NLOC.GT.0) GO TO 567                                           SF3661
      SUMC(JNX,JNY)=SUMC(JNX,JNY)-CONC(JNX,JNY)                          SF3663
      NPCELL(JNX,JNY)=NPCELL(JNX,JNY)-1                                  SF3665
      NPOLD(JNX,JNY)=NPOLD(JNX,JNY)-1                                    SF3667
  567 CONTINUE                                                           SF3668
      DO 570 ID=1,500                                                     F3670
      IF (LIMBO(ID).GT.0) GO TO 570                                       F3680
      LIMBO(ID)=IN                                                        F3690
      GO TO 590                                                           F3700
  570 CONTINUE                                                            F3710
C                                                                         F3720
  580 IF (IFLAG.LT.0) PART(2,IN)=-TEMPY                                   F3730
  590 CONTINUE                                                            F3750
C     ---END OF LOOP---                                                   F3760
C     ****************************************************************    F3762
C                                                                         F3763
C        ---INSERT TEMPORARY PARTICLES INTO LIMBO LOCATIONS---            F3764
      IF(NPTM.EQ.NP) GO TO 620                                            F3765
      IN=NPTM                                                             F3766
      DO 595 IL=1,500                                                     F3767
      IP=LIMBO(IL)                                                        F3768
      IF(IP.EQ.0) GO TO 595                                               F3769
      PART(1,IP)=PART(1,IN)                                               F3771
      PART(1,IN)=0.0                                                      F3772
      PART(2,IP)=PART(2,IN)                                               F3773
      PART(2,IN)=0.0                                                      F3774
      PART(3,IP)=PART(3,IN)                                               F3775
      PART(3,IN)=0.0                                                      F3776
      PTID(IP)=PTID(IN)                                                   F3777
      PTID(IN)=0                                                          F3778
      LIMBO(IL)=0                                                         F3779
      IN=IN-1                                                             F3781
      IF(IN.LE.NP) GO TO 596                                              F3782
  595 CONTINUE                                                            F3783
  596 NPTM=IN                                                             F3784
      GO TO 620                                                           F3785
C        ---RESTART MOVE IF PT. LIMIT EXCEEDED---                         F3790
  600 WRITE (6,700) IMOV,IN                                               F3800
      TEST=100.0                                                          F3810
      CALL GENPT                                                          F3820
      DO 610 IX=1,NMX                                                    SF3831
      DO 610 IY=1,NMY                                                    SF3841
      SUMC(IX,IY)=0.0                                                     F3850
  610 NPCELL(IX,IY)=0                                                     F3860
      TEST=0.0                                                            F3870
      GO TO 10                                                            F3880
C     ****************************************************************    F3890
  620 SUMTCH=SUMTCH+TIMV                                                  F3900
C        ---ADJUST NUMBER OF PARTICLES---                                 F3910
      NP=NPTM                                                             F3920
      WRITE (6,670) NP,IMOV                                               F3930
C     ****************************************************************    F3940
      CALL CNCON                                                          F3950
C     ****************************************************************    F3960
C        ---STORE OBS. WELL DATA FOR STEADY FLOW PROBLEMS---              F3970
      IF (S.GT.0.0) GO TO 640                                             F3980
      IF (NUMOBS.LE.0) GO TO 640                                          F3990
      J=MOD(IMOV,50)                                                      F4000
      IF (J.EQ.0) J=50                                                    F4010
      TMOBS(J)=SUMTCH                                                     F4020
      DO 630 I=1,NUMOBS                                                   F4030
      TMWL(I,J)=HK(IXOBS(I),IYOBS(I))                                     F4040
      IX=IXOBS(I)-MX+1                                                   SF4042
      IF (IXOBS(I).LT.MX.OR.IXOBS(I).GT.MMX) GO TO 630                   SF4043
      IY=IYOBS(I)-MY+1                                                   SF4044
      IF (IYOBS(I).LT.MY.OR.IYOBS(I).GT.MMY) GO TO 630                   SF4045
      TMCN(I,J)=CONC(IX,IY)                                              SF4051
  630 CONTINUE                                                           SF4053
C        ---PRINT CHEMICAL OUTPUT---                                      F4060
      IF (MOD(IMOV,50).EQ.0) IPRNT=1                                      F4065
  640 IF (IMOV.GE.NMOV) GO TO 660                                         F4070
      IF (MOD(IMOV,NPNTMV).EQ.0) IPRNT=-1                                 F4082
  650 IF (IPRNT.NE.0) CALL CHMOT                                          F4085
C     ****************************************************************    F4090
  660 RETURN                                                              F4100
C     ****************************************************************    F4110
C                                                                         F4120
C                                                                         F4130
C                                                                         F4140
  670 FORMAT (1H0,2X,2HNP,7X,2H= ,8X,I4,10X,11HIMOV     = ,8X,I4)         F4150
  680 FORMAT (1H0,10X,61HNO. OF PARTICLE MOVES REQUIRED TO COMPLETE THIS  F4160
     1 TIME STEP  = ,I4//)                                                F4170
  685 FORMAT (1H0,5X,51H*** CAUTION ***  DECAY HALF-LIFE IS LESS THAN TI  F4174R
     1MV;/23X,24HACCURACY MAY BE AFFECTED/23X,34H(REDUCE TIMV BY DECREAS  F4175R
     2ING CELDIS))                                                        F4176R
  690 FORMAT (1H0,5X,53H*** WARNING ***     QUADRANT NOT LOCATED FOR PT.  F4180
     1 NO. ,I5,11H , IN CELL ,2I4)                                        F4190
  700 FORMAT (1H0,5X,17H ***   NOTE   ***,10X,23HNPTM.EQ.NPMAX --- IMOV=  F4200
     1,I4,5X,8HPT. NO.=,I4,5X,10HCALL GENPT/)                             F4210
      END                                                                 F4220-
      SUBROUTINE CNCON                                                    G  10
      DOUBLE PRECISION TMRX,VPRM,HI,HR,HC,HK,WT,REC,RECH,TIM,AOPT,TITLE   G  20
      DOUBLE PRECISION XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,     G  30
     1                 HMIN,PYR                                           G  31
      DOUBLE PRECISION DXINV,DYINV,ARINV,PORINV                           G  35
      DOUBLE PRECISION FLW                                                G  40
      DOUBLE PRECISION DCYFCT,DCYT,DCYT2                                  G  42
      COMMON /PRMJ/ NTIM,NPMP,NPNT,NITP,N,NX,NY,NP,NREC,INT,NNX,NNY,NUMO  G  50
     1BS,NMOV,IMOV,NPMAX,ITMAX,NZCRIT,IPRNT,NPTPND,NPNTMV,NPNTVL,NPNTD,N  G  60
     2PNCHV,NPDELC,ICHK                                                   G  71
      COMMON /PRMC/ NODEID(40,40),NPCELL(20,20),NPOLD(20,20),LIMBO(500), SG  75
     1IXOBS(5),IYOBS(5)                                                   G  76
      COMMON /HEDA/ THCK(40,40),PERM(40,40),TMWL(5,50),TMOBS(50),ANFCTR  SG 100
      COMMON /HEDB/ TMRX(40,40,2),VPRM(40,40),HI(40,40),HR(40,40),HC(40, SG 110
     140),HK(40,40),WT(40,40),REC(40,40),RECH(40,40),TIM(100),AOPT(20),  SG 120
     2TITLE(10),XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,HMIN,PYR    G 130
      COMMON /HEDC/ MX,MY,MMX,MMY,NMX,NMY,MCHK                           SG 135
      COMMON /XINV/ DXINV,DYINV,ARINV,PORINV                              G 140
      COMMON /CHMA/ PART(3,6400),CONC(20,20),TMCN(5,50),VX(40,40),VY(40, SG 151
     140),CONINT(20,20),CNRECH(20,20),POROS,SUMTCH,BETA,TIMV,STORM,STORM SG 160
     2I,CMSIN,CMSOUT,FLMIN,FLMOT,SUMIO,CELDIS,DLTRAT,CSTORM               G 170
      COMMON /DIFUS/ DISP(20,20,4)                                        G 180
      COMMON /CHMC/ SUMC(20,20),VXBDY(40,40),VYBDY(40,40),PERM2(40,40)   SG 191
      COMMON /CHMR/ RF,DK,RHOB,THALF,DECAY,ADSORB,SORBI,DMASS1,CSTM2      G 195R
      DIMENSION CNCNC(20,20),CNOLD(20,20),CAVG(20,20)                     G 201R
C     ***************************************************************     G 210
      DO 10 IX=1,NMX                                                     SG 232
      DO 10 IY=1,NMY                                                     SG 242
      CNOLD(IX,IY)=CONC(IX,IY)                                            G 250
      CAVG(IX,IY)=CONC(IX,IY)                                             G 255R
   10 CNCNC(IX,IY)=0.0                                                    G 260
      APC=0.0                                                             G 270
      NZERO=0                                                             G 280
      TVA=AREA*TIMV                                                       G 290
      TVA2=TVA*0.5                                                        G 295R
      TMCHK=TIMV*10.0                                                     G 296R
      SRCDCY=0.0                                                          G 297R
      ARPOR=AREA*POROS                                                    G 300
      RFPOR=RF*POROS                                                      G 305R
      RFPORA=RFPOR*AREA                                                   G 306R
C     ***************************************************************     G 310
C     ---CONC. CHANGE DUE TO:                                             G 321R
C           MIXING AT SOURCE CELLS...                                     G 322R
C           ...WITH DECAY OF RECHARGE DURING TIME INCREMENT               G 323R
      CONST=TIMV                                                          G 341R
      DCYFCT=TIMV*DECAY                                                   G 342R
      RFFCT=CONST/RF                                                      G 344R
      DCYT=1.D0                                                           G 346A
      DCYT2=1.D0                                                          G 346B
      IF(DECAY.NE.0.0) THEN                                               G 346C
       DCYT=DEXP(-DCYFCT)                                                 G 346D
       DCYT2=DEXP(-DCYFCT*0.5D0)                                          G 346E
      END IF                                                              G 346F
      GO TO 70                                                            G 347R
   20 DO 60 IX=1,NMX                                                     SG 352
      DO 60 IY=1,NMY                                                     SG 362
      JX=IX+MX-1                                                         SG 366
      JY=IY+MY-1                                                         SG 367
      IF (THCK(JX,JY).EQ.0.0) GO TO 60                                   SG 371
      EQFCT1=RFFCT/THCK(JX,JY)                                           SG 382R
      EQFCT2=EQFCT1/POROS                                                 G 390
      C1=CAVG(IX,IY)                                                      G 401R
      IF (ABS(C1).LT.1.0E-20) C1=0.0                                     MULTICS
      C2=C1                                                               G 403R
      CLKCN=C1                                                            G 411R
      SLEAK=(HK(JX,JY)-WT(JX,JY))*VPRM(JX,JY)                            SG 421
      IF (SLEAK.GE.0.0) GO TO 25                                          G 425R
      CLKCN=CNRECH(IX,IY)*DCYT2                                           G 431A
      SRCDCY=SRCDCY+(CNRECH(IX,IY)-CLKCN)*SLEAK*TVA                       G 432R
   25 CONTINUE                                                            G 435R
      CNREC=C1                                                            G 450
      CNREC2=C1                                                           G 455
      RATE=REC(JX,JY)*ARINV                                              SG 461
      IF (RATE.GE.0.0) GO TO 27                                           G 471R
      CNREC=CNRECH(IX,IY)*DCYT2                                           G 472A
      SRCDCY=SRCDCY+(CNRECH(IX,IY)-CNREC)*RATE*TVA                        G 473R
   27 IF (RECH(JX,JY).GE.0.0) GO TO 29                                   SG 474R
      CNREC2=CNRECH(IX,IY)*DCYT2                                          G 476A
      SRCDCY=SRCDCY+(CNRECH(IX,IY)-CNREC2)*RECH(JX,JY)*TVA               SG 477R
   29 CONTINUE                                                            G 478R
      DIV=RATE+SLEAK+RECH(JX,JY)                                         SG 481
      IF (THALF.GT.TMCHK.OR.THALF.EQ.0.0) GO TO 37                        G 482R
      IF (DIV.GE.0.0) GO TO 37                                            G 484R
      IF (CNOLD(IX,IY).LE.0.0.OR.CONC(IX,IY).LE.0.0) GO TO 37             G 486R
      C1=EXP((ALOG(CNOLD(IX,IY))+ALOG(CONC(IX,IY)))*0.5)                  G 488R
C     IF (NPCELL(IX,IY).LE.0) C1=CNOLD(IX,IY)*DCYT2                       G 489A
   37 CONTINUE                                                            G 592R
      DELC=EQFCT2*(C1*DIV-RATE*CNREC-SLEAK*CLKCN-RECH(JX,JY)*CNREC2)     SG 596
   40 CNCNC(IX,IY)=CNCNC(IX,IY)+DELC                                      G 610
C     ---CONC. CHANGE DUE TO DISPERSION FOR TIMV---                       G 621R
C        ---DISPERSION WITH TENSOR COEFFICIENTS---                        G 630
      IF (BETA.EQ.0.0) GO TO 60                                           G 635
      IF (MCHK.EQ.1) GO TO 39                                            SG 638
      X1=DISP(IX,IY,1)*(CAVG(IX+1,IY)-C2)                                 G 651R
      X2=DISP(IX-1,IY,1)*(CAVG(IX-1,IY)-C2)                               G 661R
      Y1=DISP(IX,IY,2)*(CAVG(IX,IY+1)-C2)                                 G 671R
      Y2=DISP(IX,IY-1,2)*(CAVG(IX,IY-1)-C2)                               G 681R
      XX1=DISP(IX,IY,3)*(CAVG(IX,IY+1)+CAVG(IX+1,IY+1)-CAVG(IX,IY-1)-CAV  G 691R
     1G(IX+1,IY-1))                                                       G 701R
      XX2=DISP(IX-1,IY,3)*(CAVG(IX,IY+1)+CAVG(IX-1,IY+1)-CAVG(IX,IY-1)-C  G 711R
     1AVG(IX-1,IY-1))                                                     G 721R
      YY1=DISP(IX,IY,4)*(CAVG(IX+1,IY)+CAVG(IX+1,IY+1)-CAVG(IX-1,IY)-CAV  G 731R
     1G(IX-1,IY+1))                                                       G 741R
      YY2=DISP(IX,IY-1,4)*(CAVG(IX+1,IY)+CAVG(IX+1,IY-1)-CAVG(IX-1,IY)-C  G 751R
     1AVG(IX-1,IY-1))                                                     G 761R
      GO TO 48                                                           SG 762
   39 X1=0.0                                                             SG 763A
      X2=0.0                                                             SG 763B
      Y1=0.0                                                             SG 763C
      Y2=0.0                                                             SG 763D
      XX1=0.0                                                            SG 763E
      XX2=0.0                                                            SG 763F
      YY1=0.0                                                            SG 763G
      YY2=0.0                                                            SG 763H
      IF (IX.GE.NMX) GO TO 41                                            SG 763I
      X1=DISP(IX,IY,1)*(CAVG(IX+1,IY)-C2)                                SG 763J
   41 IF (IX.LE.1) GO TO 42                                              SG 763K
      X2=DISP(IX-1,IY,1)*(CAVG(IX-1,IY)-C2)                              SG 763L
   42 IF (IY.GE.NMY) GO TO 43                                            SG 763M
      Y1=DISP(IX,IY,2)*(CAVG(IX,IY+1)-C2)                                SG 763N
   43 IF (IY.LE.1) GO TO 46                                              SG 763O
      Y2=DISP(IX,IY-1,2)*(CAVG(IX,IY-1)-C2)                              SG 763P
   44 IF (IY.GE.NMY) GO TO 47                                            SG 764
      IF (IX.GE.NMX) GO TO 45                                            SG 765
      XX1=DISP(IX,IY,3)*(CAVG(IX,IY+1)+CAVG(IX+1,IY+1)-CAVG(IX,IY-1)-CAV SG 766A
     1G(IX+1,IY-1))                                                      SG 766B
   45 IF (IX.LE.1) GO TO 48                                              SG 766C
      XX2=DISP(IX-1,IY,3)*(CAVG(IX,IY+1)+CAVG(IX-1,IY+1)-CAVG(IX,IY-1)-C SG 766D
     1AVG(IX-1,IY-1))                                                    SG 766E
   46 IF (IX.LE.1) GO TO 48                                              SG 766F
      IF (IY.GE.NMY) GO TO 47                                            SG 766G
      IF (IX.GE.NMX) GO TO 48                                            SG 766H
      YY1=DISP(IX,IY,4)*(CAVG(IX+1,IY)+CAVG(IX+1,IY+1)-CAVG(IX-1,IY)-CAV SG 766I
     1G(IX-1,IY+1))                                                      SG 766J
   47 IF (IY.LE.1) GO TO 48                                              SG 766K
      IF (IX.LE.1.OR.IX.GE.NMX) GO TO 48                                 SG 766L
      YY2=DISP(IX,IY-1,4)*(CAVG(IX+1,IY)+CAVG(IX+1,IY-1)-CAVG(IX-1,IY)-C SG 766M
     1AVG(IX-1,IY-1))                                                    SG 766O
   48 CONTINUE                                                           SG 768
      CNCNC(IX,IY)=CNCNC(IX,IY)+EQFCT1*(X1+X2+Y1+Y2+XX1-XX2+YY1-YY2)      G 771
   60 CONTINUE                                                            G 780
C     ****************************************************************    G 790
      GO TO 110                                                           G 820
C     ****************************************************************    G 830
C     ---CONC. CHANGE AT NODES DUE TO CONVECTION---                       G 840
C                                  ... AND DECAY---                       G 841R
   70 DO 90 IX=1,NMX                                                     SG 852
      DO 90 IY=1,NMY                                                     SG 862
      JX=IX+MX-1                                                         SG 866
      JY=IY+MY-1                                                         SG 867
      IF (THCK(JX,JY).EQ.0.0) GO TO 90                                   SG 871
      APC=NPCELL(IX,IY)                                                   G 880
      IF (APC.GT.0.0) GO TO 80                                            G 890
      IF (REC(JX,JY).NE.0.0.OR.VPRM(JX,JY).GT.0.09) GO TO 90             SG 901
      NZERO=NZERO+1                                                       G 910
      GO TO 90                                                            G 920
   80 CONC(IX,IY)=SUMC(IX,IY)/APC                                         G 930
      CAVG(IX,IY)=(CONC(IX,IY)+CNOLD(IX,IY))*0.5                          G 935R
   90 CONTINUE                                                            G 940
C        ---CHECK NUMBER OF CELLS VOID OF PTS.---                         G 950
      IF (NZERO.GT.0) WRITE (6,290) NZERO,IMOV                            G 960
      IF (NZERO.LE.NZCRIT) GO TO 20                                       G 970
      TEST=99.0                                                           G 980
      WRITE (6,300)                                                       G 990
      WRITE (6,320)                                                       G1000
      DO 100 IY=1,NMY                                                    SG1011
  100 WRITE (6,330) (NPCELL(IX,IY),IX=1,NMX)                             SG1021
      GO TO 20                                                            G1030
C     ****************************************************************    G1040
C     ---CHANGE CONCENTRATIONS AT NODES---                                G1050
  110 DO 130 IX=1,NMX                                                    SG1062
      DO 130 IY=1,NMY                                                    SG1072
      JX=IX+MX-1                                                         SG1076
      JY=IY+MY-1                                                         SG1077
      IF (THCK(JX,JY).EQ.0.0) GO TO 120                                  SG1081
      CNCPCT=0.0                                                          G1084G
      IF (CONC(IX,IY).GT.0.0) CNCPCT=CNCNC(IX,IY)/CONC(IX,IY)             G1085G
      CONC(IX,IY)=CONC(IX,IY)+CNCNC(IX,IY)                                G1090
      SUMC(IX,IY)=0.0                                                     G1110
      IF (CNCPCT.LT.0.0) SUMC(IX,IY)=CNCPCT                               G1125G
      GO TO 130                                                           G1150
  120 IF (CONC(IX,IY).GT.0.0) WRITE (6,310) JX,JY,CONC(IX,IY)            SG1161
      CONC(IX,IY)=0.0                                                     G1170
  130 CONTINUE                                                            G1180
C     ****************************************************************    G1190
C     ---CHANGE CONCENTRATION OF PARTICLES---                             G1200
      DO 180 IN=1,NP                                                      G1210
      IF (PART(1,IN).EQ.0.0) GO TO 180                                    G1220
      INX=ABS(PART(1,IN))+0.5                                             G1230
      INY=ABS(PART(2,IN))+0.5                                             G1240
      JNX=INX-MX+1                                                       SG1244
      JNY=INY-MY+1                                                       SG1245
C     ---UPDATE CONC. OF PTS. IN SINK/SOURCE CELLS---                     G1250
      IF (REC(INX,INY).NE.0.0) GO TO 140                                  G1260
      IF (VPRM(INX,INY).LE.0.09) GO TO 150                                G1270
  140 PART(3,IN)=CONC(JNX,JNY)                                           SG1281
      GO TO 180                                                           G1290
  150 IF (CNCNC(JNX,JNY).LT.0.0) GO TO 170                               SG1301
  160 PART(3,IN)=PART(3,IN)+CNCNC(JNX,JNY)                               SG1311
      GO TO 180                                                           G1320
  170 IF (CONC(JNX,JNY).LE.0.0) GO TO 160                                SG1331
      IF (SUMC(JNX,JNY).LT.-1.0) GO TO 160                               SG1341
      PART(3,IN)=PART(3,IN)+PART(3,IN)*SUMC(JNX,JNY)                     SG1351
  180 CONTINUE                                                            G1360
      WRITE (6,280) TIM(N),TIMV,SUMTCH                                    G1370
C     ****************************************************************    G1380
C     ---COMPUTE MASS BALANCE FOR SOLUTE---                               G1390
      CSTORM=0.0                                                          G1400
      STORM=0.0                                                           G1410
      ADSORB=0.0                                                          G1411R
      DO 270 IX=1,NMX                                                    SG1422
      DO 270 IY=1,NMY                                                    SG1432
      JX=IX+MX-1                                                         SG1435
      JY=IY+MY-1                                                         SG1436
      IF (THCK(JX,JY).EQ.0.0) GO TO 270                                  SG1441
      C1=CONC(IX,IY)                                                      G1442R
      C1B=C1*THCK(JX,JY)                                                 SG1443R
      DELDCY=CNOLD(IX,IY)-CNOLD(IX,IY)*DCYT                               G1444A
      DMASS1=DMASS1-DELDCY*THCK(JX,JY)*RFPORA                            SG1446R
      SUMC(IX,IY)=0.0                                                     G1450
      WTFCT=0.0                                                           G1455
C        ---COMPUTE MASS OF SOLUTE IN STORAGE---                          G1460
      STORM=STORM+C1B*ARPOR                                               G1471R
C        ---COMPUTE MASS ADSORBED---                                      G1472R
      ADSORB=C1B*DK*RHOB*AREA+ADSORB                                      G1473R
C        ---ACCOUNT FOR MASS PUMPED IN, OUT, RECHARGED, & DISCHARGED---   G1480
      IF (REC(JX,JY)) 200,210,190                                        SG1491
  190 FCT1=NPOLD(IX,IY)                                                   G1492
      FCT2=NPCELL(IX,IY)                                                  G1494
      IF (FCT2.GT.0.0) WTFCT=FCT1/FCT2                                    G1496
      CMSOUT=CMSOUT+REC(JX,JY)*TIMV*((1.0-WTFCT)*CNOLD(IX,IY)+WTFCT*C1)  SG1498R
      GO TO 210                                                           G1510
  200 CMSIN=CMSIN+REC(JX,JY)*CNRECH(IX,IY)*TIMV                          SG1521
  210 IF (RECH(JX,JY)) 230,240,220                                       SG1531
  220 CMSOUT=CMSOUT+RECH(JX,JY)*TVA2*(CNOLD(IX,IY)+C1)                   SG1538R
      GO TO 240                                                           G1550
  230 CMSIN=CMSIN+RECH(JX,JY)*CNRECH(IX,IY)*TVA                          SG1561
C     ****************************************************************    G1570
C        ---ACCOUNT FOR BOUNDARY FLOW---                                  G1580
  240 IF (VPRM(JX,JY).EQ.0.0) GO TO 265                                  SG1586
      FLW=VPRM(JX,JY)*(WT(JX,JY)-HK(JX,JY))                              SG1601
      IF (FLW.GT.0.0) GO TO 250                                           G1610
      IF (FLW.LT.0.0) GO TO 260                                           G1620
      GO TO 265                                                           G1625
C           ---MASS IN BOUNDARY DURING TIME STEP---                       G1640
  250 FLMIN=FLMIN+FLW*CNRECH(IX,IY)*TVA                                   G1650
      GO TO 265                                                           G1655
C           ---MASS OUT DURING TIME STEP---                               G1670
  260 FCT1=NPOLD(IX,IY)                                                   G1672
      FCT2=NPCELL(IX,IY)                                                  G1674
      IF (FCT2.GT.0.0) WTFCT=FCT1/FCT2                                    G1676
      FLMOT=FLMOT+FLW*TVA*((1.0-WTFCT)*CNOLD(IX,IY)+WTFCT*C1)             G1679R
  265 NPOLD(IX,IY)=NPCELL(IX,IY)                                          G1682
      NPCELL(IX,IY)=0                                                     G1684
  270 CONTINUE                                                            G1690
C           ---SUBGRID BOUNDARIES---                                     SG1692A
      IF (MCHK.EQ.0) GO TO 275                                           SG1692B
      YT=YDEL*TIMV                                                       SG1692C
      XT=XDEL*TIMV                                                       SG1692D
      DO 272 IY=MY,MMY                                                   SG1694A
      IX=MX                                                              SG1694B
      JY=IY-MY+1                                                         SG1694C
      FLW=TMRX(IX-1,IY,1)*(HK(IX-1,IY)-HK(IX,IY))                        SG1694D
      IF (FLW.GT.0.0) FLMIN=FLMIN+FLW*YT*CNOLD(1,JY)*DCYT                SG1694G
      IF (FLW.LT.0.0) FLMOT=FLMOT+FLW*YT*CNOLD(1,JY)*DCYT                SG1694H
      IX=MMX                                                             SG1695A
      FLW=TMRX(IX,IY,1)*(HK(IX+1,IY)-HK(IX,IY))                          SG1695B
      IF (FLW.GT.0.0) FLMIN=FLMIN+FLW*YT*CNOLD(NMX,JY)*DCYT              SG1695E
  272 IF (FLW.LT.0.0) FLMOT=FLMOT+FLW*YT*CNOLD(NMX,JY)*DCYT              SG1695F
      DO 274 IX=MX,MMX                                                   SG1696A
      IY=MY                                                              SG1696B
      JX=IX-MX+1                                                         SG1696C
      FLW=TMRX(IX,IY-1,2)*(HK(IX,IY-1)-HK(IX,IY))                        SG1696D
      IF (FLW.GT.0.0) FLMIN=FLMIN+FLW*XT*CNOLD(JX,1)*DCYT                SG1696G
      IF (FLW.LT.0.0) FLMOT=FLMOT+FLW*XT*CNOLD(JX,1)*DCYT                SG1696H
      IY=MMY                                                             SG1697A
      FLW=TMRX(IX,IY,2)*(HK(IX,IY+1)-HK(IX,IY))                          SG1697B
      IF (FLW.GT.0.0) FLMIN=FLMIN+FLW*XT*CNOLD(JX,NMY)*DCYT              SG1697E
  274 IF (FLW.LT.0.0) FLMOT=FLMOT+FLW*XT*CNOLD(JX,NMY)*DCYT              SG1697F
  275 CONTINUE                                                           SG1698A
C     ****************************************************************    G1700
C        ---COMPUTE CHANGE IN MASS OF SOLUTE STORED---                    G1710
C        ---COMPUTE MASS LOST BY DECAY---                                 G1715R
      CSTORM=STORM-STORMI                                                 G1720
      CSTM2=CSTORM+ADSORB-SORBI                                           G1725R
      DMASS1=DMASS1+SRCDCY                                                G1726R
      SUMIO=FLMIN+FLMOT-CMSIN-CMSOUT                                      G1730
C     ****************************************************************    G1740
C     ---REGENERATE PARTICLES IF 'NZCRIT' EXCEEDED---                     G1750
      IF (TEST.GT.98.0) CALL GENPT                                        G1760
      TEST=0.0                                                            G1770
C     ****************************************************************    G1780
      RETURN                                                              G1790
C     ****************************************************************    G1800
C                                                                         G1810
C                                                                         G1820
C                                                                         G1830
  280 FORMAT (3H   ,11HTIM(N)   = ,1G12.5,10X,11HTIMV     = ,1G12.5,10X,  G1840
     19HSUMTCH = ,G12.5)                                                  G1850
  290 FORMAT (1H0,5X,40HNUMBER OF CELLS WITH ZERO PARTICLES  =  ,I4,5X,9  G1860
     1HIMOV  =  ,I4/)                                                     G1870
  300 FORMAT (1H0,5X,44H***   NZCRIT EXCEEDED  ---  CALL GENPT   ***/)    G1880
  310 FORMAT (1H ,5X,37H***CONC.GT.0.AND.THCK.EQ.0 AT NODE = ,2I4,4X,7HC  G1890
     1ONC = ,G10.4,4H ***)                                                G1900
  320 FORMAT (1H0,2X,6HNPCELL/)                                           G1910
  330 FORMAT (1H ,4X,20I3)                                                G1920
      END                                                                 G1930-
      SUBROUTINE OUTPT                                                    H  10
      DOUBLE PRECISION TMRX,VPRM,HI,HR,HC,HK,WT,REC,RECH,TIM,AOPT,TITLE   H  20
      DOUBLE PRECISION XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,     H  30
     1                 HMIN,PYR                                           H  31
      DOUBLE PRECISION TOTLQ,TOTLQI,TPIN,TPOUT,PQIN,PQOUT,DELQ,QOUT,QIN,  H  32
     1                 QNET,DDRW                                          H  32A
      DOUBLE PRECISION QSTR,TPUM,PUMP,TOTLQN,SRCS,SINKS,ERRMB,DENOM,      H  33
     1                 PCTERR                                             H  34
      COMMON /PRMJ/ NTIM,NPMP,NPNT,NITP,N,NX,NY,NP,NREC,INT,NNX,NNY,NUMO  H  40
     1BS,NMOV,IMOV,NPMAX,ITMAX,NZCRIT,IPRNT,NPTPND,NPNTMV,NPNTVL,NPNTD,N  H  50
     2PNCHV,NPDELC,ICHK                                                   H  61
      COMMON /HEDA/ THCK(40,40),PERM(40,40),TMWL(5,50),TMOBS(50),ANFCTR  SH  90
      COMMON /HEDB/ TMRX(40,40,2),VPRM(40,40),HI(40,40),HR(40,40),HC(40, SH 100
     140),HK(40,40),WT(40,40),REC(40,40),RECH(40,40),TIM(100),AOPT(20),  SH 110
     2TITLE(10),XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,HMIN,PYR    H 120
      COMMON /BALM/ TOTLQ,TOTLQI,TPIN,TPOUT                               H 126
      DIMENSION IH(40)                                                   SH 141
C     ***************************************************************     H 150
      TIMD=SUMT/86400.                                                    H 160
      TIMY=SUMT/(86400.0*365.25)                                          H 170
C     ---PRINT HEAD VALUES---                                             H 180
      WRITE (6,120)                                                       H 190
      WRITE (6,130) N                                                     H 200
      WRITE (6,140) SUMT                                                  H 210
      WRITE (6,150) TIMD                                                  H 220
      WRITE (6,160) TIMY                                                  H 230
      WRITE (6,170)                                                       H 240
      DO 10 IY=1,NY                                                       H 250
   10 WRITE (6,180) (HK(IX,IY),IX=1,NX)                                   H 260
      IF (N.EQ.0) GO TO 110                                               H 270
C     ***************************************************************     H 280
C     ---PRINT HEAD MAP---                                                H 290
      WRITE (6,120)                                                       H 300
      WRITE (6,130) N                                                     H 310
      WRITE (6,140) SUMT                                                  H 320
      WRITE (6,150) TIMD                                                  H 330
      WRITE (6,160) TIMY                                                  H 340
      WRITE (6,170)                                                       H 350
      DO 30 IY=1,NY                                                       H 360
      DO 20 IX=1,NX                                                       H 370
   20 IH(IX)=HK(IX,IY)+0.5                                                H 380
   30 WRITE (6,190) (IH(ID),ID=1,NX)                                      H 390
C     ***************************************************************     H 400
C     ---COMPUTE WATER BALANCE AND DRAWDOWN---                            H 410
      QSTR=0.0                                                            H 420
      PUMP=0.0                                                            H 430
      PQIN=0.0                                                            H 432
      PQOUT=0.0                                                           H 434
      TPUM=0.0                                                            H 440
      QIN=0.0                                                             H 450
      QOUT=0.0                                                            H 460
      QNET=0.0                                                            H 470
      DELQ=0.0                                                            H 480
      PCTERR=0.0                                                          H 500
      WRITE (6,290)                                                       H 510
C                                                                         H 520
      DO 80 IY=1,NY                                                       H 530
      DO 70 IX=1,NX                                                       H 540
      IH(IX)=0.0                                                          H 550
      IF (THCK(IX,IY).EQ.0.0) GO TO 70                                    H 560
      IF (REC(IX,IY).GT.0.0) GO TO 32                                     H 562
      PQIN=PQIN+REC(IX,IY)                                                H 564
      GO TO 34                                                            H 566
   32 PQOUT=PQOUT+REC(IX,IY)                                              H 568
   34 IF (RECH(IX,IY).GT.0.0) GO TO 36                                    H 572
      PQIN=PQIN+RECH(IX,IY)*AREA                                          H 574
      GO TO 38                                                            H 576
   36 PQOUT=PQOUT+RECH(IX,IY)*AREA                                        H 578
   38 IF (VPRM(IX,IY).EQ.0.0) GO TO 60                                    H 582
      DELQ=VPRM(IX,IY)*AREA*(WT(IX,IY)-HK(IX,IY))                         H 590
      IF (DELQ.GT.0.0) GO TO 40                                           H 600
      QOUT=QOUT+DELQ                                                      H 610
      GO TO 50                                                            H 620
   40 QIN=QIN+DELQ                                                        H 630
   50 QNET=QNET+DELQ                                                      H 640
   60 DDRW=HI(IX,IY)-HK(IX,IY)                                            H 650
      IH(IX)=DDRW+0.5                                                     H 660
      QSTR=QSTR+DDRW*AREA*S                                               H 670
   70 CONTINUE                                                            H 680
C        ---PRINT DRAWDOWN MAP---                                         H 690
      WRITE (6,300) (IH(IX),IX=1,NX)                                      H 700
   80 CONTINUE                                                            H 710
      TPUM=PQIN+PQOUT                                                     H 716
      PUMP=TPIN+TPOUT                                                     H 721
      TOTLQN=TOTLQ+TOTLQI                                                 H 745
      SRCS=QSTR-TPIN+TOTLQI                                               H 755
      SINKS=TPOUT-TOTLQ                                                   H 765
      ERRMB=SRCS-SINKS                                                    H 775
      DENOM=(SRCS+SINKS)*0.5                                              H 785
      IF (DENOM.EQ.0.0) GO TO 100                                         H 795
      PCTERR=ERRMB*100.0/DENOM                                            H 805
C        ---PRINT MASS BALANCE DATA FOR FLOW MODEL---                     H 830
  100 WRITE (6,240)                                                       H 840
      WRITE (6,211) TPIN                                                  H 842
      WRITE (6,212) TPOUT                                                 H 844
      WRITE (6,250) PUMP                                                  H 850
      WRITE (6,230) QSTR                                                  H 860
      WRITE (6,202) TOTLQI                                                H 862
      WRITE (6,203) TOTLQ                                                 H 864
      WRITE (6,260) TOTLQN                                                H 866
      WRITE (6,270) ERRMB                                                 H 880
      WRITE (6,280) PCTERR                                                H 883
      WRITE (6,201)                                                       H 886
      WRITE (6,202) QIN                                                   H 889
      WRITE (6,203) QOUT                                                  H 892
      WRITE (6,204) QNET                                                  H 895
      WRITE (6,211) PQIN                                                  H 898
      WRITE (6,212) PQOUT                                                 H 901
      WRITE (6,210) TPUM                                                  H 910
C     ***************************************************************     H 940
  110 RETURN                                                              H 950
C     ***************************************************************     H 960
C                                                                         H 970
C                                                                         H 980
C                                                                         H 990
  120 FORMAT (1H1,23HHEAD DISTRIBUTION - ROW)                             H1000
  130 FORMAT (1X,23HNUMBER OF TIME STEPS = ,1I5)                          H1010
  140 FORMAT (8X,16HTIME(SECONDS) = ,1G12.5)                              H1020
  150 FORMAT (8X,16HTIME(DAYS)    = ,1E12.5)                              H1030
  160 FORMAT (8X,16HTIME(YEARS)   = ,1E12.5)                              H1040
  170 FORMAT (1H )                                                        H1050
  180 FORMAT (1H0,10F12.7)                                                H1055
  190 FORMAT (1H0,30I4)                                                  SH1071
  201 FORMAT (1H0,2X,33HRATE MASS BALANCE -- (IN C.F.S.) //)              H1073
  202 FORMAT (4X,29HLEAKAGE INTO AQUIFER       = ,E12.5)                  H1076
  203 FORMAT (4X,29HLEAKAGE OUT OF AQUIFER     = ,E12.5)                  H1083
  204 FORMAT (4X,29HNET LEAKAGE    (QNET)      = ,E12.5)                  H1086
  210 FORMAT (4X,29HNET WITHDRAWAL   (TPUM)    = ,E12.5)                  H1093
  211 FORMAT (4X,29HRECHARGE AND INJECTION     = ,E12.5)                  H1096
  212 FORMAT (4X,29HPUMPAGE AND E-T WITHDRAWAL = ,E12.5)                  H1103
  230 FORMAT (4X,29HWATER RELEASE FROM STORAGE = ,1E12.5)                 H1120
  240 FORMAT (1H0,2X,38HCUMULATIVE MASS BALANCE -- (IN FT**3) //)         H1125
  250 FORMAT (4X,29HCUMULATIVE  NET  PUMPAGE   = ,1E12.5)                 H1140
  260 FORMAT (4X,29HCUMULATIVE  NET  LEAKAGE   = ,1E12.5)                 H1150
  270 FORMAT (1H0,7X,25HMASS BALANCE RESIDUAL  = ,G12.5)                  H1160
  280 FORMAT (1H ,7X,25HERROR  (AS PERCENT)    = ,G12.5/)                 H1170
  290 FORMAT (1H1,8HDRAWDOWN)                                             H1180
  300 FORMAT (3H   ,25I5)                                                SH1191
      END                                                                 H1200-
      SUBROUTINE CHMOT                                                    I  10
      DOUBLE PRECISION TMRX,VPRM,HI,HR,HC,HK,WT,REC,RECH,TIM,AOPT,TITLE   I  20
      DOUBLE PRECISION XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,     I  30
     1                 HMIN,PYR                                           I  31
      COMMON /PRMJ/ NTIM,NPMP,NPNT,NITP,N,NX,NY,NP,NREC,INT,NNX,NNY,NUMO  I  40
     1BS,NMOV,IMOV,NPMAX,ITMAX,NZCRIT,IPRNT,NPTPND,NPNTMV,NPNTVL,NPNTD,N  I  50
     2PNCHV,NPDELC,ICHK                                                   I  61
      COMMON /PRMC/ NODEID(40,40),NPCELL(20,20),NPOLD(20,20),LIMBO(500), SI  65
     1IXOBS(5),IYOBS(5)                                                   I  66
      COMMON /HEDA/ THCK(40,40),PERM(40,40),TMWL(5,50),TMOBS(50),ANFCTR  SI  90
      COMMON /HEDB/ TMRX(40,40,2),VPRM(40,40),HI(40,40),HR(40,40),HC(40, SI 100
     140),HK(40,40),WT(40,40),REC(40,40),RECH(40,40),TIM(100),AOPT(20),  SI 110
     2TITLE(10),XDEL,YDEL,S,AREA,SUMT,RHO,PARAM,TEST,TOL,PINT,HMIN,PYR    I 120
      COMMON /HEDC/ MX,MY,MMX,MMY,NMX,NMY,MCHK                           SI 125
      COMMON /CHMA/ PART(3,6400),CONC(20,20),TMCN(5,50),VX(40,40),VY(40, SI 131
     140),CONINT(20,20),CNRECH(20,20),POROS,SUMTCH,BETA,TIMV,STORM,STORM SI 140
     2I,CMSIN,CMSOUT,FLMIN,FLMOT,SUMIO,CELDIS,DLTRAT,CSTORM               I 150
      COMMON /CHMR/ RF,DK,RHOB,THALF,DECAY,ADSORB,SORBI,DMASS1,CSTM2      I 155R
      DIMENSION IC(20)                                                    I 160
C     ***************************************************************     I 170
      TMFY=86400.0*365.25                                                 I 180
      TMYR=SUMT/TMFY                                                      I 190
      TCHD=SUMTCH/86400.0                                                 I 200
      TCHYR=SUMTCH/TMFY                                                   I 210
      ERR1=0.0                                                            I 212
      ERR3=0.0                                                            I 214
      IF (IPRNT.GT.0) GO TO 100                                           I 220
C     ***************************************************************     I 230
C     ---PRINT CONCENTRATIONS---                                          I 240
      WRITE (6,160)                                                       I 250
      WRITE (6,170) N                                                     I 260
      IF (N.GT.0) WRITE (6,180) TIM(N)                                    I 270
      WRITE (6,190) SUMT                                                  I 280
      WRITE (6,450) SUMTCH                                                I 290
      WRITE (6,200) TCHD                                                  I 300
      WRITE (6,210) TMYR                                                  I 310
      WRITE (6,460) TCHYR                                                 I 320
      WRITE (6,380) IMOV                                                  I 330
      WRITE (6,220)                                                       I 340
      DO 20 IY=1,NMY                                                     SI 351
      DO 10 IX=1,NMX                                                     SI 361
   10 IC(IX)=CONC(IX,IY)+0.5                                              I 370
   20 WRITE (6,240) (IC(IX),IX=1,NMX)                                    SI 381
C     ***************************************************************     I 390
      IF (N.EQ.0) GO TO 150                                               I 400
      IF (NPDELC.EQ.0) GO TO 50                                           I 410
C                                                                         I 420
C     ---PRINT CHANGES IN CONCENTRATION---                                I 430
      WRITE (6,230)                                                       I 440
      WRITE (6,170) N                                                     I 450
      WRITE (6,180) TIM(N)                                                I 460
      WRITE (6,190) SUMT                                                  I 470
      WRITE (6,450) SUMTCH                                                I 480
      WRITE (6,200) TCHD                                                  I 490
      WRITE (6,210) TMYR                                                  I 500
      WRITE (6,460) TCHYR                                                 I 510
      WRITE (6,380) IMOV                                                  I 520
      WRITE (6,220)                                                       I 530
      DO 40 IY=1,NMY                                                     SI 541
      DO 30 IX=1,NMX                                                     SI 551
      CNG=CONC(IX,IY)-CONINT(IX,IY)                                       I 560
   30 IC(IX)=CNG                                                          I 570
   40 WRITE (6,240) (IC(IX),IX=1,NMX)                                    SI 581
C     ***************************************************************     I 590
C     ---PRINT MASS BALANCE DATA FOR SOLUTE---                            I 600
   50 RESID=SUMIO-CSTM2+DMASS1                                            I 611R
      SUMIN=FLMIN-CMSIN                                                   I 615
      IF (SUMIN.EQ.0.0) GO TO 60                                          I 625
      ERR1=RESID*100.0/SUMIN                                              I 635
      IF (SUMIO.GT.(STORMI+SORBI)) GO TO 70                               I 637R
   60 IF (STORMI.EQ.0.0) GO TO 70                                         I 650
      ERR3=-100.0*(RESID)/(STORMI+SORBI-SUMIO)                            I 661R
   70 WRITE (6,220)                                                       I 670
      WRITE (6,250)                                                       I 680
      WRITE (6,220)                                                       I 690
      WRITE (6,260) FLMIN                                                 I 700
      WRITE (6,270) FLMOT                                                 I 710
      RECIN=-CMSIN                                                        I 720
      RECOUT=-CMSOUT                                                      I 730
      WRITE (6,290) RECIN                                                 I 740
      WRITE (6,280) RECOUT                                                I 750
      WRITE (6,295) DMASS1                                                I 751R
      WRITE (6,296) ADSORB                                                I 752R
      WRITE (6,298) SORBI                                                 I 753R
      WRITE (6,300) SUMIO                                                 I 760
      WRITE (6,310) STORMI                                                I 770
      WRITE (6,320) STORM                                                 I 780
      WRITE (6,330) CSTORM                                                I 790
      WRITE (6,332) CSTM2                                                 I 795R
      WRITE (6,340)                                                       I 810
      WRITE (6,350) RESID                                                 I 820
      WRITE (6,360) ERR1                                                  I 830
   80 IF (STORMI.EQ.0.0) GO TO 90                                         I 840
      IF (SUMIN.NE.0.0.AND.SUMIO.GT.(STORMI+SORBI)) GO TO 90              I 845R
      WRITE (6,370)                                                       I 850
      WRITE (6,360) ERR3                                                  I 860
C     ***************************************************************     I 870
C     ---PRINT HYDROGRAPHS AFTER 50 STEPS OR END OF SIMULATION---         I 880
   90 IF (MOD(IMOV,50).EQ.0.AND.S.EQ.0.0) GO TO 100                       I 890
      IF (MOD(N,50).EQ.0.AND.S.GT.0.0) GO TO 100                          I 900
      IF (S.EQ.0.0.AND.N.LT.NTIM.AND.INT.GT.0) GO TO 100                  I 905
      GO TO 150                                                           I 910
  100 WRITE (6,390) TITLE                                                 I 920
      IF (NUMOBS.LE.0) GO TO 150                                          I 930
      WRITE (6,400) INT                                                   I 940
      IF (S.GT.0.0) WRITE (6,410)                                         I 950
      IF (S.EQ.0.0) WRITE (6,420)                                         I 960
C        ---TABULATE HYDROGRAPH DATA---                                   I 970
      MOZ=0                                                               I 980
      IF (S.GT.0.0) GO TO 110                                             I 990
      NTO=NMOV                                                            I1000
      IF (NMOV.GT.50) NTO=MOD(IMOV,50)                                    I1010
      GO TO 120                                                           I1020
  110 NTO=N                                                               I1031
      IF (N.GT.50) NTO=MOD(N,50)                                          I1041
  120 IF (NTO.EQ.0) NTO=50                                                I1050
      DO 140 J=1,NUMOBS                                                   I1060
      TMYR=0.0                                                            I1070
      JX=IXOBS(J)                                                        SI1074
      JY=IYOBS(J)                                                        SI1076
      WRITE (6,430) J,IXOBS(J),IYOBS(J)                                   I1080
      C1INT=0.0                                                          SI1082
      IX=JX-MX+1                                                         SI1083
      IY=JY-MY+1                                                         SI1084
      IF (JX.LT.MX.OR.JY.LT.MY.OR.JX.GT.MMX.OR.JY.GT.MMY) GO TO 125      SI1085
      C1INT=CONINT(IX,IY)                                                SI1086
      GO TO 127                                                          SI1087
  125 WRITE (6,435)                                                      SI1088
  127 WRITE (6,440) MOZ,WT(JX,JY),C1INT,                                 SI1091
     1TMYR                                                                I1100
      DO 130 M=1,NTO                                                      I1110
      TMYR=TMOBS(M)/TMFY                                                  I1120
  130 WRITE (6,440) M,TMWL(J,M),TMCN(J,M),TMYR                            I1130
  140 CONTINUE                                                            I1140
C     ****************************************************************    I1150
  150 IPRNT=0                                                             I1155
      RETURN                                                              I1165
C     ****************************************************************    I1170
C                                                                         I1180
C                                                                         I1190
C                                                                         I1200
  160 FORMAT (1H1,13HCONCENTRATION/)                                      I1210
  170 FORMAT (1X,23HNUMBER OF TIME STEPS = ,1I5)                          I1220
  180 FORMAT (8X,16HDELTA T       = ,1G12.5)                              I1230
  190 FORMAT (8X,16HTIME(SECONDS) = ,1G12.5)                              I1240
  200 FORMAT (3X,21HCHEM.TIME(DAYS)    = ,1E12.5)                         I1250
  210 FORMAT (8X,16HTIME(YEARS)   = ,1E12.5)                              I1260
  220 FORMAT (1H )                                                        I1270
  230 FORMAT (1H1,23HCHANGE IN CONCENTRATION/)                            I1280
  240 FORMAT (1H0,20I5)                                                   I1290
  250 FORMAT (1H ,21HCHEMICAL MASS BALANCE)                               I1300
  260 FORMAT (8X,25HMASS IN BOUNDARIES     = ,1E12.5)                     I1310
  270 FORMAT (8X,25HMASS OUT BOUNDARIES    = ,1E12.5)                     I1320
  280 FORMAT (8X,25HMASS PUMPED OUT        = ,1E12.5)                     I1330
  290 FORMAT (8X,25HMASS PUMPED IN         = ,1E12.5)                     I1340
  295 FORMAT (8X,25HMASS LOST BY DECAY     = ,1E12.5)                     I1341R
  296 FORMAT (8X,25HMASS ADSORBED ON SOLIDS= ,1E12.5)                     I1342R
  298 FORMAT (8X,25HINITIAL MASS ADSORBED  = ,1E12.5)                     I1343R
  300 FORMAT (8X,25HINFLOW MINUS OUTFLOW   = ,1E12.5)                     I1350
  310 FORMAT (8X,25HINITIAL MASS DISSOLVED = ,1E12.5)                     I1361R
  320 FORMAT (8X,25HPRESENT MASS DISSOLVED = ,1E12.5)                     I1371R
  330 FORMAT (8X,25HCHANGE MASS DISSOLVED  = ,1E12.5)                     I1381R
  332 FORMAT (8X,25HCHANGE TOTL.MASS STORED= ,1E12.5)                     I1382R
  340 FORMAT (1H ,5X,53HCOMPARE RESIDUAL WITH NET FLUX AND MASS ACCUMULA  I1390
     1TION:)                                                              I1400
  350 FORMAT (8X,25HMASS BALANCE RESIDUAL  = ,1E12.5)                     I1410
  360 FORMAT (8X,25HERROR  (AS PERCENT)    = ,1E12.5)                     I1420
  370 FORMAT (1H ,5X,55HCOMPARE INITIAL MASS STORED WITH CHANGE IN MASS   I1430
     1STORED:)                                                            I1440
  380 FORMAT (1X,23H NO. MOVES COMPLETED = ,1I5)                          I1450
  390 FORMAT (1H1,10A8//)                                                 I1460
  400 FORMAT (1H0,5X,65HTIME VERSUS HEAD AND CONCENTRATION AT SELECTED O  I1470
     1BSERVATION POINTS//15X,19HPUMPING PERIOD NO. ,I4////)               I1480
  410 FORMAT (1H0,16X,19HTRANSIENT  SOLUTION////)                         I1490
  420 FORMAT (1H0,15X,21HSTEADY-STATE SOLUTION////)                       I1500
  430 FORMAT (1H0,20X,22HOBS.WELL NO.    X    Y,17X,1HN,6X,40HHEAD (FT)   I1510
     1   CONC.(MG/L)    TIME (YEARS)//24X,I3,9X,I2,3X,I2//)               I1520
  435 FORMAT (1H ,3X,45H** NOTE **  THIS OBS. WELL IS LOCATED OUTSIDE/16 SI1525
     1X,24HOF THE TRANSPORT SUBGRID)                                     SI1526
  440 FORMAT (1H ,58X,I2,6X,F7.1,8X,F7.1,8X,F7.3)                         I1531
  450 FORMAT (1H ,2X,21HCHEM.TIME(SECONDS) = ,E12.5)                      I1540
  460 FORMAT (1H ,2X,21HCHEM.TIME(YEARS)   = ,E12.5)                      I1550
      END                                                                 I1560-
