Utilities

UPDATE A WORLDFILE FOR COMPATIBILITY WITH RHESSYS VERSION 5.14

Rhessys versions 5.14 + require an additional state variable for the fraction of litter cover. To add this state variable line to an existing worldfile so you can use it with the newer version, you must apply the following awk script. Create a text file with the following commands, call it add_litfrac_lines.awk:

{a=0;}

($2 == "snowpack_energy_deficit") {printf("%f %s\n%f %s\n",$1,$2,1.0,"litter.cover_fraction"); a=1;}

(a == 0) {printf("%s %s\n",$1,$2);}

To apply to an existing worldfile and create a new worldfile compatible with rhessys5.14, type the following command:

awk -f add_litfrac_lines.awk < worldfile_to_change > new_worldfile_name

BASIN DAILY WATER BALANCE FUNCTION

PASTE THE COMMANDS BELOW IN R TO CREATE THE WATER BALANCE FUNCTION

watbal = function(rb) {

rb$watbal.tmp=with(rb,precip-streamflow-trans-evap)

rb$sd=with(rb,sat_def-rz_storage-unsat_stor)

tmp=diff(rb$sd)

tmp=c(0,tmp)

rb$sddiff=tmp

tmp=diff(rb$snowpack)

tmp=c(0,tmp)

rb$snodiff=tmp

tmp=diff(rb$detention_store)

tmp=c(0,tmp)

rb$detdiff=tmp

tmp=diff(rb$litter_store)

tmp=c(0,tmp)

rb$litdiff=tmp

tmp=diff(rb$canopy_store)

tmp=c(0,tmp)

rb$candiff=tmp

tmp=diff(rb$gw.storage)

tmp=c(0,tmp)

rb$gwdiff=tmp

rb$watbal=with(rb,watbal.tmp+sddiff-snodiff-detdiff-litdiff-candiff-gwdiff)

rb

}

To use this function in an R data set:

If your basin daily output is called (for example) "bd" type "bd=watbal(bd)".

This will add a number of fields to your basin daily output file, the last of which will be called "watbal". This is your water balance error. You should see a minor numerical error here even if your water balances (on the order of 10^-6). If your watbal values are negative then water inputs are less than water outputs, and vice versa.

PATCH DAILY WATER BALANCE FUNCTIONS

NOTES:

(1) Before you run the water balance script below, you must have added a "rain" column to your patch daily output data that contains the patch-specific precipitation. Be careful to add the right precip values if your worldfile uses multiple base stations, isohyets, or a precip lapse rate.

(2) These scripts will not work if you have multiple stratums.

TO RUN:

(1) In order to run this routine, you need to have brought into R a patch daily output file (pd) and have added a rain field and also have brought in a canopy (stratum) daily output file (sd).

(2) Copy and paste into R the water balance and multipatch run functions below (i.e., all the R commands below)

(3) If you only have a single patch and a single canopy, you can run the watbal function as-is. For example, if your patch output is called "pd" and your canopy stratum output is called "sd" then type:

rp=pwatbal(pd,sd)

This will add a number of fields to your patch output file, the last of which will be called "watbal". This is your water balance error. You should see a minor numerical error here even if your water balances (on the order of 10^-6). If your watbal values are negative then water inputs are less than water outputs, and vice versa.

(4) If you have multiple patches with only a single stratum per patch, you need to run a different script to loop through patches individually, call the water balance function, then write the output. For example, if your patch output is called "pd" and your canopy stratum output is called "sd", and you want the water balance output written to a table called "out", then you would type:

out=pwatbalmult(rp,rc)

This creates a new output table, called "out" in this example, that contains date fields and a single column for each patch with the water balance error for that patch (called "P345","P578","P900", etc. where the field label number matches the patch ID). You might see a minor numerical error here even if your water balances (on the order of 10^-6). If your watbal values are negative then water inputs are less than water outputs and vice versa.

Note that this may take a long time to run if you have many patches; better to run on a single patch or a single hillslope than a full basin.

PASTE THE FOLLOWING COMMANDS INTO R FOR THE PATCH WATER BAL FUNCTION

pwatbal = function(qp,qc) {

qp$watbal.tmp=with(qp,rain+Qin-Qout-trans_sat-trans_unsat-evap-evap_surface-soil_evap)

qp$sd=with(qp,sat_def-rz_storage-unsat_stor)

tmp=diff(qp$sd)

tmp=c(0,tmp)

qp$sddiff=tmp

tmp=diff(qp$snow)

tmp=c(0,tmp)

qp$snodiff=tmp

tmp=diff(qp$detention_store)

tmp=c(0,tmp)

qp$detdiff=tmp

tmp=diff(qp$litter.rain_stor)

tmp=c(0,tmp)

qp$litdiff=tmp

tmp=diff(qc$rain_stored+qc$snow_stored)

tmp=c(0,tmp)

qp$candiff=tmp

qp$watbal=with(qp,watbal.tmp+sddiff-snodiff-detdiff-litdiff)

qp

}

PASTE THE FOLLOWING COMMANDS INTO R FOR THE MULTIPATCH RUN FUNCTION

pwatbalmult = function(rp,rc) {

pids=unique(rp$patchID)

tmp=subset(rp,rp$patchID==pids[1])

wbp=tmp[c("day","month","year")]

wbp=mkdate(wbp)

n=ncol(wbp)+1

for (i in pids) {

tmpp=subset(rp,rp$patchID==i)

tmpc=subset(rc,rc$patchID==i)

tmpp=pwatbal(tmpp,tmpc)

wbp[,n]=tmpp$watbal

names(wbp)[c(n)]=paste("P",i,sep="")

n=n+1

}

wbp

}

RESAMPLE EXISTING RHESSys METEOROLOGICAL DATA

1. read RHESSys climate data (you can use the read_rhessys_met script) into an R data file, for this example we call this new R object "clim"
2. you will also need to run the R script mkdate on the R object "clim" so all of the needed columns are included in the R object

clim = read_rhessys_met("../../clim/sage")

clim=mkdate(clim)

1. note: you must have complete calendar years or water years (depending on which you will use) for selection
2. if years are not complete - get rid of incomplete years (you can use subset, i.e.: clim = subset(clim, clim$wy > 1953), this would get rid of first year 1953 if it didn't include a full wateryear
3. this may occur at beginning and end of the record

orig = unique(clim$wy)

nyrs = length(orig)

1. determine when leap years are

dyrs = aggregate(clim$wyd, by=list(clim$wy), max)

colnames(dyrs) = c("wyd","nday")

leaps = subset(dyrs, dyrs$nday>365)

1. initialize a vector called "new" of years that you want climate for
2. you can make this however you want
3. here we show how to randomly sample from original years

new = sample(orig, replace=T, size=nyrs)

1. another possibility you could repeat the same year n times
2. my new clim sequence was a data frame, so i adjust the for loop below

new = rep(2000, times=nyrs)

1. or repeat the 5 wettest years

clim.wy = aggregate(clim, by=list(clim$wy), mean)

idx = order(clim.wy$rain)

tmp = clim.wy$wy[idx[1:5]]

new = rep(tmp, times=ceiling(nyrs/5))

1. initialize clim.new which will be your new climate sequence

clim.new=clim

clim.new[,]=0

clim$place = seq(from=1, to=length(clim$year))

startr=1

for (i in 1:(length(dyrs$wy)) ) {

cyr = new[i] # checks what yr for var seq(?)

tmpn = subset(clim,clim$wy==cyr) # extracts that var/new yr f/orig clim

lnn = nrow(tmpn) # n days in extracted yr (leap?)

oyr = orig[i] # checks what year for orig seq

tmpo = subset(clim, clim$wy==oyr) # extracts that "in seq" yr f/orig clim

lno = nrow(tmpo) # n days in orig yr (leap?)

if (lno > lnn) # compare nday in orig clim to extr clim

tmpn = rbind(tmpn, tmpo[lnn,]) # add last day (9/30) of old to new if leap

if (lnn > lno)

tmpn = tmpn[1:lno,]

endr=startr+lno-1 # which row to insert year in clim.new

clim.new[startr:endr,] = tmpn # insert extracted year into clim.new

startr=endr+1 # tally row count by number of days

done = i # useful for debugging

}

1. now lets write out the met variables in RHESSys format
2. this could be extended if there are optional met inputs
3. we use dates from the ORIGINAL clim sequence here since we want the start date of the artificial sequence to be the same

header = sprintf("%d %d %d %d", clim$year[1], clim$month[1], clim$day[1], 1)

write(header, "../clim/new.tmax")

write.table(clim.new$tmax, file="../clim/new.tmax", append=T, col.names=F, row.names=F, quote=F)

write(header, "../clim/new.tmin")

write.table(clim.new$tmin, file="../clim/new.tmin", append=T, col.names=F, row.names=F, quote=F)

1. mkdate changes precip to mm so change back by /1000

write(header, "../clim/new.rain")

write.table(clim.new$rain/1000, file="../clim/new.rain", append=T, col.names=F, row.names=F, quote=F)