water.update <- function(WAT0, RAIN, ETo){

            S = 25400/CN - 254;  IA = 0.2*S

            if(RAIN > IA){RO = (RAIN - 0.2 * S)^2/(RAIN + 0.8 * S)
            } else {
            RO = 0
            }

            if(WAT0 + RAIN - RO > FC) {DR = DC * (WAT0 + RAIN - RO - FC)
            } else {
            DR = 0
            }
            dWAT = RAIN - RO - DR - ETo
            WAT1 = WAT0 + dWAT
            return(c(WAT1,RO,DR))
          }

CN <- 60;FC <- 42;DC <- 0.02
water.update(WAT0 = 23, RAIN = 5, ETo = 2)
# 26, 0, 0

weather <- data.frame(day = 1:10 ,rain = sample(1:100, 10, replace = T), ETo = sample(1:10, 10, replace = T))

water.model <- function(weather, FC, CN, WAT0){

    WAT <- data.frame(matrix(NA, nrow =  nrow(weather), ncol = 3))
    WAT[1,1] <- WAT0 # WAT0 is a constant

    for(day in 1:(nrow(weather)-1)){

      WAT[day + 1,] = water.update(WAT[day,1],weather$rain[day],weather$ETo[day])
    }
    return(WAT)
    }

WAT0 <- 20
water.model(weather = weather, FC = FC, CN = CN, WAT0 = WAT0)

big.data <- data.frame(loc.id = rep(1:3, each = 10*3),
                            year = rep(rep(1981:1983, each = 10),times = 3),
                            day = rep(1:10, times = 3*3),
                            CN = rep(c(50,55,58), each = 10*3), # each location has a contant CN, FC and DC
                            FC = rep(c(72,76,80),each = 10*3),
                            DC = rep(c(0.02,0.5,0.8), each = 10*3),
                            WAT0 = rep(c(20,22,26), each = 10*3),
                            rain = sample(1:100,90, replace = T),
                            eto = sample(1:10,90, replace = T))

big.data %>% group_by(loc.id, year) %>% do??