C**************************************************************
      PROGRAM POR2D
c
c-------------------------------------------------------------
c POR2D calculates the pore pressure distribution within an evolving
c sedimentary basin based on a vertical cross section.
c Sediment layers may have varying thickness along dip and varying
c lithology (sand or clay).
c Porosity is a function of effective vertical stress.
c Permeability is a log-linear function of porosity.
c Uppermost layer is assumed to be at constant pressure.
c Lower boundary is no-flow boundary.
c Lateral boundaries are constant gradient boundaries.
c
c Author: Johannes Wendebourg 12/92
c
c-------------------------------------------------------------
c Input file description:
c
c 1. line: number of columns in cross section, number of sediment layers,
c          grid spacing [m]
c 2. line: density of water [kg/m3], density of sediment [kg/m3],
c          viscosity of water [cP]
c 3. line: initial basement elevation for each column [m]
c
c Following are 4 lines for each sediment layer:
c
c 1. line: length of time interval covered by layer [MMYrs]
c 2. line: thickness for each column [m]
c          (layers must have thickness > 0 everywhere!)
c 3. line: water depth for each column [m]
c 4. line: lithology indicator for each column
c          [0 = sand, 1 = clay]
c-------------------------------------------------------------
c Output is the evolving grid through time for:
c
c Elevations, porosities, permeabilities, fluid pressures,
c lithostatic pressures, flow velocities and directions
c
c-------------------------------------------------------------
      implicit none
      integer*2 mxcol,mxlay,nmax
c-------------------------------------------------------------
c Change these parameters to increase grid size:
      parameter(mxcol=50,mxlay=20,nmax=1000)
c mxcol   maximum number of columns in grid
c mxlay   maximum number of layers in grid
c nmax    = mxcol*mxlay = maximum number of equations to solve
c-------------------------------------------------------------
      integer*2 i,j,k,m,n,nrow,ncol,nlay,
     &          indx(nmax),nwell
      real*8    bas(mxcol),top(mxcol),dtlay(mxlay),
     &          thk(mxlay,mxcol),elv(mxlay,mxcol),
     &          poro(mxlay,mxcol),sr(mxlay,mxcol),
     &          kfx(mxlay,mxcol),kfz(mxlay,mxcol),
     &          dz(mxlay,mxcol),p(mxlay,mxcol),wdp(mxlay,mxcol),
     &          pbdry(mxcol),plith(mxlay,mxcol),
     &          vel(mxlay,mxcol),ang(mxlay,mxcol),
     &          amat(nmax,nmax),rhs(nmax),a(5),
     &          dx,dx2,rhow,rhos,visc,gams,gamw,subs,
     &          time,dt,elvat,d,p1,p2,dh,vx,vz,vscal,fc,
     &          dpotx,dpotz,xdia,vcomp,solid,dphi,vpore,
     &          porfunc,perm,dpor,pi,grav,secyear
      DATA pi,grav,secyear /3.1415927,9.81,31560000.0/
c
c-------------------------------------------------------------
c read data
c
c open input file
      OPEN(20, FILE='porflow.input')
      REWIND (20)
c read number of columns, layers, grid spacing
      READ (20, *) NCOL, NLAY, DX
c read densities, viscosity
      READ (20, *) RHOW, RHOS, VISC
c read initial basement
      READ (20, *) (BAS(J), J = 1, NCOL)
c read properties of layers
      DO I = 1, NLAY
         READ (20, *) DTLAY(I)
         READ (20, *) (THK(I,J), J = 1, NCOL)
         READ (20, *) (WDP(I,J), J = 1, NCOL)
         READ (20, *) (SR(I,J), J = 1, NCOL)
      END DO
c read vertical exaggeration, velocity scaling factor,well column
      READ (20, *) FC, VSCAL, NWELL
c close input file
      CLOSE(20)
c convert viscosity from cP to Pa.s
      VISC = VISC*10.0
c gravities of sediment and water
      GAMS = GRAV*RHOS
      GAMW = GRAV*RHOW
c dx squared
      DX2 = DX*DX
c
c-------------------------------------------------------------
c begin of time loop
c
      TIME = 0.0
      DO K = 1, NLAY
         NROW = K
c length of time step is duration of sedimentation
c new layer
         DT = DTLAY(K)*SECYEAR*1.0E6
c time is age of top of layer stack
         TIME = TIME + DTLAY(K)
c porosity and thickness of uppermost (loading) layer
         DO J = 1, NCOL
            PORO(K,J) = PORFUNC(SR(K,J),0.0)
            THK(K,J) = THK(K,J)/(1.0+PORO(K,J))
         END DO
c
c-------------------------------------------------------------
c calculate lithostatic pressure for each grid cell
c
         DO J = 1, NCOL
            PLITH(K,J) = 0.0
            IF (K .GT. 1) THEN
               DO I = K - 1, 1, -1
                  PLITH(I,J) = PLITH(I+1,J) + GAMS*THK(I+1,J)
               END DO
            ENDIF
         END DO
c
c-------------------------------------------------------------
c determine porosities, thicknesses and permeabilities
c
c calculate porosity from effective pressure:
c effective pressure = lithostatic - excess pore pressure
c note that p includes hydrostatic as well as excess fluid pressure
         DO J = 1, NCOL
            DO I = 1, NROW
               IF (K.GT.1 .AND. I.LT.NROW) PORO(I,J) = PORFUNC(SR(I,J),
     #            PLITH(I,J)-P(I,J)+PBDRY(J))
               DZ(I,J) = THK(I,J)*(1.0+PORO(I,J))
               KFX(I,J) = PERM(1,SR(I,J),PORO(I,J))
               KFZ(I,J) = PERM(2,SR(I,J),PORO(I,J))
            END DO
         END DO
c
c-------------------------------------------------------------
c determine new elevations and boundary pressures
c
         DO J = 1, NCOL
c subside basement by thickness of loading layer
c and by change in water depth which may be a combination
c of eustatic sea level change and tectonic subsidence
            SUBS = DZ(K,J)
            IF (WDP(K,J) .GT. 0.0) THEN
               IF (K .EQ. 1) THEN
                  SUBS = SUBS + WDP(K,J)
               ELSE
                  SUBS = SUBS + (WDP(K,J)-WDP(K-1,J))
               ENDIF
            ENDIF
            BAS(J) = BAS(J) - SUBS
            ELVAT = BAS(J)
            DO I = 1, NROW
c elevation is measured at mid-point of layer
               ELV(I,J) = ELVAT + DZ(I,J)*0.5
               ELVAT = ELVAT + DZ(I,J)
            END DO
            TOP(J) = ELVAT
         END DO
c determine boundary pressure
         DO J = 1, NCOL
            PBDRY(J) = 0.0
c Simulate topographic driven flow by setting hydrostatic pressure
c at topography above sea level assuming that topography = water table
            IF (TOP(J) .GT. 0.0) PBDRY(J) = GAMW*TOP(J)
c set upper layer pressure = hydrostatic
            P(NROW,J) = PBDRY(J)
         END DO
c
c-------------------------------------------------------------
c initialize coefficient matrix and right hand side vector
c n is number of equations
         N = NCOL*NROW
         DO I = 1, N
            RHS(I) = 0.0
            DO J = 1, N
               AMAT(I,J) = 0.0
            END DO
         END DO
c
c-------------------------------------------------------------
c determine banded n*n matrix (amat) and right hand side vector (rhs)
c this matrix is pentadiagonal
c the coefficients of a row are ordered in the following manner:
c .... a(1) (nrow-2)*0 a(2) a(5) a(3) (nrow-2)*0 a(4) ....
         N = 0
         DO J = 1, NCOL
            DO I = 1, NROW
               N = N + 1
c current pore volume of grid cell
               VPORE = DZ(I,J)*DX2*PORO(I,J)
c solid fraction of grid cell
               SOLID = 1.0 - PORO(I,J)
c effective pressure derivative
               DPHI = DPOR(SR(I,J),PORO(I,J))
c initialize temporary coefficient arrays
               DO M = 1, 5
                  A(M) = 0.0
               END DO
c determine off-diagonal coefficients
               IF (J .GT. 1) THEN
                  D = 0.5*(DZ(I,J)+DZ(I,J-1))
                  A(1)=2.0/(VISC*(1.0/(KFX(I,J)*D)+1.0/(KFX(I,J-1)*D)))
               ELSE
                  A(1) = KFX(I,J)*DZ(I,J)/VISC
               ENDIF
               IF (I .GT. 1) THEN
                  D = DX2/(VISC*(DZ(I,J)/KFZ(I,J)+DZ(I-1,J)/KFZ(I-1,J)))
                  A(2) = 2.0*D
               ENDIF
               IF (I .LT. NROW) THEN
                  D = DX2/(VISC*(DZ(I,J)/KFZ(I,J)+DZ(I+1,J)/KFZ(I+1,J)))
                  A(3) = 2.0*D
               ELSE
                  A(3) = (KFZ(I,J)*DX2)/(DZ(I,J)*VISC)
               ENDIF
               IF (J .LT. NCOL) THEN
                  D = 0.5*(DZ(I,J)+DZ(I,J+1))
                  A(4)=2.0/(VISC*(1.0/(KFX(I,J)*D)+1.0/(KFX(I,J+1)*D)))
               ELSE
                  A(4) = KFX(I,J)*DZ(I,J)/VISC
               ENDIF
c main diagonal coefficient is negative sum of off-diagonal coefficients
               DO M = 1, 4
                  A(5) = A(5) - A(M)
               END DO
c add compaction term to main diagonal
               VCOMP = (VPORE/SOLID)*DPHI/DT
               A(5) = A(5) + VCOMP
c determine right hand side vector:
               RHS(N) = 0.0
c compaction term
               RHS(N) = RHS(N) + VCOMP*(P(I,J)+THK(K,J)*GAMS)
c flux at upper boundary
               IF (I .EQ. NROW) RHS(N) = RHS(N) - A(3)*PBDRY(J)
c flux term at lateral boundaries
               IF(J.EQ.1)RHS(N)=RHS(N)-A(1)*(2.0*P(I,J+1)-P(I,J+2))
               IF(J.EQ.NCOL)RHS(N)=RHS(N)-A(4)*(2.0*P(I,J-1)-P(I,J-2))
c fill in matrix and normalize coefficients such that main diagonal = 1
               AMAT(N,N) = 1.0
               IF (J .GT. 1) AMAT(N,N-NROW) = A(1)/A(5)
               IF (I .GT. 1) AMAT(N,N-1) = A(2)/A(5)
               IF (I .LT. NROW) AMAT(N,N+1) = A(3)/A(5)
               IF (J .LT. NCOL) AMAT(N,N+NROW) = A(4)/A(5)
c normalize also right hand side vector
               RHS(N) = RHS(N)/A(5)
            END DO
         END DO
c
c-------------------------------------------------------------
c solve system of equations by LU decomposition
c
c solution vector replaces right hand side vector
         CALL LUDCMP (AMAT, N, 1000, INDX, D)
         CALL LUBKSB (AMAT, N, 1000, INDX, RHS)
c save fluid pressure
         N = 0
         DO J = 1, NCOL
            DO I = 1, NROW
               N = N + 1
               P(I,J) = RHS(N)
            END DO
         END DO
c
c-------------------------------------------------------------
c determine pore water velocities and directions
c
         DO J = 1, NCOL
            DO I = 1, NROW
c calculate lateral pressure gradient and lateral slope
               IF (J.GT.1 .AND. J.LT.NCOL) THEN
                  DPOTX = 0.5*(P(I,J+1)-P(I,J-1))
                  DH = 0.5*(ELV(I,J+1)-ELV(I,J-1))
               ELSE IF (J .EQ. 1) THEN
                  DPOTX = P(I,2) - P(I,1)
                  DH = ELV(I,2) - ELV(I,1)
               ELSE IF (J .EQ. NCOL) THEN
                  DPOTX = P(I,NCOL) - P(I,NCOL-1)
                  DH = ELV(I,NCOL) - ELV(I,NCOL-1)
               ENDIF
c calculate vertical pressure gradient
               IF (I.GT.1 .AND. I.LT.NROW) THEN
                  D = 0.5*DZ(I+1,J)/(0.5*(DZ(I+1,J)+DZ(I,J)))
                  P1 = D*P(I,J) + (1.0-D)*P(I+1,J)
                  D = 0.5*DZ(I,J)/(0.5*(DZ(I-1,J)+DZ(I,J)))
                  P2 = D*P(I-1,J) + (1.0-D)*P(I,J)
                  DPOTZ = P1 - P2
               ELSE IF (I.EQ.1 .AND. K.GT.1) THEN
                  D = 0.5*DZ(2,J)/(0.5*(DZ(2,J)+DZ(1,J)))
                  P1 = D*P(1,J) + (1.0-D)*P(2,J)
                  DPOTZ = 2.0*(P1-P(1,J))
               ELSE IF (I .EQ. NROW) THEN
                  DPOTZ = PBDRY(J) - P(NROW,J)
               ENDIF
c calculate magnitude of Darcy velocities
               VX = -DPOTX*KFX(I,J)/(VISC*DX)
               VZ = -DPOTZ*KFZ(I,J)/(VISC*DZ(I,J))
c true fluid velocity is Darcy velocity / porosity
               VEL(I,J) = DSQRT(VX*VX+VZ*VZ)/PORO(I,J)
c calculate direction of velocity
               ANG(I,J) = 0.0
               IF (VEL(I,J) .GT. 0.0) THEN
                  XDIA = VX/VEL(I,J)
                  XDIA = AMIN1(1.0,XDIA)
                  XDIA = AMAX1(-1.0,XDIA)
                  IF (VZ .GE. 0.0) THEN
                     ANG(I,J) = DACOS(XDIA)*180./PI
                  ELSE
                     ANG(I,J) = 270. + DASIN(XDIA)*180./PI
                  ENDIF
c correct for local slope and vertical exaggeration
                  ANG(I,J) = ANG(I,J) + DATAN(DH*FC/DX)*180.0/PI
                  IF (ANG(I,J) .LT. 0.0) ANG(I,J) = 360.0 + ANG(I,J)
                  IF (ANG(I,J) .GT. 360.0) ANG(I,J) = ANG(I,J) - 360.0
               ENDIF
            END DO
         END DO
c
c-------------------------------------------------------------
c write out resulting grids:
c
         WRITE (6, '(///A,F10.3)') ' Time [MMYrs] = ', TIME
c elevation
         WRITE (6, '(/A/)') ' Elevation [m] '
         WRITE (6, '(10F8.1)') (TOP(J), J = 1, NCOL)
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.1)') (ELV(I,J), J = 1, NCOL)
         END DO
         WRITE (6, '(10F8.1)') (BAS(J), J = 1, NCOL)
c lithology
         WRITE (6, '(/A/)') ' Lithology [0=sand, 1=shale] '
         DO I = NROW, 1, -1
            WRITE (6, '(I4,9I8)') (INT(SR(I,J)), J = 1, NCOL)
         END DO
c thickness
         WRITE (6, '(/A/)') ' Thickness [m] '
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.1)') (DZ(I,J), J = 1, NCOL)
         END DO
c porosity
         WRITE (6, '(/A/)') ' Porosity '
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.4)') (PORO(I,J), J = 1, NCOL)
         END DO
c vertical permeability
         WRITE (6, '(/A/)') ' log Permeability z [md]'
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.3)') (DLOG10(KFZ(I,J)*9.87E16), J = 1, NCOL
     #         )
         END DO
c excess fluid pressure
         WRITE (6, '(/A/)') ' Excess fluid pressure [MPa] '
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.3)') ((P(I,J)-PBDRY(J))*1.0E-6, J = 1, NCOL
     #         )
         END DO
c boundary pressure
         WRITE (6, '(/A/)') ' Boundary pressure [MPa] '
         WRITE (6, '(10F8.3)') (PBDRY(J)*1.0E-6, J = 1, NCOL)
c sediment load
         WRITE (6, '(/A/)') ' Lithostatic pressure [MPa] '
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.3)') (PLITH(I,J)*1.0E-6, J = 1, NCOL)
         END DO
c flow velocities
         WRITE (6, '(/A/)') ' Flow velocities [mm/a] '
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.3)') (VEL(I,J)*SECYEAR*1000.0, J = 1, NCOL)
         END DO
c flow directions
         WRITE (6, '(/A/)') ' Flow directions '
         DO I = NROW, 1, -1
            WRITE (6, '(10F8.0)') (ANG(I,J), J = 1, NCOL)
         END DO
c
c-------------------------------------------------------------
c end of time loop
      END DO
c
      END
c
c***********************************************************************
c
c Porosity is a function of effective vertical stress sigma
c Porosity                phi = phi(zero) * exp( -sigma / sigma(zero) )
c Derivative of porosity  d phi/d sigma = - phi / sigma(zero)
c
      function porfunc(sr,p)
      implicit none
      real*8     sr,p,porfunc,phi0sd,phi0sh,sig0sd,sig0sh
      DATA phi0sd,phi0sh,sig0sd,sig0sh /0.43,0.58,27.0,18.0/
c
      IF (SR .EQ. 0.) PORFUNC = PHI0SD*DEXP((-P*1.0E-6/SIG0SD))
      IF (SR .EQ. 1.) PORFUNC = PHI0SH*DEXP((-P*1.0E-6/SIG0SH))
c
      RETURN 
      END
c
      function dpor(sr,phi)
      implicit none
C...Translated by FPP 5.0 (3.03Q7) 03/04/93  15:57:19  
C...Switches: renumb=1:1,format=100:1,indal=3,conchr=#
      real*8     sr,phi,dpor,sig0sd,sig0sh
      DATA sig0sd,sig0sh /27.0,18.0/
c
      IF (SR .EQ. 0.) DPOR = -(1.0/(SIG0SD*1.0E6))*PHI
      IF (SR .EQ. 1.) DPOR = -(1.0/(SIG0SH*1.0E6))*PHI
      RETURN 
      END
c
c***********************************************************************
c
c Permeability is a function of porosity
c a log-linear function  log(k) = a + b * phi
c f is an anisotropy factor for vertical permeability kz = 1/f * kx
c
      function perm(i,sr,poro)
      implicit none
      integer*2 i
      real*8    sr,poro,perm,fsd,fsh,asd,ash,bsd,bsh
      DATA fsd,fsh,asd,ash,bsd,bsh /2.5,10.0,-13.0,-19.0,2.0,8.0/
c
      IF (SR .EQ. 0.) THEN
         PERM = 10.0**(ASD+BSD*PORO)
         IF (I .EQ. 2) PERM = (1.0/FSD)*PERM
      ELSE IF (SR .EQ. 1.) THEN
         PERM = 10.0**(ASH+BSH*PORO)
         IF (I .EQ. 2) PERM = (1.0/FSH)*PERM
      ENDIF
c
      RETURN 
      END
************************************************************************
*       Given an NxN matrix A, with physical dimenion NP, this routine *
*       replaces it by the LU decomposition of a rowwise permutation of*
*       itself(Numerical Recipies, 1986; Copyright)                    *
************************************************************************
************************************************************************
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      Implicit None
      Integer*2 NMAX
      Real*8 TINY
      PARAMETER (NMAX=1000,TINY=1.0E-20)
      Real*8 A(NP,NP),VV(NMAX),D,aamax,sum,dum
      Integer*2 indx(n),np,n,j,i,k,imax
      D = 1.
      DO 2 I = 1, N
         AAMAX = 0.
         DO 1 J = 1, N
            AAMAX = AMAX1(ABS(A(I,J)),AAMAX)
    1    CONTINUE
         IF (AAMAX .EQ. 0.) THEN
            PAUSE 'Singular matrix.'
         ENDIF
         VV(I) = 1./AAMAX
    2 CONTINUE
      DO 9 J = 1, N
         IF (J .GT. 1) THEN
            DO 4 I = 1, J - 1
               SUM = A(I,J)
               IF (I .GT. 1) THEN
                  DO 3 K = 1, I - 1
                     SUM = SUM - A(I,K)*A(K,J)
    3             CONTINUE
                  A(I,J) = SUM
               ENDIF
    4       CONTINUE
         ENDIF
         AAMAX = 0.
         DO 6 I = J, N
            SUM = A(I,J)
            IF (J .GT. 1) THEN
               DO 5 K = 1, J - 1
                  SUM = SUM - A(I,K)*A(K,J)
    5          CONTINUE
               A(I,J) = SUM
            ENDIF
            DUM = VV(I)*ABS(SUM)
            IF (DUM .GE. AAMAX) THEN
               IMAX = I
               AAMAX = DUM
            ENDIF
    6    CONTINUE
         IF (J .NE. IMAX) THEN
            DO 7 K = 1, N
               DUM = A(IMAX,K)
               A(IMAX,K) = A(J,K)
               A(J,K) = DUM
    7       CONTINUE
            D = -D
            VV(IMAX) = VV(J)
         ENDIF
         INDX(J) = IMAX
         IF (J .NE. N) THEN
            IF (A(J,J) .EQ. 0.) A(J,J) = 1.0E-20
            DUM = 1./A(J,J)
            DO 8 I = J + 1, N
               A(I,J) = A(I,J)*DUM
    8       CONTINUE
         ENDIF
    9 CONTINUE
      IF (A(N,N) .EQ. 0.) A(N,N) = 1.0E-20
      RETURN 
      END
************************************************************************
*       Solves the set of n linear equations AX=B where A is from      *
*       Subroutine LUDCMP                                              *
*       (Numerical Recipies, 1986; Copyright)                          *
************************************************************************
************************************************************************
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
      Implicit None
      real*8 A(NP,NP),B(N),SUM
      Integer*2 INDX(N),II,I,N,LL,J,NP
      II = 0
      DO 2 I = 1, N
         LL = INDX(I)
         SUM = B(LL)
         B(LL) = B(I)
         IF (II .NE. 0) THEN
            DO 1 J = II, I - 1
               SUM = SUM - A(I,J)*B(J)
    1       CONTINUE
         ELSE IF (SUM .NE. 0.) THEN
            II = I
         ENDIF
         B(I) = SUM
    2 CONTINUE
      DO 4 I = N, 1, -1
         SUM = B(I)
         IF (I .LT. N) THEN
            DO 3 J = I + 1, N
               SUM = SUM - A(I,J)*B(J)
    3       CONTINUE
         ENDIF
         B(I) = SUM/A(I,I)
    4 CONTINUE
      RETURN 
      END
