The problem with using a simple arima approach is that model identification is done when it is based upon the tacit assumption that there is no deterministic structure in the data and model parameters are invariant over time.
Your rainfall data (216 monthly values) has two significant seasonal dummies ( October and November) and a few pulses ( 6 of them all in the 4th quarter of the year ) .. one time anomalies and an arima structure (0,0,0)(1,0,0,). You need to integrate both memory (arima) and latent deterministic structure to get a useful/thorough model which has decomposed the observed data to signal and white noise i.e. an error term free of identifiable structure .
Your data 
has an ACF 
. A useful model is here 
and here 
with model summary stats here 
with forecasts here for the next 36 months 
with an ACF of here suggesting model sufficiency 
The model's residuals are here 
The Actual/Fit and Forecast is here 
You can probably come close to this result by using the TSOUTLIER package in R . I would closely examine the resultant error term to assess it's randomness i.e. the sufficiency of your model.
Hope this helps . Sorry , I do not write R code .
Note that the prediction limits are based upon bootstrapping the model residuals and you can safely truncate the lower values to be => 0.0 .
You asked for 95% prediction limits ..
REVISION OF FORECAST PREDICTION LIMITS:
Using standard computation assuming normality of the errors and no future anomalies this is the Actual/Fit and Forecast 
providing more realistic limits ... again simply change the lower limit if it is less than zero to zero.
EDITED IN RESPONSE TO YOUR COMMENT AS TO WHY AUTO.ARIMA developed it's model:
auto.arima can work if there are no latent deterministic structures . Your data has them in a major way. Simply trying a brute-force AIC optimization fails because it uses the ACF while "The correlogram should be calculated from residuals using a model that controls for intervention administration, otherwise the intervention effects are taken to be Gaussian noise, underestimating the actual autoregressive effect." see @Adam0's comment on Interrupted Time Series Analysis - ARIMAX for High Frequency Biological Data?.
auto.arima incorrectly uses seasonally differences and a seasonal ma term while the correct fix for seasonal non-stationarity (in this case) is to identify and incorporate fixed effects for the two months (Oct and Nov) . Being unable to identify large pulses and there are a few, it used an ma(1) model memory rather than fixed pulse effects.
Finally model identification is an iterative , multi-state approach like peeling an onion https://autobox.com/pdfs/ARIMA%20FLOW%20CHART.pdf .
