      MODULE DE4MODULE
        INTEGER,SAVE,POINTER  ::MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR
        INTEGER,SAVE,POINTER  ::NITERDE4,IFREQ,IPRD4,MUTD4,ID4DIM
        INTEGER,SAVE,POINTER  ::NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
        REAL   ,SAVE,POINTER  ::ACCLDE4,HCLOSEDE4,DELTL
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   ::IUPPNT
        INTEGER,          SAVE, POINTER, DIMENSION(:,:,:) ::IEQPNT
        REAL,             SAVE, POINTER, DIMENSION(:,:)   ::AU
        REAL,             SAVE, POINTER, DIMENSION(:,:)   ::AL
        REAL,             SAVE, POINTER, DIMENSION(:)     ::D4B
        REAL,             SAVE, POINTER, DIMENSION(:)     ::HDCGDE4
        INTEGER,          SAVE, POINTER, DIMENSION(:,:)   ::LRCHDE4
      TYPE DE4TYPE
        INTEGER,POINTER  ::MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR
        INTEGER,POINTER  ::NITERDE4,IFREQ,IPRD4,MUTD4,ID4DIM
        INTEGER,POINTER  ::NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW
        REAL   ,POINTER  ::ACCLDE4,HCLOSEDE4,DELTL
        INTEGER,           POINTER, DIMENSION(:,:)   ::IUPPNT
        INTEGER,           POINTER, DIMENSION(:,:,:) ::IEQPNT
        REAL,              POINTER, DIMENSION(:,:)   ::AU
        REAL,              POINTER, DIMENSION(:,:)   ::AL
        REAL,              POINTER, DIMENSION(:)     ::D4B
        REAL,              POINTER, DIMENSION(:)     ::HDCGDE4
        INTEGER,           POINTER, DIMENSION(:,:)   ::LRCHDE4
      END TYPE
      TYPE(DE4TYPE), SAVE ::DE4DAT(10)
      END MODULE DE4MODULE

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

      SUBROUTINE DE47AR(INDE4,MXITER,IGRID)
C     ******************************************************************
C     ALLOCATE STORAGE IN X ARRAY FOR D4 ARRAYS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,   ONLY:IOUT,NCOL,NROW,NLAY
      USE DE4MODULE
      use utl7module, only: URDCOM, URWORD
      use SimPHMFModule, only: USTOP
      CHARACTER*200 LINE
      double precision :: r
C     ------------------------------------------------------------------
      ALLOCATE(MXUP,MXLOW,MXEQ,MXBW,ITMX,ID4DIR)
      ALLOCATE(NITERDE4,IFREQ,IPRD4,MUTD4,ID4DIM)
      ALLOCATE(NBWL,NUPL,NLOWL,NLOW,NEQ,NUP,NBW)
      ALLOCATE(ACCLDE4,HCLOSEDE4,DELTL)
C
C1------PRINT A MESSAGE IDENTIFYING DE4 PACKAGE.
      WRITE(IOUT,1) INDE4
    1 FORMAT(1X,/1X,'DE4 -- D4 DIRECT SOLUTION PACKAGE',
     1    ', VERSION 7, 5/2/2005  INPUT READ FROM UNIT ',I4)
C
C2------CALCULATE DEFAULT VALUES FOR MXUP, MXLOW, AND MXBW.
C2------ALSO SET VALUES FOR ID4DIR (DIRECTION OF EQUATION ORDERING) AND
C2------ID4DIM (MAXIMUM NUMBER OF HEAD COEFFICIENTS FOR ONE EQUATION).
C2------ID4DIM=5 for a 2-D problem; ID4DIM=7 for 3-D.
      NODES=NCOL*NROW*NLAY
      NHALFU=(NODES-1)/2 +1
      NHALFL=NODES-NHALFU
      ID4DIM=7
      DELTL = 0.
      NBWL = 0
      NUPL = 0
      NLOWL = 0
      IF(NLAY.LE.NCOL .AND. NLAY.LE.NROW) THEN
         IF(NLAY.EQ.1) ID4DIM=5
         IF(NCOL.GE.NROW) THEN
            ID4DIR=1
            NBWGRD=NROW*NLAY+1
         ELSE
            ID4DIR=2
            NBWGRD=NCOL*NLAY+1
         END IF
      ELSE IF(NROW.LE.NCOL .AND. NROW.LE.NLAY) THEN
         IF(NROW.EQ.1) ID4DIM=5
         IF(NCOL.GE.NLAY) THEN
            ID4DIR=3
            NBWGRD=NROW*NLAY+1
         ELSE
            ID4DIR=4
            NBWGRD=NROW*NCOL+1
         END IF
      ELSE
         IF(NCOL.EQ.1) ID4DIM=5
         IF(NROW.GE.NLAY) THEN
            ID4DIR=5
            NBWGRD=NCOL*NLAY+1
         ELSE
            ID4DIR=6
            NBWGRD=NCOL*NROW+1
         END IF
      END IF
C
C3------READ AND PRINT COMMENTS, ITMX, MXUP, MXLOW, MXBW.  FOR ANY
C3------ZERO OR NEGATIVE VALUES, SUBSTITUE THE DEFAULT VALUE.
      CALL URDCOM(INDE4,IOUT,LINE)
      LLOC=1
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ITMX,R,IOUT,INDE4)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXUP,R,IOUT,INDE4)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXLOW,R,IOUT,INDE4)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXBW,R,IOUT,INDE4)
      IF(ITMX.LT.1) ITMX=1
      WRITE(IOUT,3) ITMX
    3 FORMAT(1X,'MAXIMUM ITERATIONS (EXTERNAL OR INTERNAL) =',I3)
      IF(MXUP.LT.1) MXUP=NHALFU
      IF(MXLOW.LT.1) MXLOW=NHALFL
      MXEQ=MXUP+MXLOW
      IF(MXBW.LT.1) MXBW=NBWGRD
      WRITE(IOUT,4) MXUP,MXLOW,MXBW
    4 FORMAT(1X,'MAXIMUM EQUATIONS IN UPPER PART OF [A]:',I7,/
     1       1X,'MAXIMUM EQUATIONS IN LOWER PART OF [A]:',I7,/
     2       1X,'MAXIMUM BAND WIDTH OF [AL] PLUS 1:',I5)
C
C4------ALLOCATE SPACE FOR THE DE4 ARRAYS.
      ALLOCATE (AU(ID4DIM,MXUP))
      ALLOCATE (IUPPNT(ID4DIM,MXUP))
      ALLOCATE (AL(MXBW,MXLOW))
      ALLOCATE (IEQPNT(NCOL,NROW,NLAY))
      ALLOCATE (D4B(MXEQ))
      ALLOCATE (LRCHDE4(3,ITMX))
      ALLOCATE (HDCGDE4(ITMX))
C
C1------READ IFREQ, MUTD4, ACCL, HCLOSE, AND IPRD4
      ONE=1.
      ZERO=0.
      READ(INDE4,*) IFREQ,MUTD4,ACCLDE4,HCLOSEDE4,IPRD4
      IF(ACCLDE4.LE.ZERO) ACCLDE4=ONE
      IF(IPRD4.LE.0) IPRD4=999
      IF(MUTD4.LT.0 .OR. MUTD4.GT.3) MUTD4=0
      IF(IFREQ.LT.1 .OR. IFREQ.GT.3) THEN
         WRITE(IOUT,11) IFREQ
   11    FORMAT(1X,/1X,'INVALID VALUE FOR IFREQ PARAMETER:',I8)
         CALL USTOP(' ')
      END IF
C
C2------CHECK TO SEE IF THERE IS ITERATION (ITMX>1).
      IF(ITMX.GT.1) THEN
C
C3------THERE IS ITERATION -- DETERMINE TYPE OF ITERATION BASED ON IFREG
C3------VALUE.
         IF(IFREQ.EQ.3) THEN
C
C3A-----EXTERNAL ITERATION -- SET ITERATION VARIABLES AND PRINT A
C3A-----MESSAGE.
            MXITER=ITMX
            NITERDE4=1
            WRITE(IOUT,51)
   51       FORMAT(1X,/21X,'SOLUTION BY D4 DIRECT SOLVER WITH EXTERNAL',
     1       ' ITERATION',/21X,52('-'),/)
C
C3B-----INTERNAL ITERATION -- SET ITERATION VARIABLES AND PRINT A
C3B-----MESSAGE.
         ELSE
            MXITER=1
            NITERDE4=ITMX
            WRITE(IOUT,81)
   81       FORMAT(1X,/21X,'SOLUTION BY D4 DIRECT SOLVER WITH INTERNAL',
     1       ' ITERATION'/21X,52('-')/)
         END IF
C
C3C-----PRINT ITERATION INFORMATION.
         WRITE(IOUT,91) ITMX,ACCLDE4,HCLOSEDE4,IPRD4
   91    FORMAT(1X,'MAXIMUM ITERATIONS =',I3/
     1       1X,'RELAXATION-ACCELERATION PARAMETER =',F10.6/
     2       1X,'HEAD CHANGE CRITERION FOR CLOSURE =',G15.5/
     3       1X,'D4 PRINTOUT INTERVAL =',I4)
         IF(MUTD4.EQ.1) WRITE(IOUT,92)
   92    FORMAT(1X,
     1   'CONVERGENCE PRINTOUT WILL SHOW ONLY THE NUMBER OF ITERATIONS')
         IF(MUTD4.EQ.2) WRITE(IOUT,93)
   93    FORMAT(1X,'CONVERGENCE PRINTOUT WILL BE SUPPRESSED')
C
C4------NO ITERATION -- SET ITERATION VARIABLES AND PRINT A MESSAGE.
      ELSE
         MXITER=1
         NITERDE4=1
         ACCLDE4=ONE
         WRITE(IOUT,94)
   94    FORMAT(1X,/21X,
     1   'SOLUTION BY D4 DIRECT SOLVER WITH NO ITERATION',/21X,46('-')/)
         IF(MUTD4.EQ.2) WRITE(IOUT,95)
   95    FORMAT(1X,'PRINTOUT OF MAXIMUM HEAD CHANGE WILL BE SUPPRESSED')
      END IF
C
C5------PRINT MESSAGE ABOUT FREQUENCY AT WHICH [A] IS ELIMINATED.
      IF(IFREQ.EQ.3) WRITE(IOUT,102)
  102 FORMAT(1X,'NON-LINEAR PROBLEM -- [A] MATRIX ELIMINATED',
     1   ' EVERY TIME EQUATIONS ARE REFORMULATED')
      IF(IFREQ.NE.3) WRITE(IOUT,103) IFREQ
  103 FORMAT(1X,'LINEAR PROBLEM -- [A] MATRIX ELIMINATED ONLY',
     1   ' WHEN IT CHANGES -- IFREQ=',I1)
C
C6------RETURN.
      CALL DE47PSV(IGRID)
      RETURN
      END

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

      SUBROUTINE DE47DA(IGRID)
C  Deallocate DE4 data
      USE DE4MODULE
C
      DEALLOCATE(DE4DAT(IGRID)%MXUP)
      DEALLOCATE(DE4DAT(IGRID)%MXLOW)
      DEALLOCATE(DE4DAT(IGRID)%MXEQ)
      DEALLOCATE(DE4DAT(IGRID)%MXBW)
      DEALLOCATE(DE4DAT(IGRID)%ITMX)
      DEALLOCATE(DE4DAT(IGRID)%ID4DIR)
      DEALLOCATE(DE4DAT(IGRID)%NITERDE4)
      DEALLOCATE(DE4DAT(IGRID)%IFREQ)
      DEALLOCATE(DE4DAT(IGRID)%IPRD4)
      DEALLOCATE(DE4DAT(IGRID)%MUTD4)
      DEALLOCATE(DE4DAT(IGRID)%ID4DIM)
      DEALLOCATE(DE4DAT(IGRID)%NBWL)
      DEALLOCATE(DE4DAT(IGRID)%NUPL)
      DEALLOCATE(DE4DAT(IGRID)%NLOWL)
      DEALLOCATE(DE4DAT(IGRID)%NLOW)
      DEALLOCATE(DE4DAT(IGRID)%NEQ)
      DEALLOCATE(DE4DAT(IGRID)%NUP)
      DEALLOCATE(DE4DAT(IGRID)%NBW)
      DEALLOCATE(DE4DAT(IGRID)%DELTL)
      DEALLOCATE(DE4DAT(IGRID)%ACCLDE4)
      DEALLOCATE(DE4DAT(IGRID)%HCLOSEDE4)
      DEALLOCATE(DE4DAT(IGRID)%IUPPNT)
      DEALLOCATE(DE4DAT(IGRID)%IEQPNT)
      DEALLOCATE(DE4DAT(IGRID)%AU)
      DEALLOCATE(DE4DAT(IGRID)%AL)
      DEALLOCATE(DE4DAT(IGRID)%D4B)
      DEALLOCATE(DE4DAT(IGRID)%HDCGDE4)
      DEALLOCATE(DE4DAT(IGRID)%LRCHDE4)
C
      RETURN
      END

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

      SUBROUTINE DE47PNT(IGRID)
C  Set pointers to DE4 data for grid.
      USE DE4MODULE
C
      MXUP=>DE4DAT(IGRID)%MXUP
      MXLOW=>DE4DAT(IGRID)%MXLOW
      MXEQ=>DE4DAT(IGRID)%MXEQ
      MXBW=>DE4DAT(IGRID)%MXBW
      ITMX=>DE4DAT(IGRID)%ITMX
      ID4DIR=>DE4DAT(IGRID)%ID4DIR
      NITERDE4=>DE4DAT(IGRID)%NITERDE4
      IFREQ=>DE4DAT(IGRID)%IFREQ
      IPRD4=>DE4DAT(IGRID)%IPRD4
      MUTD4=>DE4DAT(IGRID)%MUTD4
      ID4DIM=>DE4DAT(IGRID)%ID4DIM
      NBWL=>DE4DAT(IGRID)%NBWL
      NUPL=>DE4DAT(IGRID)%NUPL
      NLOWL=>DE4DAT(IGRID)%NLOWL
      NLOW=>DE4DAT(IGRID)%NLOW
      NEQ=>DE4DAT(IGRID)%NEQ
      NUP=>DE4DAT(IGRID)%NUP
      NBW=>DE4DAT(IGRID)%NBW
      DELTL=>DE4DAT(IGRID)%DELTL
      ACCLDE4=>DE4DAT(IGRID)%ACCLDE4
      HCLOSEDE4=>DE4DAT(IGRID)%HCLOSEDE4
      IUPPNT=>DE4DAT(IGRID)%IUPPNT
      IEQPNT=>DE4DAT(IGRID)%IEQPNT
      AU=>DE4DAT(IGRID)%AU
      AL=>DE4DAT(IGRID)%AL
      D4B=>DE4DAT(IGRID)%D4B
      HDCGDE4=>DE4DAT(IGRID)%HDCGDE4
      LRCHDE4=>DE4DAT(IGRID)%LRCHDE4
C
      RETURN
      END

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

      SUBROUTINE DE47PSV(IGRID)
C  Save pointers to DE4 data.
      USE DE4MODULE
C
      DE4DAT(IGRID)%MXUP=>MXUP
      DE4DAT(IGRID)%MXLOW=>MXLOW
      DE4DAT(IGRID)%MXEQ=>MXEQ
      DE4DAT(IGRID)%MXBW=>MXBW
      DE4DAT(IGRID)%ITMX=>ITMX
      DE4DAT(IGRID)%ID4DIR=>ID4DIR
      DE4DAT(IGRID)%NITERDE4=>NITERDE4
      DE4DAT(IGRID)%IFREQ=>IFREQ
      DE4DAT(IGRID)%IPRD4=>IPRD4
      DE4DAT(IGRID)%MUTD4=>MUTD4
      DE4DAT(IGRID)%ID4DIM=>ID4DIM
      DE4DAT(IGRID)%NBWL=>NBWL
      DE4DAT(IGRID)%NUPL=>NUPL
      DE4DAT(IGRID)%NLOWL=>NLOWL
      DE4DAT(IGRID)%NLOW=>NLOW
      DE4DAT(IGRID)%NEQ=>NEQ
      DE4DAT(IGRID)%NUP=>NUP
      DE4DAT(IGRID)%NBW=>NBW
      DE4DAT(IGRID)%DELTL=>DELTL
      DE4DAT(IGRID)%ACCLDE4=>ACCLDE4
      DE4DAT(IGRID)%HCLOSEDE4=>HCLOSEDE4
      DE4DAT(IGRID)%IUPPNT=>IUPPNT
      DE4DAT(IGRID)%IEQPNT=>IEQPNT
      DE4DAT(IGRID)%AU=>AU
      DE4DAT(IGRID)%AL=>AL
      DE4DAT(IGRID)%D4B=>D4B
      DE4DAT(IGRID)%HDCGDE4=>HDCGDE4
      DE4DAT(IGRID)%LRCHDE4=>LRCHDE4
C
      RETURN
      END

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