A.2 Codes of Groundwater Flow

!****************************************!  HEADS.f90

!  By Samuel Sangwon Lee, 1999

!****************************************  

   IMPLICIT NONE

   INTEGER*4 MPX,MPY,MPX1,MPY1 

! Number of pixels in x and y directions

   REAL*8 DX,DY  ! Pixel size

   PARAMETER (MPX=90,MPY=56,MPX1=MPX+1,MPY1=MPY+1,DX=10.0,DY=1.0)

   REAL*8 H(-1:MPX1,0:MPY1),X(-1:MPX1),Y(0:MPY1),DEPTH1(0:MPX),DEPTH2(0:MPX) &,HBEDI,ERMAX,OLDH,ERR

   INTEGER*4 NDEP(0:MPX),NDEPI,I,J,NUMIT

   CHARACTER DUMMY*1

! INITIALIZE ALL H(I,J) VALUES TO BE 0

   DO J=0,MPY1                DO I=-1,MPX1

   H(I,J) = 0.0    END DO       END DO

! WATER TABLE BOUNDARY

   OPEN(1,FILE='S600.tbl',STATUS='OLD')

   READ(1,'(A)')DUMMY 

! Reading variable labels (x, h, DEPTH)

   DO I=0,MPX

READ(1,*)X(I),H(I,0),DEPTH1(I),DEPTH2(I)

   END DO               CLOSE(1)

   X(-1) = X(0) - DX 

! Defining x and y values

   X(MPX1) = X(MPX) + DX

   DO J=0,MPY1

Y(J) = -DY*J     END DO

   DO I=0,MPX

  NDEP(I) = -DEPTH2(I)/DY

   END DO

! KEEP TRACK OF NUMBER OF ITERATIONS AND !OF LARGEST ERROR NO-FLOW BOUNDARIES NEED !TO BE RESET WITHIN EACH ITERATION LOOP

   NUMIT = 0

10 ERMAX = 0.0

   NUMIT = NUMIT + 1

! KEEPING NO-FLOW BOUNDARY CONDITION

   DO J=0,MPY  ! Left & right boundaries

   H(-1,J) = H(1,J)

   H(MPX1,J) = H(MPX-1,J)     END DO

   DO I=0,MPX 

! Bottom boundary at bedrock

     NDEPI = NDEP(I)

     H(I,NDEPI+1) = H(I,NDEPI-1)

     HBEDI = H(I,NDEPI+1)

     DO J=NDEPI+2,MPY1

       H(I,J) = HBEDI

     END DO             END DO

! SWEEP INTERIOR POINTS WITH 5-POINT !OPERATOR

   DO J=1,MPY           DO I=0,MPX

     IF (J.LE.NDEP(I)) THEN 

! Stop at the bedrock surface

    OLDH = H(I,J)

    H(I,J) = (H(I-1,J) + H(I+1,J) + H(I,J-1) + H(I,J+1))/4.0

    ERR = ABS(H(I,J) - OLDH)

    IF (ERR.GT.ERMAX) ERMAX=ERR

  END IF       END DO        END DO

! DO ANOTHER ITERATION IF LARGEST ERROR !AFFECTS 7TH DECIMAL PLACE

   IF(ERMAX.GT.0.000001) GO TO 10

! WE ARE DONE.  WRITE ONTO THE OUTPUT !FILE

   OPEN (1,FILE='S600.hed')

   WRITE (1,100)NUMIT

100 FORMAT(' NUMBER OF ITERATIONS: ',I5/'  X(ft)'4X'Y(ft)'6X'H(ft)')

   DO I=-1,MPX1   DO J=0,MPY1

     WRITE(1,110)X(I),Y(J),H(I,J)

   END DO               END DO

110 FORMAT(1X,F7.1,2X,F7.1,2X,F9.6)

   CLOSE(1)      STOP             END

!****************************************!       FLOWS.f90

!**************************************** 

   IMPLICIT NONE

   INTEGER*4 MPX,MPY,MPX1,MPY1 

! Number of pixels in x and y directions

   REAL*8 DX,DY  ! Pixel size

   PARAMETER (MPX=90,MPY=56,MPX1=MPX+1,MPY1=MPY+1,DX=10.0,DY=1.0)

   REAL*8 H(-1:MPX1,0:MPY1),X(-1:MPX1),Y(0:MPY1),DEPTH1(0:MPX),DEPTH2(0:MPX) & ,FX(-1:MPX,0:MPY),FY(-1:MPX,0:MPY),K(-1:MPX,0:MPY) &,K1,K2,HX1, HX2,HY1,HY2

   INTEGER*4 NDEP1,NDEP2,I,J

 CHARACTER FTBL*40,FIN*40,FOUT*40,DUMMY*1

! READ DATA FILE NAMES

   FTBL = 'S600.tbl'

   WRITE(*,'(A$)')' Input head value (H) file: '

   READ(*,'(A)')FIN

   FIN = 'S600.hed'

   WRITE(*,'(A$)')' Output flow rate (F) file: '

   READ(*,'(A)')FOUT

   FOUT = 'S600.flw'

! READ HYDRAULIC CONDUCTIVITY

   WRITE(*,'(A$)')' Hydraulic conductivity (1st layer) : '

   READ(*,*)K1

   WRITE(*,'(A$)')' Hydraulic conductivity (2nd layer) : '

   READ(*,*)K2

! READ WATER TABLE & LAYER DEPTH DATA

   OPEN(1,FILE=FTBL,STATUS='OLD')

   READ(1,'(A)')DUMMY 

! Reading variable labels (x, h, DEPTH1, !DEPTH2)

   DO I=0,MPX

 READ(1,*)X(I),H(I,0),DEPTH1(I),DEPTH2(I)

   END DO               CLOSE(1)

! READ H(I,J) VALUES FROM EXISTING FILE

   OPEN(1,FILE=FIN)

   READ(1,100)DUMMY  ! Skip two lines

   READ(1,100)DUMMY

100 FORMAT(A1)

   DO I=-1,MPX1 DO J=0,MPY1

     READ(1,*)X(I),Y(J),H(I,J)

   END DO       END DO    CLOSE(1)

! DEFINE THE FLOW RATE PIXELS SHIFTED BY !HALF A PIXEL

   DO I=-1,MPX

     X(I) = (X(I) + X(I+1))/2

   END DO

   DO J=0,MPY

     Y(J) = (Y(J) + Y(J+1))/2

   END DO

! SET HYDRAULIC CONDUCTIVITY K(I,J) !VALUES

   DEPTH1(-1) = DEPTH1(0)

   DEPTH2(-1) = DEPTH2(0)

   DO I=-1,MPX

     NDEP1 = -DEPTH1(I)/DY

     NDEP2 = -DEPTH2(I)/DY

     DO J=0,NDEP1-1

       K(I,J) = K1  ! First layer

     END DO

     DO J=NDEP1,NDEP2-1

       K(I,J) = K2  ! Second layer

     END DO

     DO J=NDEP2,MPY

       K(I,J) = 0   ! Bedrock

     END DO             END DO

! CALCULATE THE FLOW RATES

   DO J=0,MPY           DO I=-1,MPX

   HX1 = (H(I,J) + H(I,J+1)) / 2.0

   HX2 = (H(I+1,J) + H(I+1,J+1)) / 2.0

   FX(I,J) = -K(I,J) * (HX2 - HX1) / DX

   HY1 = (H(I,J) + H(I+1,J)) / 2.0

   HY2 = (H(I,J+1) + H(I+1,J+1)) / 2.0

   FY(I,J) = -K(I,J) * (HY1 - HY2) / DY

   END DO               END DO

! WRITE ONTO A FILE

   OPEN (1,FILE=FOUT)

   WRITE (1,110)K1, K2

110 FORMAT(' Hydraulic conductivity: 'F7.1' (1st Layer), 'F7.1' (2nd Layer)')

   WRITE(1,115)

115FORMAT('X(ft)'4X'Y(ft)'9X'Fx'11X' Fy')

   DO J=0,MPY           DO I=-1,MPX

   WRITE(1,120)X(I),Y(J),FX(I,J),FY(I,J)

   END DO               END DO

120 FORMAT(1X,F7.1,2X,F7.1,1X,F12.7, 1X,F12.6)

   CLOSE(1)      STOP             END

!***************************************

!       VFLOW.f90

!***************************************

   IMPLICIT NONE

   INTEGER*4 MPX,MPY,MPX1,MPY1 

! Number of pixels in x and y directions

   REAL*8 DX,DY  ! Pixel size

   PARAMETER (MPX=90,MPY=56,MPX1=MPX+1,MPY1=MPY+1,DX=10.0,DY=1.0)

   REAL*8 X(-1:MPX),Y(0:MPY), V(0:MPX,0:MPY),FX(-1:MPX,0:MPY),FY(-1:MPX,0:MPY) &  ,VIJ1,VIJ2,VIJ3,VIJ4, VOLD,ERR,ERMAX

   INTEGER*4 I,J,NUMIT

   CHARACTER FIN*40,FOUT*40,DUMMY*1

! READ FILE NAMES

   WRITE(*,'(A$)')' Input flow rate (F) file : '

   READ(*,'(A)')FIN

   FIN = 'S600.flw'

   WRITE(*,'(A$)')' Output flow function (V) file : '

   READ(*,'(A)')FOUT

   FOUT = 'S600.vxy'

! READ FX(I,J) & FY(I,J) VALUES FROM !EXISTING FILE

   OPEN(1,FILE=FIN)

   READ(1,100)DUMMY  ! Skip labels

   READ(1,100)DUMMY

100 FORMAT(A1)

   DO J=0,MPY           DO I=-1,MPX

   READ(1,*)X(I),Y(J),FX(I,J),FY(I,J)

   END DO      END DO        CLOSE(1)

! DEFINE THE FLOW FUNCTION PIXELS SHIFTED !BY HALF A PIXEL

   DO I=MPX,0,-1

     X(I) = ( X(I) + X(I-1) ) / 2

   END DO

      DO J=MPY,1,-1

     Y(J) = ( Y(J) + Y(J-1) ) / 2

   END DO

   Y(0) = 0

! REDEFINE FX & FY

   DO I=-1,MPX            DO J=MPY,1,-1

   FY(I,J) = ( FY(I,J) + FY(I,J-1) ) / 2

   END DO               END DO

   DO J=0,MPY           DO I=0,MPX

     FX(I,J)= ( FX(I,J) + FX(I-1,J) ) / 2

   END DO               END DO

! INITIALIZE THE FLOW FUNCTION V(X,Y)

   DO I=0,MPX           DO J=0,MPY

     V(I,J) = 0

   END DO              END DO

! CALCULATE THE FLOW FUNCTION V(X,Y)

! V(X,Y)=0 AT THE RIGHT, LEFT, BOTTOM !BOUNDARIES

  NUMIT = 0  ! Number of iterations set 0

10  ERMAX = 0 

! Initialize error maximum to 0

  NUMIT = NUMIT + 1

  DO I=1,MPX-1

    VOLD = V(I,0)

    V(I,0) = V(I,1) + FX(I,0)*DY

    ERR = ABS( V(I,0) - VOLD )

    IF (ERR.GT.ERMAX) ERMAX = ERR

  DO J=1,MPY-1

    VOLD = V(I,J)

    VIJ1 = V(I,J+1) + FX(I,J)*DY

    VIJ2 = V(I,J-1) + FX(I,J-1)*(-DY)

    VIJ3 = V(I-1,J) - FY(I-1,J)*DX

    VIJ4 = V(I+1,J) - FY(I,J)*(-DX)

    V(I,J)=(VIJ1 + VIJ2 + VIJ3 + VIJ4)/4

       ERR = ABS( V(I,J) - VOLD )

       IF (ERR.GT.ERMAX) ERMAX = ERR

     END DO             END DO

   IF (ERMAX .GT. 0.001) GOTO 10 

! Repeat until error < 0.001

! WRITE ONTO A FILE

   OPEN (1,FILE=FOUT)

   WRITE(1,105)'Number of iterations : ',NUMIT

105 FORMAT(1X,A,I5)

   WRITE(1,110)'X(ft)','Y(ft)','V(X,Y)'

110 FORMAT(1X,2X,A,6X,A,7X,A)

   DO J=0,MPY           DO I=0,MPX

     WRITE(1,120)X(I),Y(J),V(I,J)

   END DO               END DO

120 FORMAT(1X,F7.1,2X,F7.1,2X,F12.3)

   CLOSE(1)       STOP      END