explorations/extractFuns/example/FAO56_dualcropcoeff.R

In duncantl/CodeAnalysis: Tools for Static Analysis of R code

#TO-DO:
  #1.  Re-run all 30% and 80% allowable depletion scenarios because of correction of KsCalc function.
  #2.  Correct Ei and Ep functions such that TEW cannot be exceeded by Dei or Dep, respectively, in the upper layer [DONE].  
  #3.  Re-run all almond and walnut scenarios, since model codes were changed when tomatoes, etc. were added to the mix [DONE]
  #5.  Modify winter Kcb calculation for alfalfa: allow for fall growth in intermountain region before frost termination; allow for higher peak values in Fall then decline to 0.8 before rising again in Feb [DONE]
  #6.  Modify Fc equation for alfalfa, so that it is directly based off of Kcb values [DONE]
  #7.  Allow for winter residual "mulch" in intermountain alfalfa production (for every 10% effective surface coverage results in 5% reduction in TEW)
  #8.  Modify irrigation decision so that irrigations can't occur before 2/?? in Central Valley, even with drought stress but ensure this doesn't adversely affect results calculations that depend on definitions of Jdev [DONE]
  #9. Add alfalfa zone to model scaffold [DONE]
  #10. Add grape zone to model scaffold [DONE]
#modified KeiCalc and KepCalc on 9/11/17 to correct for overestimation of evaporation from sandy soils with very low TEW (i.e. <12 mm)
  #11. Simple tests or monte-Carlo simulation of effect of following effects (1) climatic adjustment of Kcb values (2) bloom and leaf-drop assumptions (3) Kcb values themselves
  #12. Test effect of varying alfalfa height on Kcb calcs
#changed order of DPei and DPep on 9/11/17  
# changed order of Ir and Ks calculation on 8/23/17
# changed Ir decision function on 8/23/17 to accomodate different irrigation decisions for wine grapes
# concept for wine irrigation decisions is to set a min Ks threshold and irrigate to allowable depletion when that threshold is crossed, as opposed to irrigating to field capacity when allowable depletion is crossed
#script to implement the FAO56 dual crop coefficient ET routine
modelscaffoldDir <- 'C:/Users/smdevine/Desktop/Allowable_Depletion/model_scaffold/run_model/Oct2017' #location of input data; copied from Sep2017 on 10/18/17
resultsDir <- 'C:/Users/smdevine/Desktop/Allowable_Depletion/results/Oct2017'
rounding_digits <- 3
setwd(modelscaffoldDir)
irrigation.parameters <- read.csv('irrigation_parameters.csv', stringsAsFactors = F)
crop.parameters.df <- read.csv('crop_parameters.csv', stringsAsFactors = F) #necessary parameters by crop to run the FAO56 dual crop coefficient ET model
P.df <- read.csv('PRISM.precip.data.updated9.13.17.csv', stringsAsFactors = F) #this is a daily summary of precip from 10/1/2003-6/25/17 from 'free' daily PRISM 4km resolution for cells of interest in California, created in download_PRISM.R script (from 6/26/17 download); blanks checked for in 'data_QA_QC.R'
U2.df <- read.csv('SpatialCIMIS.U2.updated9.13.17.csv', stringsAsFactors = F) #this is a daily summary of wind data from download of spatial CIMIS data, created in spatialCIMIS.R script.  No missing data except for cell 148533
RHmin.df <- read.csv('SpatialCIMIS.RHmin.updated9.13.17.csv', stringsAsFactors = F) #this is a daily summary of minimum relative humidity, estimated from download of spatial CIMIS Tdew and Tmax data, created in spatialCIMIS.R script.  Blanks filled on "12_08_2011" in data_QA_QC.R.  Now, no missing data except for cell 148533
ETo.df <- read.csv('SpatialCIMIS.ETo.updated9.13.17.csv', stringsAsFactors = F) #this is a daily summary of reference ET from download of spatial CIMIS data, created in spatialCIMIS.R script.  Blanks filled on multiple days in data_QA_QC.R.  Now, no missing data except for cell 148533
model.scaffold <- read.csv('model_scaffold_majcomps.v2.csv', stringsAsFactors = F)
model.scaffold <- model.scaffold[order(model.scaffold$unique_model_code), ]
model.scaffold$longitude_AEA <- NULL
model.scaffold$latitude_AEA <- NULL
model.scaffold$grape.zone[which(model.scaffold$grape.zone=='Sonoran Basin and Range')] <- 'Central California Valley' #this zone will be treated like 'Central California Valley,' that is, table grape, raisin, or low-quality wine grape
cropscape_legend <- read.csv('cropscape_legend.txt', stringsAsFactors = FALSE)
# sum(model.scaffold$crop_code==grape_code) #80312
# sum(model.scaffold$grape.zone=='Central California Valley' | model.scaffold$grape.zone=='Central California Foothills and Coastal Mountains', na.rm=TRUE) #79767
# dim(model.scaffold) #387970 rows

#now merge SpCIMIS data updates and write to disk when necessary
# setwd(file.path(modelscaffoldDir, 'SpCIMIS'))
# U2.update <- read.csv('SpatialCIMIS.U2update.rounded.csv')
# ETo.update <- read.csv('SpatialCIMIS.EToupdate.rounded.csv')
# RHmin.update <- read.csv('SpatialCIMIS.minRHupdate.rounded.csv', stringsAsFactors = F)
# setwd(modelscaffoldDir)
# P.update <- read.csv('PRISM_precip_data_update.csv', stringsAsFactors = FALSE)
# dim(U2.update)
# dim(U2.df)
# U2.df <- rbind(U2.df, U2.update)
# dim(U2.df)
# dim(ETo.df)
# dim(ETo.update)
# ETo.df <- rbind(ETo.df, ETo.update)
# dim(RHmin.df)
# dim(RHmin.update)
# RHmin.df <- rbind(RHmin.df, RHmin.update)
# dim(P.df)
# P.df[nrow(P.df),1:20] #last date was 6/25/17 from June 26 2017 download
# P.df <- P.df[-(which(P.df$dates=='12_01_2016'):nrow(P.df)), ]
# dim(P.df)
# P.df[nrow(P.df),1:20]
# P.df <- rbind(P.df, P.update)
# write.csv(P.df, 'PRISM.precip.data.updated9.13.17.csv', row.names=FALSE)
# write.csv(U2.df, 'SpatialCIMIS.U2.updated9.13.17.csv', row.names=FALSE)
# write.csv(ETo.df, 'SpatialCIMIS.ETo.updated9.13.17.csv', row.names=FALSE)
# write.csv(RHmin.df, 'SpatialCIMIS.RHmin.updated9.13.17.csv', row.names=FALSE)

#define functions implement FA56 dual crop coefficients
#includes subroutine that separates evaporable water as P vs. Irr sourced
#some results abbreviations
#E=evaporation
#T=transpiration
#ET=evapotranspriation
#Irr.1=first irrigation of year
#Irr.Last=last irrigation of year
#GW.ET.to.Irr1=ET sourced from green water from flowering to first irrigation (should subtract previous year's carryover irrigation storage)
#ET.growing=total ET from flowering (Jdev) to leaf-drop (Jharv)
#ET.WY=total annual ET on water-year basis
#deep.perc is annual deep percolation ()
#GW.capture.net is net change in soil root zone depletion from Jharv (leaf drop) to Jdev (flowering and development)
#end.season.Dr is soil root zone depletion at Jharv (leaf drop)

#temporary arguments for working inside the function when necessary
# cropname <- 'alfalfa.intermountain'
# cropcode <- alfalfa_code
# AD.percentage <- 30
# root_depth <- '1.0m'
# irr.type <- 'Border'
# results_file <- 'new'
# row_start <- 1
# RDI.min <- NA
# alfalfa.zone <- 'Intermountain'
# alfalfa_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Alfalfa'] #75380 total
# grape_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Grapes']
# almond_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Almonds']
# walnut_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Walnuts']
# pistachio_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Pistachios']
FAO56DualCropCalc <- function(cropname, cropcode, AD.percentage, root_depth, irr.type, crop.parameters.df, model.scaffold, U2.df, P.df, ETo.df, RHmin.df, results_file, row_start, RDI.min, alfalfa.zone, grape.zone) {
  alfalfa_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Alfalfa']
  grape_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Grapes']
  almond_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Almonds']
  walnut_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Walnuts']
  pistachio_code <- cropscape_legend$VALUE[cropscape_legend$CLASS_NAME=='Pistachios']
  CropParametersDefine <- function(crop.parameters, cropname) {
    crop.parameters <- crop.parameters[which(crop.parameters$crop==cropname), ]  
    bloom.date <- strptime(paste0(as.character(crop.parameters$bloom.mo), '/', as.character(crop.parameters$bloom.day)), '%m/%d')
    if (cropname=='alfalfa.intermountain' | cropname=='alfalfa.CV' | cropname=='alfalfa.imperial') {
      crop.parameters$Jini <- as.integer(format.Date(bloom.date, '%j')) #this covers the first cutting
      crop.parameters$Jdev <- crop.parameters$Jini + crop.parameters$Lini
      crop.parameters$Jmid <- crop.parameters$Jdev + crop.parameters$Ldev
      crop.parameters$Jlate <- crop.parameters$Jmid + crop.parameters$Lmid
      crop.parameters$Jharv <- crop.parameters$Jlate + crop.parameters$Llate
      return(crop.parameters)
    } else { #these cover the orchard and vine crops
      crop.parameters$Jdev <- as.integer(format.Date(bloom.date, '%j'))
      crop.parameters$Jmid <- crop.parameters$Jdev + crop.parameters$Ldev
      crop.parameters$Jlate <- crop.parameters$Jmid + crop.parameters$Lmid
      crop.parameters$Jharv <- crop.parameters$Jlate + crop.parameters$Llate
      return(crop.parameters)
    }
  }
  KcbDefine <- function(doy, crop) { #find standard value of Kcb by day of year relative to three reference Kcb points defined by crop parameters
    parameters <- crop.parameters[which(crop.parameters$crop==crop), ]
    if (crop=='alfalfa.intermountain' | crop=='alfalfa.CV' | crop=='alfalfa.imperial') {
      cycle.1.length <- sum(parameters[['Lini']], parameters[['Ldev']], parameters[['Lmid']], parameters[['Llate']])
      alfalfa.Kcb.cycle1 <- function(Kcb.winter, final.else) {
        if (crop=='alfalfa.intermountain') {
          ifelse(doy < parameters[['Jdev']], Kcb.winter, ifelse(doy < parameters[['Jmid']], Kcb.winter + (doy - parameters[['Jdev']]) / parameters[['Ldev']] * (parameters[['Kcb.mid']] - Kcb.winter), ifelse(doy < parameters[['Jlate']] , parameters[['Kcb.mid']], ifelse(doy < parameters[['Jharv']], parameters[['Kcb.mid']] + (doy - parameters[['Jlate']])/ parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), final.else))))
        }
        else if (crop=='alfalfa.CV') {
          ifelse(doy < parameters[['Jini']], Kcb.winter, ifelse(doy < parameters[['Jdev']] , Kcb.winter, ifelse(doy < parameters[['Jmid']], Kcb.winter + (doy - parameters[['Jdev']]) / parameters[['Ldev']] * (parameters[['Kcb.mid']] - Kcb.winter), ifelse(doy < parameters[['Jlate']] , parameters[['Kcb.mid']], ifelse(doy < parameters[['Jharv']], parameters[['Kcb.mid']] + (doy - parameters[['Jlate']])/ parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), final.else)))))
        } else { #then covers Imperial valley production just before and after Jan cutting
            ifelse(doy < parameters[['Jini']] - parameters[['Llate']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jini']], parameters[['Kcb.mid']] + (doy - parameters[['Jini']] + parameters[['Llate']])/ parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), ifelse(doy < parameters[['Jdev']] , parameters[['Kcb.ini']], ifelse(doy < parameters[['Jmid']], parameters[['Kcb.ini']] + (doy - parameters[['Jdev']]) / parameters[['Ldev']] * (parameters[['Kcb.mid']] - parameters[['Kcb.ini']]), ifelse(doy < parameters[['Jlate']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jharv']], parameters[['Kcb.mid']] + (doy - parameters[['Jlate']])/ parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), final.else))))))
        }
      }
      alfalfa.Kcb.cycles <- function(cycle.no, cycle.length, final.else) {
        ifelse(doy < parameters[['Jini']] + cycle.1.length + cycle.length*(cycle.no - 2) + parameters[['Lini.cycle']], parameters[['Kcb.ini']], ifelse(doy < parameters[['Jini']] + cycle.1.length + cycle.length*(cycle.no - 2) + parameters[['Lini.cycle']] + parameters[['Ldev.cycle']], parameters[['Kcb.ini']] + (doy - parameters[['Jini']] - cycle.1.length - cycle.length*(cycle.no - 2) - parameters[['Lini.cycle']]) / parameters[['Ldev.cycle']] * (parameters[['Kcb.mid']] - parameters[['Kcb.ini']]), ifelse(doy < parameters[['Jini']] + cycle.1.length + cycle.length*(cycle.no - 2) + parameters[['Lini.cycle']] + parameters[['Ldev.cycle']] + parameters[['Lmid.cycle']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jini']] + cycle.1.length + cycle.length*(cycle.no - 2) + parameters[['Lini.cycle']] + parameters[['Ldev.cycle']] + parameters[['Lmid.cycle']] + parameters[['Llate.cycle']], parameters[['Kcb.mid']] + (doy - parameters[['Jini']] - cycle.1.length - cycle.length*(cycle.no - 2) - parameters[['Lini.cycle']] - parameters[['Ldev.cycle']] - parameters[['Lmid.cycle']]) / parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), final.else))))
      }
      CV.alfalfa.last <- function(final.else) {
        ifelse(doy < parameters[['Jini']] + cycle.1.length + parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) + parameters[['Lini']], parameters[['Kcb.ini']], ifelse(doy < parameters[['Jini']] + cycle.1.length + parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) + parameters[['Lini']] + parameters[['Ldev']], parameters[['Kcb.ini']] + (doy - parameters[['Jini']] - cycle.1.length - parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) - parameters[['Lini']]) / parameters[['Ldev']] * (parameters[['Kcb.mid']] - parameters[['Kcb.ini']]), ifelse(doy < parameters[['Jini']] + cycle.1.length + parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jini']] + cycle.1.length + parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']] + parameters[['Llate']], parameters[['Kcb.mid']] + (doy - parameters[['Jini']] - cycle.1.length - parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) - parameters[['Lini']] - parameters[['Ldev']] - parameters[['Lmid']]) / parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), final.else))))
      }
      winter.Kcb.alfalfa <- function(cycle.no=1, final.else=NA) { #could add dormancy reduction factor here
        if (crop=='alfalfa.CV') {
          ifelse(doy < parameters[['Jharv']] + cycle.1.length + (parameters[['cycle.no']] - 2) * parameters[['cycle.length']] + parameters[['Lini']], parameters[['Kcb.ini']], ifelse(doy < parameters[['Jharv']] + cycle.1.length + (parameters[['cycle.no']] - 2) * parameters[['cycle.length']] + parameters[['Lini']] + parameters[['Ldev']], parameters[['Kcb.ini']] + (doy - parameters[['Jharv']] - cycle.1.length - (parameters[['cycle.no']] - 2) * parameters[['cycle.length']] - parameters[['Lini']]) / parameters[['Ldev']] * (parameters[['Kcb.winter']] - parameters[['Kcb.ini']]), ifelse(doy < parameters[['Jharv']] + cycle.1.length + (parameters[['cycle.no']] - 2) * parameters[['cycle.length']] + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jharv']] + cycle.1.length + (parameters[['cycle.no']] - 2) * parameters[['cycle.length']] + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']] + parameters[['Llate']] + parameters[['Llate']], parameters[['Kcb.mid']] + (doy - parameters[['Jharv']] - cycle.1.length - parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) - parameters[['Lini']] - parameters[['Ldev']] - parameters[['Lmid']]) / (parameters[['Llate']] + parameters[['Llate']]) * (parameters[['Kcb.winter']] - parameters[['Kcb.mid']]), parameters[['Kcb.winter']]))))
        }
        else if (crop=='alfalfa.intermountain') {
          ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 1) * parameters[['cycle.length']] + parameters[['Lini']], parameters[['Kcb.ini']], ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 1) * parameters[['cycle.length']] + parameters[['Lini']] + parameters[['Ldev']], parameters[['Kcb.ini']] + (doy - parameters[['Jharv']] - (parameters[['cycle.no']] - 1) * parameters[['cycle.length']] - parameters[['Lini']]) / parameters[['Ldev']] * (parameters[['Kcb.mid']] - parameters[['Kcb.ini']]), ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 1) * parameters[['cycle.length']] + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 1) * parameters[['cycle.length']] + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']] + parameters[['Lfrost.kill']], parameters[['Kcb.mid']] + (doy - parameters[['Jharv']] - (parameters[['cycle.no']] - 1) * parameters[['cycle.length']] - parameters[['Lini']] - parameters[['Ldev']] - parameters[['Lmid']]) / parameters[['Lfrost.kill']] * (parameters[['Kcb.dorm']] - parameters[['Kcb.mid']]), parameters[['Kcb.dorm']]))))
        } else { #this covers Imperial Valley
          ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 3) * parameters[['cycle.length']] + (cycle.no - parameters[['cycle.no']] + 1) * cycle.1.length + parameters[['Lini']], parameters[['Kcb.ini']], ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 3) * parameters[['cycle.length']] + (cycle.no - parameters[['cycle.no']] + 1) * cycle.1.length + parameters[['Lini']] + parameters[['Ldev']], parameters[['Kcb.ini']] + (doy - parameters[['Jharv']] - (parameters[['cycle.no']] - 3) * parameters[['cycle.length']] - (cycle.no - parameters[['cycle.no']] + 1) * cycle.1.length - parameters[['Lini']]) / parameters[['Ldev']] * (parameters[['Kcb.mid']] - parameters[['Kcb.ini']]), ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 3) * parameters[['cycle.length']] + (cycle.no - parameters[['cycle.no']] + 1) * cycle.1.length + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 3) * parameters[['cycle.length']] + (cycle.no - parameters[['cycle.no']] + 1) * cycle.1.length + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']] + parameters[['Llate']], parameters[['Kcb.mid']] + (doy - parameters[['Jharv']] - (parameters[['cycle.no']] - 3) * parameters[['cycle.length']] - (cycle.no - parameters[['cycle.no']] + 1) * cycle.1.length - parameters[['Lini']] - parameters[['Ldev']] - parameters[['Lmid']]) / parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), final.else))))
        }
      }
      if (crop=='alfalfa.intermountain') {
        alfalfa.Kcb.cycle1(parameters[['Kcb.dorm']], alfalfa.Kcb.cycles(2, parameters[['cycle.length']], alfalfa.Kcb.cycles(3, parameters[['cycle.length']], winter.Kcb.alfalfa())))
      }
      else if (crop=='alfalfa.CV') {
        alfalfa.Kcb.cycle1(parameters[['Kcb.winter']], alfalfa.Kcb.cycles(2, parameters[['cycle.length']], alfalfa.Kcb.cycles(3, parameters[['cycle.length']], alfalfa.Kcb.cycles(4, parameters[['cycle.length']], alfalfa.Kcb.cycles(5, parameters[['cycle.length']], alfalfa.Kcb.cycles(6, parameters[['cycle.length']], CV.alfalfa.last(winter.Kcb.alfalfa())))))))
      } else { #crop is alfalfa.imperial
        alfalfa.Kcb.cycle1(Kcb.winter=NA, alfalfa.Kcb.cycles(2, parameters[['cycle.length']], alfalfa.Kcb.cycles(3, parameters[['cycle.length']], alfalfa.Kcb.cycles(4, parameters[['cycle.length']], alfalfa.Kcb.cycles(5, parameters[['cycle.length']], alfalfa.Kcb.cycles(6, parameters[['cycle.length']], alfalfa.Kcb.cycles(7, parameters[['cycle.length']], alfalfa.Kcb.cycles(8, parameters[['cycle.length']], alfalfa.Kcb.cycles(9, parameters[['cycle.length']], winter.Kcb.alfalfa(cycle.no = 10, final.else = winter.Kcb.alfalfa(cycle.no = 11, final.else=NA)))))))))))
      }
    } else {
        ifelse(doy < parameters[['Jdev']], parameters[['Kcb.dorm']], ifelse(doy < parameters[['Jmid']], parameters[['Kcb.ini']] + (doy - parameters[['Jdev']]) / parameters[['Ldev']] * (parameters[['Kcb.mid']] - parameters[['Kcb.ini']]), ifelse(doy < parameters[['Jlate']], parameters[['Kcb.mid']], ifelse(doy < parameters[['Jharv']], parameters[['Kcb.mid']] + (doy - parameters[['Jlate']]) / parameters[['Llate']] * (parameters[['Kcb.end']] - parameters[['Kcb.mid']]), parameters[['Kcb.dorm']]))))
    }
  }
  #adjust standard Kcb values by daily climate according to equation #5 in Allen et al. 2005 and calculate Kc max at the same time
  KcbAdj <- function(Kcb.daily, crop.parameters, crop, U2, RHmin){
    height <- crop.parameters$height[which(crop.parameters$crop==crop)]
    Kcb.climate.adj <- Kcb.daily + (0.04 * (U2 - 2) - 0.004 * (RHmin - 45)) * (height / 3) ^ 0.3
    Kcb.climate.adj[which(Kcb.climate.adj < 0)] <- 0
    Kc.max <- pmax(1.2 + (0.04 * (U2 - 2) - 0.004 * (RHmin - 45)) * (height / 3) ^ 0.3, Kcb.climate.adj + 0.05)
    return(as.data.frame(cbind(Kcb.climate.adj, Kc.max)))
  }
  #calculate fraction of cover ('fc') for specific day of year relative to three points defined in crop parameters
  fcCalc <- function(doy, parameters, crop) {
    parameters <- crop.parameters[which(crop.parameters$crop==crop), ]
    if (crop == 'alfalfa.CV' | crop == 'alfalfa.imperial') {
      Kcb.std/parameters[['Kcb.mid']]
    }
    else if (crop == 'alfalfa.intermountain') {
      ifelse(doy < parameters[['Jdev']], 0, ifelse(doy < parameters[['Jharv']] + (parameters[['cycle.no']] - 1) * parameters[['cycle.length']] + parameters[['Lini']] + parameters[['Ldev']] + parameters[['Lmid']] + parameters[['Lfrost.kill']], Kcb.std/parameters[['Kcb.mid']], 0))
    } else {
        ifelse(doy < parameters[['Jdev']], parameters[['fc.ini']], ifelse(doy < parameters[['Jmid']], parameters[['fc.ini']] + (doy - parameters[['Jdev']]) / parameters[['Ldev']] * (parameters[['fc.mid']] - parameters[['fc.ini']]), ifelse(doy < parameters[['Jlate']], parameters[['fc.mid']], ifelse(doy < parameters[['Jharv']], parameters[['fc.mid']] + (doy-parameters[['Jlate']]) / parameters[['Llate']] * (parameters[['fc.end']] - parameters[['fc.mid']]), parameters[['fc.ini']]))))
    }
  }
  #calculate fraction of soil wetted by both irrigation and precipitation and exposed to rapid drying (fewi) and fraction of soil exposed to rapid drying and is wetted by precipitation only (fewp).  These are dependent upon fraction of cover calculation above
  fwSelect <- function(irr.parameters, irr.type) {
    irr.parameters$fw[which(irr.parameters$irrigation.type==irr.type)]
  }
  fewiCalc <- function(fc, fw) { #see p.147 of FAO-56 Chp.7 re: drip
    fewi.temp <- pmin(1 - fc, fw)
    fewi.temp[which(fewi.temp < 0.001)] <- 0.001 #lower limit on fewi for numeric stability
    fewi.temp[which(fewi.temp > 1)] <- 1 #upper limit on fewi for numeric stability
    return(fewi.temp)
  }
  fewpCalc <- function(fc, fewi) {
    fewp.temp <- 1 - fc - fewi
    fewp.temp[which(fewp.temp < 0.001)] <- 0.001 #lower limit on fewp for numeric stability
    fewp.temp[which(fewp.temp > 1)] <- 1 #upper limit on fewp for numeric stability
    return(fewp.temp)
  }
  # TEWCalc <- function(ThetaFC, ThetaWP, REW, Ze=0.10) { #Ze is depth of upper soil 
  #   #layer where evaporation occurs in meters
  #   result <- 1000 * (ThetaFC - 0.5 * ThetaWP) * Ze
  #   if (result < REW) {
  #     stop(print('Hay una problema con TEW')) #could change to next and print or 
  #     #save model code
  #   } 
  #   return(result)
  # }
  DepInitialCalc <- function(Dep.end, P, i) {
    max(Dep.end[i-1] - P[i], 0) # DON'T RUN THIS IF i=1
  }
  DeiInitialCalc <- function(Dei.end, P, Ir, fw, i) { #DON'T RUN THIS WHEN i=1
    max(Dei.end[i - 1] - P[i] - Ir[i - 1] / fw, 0) #have to use irrigation from 
    #previous day because current day irrigation decision is dependent on this calc
  }
  KrCalc <- function(TEW, REW, De.initial, i) { #can be used to calculate Kri or Krp
    max(0, if(De.initial[i] < REW) { #could initialize this in vector above
      1
    } else if(TEW==REW) {
      0
    } else {
      (TEW - De.initial[i]) / (TEW - REW)
    })
  }
  WCalc <- function(TEW, Dei.initial, Dep.initial, fewp, fewi, i) {
    1 / (1 + (fewp[i] / fewi[i]) * max(TEW - Dep.initial[i], 0.001) / max((TEW - Dei.initial[i]), 0.001))
  }
  KeiCalc <- function(Kri, W, Kcmax, Kcb, fewi, TEW, Dei.end, DPei, P, ETo, i) {
    Kei.est <- min(Kri[i] * W[i] * (Kcmax[i] - Kcb[i]), fewi[i] * Kcmax[i])
    #print(Kei.est)
    TEW.check <- max(Dei.end[i - 1] - P[i] - Ir[i - 1] / fw + Kei.est * ETo[i] / fewi[i] + DPei[i], 0)
    #print(TEW.check)
    if (TEW.check > TEW) {
      return(fewi[i] * (TEW - Dei.end[i - 1] + P[i] + Ir[i - 1] / fw - DPei[i]) / ETo[i]) #as revised Kep estimate
    } else {
      return(Kei.est)
    }
  }
  KepCalc <- function(Krp, W, Kcmax, Kcb, fewp, TEW, Dep.end, DPep, P, ETo, i) {
    Kep.est <- min(Krp[i] * (1 - W[i]) * (Kcmax[i] - Kcb[i]), fewp[i] * Kcmax[i])
    #print(Kep.est)
    TEW.check <- max(Dep.end[i - 1] - P[i] + (Kep.est * ETo[i]) / fewp[i] + DPep[i], 0)
    #print(TEW.check)
    if (TEW.check > TEW) {
      return(fewp[i] * (TEW - Dep.end[i - 1] + P[i] - DPep[i]) / ETo[i]) #as revised Kep estimate
    } else {
        return(Kep.est)
    }
  }
  EpCalc <- function(ETo, Kep, i) {
    ETo[i] * Kep[i]
  }
  EiCalc <- function(ETo, Kei, i) {
    ETo[i] * Kei[i]
  }
  DPepCalc <- function(P, Dep.end, i) { #DON'T RUN THIS FOR i=1
    max(P[i] - Dep.end[i-1], 0)
  }
  DepEndCalc <- function(Dep.end, P, Ep, fewp, DPep, i) { #DON'T RUN THIS FOR i=1
    Dep.end.est <- max(Dep.end[i - 1] - P[i] + Ep[i] / fewp[i] + DPep[i], 0)
  } #ignores transpiration and runoff from upper layer
  DPeiCalc <- function(P, Ir, fw, Dei.end, i) { #DON'T RUN THIS FOR i=1
    max(0, P[i] + Ir[i - 1] / fw - Dei.end[i-1])
  }
  DeiEndCalc <- function(Dei.end, P, Ir, fw, fewi, Ei, DPei, i) { #DON'T RUN THIS FOR i=1
    max(Dei.end[i - 1] - P[i] - Ir[i - 1] / fw + Ei[i] / fewi[i] + DPei[i], 0)
  } #again, ignoring tranpiration from upper 10 cm and runoff, eqn. 21)
  KcnsCalc <- function(Kcb, Kei, Kep, i) {
    Kcb[i] + Kei[i] + Kep[i]
  }
  ETcnsCalc <- function(Kc.ns, ETo, i) {
    Kc.ns[i] * ETo[i]
  }
  DrInitialCalc <- function(Dr.end, ETc.ns, P, Ir, i) { #DON'T RUN THIS FOR i=1
    Dr.end[i - 1] - P[i] - Ir[i-1] + ETc.ns[i]
  }
  #need to refine days.no.irr based on ETo data for a given location
  DaysNoIrr <- function(P, ETo, Kcb.adjusted, AD, doys.model, years, Jmid, Jharv, buffer.days=if(cropname=='alfalfa.intermountain') {53} else {30}) { #Jharv(cycle 1) + (cycle.no - 1) * cycle.length + Lini + Ldev + Lmid + L.frost.kill * 0.444 = 298, in ref. to alfalfa.intermountain; doy 322 for alfalfa.CV from 292 [last day in harvest.days] + 30 [cycle.length]
    if (cropname=='alfalfa.imperial') {
      return(0)
    }
    else if (cropname == 'alfalfa.CV' | cropname == 'alfalfa.intermountain') {
      last.irr.possible <- harvest.days[length(harvest.days)] + buffer.days
      df <- data.frame(cbind(Kcb.adjusted, ETo, P, doys.model, years))
      df <- df[which(df$doys.model >= 200 & df$doys.model <= last.irr.possible), ]
      df$ETcb <- df$Kcb.adjusted * df$ETo
      daily.ETcb <- tapply(df$ETcb, df$doys.model, mean)
      daily.P <- tapply(df$P, df$doys.model, mean)
      daily.WB <- daily.ETcb - daily.P
      for (i in 1:length(daily.WB)) {
        if (sum(daily.WB) < AD) {
          last.irr <- as.integer(names(daily.WB[1]))
          while(last.irr %in% harvest.days) {
            last.irr <- last.irr + 1
          }
          return(last.irr.possible - last.irr)
        } else {
            daily.WB <- daily.WB[-1]
        }
      }
    }
    else if (cropname=='grapes.wine') {
      return(30)
    } else {
      df <- data.frame(cbind(Kcb.adjusted, ETo, P, doys.model, years))
      df <- df[which(df$doys.model >= Jmid & df$doys.model <= Jharv), ]
      df$ETcb <- df$Kcb.adjusted * df$ETo
      daily.ETcb <- tapply(df$ETcb, df$doys.model, mean)
      daily.P <- tapply(df$P, df$doys.model, mean)
      daily.WB <- daily.ETcb - daily.P
      for (i in 1:length(daily.WB)) {
        if (sum(daily.WB) < AD) {
          return(Jharv - as.integer(names(daily.WB[1])))
        } else {
          daily.WB <- daily.WB[-1]
        }
      }
    }
  }
  GrapeLastIrr <- function() {
    if (cropname == 'grapes.wine') {
      df <- data.frame(cbind(Kcb.adjusted, ETo, P, doys.model))
      df <- df[which(df$doys.model >= Jharv - days.no.irr & df$doys.model <= Jharv), ]
      df$ETcb <- df$Kcb.adjusted * df$ETo
      daily.ETcb <- tapply(df$ETcb, df$doys.model, mean)
      daily.P <- tapply(df$P, df$doys.model, mean)
      last.irr <- Dr.initial[i] + sum(daily.ETcb) - sum(daily.P) - AD
      if (last.irr > Dr.initial[i]) {
        return(Dr.initial[i])
      } else {
        return(max(last.irr, 0))
      }
    }
  }
  HarvestDays <- function(crop) { #assume 4 day buffer on either end of each harvest
    if (crop == 'alfalfa.imperial') {
      parameters <- crop.parameters[which(crop.parameters$crop==crop), ]
      harvest.dates <- c(parameters[['Jini']], seq(from=parameters[['Jharv']], to=parameters[['Jharv']] + parameters[['cycle.length']]*(parameters[['cycle.no']]-3), by=parameters[['cycle.length']]), parameters[['Jharv']] + parameters[['cycle.length']]*(parameters[['cycle.no']]-3) + (parameters[['Jharv']] - parameters[['Jini']]))
      harvest.dates.all <- c(harvest.dates, harvest.dates-1, harvest.dates-2, harvest.dates-3, harvest.dates-4, harvest.dates+1, harvest.dates+2, harvest.dates+3, harvest.dates+4)
      harvest.dates.all[order(harvest.dates.all)]
    }
    else if (crop == 'alfalfa.CV') {
      parameters <- crop.parameters[which(crop.parameters$crop==crop), ]
      harvest.dates <- c(seq(from=parameters[['Jharv']], to=parameters[['Jharv']] + parameters[['cycle.length']]*(parameters[['cycle.no']] - 2), by=parameters[['cycle.length']]), parameters[['Jharv']] + parameters[['cycle.length']]*(parameters[['cycle.no']] - 2) + (parameters[['Jharv']] - parameters[['Jini']]))
      harvest.dates.all <- c(harvest.dates, harvest.dates - 1, harvest.dates - 2, harvest.dates - 3, harvest.dates - 4, harvest.dates + 1, harvest.dates + 2, harvest.dates + 3, harvest.dates + 4)
      harvest.dates.all[order(harvest.dates.all)]
    } else { #this is for intermountain alfalfa
      parameters <- crop.parameters[which(crop.parameters$crop==crop), ]
      harvest.dates <- seq(from=parameters[['Jharv']], to=parameters[['Jharv']] + (parameters[['cycle.no']] - 1) * parameters[['cycle.length']], by=parameters[['cycle.length']])
      harvest.dates.all <- c(harvest.dates, harvest.dates - 1, harvest.dates - 2, harvest.dates - 3, harvest.dates - 4, harvest.dates + 1, harvest.dates + 2, harvest.dates + 3, harvest.dates + 4)
      harvest.dates.all[order(harvest.dates.all)]
    }
  }
  IrCalc <- function(Ks, RDI.min, AD, Dr.initial, doys.model, Jdev, Jharv, days.no.irr, i, buffer.days=if(cropname=='alfalfa.intermountain') {53} else {30}) { #see notes defining these above in days.no.irr calc
    if (cropname == 'alfalfa.imperial') {
      if (Dr.initial[i] > AD & !(doys.model[i] %in% harvest.days)) {
        return(Dr.initial[i])
      } else {
          return(0)
        }
    } else if (cropname == 'alfalfa.CV' | cropname == 'alfalfa.intermountain') {
        if (doys.model[i] < crop.parameters$Jini[which(crop.parameters$crop==cropname)] | doys.model[i] %in% harvest.days | doys.model[i] > harvest.days[length(harvest.days)] + buffer.days - days.no.irr) {
          return(0)
        } else if(Dr.initial[i] > AD | doys.model[i] == harvest.days[length(harvest.days)] + buffer.days - days.no.irr) {
        return(Dr.initial[i])
      } else {
          return(0)
      }
    } else if (doys.model[i] < Jdev | doys.model[i] > Jharv - days.no.irr) {
      return(0)
    } else if (cropname == 'grapes.wine' & doys.model[i] < Jharv - days.no.irr & Ks[i] < RDI.min) {
      return(Dr.initial[i] - AD)
    } else if (Dr.initial[i] > AD & cropname != 'grapes.wine') {
      return(Dr.initial[i])
    } else if (doys.model[i] == Jharv - days.no.irr & cropname=='grapes.wine') {
      return(GrapeLastIrr())
    } else if (doys.model[i] == Jharv - days.no.irr) {
      return(Dr.initial[i])
    } else {
      return(0)
    }
  }
  DPrCalc <- function(P, Ir, ETc.ns, Dr.end, i) { #DON'T RUN THIS FOR i=1
    max(P[i] + Ir[i-1] - ETc.ns[i] - Dr.end[i - 1], 0)
  }
  KsCalc <- function(Dr.initial, PAW, stress.point=0.5*PAW, i) {
    if (Dr.initial[i] > stress.point) {
      (PAW - Dr.initial[i])/(PAW - stress.point)
    } else { 
      return(1)
    }
  }
  KcactCalc <- function(Ks, Kcb, Kei, Kep, i) {#equation 4 in Allen et al. 2005
    Ks[i] * Kcb[i] + Kei[i] + Kep[i]
  }
  ETcactCalc <- function(Kc.act, ETo, i) {#equation 3 in Allen et al. 2005
    Kc.act[i] * ETo[i]
  }
  DrEndCalc <- function(Dr.end, P, Ir, Kc.act, ETo, DPr, i) { #NOT TO BE USED for i=1
    Dr.end[i - 1] - P[i] - Ir[i-1] + Kc.act[i] * ETo[i] + DPr[i]
  }
  
  #results function to summarize each model's results
  #find time to first and last irrigation
  #Jdev is 1 and Jharv is 365 for alfalfa.CV and alfalfa.imperial; Jdev is Jini for alfalfa.intermountain and Jharv is doy 327 (complete dormancy), even though actual irrigations are shortened to <doy293 for alfalfa.intermountain and <doy322 for alfalfa.CV
  IrDateCalc <- function(df) { #as.Date('1900-01-01) is a proxy for NA
    first.irr.index <- if (length(which(df$Ir >0)) > 1 & df$doys.model[1] <= Jdev) {
      head(which(df$Ir > 0), 1)} else {
        if (length(which(df$Ir >0)) == 1 & df$doys.model[1] <= Jdev) {
          which(df$Ir > 0)} else {vector()}
      }
    last.irr.index <- if (length(which(df$Ir >0)) > 1 & df$doys.model[nrow(df)] >= Jharv - days.no.irr) {
      tail(which(df$Ir > 0), 1)} else {
        if (length(which(df$Ir >0)) == 1 & df$doys.model[nrow(df)] >= Jharv - days.no.irr) {
          which(df$Ir > 0)} else {vector()}
      }
    if (length(first.irr.index) == 0 & length(last.irr.index) == 0) {
      data.frame(Irr.1=as.Date('1900-01-01'), Irr.Last=as.Date('1900-01-01'))
    }
    else if (length(first.irr.index) == 0) {
      data.frame(Irr.1=as.Date('1900-01-01'), Irr.Last=df$dates[last.irr.index])
    }
    else if (length(last.irr.index) == 0) {
      data.frame(Irr.1=df$dates[first.irr.index], Irr.Last=as.Date('1900-01-01'))
    } else {
      data.frame(Irr.1=df$dates[first.irr.index], Irr.Last=df$dates[last.irr.index])
    }
  }
  #do.call(rbind, lapply(split(model.result, model.result$years), IrDateCalc))
  
  #calculate Green Water utilized
  #this asssumes that residual irrigation water storage from previous fall will not contribute to the following year's growing season ET for the purpose of calculating green water utilization.  However a correction is applied to the growing season ET for residual irrigation water storage to correctly estimate green water utilization within the same year.
  WaterBalanceCalc <- function(df) { #concept is to run by year on a data.frame trimmed to Jdev-Jharv each year
    #df <- model.result[which(model.result$years==2004),]
    last.irr.index <- if (length(which(df$Ir >0)) > 1 & df$doys.model[nrow(df)] >= Jharv - days.no.irr) {
      tail(which(df$Ir > 0), 1)} else {
        if (length(which(df$Ir >0)) == 1 & df$doys.model[nrow(df)] >= Jharv - days.no.irr) {
          which(df$Ir > 0)} else {vector()}
      }
    jharv_index <- which(df$doys.model==Jharv)
    jdev_index <- which(df$doys.model==Jdev)
    if (df$doys.model[1] > Jdev | df$doys.model[nrow(df)] < Jharv) { #if there is not a complete growing season, don't report any data
      if (length(jharv_index)==0) { #data begins after leaf-drop or ends before leaf-drop; either way, can't get entire season or end season data
        data.frame(RAW.end.season = NA, PAW.end.season = NA, Dr.end.season = NA, P.end.season=NA, Irr.end.storage = NA, GW.ET.growing = NA, Irr.app.total = NA, Irr.app.last = NA, ET.growing = NA, E.growing = NA, T.growing = NA, H2O.stress=NA)
      }
      else if (length(last.irr.index)==0) { #implies no data when last irrigation occurred but there is data at Jharv; thus can get end of season storage data
        data.frame(RAW.end.season = max(AD - df$Dr.end[jharv_index], 0), PAW.end.season = max(PAW - df$Dr.end[jharv_index], 0), Dr.end.season = df$Dr.end[jharv_index], P.end.season = NA, Irr.end.storage = NA, GW.ET.growing = NA, Irr.app.total = NA, Irr.app.last = NA, ET.growing = NA, E.growing = NA, T.growing = NA, H2O.stress=NA)
      } else { #implies data exists from at least when last irrigation occurred to Jharv, so can get everything but entire growing season data
        data.frame(RAW.end.season = max(AD - df$Dr.end[jharv_index], 0), PAW.end.season = max(PAW - df$Dr.end[jharv_index], 0), Dr.end.season = df$Dr.end[jharv_index], P.end.season = sum(df$P[last.irr.index:jharv_index]), Irr.end.storage = max(AD - df$Dr.end[jharv_index] - sum(df$P[last.irr.index:jharv_index]), 0), GW.ET.growing = NA, Irr.app.total = NA, Irr.app.last = df$Ir[last.irr.index], ET.growing = NA, E.growing = NA, T.growing = NA, H2O.stress=NA)
      }
    } else { #entire season's coverage available
      if (length(last.irr.index)==0) { #but no final irrigation occurred
        data.frame(RAW.end.season = max(AD - df$Dr.end[jharv_index], 0), PAW.end.season = max(PAW - df$Dr.end[jharv_index], 0), Dr.end.season = df$Dr.end[jharv_index], P.end.season = NA, Irr.end.storage = NA, GW.ET.growing = sum(df$ETc.act[jdev_index:jharv_index] - df$Ir[jdev_index:jharv_index]), Irr.app.total = sum(df$Ir), Irr.app.last = NA, ET.growing = sum(df$ETc.act[jdev_index:jharv_index]), E.growing = sum(df$Ei[jdev_index:jharv_index], df$Ep[jdev_index:jharv_index]), T.growing = sum(df$ETc.act[jdev_index:jharv_index] - df$Ei[jdev_index:jharv_index] - df$Ep[jdev_index:jharv_index]), H2O.stress=sum(df$ETc.ns - df$ETc.act))
      } else {
        irr_storage <- max(AD - df$Dr.end[jharv_index] - sum(df$P[last.irr.index:jharv_index]), 0)
      data.frame(RAW.end.season = max(AD - df$Dr.end[jharv_index], 0), PAW.end.season = max(PAW - df$Dr.end[jharv_index], 0), Dr.end.season = df$Dr.end[jharv_index], P.end.season = sum(df$P[last.irr.index:jharv_index]), Irr.end.storage = irr_storage, GW.ET.growing = sum(df$ETc.act[jdev_index:jharv_index] - df$Ir[jdev_index:jharv_index]) + irr_storage, Irr.app.total = sum(df$Ir), Irr.app.last = df$Ir[last.irr.index], ET.growing = sum(df$ETc.act[jdev_index:jharv_index]), E.growing = sum(df$Ei[jdev_index:jharv_index], df$Ep[jdev_index:jharv_index]), T.growing = sum(df$ETc.act[jdev_index:jharv_index] - df$Ei[jdev_index:jharv_index] - df$Ep[jdev_index:jharv_index]), H2O.stress=sum(df$ETc.ns - df$ETc.act))
      }
    }
  }
  #do.call(rbind, lapply(split(model.result, model.result$years), WaterBalanceCalc))
  
  ##this splits the overall results data.frame into subsets by year and then runs the GreenWaterCalc on each subset via lapply.  The result from each year is then bound together via rbind called by do.call; AD minus Dr.end at leaf-drop is the readily available water remaining in storage at leaf-drop; the source of this readily available water is minus precip since the last irrigation is Irr.End.Storage
  GreenWaterIrr1Calc <- function(df) { #works on a data.frame split by year
    #df <- model.result[model.result$years==2004,]
    first.irr.index <- if (length(which(df$Ir >0)) > 1 & df$doys.model[1] <= Jdev) {
      head(which(df$Ir > 0), 1)} else {
        if (length(which(df$Ir >0)) == 1 & df$doys.model[1] <= Jdev) {
          which(df$Ir > 0)} else {vector()}
      }
    jdev_index <- which(df$doys.model==Jdev)
    #print(df$years[1])
    if (df$doys.model[1] > Jdev | length(first.irr.index)==0) {
      data.frame(GW.ET.to.Irr1=NA, GW.E.to.Irr1=NA, GW.T.to.Irr1=NA)
    } else {
      data.frame(GW.ET.to.Irr1 = sum(df$ETc.act[jdev_index:first.irr.index]), GW.E.to.Irr1 = sum(df$Ei[jdev_index:first.irr.index] + df$Ep[jdev_index:first.irr.index]), GW.T.to.Irr1 = sum(df$Kcb.adjusted[jdev_index:first.irr.index] * df$Ks[jdev_index:first.irr.index] * df$ETo[jdev_index:first.irr.index]))
    }
  }
  #do.call(rbind, lapply(split(model.result, model.result$years), GreenWaterIrr1Calc))
  
  #determine deep percolation and annual water balance using subsetting by year
  DeepPercCalc <- function(df) { #assumes Jharv index is after 10/1
    #df <- model.result[which(model.result$years==2017), ]
    #print(df$years[1])
    jan.1.index <- which(df$dates==as.Date(paste0(as.character(df$years[1]), '-01-01')))
    first.irr.index <- if (length(which(df$Ir >0)) > 1 & df$doys.model[1] <= Jdev) {
      head(which(df$Ir > 0), 1)} else {
        if (length(which(df$Ir >0)) == 1 & df$doys.model[1] <= Jdev) {
          which(df$Ir > 0)} else {vector()}
      }
    last.irr.index <- if (length(which(df$Ir >0)) > 1 & df$doys.model[nrow(df)] >= Jharv - days.no.irr) {
      tail(which(df$Ir > 0), 1)} else {
        if (length(which(df$Ir >0)) == 1 & df$doys.model[nrow(df)] >= Jharv - days.no.irr) {
          which(df$Ir > 0)} else {vector()}
      }
    jul.1.index <- which(df$dates==as.Date(paste0(as.character(df$years[1]), '-07-01')))
    dec.31.index <- which(df$dates==as.Date(paste0(as.character(df$years[1]), '-12-31')))
    #print(df$years[1])
    if (df$dates[1] > as.Date(paste0(as.character(df$years[1]), '-01-01')) | df$dates[nrow(df)] < as.Date(paste0(as.character(df$years[nrow(df)]), '-12-31'))) { #this means they ain't a whole year's data
      if (length(last.irr.index) == 0 & length(first.irr.index) == 0) { #not sufficient coverage to calc anything
        data.frame(ET.annual = NA, E.annual = NA, T.annual = NA, deep.perc.annual = NA, winter.deep.perc = NA, post.Irr1.deep.perc=NA, fall.deep.perc = NA)
      }
      else if (length(jan.1.index) != 0 & length(first.irr.index) != 0 & length(jul.1.index) == 0) { #winter deep perc only
        data.frame(ET.annual = NA, E.annual = NA, T.annual = NA, deep.perc.annual = NA, winter.deep.perc = sum(df$DPr[jan.1.index:first.irr.index]), post.Irr1.deep.perc=NA, fall.deep.perc = NA)
      }
      else if (length(jan.1.index) != 0 & length(jul.1.index) != 0 & length(dec.31.index) == 0) { #winter and post Irr deep perc
        data.frame(ET.annual = NA, E.annual = NA, T.annual = NA, deep.perc.annual = NA, winter.deep.perc = sum(df$DPr[jan.1.index:first.irr.index]), post.Irr1.deep.perc=if(first.irr.index < jul.1.index) {sum(df$DPr[(first.irr.index+1):jul.1.index])}else{NA}, fall.deep.perc = NA)
      }
      else if (length(jan.1.index) == 0 & length(first.irr.index) != 0 & length(jul.1.index) != 0 & length(last.irr.index) == 0) { #only Spring deep perc available
        data.frame(ET.annual = NA, E.annual = NA, T.annual = NA, deep.perc.annual = NA, winter.deep.perc = NA, post.Irr1.deep.perc=if(first.irr.index < jul.1.index) {sum(df$DPr[(first.irr.index+1):jul.1.index])}else{NA}, fall.deep.perc = NA)
      }
      else if (length(jan.1.index) == 0 & length(first.irr.index) != 0 & length(dec.31.index) != 0) { #spring and fall deep perc
        data.frame(ET.annual = NA, E.annual = NA, T.annual = NA, deep.perc.annual = NA, winter.deep.perc = NA, post.Irr1.deep.perc=if(first.irr.index < jul.1.index) {sum(df$DPr[(first.irr.index+1):jul.1.index])}else{NA}, fall.deep.perc =  sum(df$DPr[last.irr.index:dec.31.index]))
      } else {#fall perc only
        data.frame(ET.annual = NA, E.annual = NA, T.annual = NA, deep.perc.annual = NA, winter.deep.perc = NA, post.Irr1.deep.perc=NA,  fall.deep.perc =  sum(df$DPr[last.irr.index:(dec.31.index-1)]))
      }
    } else { #entire annual coverage available
      if (length(first.irr.index)==0 & length(last.irr.index)==0) { #but no irrigations
        data.frame(ET.annual = sum(df$ETc.act), E.annual = sum(df$Ei, df$Ep), T.annual = sum(df$ETc.act, -df$Ei, -df$Ep), deep.perc.annual = sum(df$DPr), winter.deep.perc = NA, post.Irr1.deep.perc = NA, fall.deep.perc = NA)
      } else {
        data.frame(ET.annual = sum(df$ETc.act), E.annual = sum(df$Ei, df$Ep), T.annual = sum(df$ETc.act, -df$Ei, -df$Ep), deep.perc.annual = sum(df$DPr), winter.deep.perc = sum(df$DPr[jan.1.index:first.irr.index]), post.Irr1.deep.perc = if(first.irr.index < jul.1.index) {sum(df$DPr[(first.irr.index+1):jul.1.index])}else{NA}, fall.deep.perc = sum(df$DPr[last.irr.index:(dec.31.index)]))
      }
    }
  }
  #do.call(rbind, lapply(split(model.result, model.result$years), DeepPercCalc))
  
  #determine green water capture from leaf-drop to flowering
  GreenWaterCaptureCalc <- function(df) {
    if (cropname=='alfalfa.imperial') {
      return(data.frame(GW.capture.net = NA))
    } else if (df$doys.model[1] > Jharv | df$doys.model[nrow(df)] < Jmid) {
        return(data.frame(GW.capture.net = NA))
    } else {
        return(data.frame(GW.capture.net = df$Dr.end[which(df$doys.model == Jharv)] - df$Dr.end[which(df$doys.model == Jdev)]))
    }
  }
  #do.call(rbind, lapply(split(model.result, model.result$water.year), GreenWaterCaptureCalc))
  if (length(U2.df$DOY)==length(RHmin.df$DOY) & length(U2.df$DOY)==length(ETo.df$DOY)) {
    doys.model <- U2.df$DOY
    dates <- as.Date(U2.df$dates, format='%m_%d_%Y')
    days <- U2.df$day
    months <- U2.df$month
    years <- U2.df$year
    water.year <- years
    water.year[which(months >= 10)] <- years[which(months >= 10)] + 1
    print('Temporal coverages match in Spatial CIMIS.')
  } else {
    print('There are differing temporal coverages in the Spatial CIMIS data.')
  }
  model.length <- nrow(ETo.df)
  #get precip and spatial CIMIS data to same temporal end point
  last.date <- ETo.df$dates[nrow(ETo.df)]
  P.df <- P.df[1:which(P.df$dates==last.date), ]
  cropname.dir <-  paste0(cropname, '_majcomps')
  scenario.name <- if (cropname == 'grapes.wine') { #this was modified after the fact so that AD or RDI.min precedes root depth
    paste0(cropname.dir, '/scenario_', root_depth, as.character(RDI.min), 'RDI.min')} else {
        paste0(cropname.dir, '/scenario_', root_depth, as.character(AD.percentage), 'AD')
      } #need to add irr.type here
  if (cropname=='alfalfa.intermountain' | cropname=='alfalfa.CV' | cropname == 'alfalfa.imperial') {
    paw.var <- paste0('z', root_depth, '_cmH2O_unmodified_comp') 
  } else {
      paw.var <- paste0('z', root_depth, '_cmH2O_modified_comp')
  }
  if (!dir.exists(file.path(resultsDir, cropname.dir))) {
    dir.create(file.path(resultsDir, cropname.dir))
  }
  if (!dir.exists(file.path(resultsDir, scenario.name))) {
    dir.create(file.path(resultsDir, scenario.name))
  }
#initialization assumptions
  first_irrigation <- 0
  TEW.fraction <- 0.5
#limit model.scaffold to cropname
  if (cropname=='alfalfa.intermountain' | cropname=='alfalfa.CV' | cropname=='alfalfa.imperial') {
    model.scaffold.crop <- model.scaffold[which(model.scaffold$crop_code==cropcode & model.scaffold$alfalfa.zone==alfalfa.zone), ] #further refine model.scaffold.crop if cropname == alfalfa according to geography
  } else if (cropname=='grapes.table' | cropname=='grapes.wine') {
    model.scaffold.crop <- model.scaffold[grepl(grape.zone, model.scaffold$grape.zone), ] #further refine model.scaffold.crop if cropname == grape according to geography; wine.grapes will be run with 'Central California Foothills and Coastal Mountains'; grapes.table will be run for 'Central California Valley' and 'Sonora Basin and Range', the latter is forced to 'Central California Valley' after reading in model scaffold for simplicity's sake
  } else {
      model.scaffold.crop <- model.scaffold[which(model.scaffold$crop_code==cropcode), ] #80,401 unique crop, soil, and climate combinations for almond (spatially, this is the equivalent of only 7,236 ha)
  }
#make a results data.frame
  if (results_file == 'new') {
    model.length.yrs <- max(ETo.df$year) - min(ETo.df$year) + 1 #data starts 10/2003
    paw.vector <- model.scaffold.crop[[paw.var]]
  #print(head(paw.vector))
    #mukey	crop_code	PRISMcellnumber	CIMIScellnumber	unique_model_code	full_matrix_rownum	n_compkeys	cokey	compname	comppct_r	majcompflag	TEW	REW	surface.depth	z2.0m_cmH2O_modified_comp	Model.Year
    model.scaffold2 <- model.scaffold.crop[ ,c(-11:-22, -26:-27)] #takes out paw data, grape.zone, and alfalfa.zone from model scaffold for pasting results later
    model.scaffold2$paw <- paw.vector
    colnames(model.scaffold2)[14] <- paw.var
    model.scaffold.results <- model.scaffold2[rep(seq.int(1, nrow(model.scaffold2)), model.length.yrs), 1:ncol(model.scaffold2)] #makes a new data.frame with each row repeated model.length.yrs number of times
    model.scaffold.results <- model.scaffold.results[order(model.scaffold.results$unique_model_code, model.scaffold.results$cokey), ] #for some reason it is produced out of order
    model.scaffold.results <- data.frame(model.scaffold.results, Model.Year=rep(seq(from=min(ETo.df$year), to=max(ETo.df$year), by=1), times=nrow(model.scaffold.crop)), Irr.1=as.Date('1900-01-01'), Irr.Last=as.Date('1900-01-01'), RAW.end.season=NA, PAW.end.season=NA, Dr.end.season=NA, P.end.season=NA, Irr.end.storage=NA, GW.ET.growing=NA, Irr.app.total=NA, Irr.app.last=NA, ET.growing=NA, E.growing=NA, T.growing=NA, H2O.stress=NA, GW.ET.to.Irr1=NA, GW.E.to.Irr1=NA, GW.T.to.Irr1=NA, ET.annual=NA, E.annual=NA, T.annual=NA, deep.perc.annual=NA, winter.deep.perc=NA, post.Irr1.deep.perc=NA, fall.deep.perc=NA, GW.capture.net=NA)
    model.scaffold.results$unique_model_code_final <- paste0(as.character(model.scaffold.results$unique_model_code), as.character(model.scaffold.results$cokey)) #need to use as.character to preserve integrity of long integers
#model.scaffold.results[which(model.scaffold.results$unique_model_code==100058 & model.scaffold.results$cokey==13094564), ]
    rm(model.scaffold2)
  } else {
    setwd(file.path(resultsDir, scenario.name))
    #fname <- list.files(pattern = glob2rx('*_FAO56results.csv'))
    #model.scaffold.results <- read.csv(fname, stringsAsFactors = FALSE)
    model.scaffold.results <- read.csv(results_file, stringsAsFactors = FALSE)
    model.scaffold.results$Irr.1 <- as.Date(model.scaffold.results$Irr.1)
    model.scaffold.results$Irr.Last <- as.Date(model.scaffold.results$Irr.Last)
  }
  crop.parameters <- CropParametersDefine(crop.parameters.df, cropname)
  Kcb.std <- KcbDefine(doys.model, cropname) #this will be substituted with a crop code
  fc <- fcCalc(doys.model, crop.parameters, cropname) #TO-DO: implement alternative fc calculation in accordance with Eq. 11 from Allen et al. 2005: ((Kcb-Kcmin)/(Kcmax-Kcmin))^(1+0.5*h).  However, this produced a strange result in spreadsheet model for almonds, where increasing h decreases fc.
  if (cropname == 'alfalfa.intermountain' | cropname == 'alfalfa.CV' | cropname == 'alfalfa.imperial') {
    harvest.days <- HarvestDays(cropname)
  }
  if (cropname=='alfalfa.intermountain') {
    Jdev <- crop.parameters$Jini[crop.parameters$crop==cropname] #Jdev is intended to represent start of growing season as DOY
  } else if (cropname=='alfalfa.CV' | cropname == 'alfalfa.imperial') {
    Jdev <- 1 #day 1 is start of green alfalfa growth, though supressed by cold
  } else {
      Jdev <- crop.parameters$Jdev[crop.parameters$crop==cropname]
  }
  Jmid <- crop.parameters$Jmid[crop.parameters$crop==cropname]
  Jlate <- crop.parameters$Jlate[crop.parameters$crop==cropname]
  if (cropname=='alfalfa.intermountain') {
    Jharv <- harvest.days[length(harvest.days)] + crop.parameters$Lini.cycle[crop.parameters$crop==cropname] + crop.parameters$Ldev.cycle[crop.parameters$crop==cropname] + crop.parameters$Lmid.cycle[crop.parameters$crop==cropname] + crop.parameters$Lfrost.kill[crop.parameters$crop==cropname] #Jharv is intended to represent end of growing season where alfalfa is dormant in this region, approx. 11/23 that would be best varied as function of TMIN
  } else if (cropname == 'alfalfa.imperial' | cropname=='alfalfa.CV') {
    Jharv <- 365 #what about leap years?
  } else {
      Jharv <- crop.parameters$Jharv[crop.parameters$crop==cropname]
  }
  fw <- fwSelect(irrigation.parameters, irr.type)
  fewi <- fewiCalc(fc, fw)
  fewp <- fewpCalc(fc, fewi)
#loop through all rows of model scaffold but only do these operations once for each model.scaffold.crop row
  set.seed(461980)
  rows.to.sample <- sample(1:nrow(model.scaffold.crop), 0.005*nrow(model.scaffold.crop))
  if (nrow(model.scaffold.crop) > 10000) {
    save.times <- seq(from=10000, to=nrow(model.scaffold.crop), by=10000)
  } else {save.times <- 1000}
  for (n in row_start:nrow(model.scaffold.crop)) {
    model.code <- model.scaffold.crop$unique_model_code[n]
    PAW <- model.scaffold.crop[[paw.var]][n]*10
    AD <- (AD.percentage/100)*PAW
    cokey <- model.scaffold.crop$cokey[n]
    REW.parameter <- model.scaffold.crop$REW[n]
    TEW.parameter <- model.scaffold.crop$TEW[n]
    if (is.na(AD) | is.na(REW.parameter) | is.na(TEW.parameter)) {
      next(print(paste('Soils data is missing for cokey ', as.character(cokey)))) #TO-DO: write this result to separate file of NAs
    }
    if (AD==0 | TEW.parameter==0 | REW.parameter==0) {
      print(paste('AD, TEW, or REW is 0 for cokey ', as.character(cokey)))
      next
    }
    if (REW.parameter > TEW.parameter) { #there are several hundered instances where this was the result of the SSURGO aggregation work
      REW.parameter <- TEW.parameter
    }
  #identify climatic and soil parameters
    spCIMIScell <- model.scaffold.crop$CIMIScellnumber[n]
    PRISMcell <- model.scaffold.crop$PRISMcellnumber[n]
    comppctr <- model.scaffold.crop$comppct_r[n]
    mukey <- model.scaffold.crop$mukey[n]
    P <- P.df[ , which(colnames(P.df)==paste0('cell_', as.character(PRISMcell)))]
    ETo <- ETo.df[ , which(colnames(ETo.df)==paste0('cell_', as.character(spCIMIScell)))]
    U2 <- U2.df[ ,which(colnames(U2.df)==paste0('cell_', as.character(spCIMIScell)))]
    RHmin <- RHmin.df[ ,which(colnames(RHmin.df)==paste0('cell_', as.character(spCIMIScell)))]
    if (all(is.na(P)) | all(is.na(ETo))) {
      print(paste('Climate data is missing for scenario number', as.character(n))) #TO-DO: write this result to separate file of NAs
      next
    }
    Kcb.df <- KcbAdj(Kcb.std, crop.parameters=crop.parameters, cropname, U2, RHmin) #object 'crop.parameters' not found
    Kcb.adjusted <- Kcb.df$Kcb.climate.adj
    days.no.irr <- DaysNoIrr(P, ETo, Kcb.adjusted, AD, doys.model, years, Jmid, Jharv)
    if (is.null(days.no.irr)) {
      print(paste0('Null value returned for days.no.irr for model.scaffold.crop n = ', as.character(n), ', likely a result of a low PAW: ', as.character(PAW), ' mm.'))
      next
    }
    Kcmax <- Kcb.df$Kc.max
    Dei.initial <- numeric(length = model.length)
    DPei <- numeric(length = model.length)
    Dep.initial <- numeric(length = model.length)
    Kri <- numeric(length = model.length)
    Kei <- numeric(length = model.length)
    Ei <- numeric(length = model.length)
    Dei.end <- numeric(length = model.length)
    Krp <- numeric(length = model.length)
    Kep <- numeric(length = model.length)
    Ep <- numeric(length = model.length)
    DPep <- numeric(length = model.length)
    Dep.end <- numeric(length = model.length)
    W <- numeric(length = model.length)
    Kc.ns <- numeric(length = model.length)
    ETc.ns <- numeric(length = model.length)
    Dr.initial <- numeric(length = model.length)
#assumes depletion is zero for root zone except for half of TEW in upper layer
    Ir <- numeric(length = model.length)
    DPr <- numeric(length = model.length)
    Ks <- numeric(length = model.length)
    Kc.act <- numeric(length = model.length)
    Dr.end <- numeric(length = model.length)
    ETc.act <- numeric(length = model.length)
  #now for i=1, which will have slighly modified execution as it is the model initialization
    Dei.initial[1] <- max(TEW.parameter * TEW.fraction - P[1], 0) #this is an initial estimate of 
  #the water balance for the exposed surface soil where irrigation is not applied 
  #that assumes all daily precip occurs early in the morning so as to estimate Ke 
  #and Kr.  Same done for Dep.initial
    Dep.initial[1] <- max(TEW.parameter * TEW.fraction - P[1], 0)
    Kri[1] <- KrCalc(TEW.parameter, REW.parameter, Dei.initial, 1)
    Krp[1] <- KrCalc(TEW.parameter, REW.parameter, Dep.initial, 1)
    W[1] <- WCalc(TEW.parameter, Dei.initial, Dep.initial, fewp, fewi, 1)
    DPei[1] <- max(P[1] - TEW.parameter * TEW.fraction, 0) #initial estimate assumes irrigation is zero on previous day
    DPep[1] <- DPepCalc(P, Dep.initial, 1)
    Kei[1] <- min(Kri[1] * W[1] * (Kcmax[1] - Kcb.adjusted[1]), fewi[1] * Kcmax[1])
    Kep[1] <- min(Krp[1] * (1 - W[1]) * (Kcmax[1] - Kcb.adjusted[1]), fewp[1] * Kcmax[1])
    Ep[1] <- EpCalc(ETo, Kep, 1)
    Ei[1] <- EiCalc(ETo, Kei, 1)
    Dep.end[1] <- min(TEW.parameter, max(TEW.parameter * TEW.fraction - P[1] + Ep[1] / fewp[1] + DPep[1], 0)) #replaces Dep.end[i-1] with TEW.parameter * TEW.fraction
    Dei.end[1] <- min(TEW.parameter, max(TEW.parameter * TEW.fraction - P[1] + Ei[1] / fewi[1] + DPei[1], 0)) #replaces Dei.end[i-1] with Dei.intial[i]
    Kc.ns[1] <- KcnsCalc(Kcb.adjusted, Kei, Kep, 1)
    ETc.ns[1] <- ETcnsCalc(Kc.ns, ETo, 1)
    Dr.initial[1] <- max(TEW.parameter * TEW.fraction - P[1] + ETc.ns[1], 0) #initial calc
    Ks[1] <- KsCalc(Dr.initial=Dr.initial, PAW=PAW, i=1)
    Ir[1] <- IrCalc(Ks, RDI.min, AD, Dr.initial, doys.model, Jdev, Jharv, days.no.irr, 1)
    DPr[1] <- max(max(P[1] + Ir[1] - TEW.parameter * TEW.fraction - ETc.ns[1], 0)) #initial calc
    Kc.act[1] <- KcactCalc(Ks, Kcb.adjusted, Kei, Kep, 1)
    ETc.act[1] <- ETcactCalc(Kc.act, ETo, 1)
    Dr.end[1] <- TEW.parameter * TEW.fraction - P[1] + Kc.act[1] * ETo[1] + DPr[1] #initial calc
    for (i in 2:model.length) { #now for days 2...model.length after initialization
      Dei.initial[i] <- DeiInitialCalc(Dei.end, P, Ir, fw, i)
      Dep.initial[i] <- DepInitialCalc(Dep.end, P, i)
      Kri[i] <- KrCalc(TEW.parameter, REW.parameter, Dei.initial, i)
      Krp[i] <- KrCalc(TEW.parameter, REW.parameter, Dep.initial, i)
      W[i] <- WCalc(TEW.parameter, Dei.initial, Dep.initial, fewp, fewi, i)
      DPep[i] <- DPepCalc(P, Dep.initial, i)
      DPei[i] <- DPeiCalc(P, Ir, fw, Dei.initial, i)
      Kei[i] <- KeiCalc(Kri, W, Kcmax, Kcb.adjusted, fewi, TEW.parameter, Dei.end, DPei, P, ETo, i)
      Kep[i] <- KepCalc(Krp, W, Kcmax, Kcb.adjusted, fewp, TEW.parameter, Dep.end, DPep, P, ETo, i)
      Ep[i] <- EpCalc(ETo, Kep, i)
      Ei[i] <- EiCalc(ETo, Kei, i)
      Dep.end[i] <- DepEndCalc(Dep.end, P, Ep, fewp, DPep, i)
      Dei.end[i] <- DeiEndCalc(Dei.end, P, Ir, fw, fewi, Ei, DPei, i)
      Kc.ns[i] <- KcnsCalc(Kcb.adjusted, Kei, Kep, i)
      ETc.ns[i] <- ETcnsCalc(Kc.ns, ETo, i)
      Dr.initial[i] <- DrInitialCalc(Dr.end, ETc.ns, P, Ir, i)
      Ks[i] <- KsCalc(Dr.initial=Dr.initial, PAW=PAW, i=i) #corrected on 10/18/17
      Ir[i] <- IrCalc(Ks, RDI.min, AD, Dr.initial, doys.model, Jdev, Jharv, days.no.irr, i)
      DPr[i] <- DPrCalc(P, Ir, ETc.ns, Dr.end, i)
      Kc.act[i] <- KcactCalc(Ks, Kcb.adjusted, Kei, Kep, i)
      Dr.end[i] <- DrEndCalc(Dr.end, P, Ir, Kc.act, ETo, DPr, i)
      ETc.act[i] <- ETcactCalc(Kc.act, ETo, i) #could take this out of loop
    }
    model.result <- data.frame(dates, months, days, years, water.year, doys.model, P, ETo, RHmin, U2, lapply(X=list(Kcb.std=Kcb.std, Kcb.adjusted=Kcb.adjusted, Kcmax=Kcmax, fceff=fc, fw=fw, fewi=fewi, fewp=fewp, Dei.initial=Dei.initial, Dep.initial=Dep.initial, Kri=Kri, Krp=Krp, W=W, Kei=Kei, Kep=Kep, Ei=Ei, Ep=Ep, Dpei=DPei, DPep=DPep, Dei.end=Dei.end, Dep.end=Dep.end, Kc.ns=Kc.ns, ETc.ns=ETc.ns, Dr.initial=Dr.initial, Ir=Ir, DPr=DPr, Ks=Ks, Kc.act=Kc.act, ETc.act=ETc.act, Dr.end=Dr.end), round, digits=rounding_digits))
    model.scaffold.results[which(model.scaffold.results$unique_model_code==model.code & model.scaffold.results$cokey == cokey), which(colnames(model.scaffold.results)=='Irr.1'):(which(colnames(model.scaffold.results)=='GW.capture.net'))] <- merge(cbind(do.call(rbind, lapply(split(model.result, model.result$years), IrDateCalc)), do.call(rbind, lapply(split(model.result, model.result$years), WaterBalanceCalc)), do.call(rbind, lapply(split(model.result, model.result$years), GreenWaterIrr1Calc)), do.call(rbind, lapply(split(model.result, model.result$years), DeepPercCalc))), do.call(rbind, lapply(split(model.result, model.result$water.year), GreenWaterCaptureCalc)), by="row.names", all=TRUE)[ ,2:26]
    print(paste(scenario.name, as.character(n)))
    if (n==1 | n %in% rows.to.sample) {
      setwd(file.path(resultsDir, scenario.name))
      if (cropname=='grapes.wine') {
        write.csv(model.result, paste0(cropname, root_depth, 'RDI.min', as.character(RDI.min), '_', as.character(model.code), '_', as.character(cokey), '_', Sys.Date(), '.csv'), row.names=FALSE)
      } else {
      write.csv(model.result, paste0(cropname, root_depth, 'AD', as.character(AD.percentage), '_', as.character(model.code), '_', as.character(cokey), '_', Sys.Date(), '.csv'), row.names=FALSE)
      }
    }
    if (n %in% save.times) { #was (n==100 | n %in% save.times)
      setwd(file.path(resultsDir, scenario.name))
      if (cropname=='grapes.wine') {
        write.csv(model.scaffold.results, paste0(cropname, root_depth, 'RDI.min', as.character(RDI.min), '_FAO56results.csv'), row.names=FALSE)
      } else {
        write.csv(model.scaffold.results, paste0(cropname, root_depth, 'AD', as.character(AD.percentage), '_FAO56results.csv'), row.names=FALSE)
      }
    } else {next}
  }
  setwd(file.path(resultsDir, scenario.name))
  if (cropname=='grapes.wine') {
    write.csv(model.scaffold.results, paste0(cropname, root_depth, 'RDI.min', as.character(RDI.min), '_FAO56results.csv'), row.names=FALSE)
  } else {
      write.csv(model.scaffold.results, paste0(cropname, root_depth, 'AD', as.character(AD.percentage), '_FAO56results.csv'), row.names=FALSE)
  }
  metadata <- cbind(data.frame(date.run=Sys.Date(), crop=cropname, cropscape.code=cropcode, AD.percentage=AD.percentage, rooting.depth=root_depth, irrigation.type=irr.type, paw.varname = paw.var, model.days=model.length, first.day=dates[1], last.day=dates[length(dates)], n.models=nrow(model.scaffold.crop)), crop.parameters[which(crop.parameters$crop==cropname), 2:ncol(crop.parameters)], irrigation.parameters[which(irrigation.parameters$irrigation.type==irr.type), 'fw'])
  colnames(metadata)[ncol(metadata)] <- 'fw'
  if (cropname == 'grapes.wine') {
    write.csv(metadata, paste0(cropname, root_depth, 'RDI.min', as.character(RDI.min), '_model_metadata.csv'), row.names = FALSE)
  } else {
      write.csv(metadata, paste0(cropname, root_depth, 'AD', as.character(AD.percentage), '_model_metadata.csv'), row.names = FALSE)
  }
}


#legend for FAO56 abbreviations
#De,j=De,j-1 - P,j - Ij/fw + Ej/fewi + DPei,j (again, ignoring tranpiration from upper 10 cm and runoff, eqn. 21) 
#pseudo-code outline of 'separate prediction of evaporation from soil wetted by precipitation only' following Allen et al. 2005, except ignoring runoff, essentially assuming that runoff will really only occur when soils are near field capacity, so partitioning this as 'deep percolation' is acceptable and is consciously preferred over introduced errors from the curve number approach
#Ke = Kei + Kep (eqn. 14)
#where Kei=evaporation coefficient for the exposed fraction of the soil wetted by both irrigation and by precipitation and Kep=evaporation coefficient for the exposed fraction of the soil wetted by precipitation only
#Kei = Kri*W*(Kcmax-Kcb) <= fewi*Kcmax (eqn. 15)
#Kep = Krp*(1-W)*(Kcmax-Kcb) <= fewp*Kcmax (eqn. 16)
#where fewi=fraction of soil wetted by both irrigation and precipitation and is exposed to rapid drying due to exposure to solar radiation and/or ventilation and is calculated as: min(1-fc, fw) (eqn. 18);
#where fewp=fraction of soil exposed to rapid drying and is wetted by precipitation only and is calculated as: 1-fc-fewi (eqn. 17)
#where Kri and Krp=evaporation reduction coefficients for the fewi and fewp fractions, respectively
#where W=weighting coefficient for partitioning the energy available for evaporation in the fewi and fewp soil fractions, depending on water availability, and is calculated as: 1/(1+(fewp/fewi)*(TEW-Dep)/(TEW-De)) (eqn. 19);
#where De=cumulative depletion depth (mm) from the evaporating layer for the fewi fraction of soil
#where Dep=cumulative depletion depth (mm) from the evaporating layer for the fewp fraction of soil
#finally, the water balance formation for fraction of soil only wetted by precipitation:
#Dep,j=Dep,j-1 - Pj + Ep,j/fewp + DPep,j (ignoring transpiration from upper 10 cm and runoff, eqn. 20)
#where Dep,j-1 and Dep,j are the cumulative depletion depths at the ends of days j-1 and j in the fewp fraction of the surface (mm)
#where Ep,j=evaporation from kewp fraction on day j and is calculated as: Kep*ETo in mm
#where DPep,j=deep percolation from the fewp fraction of the evaporation layer on day j if soil water content exceeds field capacity
#where 0 <= Dep, j <= TEW
#TEW=total evaporable water in upper layer and is calculated as: 1000*(field_capacity-0.5*wilting_point)*Ze
#where Ze=effective depth of wetting and field capacity and wilting point is assumed to be moisture content at 1/3 and 15 bars, respectively, determined from the SSURGO database
#and the water balance formulation for the fraction of soil wetted by both precipitation and irrigation
#De,j=De,j-1 - P,j - Ij/fw + Ej/fewi + DPei,j (again, ignoring tranpiration from upper 10 cm and runoff, eqn. 21)
#where 0 <= De,j <= TEW
#where Kri=(TEW-De,j-1)/(TEW-REW) (eqn. 22)
#where Krp=(TEW-Dep,j-1)/(TEW-REW) (eqn. 23)
#where Dpei,j=Pj + Ij/fw - Dei,j-1 >= 0
#where DPep,j=Pj-Dep,j-1 >= 0
#Kr depends on daily water balance of surface layer (upper 10-15 cm); Few will also depend on daily water balance if separate calculations are done for soil wetted by precipitation vs. by both precip and irrigation (see "Separate Prediction of Evaporation from Soil Wetted by Precipitation Only" in Allen et al. 2005)