!23456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012
! TheHeatIsOn generates a temperature record with daily values of the mean, minimum, and maximum temperature and a record of daily
! extraterrestrial radiation and potential evapotranspiration. 
 MODULE DoublePrecision
!  This module defines the precision of Kind DP. 
   IMPLICIT NONE
   INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(20, 300)  
 END MODULE DoublePrecision
!
 MODULE Constants 
!  This module contains the numerical value of constants encountered in the calculations.
   USE DoublePrecision
   IMPLICIT NONE
   REAL (KIND = DP), PARAMETER :: TwoPi = 6.283185307179586476925287_DP, Gsc = 1360.0_DP
!  TwoPi: two times pi
!  Gsc: solar constant (Watts per square meter)
 END MODULE Constants
!
 MODULE TemperatureParameters
!  This module contains the parameters needed to generate the temperature record.
   USE DoublePrecision
   IMPLICIT NONE
   REAL (KIND = DP):: Ac, Ao, AR1var, Meanf, SDa, SDc, SDm, SDo, SDt, Shift, Tavg
!  Ac: amplitude of the annual sinusoidal temperature signal for clear days (deg Celcius)
!  Ao: amplitude of the annual temperature signal for overcast days (deg Celcius)
!  AR1var: the autoregressive variable representing the persistence of the earlier signal (0 <= AR1var < 1) 
!  Meanf: the mean of the normal distribution that is made the exponent of e to arrive at the lognormal distribution of the half
!         the daily temperature range during a day
!  SDa: the standard deviation of the variation of the annual amplitude of the daily mean temperature (deg Celcius)
!  SDc: standard deviation of the normal distribution that is made the exponent of e to arrive at the lognormal distribution of
!       half the daily temperature range during clear days
!  SDm: standard deviation of the normal distribution of the variation around the annual trend of the daily mean temperature
!  SDo: like SDc, for overcast days
!  SDt: standard deviation of the normal distribution of the variation around the mean Tavg of the annual mean temperature 
!       (deg Celcius)
!  Shift: shift of the annual temperature signal (d)
!  Tavg: annual mean temperature (deg Celcius)
 END MODULE TemperatureParameters
!
 MODULE ETParameter
!  This module contains the parameter needed to generate the evapotranspiration record.
   USE DoublePrecision
   IMPLICIT NONE
   REAL (KIND = DP):: Latitude 
!  Latitude: The latitude of the location of interest in radians
 END MODULE ETParameter
!
 MODULE DailyWeather
!  This module contains the arrays with daily weather record (rainfall, temperatures, evapotranspiration, extraterrestrial
!  radiation, and clear/cloudy days. It also contains the length of the weather record.
   USE DoublePrecision
   IMPLICIT NONE
   INTEGER :: NrOfDays, NrOfOvercastDays
!  NrOfDays: length of the temperature and evapotranspiration records
!  NrOfOvercastDays: total number of overcast days in the record
   INTEGER, ALLOCATABLE :: OvercastRecord(:)
!  OvercastRecord: for each day, the entry 0 indicates the day was clear, the entry 1 represents an overcast day
   REAL :: FractionOfOvercastDays
!  FractionOfOvercastDays: the fraction of days that wer cloudy in the entire simulation record.
   REAL (KIND = DP), ALLOCATABLE :: DailyEpotRadiation(:,:), DailyRain(:), DailyTemp(:,:), RainfallSum(:)
!  DailyEpotRadiation: array with daily potential evapotranspiration according to two equations, and the daily extraterrestrial
!                      radiation.
!  DailyRain: array with daily rainfall. the day number is given by the location in the array.
!  DailyTemp: array with daily average, minimum, and maximum temperature.
!  RainfallSum: 30-day sum of the rainfall, centered on the end of the day of interest.
 END MODULE DailyWeather
! 
 MODULE Staircase 
!  This module contains the variable of the staircase function that links the probability of an overcast day to that day’s rainfall.
   USE DoublePrecision 
   IMPLICIT NONE
   REAL (KIND = DP) :: Plow, Phigh, f1, f2, f3 
!  Plow and Phigh: thresholds of the daily amount of rainfall (mm). f1, f2, f3: probability that a day is overcast when the rainfall
!  is smaller than Plow, between Plow and Phigh, and above Phigh, respectively.
 END MODULE Staircase
!
 PROGRAM TheHeatIsOn
!
! TheHeatIsOn generates a temperature record with daily values of the mean, minimum, and maximum temperature. The daily trend is
! composed of a long-term average with annual normally distributed fluctuations and a sinusoidal annual trend around the annual
! average. The amplitude of the annual trend has normally distributed variations around its mean. Superimposed upon this trend
! are first-order autoregressive variations of the detrended daily mean temperature, with the shocks normally distributed with zero
! mean.
! For cloudy days, the amplitude of the sinusoidal trend is reduced. The minimum and maximum temperatures are computed by 
! superimposing lognormally distributed fluctuations. The standard deviation of these is reduced for cloudy days as well.
! The probability that a day is cloudy (overcast) depends on the amount of rainfall of that day through a three-step staircase
! function.
! After the temperature record is completed, TheHeatGoesOn generates a daily potential evapotranspiration record based on the
! original (Hargreaves, 1994) and modified version (Droogers and Allen, 2002) of the Hargreaves equation. The expression for the
! extraterrestrial radiation appearing in both equations is taken form the HYDRUS-1D manual (Simunek et al., 2013, p. 42-43).
! The conversion of the extraterrestrial radiation to equivalent mm of evapotranspiration is carried out using Evett's (2000)
! relationship between the heat of vaporization and temperature (Evett's Eq. 5.45).
!
! Full details are in de Rooij (2018).
!
! Interface: the length of the time period of the temperature record and the daily rainfall are obtained from file
! Daily_Rainfall.OUT. The parameters of the temperature signal and of the staircase function, as well as the seed for the
! random generator are read from file TempEpot.In. This file also provides the latitude (in radians) of the location of
! interest. Monthly and annual rainfall sums are obtained from files Monthly_Rainfall.OUT and Annual_Rainfall.OUT.
! The files with daily, monthly, and annual rainfall are produced by the rainfall generator WhollStopTheRain3.EXE.
!
! The format of input file TempEpot.IN is
!
! Line 1: Will not be read, can be used for labels
! Line 2: Tavg  SDt  Shift  Ao  Ac  SDa  AR1var  SDm  Meanf  SDo  SDc
! Line 3: Will not be read, can be used for labels 
! Line 4: Plow  Phigh  f1  f2  f3
! Line 5: Will not be read, can be used for labels
! Line 6: Latitude
! Line 7: Will not be read, can be used for labels
! Line 8: Seed1
!
! Tavg: long-term average temperature (deg Celcius)
! Shift: time-shift of the sinusoidal annual temperature signal (days)
! Ao: average amplitude of the annual sinusoidal temperature signal for overcast days
! Ac: average amplitude of the annual sinusoidal temperature signal for clear days
! SDa: standard deviation of the normal distribution that governs the interannual variation of Ac (the mean is Tavg)
! AR1var: 1st-order auto-regressive coefficient by which yesterday's detrended temperature is multiplied to provide a memory effect
!         in the day-to-day temperature signal. 0 <= AR1var < 1
! SDm: standard deviation of the normally distributed variations of the fluctuations of the mean daily temperature around the 
!      sinusoidal yearly trend. The mean equals zero.
! Meanf: the mean of the normal distribution of the natural logarithm of the magnitude of the difference between the maximum
!        temperature and the daily mean.
! SDo: the standard deviation of the normal distribution of the natural logarithm of the magnitude of the difference between the
!      maximum temperature and the daily mean. The value is valid for overcast days. 
! SDc: the standard deviation of the normal distribution of the natural logarithm of the magnitude of the difference between the
!      maximum temperature and the daily mean. The value is valid for clear days. 
! Plow: the daily precipitation amount below which the probability of a cloudy day equals f1. 
! Phigh: the daily precipitation amount above which the probability of a cloudy day equals f3. For daily rainfall amounts between
!        Plow and Phigh, the probability of a cloudy day equals f2. 
! f1, f2, f3: see definitions of Plow and Phigh. The value of these parameters must lie between 0 and 1 (inclusive)
! Latitude: geographical latitude of the location of interest in radians. If the latitude is expressed in degrees that value needs
!           to be multiplied by pi and divided by 180. 
! Seed1: a large integer that is used to initiate the random generator
! 
! All input parameters except Seed1 are REAL (double precision). Seed1 is an integer.
! 
! The daily record of temperatures and clear/cloudy days is written to file Temperature_Clouds.OUT. The daily record of 
! potential evapotranspiration and extraterrestrial radiation is written to file Epot.OUT.
! File Monthly_Temp.OUT contains the monthly averages of the daily mean, minimum, and maximum temperatures, and the number of
! cloudy days in each month. File Annual_Temp.OUT has the same data for each year. File Monthly_Epot.OUT stores the monthly sums
! of potential evapotranspiration according to Hargreaves (1994) and Droogers and Allen (2002), as well as the sum of 
! extraterrestrial radiation, all converted to mm water column. The evapotranspiration deficit for both these equations is also
! written to that file. 
! File Annual_Epot.OUT has the same data for each year.
! Annual and monthly averages of key weather characteristics are written to file Weather_Statistics.OUT.
!
! References
! de Rooij, G.H. 2018. A simple weather generator for applications with limited data availability. In preparation.
! Droogers, P., and Allen, R. G.: Estimating reference evapotranspiration under inaccurate data conditions, Irrigation and Drainage
!   Systems, 16, 33-45, 2002.
! Evett, S. R.: Energy and water balances at soil-plant-atmosphere interfaces, In: Sumner. M. E. (ed.) Handbook of soil science,
!   A-129 – A-182, CRC Press, Boca Raton, Fla, U.S.A., 2000.
! Hargreaves, G. H.: Defining and using reference evapotranspiration, J. Irrigation and Drainage Engineering, 120, 1132-1139, 1994.
! Simunek, J. Sejna, M., Saito, H., Sakai, M., and van Genuchten, M. Th.: The HYDRUS-1D software package for simulating the
!   one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, Version 4.17, Dept. Env. Sci.,
!   Univ. Calif., Riverside, California, U.S.A., 2013.
!
! Gerrit H. de Rooij, March 2018
!
 USE DoublePrecision
 USE TemperatureParameters
 USE ETParameter
 USE Staircase
 USE DailyWeather, ONLY: NrOfDays
 IMPLICIT NONE
 INTEGER :: I,Seed1, SeedSize
! I: Counter
! Seed1: the seed that provides Seed(1) of the array that initializes the random generator
! SeedSize: the size of array Seed that initializes the random generator
 INTEGER, ALLOCATABLE :: Seed(:)
! Seed: an array with preferably large integers to initialize the random generator.
!
! Read input.
 OPEN (20, FILE = 'TempEpot.IN', STATUS = 'OLD', ACTION = 'READ')
 READ (20, *)
 READ (20, *) Tavg, SDt, Shift, Ao, Ac, SDa, AR1var, SDm, Meanf, SDo, SDc
 READ (20, *)
 READ (20, *) Plow, Phigh, f1, f2, f3
 READ (20, *)
 READ (20, *) Latitude
 READ (20, *)
 READ (20, *) Seed1
 CLOSE (20)
! Count the number of records in the output file of the rainfall generator with the daily rainfall sums (Daily_Rainfall.OUT).
 OPEN (21, FILE = 'Daily_Rainfall.OUT', STATUS = 'OLD', ACTION = 'READ')
 NrOfDays = 0 
 DO
   READ (21, *, END = 1001) 
   NrOfDays = NrOfDays + 1
 END DO
1001 CLOSE (21)
! Initialize the random generator.
 CALL Random_seed (SIZE = SeedSize)
 Allocate (Seed(1:SeedSize))
 Seed(1) = Seed1
 DO I = 2, SeedSize, 1
   Seed(I) = Seed(1) + (I * 89) ! 89 is the second largest prime number smaller than 100
 END DO
 CALL Random_seed (PUT = Seed)
 Deallocate (Seed)
! Generate the temperature record.
 CALL SummerTime
! Generate the evapotranspiration record.
 CALL WhatGoesUp
! Write the records to output.
 CALL OutputWriter
! Subroutine declarations
 CONTAINS 
!
 SUBROUTINE OutputWriter
!  OutputWriter writes the temperature record, the record of clear and cloudy days, the extraterrestrial radiation record, and the
!  record of potential evapotranspiration to output files. It als generates montnhly and annual records of the average of the daily
!  mean, minimum, and maximum temperature, of the extraterrestrial radiation, of the number of cloudy days, and monthly and annual
!  sums of potential evapotranspiration and extraterrestrial radiation, both expressed in equivalent mm of water column.
!
!  Interface: the daily records (and their length) of temperature, cloud cover, radiation, and potential evapotranspiration are
!  taken from Module DailyWeather. The array with 30-day rainfall sums is also from that Module.
!  The daily record of temperatures and clear/cloudy days is written to file Temperature_Clouds.OUT. The daily record of 
!  potential evapotranspiration and extraterrestrial radiation is written to file Epot.OUT.
!  File Monthly_Temp.OUT contains the monthly averages of the daily mean, minimum, and maximum temperatures, and the number of
!  cloudy days in each month. File Annual_Temp.OUT has the same data for each year. File Monthly_Epot.OUT stores the monthly sums
!  of potential evapotranspiration according to Hargreaves (1994) and Droogers and Allen (2002), as well as the sum of 
!  extraterrestrial radiation, all converted to mm water column. The evapotranspiration deficit for both these equations is also
!  written to that file. File Annual_Epot.OUT has the same data for each year. Negative values of the evapotranspiration deficit
!  indicate that the rainfall is not sufficient to offset potential evapotranspiration.
!  Annual and monthly averages of key weather characteristics are written to file Weather_Statistics.OUT.
! 
!  Gerrit H. de Rooij, March 2018
!
   USE DoublePrecision
   USE DailyWeather, Only: NrOfDays, NrOfOvercastDays, OvercastRecord, FractionOfOvercastDays, DailyEpotRadiation, DailyRain, &
                         & DailyTemp, RainfallSum
!
   IMPLICIT NONE
   INTEGER :: CyclicMonth, CyclicYear, EndDay, FourYearPeriodNr, I, J, NrOfYears, StartDay
!  CyclicMonth, CyclicYear: number of the month/year of interest within a four-year cycle.
!  EndDay: final day of the month or year of interest
!  FourYearPeriodNr: the number of the four-year period in wich the month/year of interest occurs, or the one preceding it
!  I, J: counters
!  NrOfYears: the number of years in the simulation record
!  StartDay: first day of the month or year of interest
   INTEGER, Allocatable :: CloudyDaysPerYear(:), CloudyDaysPerMonth(:)
   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 /) 
!  MonthStartTimes: 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.
   REAL (KIND = DP) :: AnnualMeanCloudyDays, Dummy
   REAL (KIND = DP), Dimension(1:2) :: AnnualMeanEpotDeficit
   REAL (KIND = DP), Dimension(1:3) :: AnnualMeanEpotRa, AnnualMeanTemp
   REAL (KIND = DP), Dimension(1:12) :: MonthlyMeanCloudyDays
   REAL (KIND = DP), Dimension(1:12,1:2) :: MonthlyMeanEpotDeficit 
   REAL (KIND = DP), Dimension(1:12,1:3) :: MonthlyMeanEpotRa, MonthlyMeanTemp
   REAL (KIND = DP), Allocatable :: AnnualEpotDeficit(:,:), AnnualEpotRadiation(:,:), AnnualRain(:), AnnualTemp(:,:), &
                                  & EpotDeficit(:,:), MonthlyEpotDeficit(:,:), MonthlyEpotRadiation(:,:), MonthlyRain(:),&
                                  & MonthlyTemp(:,:)
!  
!  Write the temperature record (including the nr of cloudy days and which days are overcast) to the output file
!  Temperature_Clouds.OUT.
   OPEN (30, FILE = 'Temperature_Clouds.OUT', STATUS = 'UNKNOWN')
   WRITE (30, 100) NrOfOvercastDays, FractionOfOvercastDays
   WRITE (30, 110) 
   DO I = 1, NrOfDays
     WRITE (30, 120) I, DailyTemp(I,1), DailyTemp(I,2), DailyTemp(I,3), OvercastRecord(I)
   END DO
100 FORMAT (' Number of days with cloud cover: ',I7,'; Fraction of the total time period: ',F6.3)
110 FORMAT (' Day number  Average temp  Minimum temp  Maximum temp  Overcast(1)/Clear(0)')
120 FORMAT (1X, I11, 3(2X, ES12.4), 10X, I1)
   CLOSE (30)
!  Write the record of potential evapotranspiration and solar radiation to output file Epot.OUT
   OPEN (31, FILE = 'Epot.OUT', STATUS = 'UNKNOWN')
   WRITE (31, 130) 
   DO I = 1, NrOfDays
     WRITE (31, 140) I, DailyEpotRadiation(I,1), DailyEpotRadiation(I,2), DailyEpotRadiation(I,3), RainfallSum(I)
   END DO
   CLOSE (31)
130 FORMAT (' Day number      Epot Hargreaves (1994)  Epot Droogers and Allen (2002)  Extraterrestrial radiation  30-day rain sum'/&
          & '                 mm per day              mm per day                      equivalent mm per day       mm')
140 FORMAT (1X, I9, 2X, ES15.4, 9X, ES15.4, 17X, ES15.4, 13X, ES15.4)
!
!  Create files with daily, monthly, and annual sums of potential evapotranspiration and averages of the daily
!  mean and extreme temperatures.
   OPEN(32, FILE = 'Monthly_Temp.OUT', STATUS = 'UNKNOWN')
   OPEN(33, FILE = 'Annual_Temp.OUT', STATUS = 'UNKNOWN')
   OPEN(34, FILE = 'Monthly_Epot.OUT', STATUS = 'UNKNOWN')
   OPEN(35, FILE = 'Annual_Epot.OUT', STATUS = 'UNKNOWN')
!  Run through all days of the simulated period (1461 days in each four-year period). 
   NrOfYears = 4 * NrOfDays / 1461
   Allocate (AnnualTemp(1:NrOfYears,1:3))
   Allocate (MonthlyTemp(1: 12 * NrOfYears, 1:3))
   Allocate (AnnualEpotRadiation(1:NrOfYears,1:3))
   Allocate (MonthlyEpotRadiation(1: 12 * NrOfYears, 1:3))
   Allocate (CloudyDaysPerYear(1:NrOfYears))
   Allocate (CloudyDaysPerMonth(1: 12 * NrOfYears))
   Allocate (MonthlyRain(1: 12 * NrOfYears))
   Allocate (AnnualRain(1:NrOfYears))
   Allocate (EpotDeficit(1:NrOfDays, 1:2))
   Allocate (AnnualEpotDeficit(1:NrOfYears, 1:2))
   Allocate (MonthlyEpotDeficit(1: 12 * NrOfYears, 1:2))
!  Calculate the daily evapotranspiration deficit (Rainfall - Epot) for both equations for the potential evapotranspiration. 
   DO I = 1, NrOfDays
     EpotDeficit(I,1) = DailyRain(I) - DailyEpotRadiation(I,1)
     EpotDeficit(I,2) = DailyRain(I) - DailyEpotRadiation(I,2)
   END DO 
   WRITE(32, 150)
   WRITE(34, 160)
!  Read the monthly rainfall record.
   OPEN(36, FILE = 'Monthly_Rainfall.OUT', STATUS = 'OLD', ACTION = 'READ')
   DO I = 1, 12 * NrOfYears
     READ (36, *) Dummy, MonthlyRain(I)
   END DO
   CLOSE (36)
!  Read the annual rainfall record. 
   OPEN(37, FILE = 'Annual_Rainfall.OUT', STATUS = 'OLD', ACTION = 'READ')
   DO I = 1, NrOfYears
     READ (37, *) Dummy, AnnualRain(I)
   END DO
   CLOSE (37)
!  Determine the monthly sums and averages 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)
!    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) + 1431
       EndDay = ((FourYearPeriodNr - 1) * 1461) + 1461
     END IF
     CloudyDaysPerMonth(I) = SUM(OvercastRecord(StartDay:EndDay))
     DO J = 1, 3
       MonthlyTemp(I,J) = SUM(DailyTemp(StartDay:EndDay, J)) / REAL(EndDay - StartDay + 1, DP)
       MonthlyEpotRadiation(I,J) = SUM(DailyEpotRadiation(StartDay:EndDay, J))
     END DO
     DO J = 1, 2
       MonthlyEpotDeficit(I,J) = SUM(EpotDeficit(StartDay:EndDay, J))
     END DO
     WRITE(32, 170) I, MonthlyTemp(I,1), MonthlyTemp(I,2), MonthlyTemp(I,3), CloudyDaysPerMonth(I)
     WRITE(34, 180) I, MonthlyRain(I), MonthlyEpotRadiation(I,1), MonthlyEpotRadiation(I,2), MonthlyEpotRadiation(I,3), &
                  & MonthlyEpotDeficit(I,1), MonthlyEpotDeficit(I,2)
   END DO
150 FORMAT (' Month number  Average temp  Avg min temp  Avg max temp  Nr of overcast days')
160 FORMAT (' Month number     Rainfall  Epot Hargreaves  Epot Droogers &  Extraterr. rad.  ET Deficit Har.   ET Deficit D&A'/&
          & '                        mm        (1994) mm  Allen (2002) mm    equivalent mm       (1994)  mm        (2002) mm') 
170 FORMAT (2X, I11, 3(2X, ES12.4), 18X, I3)
180 FORMAT (1X, I11, 2X, ES12.4, 5(2X, ES15.4))
!
   WRITE(33, 190)
   WRITE(35, 200)
!  Determine the annual sums and averages and write them to output.
   DO I = 1, NrOfYears
     CyclicYear = MOD(I,4)
     IF (CyclicYear == 0) CyclicYear = 4
!    Determine in which four-year period year nr I occurs.
     FourYearPeriodNr = INT(REAL(I) / 4.0)
!    Determine the first and last day of the year of interest.
     IF (CyclicYear < 4) THEN
       StartDay = (FourYearPeriodNr * 1461) + ((CyclicYear - 1) * 365) + 1    
       EndDay = StartDay + 364
     ELSE
       StartDay = ((FourYearPeriodNr - 1) * 1461) + (3 * 365) + 1 
      EndDay = StartDay + 365
     END IF
     CloudyDaysPerYear(I) = SUM(OvercastRecord(StartDay:EndDay))
     DO J = 1, 3
       AnnualTemp(I, J) = SUM(DailyTemp(StartDay:EndDay, J)) / REAL(EndDay - StartDay + 1, DP)
       AnnualEpotRadiation(I, J) = SUM(DailyEpotRadiation(StartDay:EndDay, J))
     END DO
     DO J = 1, 2
       AnnualEpotDeficit(I,J) = SUM(EpotDeficit(StartDay:EndDay, J))
     END DO
     WRITE(33, 170) I, AnnualTemp(I,1), AnnualTemp(I,2), AnnualTemp(I,3), CloudyDaysPerYear(I)
     WRITE(35, 180) I, AnnualRain(I), AnnualEpotRadiation(I,1), AnnualEpotRadiation(I,2), AnnualEpotRadiation(I,3), &
                  & AnnualEpotDeficit(I,1), AnnualEpotDeficit(I,2) 
   END DO
190 FORMAT ('  Year number  Average temp  Avg min temp  Avg max temp  Nr of overcast days') 
200 FORMAT ('  Year number     Rainfall  Epot Hargreaves  Epot Droogers &  Extraterr. rad.  ET Deficit Har.   ET Deficit D&A'/&
          & '                        mm        (1994) mm  Allen (2002) mm    equivalent mm       (1994)  mm        (2002) mm') 
!
   CLOSE (32)
   CLOSE (33)
   CLOSE (34)
   CLOSE (35)
!
!  Create a climate summary with the average values of key statistics for every month.
!  Compute the mean of the monthly and annual key values.
   MonthlyMeanTemp = 0.0_DP
   AnnualMeanTemp = 0.0_DP
   MonthlyMeanCloudyDays = 0.0_DP
   AnnualMeanCloudyDays = 0.0_DP
   MonthlyMeanEpotRa = 0.0_DP
   AnnualMeanEpotRa = 0.0_DP
   MonthlyMeanEpotDeficit = 0.0_DP
   AnnualMeanEpotDeficit = 0.0_DP
   DO I = 0, NrOfYears - 1
     AnnualMeanTemp = AnnualMeanTemp + AnnualTemp(I + 1, 1:3)
     AnnualMeanEpotRa = AnnualMeanEpotRa + AnnualEpotRadiation(I + 1, 1:3)
     AnnualMeanEpotDeficit = AnnualMeanEpotDeficit + AnnualEpotDeficit(I + 1, 1:2)
     AnnualMeanCloudyDays = AnnualMeanCloudyDays + REAL(CloudyDaysPerYear(I + 1), DP)
     DO J = 1, 12
       MonthlyMeanTemp(J, 1:3) = MonthlyMeanTemp(J, 1:3) + MonthlyTemp((12 * I) + J, 1:3)
       MonthlyMeanEpotRa(J, 1:3) = MonthlyMeanEpotRa(J, 1:3) + MonthlyEpotRadiation((12 * I) + J, 1:3)
       MonthlyMeanEpotDeficit(J, 1:2) = MonthlyMeanEpotDeficit(J, 1:2) + MonthlyEpotDeficit((12 * I) + J, 1:2)
       MonthlyMeanCloudyDays(J) = MonthlyMeanCloudyDays(J) + REAL(CloudyDaysPerMonth((12 * I) + J), DP)
     END DO
   END DO
   AnnualMeanTemp = AnnualMeanTemp / REAL(NrOfYears, DP)
   AnnualMeanEpotRa = AnnualMeanEpotRa / REAL(NrOfYears, DP) 
   AnnualMeanEpotDeficit = AnnualMeanEpotDeficit / REAL(NrOfYears, DP) 
   AnnualMeanCloudyDays = AnnualMeanCloudyDays / REAL(NrOfYears, DP)
   MonthlyMeanTemp = MonthlyMeanTemp / REAL(NrOfYears, DP)
   MonthlyMeanEpotRa = MonthlyMeanEpotRa / REAL(NrOfYears, DP) 
   MonthlyMeanEpotDeficit = MonthlyMeanEpotDeficit / REAL(NrOfYears, DP) 
   MonthlyMeanCloudyDays = MonthlyMeanCloudyDays / REAL(NrOfYears, DP)
!  Write the descriptive statistics to output.
   OPEN (38, FILE = 'Weather_Statistics.OUT', STATUS = 'UNKNOWN')
   WRITE(38, 210) AnnualMeanTemp(1), AnnualMeanTemp(2), AnnualMeanTemp(3), AnnualMeanCloudyDays, AnnualMeanEpotRa(1), &
                & AnnualMeanEpotRa(2), AnnualMeanEpotRa(3), AnnualMeanEpotDeficit(1), AnnualMeanEpotDeficit(2)
210 FORMAT(' Annual mean values of',/,&
          &'  Temperature    Daily min.    Daily max.   Cloudy days   Epot Hargr.  Epot Dro&All  Exttrs. rad.    ET Deficit&
          &   ET Deficit D&A',/,&
          &'  deg Celcius   temp. deg C   temp. deg C                   (1994) mm     (2002) mm        eq. mm   H (1994) mm&
          &    D&A (2002) mm',/, 1X,ES12.4,7(2X,ES12.4),5X,ES12.4)
   WRITE(38, 220)
220 FORMAT(/,' Monthly mean values of',/,&
        &'            Temperature    Daily min.    Daily max.   Cloudy days   Epot Hargr.  Epot Dro&All  Exttrs. rad.    ET Deficit&
        &     ET Deficit',/,&
        &' Month nr   deg Celcius   temp. deg C   temp. deg C                   (1994) mm     (2002) mm        eq. mm   H (1994) mm&
        &  D&A (2002) mm',/)
   DO I = 1, 12
     WRITE(38, 230) I, (MonthlyMeanTemp(I,J), J=1,3), MonthlyMeanCloudyDays(I), (MonthlyMeanEpotRa(I,J), J=1,3),&
                 & (MonthlyMeanEpotDeficit(I,J), J=1,2)
   END DO 
230 FORMAT(6X,I3,8(2X,ES12.4),3X,ES12.4)
   CLOSE (38)
   Deallocate (DailyTemp)
   Deallocate (OvercastRecord)
   Deallocate (RainfallSum)
   Deallocate (DailyEpotRadiation)
   Deallocate (AnnualTemp)
   Deallocate (MonthlyTemp)
   Deallocate (AnnualEpotRadiation)
   Deallocate (MonthlyEpotRadiation)
   Deallocate (CloudyDaysPerYear)
   Deallocate (CloudyDaysPerMonth)
   Deallocate (DailyRain)
   Deallocate (MonthlyRain)
   Deallocate (AnnualRain)
   Deallocate (EpotDeficit)
   Deallocate (AnnualEpotDeficit)
   Deallocate (MonthlyEpotDeficit) 
! 
 END SUBROUTINE OutputWriter
!
 SUBROUTINE SummerTime
!  SummerTime generates a temperature record with daily values of the mean, minimum, and maximum temperature. The daily trend is
!  composed of a long-term average and a sinusoidal annual trend around the average. The amplitude of the annual trend has
!  normally distributed variations around its mean. Superimposed upon this trend are first-order autoregressive variations of the 
!  detrended daily mean temperature, with the shocks normally distributed with zero mean.
!  For cloudy days, the amplitude of the sinusoidal trend is reduced. The minimum and maximum temperatures are computed by 
!  superimposing lognormally distributed fluctuations. The standard deviation of these is reduced for cloudy days as well.
!  The probability that a day is cloudy (overcast) depends on the amount of rainfall of that day through a three-step staircase
!  function.
!
!  Full details are in de Rooij (2018).
!
!  Interface: the length of the time period of the temperature record and the daily rainfall are obtained from file
!  Daily_Rainfall.OUT, one of the output files of the rainfall generator WhollStopTheRain3. The parameters of the temperature signal
!  and of the staircase function, as well as the seed for the random generator are read from file Temperature.In. This file also
!  provides the latitude (in radians) of the location of interest. The temperature and overcast record is written to file
!  Temperature_Clouds.OUT. The daily temperature and rainfall records are stored in Module DailyWeather.
!
!  Reference
!  de Rooij, G.H. 2018. A simple weather generator for applications with limited data availability. In preparation.
!
!  Gerrit H. de Rooij, March 2018
!
   USE DoublePrecision
   USE TemperatureParameters
   USE DailyWeather, ONLY: NrOfDays, NrOfOvercastDays, FractionOfOvercastDays, OvercastRecord, DailyRain, DailyTemp
   USE Constants, ONLY: TwoPi
   USE Staircase 
   IMPLICIT NONE
   INTEGER :: I
!  I: counter
   REAL (KIND = DP):: DetrendedTempYesterday, Dummy, HalfRange, SDf, TempTrend
!  An: amplitude of the annual sinusoidal temperature signal (set to Ac oder Ao) 
!  Dummy: Dummy variable that helps skip a column of an input file that is not required.
!  HalfRange: maximum temperature of day I minus average temperature of day I
!  SDf: standard deviation of the normal distribution that is made the exponent of e to arrive at the lognormal distribution of
!       half the daily temperature range (set to SDc or SDo)
!  TempTrend: the mean daily temperature for a given day according to the annual trend (deg Celcius)
   REAL (KIND = DP), ALLOCATABLE :: Amplitude(:), AnnualAverageT(:)
!  Amplitude: array with the amplitude of the annual trend for every year, given day by day taking into account if the day is 
!             cloudy or clear.
!  Tavg: array with the annual values of the long-term average temperature. The value changes for each year.
   Allocate (DailyRain (1:NrOfDays))
   Allocate (DailyTemp (1:NrOfDays, 1:3))
   Allocate (OvercastRecord (1:NrOfDays))
   Allocate (AnnualAverageT (1:NrOfDays))
   Allocate (Amplitude (1:NrOfDays))
!  Read the rainfall data.
   OPEN (40, FILE = 'Daily_Rainfall.OUT', STATUS = 'OLD', ACTION = 'READ')
   DO I = 1, NrOfDays
     READ (40, *) Dummy, DailyRain(I)
   END DO
   CLOSE (40)
!  Determine a time line of years and annual amplitudes. Initialize array OvercastRecord.
   OvercastRecord = 0
   CALL WhatYearIsIt(AnnualAverageT, Amplitude, OvercastRecord)
   NrOfOvercastDays = SUM(OvercastRecord)
!  Generate the temperature record.
   DetrendedTempYesterday = NormalVariate(0.0_DP, SDm)
   DO I = 1, NrOfDays
!    Set the standard deviation of the fluctuations of the daily extremes depending on whether the day is overcast or clear.
     IF (OvercastRecord(I) == 0) THEN
       SDf = SDc
     ELSE
       SDf = SDo
     END IF
!    Calculate the various temperatures.
     TempTrend = AnnualAverageT(I) + (Amplitude(I) * sin(TwoPi * (Shift + REAL(JulianDay(I), DP) - 0.5_DP) / 365.0_DP))
     DailyTemp(I,1) = TempTrend + (AR1var * DetrendedTempYesterday) + NormalVariate(0.0_DP, SDm) 
     DetrendedTempYesterday = DailyTemp(I,1) - TempTrend 
     HalfRange = EXP(NormalVariate(Meanf, SDf))
     DailyTemp(I,2) = DailyTemp(I,1) - HalfRange
     DailyTemp(I,3) = DailyTemp(I,1) + HalfRange
   END DO
   FractionOfOvercastDays = REAL(NrOfOvercastDays) / REAL(NrOfDays)
!  Deallocate the allocatable arrays that are not needed outside the routine.
   Deallocate (AnnualAverageT)
   Deallocate (Amplitude)
!
 END SUBROUTINE SummerTime
!
 SUBROUTINE WhatGoesUp
!  WhatGoesUp generates a daily potential evapotranspiration record based on the original (Hargreaves, 1994) and modified version
!  (Droogers and Allen, 2002) of the Hargreaves equation. The expression for the extraterrestrial radiation appearing in both 
!  equations is taken form the HYDRUS-1D manual (Simunek et al., 2013, p. 42-43). The conversion of the extraterrestrial radiation
!  to equivalent mm of evapotranspiration is carried out using Evett's (2000) relationship between the heat of vaporization and 
!  temperature (his Eq. 5.45).
!
!  Interface: the length of the time period (NrOfDays) of the temperature record and the daily rainfall as well as these records
!  themselves (DailyRain and DailyTemp) are obtained from Module DailyWeather.
!  The latitude of the location of interest (rad), needed to compute the extraterrestrial radiation, is taken from Module
!  ETParameter.
!  The evapotranspiration record is written to output file Epot.OUT.
!
!  References
!  Droogers, P., and Allen, R. G.: Estimating reference evapotranspiration under inaccurate data conditions, Irrigation and 
!    Drainage Systems, 16, 33-45, 2002.
!  Evett, S. R.: Energy and water balances at soil-plant-atmosphere interfaces, In: Sumner. M. E. (ed.) Handbook of soil science,
!    A-129 – A-182, CRC Press, Boca Raton, Fla, U.S.A., 2000.
!  Hargreaves, G. H.: Defining and using reference evapotranspiration, J. Irrigation and Drainage Engineering, 120, 1132-1139, 1994.
!  Simunek, J., Sejna, M., Saito, H., Sakai, M., and van Genuchten, M. Th.: The HYDRUS-1D software package for simulating the
!    one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, Version 4.17, Dept. Env. Sci.,
!    Univ. Calif., Riverside, California, U.S.A., 2013.
! 
!  Gerrit de Rooij, March 2018
!
   USE DoublePrecision
   USE DailyWeather, ONLY: NrOfDays, DailyEpotRadiation, DailyRain, DailyTemp, RainfallSum 
   USE Constants
   USE ETParameter
!
   IMPLICIT NONE
!
   INTEGER :: I
!    I: counter
   REAL (KIND = DP) :: CosLat, GonioArgPart, GonioArgFull, RaLead, RelativeSunDistance, SinLat, &
                     & SolarDeclination, SunsetHourAngle, TanLat, TemperatureRange
!  CosLat: cosine of the latitude of the location
!  GonioArgPart: constant part of the argument of goniometric functions involving the Julian day number
!  GonioArgFull: full argument of goniometric functions that involves the Julian day number 
!  RaLead: leading term of the expression for the extraterrestrial radiation (MJ per square meter per day)
!  RelativeSunDistance: relative distance of the Earth to the Sun (mean distance = 1)
!  SinLat: sine of the latitude of the location
!  TanLat: tangent of the latitude of the location
!  TemperatureRange: the difference between maximum and minimum temperature during the day of interest
!
   Allocate (RainfallSum (1:NrOfDays))
   Allocate (DailyEpotRadiation (1:NrOfDays, 1:3))
!
   OPEN (50, FILE = 'WhatgoesupDiagnostics.OUT', STATUS = 'UNKNOWN')
!
!  Compute the 30-day rainfall sum for the first 14 days. Pad up the missing days at the start of the record with those at the end.
!  N.B. When data are read from a file, round-off errors can cause minute negative rainfall sums (< 10E-30). Make sure these are 
!  avoided and write a warning if significant negative values occur.
   DO I = 1, 14
     RainfallSum(I) = SUM(DailyRain(NrOfDays-14+I:NrOfDays)) + SUM(DailyRain(1:I+15))
     IF (RainfallSum(I) < 0.0000000000000000000_DP) THEN
       IF (RainfallSum(I) < -0.00000099999999999999_DP) WRITE(50, 300) I, RainfallSum(I)
       RainfallSum(I) = 0.0000000000000000000_DP
     END IF
   END DO
!  Calculate the 30-day rainfall sum for the final 15 days. Take the missing days at the end of the record from the beginning.
   DO I = 1, 15
     RainfallSum(NrOfDays-15+I)= SUM(DailyRain(1:I)) + SUM(DailyRain(NrOfDays-29+I:NrOfDays))
     IF (RainfallSum(I) < 0.0000000000000000000_DP) THEN
       IF (RainfallSum(I) < -0.00000099999999999999_DP) WRITE(50, 300) I, RainfallSum(I)
       RainfallSum(I) = 0.0000000000000000000_DP
     END IF
   END DO
!  Calculate the intermediate 30-day sums.
   DO I = 15, NrOfDays - 15
     RainfallSum(I) = RainfallSum(I-1) - DailyRain(I-15) + DailyRain(I+15)
     IF (RainfallSum(I) < 0.0000000000000000000_DP) THEN
       IF (RainfallSum(I) < -0.00000099999999999999_DP) WRITE(50, 300) I, RainfallSum(I)
       RainfallSum(I) = 0.0000000000000000000_DP
     END IF
   END DO
300 FORMAT(' The 30-day rainfall sum of day ',I9,' equals ',ES15.4,'.')
!  Calculate various terms in the equations that are independent of time. Start with the leading term of the extraterrestrial
!  radiation.
   RaLead = 0.1728_DP * Gsc / TwoPi
!  Terms that involve the latitude of the location of interest.
   SinLat = SIN(Latitude)
   CosLat = COS(Latitude)
   TanLat = TAN(Latitude)
!  The constant part of the function argument involving the Julian day number.
   GonioArgPart = TwoPi / 365.0_DP
!  Generate the potential evapotranspiration record
   DO I = 1, NrOfDays
     ! Calculate the extraterrestrial radiation in MJ per square meter per day. Start with the solar declination,
     ! then the sunset hour angle, then the relative distance to the Sun.
     GonioArgFull = GonioArgPart * JulianDay(I)
     SolarDeclination = 0.409_DP * SIN(GonioArgFull - 1.39_DP)
     SunsetHourAngle = ACOS(-TanLat * TAN(SolarDeclination))
     RelativeSunDistance = 1.0_DP + (0.033_DP * COS(GonioArgFull))
     DailyEpotRadiation(I,3) = (SunsetHourAngle*SinLat*SIN(SolarDeclination)) + (SIN(SunsetHourAngle)*CosLat*COS(SolarDeclination))
     ! The expression below gives the extraterrestrial radiation in MegaJoules per square meter per day.
     DailyEpotRadiation(I,3)  = DailyEpotRadiation(I,3)  * RelativeSunDistance * RaLead
     ! Convert the extraterrestrial radiation to equivalent mm of evapotranspiration depending on the mean temperature of the day,
     ! based on Eq. 5.45 of Evett (2000).
     DailyEpotRadiation(I,3)  = DailyEpotRadiation(I,3) / (2.501_DP - (0.002370_DP * DailyTemp(I,1)))
     ! Calculate the potential evaporation according to the original Hargreaves equation.
     TemperatureRange = DailyTemp(I, 3) - DailyTemp(I, 2)
     ! Avoid negative values for very low average temperatures DailyTemp(I,1).
     DailyEpotRadiation(I,1) = 0.0023_DP * DailyEpotRadiation(I,3) * MAX((DailyTemp(I,1) + 17.8_DP), 0.0000000000000000000_DP) * &
                             & SQRT(TemperatureRange)
     ! Calculate the potential evaporation according to the modified Hargreaves equation. Avoid negative values during very
     ! wet periods.
     DailyEpotRadiation(I,2) = (MAX((TemperatureRange - (0.0123_DP * RainfallSum(I))),0.0000000000000000000_DP))**(0.76_DP)
     ! Avoid negative values for very low average temperatures DailyTemp(I, 1).
     DailyEpotRadiation(I,2) = 0.0013_DP * DailyEpotRadiation(I,3) * MAX((DailyTemp(I,1) + 17.0_DP), 0.0000000000000000000_DP) * &
                             & DailyEpotRadiation(I,2)
   END DO
   CLOSE (50)
!
 END SUBROUTINE WhatGoesUp   
!
 SUBROUTINE WhatYearIsIt(AnnualAverageT, Amplitude, OvercastRecord)
!  WhatYearIsIt determines the year of the weather record for every day of that record, and determines the average temperature 
!  for that year. It also determines whether a day is cloudy or clear. With this information it determines the value of the
!  amplitude of the annual signal for each day. The amplitude of the sinuoidal trend varies from year to year with normally
!  distributed variations around its annual mean, which itself is normally distributed around its long-term mean. For overcast
!  days the amplitudenis reduced. For every year a random fluctuation is determined for the annual mean temperature and the
!  amplitude of the annual fluctuation around that mean and applied. The latter is pllied to the amplitude values of clear and
!  overcast days. 
!
!  Interface: the number of days (NrOfDays) and the daily rain (DailyRain) are obtained from Module DailyWeather. The long-term
!  mean temperature (Tavg) and the standard deviation of the yearly fluctuations (SDt) are obtained from Module
!  TemperatureParameters. The amplitude of the annual signal for clear (Ac) and overcast (Ao) days, and the standard deviation
!  around these means for clear days (SDa) are also received through from Module TemperatureParameters. 
!  The realizations of the cloudiness (0: clear, 1: cloudy) are returned through array OvercastRecord. The daily values of the
!  annual mean temperature and the amplitude are returned through arrays AnnualAverageT and Amplitude.
! 
!  Gerrit H. de Rooij, March 2018
!
   USE DoublePrecision
   USE TemperatureParameters, ONLY: Ac, Ao, Tavg, SDa, SDt
   USE DailyWeather, ONLY: NrOfDays, DailyRain
!
   IMPLICIT NONE
!
   INTEGER :: I, QuatDay
   INTEGER, Dimension(:) :: OvercastRecord
   REAL (KIND = DP) :: OneAmplitude, OneTavg, U
   REAL (KIND = DP), Dimension(:) :: Amplitude, AnnualAverageT
!
   DO I = 1, NrOfDays
!    Determine which day number it is within a four-year cycle.
     QuatDay = MODULO(I, 1461)
     IF (QuatDay == 0) QuatDay = 1461
!    Reset the annual mean temperature and the amplitude of the sinusoidal temperature signal for clear days on January 1st.
     IF ((QuatDay == 1) .OR. (QuatDay == 366) .OR. (QuatDay == 731) .OR. (QuatDay == 1096)) THEN
       OneTavg = Normalvariate(Tavg, SDt)
       OneAmplitude = NormalVariate(Ac, SDa)
     END IF
     AnnualAverageT(I) = OneTavg
!    Determine if the day is cloudy.
     CALL Random_Number(U)
     IF (U < ProbOvercast(DailyRain(I))) THEN 
!      Change OvercastRecord(I) from zero to one and set the amplitude to its overcast value.
       OvercastRecord(I) = 1
       Amplitude(I) = OneAmplitude * (Ao / Ac)
     ELSE ! Set the amplitude to clear-sky values.
       Amplitude(I) = OneAmplitude
     END IF   
   END DO 
! 
 END SUBROUTINE WhatYearIsIt
!
 FUNCTION JulianDay(I)
!  Function JulianDay determines the Julian day number (from 1 on January 1st to 365 on December 31st) of day nr I of the weather  
!  record. The function assumes that day 1 of the record is January first of a year following a leap year. February 29 receives
!  the same Julian day number as February 28.
! 
!  Interface: The day number of interest is received through the argument list. The Julian day number is returned through the
!  function value. 
!
!  Gerrit de Rooij, January 2018
!
   IMPLICIT NONE
   INTEGER :: JulianDay, QuatDay
!  JulianDay: Julian day number (from 1 of January 1st to 365 on December 31st)
!  QuatDay: similar to JulianDay, but for a four-year period of which the last year is a leap year (running from 1 to 1461)
   INTEGER, INTENT(IN) :: I
!  I: day number of interest in the time record
!
!  Determine the Julian day number. Start by determining if we are before, on, or after February 29.
   QuatDay = MODULO(I, 1461) ! The day number within a four-year period.
   IF (QuatDay == 0) QuatDay = 1461
   IF (QuatDay < 1551) THEN ! Before February 29.
     JulianDay = MODULO(QuatDay, 365)
   ELSE IF (QuatDay == 1551) THEN ! February 29 receives the same Julian day number as February 28.
     JulianDay = 59
   ELSE ! The remaining leap year days receive Julian day numbers as if they were regular days.
     JulianDay = MODULO(QuatDay - 1, 365)
   END IF
   IF (JulianDay == 0) JulianDay = 365
 END FUNCTION JulianDay
!
 FUNCTION NormalVariate(Mean, Standarddeviation)
!  Function NormalVariate returns a random variate with Normal distribution with the mean and standard deviation provided in
!  the argument list.
!  The expression reguires two variates with a standard uniform distribution (see Press et al., 1992, p. 279)
!
!  Interface: the mean and standard deviation are received through the argument list. The random value is returned through the
!  function value.
!  
!  Reference:
!  Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical Recipes in FORTRAN. The art of scientific
!  computing, second ed., Cambridge university Press, Cambridge, United Kingdom, 963 p, 1992.
!
!  Gerrit de Rooij, February 2018
!
   USE DoublePrecision
   USE Constants, ONLY: TwoPi
   IMPLICIT NONE
   REAL (KIND = DP), INTENT(IN) :: Mean, StandardDeviation 
   REAL (KIND = DP) :: NormalVariate, U1, U2
!  NormalVariate: a random variable with a normal distribution
!  U1, U2: random variables with a standard uniform distribution
!
   IF (Standarddeviation <= 0.0_DP) THEN
     NormalVariate = Mean
   ELSE
     CALL Random_Number(U1)
     CALL Random_Number(U2)
     NormalVariate = SQRT(-2.0_DP * LOG(U1)) * COS(TwoPi * U2)
     NormalVariate = Mean + (StandardDeviation * NormalVariate)
   END IF
 END FUNCTION NormalVariate
!
 FUNCTION ProbOvercast(AmountOfRain) 
!  Function ProbOvercast determines the probability that a day is overcast based on the amount of rainfall during that day. 
!  It adopts a three-step staircase function.
!
!  Interface: The parameters of the staircase function are taken from Module Staircase. The daily rainfall is received through
!  argument AmountOfRain in the heading. The probability that the day is overcast is returned to the calling program unit
!  through the function value.
!
!  Gerrit de Rooij, January 2018
!
   USE DoublePrecision
   USE Staircase
   IMPLICIT NONE
   REAL (KIND = DP) :: ProbOvercast
   REAL (KIND = DP), INTENT(IN) :: AmountOfRain
   IF (AmountOfRain < Plow) THEN
     ProbOvercast = f1
   ELSE IF (AmountOfRain < Phigh) THEN
     ProbOvercast = f2
   ELSE
     ProbOvercast = f3
   END IF
 END FUNCTION ProbOvercast
!
 END PROGRAM TheHeatIsOn
