microclimate-model-tutorial.Rmd

---
title: "Introduction to the NicheMapR microclimate model"
author: "Michael Kearney"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Microclimate Model Tutorial}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, echo = FALSE}
knitr::opts_chunk$set(
 eval = TRUE
)
```
This vignette provides a detailed tutorial of the Niche Mapper microclimate model as implemented in the package ```NicheMapR```. See the [microclimate-model-theory-equations](microclimate-model-theory-equations.Rmd) vignette for details on the underlying theory and equations.

## Getting started: working with the micro_global() function

The ```NicheMapR``` package includes data and functions for using a global monthly climate database for computing microclimates. This tutorial illustrates the different capabilities of the microclimate model using this database. The vignettes [microclimate-IO](microclimate-IO.Rmd) and [microclimate-hourly-example](microclimate-hourly-example.Rmd) illustrate how to customise the model to work with other sources of data.

### Basic operation: modelling microclimates for the average day of each month

This first example involves the most basic case of running the model to produce 12 days of output, one for each month of the year, without invoking the soil moisture or snow subroutines. In this mode of operation, the model does three iterations of each day, starting with a uniform soil temperature profile and the first two simuulations acting as 'burn-in' runs prior to the third and final run which is provided as the final output.

Try running the model for a place of interest to you, as follows:
```{r, message=FALSE, warnings=FALSE}
library(NicheMapR) 
longlat <- c(-89.40123, 43.07305) # Madison, Wisconsin, USA
micro <- micro_global(loc = longlat)
```

The results are stored as a list in the variable ```micro```. The main ones to focus on for now are

* ```metout``` The above ground micrometeorological conditions under the minimum specified shade
* ```shadmet``` The above ground micrometeorological conditions under the maximum specified shade
* ```soil``` Hourly predictions of the soil temperatures under the minimum specified shade
* ```shadsoil``` Hourly predictions of the soil temperatures under the maximum specified shade

The first two rows of the output tables ```metout``` and ```shadmet``` look like this

```{r, echo=FALSE, results='asis', message=FALSE, warnings=FALSE}
knitr::kable(head(micro$metout[, 1:9], 2), digits = 2)
knitr::kable(head(micro$metout[, 10:18], 2), digits = 2)
```

and show, for each day of the year and for each hour of the day, a series of aboveground microclimatic conditions as follows:

1. DOY - day of year
1. TIME - time of day (mins)
1. TALOC - air temperature (&deg;C) at local height (specified by 'Usrhyt' variable)
1. TAREF - air temperature (&deg;C) at reference height (1.2m)
1. RHLOC - relative humidity (%) at local height (specified by 'Usrhyt' variable)
1. RH  - relative humidity (%) at reference height (1.2m)
1. VLOC - wind speed (m/s) at local height (specified by 'Usrhyt' variable)
1. VREF - wind speed (m/s) at reference height (1.2m)
1. SNOWMELT - snowmelt (mm)
1. POOLDEP - water pooling on surface (mm)
1. PCTWET - soil surface wetness (%)
1. ZEN - zenith angle of sun (degrees - 90&deg; = below the horizon)
1. SOLR - solar radiation (W/m2)
1. TSKYC - sky radiant temperature (&deg;C)
1. DEW - dew formation (mm/h)
1. FROST - frost formation (mm/h)
1. SNOWFALL - snow predicted to have fallen (mm)
1. SNOWDEP - predicted snow depth (cm)

Note that for air temperature, relative humidity and wind speed, there's a 'reference height' value and a 'local height' value. The reference height value is the height at which the meteorological observations were made that were given to the model as input. The local height values are the microclimatic values at a height determined by the variable ```Usrhyt```, which you can set to be relevant to the organism you are interested in (by default it is 1cm). It should be at the midpoint of the animal or plant's height.

Check what you got with:
```
metout<-micro$metout # put the metout result into a variable 'metout'
head(metout,24) # show the first 24 rows, i.e. the first day

shadmet<-micro$shadmet # put the shadmet result into a variable 'shadmet'
head(shadmet,24) # show the first 24 rows, i.e. the first day
```

Because the snow model wasn't flagged to run, there will be zeros in all the snow-related columns.

The first two rows of the output tables ```soil``` and ```shadsoil``` look like this 

```{r, echo=FALSE, results='asis', message=FALSE, warnings=FALSE}
knitr::kable(head(micro$soil[, 1:12], 2), digits = 2)
```

and show, for each day of the year and for each hour of the day, the soil temperatures as follows:

1. DOY - day of year
1. TIME - time of day (mins)
1. Columns 3-12 D0cm ... - soil temperatures at each of the 10 specified depths

Have a look at yours:
```
soil<-micro$soil # put the soil result into a variable 'soil'
head(soil,24) # show the first 24 rows, i.e. the first day

shadsoil<-micro$shadsoil # put the metout result into a variable 'shadsoil'
head(shadsoil,24) # show the first 24 rows, i.e. the first day
```
The 10 different depths used by the model are specified with ```DEP``` and by default are `DEP=c(0., 2.5,  5.,  10.,  15.,  20.,  30.,  50.,  100.,  200.)`.

There are some example plots of all the variables in the help for the ```micro_global``` function. Here is a plot for all the soil temperatures in Madison, Wisconsin in July and December:

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
require(lattice, quietly = TRUE)
soil <- as.data.frame(micro$soil) # get the soil data
minshade <- micro$minshade[1] # get the value for minimum shade
with(subset(soil, DOY==196 | DOY==349), {xyplot(D0cm + D2.5cm + D5cm + D10cm + D15cm + D20cm + D30cm + D50cm + D100cm + D200cm ~ TIME | as.factor(DOY), ylim = c(-20, 70), xlab = "Time of Day (min)", ylab = "Soil Temperature (deg C)", auto.key = list(columns = 5), as.table = TRUE, type = "b", main = paste(minshade,"% shade"))})
```

Notice how the deeper one goes into the soil, the lower the amplitude of the fluctuations but also that the maximum temperature is reached later in the day. The daily fluctuations have largely disappeared by 50cm and at 200cm the temperature is stable year round. In fact the model assumes the deep soil temperature is the annual mean of the air temperatures, and this is a 'boundary condition' in the solution of the soil heat budget.

### Simulating shading by vegetation

By default the model runs twice, for two different levels of shading by vegetation. The particular shade levels can be selected with the parameters ```maxshade``` and ```minshade```. Also, the second shade level simulation can be turned off if only one shade level is required with the parameter ```runshade```. The code below runs the model for one shade level only, 50%:

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
micro <- micro_global(loc = longlat, runshade = 0, minshade = 50)
soil <- as.data.frame(micro$soil) # get the soil data
minshade<-micro$minshade[1] # get the value for minimum shade
with(subset(soil, DOY==196 | DOY==349),{xyplot(D0cm + D2.5cm + D5cm + D10cm + D15cm + D20cm + D30cm + D50cm + D100cm + D200cm ~ TIME | as.factor(DOY), ylim = c(-20, 70), xlab = "Time of Day (min)", ylab = "Soil Temperature (deg C)", auto.key=list(columns = 5), as.table = TRUE, type = "b", main = paste(minshade, "% shade"))})
```

Notice how, the effect of shade is to not only reduce the daytime soil temperatures, but also to increase the nighttime temperatures. This is due to the increased levels of infrared radiation reaching the ground surface because the shading vegetation is assumed to be at the reference height air temperature while the 'sky temperature' is typically much colder, depending on the level of cloud cover (have a look at the variable ```Tsky``` in the ```metout``` and ```shadmet``` outputs).

### Simulating terrain effects - slope, aspect, hillshade

Slope and aspect have strong effects on microclimates through the effects of direct solar radiation. The default setting is flat ground in the ```NicheMapR``` microclimate model. Here is the effect of a 45 degree southerly slope on the 0% shade simulation:

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
micro <- micro_global(loc = longlat, runshade = 0, minshade = 0, slope = 45, aspect = 180)
soil <- as.data.frame(micro$soil) # get the soil data
minshade <- micro$minshade[1] # get the value for minimum shade
with(subset(soil, DOY==196 | DOY==349), {xyplot(D0cm + D2.5cm + D5cm + D10cm + D15cm + D20cm + D30cm + D50cm + D100cm + D200cm ~ TIME | as.factor(DOY), ylim = c(-20, 70), xlab = "Time of Day (min)", ylab = "Soil Temperature (deg C)", auto.key=list(columns = 5), as.table = TRUE, type = "b", main = paste(minshade, "% shade, 45 degree slope, 180 degrees aspect"))})
```

The strongest effect is on the December soil temperatures.

Hillshade also has a strong effect on microclimates. Gullies and gorges receive a shorter window of direct solar radiation and have a smaller 'view' of the sky, and so don't cool down as much at night. Horizon angles can be specified to simulated these kinds of effects by supplying in the parameter ```hori``` a vector of 24 horizon angles (angle to nearest hill/cliff) starting due north and continuing clockwise in 15 degree intervals. Here is an example for a steep-sided gully running north-south:

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
micro <- micro_global(loc = longlat, runshade = 0, minshade = 0, hori = c(0, 0, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65, 0, 0, 65, 65, 65, 65, 65, 65, 65, 65, 65, 65))
soil <- as.data.frame(micro$soil) # get the soil data
minshade <- micro$minshade[1] # get the value for minimum shade
with(subset(soil,DOY==196 | DOY==349),{xyplot(D0cm + D2.5cm + D5cm + D10cm + D15cm + D20cm + D30cm + D50cm + D100cm + D200cm ~ TIME | as.factor(DOY), ylim = c(-20, 70), xlab = "Time of Day (min)", ylab = "Soil Temperature (deg C)", auto.key = list(columns = 5), as.table = TRUE, type = "b", main=paste(minshade,"% shade, north-south gully"))})
```

### Changing the soil properties

Another important influence on microclimatic conditions is the substrate type. The substrate solar reflectivity can be changed with the parameter ```REFL`` and the emissivity with the parameter ```SLE``` (both range from 0-1).

Soil thermal properties are also importance. Soil comprises a mixture of minerals, air and water. The density (Mg/m3), thermal conductivity (W/mK) and specific heat capacity (J/kg-K) of the mineral component are specified by the parameters ```Density```, ```Thcond``` and ```SpecHeat```, respectively. The bulk density determines the proportion of the soil volume that is air, parameter ```BulkDensity``` (Mg/m3). The overall soil conductivity, specific heat and density values are then adjusted for bulk density and soil moisture (see below), as well as for soil temperature, using themethods described in Campbell et al. (1994, equations 8, 9) and Campbell and Norman (Campbell & Norman 1998, equations 8.13, 8.17, 8.20 and data in Tables 8.2 and 9.1). All these properties can be specified separately for each soil node in the model. 

The default values are typical for soil minerals as described in Campbell and Norman (1998) Table 8.2, but  an 'A-horizon' layer for the first node, reflecting properties of organic matter, is also simulated (see Kearney et al. 2014). This can be turned off by specifying `cap = 0`. 

For simulating rock, ```BulkDensity``` should be set the same as ```Density``` and this can be easily set in ```micro_global``` by setting ```soiltype``` to value 0 (which also sets cap = 0):

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
micro <- micro_global(runshade = 0, minshade = 0, soiltype = 0)
soil <- as.data.frame(micro$soil) # get the soil data
minshade <- micro$minshade[1] # get the value for minimum shade
with(subset(soil, DOY==196 | DOY==349), {xyplot(D0cm + D2.5cm + D5cm + D10cm + D15cm + D20cm + D30cm + D50cm + D100cm + D200cm ~ TIME | as.factor(DOY), ylim = c(-20, 70), xlab = "Time of Day (min)", ylab = "Soil Temperature (deg C)",  auto.key = list(columns = 5), as.table = TRUE, type = "b", main = paste(minshade, "% shade, rock substrate"))})
```

Notice how the heat penetrates more deeply into the rock, producing a gentler temperature gradient with depth, cooler near-surface temperatures during the day and warmer near-surface temperatures at night. 

The texture of the soil, i.e. the proportion of silt, sand and clay and associated particle size distributions, has a powerful influence on soil moisture and is discussed further in the next section.

### Soil moisture

The microclimate model has the option of either specifing soil moisture or modelling it from first principles. The ```micro_global``` function by default monthly soil moisture estimates obtained from the Climate Prediction Center(NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, <http://www.esrl.noaa.gov/psd/>).

The model will compute soil moisture when the ```runmoist``` parameter = 1. This invokes an implementation of Campbell’s (1985) Program 11.1 for simulating depth-specific soil moisture, water potential and humidity gradients in the presence of vegetation. It uses matric potential as the dependent variable (rather than matric flux potential) and as a consequence is most computationally efficient when soil is relatively dry (it may significantly extend run times in wet environments). 

The soil moisture calculations are made at the 10 user-specified nodes (parameter ```DEP```) and also at the midpoints between them, but are only reported for the 10 nodes in ```DEP```. The user input variables specific to this module include, for each of these 19 nodes, the Campbell (or Clapp and Hornberger) b exponent (```BB```), the air entry potential ```PE``` (J/kg), the hydraulic conductivity ```KS``` (kg s m-3), the bulk density ```BD``` (Mg m-3), the density ```DD``` (Mg m-3) the root density ```L`` (m m-3) and the leaf area index LAI. All other model parameters are fixed at those suggested by Campbell (1985). The model has been tested against detailed observations across Australia with favourable results (Kearney et al. in prep).

Try running the model for your site with the soil moisture model turned on:
```{r, message=FALSE, warnings=FALSE}
micro <- micro_global(loc = longlat, runmoist = 1)
```

The output tables ```soilmoist``` and ```shadmoist```, which previously were zero, will now have soil moisture values (m3/m3) in them: 
```{r, message=FALSE, warnings=FALSE, echo=FALSE, results='asis'}
knitr::kable(head(micro$soilmoist[, 1:9], 2))
knitr::kable(head(micro$shadmoist[, 10:12], 2))
```

The model also gives you the soil water potential (the water-equivalent to temperature and the driving force for water flow, e.g. into or out of water-permeable eggs or skin) in the outputs ```soilpot``` and ```shadpot```:
```{r, message=FALSE, warnings=FALSE, echo=FALSE, results='asis'}
knitr::kable(head(micro$soilpot[, 1:9], 2))
knitr::kable(head(micro$shadpot[, 10:12], 2))
```

Finally, the soil humidity is also computed (outputs ```humid``` and ```shadhumid```):
```{r, message=FALSE, warnings=FALSE, echo=FALSE, results='asis'}
knitr::kable(head(micro$humid[, 1:9], 2))
knitr::kable(head(micro$shadhumid[, 10:12], 2))
```

Note that soil humidity is typically close to saturation, even in 'dry' soils. The situation in a burrow may lower than this but still much higher than aboveground conditions (see Schmidt-Nielsen & Schmidt-Nielsen. 1950).

Plotting the results for soil moisture through time across the year:
```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
soilmoist <- as.data.frame(micro$soilmoist) # get the minimum shade soil moisture output
minshade <- micro$minshade[1] # get the value for minimum shade
# append dates
days <- rep(seq(1, 12), 24)
days <- days[order(days)]
dates <- days + soilmoist$TIME / 60 / 24 - 1 # dates for hourly output
soilmoist<-cbind(dates, soilmoist)
for(i in 1:10){
  if(i == 1){
    plot(soilmoist[, i + 3] ~ soilmoist[, 1], ylim = c(0, 0.5), xlab = "Month",  ylab = "Soil Moisture (% vol)", col = i, type = "l", main = paste("soil moisture, ", minshade, "% shade"))
  }else{
    points(soilmoist[, i + 3] ~ soilmoist[, 1], col = i, type = "l")
  }
}
```

By default the model uses ```soiltype``` = 4, where Rock = 0, sand = 1, loamy sand = 2, sandy loam = 3, loam = 4, silt loam = 5, sandy clay loam = 6, clay loam = 7, silt clay loam = 8, sandy clay = 9, silty clay = 10, clay = 11, based on Campbell and Norman 1990 Table 9.1. One can also manually set the soil properites by setting ```soiltype = 12```, upon which the model will use the user-specified values of ```PE```, ```KS```, ```BB```, ```DD``` and ```DD```.

Rerunning with clay as the soil type:
```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
micro<-micro_global(loc = longlat, runmoist = 1, soiltype = 11)
soilmoist<-as.data.frame(micro$soilmoist) # get the minimum shade soil moisture output
minshade<-micro$minshade[1] # get the value for minimum shade
# append dates
days<-rep(seq(1,12),24)
days<-days[order(days)]
dates<-days+soilmoist$TIME/60/24-1 # dates for hourly output
soilmoist<-cbind(dates,soilmoist)
for(i in 1:10){
  if(i==1){
    plot(soilmoist[,i+3]~soilmoist[,1], ylim=c(0,0.5),xlab = "Month", ylab = "Soil Moisture (%
      vol)",col=i,type = "l",main=paste("soil moisture, ",minshade,"% shade"))
  }else{
    points(soilmoist[,i+3]~soilmoist[,1] ,col=i,type = "l")
  }
}
```

we see strong changes in the soil moisture values, with the fine-grained clay particles holding on more strongly to the water resulting in higher soil water content. Also very apparent in these plots is that the model was not reaching steady state soil moisture values with only 3 iterations of each day. This can be resolved by running the model for longer periods, as discussed in the next section.

### Running at finer time intervals and for longer time periods

It is  possible to run the model for multiple years by changing the ```nyears``` parameter. This can be useful when running the soil moisture and snow subroutines to allow annual cycles to stabilise. To illustrate, if we run the last example in the soil moisture section for five years, the soil moisture at all but the deepest levels have started to reach steady state:

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
nyears<-5
micro<-micro_global(loc = longlat, runmoist = 1, soiltype = 11, nyears = nyears)
soilmoist<-as.data.frame(micro$soilmoist) # get the minimum shade soil moisture output
minshade<-micro$minshade[1] # get the value for minimum shade
# append dates
days<-rep(seq(1,12*nyears)/12,24)
days<-days[order(days)]
dates<-days+soilmoist$TIME/60/24/(12) # dates for hourly output
soilmoist<-cbind(dates,soilmoist)
for(i in 1:10){
  if(i==1){
    plot(soilmoist[,i+3]~soilmoist[,1], ylim=c(0,0.5),xlab = "Year", ylab = "Soil Moisture (%
      vol)",col=i,type = "l",main=paste("soil moisture, ",minshade,"% shade"))
  }else{
    points(soilmoist[,i+3]~soilmoist[,1],col=i,type = "l")
  }
}
```

It is also possible to run the model for more frequent intervals than once per month by changing the parameter ```timeinterval```. If ```timeinterval``` is set to something less than 365, then the model runs as described in the previous section with three iterations per day, and uses the 'spline' function to interpolate the climate data from monthly to the chosen time interval. If 'timeinterval' is set to the maximum value of 365, it will apply the three iterations to the first day of the simulation to get started and then, from that day on, use the previous day's conditions as the initial state for the next day's conditions.

Here is an example running the model for 365 days and plotting soil moisture:
```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
timeinterval <- 365
micro <- micro_global(loc = longlat, runmoist = 1, soiltype = 11, timeinterval = timeinterval)
soilmoist <- as.data.frame(micro$soilmoist)  # get the minimum shade soil moisture output
minshade <- micro$minshade[1]  # get the value for minimum shade
# append dates
days <- rep(seq(1, timeinterval), 24)
days <- days[order(days)]
dates <- days + soilmoist$TIME/60/24 - 1  # dates for hourly output
soilmoist <- cbind(dates, soilmoist)
for (i in 1:10) {
    if (i == 1) {
        plot(soilmoist[, i + 3] ~ soilmoist[, 1], ylim = c(0, 0.5), xlab = "Day of Year", ylab = "Soil Moisture (% vol)", col = i, type = "l", main = paste("soil moisture, ", minshade, "% shade"))
    } else {
        points(soilmoist[, i + 3] ~ soilmoist[, 1], col = i, type = "l")
    }
}
```

Here are the associated soil temperatures:
```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
soil <- as.data.frame(micro$soil)  # get the minimum shade soil temperature output
minshade <- micro$minshade[1]  # get the value for minimum shade
# append dates
days <- rep(seq(1, timeinterval), 24)
days <- days[order(days)]
dates <- days + soil$TIME/60/24 - 1  # dates for hourly output
soil <- cbind(dates, soil)
for (i in 1:10) {
    if (i == 1) {
        plot(soil[, i + 3] ~ soil[, 1], ylim = c(-20, 70), xlab = "Day of Year", ylab = "Soil Temperature (deg C)", col = i, type = "l", main = paste("soil temperature, ", minshade, "% shade"))
    } else {
        points(soil[, i + 3] ~ soil[, 1], col = i, type = "l")
    }
}
```


### Snow

Snow has a powerful influence on soil temperatures and can be simulated in the NicheMapR microclimate model by setting the parameter ```snowmodel = 1```. This invokes an algorithm that builds up to 8 nodes of snow above the surface as a function of incoming precipitation and losses due to melting. The model incorporates a declining albedo with time since snow and has five user-controlled parameters: ```snowtemp``` - the air temperature (&deg;C) at which rain falls as snow, ```snowdens``` (Mg/m3) - the density of the snow, ```snowmelt``` (-) - the proportion of calculated snow melt that doesn't refreeze into the snowpack,  ```rainmelt``` (-) a parameter in an equation from Anderson's SNOW17 model (<http://www.nws.noaa.gov/oh/hrl/nwsrfs/users_manual/part2/_pdf/22snow17.pdf>) that melts snow with rainfall as a function of air temperature. The snow density can also be made to vary as a linear function of day of year by providing a slope and intercept as the variable ```densfun```, which is only used if non-zero. There is also a parameter ```rainfrac``` which controls the fraction of rain that falls on the first day of the month, with 0 meaning rain falls evenly. This parameter allows something other than an even intensity of rainfall when interpolating the montly rainfall data. This snow model has been tested with favourable results in the Australian Alps and through the USA (Kearney et al. in prep).

Here is the consequence of invoking the snow model for soil temperature predictions at Madison, Wisconsin as well as the snow depth and snowfall predictions (running for 2 years to get a steady periodic and plotting only the 2nd year):

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
timeinterval <- 365
nyears <- 2  # running two years, first one acting as a 'burn in' year and discarded
micro <- micro_global(loc = longlat, runmoist = 1, snowmodel = 1, timeinterval = timeinterval, nyears = 2)
soil <- as.data.frame(micro$soil)[(365 * 24 + 1):(365 * 24 * nyears), ]  # get the minimum shade soil temperature output, discarding the first year
metout <- as.data.frame(micro$metout)[(365 * 24 + 1):(365 * 24 * nyears), ]  # get the minimum shade above ground conditions, discarding the first year
minshade <- micro$minshade[1]  # get the value for minimum shade
# append dates
days <- rep(seq(366, timeinterval * nyears), 24)
days <- days[order(days)]
dates <- days + soil$TIME / 60 / 24 - 1  # dates for hourly output
soil <- cbind(dates, soil)
metout <- cbind(dates, metout)
for (i in 1:10) {
    if (i == 1) {
        plot(soil[, i + 3] ~ soil[, 1], ylim = c(-20, 70), xlab = "Day of Year", ylab = "Soil Temperature (deg C)", col = i, type = "l", main = paste("soil  temperature, ", minshade, "% shade"))
    } else {
        points(soil[, i + 3] ~ soil[, 1], col = i, type = "l")
    }
}
plot(metout$SNOWDEP ~ metout$dates, xlab = "Time of Day (min)", ylab = "snow depth, cm / snow fall, mm", type = "h", main = paste("snow depth (cm) and snow fall (mm) ", minshade, "% shade", sep = ""), col = "light blue")
points(metout$SNOWFALL ~ metout$dates, xlab = "Time of Day (min)", type = "h", col = "blue")
```


### Other outputs of interest

Plant transpiration and water potential 
```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
plant <- as.data.frame(micro$plant)[(365 * 24 + 1):(365 * 24 * nyears), ]  # get the minimum shade plant output, discarding the first year
soilpot <- as.data.frame(micro$soilpot)[(365 * 24 + 1):(365 * 24 * nyears), ]  # get the minimum shade soil water potential output, discarding the first year
plot(plant$TRANS ~ dates, xlab = "Time of Day (min)", ylab = "transpiration, g/h/m^2", type = "h", main = paste("plant transpiration rate ", minshade, "% shade", sep = ""), col = "light blue")
plot(plant$LPT ~ dates, xlab = "Time of Day (min)", ylab = "leaf water potential, J/kg", type = "h", main = paste("plant water potential ", minshade, "% shade", sep = ""), col = "darkgreen")
for (i in 2:10) {
        points(plant[, i + 4] ~ dates, col = i, type = "l")
}

for (i in 2:10) {
    if (i == 2) {
        plot(soilpot[, i + 2] ~ dates, ylim = c(-250, 0), xlab = "Day of Year", ylab = "soil water potential (J/kg)", col = i, type = "l", main = paste("soil water potential ", minshade, "% shade"))
    } else {
        points(soilpot[, i + 2] ~ dates, col = i, type = "l")
    }
}
```

Solar spectra

```{r, message=FALSE, warnings=FALSE, fig.width=7, fig.height=6}
micro <- micro_global(loc = longlat, solonly = 1, clearsky = 1, lamb = 1, IUV = 1, microclima = 1)

metout <- as.data.frame(micro$metout)
drlam <- as.data.frame(micro$drlam) # direct solar
srlam <- as.data.frame(micro$srlam) # scattered solar
drrlam <- as.data.frame(micro$drrlam) # direct Rayleigh solar
diffuse_frac <- micro$diffuse_frac

par(mfrow = c(3, 1))
# solstice (mid December)
lambdas <- as.numeric(colnames(drlam[3:113])) # get wavelengths
plot(lambdas, drrlam[277, 3:113], type = 'l', ylab = 'irradiance, W/m2', xlab = 'wavelength, nm', ylim = c(0, 1.8), main = 'solar spectra, mid December')
points(lambdas, drlam[277, 3:113], type = 'l', col = 2)
points(lambdas, srlam[277, 3:113], type = 'l', col = 3)

# equinox (mid March)
plot(lambdas, drrlam[61, 3:113], type = 'l', ylab = 'irradiance, W/m2', xlab = 'wavelength, nm', ylim = c(0, 1.8), main = 'solar spectra, mid March')
points(lambdas, drlam[61, 3:113], type = 'l', col = 2)
points(lambdas, srlam[61, 3:113], type = 'l', col = 3)

# solstice (mid June)
plot(lambdas, drrlam[133, 3:113], type = 'l', ylab = 'irradiance, W/m2', xlab = 'wavelength, nm', ylim = c(0, 1.8), main = 'solar spectra, mid June')
points(lambdas, drlam[133, 3:113], type = 'l', col = 2)
points(lambdas, srlam[133, 3:113], type = 'l', col = 3)

library(pracma)
row2calc <- 61
SOLR_dir <- trapz(lambdas, t(drlam[row2calc, 3:113]))
SOLR_dif <- trapz(lambdas, t(srlam[row2calc, 3:113]))

metout$SOLR[row2calc] # horizontal plane total solar, W/m^2
SOLR_dir + SOLR_dif # sum the direct and diffuse integrations
diffuse_frac[row2calc] # microclima predicted diffuse fraction
SOLR_dif / (SOLR_dir + SOLR_dif) # NMR predicted diffuse fraction

```

### Other parameters

In the help for the ```micro_global``` function you'll find information on a number of other parameters that may be of use for your application. These allow control of the eccentricity of the earth's orbit ```EC``` (if you want to go way back into the past), changes in the timing of daily maxima and minima for air temperature, cloud cover, wind speed and relative humidity ```TIMINS``` and ```TIMAXS```, the optical effects of the amount of water in the upper atmosphere ```CMH2O```, whether to apply the Global Aerosol Data Set (package ```GADS```) to determine solar attenuation due to dust etc. in the atmosphere ```rungads```, whether the model outputs csv files ```writecsv``` or csv files of all inputs ```writeinput``` (mainly used for the maintainers of the program when testing the Fortran code).

## References

Campbell, G. S. 1985. Soil Physics with Basic: Transport Models for Soil-Plant Systems. Elsevier, Amsterdam.

Campbell, G. S., J. D. J. Jungbauer, W. R. Bidlake, and R. D. Hungerford. 1994. Predicting the effect of temperature on soil thermal conductivity. Soil Science 158:307-313.

Campbell, G. S., and J. M. Norman. 1998. Environmental Biophysics. Springer, New York.

Kearney, M. R., A. Shamakhy, R. Tingley, D. J. Karoly, A. A. Hoffmann, P. R. Briggs, and W. P. Porter. 2014. Microclimate modelling at macro scales: a test of a general microclimate model integrated with gridded continental-scale soil and weather data. Methods in Ecology and Evolution 5:273-286.

Schmidt-Nielsen, B., and K. Schmidt-Nielsen. 1950. Evaporative Water Loss in Desert Rodents in Their Natural Habitat. Ecology 31:75-85.