! Lateral groundwater flow according to the Dupuit assumptions.
! This program computes non-steady saturated flow in a strip aquifer of constant thickness D (L). The saturated hydraulic 
! conductivity K (L/T) and the storage coefficient mu are uniform. The aquifer may be leaky, with alpha (< 0; 1/T) and b
! (L/T) characterizing the exchange with the aquifer below. If the lower aquifer has a constant hydraulic head H2 (L),
! -1/alpha (T) is the resistance of the aquitard between the aquifers, and b = -alpha * H2.
! The aquifer has a no-flow boundary at x1 = 0 and a prescribed head at x1 = L (surface water with prescribed head).
! The flow obeys Dupuit's assumptions: strictly horizontal flow, no vertical gradients in hydraulic head H (L), hence 
! any recharge or extraction is instantly divided over the entire aquifer thickness. In case of phreatic flow, variations
! in the groundwater level are assumed to have a negligible effect on aquifer thickness. The open water bordering the
! aquifer at x1 = L extends over the full thickness of the aquifer.
! The water level in this open water body and the recharge R are arbitrary functions of time. The model assumes R to be
! distributed uniformly over the aquifer, and to be constant during the interval between any two observations. The initial
! hydraulic head across the aquifer strip is flexible and specified on input. In case of steady flow, this is irrelevant
! and does not need to be given. In case of a flat groundwater level, H should be specified only at x1 = 0 and x1 = L for
! time zero.
! Open water levels H(L,t) (L), denoted H1 must be specified at t = zero and the desired end of the time period to be
! simulated, and at an arbitrary number of times in between. The forcings (R(t), H(L,t) and H(x1,0) can be specified at
! non-equidistant spacing in time or space, without limitations to the number of observations. 
!
MODULE DoublePrecision
! This module defines the precision of Kind DP. 
  IMPLICIT NONE
  INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(20, 300)  
END MODULE DoublePrecision
!
PROGRAM DupuitFlow
! 
!2345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234
! Interface: All input is list directed. Variables should be separated by arbitrary numbers of spaces. Care must be taken 
! to ensure that all input variables appear in their designated line in the input files. Two character strings that are part
! of the input each appear on a separate line. 
! The units of the variables can be freely chosen, but must be consistent throughout. The input and output files use 
! dimensional values of spatial coordinate x1 and time t.
!
! Input file Dupuitflow.IN
! Line 1: Problem label (<= 60 characters)
! Line 2: will not be read, can be used for labels
! Line 3: K, D, L, mu, alpha (designated a in the referenced literature), b
! Line 4: will not be read, can be used for labels
! Line 5: Character string designated steady flow ('Steady', 'steady', or 'STEADY') or transient flow (any other string).
! Line 6: will not be read, can be used for labels 
! Line 7: Number of locations on the x1-axis for which output is desired (NrOutputLocations), Number of times for which
!         output is desired (NrOutputTimes), the Maximum value of n that will be considered in the infinite sums (Maxn),
!         The number of terms of the infinite sums that have to be calculated between convergernce checks (Interval) and 
!         the tolereance criterion that will determine if the sum has converged satisfactorily and no further terms need
!         to be computed(Tolerance). 
!         Note that NrOutputTimes, Maxn, Interval, and Tolerance are irrelevant when the flow is steady and may be
!         omitted. For transient flow, NrOutputTimes must be at least 2.
! Line 8: will not be read, can be used for labels
! For steady flow:
!   Line 9: The hydraulic head in the surface water (H(L); H1(1,2) in the code), the recharge rate (R; R(1,2)in the code).
!   End of file
! For transient flow:
!   Line 9: The number of locations on the x1 axis where the initial hydraulic head H(x,0) is given (NrLocationsHinitial),
!           the number of times the hydraulic head of the surface water (H(L,t) is given (NrTimesH1), the number of 
!           times the recharge rate R(t) is given (NrTimesR). Note that a given value of R is valid for the time interval
!           between the previous time R was given (or t = 0) and the time for which the current value is specified.
!   End of file.
!
! Input file Hinitial.IN (only required for transient flow)
! The number of lines must be equal to NrLocationsHinitial. Each line contains two real numbers. The first is a location 
! on the x1-axis, the second is the value of the hydraulic head at that location at t = 0. (x1, H(x1,0))
!
! Input file H1.IN (only required for transient flow)
! The number of lines must be equal to NrTimesH1. Each line contains two real numbers. The first is a time, the second is
! the value of the hydraulic head in the surface water at that time. (t, H(L,t)).
!
! Input file R.IN (only required for transient flow)
! The number of lines must be equal to NrTimesR. Each line contains two real numbers. The first is a time, the second is
! the value of the recharge rate for the time period between the times indicated on the line and the line above (t, R(t)).
! For the first line the infiltration rate is valid from time zero to the indicated time.
!
! Input file OutputTimes.IN (only required for transient flow)
! The number of lines must be equal to NrOutputTimes. Each line contains a single real number: a time for which output is 
! desired. 
!
! Input file OutputLocations.IN
! The number of lines must be equal to NrOutputLocations. Each line contains a single real number: a location along the
! x1-axis for which output is desired. 
!
! On output, selected input data are echoed to file DupuitFlow.OUT (sometimes after corrections or slight processing), and
! some information regarding the computations is written to file 'Diagnostics.OUT'. The results of the calculations are the 
! hydraulic head H, the flux to/from the aquifer Q, and other aquifer-scale properties: average head Havg, aquifer-scale
! hydraulic conductivity Kup, and a perturbation measure P. H is given in file H.OUT at the desired times and locations.
! For the same times, the flux per bank length Q to (+) or from (-)the open water, Havg, and the resulting
! upscaled hydraulic conductivity Kup and perturbation measure P are written to file AquiferScale.OUT. Perturbation measure P
! indicates how far the aquifer is from a relaxed state is also part of the output. This value is the normalized difference
! between the actual and asymptotic value of the upscaled hydraulic conductivity appropriate for the conditions at any
! given time.
! An aquifer can have three asymptotic aquifer-scale hydraulic conductivities and three associated characteristic times. All
! of these are given in file Dupuitflow.OUT - the user needs to select the appropriate values for further analysis or
! interpretation. 
!
! Of the output files, only Diagnostics.OUT will be overwritten by a subsequent run of Dupuitflow (Status = 'Unknown').
! The presence of the other files in the same folder as the executable will cause the program to crash (Status = 'New)'.
!
! References
! de Rooij, G.H. 2012. Transient flow between aquifers and surface water: analytically derived field-scale hydraulic heads
! and fluxes. Hydrol. Earth Syst. Sci., 16: 649-669. doi:10.5194/hess-16-649-2012.
!
! Author: Gerrit H. de Rooij
! June 2013, updated February 2016
!
USE DoublePrecision
! Variable declarations
IMPLICIT NONE
CHARACTER (LEN = 60) :: ProblemLabel
CHARACTER (LEN = 6) :: FlowCase
INTEGER :: I, J, n, Counter, CounterH1, CounterR, NrLocationsHinitial, NrtimesH1, NrtimesR, NumChanges, Maxn, &
           Interval, NrOutputTimes, NrOutputLocations, Start
INTEGER, Dimension (:), Allocatable :: nmax
REAL(KIND = DP) :: K, D, mu, L, alpha, A, b, BCInt, t, Term1, Term2, Term3, Help1, Help2, Plusmin, Pin, nreal, Tolerance, &
          First, KupRconstant, KupZeroRecharge, KupLeakyAquifer,CharTimeRconstant,CharTimeZeroRecharge, CharTimeLeakyAquifer
REAL(KIND = DP), PARAMETER :: PI = 3.14159265358979323846
REAL(KIND = DP), Dimension (:), Allocatable :: ICInt, OutputTimes, DimensionlessOutputTimes, OutputLocations, Havg, &
                                               H1OutputTimes, Q, Kup, P                                                
REAL(KIND = DP), Dimension (:,:), Allocatable :: Hinitial, Hinitialspare, H1, R, ChangeDrivers, H, Mn
!
! Provide an initial value for FlowCase. If it is missing from the input file, the segment of code determining FlowCase still 
! has a starting value that allows the program to execute.
FlowCase = 'Generl'
! Read input
! Data that need to be provided in the input files:
! An indicator that the flow is steady by a character string 'Steady'; in this case an initial condition needs not be 
! specified. N.B. This is only allowed if alpha = 0 because the groundwater level is assumed to be parabolic.
! K, D, L, a (alpha), b, mu; Array with H(x,0) and its observation locations; array with observations of R and the times 
! of observation (R); array with observations of H(1,t) and the times of observation (H1); array with times at which output
! is desired. For all these arrays, the number of rows also needs to be explicitly given.
OPEN (30, FILE = 'Diagnostics.OUT', STATUS = 'UNKNOWN')
OPEN (21, FILE = 'DupuitFlow.IN', STATUS = 'OLD', ACTION = 'READ')
READ (21, "(A)") ProblemLabel
READ (21, *)
READ (21, *) K, D, L, mu, alpha, b
READ (21, "(/,A)") FlowCase
IF ((Flowcase == 'Steady') .OR. (FlowCase == 'STEADY') .OR. (FlowCase == 'steady')) FlowCase = 'Steady' 
READ (21, *)
IF (Flowcase == 'Steady') THEN
  READ (21, *) NrOutputLocations
ELSE
  READ (21, *) NrOutputLocations, NrOutputTimes, Maxn, Interval, Tolerance
END IF
READ (21, *)
WRITE (30, 600) ProblemLabel
600 FORMAT (1X, A)
! Read the required array sizes and allocate the arrays with external forcings accordingly.
! When Flowcase = 'Steady', the values of R and H1 must be specified. The files with records of surface water levels and recharge
! will not be required then, and neither is the file with desired output times.
IF (FlowCase == 'Steady') THEN
  Allocate (H1(1:1, 1:2))
  Allocate (R(1:1, 1:2))
  Allocate (OutputTimes(1:1))
  Allocate (DimensionlessOutputTimes(1:1))
  READ (21, *) H1(1,2), R(1,2)
  H1(1,1) = 0.0_DP
  R(1,1) = 0.0_DP
  NrTimesH1 = 1
  NrTimesR = 1
  NrOutputTimes = 1
  OutputTimes(1) = 0.0_DP
  WRITE (30, 605) FlowCase 
ELSE
  READ (21, *) NrLocationsHinitial, NrTimesH1, NrTimesR
  Allocate (Hinitial(1:NrLocationsHinitial, 1:2))
  Allocate (H1(1:NrTimesH1, 1:2))
  Allocate (R(1:NrTimesR, 1:2))
  Allocate (OutputTimes(1:NrOutputTimes))
  Allocate (DimensionlessOutputTimes(1:NrOutputTimes))
  Allocate (nmax(1:NrOutputTimes))
END IF
605 FORMAT (1X,' FlowCase: ', A)
CLOSE (21)
! The arrays with recharge, surface water levels, the initial groundwater level, and the desired output times and locations 
! are given in separate files.

! Allocate the remaining input and output arrays
Allocate (OutputLocations(1:NrOutputLocations))
Allocate (ChangeDrivers(1:NrtimesH1 + NrtimesR, 1:3)) 
Allocate (ICInt(0:Maxn))
Allocate (Mn(1:NrOutputTimes, 0:Maxn))
Allocate (H(1:NrOutputLocations, 1:NrOutputTimes))
Allocate (H1OutputTimes(1:NrOutputTimes))
Allocate (Havg(1:NrOutputTimes))
Allocate (Q(1:NrOutputTimes))
Allocate (Kup(1:NrOutputTimes))
Allocate (P(1:NrOutputTimes))
! Read input files with arrays of observations, and desired times and locations for output.
IF (FlowCase /= 'Steady') THEN
  OPEN (21, FILE = 'Hinitial.IN', STATUS = 'OLD', ACTION = 'READ')
    DO I = 1, NrLocationsHinitial
      READ (21, *) Hinitial(I,1), Hinitial(I,2)
    END DO
  CLOSE (21)
  OPEN (21, FILE = 'H1.IN', STATUS = 'OLD', ACTION = 'READ')
    DO I = 1, NrTimesH1
      READ (21, *) H1(I,1), H1(I,2)
    END DO
  CLOSE (21)
  OPEN (21, FILE = 'R.IN', STATUS = 'OLD', ACTION = 'READ')
    DO I = 1, NrTimesR
      READ (21, *) R(I,1), R(I,2)
    END DO
  CLOSE (21)
  OPEN (21, FILE = 'OutputTimes.IN', STATUS = 'OLD', ACTION = 'READ')
    DO I = 1, NrOutputTimes
      READ (21, *) OutputTimes(I)
    END DO
  CLOSE (21)
END IF
OPEN (21, FILE = 'OutputLocations.IN', STATUS = 'OLD', ACTION = 'READ')
DO I = 1, NrOutputLocations
  READ (21, *) OutputLocations(I)
END DO
CLOSE (21)
! Ensure that the initial hydraulic head is read for x = 0 and x = 1 (dimensionless). If this is not the case, add
! these observation points and assign to them the value of the nearest observation provided in the input. Next, shift
! observation points to x = 0 and x = 1 in case they are within the domain instead of exactly at the domain limits.
Nonsteady1: IF (FlowCase /= 'Steady') THEN
  Hi1: IF ((Hinitial(1,1) > 0.0_DP) .OR. (Hinitial(NrLocationsHinitial,1) < L)) THEN
    Allocate (Hinitialspare(1:NrLocationsHinitial, 1:2))
    Hinitialspare = Hinitial
    Deallocate (Hinitial)
    Hi2: IF ((Hinitialspare(1,1) > 0.0_DP) .AND. (Hinitialspare(NrLocationsHinitial,1) < L)) THEN
      Allocate (Hinitial(1:NrLocationsHinitial+2, 1:2))
      Hinitial(1,1) = 0.0_DP
      Hinitial(NrLocationsHinitial+2,1) = L
      Hinitial(1,2) = Hinitialspare(1,2)
      Hinitial(NrLocationsHinitial+2,2) = Hinitialspare(NrLocationsHinitial,2)
      Hinitial(2:NrLocationsHinitial+1,1:2) = Hinitialspare(1:NrLocationsHinitial,1:2)
      NrLocationsHinitial = NrLocationsHinitial + 2 
      WRITE (30, 610) NrLocationsHinitial
    ELSE
      Allocate (Hinitial(1:NrLocationsHinitial+1, 1:2))
      Hi3: IF (Hinitialspare(1,1) > 0.0_DP) THEN
        Hinitial(1,1) = 0.0_DP
        Hinitial(1,2) = Hinitialspare(1,2)
        Hinitial(2:NrLocationsHinitial+1,1:2) = Hinitialspare(1:NrLocationsHinitial,1:2)
        NrLocationsHinitial = NrLocationsHinitial + 1
        WRITE (30, 615) NrLocationsHinitial
      ELSE
        Hinitial(NrLocationsHinitial+1,1) = L
        Hinitial(NrLocationsHinitial+1,2) = Hinitialspare(NrLocationsHinitial,2)
        Hinitial(1:NrLocationsHinitial,1:2) = Hinitialspare(1:NrLocationsHinitial,1:2)
        NrLocationsHinitial = NrLocationsHinitial + 1
        WRITE (30, 620) NrLocationsHinitial
      END IF Hi3
    END IF Hi2
    Deallocate (Hinitialspare)
  END IF Hi1
! The next two statements shift observation points outside the domain to x = 0 and x = L (e.g., in case of roundoff
! error.)
  IF (Hinitial(1,1) < 0.0_DP) THEN
    Hinitial(1,1) = 0.0_DP
    WRITE (30, 625) 
  END IF
  IF (Hinitial(NrLocationsHinitial,1) > L) THEN
    Hinitial(NrLocationsHinitial,1) = L
    WRITE (30, 630)
  END IF
END IF Nonsteady1
610 FORMAT (' H(x,0) was given for neither x = 0 nor for x = L; added these values and updated NrLocationsHinitial to',I5)
615 FORMAT (' H(x,0) was not given for x = 0; added its value and updated NrLocationsHinitial to',I5)
620 FORMAT (' H(x,0) was not given for x = L; added its value and updated NrLocationsHinitial to',I5)
625 FORMAT (' Initial value for x = 0 out of domain (x < 0); shifted to x = 0.')
630 FORMAT (' Initial value for x = L out of domain (x > L); shifted to x = L.')
! Make sure alpha is negative (physical requirement), and calculate A. 
alpha = -ABS(alpha)
A = alpha * L * L /(mu * K * D)
! Make sure the maximum number of terms to be evaluated is at least equal to Interval. This ensures proper execution of the 
! convergence checks, which operate with intervals with the number of terms set by Interval.
IF (Interval < 1) Interval = 1
IF (Maxn <= Interval - 1) Maxn = Interval - 1
OPEN (21, FILE = 'DupuitFlow.OUT', STATUS = 'NEW', ACTION = 'WRITE')
! Check if the flow is steady. If not, check if the initial groundwater is level. Based on the outcome, define variable
! FlowCase.
IF ((FlowCase == 'Steady') .AND. (alpha /= 0.0_DP)) THEN
! Write an error message and exit the program.
  WRITE (21, 300) ProblemLabel
  GO TO 100
END IF
IF (FlowCase /= 'Steady') THEN
! Check if the groundwater is initially level
  IF ((NrLocationsHinitial == 2) .AND. (Hinitial(1,2) == Hinitial(2,2))) THEN 
    FlowCase = 'FlatGL'
  ELSE
    FlowCase = 'Generl'  
  END IF 
END IF

! Echo input to the first output file. 
WRITE (21, 305) ProblemLabel
IF (FlowCase == 'Steady') THEN
  WRITE (21, 310) K, D, L, mu, alpha, b, A
ELSE
  WRITE (21, 310) K, D, L, mu, alpha, b, A, Maxn + 1, Interval, Tolerance
END IF
WRITE (21, 315) FlowCase
IF (FlowCase == 'Steady') THEN
  WRITE (21, 320) OutputLocations(1), OutputLocations(NrOutputLocations)
ELSE
  WRITE (21, 325) NrLocationsHinitial, Hinitial(1,1), Hinitial(NrLocationsHinitial,1)
  WRITE (21, 330) NrTimesH1, H1(1,1), H1(NrTimesH1,1)
  WRITE (21, 335) NrTimesR, R(1,1), R(NrTimesR,1)
  WRITE (21, 320) OutputLocations(1), OutputLocations(NrOutputLocations)
  WRITE (21, 340) OutputTimes(1), OutputTimes(NrOutputTimes)
END IF
300 FORMAT(1X, A/, ' Steady flow but aquifer leaky.')
305 FORMAT (' Output for ', A)
310 FORMAT(/,'          K          D          L         mu      alpha          b          A',/,&
           1X, F10.4,' ',F10.4,' ',F10.4,' ',F10.4,' ',F10.4,' ',F10.4,' ',F10.4/&
           ' Evaluated sum terms  Interval  Tolerance',/, 15X, I5, 5X, I5,' ',F10.4) 
315 FORMAT(/,' General flow problem, initially flat groundwater level, or steady flow: ', A)
320 FORMAT(/,' First and last location for which output will be generated: ', F10.4, ', ', F10.4)
325 FORMAT(/,' Number of locations for which H(x,0) is given:       ',I5,/,&
             ' First and last of those locations: ',F10.4,', ',F10.4) 
330 FORMAT(/,' Number of times for which H(1,t) is given:           ',I5,/,&
             ' First and last of those times:     ',F10.4,', ',F10.4)
335 FORMAT(/,' Number of times for which R(t) is given:             ',I5,/,&
             ' First and last of those times:     ',F10.4,', ',F10.4) 
340 FORMAT(' First and last time for which output will be generated:     ', F10.4,', ', F10.4)       
! Segment of the code that creates an array of observation times for recharge R(t) and surface water level H(1,t), along
! with the values at these times.
! Initialize various counters.
CounterH1 = 1
CounterR = 1
I = 1
! Convert the observation times in the arrays containing them to dimensionless times. The time coordinate is set to zero at
! the time of the first observation of H1(t).
R(:,1) = K * D * (R(:,1) - H1(1,1))/(L * L)
DimensionlessOutputTimes = K * D * (OutputTimes - H1(1,1)) / (L*L)
H1(:,1) = H1(:,1) - H1(1,1)
H1(:,1) = K * D * H1(:,1) / (L * L)
! Start an unbounded DO loop to sort the times at which the external forcings are observed.
SortForc: DO
  StopLoop: IF((CounterH1 == NrtimesH1) .AND. (CounterR == NrtimesR)) THEN
!   Complete writing times and values to ChangeDrivers, then exit the DO-loop. The stopping criteria allow at most
!   one observation of R after t = H1(1,NrtimesH1).
    SL2: IF (H1(CounterH1,1) == R(CounterR,1)) THEN !H1 and R observed at the same time.
      ChangeDrivers(I,1) = H1(CounterH1,1)
      ChangeDrivers(I,2) = H1(CounterH1,2)
      ChangeDrivers(I,3) = R(CounterR,2)
      EXIT
    ELSE
      SL3: IF (H1(CounterH1,1) < R(CounterR,1)) THEN !The final value of H1 was observed before that of R.
        ChangeDrivers(I,1) = H1(CounterH1,1)
        ChangeDrivers(I,2) = H1(CounterH1,2)
        ChangeDrivers(I,3) = R(CounterR,2) !An intermediate value of R is equal to the observed value at the
!                                           end of its time interval.
        EXIT
!       Note that it is not necessary to examine what happens after the time the final H1 value was observed,
!       since this lies beyond the time period for which the calculations are performed.

      ELSE !The final observation of R occurs before t=H1(1,NrtimesH1). This is undesirable, but if it occurs,
!           assume R=0 between the final observation and t=H1(1,NrtimesH1).
        ChangeDrivers(I,1) = R(CounterR,1)
!       Interpolate H1 between the measurements before and after the measurement of R.
        ChangeDrivers(I,2) = H1(CounterH1 - 1,2)+((H1(CounterH1,2) - H1(CounterH1 - 1,2))* &
                             (ChangeDrivers(I,1) - H1(CounterH1 - 1,1))/(H1(CounterH1,1) - H1(CounterH1 - 1,1)))
        ChangeDrivers(I,3) = R(CounterR,2)
        ChangeDrivers(I+1,1) = H1(CounterH1,1)
        ChangeDrivers(I+1,2) = H1(CounterH1,2)
        ChangeDrivers(I+1,3) = 0.0_DP
        I = I + 1
        EXIT 
      END IF SL3
    END IF SL2
  END IF StopLoop
! R may be observed after t=0; the value is then assumed to apply to the period between t=0 and the observation time.
! The following IF blocks fill ChangeDrivers one observation time at a time, starting with the first observation at t = 0.
  SortObs: IF (H1(CounterH1,1) == R(CounterR,1)) THEN !H1 and R observed at the same time.            
    ChangeDrivers(I,1) = H1(CounterH1,1)
    ChangeDrivers(I,2) = H1(CounterH1,2)
    ChangeDrivers(I,3) = R(CounterR,2)
    CounterH1 = CounterH1 + 1
    CounterR = CounterR + 1
    I=I+1
  ELSE
    SO2:IF (H1(CounterH1,1) < R(CounterR,1)) THEN !H1 was observed before R.
      ChangeDrivers(I,1) = H1(CounterH1,1)
      ChangeDrivers(I,2) = H1(CounterH1,2)
      ChangeDrivers(I,3) = R(CounterR,2) !This is the next observed value of R.
      CounterH1 = CounterH1 + 1
      I=I+1
    ELSE !R was observed before H1
      ChangeDrivers(I,1) = R(CounterR,1)
!     Interpolate H1 between the previous and next observation of H1 to estimate the value at the time R was
!     observed.
      ChangeDrivers(I,2) = H1(CounterH1 - 1,2)+((H1(CounterH1,2)-H1(CounterH1 - 1,2))* &
                           (Changedrivers(I,1)-H1(CounterH1 - 1,1))/(H1(CounterH1,1)-H1(CounterH1 - 1,1)))
      ChangeDrivers(I,3) = R(CounterR,2) 
      CounterR = CounterR + 1      
      I=I+1
    END IF SO2
  END IF SortObs          
END DO SortForc
NumChanges = I
WRITE (30, 635) I
635 FORMAT (/, ' Changedrivers completed; Number of times H(L,t) and/or R change: ',I5)

! The full solution contains two integrals. 
! The next section calculates the integral that incorporates the effect of the initial phreatic surface. Depending on the
! value of FlowCase, the general formulation is chosen, or the special case of steady flow in a non-leaky aquifer, or the
! case of an initially uniform hydraulic head. The results are stored in array ICInt.
IF (FlowCase == 'Steady') THEN
! The integral involving the initial condition can be calculated for steady flow, but is not required for the solution as H(x),
! the average H, and Q can be readily calculated in explicit form without resorting to infinite sums. The code for calculating
! the sum terms is given for completeness, but commented out.
!
! Plusmin = -1.0_DP
! ICInt = R(1,2) * L * L / (2.0_DP * K * D)
! DO n = 0, Maxn
!   nreal = REAL(n, DP)
!   Plusmin = -1.0_DP * Plusmin
!   Pin = (nreal + 0.5_DP) * PI
!   ICInt(n) = Plusmin * ICInt(n) * (1.0_DP - (((Pin*Pin) - 2.0_DP) / (Pin*Pin))) / Pin
! END DO
ELSE IF (FlowCase == 'FlatGL') THEN
  Plusmin = -1.0_DP
  Help1 = (Hinitial(1,2) - H1(1,2)) / PI
  DO n = 0, Maxn
    nreal = REAL(n, DP)
    Plusmin = -1.0_DP * Plusmin
    ICInt(n) = Plusmin * Help1 / (nreal + 0.5_DP)
  END DO
ELSE
! Implement the general solution of the integral over x. Start with making the spatial coordinate dimensionless.
  Plusmin = -1.0_DP
  Hinitial(:,1) = Hinitial(:,1)/L
  DO n = 0, Maxn
    nreal = REAL(n, DP)
    Plusmin = -1.0_DP * Plusmin
    Pin = (nreal + 0.5_DP) * PI
    Term2 = 0.0_DP
    Term3 = 0.0_DP
    DO I = 1, NrLocationsHinitial - 1
      Term2 = Term2 + (Hinitial(I,2) - (Hinitial(I,1) * (Hinitial(I+1,2)-Hinitial(I,2))/(Hinitial(I+1,1)-Hinitial(I,1))))* &
              (sin(Pin * Hinitial(I+1,1)) - sin(Pin * Hinitial(I,1)))
      Help1 = (Hinitial(I+1,1) * sin(Pin * Hinitial(I+1,1))) - (Hinitial(I,1) * sin(Pin * Hinitial(I,1)))
      Help2 = (cos(Pin * Hinitial(I+1,1)) - cos(Pin * Hinitial(I,1))) / Pin
      Term3 = Term3 + ((Hinitial(I+1,2)-Hinitial(I,2))/(Hinitial(I+1,1)-Hinitial(I,1))) * (Help1 + Help2) / Pin
    END DO
    ICInt(n) = (Term2 / Pin) + Term3 - (H1(1,2) * Plusmin / Pin)
  END DO
END IF

Nonsteady2: IF (FlowCase /= 'Steady') THEN
! Calculate Mn and Q for all output times. The infinite sum converges slowest for Q, therefore it is calculated before H.
! The number of terms required to reach convergence is retained. The same number of terms is later used to calculate the 
! hydraulic head in space and time (H), and the spatial average of H (Havg).
!
! Initialize Q.
  Q = 0.0_DP
! Calculate Mn and Q. For t = 0, the calculation of Mn and Q is easier and carried out separately.
  Start = 1
  tzero: IF (DimensionlessOutputTimes(1) <= 0.0_DP) THEN
    Start = 2
    Mn(1,:) = PI * ICInt(:)
    Plusmin = -1.0_DP
    DO n = 0, Maxn
      nreal = REAL(n, DP)
      Plusmin = -1.0_DP * Plusmin
      Q(1) = Q(1) + ((2.0_DP * K * D / L) * (nreal + 0.5_DP) * Mn(1,n) / Plusmin)
!     Convergence check. Avoid dividing by zero.
      Converge1: IF (Q(1) /= 0.0_DP) THEN
        IF (n == 0) THEN
          First = Q(1)
          J = -1
        ELSE
          IF ((n - J) == Interval) THEN
            IF (ABS((Q(1) - First) / Q(1)) > Tolerance) THEN ! No convergence yet. Continue the summation.
              First = Q(1)
              J = n
            ELSE ! Convergence achieved. Store the current value of n and leave the loop.
              nmax(1) = n
              EXIT
            END IF
          END IF
        END IF
      END IF Converge1
    END DO
    IF (n == Maxn + 1) THEN ! The sum did not converge, write an error message.
      WRITE (30, 640) DimensionlessOutputTimes(1) 
    END IF
  END IF tzero
640 FORMAT (' The infinite sum for Q did not converge for dimensionless time ',F10.4)
! Calculate Mn and Q for t > 0.
  Outputtime: DO I = Start, NrOutputTimes
!   Make the time coordinate dimensionless
    t = DimensionlessOutputTimes(I)
    Help1 = EXP(alpha * L * L * t / (mu * K * D))
    Plusmin = -1.0_DP
    Sum: DO n = 0, Maxn
      nreal = REAL(n, DP)
      Plusmin = -1.0_DP * Plusmin
      Pin = (nreal + 0.5_DP) * PI
      Term1 = PI * Help1 * EXP(-Pin * Pin * t / mu) * ICInt(n)
      CALL BCIntegral(BCInt, t, n, ChangeDrivers, NumChanges, H1, K, D, L, mu, A, b, PI)
      Mn(I,n) = Term1 + (Help1 * Plusmin * BCInt / (nreal + 0.5_DP))
      Q(I) = Q(I) + ((2.0_DP * K * D / L) * (nreal + 0.5_DP) * Mn(I,n) / Plusmin)
!     Convergence check. Avoid dividing by zero.
      Converge2: IF (Q(I) /= 0.0_DP) THEN
        IF (n == 0) THEN
          First = Q(I)
          J = -1
        ELSE
          IF ((n - J) == Interval) THEN
            IF (ABS((Q(I) - First) / Q(I)) > Tolerance) THEN ! No convergence yet. Continue the summation.
              First = Q(I)
              J = n
            ELSE ! Convergence achieved. Store the current value of n and leave the loop.
              nmax(I) = n
              EXIT Sum
            END IF
          END IF
        END IF
      END IF Converge2
    END DO Sum
    IF (n == Maxn + 1) THEN ! The sum did not converge, write an error message.
      WRITE (30, 640) t 
    END IF
  END DO Outputtime
END IF Nonsteady2

! Calculate H(x,t)
Steady1: IF (FlowCase == 'Steady') THEN
  DO J = 1, NrOutputLocations
    H(J,1) = H1(1,2) + (R(1,2) * L * L * (1.0_DP - ((OutputLocations(J) / L)**2.0_DP)) / (2.0_DP * K * D))
  END DO
ELSE
  Hxt: DO I = 1, NrOutputTimes
!   Initialize H and Havg with the hydraulic head in the surface water at the current time.
    Hxt2a: DO J = 1, NrTimesH1 - 1
      IF ((DimensionlessOutputTimes(I) >= H1(J,1)) .AND. (DimensionlessOutputTimes(I) <= H1(J+1,1))) THEN
!       Interpolate between the two values of the surface water level.
        H1OutputTimes(I) = H1(J,2) + ((H1(J+1,2) - H1(J,2)) * (DimensionlessOutputTimes(I) - H1(J,1)) / (H1(J+1,1) - H1(J,1)))
        Havg(I) = H1OutputTimes(I)
        H(:,I) = H1OutputTimes(I)
        EXIT
      END IF      
    END DO Hxt2a
!   Calculate H(x,t) by adding the infinite sum in its expression to H1(t)
    Hxt2b: DO J = 1, NrOutputlocations
      Hxt3: DO n = 0, nmax(I)
        nreal = REAL(n, DP)
        H(J,I) = H(J,I) + ((2.0_DP / PI) * cos((nreal + 0.5_DP) * PI * OutputLocations(J) / L) * Mn(I,n))
      END DO Hxt3
    END DO Hxt2b  
  END DO Hxt
END IF Steady1
! Calculate Havg and Kup, as well as Q in case of steady flow.
DO I = 1, NrOutputTimes
  IF (FlowCase == 'Steady') THEN
    Havg(I) = (R(1,2) * L * L / (3.0_DP * K * D)) + H1(1,2)
    Q(I) = R(1,2) * L
    Kup(I) = Q(I) / (Havg(I) - H1(1,2))
  ELSE
    Plusmin = -1.0_DP
    DO n = 0, nmax(I)
      nreal = 1.0_DP * n
      Plusmin = -1.0_DP * Plusmin
      Havg(I) = Havg(I) + (2.0_DP * Plusmin * Mn(I,n) / ((nreal + 0.5_DP) * PI * PI))
    END DO  
    Kup(I) = Q(I) / (Havg(I) - H1OutputTimes(I)) 
  END IF
END DO
! Calculate applicable characteristic times for the linear reservoir approximation: the asymptotic values of Kup and the
! characteristic times. These are calculated for all configurations (leaky vs. non leaky, and for non-leaky: constant R 
! and zero R). The values that do not apply to the conditions given should be used for comparison only as they do not characterize
! the aquifer.
KupRconstant = 3.0_DP * K * D / L
KupZeroRecharge = PI * PI * K * D / (4.0_DP * L)
Help1 = 0.0_DP
Help2 = 0.0_DP
Kaquifer: IF (alpha == 0.0_DP) THEN
  ! For alpha = 0, the series in the expression for the asymptotic upscaled conductivity converge to make the expression equivalent
  ! to that for constant,non-zero recharge.
  KupLeakyAquifer = KupRconstant
ELSE
  DO n = 0, Maxn
    nreal = REAL(n, DP)
    Pin = ((nreal + 0.5_DP) * PI)**2.0_DP
    Help1 = Help1 + (1.0_DP / (Pin - (A * mu)))
    Help2 = Help2 + (1.0_DP / ((Pin * ((nreal + 0.5_DP)**2.0_DP)) - (((nreal + 0.5_DP)**2.0_DP) * (A * mu))))
!   Convergence check. Avoid dividing by zero.
    Converge3: IF (Help1 /= 0.0_DP) THEN
      IF (n == 0) THEN
        First = Help1 / Help2
        J = -1
      ELSE
        IF ((n - J) == Interval) THEN
          IF (ABS(((Help1/Help2) - First) / (Help1/Help2)) > Tolerance) THEN ! No convergence yet. Continue the summation.
            First = Help1 / Help2
            J = n
          ELSE ! Convergence achieved. Leave the loop.
            EXIT
          END IF
        END IF
      END IF
    END IF Converge3
  END DO
  IF (n == Maxn + 1) WRITE (30, 650) !The sum did not converge, write an error message.
  KupLeakyAquifer = PI * PI * K * D * Help1 / (Help2 * L)
END IF Kaquifer
650 FORMAT (' The infinite sum for KupLeakyAquifer did not converge.')
CharTimeRconstant = mu / ((KupRconstant / L) - alpha)
CharTimeZeroRecharge = mu / ((KupZeroRecharge / L) - alpha)
CharTimeLeakyAquifer = mu / ((KupLeakyAquifer / L) - alpha)
! Calculate the perturbation measure P indicating how far removed from asymptpotic behavior the aquifer is.
! First determine which conditions apply.
Steady2: IF (FlowCase == 'Steady') THEN
  P(1) = (Kup(1) - KupRconstant) / KupRconstant
ELSE
  Pert1: IF (alpha /= 0.0_DP) THEN
    DO I = 1, NrOutputTimes
      P(I) = (Kup(I) - KupLeakyAquifer) / KupLeakyAquifer
    END DO
  ELSE
    Pert2: DO I = 1, NrOutputTimes
!     Determine whether or not R = 0 at the given time, and use the corresponding asymptotic value of Kup.
      Pert3: DO J = 1, NrTimesR -1
        Pert4: IF ((R(J,1) < DimensionlessOutputTimes(I)) .AND. (R(J+1,1) >= DimensionlessOutputTimes(I))) THEN
          Pert5: IF (R(J+1,2) == 0.0_DP) THEN
            P(I) = (Kup(I) - KupZeroRecharge) / KupZeroRecharge
          ELSE
            P(I) = (Kup(I) - KupRconstant) / KupRconstant
          END IF Pert5
        END IF Pert4
      END DO Pert3
    END DO Pert2
  END IF Pert1
END IF Steady2
! Write output.
OPEN (22, FILE = 'H.OUT', STATUS = 'NEW', ACTION = 'WRITE')
WRITE (22, 400) ProblemLabel
DO J = 1, NrOutputTimes
  DO I = 1, NrOutputLocations
    WRITE (22, 405) OutputTimes(J), OutputLocations (I), H(I,J)
  END DO
END DO
CLOSE (22)
OPEN (22, FILE = 'AquiferScale.OUT', STATUS = 'NEW', ACTION = 'WRITE')
WRITE (22, 410) ProblemLabel
DO I = 1, NrOutputTimes
  WRITE (22, 415) Outputtimes(I), Havg(I), Q(I), Kup(I), P(I)
END DO
CLOSE (22)
WRITE (21, 345) KupZeroRecharge, KupRconstant, KupLeakyAquifer
WRITE (21, 350) CharTimeZeroRecharge, CharTimeRconstant, CharTimeLeakyAquifer
CLOSE (21)
400 FORMAT(' Hydraulic head for ',A,/,'         time     x-coord.            H')
405 FORMAT(1X, ES12.4, 1X, ES12.4, 1X, ES12.4)
410 FORMAT(' Aquifer-scale variables for ',A,/,'         Time    Average H            Q          Kup            P')
415 FORMAT(1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4)
345 FORMAT(//' Asymptotic aquifer-scale hydraulic conductivities',/,' Zero recharge    R constant     Leaky aquifer',/,&
           1X, ES12.4, 3X, ES12.4, 3X, ES12.4)
350 FORMAT(//' Characteristic times of the aquifer',/,' Zero recharge    R constant     Leaky aquifer',/,&
           1X,ES12.4, 3X, ES12.4, 3X, ES12.4)
! Deallocate arrays
IF (FlowCase /= 'Steady') THEN 
  Deallocate (Hinitial)
  Deallocate (nmax)
END IF
100 Deallocate (H1)
Deallocate (R)
Deallocate (OutputTimes)
Deallocate (DimensionlessOutputTimes)
Deallocate (OutputLocations)
Deallocate (ChangeDrivers) 
Deallocate (ICInt)
Deallocate (Mn)
Deallocate (H)
Deallocate (H1OutputTimes)
Deallocate (Havg)
Deallocate (Q)
Deallocate (Kup)
Deallocate (P)
CLOSE (30)
! Subroutine declaration.
CONTAINS
SUBROUTINE BCIntegral(BCInt, t, n, ChangeDrivers, NumChanges, H1, K, D, L, mu, A, b, PI)
! BCIntegral calculates the integral involving the boundary conditions and external forcings for a given dimensionless time t
! and n (the index for the term in the infinite sum in which the integral appears).
!
! Interface: All required input is supplied through the dummy argument list, with the names of the variables identical to
! those in the main program. The value of the integral is returned to the main program in variable BCInt in the argument list.
USE DoublePrecision
IMPLICIT NONE
INTEGER :: J, n, NumChanges
REAL(KIND = DP), Intent(in) :: K, D, mu, L, A, b, t, PI
REAL (KIND = DP) :: Gamman, GamnA, ExpGamt, ExpGamAt, ExpGamAti, ExpGamAtipl1, t1, H1t1, dHdt1, &
                    Term1, Term2, Term3, Term4, Term5, Term6, Help1, H1i, H1ipl1, ti, tipl1, Ripl1, nreal
REAL(KIND = DP), Intent(out) :: BCInt                                                               
REAL(KIND = DP), Dimension (:,:), Allocatable, Intent(in) :: H1, ChangeDrivers
!
! For the given value for n and t, calculate the integral wiht the boundary conditions and the external forcings. For t = 0,
! its value is simply zero.
tzero: IF (t == 0.0_DP) THEN
  BCInt = 0.0_DP
ELSE
  dHdt1 = (ChangeDrivers(2,2) - ChangeDrivers(1,2)) / (ChangeDrivers(2,1) - ChangeDrivers(1,1))
  IF (t <= ChangeDrivers(2,1)) THEN
    t1 = t
!   Find H(1,t) by linear interpolation if t falls within the first observation interval in ChangeDrivers.
    H1t1 = H1(1,2) + ((ChangeDrivers(2,2) - H1(1,2)) * (t - ChangeDrivers(1,1)) / (ChangeDrivers(2,1) - ChangeDrivers(1,1)))
  ELSE
    t1 = ChangeDrivers(2,1)
    H1t1 = ChangeDrivers(2,2)
  END IF
! Determine in which interval of ChangeDrivers the current time falls.
  DO J = 1, NumChanges - 1
    IF ((t >= ChangeDrivers(J,1)) .AND. (t <= ChangeDrivers (J+1,1))) EXIT
  END DO
! Calculate the terms that only depend on n and either the time of interest (t), or the end of the first interval in
! ChangeDrivers, whichever occurs first.
  nreal = REAL(n, DP)
  Gamman = -(((nreal+0.5_DP)*PI)**2.0_DP) / mu
  GamnA = A + Gamman
  ExpGamt = EXP(Gamman*t)
  ExpGamAt = EXP((Gamman * (t - t1)) - (A * t1))
  Term1 = (ExpGamt  - ExpGamAt) / GamnA
  Term3 = Term1 * ((L * L * (b + ChangeDrivers(2,3)) / (mu * K * D)) - dHdt1)
  Term1 = Term1 * A * H1(1,2)
  Term2 = A * dHdt1 * (ExpGamt - (((GamnA * t1) + 1.0_DP)* ExpGamAt)) / (GamnA * GamnA)
! Calculate the terms that involve sums over intervals between observations of the surface water level and recharge R.
  Term4 = 0.0_DP
  Term5 = 0.0_DP
  Term6 = 0.0_DP
  IF (t > ChangeDrivers(2,1)) THEN
    ObsInt: DO Counter = 2, J
      H1i = ChangeDrivers(Counter,2)
      ti = ChangeDrivers(Counter,1)
!     For the final time interval, avoid that the values of the time and H1 overshoot to their values for the next 
!     observation time. Since R is a rate assumed constant over a time interval, its value does not need modification. 
      IF (Counter == J) THEN
        tipl1 = t
        H1ipl1 = H1i + ((ChangeDrivers(Counter+1,2) - H1i) * (t - ti) / (ChangeDrivers(Counter+1,1) - ti))
      ELSE
        tipl1 = ChangeDrivers(Counter+1,1)
        H1ipl1 = ChangeDrivers(Counter+1,2)
      END IF
      Ripl1 = ChangeDrivers(Counter+1,3)
      ExpGamAti = EXP((Gamman * (t - ti)) - (A * ti))
      ExpGamAtipl1 = EXP((Gamman * (t - tipl1)) - (A * tipl1))
      Help1 = (H1ipl1 - H1i) / (tipl1 - ti)
      Term4 = Term4 + ((H1i - (ti * Help1)) * (ExpGamAti - ExpGamAtipl1))
      Term5 = Term5 + (Help1 * ((((GamnA * ti) + 1.0_DP) * ExpGamAti) - (((GamnA * tipl1) + 1.0_DP) * ExpGamAtipl1)))
      Term6 = Term6 + (((L * L * (b + Ripl1) / (mu * K * D)) - Help1) * (ExpGamAti - ExpGamAtipl1))
    END DO ObsInt
    Term4 = A * Term4 / GamnA
    Term5 = A * Term5 / (GamnA * GamnA)
    Term6 = Term6 / GamnA
  END IF
  BCInt = Term1 + Term2 + Term3 + Term4 + Term5 + Term6
! The following IF-loop may be useful if diificulties arise, but generates considerable output.
! IF ((n == 0) .OR. (n == 2) .OR. (n == Maxn)) THEN
!   WRITE (30, 640) I, n, BCInt(I,n), Term1, Term2, Term3, Term4, Term5, Term6
! END IF
! 640 FORMAT (/,' Integral involving BC and external forcings. Output time nr =',I2,', n =',I4,', BCint(n)=',F10.4,/,&
!             ' Term1 =',F10.4,' Term2 =',F10.4,' Term3 =',F10.4,' Term4 =',F10.4,' Term5 =',F10.4,' Term6 =',F10.4)
END IF tzero
END SUBROUTINE BCIntegral
END PROGRAM DupuitFlow



