Load Parameters (You could use DoE)

tk=40000 #tank size in gallons
rf=3000  #roof size in square feet
its=1000 #number of its to run simulation
yrs=30 #number of yrs to run simulation
mo=12 #number of mo each year
t=yrs*mo #short hand for total mo
prime=10000 #initial gallons primed in water tank
constant=7.48 #gallons for each cubic feet of rain

Rainfall Data

rainfall=matrix(c(0.091, 0.328,0.033,0.342,0.221,0.2,0.053,0.396,0,0.023,0.342,0.175,
                  0.106, 0.084,0.033,0.221,0.148,0.231,0.039,0.007,0,0.038,0.376,0.026,
                  0.035, 0.232,0.128,0.197,0.359,0.321,0.097,0.396,0.124,0,0.243,0.032,
                  0.142, 0.341,0.203,0.009,0.715,0.005,0.078,0.238,0,0.069,0.231,0.018,
                  0.118, 0.386,0.108,0.013,0.218,0.304,0.192,0.654,0.038,0.136,0.084,0.161,
                  0.377, 0.069,0,0.346,0.857,0.056,0.055,0.646,0,0.014,0.07,0.036,
                  0.077, 0,0.649,0.701,0.214,0.304,0.075,1.342,0.228,0.116,0.234,0.097,
                  0, 0.52,0.09,0.079,0.393,0.163,0,0.703,0,0.128,0.013,0.01,
                  0.158, 0.323,0.573,0.662,0.12,0.074,0.443,0.182,0,0.413,0.343,0.139,
                  0.605, 0.296,0.593,0.182,0.412,0.159,0.149,0,0,0.692,0.021,0.253,
                  0, 0.328,0.198,0.023,0.625,0.082,0.022,0.031,0,0.063,0.002,0.061,
                  0.092, 0.223,0.312,0.012,0.006,0,0.217,0.018,0.003,0.16,0.083,0.137),
                byrow=TRUE, c(12,12))

Storage of Results / Seed Initialization

h20=matrix(rep(0,yrs*mo*its),nrow=yrs*mo)#initialize storage vector
set.seed(1234)

Simulate Demand & Efficiency

demand=runif(t*its, 4000,5200) #demand in gallons
dim(demand)=c(t,its)
efficiency=runif(t*its,900,980)/1000 #estimated efficiency
dim(efficiency)=c(t,its)

Rain Sample

myf=function(x) sample(rainfall[x, 1:12],its, replace=T)
rain=matrix(rep(0,t*its), nrow=t)

k=0
for (i in 1:yrs){
  for (j in 1:mo){
    k=k+1
    rain[k,]=myf(j) 
  }}

Calculations

capture=rf*rain*efficiency*constant #all rainfall captured
out_of_water=matrix(c(rep(0,yrs*mo*its)),nrow=yrs*mo) #tracks amount of water needed but not on-hand
overflow=matrix(c(rep(0,yrs*mo*its)),nrow=yrs*mo) #tranks the overflow

Initial Calculation Matrix (no overflow or out-of-water indicators)

temp=matrix(rep(0,t*its), nrow=t)
temp[1,]=prime+capture[1,1:its]-demand[1,1:its] #if temp is negative, out of water
temp[2:t,]=h20[1:t-1,]+capture[2:t,]-demand[2:t,]

Simulation Runs

for (i in 1:t){
  overflow[i,]=max(0,temp[i,]-tk) #either no overflow or the temp calculation less the tank size
  out_of_water[i,]=max(0,0-temp[i,]) #either not out  of water or the amount below zero
  h20[i,]=temp[i,]-overflow[i,]+out_of_water[i,] #subtract the overflow and add back in the out of water
}

Histogram

par(mfrow=c(2,2))
hist(h20, main='Histogram of Water in Tank', xlab='Water')
hist(h20[,1])