!23456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012
! Whollstoptherain simulates rainfall according to one of several Bartlett-Lewis-type models.
 MODULE DoublePrecision
!  This module defines the precision of Kind DP. 
   IMPLICIT NONE
   PUBLIC
   INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(20, 300)  
 END MODULE DoublePrecision
! 
 MODULE Parameters
!  This module stores the input parameters. It also gives the start times (at midnight) of all 48 months in a four-year period,
!  the last year of which is a leap year. Day 1155 is February 29. Note that the last and first entry both refer to midnight of
!  December 31st.
   USE DoublePrecision
   IMPLICIT NONE
   CHARACTER (LEN = 120) :: ProblemLabel
   CHARACTER (LEN = 11) :: ModelType, Months_Yes_or_No
   CHARACTER (LEN = 3) :: Output, FullRecord, HourlyRain
   INTEGER :: NrOfFourYearPeriods, NrOfPeriods, NrOfYears, Seed1                  
   REAL (KIND = DP), Dimension (:), Allocatable :: alpha, kappa, lambda, lower_eta, murate, mux, nu, phi, StartTimes
   REAL (KIND = DP), Dimension (1:49), Parameter :: MonthStartTimes = &
   (/ 0.0_DP,31.0_DP,59.0_DP,90.0_DP,120.0_DP,151.0_DP,181.0_DP,212.0_DP,243.0_DP,273.0_DP,304.0_DP,334.0_DP,365.0_DP,396.0_DP,&
   424.0_DP,455.0_DP,485.0_DP,516.0_DP,546.0_DP,577.0_DP,608.0_DP,638.0_DP,669.0_DP,699.0_DP,730.0_DP,761.0_DP,789.0_DP,820.0_DP,&
   850.0_DP,881.0_DP,911.0_DP,942.0_DP,973.0_DP,1003.0_DP,1034.0_DP,1064.0_DP,1095.0_DP,1126.0_DP,1155.0_DP,1186.0_DP,1216.0_DP,&
   1247.0_DP,1277.0_DP,1308.0_DP,1339.0_DP,1369.0_DP,1400.0_DP,1430.0_DP,1461.0_DP /) 
 END MODULE Parameters
!
 MODULE StormParameters
!  This module stores the parameters of individual storms, as well as the starting times of each period in the rainfall record
!  for which input parameters are asigned (AllStartTimes).
   USE DoublePrecision
   IMPLICIT NONE
   INTEGER :: NrOfStorms
   REAL (KIND = DP), Dimension (:), Allocatable :: AllStartTimes, StormStartTimesOld
 END MODULE StormParameters
!
 MODULE CellParameters
! This module stores variables that give the start time, end time, and rainfall rate of all cells, as well as the total amount 
! of rain calculated from that information.
   USE DoublePrecision
   IMPLICIT NONE
   INTEGER, Dimension (:), Allocatable :: NrOfCells
   REAL (KIND = DP) :: TotalAmountOfRain
   REAL (KIND = DP), Dimension (:,:), Allocatable :: AllCellsStartEndDepthOld
 END MODULE CellParameters
!
 MODULE RainfallTimeline
! This module stores the detailed timeline of generated rainfall: the number of intervals of varying lengths into which the 
! rainfall record is divded, and the constant rainfall rate of each of these intervals. It also stores the daily, hourly,
! monthly, and annual sums of rainfall. 
   USE DoublePrecision
   IMPLICIT NONE
   INTEGER :: NrOfIntervals
   REAL (KIND = DP), Dimension (:,:), Allocatable :: RainfallRate, DailySum, HourlySum, MonthlySum, AnnualSum
 END MODULE RainfallTimeline
!
 PROGRAM Whollstoptherain
! 
! This program simulates rainfall according to one of several Bartlett-Lewis-type models described in Pham et al. (2013)(An
! assessment of the ability of Bartlett-Lewis type of rainfall models to reproduce drought statistics, HESS 17:5167-5183,
! doi: 10.5194/hess-17-5176-2013.Of these, the most basic is left out. The others are the model described by Rodriguez-Iturbe et al.
! 1988 (A point process model for rainfall: further developments, Proc. R. Soc. Lond. A 417:283-298) or modifications thereof. 
! One of these modifications is Onof et al.'s (2013) suggestion to truncate the gamma-distribution of eta (see below). Onof et al.
! (2013) is unpublished.
! Another is to replace the exponential distribution of the rainfall rate in a rainfall cell by a gamma distribution (Onof and
! Wheater, 1994. Improvements to the modelling of British rainfall using a modified random parameter Bartlett-Lewis rectangular
! pulse model. J. Hydrol. 157:177-195.). 
! 
! This leads to four model configurations (abbreviations according to Pham et al., 2013).
! 1) MBL: Bartlett-Lewis with the parameters characterizing a storm (pdf of storm duration and pdf of cell arrival times) related to
! a gamma-distributed parameter eta. The rainfall rate in a rain cell has an exponential distribution. This is the Modified
! Bartlett-Lewis model, introduced by Rodriguez-Iturbe et al. (1988).
! 2) TBL: As MBL, but with a truncated gamma distribution for eta (Truncated Bartlett-Lewis).
! 3) MBLG: As MBL, but with the rainfall rate per cell having a gamma distribution.
! 4) TBLG: As TBL, but with the rainfall rate per cell having a gamma distribution.
! For very long rainfall series, using the untruncated gamma distribution can lead to physically unrealistic rainfall amounts for
! small values of eta because the moments of the gamma distribution are unbounded for values of alpha < 4 (Onof et al., 2103).
! Excessively small values of eta can be avoided by simply not permitting values below a threshold epsilon, which is the rationale
! for TBL and TBLG. Pham et al. (2013) found that MBL outperformed the other models, but still TBL might be a good choice.
! For that reason, MBL and MBLG are not explicitly included in the input. To run an MBL- or MBLG-type model, the parameter
! lower_eta,  which gives the lower threshold for eta, needs to be set to zero in the input file for each period for which the
! parameters are given.
!
! The code requires the model parameters for arbitrary periods of time in which a year is subdivided (usually for each month).
! Storms are generated for a given period until the starting time of a storm is beyond that period. This starting time is not
! discarded, but it serves as the starting time of the first storm in the period in which it occurs. The parameters of that
! storm are generated using the pdfs for the period in which it occurs, not the period in which it was generated. Subsequent
! storms are generated using the parameters of that period. 
! Rain cells in a storm are generated until a rain cell has a starting time that exceeds the time at which the storm ends. 
! This latest cell is discarded. Rain cells may overlap, creating irregular rainfall rates over time.
!
! The resulting rainfall from individual cells is used to create a time line of rainfall rates (including zero, for dry spells) that
! gives a continuous record of the rainfall rate as a function of time. This record is processed to provide hourly, daily, monthly,
! and annual sums. All three records are written to output files.
! 
! Interface: Input file RainPar.IN contains a problem label (ProblemLabel), a label (ModelType) that indicates whether the rainfall
! rate of a cell is gamma ('Gam') or exponentially ('Exp') distributed, and a label (Months_Yes_or_No) to indicated whether the time
! periods for which the parameters are given are months ('Yes' or 'No'). If 'No', the number of periods for which parameters are
! given (>= 1) needs to be specified (NrOfPeriods).
! Label 'Output' indicates whether or not additional information about the structure of rain storms should be written to file
! Diagnostics.OUT.
! Irrespective of the periods being months or not, the number of years for which rainfall is to be generated needs to be given
! (NrOfYears). This number is rounded up to a multiple of four. The simulated rainfall record starts on January 1st 00:00 hrs, and 
! ends on December 31, 24:00 hrs. In the simulated period, every fourth year is a leap year.
! The model assumes that the time unit is a day, and that the specified periods (if they are not months) together constitute a 
! non-leap calendar year. 
! The following parameters are required for each time period: 
! lambda [T^-1] (governs the exponentially distributed intervals between storm arrival times)
! mux ([LT^-1] when the rainfall rate per cell is exponentially distributed, dimensionless if it is Gamma-distributed - mux then
!    is the shape parameter)(governs cell depth, i.e. rainfall rate)
! murate [L^-1 T] (this is the rate parameter of the gamma distribution of the rainfall rate)
! alpha (dimensionless shape parameter of the gamma distribution of eta, with eta [T^-1] governing the exponential distribution
!    of cell durations within a storm)
! nu [T^-1] (rate parameter of that gamma distribution)
! kappa (ratio of beta [T^-1] over eta)
! phi (ratio of gamma [T^-1] over eta)
! lower_eta [T^-1](the minimum permitted value of eta). 
! 
! N.B. Eta in Pham et al. (2013) has dimension T-1, and its rate parameter nu has dimension T. Onof and Wheater (1994) gave the
! scale parameter delta the opposite units of the rainfall rate, it is therefore a rate parameter, just as nu. In the code, their
! shape parameter p is represented by mux. Their rate parameter delta is represented by murate. Variable mux holds the reciprocal
! of the parameter of the exponential distribution of the rainfall rate if it is not Gamma-distributed, which is the mean rainfall
! rate [LT^-1]. Inconsistent with this, but consistent with the notation in the cited literature, the rate parameter of the
! rainfall rate has dimensions TL^-1. The reciprocals of nu and murate must be used by the routine that generates Gamma variates.
! 
! When lower_eta is set to zero, the model simplifies to the MBL or MBLG.  Beta (not required on input) governs the exponential
! distribution of intervals between arrival times of rain cells within storms, and gamma [T^-1] (not required on input) the 
! exponentially distributed storm duration.  For a full description see Rodriguez-Iturbe et al. (1988) and Pham et al. (2013).
! The latter also gives monthly parameter values for the Uccle (Ukkel) rainfall gauge in Belgium. Only if the time periods do not
! reflect months should the input file contain the starting times (in days since the start of the year) of the periods for which
! the parameters are valid (StartTimes). This implies the first period must start at time 0.0. The starting times should be valid
! for non-leap years. For leap years, the program adds a day to the period with February 29th in it. The final period of which the
! starting time is given on input is assumed to end at the end of a non-leap year (365.0 days).
!
! Finally, RainPar.IN contains an arbitrary integer number that helps initialize the random number generator (Seed1).
!
! The input file has the following format
! Line 1: Problem label (maximum 120 characters)
! Line 2: will not be read, can be used for a label.
! Line 3: Months_Yes_or_No. Fill in YES, Yes, yes, Y, or y if the model parameters are given for months. Other input will be
!         treated as no.
! Line 4: will not be read, can be used for a label.
! Line 5: ModelType. Fill in GAMMA, Gamma, gamma, GAM, Gam, or gam if the rainfall rate of a cell is gamma distributed. Anything
!         else will result in an exponential distribution rainfall rates for individual cells.
! Line 6: will not be read, can be used for a label
! Line 7: Output. Fill in YES, YE, or Y (case insensitive) if additional informative output about storms and cells is desired. 
!         Other input will be treated as no.
! Line 8: will not be read, can be used for a label.
! Line 9: FullRecord. Fill in YES, YE, or Y (case insensitive) if a file with the continuous time line of time intervals and
!         their rainfallrate is to be generated. Other input will be treated as no.
! Line 10:will not be read, can be used for a label.
! Line 11:HourlyRain. Fill in YES, YE, or Y (case insensitive) if a file with hourly sums of rainfall is to be generated. 
!         Other input will be treated as no.
! Line 12:will not be read, can be used for labels.
! Line 13:Seed1, (NrOfPeriods, if Months_Yes_or_No equals 'Yes' or one of its equivalents), NrOfYears. 'Seed1' is a large, arbitrary
!         integer that initializes the random number generator. 'NrOfPeriods' gives the number of periods in a year for which the
!         rainfall parameters are available. This number defaults to twelve when Months_Yes_or_No equals 'Yes' or one of its
!         equivalents). NrOfYears gives the number of years that need to be covered by the rainfall record. The rainfall record
!         starts at January 1st of a year and ends at December 31st of another year. The code rounds up the number of years to the
!         next multiple of four to facilitate the inclusion of leap year. In the generated rainfall record every fourth year is a
!         leap year.
! Line 14:will not be read, can be used for labels.
! Line 15 through line (14 + NrOfPeriods) give the parameters for period I (I runs from 1 to NrOfPeriods) in the following order:
!        lambda(I)  mux(I)  (murate(I), if ModelType = Gam or any of its equivalents)  alpha(I)  nu(I)  kappa(I)  phi(I) 
!        epsilon(I) (StartTimes(I), if Months_Yes_or_No equals 'Yes' or one of its equivalents)
!        murate(I) is needed if the rainfall rate per cell is gamma distributed. In that case, mux(I) gives the shape parameter
!        of the gamma distribution, while murate gives the rate parameter. StartTimes(I) only needs to be included if the periods
!        do not reflect months. The first entry needs to be 0.0, and the times must be in ascending order. The last entry
!        (StartTimes(NrOfPeriods) should not exceed 365.0 The start times are valid for a non-leap year. The code will increase the
!        period in which February 29 falls by one day every fourth year.
! 
! N.B. The code requires that the unit of time is one day. If the parameters were determined for other time units they need to be
! converted. In the table below, the ratio day/unit is the number that results if the origin unit (unit) and the target unit (day)
! are both expressed in terms of the same unit. If the original unit is 1 hour, day/unit becomes 24.
!
! Parameter      Dimension    Value in terms of
!                             ___________________________________________________________________
!                             original unit  days            days in case original unit is 1 hour
! _______________________________________________________________________________________________
! lambda(I)      1/T          x              x * day/unit    24x
! mux(I)         L/T          x              x * day/unit    24x (ModelType /= Gamma)
! mux(I)          -           x              x               x (ModelType = Gamma)
! murate(I)      T/L          x              x * unit/day    x/24
! alpha(I)        -           x              x               x
! nu(I)           T           x              x * unit/day    x/24
! kappa(I)        -           x              x               x
! phi(I)          -           x              x               x
! epsilon(I)     1/T          x              x * day/unit    24x
! StartTimes(I)   T           x              x * unit/day    x/24
!
! Output: 
! File Diagnostics.OUT contains an echo of the input as well as information about the way the input data was processed. It also
! contains some statistics about the generated rainfall record.
! File Rainfall_Rate.OUT contains the times at which the rainfall rate changes and gives the rainfall rate for the period between 
! the time in the same record and the time in the previous record (or zero for the first record). The rainfall rates appear in the
! first column, the times separating the time intervals are in the second column.
! File Hourly_Rainfall.OUT contains the hourly sums for the rainfall of the entire simulation period. The sum in a record is valid
! for the hour preceding the time in that record. This file is only produced if HourlyRain (in Module Parameters) equals 'Yes'.
! File Daily_Rainfall.OUT contains the daily sums for the days ending at the times for which the sums are given.
! For both files, the first column gives the time (the unit is day), and the second column the total rainfall in the hour or
! day preceding the time in the same row.
! Files Monthly_Rainfall.OUT and Annual_Rainfall.OUT give the monthly (annual) sum in the second column for the month (year)
! number in the first column.
! The hourly, daily, monthly, and annual sums of rainfall are stored in arrays HourlySum, DailySum, MonthlySum, and AnnualSum in
! Module RainfallTimeline. 
!
! Gerrit H. de Rooij, October 2014. Modifications, February 2018
!
 USE DoublePrecision
 USE Parameters
 USE StormParameters, ONLY: AllStartTimes
 USE CellParameters, ONLY: NrOfCells
! Variable declarations
 IMPLICIT NONE
! INTEGER ClusterEndStorm, ClusterStartStorm, DayCount, EndCell, EndingCell K &
!       & NrStart, NrStartOld, OldestActiveCell, StartCell
! REAL (KIND = DP) :: ClusterStartTime, HourlySum&
!                   & NextCellStart, PeriodDuration, SoonestEnd, StormInterval                   
! REAL (KIND = DP), Dimension (:), Allocatable :: 
!                                               &  RainPerInterval, RelativeError,  &
!                                               &    
! REAL (KIND = DP), Dimension (:,:), Allocatable ::  DailySum, RainfallRate, RainfallRateOld, SortedArray, &
!                                                & StartEnd, UnsortedSegment                                                 
!
 INTEGER :: I, Help
 INTEGER, Dimension (:), Allocatable :: Seed
!
! Read input and open the file for diagnostic information.
 OPEN (22, FILE = 'Diagnostics.OUT', STATUS = 'UNKNOWN')
 CALL InputReader
! Initialize random number generator.
 CALL Random_seed (SIZE = Help)
 Allocate (Seed(1:Help))
 Seed(1) = Seed1
 DO I = 2, Help, 1
   Seed(I) = Seed(1) + (I * 97) ! 97 is the largest prime number smaller than 100
 END DO
 CALL Random_seed (PUT = Seed)
 Deallocate (Seed)
! Generate the starting times of the rainstorms.
 CALL StormStarts
! Generate the rainfall within each storm
 CALL RainMaker
! Sort out the timeline of the rain cells and generate the detailed rainfall record, with a complete timeline of consecutive,
! non-equidistant, time intervals and their rainfall rates.
 CALL Timeline
! Generate the hourly, daily, monthly, and annual rainfall sums and write the full and derived rainfall records to output files.
 CALL OutputWriter
!
 CLOSE (22)
!
 Deallocate (lambda)
 Deallocate (mux)
 Deallocate (alpha)
 Deallocate (nu)
 Deallocate (kappa)
 Deallocate (phi)
 Deallocate (lower_eta)
 Deallocate (NrOfCells)
 Deallocate (StartTimes)
 Deallocate (AllStartTimes)
 IF (ModelType == 'Gam') Deallocate (murate)
!
 CONTAINS
!
 SUBROUTINE DrawGammaValue (Alpha, Nu, MinimumValue, Variate)
! DrawGammaValue generates a variate from a gamma distribution with shape parameter Alpha and scale parameter Nu. The random value
! is only accepted if it is no smaller than MinimumValue.
! The algorithm is that of Xi, B., K. Ming Tan, and C. Li. 2013. Logarithmic transformation-based gamma random number generators.
! Journal of Statistical Software 55(4):1-17. http://www.jstatsoft.org/ 
!
! Interface: shape parameter Alpha, scale parameter Nu, and the smallest permitted value of the variate are obtained on input, the
! variate is passed on to the main program through variable Variate. The scale parameter and the truncation parameter MinimumValue
! have units equal to that of the reciprocal of the time coordinate in case the routine is called to calculate the paramter
! governing the storm duration (1/day in this code). When the rainfall rate is to be generated, the dimension of these parameters
! is day/length (reciprocal of the equivalent water layer per time). The shape parameter is dimensionless.
! The subroutine calls Random-Number twice for variates from the standard uniform distribution. The precision of the REAL variables
! for which a high precision is required is provided by MODULE Precision.
!
! Gerrit H. de Rooij, October 2014, modified January 2018
!
 USE DoublePrecision
 IMPLICIT NONE
 REAL (KIND = DP), Intent(in) :: Alpha, Nu, MinimumValue
 REAL (KIND = DP), Intent(out) :: Variate
 REAL (KIND = DP) :: LnAlpha, SqrRtAlpha, Bs, Bw, Bmin, Bmax, Umax, Vmin, Vmax, U, V, T, T1
! This routine generates a variate from the unit gamma distribution according to algorithm 1 or 2 in Xi, B., K. Ming Tan, and
! C. Li. 2013. This variate is stored in T1 and then transformed to the desired gamma distribution by means of the scale parameter
! Nu. When it exceeds the minimum permitted value it is passed on to the main program.
! Algorithm 1 works best for Alpha > 1, algorithm 2 only works when Alpha <= 1.
!
 AlphaBigOrSmall: IF (Alpha > 1.0_DP) THEN
!  Use algorithm 1
   LnAlpha = LOG(Alpha)
   SqrRtAlpha = SQRT(Alpha)
!  The calculations below define the envelope within which the permissible sampling region of the variate lies.
   IF (LnAlpha < 0.20931492_DP) THEN
     Bw = LOG(2.0_DP / EXP(1.0_DP)) + ((Alpha - LnAlpha) / 2.0_DP)
   ELSE IF (LnAlpha < 0.52122324_DP) THEN
     Bw = (-0.13546023_DP * LnAlpha) + 0.12789964_DP
   ELSE IF (LnAlpha < 1.76421669_DP) THEN
     Bw = (-0.08476353_DP * LnAlpha) + 0.10147534_DP
   ELSE
     Bw = -0.04806589_DP
   END IF
   IF (LnAlpha <= -3.33318991_DP) THEN
     Bs = (0.30625300_DP * LnAlpha) + 0.27127295_DP
   ELSE IF (LnAlpha <= 1.44893155_DP) THEN
     Bs = (0.12465180_DP * LnAlpha) - 0.33403833_DP
   ELSE
     Bs = -0.15342641_DP
   END IF
   Bmax = EXP(Bs)
   Bmin = -EXP(Bw)
!
!  Take two random numbers from the standard uniform distribution and transform to a variate that has a gamma distribution if it  
!  falls within the permitted region. Repeat the process until a variate is found.
   DO
     CALL Random_Number(U)
     CALL Random_Number(V)
     V = (V * (Bmax - Bmin)) + Bmin
     T = V / U
     T1 = EXP((T / SqrRtAlpha) + LnAlpha)
!    Check if the variate is in the permitted region, and if it exceeds the smallest permitted value.
     IF (((2.0_DP * LOG(U)) <= (Alpha - T1 + (SqrRtAlpha * T))) .AND. ((T1 * Nu) >= MinimumValue)) THEN
!    Transform the variate to a gamma distribution with parameters Alpha and Nu, pass on to the main program, and exit the loop.
       Variate = T1 * Nu
       EXIT
     END IF
   END DO
 ELSE
!  Use algorithm 2.
   Umax = (Alpha * EXP(-1.0_DP)) ** (Alpha / 2.0_DP)
   Vmin = 2.0_DP * EXP(-1.0_DP)
   Vmax = Vmin * Alpha / (EXP(1.0_DP) - Alpha)
   Vmin = -1.0_DP * Vmin
!  Take two random numbers from the standard uniform distribution and transform to a variate that has a gamma distribution if it
!  falls within the permitted region. Repeat the process until a variate is found.
   DO
     CALL Random_Number(U)
     CALL Random_Number(V)
     U = Umax * U
     T = ((V * (Vmax - Vmin)) + Vmin) / U
     T1 = EXP(T / Alpha)
!    Check if the variate is in the permitted region, and if it exceeds the minimum permitted value.
     IF (((2.0_DP * LOG(U)) <= (T - T1)) .AND. ((T1 * Nu) >= MinimumValue)) THEN
!    Transform the variate to have a gamma distribution with parameters Alpha and Nu, pass on to the main program, and exit the loop.
       Variate = T1 * Nu
       EXIT
     END IF
   END DO
 END IF AlphaBigOrSmall
 END SUBROUTINE DrawGammaValue
! 
 SUBROUTINE InputReader
! InputReader reads the input and processes it for use by the rest of the program.
! 
! Interface: Input file RainPar.IN contains a problem label (ProblemLabel), a label (ModelType) that indicates whether the rainfall
! rate of a cell is gamma ('Gam') or exponentially ('Exp') distributed, and a label (Months_Yes_or_No) to indicated whether the time
! periods for which the parameters are given are months ('Yes' or 'No'). If 'No', the number of periods for which parameters are
! given (>= 1) needs to be specified (NrOfPeriods).
! Parameters 'Output', 'FullRecord', and 'HourlyRain' (all Yes/No) determine if informative output about storms and rain cells
! Needs to be written to file Diagnostics.OUT (Output), if the continuous timeline of time intervals and rainfall rates needs to
! be wriiten to output (FullRecord), and if record of hourly raifall sums needs to be writen to output (HourlyRain). For very long 
! simulation periods, this can give very large output files.
! Irrespective of the periods being months or not, the number of years for which rainfall is to be generated needs to be given
! (NrOfYears). This number is rounded up to a multiple of four. The simulated rainfall record starts on January 1st 00:00 hrs, and 
! ends on December 31, 24:00 hrs. In the simulated period, every fourth year is a leap year.
! The model assumes that the time unit is a day, and that the specified periods (if they are not months) together constitute a 
! non-leap calendar year. 
! The following parameters are required for each time period: 
! lambda [T^-1] (governs the exponentially distributed intervals between storm arrival times)
! mux ([LT^-1] when the rainfall rate per cell is exponentially distributed, dimensionless if it is Gamma-distributed - mux then
!    is the shape parameter)(governs cell depth, i.e. rainfall rate)
! murate [L^-1 T] (this is the rate parameter of the gamma distribution of the rainfall rate)
! alpha (dimensionless shape parameter of the gamma distribution of eta, with eta [T^-1] governing the exponential distribution
!    of cell durations within a storm)
! nu [T^-1] (rate parameter of that gamma distribution)
! kappa (ratio of beta [T^-1] over eta)
! phi (ratio of gamma [T^-1] over eta)
! lower_eta [T^-1](the minimum permitted value of eta). 
! 
! N.B. Eta in Pham et al. (2013) has dimension T-1, and its rate parameter nu has dimension T. Onof and Wheater (1994) gave the
! scale parameter delta the opposite units of the rainfall rate, it is therefore a rate parameter, just as nu. In the code, their
! shape parameter p is represented by mux. Their rate parameter delta is represented by murate. Variable mux holds the reciprocal
! of the parameter of the exponential distribution of the rainfall rate if it is not Gamma-distributed, which is the mean rainfall
! rate [LT^-1]. Inconsistent with this, but consistent with the notation in the cited literature, the rate parameter of the
! rainfall rate has dimensions TL^-1. The reciprocals of nu and murate must be used by the routine that generates Gamma variates.
! 
! When lower_eta is set to zero, the model simplifies to the MBL or MBLG.  Beta (not required on input) governs the exponential
! distribution of intervals between arrival times of rain cells within storms, and gamma [T^-1] (not required on input) the 
! exponentially distributed storm duration.  For a full description see Rodriguez-Iturbe et al. (1988) and Pham et al. (2013).
! Only if the time periods do not
! reflect months should the input file contain the starting times (in days since the start of the year) of the periods for which
! the parameters are valid (StartTimes). This implies the first period must start at time 0.0. The starting times should be valid
! for non-leap years. For leap years, the program adds a day to the period with February 29th in it. The final period of which the
! starting time is given on input is assumed to end at the end of a non-leap year (365.0 days).
!
! Finally, RainPar.IN contains an arbitrary integer number that helps initialize the random number generator (Seed1).
!
! The parameters and arrays arising from them are passed on to other program elements through Module Parameters. Information about
! the precision of the real varaibles with higher than single precision is obtained from Module DoublePrecision.
!
! Gerrit H. de Rooij, February 2018
!
 USE DoublePrecision
 USE Parameters
! Variable declarations
 INTEGER :: Help, I, J
! 
 OPEN (21, FILE = 'RainPar.IN', STATUS = 'OLD', ACTION = 'READ') 
 READ (21, '(A)') ProblemLabel
 READ (21, '(/,A)') Months_Yes_or_No
 READ (21, '(/,A)') Output
 IF ((Output == 'Y').OR.(Output == 'YE').OR.(Output == 'YES')) THEN
   Output = 'Yes'
 ELSE
   Output = 'No'
 END IF
 READ (21, '(/,A)') FullRecord
 IF ((FullRecord == 'Y').OR.(FullRecord == 'YE').OR.(FullRecord == 'YES')) THEN
   FullRecord = 'Yes'
 ELSE
   FullRecord = 'No'
 END IF
 READ (21, '(/,A)') HourlyRain
 IF ((HourlyRain == 'Y').OR.(HourlyRain == 'YE').OR.(HourlyRain == 'YES')) THEN
   HourlyRain = 'Yes'
 ELSE
   HourlyRain = 'No'
 END IF
 READ (21, '(/,A)') ModelType
 IF ((ModelType == 'G').OR.(ModelType == 'GA').OR.(ModelType == 'GAM').OR.(ModelType == 'GAMM').OR.(ModelType == 'GAMMA')) THEN
   ModelType = 'Gam'
 ELSE
   ModelType = 'Exp'
 END IF
 IF ((Months_Yes_or_No == 'YES') .OR. (Months_Yes_or_No == 'YE') .OR. (Months_Yes_or_No == 'Y')) THEN
! Prepare for leap years: each four-year period has 48 months.
   Months_Yes_or_No = 'Yes'
! First read the seed for the random generator and the number of years for which rainfall needs to be generated, then allocate
! the arrays of the input variables.
   READ (21, *)
   READ (21, *) Seed1, NrOfYears
   NrOfPeriods = 12 
   Allocate (lambda (1:48))
   Allocate (mux (1:48))
   Allocate (alpha (1:48))
   Allocate (nu (1:48))
   Allocate (kappa (1:48))
   Allocate (phi (1:48))
   Allocate (lower_eta (1:48))
   Allocate (StartTimes (1:49))
   IF (ModelType == 'Gam') Allocate (murate (1:48))
 ELSE
   Months_Yes_or_No = 'No'
! First read the seed for the random generator, the number of periods in a year for which the input variables are valid, 
! and the number of years for which rainfall needs to be generated, then allocate the arrays of the input variables.
   READ (21, *) 
   READ (21, *) Seed1, NrOfPeriods, NrOfYears
   Allocate (lambda (1:(4 * NrOfPeriods)))
   Allocate (mux (1:(4 * NrOfPeriods)))
   Allocate (alpha (1:(4 * NrOfPeriods)))
   Allocate (nu (1:(4 * NrOfPeriods)))
   Allocate (kappa (1:(4 * NrOfPeriods)))
   Allocate (phi (1:(4 * NrOfPeriods)))
   Allocate (lower_eta (1:(4 * NrOfPeriods)))
   Allocate (StartTimes (1:(4 * NrOfPeriods) + 1))
   IF (ModelType == 'Gam') Allocate (murate (1:(4 * NrOfPeriods)))
 END IF
 READ (21, *)
! Read the input parameters for all periods. If these periods are not months, this includes the starting day of each period 
! in a non-leap year.
 DO I = 1, NrOfPeriods
   IF (Months_Yes_or_No == 'Yes') THEN
     IF (ModelType == 'Gam') THEN
       READ (21, *) lambda(I), mux(I), murate(I), alpha(I), nu(I), kappa(I), phi(I), lower_eta(I)
     ELSE
       READ (21, *) lambda(I), mux(I), alpha(I), nu(I), kappa(I), phi(I), lower_eta(I)
     END IF
!    The following DO loop assigns the input variables of January to all Januaries in the four-year period, and similarly 
!    for the other months.
     DO J = 12, 36, 12
       lambda(I + J) = lambda(I)
       mux(I + J) = mux(I)
       alpha(I + J) = alpha(I)
       nu(I + J) = nu(I)
       kappa(I + J) = kappa(I)
       phi(I + J) = phi(I)
       lower_eta(I + J) = lower_eta(I)
       IF (ModelType == 'Gam') murate(I + J) = murate(I)
     END DO
   ELSE
     IF (ModelType == 'Gam') THEN
       READ (21, *) lambda(I), mux(I), murate(I), alpha(I), nu(I), kappa(I), phi(I), lower_eta(I), StartTimes(I)
     ELSE
       READ (21, *) lambda(I), mux(I), alpha(I), nu(I), kappa(I), phi(I), lower_eta(I), StartTimes(I)
     END IF
   END IF
 END DO
 CLOSE (21)
!
 StartOfPeriods: IF (Months_Yes_or_No == 'Yes') THEN
! Below are the starting times (at midnight) of 48 months in a four-year cycle, the final one of which is a leap year. 
! The final time is the starting time of the next cycle. 
   StartTimes = MonthStartTimes
 ELSE
!  The code in this part of the IF construct finds the start and end of each period of time within a four-year cycle for which the
!  input parameters are given.
!  It includes a leap year as the final year. It is assumed that the time periods given on input together cover a single
!  non-leap year. 
!  Ensure that the first period starts at time zero.
   StartTimes(1) = 0.0_DP
!  Find the period that includes February 28.
   IF (StartTimes(2) >= 59.0_DP) THEN
!    February 28 is in period 1. Make sure Feb 29 ends up in period 1 of the fourth year.
     I = 1 + (3 * NrOfPeriods)
   ELSE
     IF (NrOfPeriods > 1) THEN
!      StartTimes(1) = 0 by definition, StartTimes(2) is the end of the the first period and was already covered above, so we 
!      can move on to the start of the second period, at StartTimes(3).
       DO J = 3, NrOfPeriods + 1
         IF (StartTimes(J) >= 59.0_DP) THEN
!          February 28 is in Period J - 1. Make sure Feb 29 ends up in that period in the fourth year, and then leave the loop.
           I = J - 1 + (3 * NrOfPeriods)
           EXIT
         END IF
       END DO
     END IF
   END IF
!  Calculate the starting times of all periods within the first four years from those of year 1. Start with year two and three.
   DO J = NrOfPeriods + 1, (2 * NrOfPeriods) 
     StartTimes(J) = StartTimes(J - NrOfPeriods) + 365.0_DP
     StartTimes(J + NrOfPeriods) = StartTimes(J - NrOfPeriods) + 730.0_DP
   END DO
!  Year 4 up to and including the period with February 29 in it (given by I).
   DO J = (3 * NrOfPeriods + 1), I 
     StartTimes(J) = StartTimes(J - (3 * NrOfPeriods)) + 1095.0_DP
   END DO
!  The rest of year 4.
   DO J = I + 1, (4 * NrOfPeriods)
     StartTimes(J) = StartTimes(J - (3 * NrOfPeriods)) + 1096.0_DP
   END DO
!  The end of the final day of the 4-year period.
   StartTimes((4 * NrOfPeriods) + 1) = 1461.0_DP
!  Now assign the correct parameter values to all periods within the four-year cycle for user-specified periods.
   DO I = NrOfPeriods, 3 * NrOfPeriods, NrOfPeriods
     DO J = 1, NrOfPeriods
       lambda(I + J) = lambda(J)
       mux(I + J) = mux(J)
       alpha(I + J) = alpha(J)
       nu(I + J) = nu(J)
       kappa(I + J) = kappa(J)
       phi(I + J) = phi(J)
       lower_eta(I + J) = lower_eta(J)
       IF (ModelType == 'Gam') murate(I + J) = murate(J)
     END DO
   END DO
 END IF StartOfPeriods
!
! Round up the years for which rainfall has to be generated to multiples of four years. Possible values of Help are 0, 
! 1, 2, and 3.
 Help = MOD(NrOfYears, 4)
 IF (Help == 0) HELP = 4
 NrOfYears = NrOfYears + 4 - Help
 NrOfFourYearPeriods = NrOfYears / 4
! 
! Echo (modified) input to the output file, including a sample of the input arrays for verification purposes.
 IF (ModelType == 'Gam') THEN
   WRITE(22, 100) ProblemLabel, Months_Yes_or_No, ModelType, NrOfPeriods, NrOfYears, StartTimes(2), lambda(2), mux(2),&
                & murate(2), alpha(2), nu(2), kappa(2), phi(2), lower_eta(2)
 ELSE
   WRITE(22, 110) ProblemLabel, Months_Yes_or_No, ModelType, NrOfPeriods, NrOfYears, StartTimes(2), lambda(2), mux(2),&
                & alpha(2), nu(2), kappa(2), phi(2), lower_eta(2)
 END IF
100 FORMAT (1X, A,/,' Are periods months? ', A,' Gamma or exp distr.? ', A,/,&
           &' Nr of periods in a year  Nr of years',/, 1X, I5, 20X, I5, /,&
           &' Period 2: Start time      lambda     mu-shape      mu-rate        alpha           nu        kappa          phi   &
           &minimum eta',/,8X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1x, ES12.4, 1X, ES12.4, 1X,&
           & ES12.4)
110 FORMAT (1X, A,' Are periods months? ', A,' Gamma or exp distr.? ', A,/,&
           &' Nr of periods in a year  Nr of years',/, 1X, I5, 20X, I5, /,&
           &' Period 2: Start time      lambda     mu-sub-x        alpha           nu        kappa          phi   minimum eta',/,&
           &8X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1X, ES12.4, 1x, ES12.4, 1X, ES12.4)
 END SUBROUTINE InputReader
!
 SUBROUTINE OutputWriter
! OutputWriter generates hourly, daily, monthly, and annual rainfall sums and writes the full rainfall record and these
! derived records of sums to various output files.
!
! Interface: The full rainfall record is taken from RainfallRate in Module RainfallTimeline. The full time line in this array is
! written to file RainfallRate.OUT. In this file, the first column gives series of rainfall rates, the second a series of times.
! For row x the first column applies to the time interval between the times in row x-1 and x. The rainfall
! rate in he first row applies to the time interval between zero and the time in the first row.
! If HourlyRain (in Module Parameters) equals 'Yes', the hourly rainfall sums are written to file Hourly_Rainfall.OUT. The daily
! rainfall sums are written to file Daily_Rainfall.OUT. Similarly, monthly and annual sums are written to file
! Monthly_Rainfall.OUT and Annual_Rainfall.OUT, respectively.
! The hourly, daily, monthly, and annual sums of rainfall are stored in arrays HourlySum, DailySum, MonthlySum, and AnnualSum, 
! respectively, in Module RainfallTimeline. 
! Descriptive statistics of the rainfall record are written to file Rain_Statistics.OUT.
! Information about the required precision of double precision real variables is obtained from Module DoublePrecision.
!
! Gerrit H. de Rooij, March 2018
!
 USE DoublePrecision
 USE Parameters, ONLY: HourlyRain, FullRecord, NrOfFourYearPeriods, NrOfYears, MonthStartTimes
 USE RainfallTimeline
! Variable declarations
 IMPLICIT NONE
 INTEGER :: CyclicMonth, DayCount, EndDay, FourYearPeriodNr, Help, Help2, I, J, StartDay
 REAL (KIND = DP) :: AnnualMean, AnnualSD, Time
 REAL (KIND = DP), Dimension(1:12) :: MonthlyMean, MonthlySD
!
! Create an output file with times and rainfall rates since the previous time for all these times. If the unsaturated zone
! is considered passive, this file can serve as input for the groundwater model DupuitFlow. In that case it needs to be
! renamed R.IN. Note that the first column (Zero rainfall at time zero) from array RainfallRate does not need to be written
! to the output file. 
 IF (FullRecord == 'Yes') THEN
 OPEN(23, FILE = 'Rainfall_Rate.OUT', STATUS = 'UNKNOWN')
   DO I = 1, NrOfIntervals
     WRITE (23, 500) RainfallRate(I, 2), RainfallRate(I, 1)
   END DO
   CLOSE (23)
 END IF
500 FORMAT (1X, ES12.4, 1X, F15.6)
!
! Also create files with hourly, daily, monthly, and annual rainfall sums.
 OPEN(24, FILE = 'Daily_Rainfall.OUT', STATUS = 'UNKNOWN')
 IF (HourlyRain == 'Yes') OPEN(25, FILE = 'Hourly_Rainfall.OUT', STATUS = 'UNKNOWN')
 OPEN(26, FILE = 'Monthly_Rainfall.OUT', STATUS = 'UNKNOWN')
 OPEN(27, FILE = 'Annual_Rainfall.OUT', STATUS = 'UNKNOWN')
! Run through all hours of the simulated period (35064 hours in each four-year period). Figure out which time interval(s) are 
! represented in each hour, and calculate the hourly rainfall sum from that. Sum these in 24-hour blocks to find daily sums of
! rainfall.
! Help stores the number of the interval with constant rainfall (including zero) stored in RainfallRate that was active at the end
! of the previous hour. Help2 serves as a counter for the hours. Every time it reaches 24, a daily rainfall amount can be written
! to output using DailySum. DailySum also stores every daily rainfall amount so it is available for processing by the unsaturated
! zone module that uses daily input to compute daily recharge rates.
 Help = 1
 Help2 = 0
 Allocate (DailySum(1: 1461 * NrOfFourYearPeriods, 1:2))
 Allocate (HourlySum(1: 35064 * NrOfFourYearPeriods, 1:2))
 Allocate (MonthlySum(1: 12 * NrOfYears, 1:2))
 Allocate (AnnualSum(1: NrOfYears, 1:2))
! DailySum needs to initialized.
 DailySum = 0.0_DP
 DayCount = 0
 EveryHour: DO I = 1, 35064 * NrOfFourYearPeriods
!  Calculate the time in days since the start of the simulation period.
   Time = REAL(I, DP) / 24.0_DP
   Help2 = Help2 + 1
!  Find the active episode(s) in the hour ending at Time. The number of intervals in RainfallRate is given by NrOfIntervals.
   IF (RainfallRate(Help, 1) > Time) THEN
!    There is only one active time interval in this hour. It continues in the next, so Help should not be updated.
     HourlySum(I, 2) = RainfallRate(Help, 2) / 24.0_DP
   ELSE
!    Find the all intervals within this hour and compute their contributions to the total rainfall of this hour. 
!    Of the first interval, only the section that is contained within this hour contributes to the hourly sum.
     HourlySum(I, 2) = RainfallRate(Help, 2) * (RainfallRate(Help, 1) - Time + (1.0_DP/ 24.0_DP)) 
     DO J = Help + 1, NrOfIntervals
       IF (RainfallRate(J, 1) <= Time) THEN
         HourlySum(I, 2) = HourlySum(I, 2) + (RainfallRate(J, 2) * (RainfallRate(J, 1) - RainfallRate(J-1, 1)))
       ELSE
!        For the last interval, only the section within the hour contributes to the sum. Help is updated so the next hour starts
!        with the contribution of this interval.
         HourlySum(I, 2) = HourlySum(I, 2) + (RainfallRate(J, 2) * (Time - RainfallRate(J-1, 1)))
         Help = J
         EXIT
       END IF
     END DO 
   END IF
!  Store the time in array HourlySum.
   HourlySum(I, 1) = Time
   IF (HourlyRain == 'Yes') WRITE (25, 510) HourlySum(I, 1), HourlySum(I, 2)
510  FORMAT (1X, F16.6, 1X, ES12.4)
!  Check if 24 hours have passed. If so, store the daily sum in its array and write it to output and reset the counter for the
!  hours (Help2) and the summed rainfall (DailySum).
   DailySum(DayCount + 1, 2) = DailySum(DayCount + 1, 2) + HourlySum(I, 2)
   IF (Help2 == 24) THEN
     DayCount = DayCount + 1
     DailySum(DayCount, 1) = Time
     WRITE(24, 520) Time, DailySum(DayCount, 2)
     Help2 = 0
   END IF
 END DO EveryHour
 CLOSE (24)
 IF (HourlyRain == 'Yes') CLOSE (25)
! Determine the monthly sums and write them to output.
 DO I = 1, 12 * NrOfYears
!  Determine which month it is in a four-year period.
   CyclicMonth = MOD(I, 48)
   IF (CyclicMonth == 0) CyclicMonth = 48
!  Determine in which four-year period month nr I occurs.
   FourYearPeriodNr = FLOOR(REAL(I) / 48.0)
   MonthlySum(I,1) = ANINT(REAL(I))
!  Determine the first and last day of the month of interest.
   IF (CyclicMonth < 48) THEN
     StartDay = (FourYearPeriodNr * 1461) + NINT(MonthStartTimes(CyclicMonth)) + 1
     EndDay = (FourYearPeriodNr * 1461) + NINT(MonthStartTimes(CyclicMonth + 1))
   ELSE
     StartDay = ((FourYearPeriodNr - 1) * 1461) + NINT(MonthStartTimes(CyclicMonth)) + 1
     EndDay = ((FourYearPeriodNr - 1) * 1461) + NINT(MonthStartTimes(CyclicMonth + 1))
   END IF
   MonthlySum(I,2) = SUM(DailySum(StartDay:EndDay, 2))
   WRITE(26, 520) MonthlySum(I,1), MonthlySum(I,2)
 END DO
! Determine the annual sums and write them to output.
 DO I = 1, NrOfYears
   AnnualSum(I,1) = ANINT(REAL(I))
!  The first month of year I is [(I - 1) * 12 + 1]. The last month of year I is (I * 12).
   AnnualSum(I, 2) = SUM(MonthlySum(((I - 1) * 12) + 1 : I * 12, 2)) 
   WRITE(27, 520) AnnualSum(I,1), AnnualSum(I,2)
 END DO
520 FORMAT (1X, F15.1, 1X, ES12.4)
 CLOSE (26)
 CLOSE (27)
! Compute the mean of the monthly and annual rainfall.
 MonthlyMean = 0.0_DP
 AnnualMean = 0.0_DP
 DO I = 0, NrOfYears - 1
   AnnualMean = AnnualMean + AnnualSum(I + 1, 2)
   DO J = 1, 12 
     MonthlyMean(J) = MonthlyMean(J) + MonthlySum((12 * I) + J, 2)
   END DO
 END DO
 AnnualMean = AnnualMean / REAL(NrOfYears, DP)
 MonthlyMean = MonthlyMean / REAL(NrOfYears, DP)
! Compute the standard deviation of monthly and annual rainfall. Use the unbiased estimator.
 MonthlySD = 0.0_DP
 AnnualSD = 0.0_DP
 DO I = 0, NrOfYears - 1
   AnnualSD = AnnualSD + ((AnnualSum(I + 1, 2) - AnnualMean)**2.0_DP)
   DO J = 1, 12 
     MonthlySD(J) = MonthlySD(J) + ((MonthlySum((12 * I) + J, 2) - MonthlyMean(J))**2.0_DP)
   END DO
 END DO 
 AnnualSD = SQRT(AnnualSD / (REAL(NrOfYears, DP) - 1.0_DP))
 DO J = 1, 12
   MonthlySD(J) = SQRT(MonthlySD(J) / (REAL(NrOfYears, DP) - 1.0_DP))
 END DO
! Write the descriptive monthly and annual statistics to output.
 OPEN (28, FILE = 'Rain_Statistics.OUT', STATUS = 'UNKNOWN')
 WRITE(28, 530) AnnualMean, AnnualSD
530 FORMAT(/,' The mean annual rainfall is: ',ES12.4,'. The standard deviation of the annual rainfall is: ',ES12.4,'.',/)
 WRITE(28, 540)
540 FORMAT(' Month     Mean rainfall  Standard deviation',/,' ___________________________________________')
 WRITE(28, 550) MonthlyMean(1), MonthlySD(1)
 WRITE(28, 551) MonthlyMean(2), MonthlySD(2)
 WRITE(28, 552) MonthlyMean(3), MonthlySD(3)
 WRITE(28, 553) MonthlyMean(4), MonthlySD(4)
 WRITE(28, 554) MonthlyMean(5), MonthlySD(5)
 WRITE(28, 555) MonthlyMean(6), MonthlySD(6)
 WRITE(28, 556) MonthlyMean(7), MonthlySD(7)
 WRITE(28, 557) MonthlyMean(8), MonthlySD(8)
 WRITE(28, 558) MonthlyMean(9), MonthlySD(9) 
 WRITE(28, 559) MonthlyMean(10), MonthlySD(10)
 WRITE(28, 560) MonthlyMean(11), MonthlySD(11)
 WRITE(28, 561) MonthlyMean(12), MonthlySD(12)  
550 FORMAT(' January    ',ES12.4,8X,ES12.4)
551 FORMAT(' February   ',ES12.4,8X,ES12.4)
552 FORMAT(' March      ',ES12.4,8X,ES12.4)
553 FORMAT(' April      ',ES12.4,8X,ES12.4)
554 FORMAT(' May        ',ES12.4,8X,ES12.4)
555 FORMAT(' June       ',ES12.4,8X,ES12.4)
556 FORMAT(' July       ',ES12.4,8X,ES12.4)
557 FORMAT(' August     ',ES12.4,8X,ES12.4)
558 FORMAT(' September  ',ES12.4,8X,ES12.4)
559 FORMAT(' October    ',ES12.4,8X,ES12.4)
560 FORMAT(' November   ',ES12.4,8X,ES12.4)
561 FORMAT(' December   ',ES12.4,8X,ES12.4)
 CLOSE (28)
 Deallocate (DailySum)
 Deallocate (HourlySum)
 Deallocate (MonthlySum)
 Deallocate (AnnualSum)
!
 END SUBROUTINE OutputWriter
!
 SUBROUTINE RainMaker
! Rainmaker determines the duration of each storm and generates the raincells in them.
!
! Interface: The number of storms (NrOfStorms) and their start times (StormStartTimesOld) are obtained from Module StormParameters.
! Module Parameters (ModelType) informs the subroutine if the rainfall rate is gamma or epxonentially distributed. It also provides
! the number of periods per year (NrOfPeriods) for which parameters are provided and the number of years (NrOfYears) for which
! rainfall needs to be generated. The arrays with the various model parameters for all periods are also obtained from Module
! Parameters.These are used to generate storm-specific parameters.
! The start time, end time, and rainfall rate of each rain cell is passed on to other program elements through Module
! CellParameters (AllCellsStartEndDepthOld). This module also contains the total amount of rain calculated from the individual
! cells (TotalAmountOfRain).
! Information about the required precision of double precision real variables is obtained from Module DoublePrecision.
!
! Gerrit H. de Rooij, February 2018
!
 USE DoublePrecision
 USE Parameters, ONLY: Output, ModelType, NrOfYears, NrOfPeriods, alpha, kappa, lower_eta, murate, mux, nu, phi
 USE StormParameters
 USE CellParameters
 IMPLICIT NONE
 INTEGER :: Help, Help2, I, J
 REAL (KIND = DP) :: alpha_p, CellInterval, Duration, Help3, lower_eta_p, kappa_p, murate_p, mux_p, nu_p, phi_p, Time
 REAL (KIND = DP), Dimension (:), Allocatable :: beta, eta, gamma, RainPerCell, StormEndTimes
! beta, eta, and gamma store parameters that vary from storm to storm.
 REAL (KIND = DP), Dimension (:,:), Allocatable ::  AllCellsStartEndDepth, CellStartEndDepth, CellStartEndDepthOld
!
! With the total number of storms available, allocate the arrays with storm-specific parameters beta, eta, and gamma, as well as
! the array that stores the times when storms end.
 Allocate (beta(1:NrOfStorms))
 Allocate (eta(1:NrOfStorms))
 Allocate (gamma(1:NrOfStorms))
 Allocate (StormEndTimes(1:NrOfStorms))
! The array NrOfCells stores the number of cells for every storm in the simulation period.  
 Allocate (NrOfCells(1: NrOfStorms))
! In the code below, the number of cells per storm will be determined. Since the total number of cells is required every time 
! a storm has been processed, i.e.,the number of cells in it has been determined, the sum operation should not be influenced by 
! array entries pertaining to later storms. Therefore, initialize this array to avoid this.
 NrOfCells = 0
! 
! With the full set of storms available, use alpha and nu for each storm to randomly draw a value for the gamma-distributed eta. 
! With the values of kappa and phi for each storm, determine beta and gamma for each storm. Use those to generate the duration 
! and the number of cells in it.
! Help stores the period in which the previous storm started to avoid having to run through all periods for each storm.
 Help = 1
 AllStorms: DO I = 1, NrOfStorms
!  Find the period in which the storm starts, and adopt the proper values of alpha and nu. Then generate a value of eta for 
!  each storm, and calculate the corresponding values of beta and gamma (see Rodriguez-Iturbe et al., 1988).
   DO J = Help, (NrOfYears * NrOfPeriods)
!    Determine the time period J in which storm I starts.  
     IF ((StormStartTimesOld(I) >= AllStartTimes(J)) .AND. (StormStartTimesOld(I) < AllStartTimes(J + 1))) THEN
!      Find the correct values of alpha and nu, use them to generate a value of eta, and use that to determine beta and gamma
!      using the appropriate values for kappa and phi.
       Help2 = MOD(J, (4 * NrOfPeriods))
       IF (Help2 == 0) Help2 = 4 * NrOfPeriods
       alpha_p = alpha(Help2)
       nu_p = nu(Help2)
       lower_eta_p = lower_eta(Help2)
!      Subroutine DrawGammaValue generates a random variate from a gamma distribution with parameters alpha and nu.
       CALL DrawGammaValue(alpha_p, 1.0_DP/nu_p, lower_eta_p, eta(I)) 
!      With eta known, the values of beta and gamma can be computed using the values of kappa and phi valid for the period 
!      in which the storm occurs.      
       kappa_p = kappa(Help2)
       phi_p = phi(Help2)
       beta(I) = kappa_p * eta(I)
       gamma(I) = phi_p * eta(I)
!      For the calculation of the cell depth, the appropriate value(s) of mux (and murate) is/are required.
       mux_p = mux(Help2)
       IF (ModelType == 'Gam') murate_p = murate(Help2)
!      Since the next storm cannot precede the current storm, there is no need to search earlier periods. Therefore, update 
!      the period in which to start the search.
       Help = J
       EXIT
     END IF
   END DO
!  Determine the duration of storm I. With this, the time the storm ends is computed and stored in array StormEndTimes.
!  Subroutine StormDuration generates a random number from an exponential distribution with parameter gamma(I) and adds that to the
!  starting time of the storm to give the time it ends.
   CALL StormDuration (gamma(I), StormStartTimesOld(I), StormEndTimes(I))
!  Determine the number of cells and their starting times in each storm. The algorithm to do so is comparable to that for the 
!  generation of storms above. In addition, generate the duration and cell depth for each cell.
   Time = 0.0_DP
!  Each storm starts with a cell, so the minimum value of all array elements of NrOfCells is 1.
   NrOfCells(I) = 1
!  Initializing allocation of the array and the back-up array that store the starting and ending times of each cell as well as the
!  cell depth (rainfall rate).
   Allocate (CellStartEndDepth (1:1, 1:3))
   Allocate (CellStartEndDepthOld (1:1, 1:3))
!  The storm starts with a rain cell.
   CellStartEndDepth(1,1) = StormStartTimesOld(I)
!  Generate the duration of the first cell.
!  Subroutine StormCell generates a random number from an exponential distribution with parameter eta(I).
   CALL StormCell (eta(I), Duration)  
   CellStartEndDepth(1,2) = StormStartTimesOld(I) + Duration
!  Generate the rainfall rate (depth) of the first cell.
   IF (ModelType == 'Gam') THEN
!    Generate the rainfall rate from a gamma distribution. Note that mux_p is the shape parameter of the Gamma-distributed 
!    rainfall rate.
     CALL DrawGammaValue(mux_p, 1.0_DP/murate_p, 0.0_DP, CellStartEndDepth(1,3))
   ELSE
!    Generate a rainfall rate from an exponential distribution.
     CALL StormCell(1.0_DP/mux_p, CellStartEndDepth(1,3))
   END IF
!  Copy this to the backup array.
   CellStartEndDepthOld = CellStartEndDepth 
!  The open DO-loop CellGenerator generates cells for storm I until the starting time of a cell exceeds the end time of a storm.  
   CellGenerator: DO
     CALL StormCell(beta(I), CellInterval) 
!    Express the starting time of a cell relative to the start of the storm to which it belongs.
     Time = Time + CellInterval
!    The generated cell only needs to be stored if its starting time does not exceed the end time of its storm.
     IF (Time < (StormEndTimes(I) - StormStartTimesOld(I))) THEN
!      Increase the size of the array for the cell properties by one, and copy the previous cells in it.
       NrOfCells(I) = NrOfCells(I) + 1
       Deallocate (CellStartEndDepth)
       Allocate (CellStartEndDepth(1:NrOfCells(I), 1:3)) 
       CellStartEndDepth(1:NrOfCells(I) - 1, 1:3) = CellStartEndDepthOld 
!      Store the starting time of the latest cell.       
       CellStartEndDepth(NrOfCells(I), 1) = StormStartTimesOld(I) + Time
!      Generate the duration and the depth of this cell.
       CALL StormCell(eta(I), Duration) 
       IF (ModelType == 'Gam') THEN
!        Generate the rainfall rate from a gamma distribution.
         CALL DrawGammaValue(mux_p, 1.0_DP/murate_p, 0.0_DP, CellStartEndDepth(NrOfCells(I), 3))
       ELSE
!        Generate a rainfall rate from an exponential distribution.
         CALL StormCell (1.0_DP/mux_p, CellStartEndDepth(NrOfCells(I), 3))
       END IF  
!      Store the end time.
       CellStartEndDepth(NrOfCells(I), 2) = CellStartEndDepth(NrOfCells(I), 1) + Duration
!      Increase the back-up array and store the cell properties in it so CellStartEndDepth can be increased in the next cycle in
!      the unbounded DO loop.
       Deallocate (CellStartEndDepthOld)
       Allocate (CellStartEndDepthOld(1:NrOfCells(I), 1:3))
       CellStartEndDepthOld = CellStartEndDepth
     ELSE
!      The number of cells and their properties have been determined for storm I. Move on to the next storm - revert control back
!      to the outer DO-Loop (AllStorms). 
       IF (Output == 'Yes') WRITE (22, 300) I, StormStartTimesOld(I), StormEndTimes(I), NrOfCells(I)
       EXIT
     END IF    
   END DO CellGenerator
!  In the next cycle of the outer DO Loop (AllStorms), the arrays with the cell properties will be reset. Store them in arrays
!  AllCellsStartEndDepth and AllCellsStartEndDepthOld.  
   Allocate (AllCellsStartEndDepth(1:SUM(NrOfCells), 1:3))
   IF (I == 1) THEN
     Allocate (AllCellsStartEndDepthOld(1:NrOfCells(1), 1:3))
     AllCellsStartEndDepth(1:NrOfCells(1), 1:3) = CellStartEndDepth
     AllCellsStartEndDepthOld = CellStartEndDepth
   ELSE
!    Copy the cell properties of the previous storms into AllCellsStartEndDepth.
     AllCellsStartEndDepth(1 : (SUM(NrOfCells) - NrOfCells(I)), 1:3) = AllCellsStartEndDepthOld    
!    Add the cell properties of the current storm (storm I).
     AllCellsStartEndDepth((SUM(NrOfCells) - NrOfCells(I) + 1) : SUM(NrOfCells), 1:3) = CellStartEndDepth
   END IF
!  With the old starting times copied into the newly increased array, the backup array can be increased in turn.
   Deallocate (AllCellsStartEndDepthOld)
   Allocate (AllCellsStartEndDepthOld(1:SUM(NrOfCells), 1:3))
!  The backup array is ready to receive the latest results. When this is done, array AllCellsStartEndDepth can be deallocated so it
!  can be reallocated after the cells for the next storm have been generated. 
   AllCellsStartEndDepthOld = AllCellsStartEndDepth
   Deallocate (AllCellsStartEndDepth)
!  Now that the cell properties have been properly stored, the arrays used in DO loop CellGenerator can be deallocated to allow them
!  to be reset for the next storm.
   Deallocate (CellStartEndDepth)
   Deallocate (CellStartEndDepthOld)
 END DO AllStorms
! As a balance check, calculate the total rainfall from all cells. For efficiency, store the total number of cells in Help.
 Help = SUM(NrOfCells)
 Allocate (RainPerCell(1:Help))
! Note that at this stage, array AllCellsStartEndDepth is deallocated, but its contents are stored in its backup array
! AllCellsStartEndDepthOld. Use that array instead.
 DO I = 1, Help
   RainPerCell(I) = (MIN((REAL((NrOfYears / 4), DP) * 1461.0_DP), AllCellsStartEndDepthOld(I,2))&
                   & - MIN((REAL((NrOfYears / 4), DP) * 1461.0_DP), AllCellsStartEndDepthOld(I,1))) * AllCellsStartEndDepthOld(I,3)
!  Debug start
   IF ((RainPerCell(I) < 0.0_DP) .OR. (RainPerCell(I) > 1000.0_DP)) THEN
     WRITE(22, 310) I, RainPerCell(I), (AllCellsStartEndDepthOld(I,J), J = 1, 3)
   END IF
!  Debug end
 END DO
 TotalAmountOfRain = SUM(RainPerCell)
 WRITE (22, 320) TotalAmountOfRain, NrOfStorms, Help !Write the total amount of rain, the number of storms, and the number of cells
!                                                     to the output file.
 Help3 = MAXVAL(RainPerCell)
 WHERE (RainPerCell == 0.0_DP) RainPerCell = Help3
 WRITE (22, 330) MINVAL(RainPerCell), Help3 ! Write the range of rainfall amounts to the output file.
300 FORMAT(' Storm ',I10,' starts at ',ES12.4,', ends at ',ES12.4,', and has ',I10,' cells.')
310 FORMAT(' Cell ',I10,' has ',ES12.4,' units of rainfall. It starts at time ',ES12.4,', ends at ',ES12.4,&
          &', and has rainfall rate ',ES12.4,'.')
320 FORMAT(/,' A total amount of rain of ', ES12.4,' fell in ',I10,' storms that had a total of ',I7,' cells.')
330 FORMAT(' The smallest non-zero cell produced ',ES12.4,' units of rainfall, the largest cell produced ',ES12.4,' units.')
 Deallocate (beta)
 Deallocate (eta)
 Deallocate (gamma)
 Deallocate (StormEndTimes)
!
 END SUBROUTINE RainMaker
!
 SUBROUTINE Sorter(ArrayToBeSorted, SortedArray)
! Sorter sorts a 1-dimensional array according to the elements in its first dimension.
! In this application, the length of the array is determined by the number of rain cells between dry intervals separating rain
! storms. The number of cells involved will be limited for most climates, and many of them will already be in the proper order.
! For those reasons, a relatively inefficient algorithm that is more efficient than more powerful  algorithms if the array is
! better sorted is implemented based on Shell's method as described in Press, W.H., Teukolsky, S.A., Vetterling, W.T., and
! Flannery, B.P. 1996. Numerical Recipes in Fortran 90, 2nd edition, Cambridge University Press, Cambridge, UK (p.1168), and in
! http://en.wikipedia.org/wiki/Shellsort.
!	
! Interface: The unsorted array is received through the argument list, and the sorted array is passed back through SortedArray
! in the argument list. The second columns of SortedArray contains the original ranking of then unsorted array. This is helpful
! if the sorted array was part of a larger one that needs to be sorted according to the column passed on to Sorter. 
! The precision of the REAL variables for which a high precision is required is provided by MODULE Precision.
!
! Gerrit H. de Rooij, January 2018
!
 USE DoublePrecision
 IMPLICIT NONE
 INTEGER :: L, M, N, NrOfElements, P, Q, Step, UpperBound
 REAL (KIND = DP), Dimension (1:2) :: New   
 REAL (KIND = DP), Intent(in), Dimension (:) :: ArrayToBeSorted
 REAL (KIND = DP), Intent(inout), Dimension (:,:) :: SortedArray
 REAL (KIND = DP), Dimension (:), Allocatable :: Increment, IncrementOld
 REAL (KIND = DP), Dimension (:,:), Allocatable :: Dummy   
!
! Start by finding the series of increments of data pairs that need to be sorted. The increment series is defined according to
! Tokuda's sequence (see the Wikipedia reference). 
 UpperBound = SIZE(ArrayToBeSorted) 
 Allocate (Increment(1))
 Allocate (IncrementOld(1))
 L = 1
 Increment(1) = 1.0_DP
 DO
   L = L + 1
   IncrementOld = Increment
   Deallocate (Increment)
   Allocate (Increment(L))
   Increment(1: L - 1) = IncrementOld
   Increment(L) = (2.25_DP * Increment(L - 1)) + 1.0_DP
   IF (INT(Increment(L)) > UpperBound) THEN
     EXIT
   END IF
   Deallocate (IncrementOld)
   Allocate (IncrementOld(L))
 END DO
 Deallocate (IncrementOld)
!Start with the sort according to the largest increment and move down to increment 1.
 Sort: DO M = L - 1, 1, -1
!  The size of the increment is determined by Step. It decreases with every passage through the DO loop.          
   Step = INT(Increment(M))
!  Find the subarrays that need to be sorted for this increment and their sizes. The number of subarrays equals the step size. 
   DO N = 1, Step
!    N.B. The division in the next statement contains integers only, and the result will itself be an integer. As the number of
!    the starting element of the subarray increases, the subarray tends to become smaller, as fewer Steps fit in the remainder
!    of the array.
     NrOfElements = ((UpperBound - N) / Step) + 1
!    Prepare the array Dummy that will be sorted.
     Allocate (Dummy(1:NrOfElements, 1:2))
     DO P = 0, NrOfElements - 1
        Dummy(P + 1, 1) = ArrayToBeSorted(N + (P * Step))  
        Dummy(P + 1, 2) = REAL(N + (P * Step), DP)  ! This stores the original position of the element that will be sorted. 
     END DO
!    Sort Dummy and pass on the sorted results to SortedArray. The sorting algorithm is straight insertion (see Press, W.H.,
!    Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P. 1992. Numerical Recipes in Fortran, 2nd edition, Cambridge 
!    University Press, Cambridge, UK (p.321). 
!    In the first step, the two left-most array elements are sorted. Then, the third is placed in the correct position with 
!    respect to these two, then the fourth with respect to the first three, and so on.
     DO P = 1, NrOfElements - 1
       New = Dummy(P + 1, 1:2) 
       DO Q = P, 1, -1
!        Run down the already sorted array elements from largest to smallest. Every time when New is smaller than an element,
!        this element is moved one position to the right. If an element is smaller than New, it is inserted in the slot 
!        immediately after that element.
         IF (New(1) < Dummy(Q,1)) THEN
           Dummy(Q + 1, 1:2) = Dummy(Q,1:2) 
           IF (Q == 1) Dummy(1, 1:2) = New   
!          New is smaller than all previously sorted elements. Place it in the first slot.
         ELSE
           Dummy(Q + 1, 1:2) = New 
           EXIT
         END IF
       END DO        
     END DO
     DO P = 0, NrOfElements - 1
       SortedArray((N + (P * Step)), 1:2) = Dummy(P + 1, 1:2)
     END DO
     Deallocate (Dummy)
   END DO
 END DO Sort
 Deallocate (Increment)
 END SUBROUTINE Sorter
!
 SUBROUTINE StormCell (PdfParameter, Variate)
! StormCell generates an exponentially distributed variate.
! 
! Interface: The parameter of the exponential distribution is obtained on input through PdfParameter. The generated 
! variate is passedon to the main program through Variate.
! The subroutine calls Random_Number once to obtain a variate from the standard uniform distribution. The unit of PdfParameter
! is  1/unit of the time coordinate,so 1/day in this code. The precision of the REAL variables for which a high precision is
! required is provided by MODULE Precision.
!
! Gerrit H. de Rooij, October 2014, modified January 2018
!
 USE DoublePrecision
 IMPLICIT NONE
 REAL (KIND = DP), Intent(in) :: PdfParameter
 REAL (KIND = DP), Intent(out) :: Variate
 REAL (KIND = DP) :: W
! Acquire a variate from the standard uniform distribution and transform it to an exponentially distributed variate
! (inverse transform sampling)
 CALL Random_Number(W)
 Variate = -LOG(W) / PdfParameter
 END SUBROUTINE StormCell
!
 SUBROUTINE StormDuration (PdfParameter, StormStartTime, StormEndTime)
! Subroutine StormDuration generates a random number from an exponential distribution with parameter gamma(I) and adds that to the
! starting time of the storm to give the time it ends.
! 
! Interface: The parameter of the exponential distribution is obtained on input through PdfParameter and passed on to routine
! StormCell. StormCell passes on the duration of the storm of interest through variable StormDuration. This duration is added to the
! starting time of the storm obtained from the main program through variable StormStartTime. The calculated end time of the storm is
! passed on to the main program through variable StormEndTime. The precision of the REAL variables for which a high precision is
! required is provided by MODULE Precision.
!
! Gerrit H. de Rooij, October 2014, modified January 2018
!
 USE DoublePrecision
 IMPLICIT NONE
 REAL (KIND = DP), Intent(in) :: PdfParameter, StormStartTime
 REAL (KIND = DP), Intent(out) :: StormEndTime
 REAL (KIND = DP) :: DurationOfStorm
!
 CALL StormCell(PdfParameter, DurationOfStorm)
 StormEndTime = StormStartTime + DurationOfStorm
 END SUBROUTINE StormDuration
!
 SUBROUTINE StormStarts
! StormStarts generates the start times of all storms.
 USE DoublePrecision
 USE Parameters, ONLY: Output, NrOfFourYearPeriods, NrOfPeriods, NrOfYears, lambda, StartTimes
 USE StormParameters
! Variable declarations
 CHARACTER (LEN = 3) :: FirstStorm
 INTEGER :: I, J
 REAL (KIND = DP) :: lambda_p, PeriodDuration, StormInterval, Time
 REAL (KIND = DP), Dimension (:), Allocatable :: StormStartTimes
!
! For each period, generate a series of storms. The total number of periods equals NrOfYears * NrOfPeriods. 
! N.B. The latest storm is defined as the first storm that has its starting time in a later period (month) than the one that is
! currently simulated. No rainfall is generated for this storm, but its start time is stored. The value of lambda is updated to
! reflect the period (month) in which it has its starting point and a new set of storms starting at that time is generated using 
! the updated value of lambda and a new set of random numbers.
! Note that this algorithm permits an entire period (month), and even multiple consecutive periods, to be without rain.
!
! Define the start times of all periods in the simulation period.
 Allocate (AllStartTimes (1 : (NrOfYears * NrOfPeriods) + 1))
 IF (NrOfFourYearPeriods == 1) THEN
   AllStartTimes = StartTimes
 ELSE
   AllStartTimes(1:4 * NrOfPeriods) = StartTimes(1:4 * NrOfPeriods)
   DO J = 2, NrOfFourYearPeriods 
     DO I = 1, 4 * NrOfPeriods
       AllStartTimes(((J - 1)* 4 * NrOfPeriods) + I) = StartTimes(I) + (1461.0_DP * (J - 1))
     END DO
   END DO
   AllStartTimes((NrOfYears * NrOfPeriods) + 1) = REAL(NrOfFourYearPeriods, DP) * 1461.0_DP
 END IF
!
! Generate storm starting times within each period using the lambda values supplied on input. Check in which later period the
! last storm generated for a given period starts.
! N.B. At the end of this DO construct, the starting times of all storms are stored in StormStartTimesOld, while StormStartTimes
! is deallocated at the time the loop is terminated.
 Time = 0.0_DP
 NrOfStorms = 0 
 AllPeriods: DO I = 1, NrOfYears * NrOfPeriods
   PeriodDuration = AllStartTimes(I + 1) - AllStartTimes(I)
!  Check if period I contains a storm. If not, move to the next period. This IF construct also terminates the DO loop AllPeriods if
!  there is no storm in the final period.
   IF (Time > PeriodDuration) THEN
     Time = Time - PeriodDuration
     CYCLE AllPeriods
   END IF 
!  Determine which value of lambda to use. J gives the number of the period within the first 4 years that corresponds to period I
!  anywhere in the simulation period.
   J = MOD(I, (4 * NrOfPeriods))
   IF (J == 0) J = 4 * NrOfPeriods !The Modulo operation only returns the value zero for the last period in each four-year cycle.
   lambda_p = lambda(J)  
!  The DO loop AllPeriods keeps generating storms in period I until the starting time of a storm exceeds the end of the period. 
!  The array StormStartTimes is increased with every generated storm and stores the starting times of all storms. Note that the 
!  last time stored is of a storm that is outside the simulation period and should be ignored. 
!  The backup array StormStartTimesOld stores the starting times generated so far while the size of StormStartTimes is increased. 
!  Once the values have been copied back into the increased array, the backup array is deallocated and increased and allocated 
!  again. The arrays keep on increasing and storing starting times with respect to the start of the simulation period until the 
!  DO loop AllPeriods cycled through all the periods (counted by I) in the total simulation period.
   FirstStorm = 'Yes' ! This variable equals 'Yes' only if DO-loop StartStorm is in its first cycle for period I.
   StartStorm: DO
!    Increase the size of the array storing the starting times of the storms by one element to hold the next storm. The final run 
!    through this unbounded DO loop will generate a storm starting time that lies outside of the simulation period, and thus will 
!    not be required later in the program. The IF construct at the start of the DO-loop AllPeriods will respond to this storm 
!    starting time by increasing I to its final value and then orderly exiting the outer DO loop. At that time, StormStartTimes is 
!    in deallocated state, and all storm starting times but the one(i.e., all that are within the simulation period) are retained 
!    in the backup array StormStartTimesOld.
     Allocate (StormStartTimes (1:NrOfStorms + 1))
!    Feed all elements of this array except the last one with the previously generated storm starting times stored in the backup
!    array.
     IF (NrOfStorms >= 1) THEN
       StormStartTimes(1:NrOfStorms) = StormStartTimesOld 
       Deallocate (StormStartTimesOld)
     END IF
!    Subroutine StormCell randomly generates a number from an exponential distribution with parameter lambda_p.
     CALL StormCell(lambda_p, StormInterval) 
     NrOfStorms = NrOfStorms + 1
!    The starting time of the newest storm is calculated from that of the previous storm and the interval between the two, and
!    stored.
     IF (NrOfStorms == 1) THEN
       StormStartTimes(1) = StormInterval
     ELSE
       StormStartTimes(NrOfStorms) = StormStartTimes(NrOfStorms - 1) + StormInterval
     END IF
!    The array is copied into a backup array so it can be deallocated to be increased by one element for the next iteration.
     Allocate (StormStartTimesOld(1:NrOfStorms))
     StormStartTimesOld = StormStartTimes
     Deallocate (StormStartTimes)
     Time = Time + StormInterval
     IF (Time > PeriodDuration) THEN
       IF ((FirstStorm == 'Yes') .AND. (Output == 'Yes')) WRITE(22, 210) I, StormInterval
!      The next statements ensure that the first storm in a later period starts at the time determined by the last call of 
!      StormCell within DO loop AllPeriods.
       Time = Time - PeriodDuration
       EXIT StartStorm
     END IF
     FirstStorm = 'No'
   END DO StartStorm
 END DO AllPeriods
! The last storm generated will always be outside the simulation period and should not be taken into account.
 NrOfStorms = NrOfStorms - 1
! If the final storm just within the simulation period starts exactly at the end of the simulation period it should be ignored
! as well.
 IF (StormStartTimesOld(NrOfStorms) >= AllStartTimes((NrOfYears * NrOfPeriods) + 1)) NrOfStorms = NrOfStorms - 1
 WRITE(22, 200) NrOfStorms
200 FORMAT (/,' Start times of ',I12,' storms generated.')
210 FORMAT (' Period ',I10,' is followed by a dry period. The time between its last storm and the next is ',ES12.4,'.')
! 
 END SUBROUTINE StormStarts
!
 SUBROUTINE Timeline
! Timeline takes the array that stores all rain cells with their start times, end times, and rainfall rates and sorts the cells
! in order of ascending start times. From the sorted arrays it derives the detailed rainfall record: every time a rain cell 
! starts or stops, a time is set and the rainfall rate for the time period that ends at that time is determined by adding the 
! rainfall rates of all acitve cells during that time interval. When that is completed, the routine adds the amounts of rainfall
! of all time intervals and compares the result with sum of the rainfall of all cells calculated by Subroutine RainMaker.
!
! Interface: The start and end times of the rain cells are taken from array AllCellsStartEndDepthOld in Module CellParameters.
! The rainfall rates of the cells are obtained from the same array. The number of cells per storm (NrOfCells) is taken from
! Module CellParameters. The number of storms (NrOfStorms) is obtained from Module StormParameters. The number of years
! (NrOfYears) of the rainfall record and the number of periods into which a year is divided (NrOfPeriods) are taken from
! Module Parameters. The total amount of rain (TotalAmountOfRain) is obtained from Module CellParameters. The required
! precision of the REAL variables id obtained from Module DoublePrecision.
! Informative output is written to file Diagnositic.OUT. The time line of rainfall (the start and end times of all intervals
! with a constant rainfall rate, together with their rainfall rates are provided to other program units through array
! RainfallRate in Module RainfallTimeline.
!
! Gerrit H. de Rooij, February 2018.
!
 USE DoublePrecision
 USE Parameters, ONLY: NrOfFourYearPeriods
 USE CellParameters
 USE StormParameters, ONLY: NrOfStorms
 USE RainfallTimeline, ONLY: NrOfIntervals, RainfallRate
!
! Variable declarations
 IMPLICIT NONE
 INTEGER :: ClusterEndStorm, ClusterStartStorm, EndCell, EndingCell, Help, I, J, NrStart, NrStartOld, OldestActiveCell,&
          & StartCell
 REAL (KIND = DP) :: ClusterStartTime, NextCellStart, SoonestEnd, RelativeError, Time
 REAL (KIND = DP), Dimension (:), Allocatable :: RainPerInterval
 REAL (KIND = DP), Dimension (:,:), Allocatable :: RainfallRateOld, SortedArray, StartEnd, UnsortedSegment
! N.B. The cells appear clustered by storm in array AllCellsStartEndDepthOld. If storms overlap (and they frequently do in humid 
! climates) this means the cells are not in chronological order, which causes considerable problems when a time line of the 
! rainfall record is constructed below. The array therefore needs to be sorted in order of ascending cell starting times.
!
! Prepare an array of the starting times of the first and last cell of each storm.
 StartCell = 0
 Allocate (StartEnd(1:NrOfStorms, 1:2))
 DO I = 1, NrOfStorms
   StartEnd(I, 1) = AllCellsStartEndDepthOld(StartCell + 1, 1)
   StartCell = StartCell + NrOfCells(I)
   StartEnd(I, 2) = AllCellsStartEndDepthOld(StartCell, 1)
 END DO
! Starting with the last storm, start looking backward for cells overlapping with this storm, and for cells overlapping with 
! those storms to find the first storm of the cluster. Once a cluster of overlapping storms has been identified, the segment of 
! the array containing the cells of those storms is sorted.
! Then, the storm preceding the first storm of the cluster is taken as the starting point for the cluster preceding it.
 ClusterEndStorm = NrOfStorms
 OverlappingStorms: DO
! Find the starting time of the storm that ends the cluster.
   ClusterStartTime = StartEnd(ClusterEndStorm, 1)
! Find the sequence of storms overlapping with one another that ends with ClusterEndStorm.
   ClusterStartStorm = ClusterEndStorm
   DO I = ClusterEndStorm - 1, 1, -1
     IF (StartEnd(I, 2) > ClusterStartTime) THEN
       ClusterStartTime = StartEnd(I, 1)
       ClusterStartStorm = I
     END IF
   END DO
   IF (ClusterStartStorm == ClusterEndStorm) THEN
!    This storm is not part of a cluster, start search for clusters from the storm preceding it.
     ClusterEndStorm = ClusterEndStorm - 1
     IF (ClusterEndStorm == 1) THEN
       EXIT
     ELSE
       CYCLE
     END IF
   ELSE
!    There is a cluster. Sort the cells of the showers that comprise the cluster. Then proceed the search for other clusters
!    taking the storm preceding this cluster as a starting point.
!    Identify the first and last cell of the cluster.
     StartCell = 1
     DO I = 1, ClusterStartStorm - 1
       StartCell = StartCell + NrOfCells(I)
     END DO
     EndCell = StartCell - 1
     DO I = ClusterStartStorm, ClusterEndStorm
       EndCell = EndCell + NrOfCells(I)
     END DO
!    Sort the segment of array AllCellsStartEndDepthOld between StartCell and Endcell (inclusive) according to the cell starting
!    times (AllCellsStartEndDepthOld(1,:,:)).
!    This involves cells nr. FirstCell through FirstCellOfNextStorm - 1.
     Allocate (SortedArray(1:EndCell - StartCell + 1, 1:2)) 
     Allocate (UnsortedSegment(1:EndCell - StartCell + 1, 1:3)) 
     DO I = 1, EndCell - StartCell + 1 
       SortedArray(I,1) = -AllCellsStartEndDepthOld(I + StartCell - 1, 1) 
       SortedArray(I,2) = -I ! After the sorting operaton below all values should be positive.
       UnsortedSegment(I, 1:3) = AllCellsStartEndDepthOld(I + StartCell - 1, 1:3) 
     END DO
     CALL Sorter(AllCellsStartEndDepthOld(StartCell:EndCell, 1), SortedArray)
     DO I = 1, EndCell - StartCell + 1 ! 
       AllCellsStartEndDepthOld(I + StartCell - 1, 1:3) = UnsortedSegment(INT(SortedArray(I,2)), 1:3) 
     END DO 
     Deallocate (SortedArray) 
     Deallocate (UnsortedSegment) 
!    Exit the loop if the first storm was reached. If not, start the search for the next cluster.
     IF (ClusterStartStorm == 1) THEN
       EXIT
     ELSE
       ClusterEndStorm = ClusterStartStorm - 1
     END IF   
   END IF
 END DO OverlappingStorms
!
! From the times in AllCellsStartEndDepthOld create a time line with all instances of a cell beginning and ending (N.B. Overlap is
! likely, identical times for multiple changes are possible.) For each interval of that time line, sum the depths of the active
! cells to arrive at the rainfall rate as a function of time.
! Array RainfallRate stores the times at which the rainfall rate changes and the rainfall rate for time periods between these
! changes. The first entries in both rows are zero (zero rainfall at the start of the simulation period). Note that the first column
! in RainFallRate is numbered zero. The rainfall rate stored in RainFallRate(I,2) is valid for the time period between
! RainfallRate(I-1,1) and RainfallRate(I,1).
!
! Allocate the backup version of the array for the time line. The version for the updates will be allocated in the DO loop below.
 Allocate (RainfallRateOld(0:1, 1:2))
! No rain at the start of the simulation period.
 RainfallRateOld(0,1) = 0.0_DP
 RainfallRateOld(0,2) = 0.0_DP
! Enter the first dry spell until the start of cell 1.
 RainfallRateOld(1,1) = AllCellsStartEndDepthOld(1,1)
 RainfallRateOld(1,2) = 0.0_DP
! Initialize the variable that stores the time at which the next cell starts producing rain with the starting time of the first
! cell.
 NextCellStart = AllCellsStartEndDepthOld(1,1)
! Initialize counters for the unbounded DO loop 'Timeline'. NrStart stores the total number of cells that have become active (and
! possibly inactive again) at the time for which the DO Loop is executed. It therefore equals the number of the latest cell that has
! become active. Counter I reflects the number of the time interval on the time line as it is being constructed.
 NrStart = 1
 NrStartOld = 1
 Help = SUM(NrOfCells)
! OldestActiveCell stores the oldest cell that is still producing rain at a given time. It limits the range of cells that need to be
! checked for rainfall, which will reduce runtime especially when very long rainfall records need to be generated. Its initial value
! is 1 for obvious reasons.
 OldestActiveCell = 1
! Counter I keeps track of the number of intervals on the time line.
 I = 1
 TimelineMaker: DO 
!  Allocate the array for the time line and copy the backup array into it. 
   Allocate (RainfallRate(0:I+1, 1:2))
   RainfallRate(0:I,1:2) = RainfallRateOld
   Deallocate (RainfallRateOld)
!  Find the first occurrence of change: an active cell ending, or a new cell starting to rain. The time from which to start looking
!  forward to find this change is determined by I: RainfallRate(I,1).
   Time = RainfallRate(I,1)
!  Find the starting time of the next cell, and store in the corresponding cell number in NrStart. If the last cell has been
!  activated, the next starting time is set equal to the end of the simulation period. This sequence should only be started when
!  the current time has not yet progressed beyond the starting time of the last cell for which it was determined.
   IF (Time >= NextCellStart) THEN
     IF (NrStart < Help) THEN 
       DO J = NrStart + 1, Help
         IF (AllCellsStartEndDepthOld(J, 1) > Time) THEN
           NextCellStart = AllCellsStartEndDepthOld(J, 1)
           NrStart = J   
           EXIT 
         END IF
       END DO
!      Check if multiple cells start at the same time and update NrStart accordingly. 
       DO J = NrStart + 1, Help
         IF (AllCellsStartEndDepthOld(J, 1) == NextCellStart) THEN
           NrStart = J
         ELSE
           EXIT
         END IF
       END DO
     ELSE
!      In the previous cycle through DO-loop TimelineMaker, NrStart has been increased to the total number of cells. Consequently,
!      no cells will become active anymore, but those that still are can stop producing rainfall before the end of the simulation
!      period. More intervals on the time line may therefore be required. To ensure a proper completion of all intervals before the
!      DO-loop TimelineMaker is terminated, NextCellStart of an imaginary cell is set to the end of the simulation period.
       NextCellStart = REAL(NrOfFourYearPeriods, DP) * 1461.0_DP 
     END IF
   END IF
!
!  Find the cells that are active immediately after the current time, and determine their combined rainfall rate. Also find the 
!  cell that will stop producing rainfall the soonest.
!
!  Initialize the rainfall rate for this period. Also set the initial ending time of the time interval to the end of the simulation
!  period. Initialize the variable that stores the cell that will become inactive the soonest. The initial value of EndingCell
!  cannot be equal to a valid cell number.
   SoonestEnd = REAL(NrOfFourYearPeriods, DP) * 1461.0_DP
   RainfallRate(I + 1,2) = 0.0_DP
   EndingCell = 0
!  The following DO-loop adds the rainfall of all active cells to calculate the total rainfall rate for the interval starting 
!  at Time. It also checks which of those cells will cease producing rainfall the soonest.
   DO J = 1, NrStartOld
     IF ((AllCellsStartEndDepthOld(J,1) <= Time) .AND. (AllCellsStartEndDepthOld(J,2) > Time)) THEN
!      Compute total rainfall in this time interval. Note that these cells will remain active at least until the end of the 
!      time interval that starts at Time.
       RainfallRate(I + 1,2) = RainfallRate(I + 1,2) + AllCellsStartEndDepthOld(J,3)
!      Store the cell that stops being active the soonest. By checking if a cell's ending time is smaller but not equal to the 
!      smallest value found so far, the lowest cell number of all cells ending the earliest (at SoonestEnd) is stored, which is 
!      the oldest of the cells that ends at SoonestEnd. 
!      N.B. During a dry spell, the condition in the above IF statement is evaluated as 'False' for all J, and the following IF 
!      statement will not be executed:
!      SoonestEnd will keep its initial value (end of the simulation period), and EndingCell will equal zero. Thus, a zero-valued
!      EndingCell is an indicator of a dry period.
       IF (AllCellsStartEndDepthOld(J,2) < SoonestEnd) THEN
         SoonestEnd = AllCellsStartEndDepthOld(J,2)
         EndingCell = J     
       END IF
     END IF
   END DO 
!  With the earliest next cell and earliest end of the currently active cells known, determine which is first. Update the 
!  array with the time line accordingly.
!  N.B. Above, the time for the next cell to start was capped at the end of the simulation time. If SoonestEnd is larger than 
!  that, the currently calculated value of the rainfall rate is valid until the end of the simulation period, and the DO loop 
!  TimelineMaker can be exited. Also update the oldest active cell if needed.
   I = I + 1
!  During dry periods, SoonestEnd is set equal to the end of the simulation period. This is also the largest possible value of
!  NextCellStart. During a dry spell, the condition in the next IF statement evaluates as 'False'. Consequently, the first 
!  block in the IF construct is only executed during rainfall.
   FixTimeInterval: IF (SoonestEnd < NextCellStart) THEN
!    The next time interval starts with a cell becoming inactive. Update the oldest active cell if the currently oldest active 
!    cell is the one reaching its end.
     RainfallRate(I,1) = SoonestEnd
     IF ((OldestActiveCell == EndingCell) .AND. (OldestActiveCell < Help)) THEN
       DO J = OldestActiveCell + 1, NrStartOld
!        Look for the oldest active cell that will end after OldestActiveCell becomes inactive. Note that Time is the start of the
!        current interval, and SoonestEnd is its end.
!        OldestActiveCell should be updated to store the oldest cell that produces rain after time SoonestEnd (the start of the next
!        time interval).
         IF ((AllCellsStartEndDepthOld(J,1) <= Time) .AND. (AllCellsStartEndDepthOld(J,2) > SoonestEnd)) THEN 
           OldestActiveCell = J
           EXIT 
         END IF
       END DO
     END IF
   ELSE
!    This section of the IF construct is executed in each of the following cases: 1) for a dry period starting at Time, 2) during a
!    rainfall period, when another cell will start contributing rain at time NextCellStart, or 3) when the end of the simulation
!    period is reached. 
!    For all three cases, the start of the next time period needs to be stored and NrStartOld needs to be updated.
     RainfallRate(I,1) = NextCellStart
     NrStartOld = NrStart
     IF (EndingCell == 0) THEN
!      A dry period starts at Time. It will end at time NextCellStart. OldestActiveCell can be set equal to the oldest cell 
!      becoming active: NrStart or an older cell starting at the same time (NextCellStart). At the start of the next interval, 
!      this cell will indeed become active and thus become the oldest active cell at that time.
       DO J = OldestActiveCell, NrStart
         IF (AllCellsStartEndDepthOld(J,1) == NextCellStart) THEN 
           OldestActiveCell = J
           EXIT 
         END IF
       END DO
     END IF
   END IF FixTimeInterval
!  Check if the end of the simulation period has been reached. If it has, RainfallRateOld needs not be allocated anymore and 
!  RainfallRate is in order.
   IF (RainfallRate(I,1) >= (REAL(NrOfFourYearPeriods, DP) * 1461.0_DP)) THEN
     RainfallRate(I,1) = REAL(NrOfFourYearPeriods, DP) * 1461.0_DP
!    The timeline is complete - leave DO-loop 'TimelineMaker'.
     EXIT
   END IF
!  Update the backup array and copy and then deallocate RainfallRate.
   Allocate (RainfallRateOld(0:I,1:2))
   RainfallRateOld = RainfallRate
   Deallocate (RainfallRate)
 END DO TimelineMaker
! N.B. Counter I stores the total number of time intervals in the rainfall series. This information needs to be written to 
! the output file as it is needed by DupuitFlow (a code for saturated flow). It als needs to be passed on to subroutine
! OutputWriter through Module RainfallTimeline.
 NrOfIntervals = I
 WRITE (22, 400) I
400 FORMAT (/, ' The rainfall time series has ',I10,' intervals.')
! Check mass balance. Store the relative mass balance error in RelativeError. Array RainPerInterval stores the amount of 
! rain [L] cumulated over each interval.
 Allocate (RainPerInterval (1:I))
 RainPerInterval = (RainfallRate(1:I,1) - RainfallRate(0:I-1,1)) * RainfallRate(1:I,2)
 RelativeError = (TotalAmountOfRain - SUM(RainPerInterval)) / TotalAmountOfRain
 WRITE (22, 410) TotalAmountOfRain, RelativeError
410 FORMAT (' The sum of the rainfall per interval in the rainfall record is: ',ES12.4, &
          & ', leading to a relative mass balance error of: ',ES12.4)
 END SUBROUTINE Timeline
!
!23456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012
 END PROGRAM Whollstoptherain
