!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