knitr::opts_chunk$set( eval = TRUE )
This vignette runs through the core inputs and outputs of the NicheMapR microclimate model. It provides a basis for developing scripts and functions that are customised to particular databases of climate/weather, terrain, vegetation and soil properties.
writecsv <- 0 # make Fortran program write output as csv files microdaily <- 0 # run microclimate where each day is iterated 3 times starting with the initial condition of uniform soil temp at mean monthly temperature (0) or with initial conditions of each day based on final conditions of previous day (1) runshade <- 1 # run the model twice, once for each shade level (1) or just for the first shade level (0)? IR <- 0 # Clear-sky longwave radiation computed using Campbell and Norman (1998) eq. 10.10 (includes humidity) (0) or Swinbank formula (1) runmoist <- 1 # run soil moisture model (0=no, 1=yes)? snowmodel <- 1 # run the snow model (0=no, 1=yes)? - note that this runs slower hourly <- 0 # use hourly input data? solonly <- 0 # Only run SOLRAD to get solar radiation? 1=yes, 0=no lamb <- 0 # Return wavelength-specific solar radiation output? IUV <- 0 # Use gamma function for scattered solar radiation? message <- 0 #allow the Fortran integrator to output warnings? (1) or not (0) fail <- 1*24*365 # how many restarts of the integrator before the Fortran program quits (avoids endless loops when solutions can't be found)
ndays <- 12 # number of time intervals to generate predictions for over a year (must be 12 <= x <=365) doy <- c(15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349) # middle day of each month idayst <- 1 # start day (legacy parameter) ida <- ndays # end day (legacy parameter) HEMIS <- 1 # chose hemisphere ALAT <- 43 # degrees latitude AMINUT <- 4.383 # minutes latitude ALONG <- 89 # degrees longitude ALMINT <- 24.074 # minutes latitude ALREF <- 89 # reference longitude for time zone EC <- 0.0167238 # Eccenricity of the earth's orbit (current value 0.0167238, ranges between 0.0034 to 0.058)
RUF <- 0.004 # Roughness height (m), , e.g. sand is 0.05, grass may be 2.0, current allowed range: 0.001 (snow) - 2.0 cm. Refhyt <- 2 # Reference height (m), reference height at which air temperature, wind speed and relative humidity input data are measured Usrhyt <- 0.01# local height (m) at which air temperature, relative humidity and wind speed calculatinos will be made ZH <- 0 # heat transfer roughness height (m) for Campbell and Norman air temperature/wind speed profile (invoked if greater than 1, 0.02 * canopy height in m if unknown) D0 <- 0 # zero plane displacement correction factor (m) for Campbell and Norman air temperature/wind speed profile (0.6 * canopy height in m if unknown) # Next four parameters are segmented velocity profiles due to bushes, rocks etc. on the surface #IF NO EXPERIMENTAL WIND PROFILE DATA SET ALL THESE TO ZERO! (then roughness height is based on the parameter RUF) Z01 <- 0 # Top (1st) segment roughness height(m) Z02 <- 0 # 2nd segment roughness height(m) ZH1 <- 0 # Top of (1st) segment, height above surface(m) ZH2 <- 0 # 2nd segment, height above surface(m)
SLE <- 0.96 # substrate longwave IR emissivity (decimal %), typically close to 1 REFL <- 0.10 # substrate solar reflectivity (decimal %) CMH2O <- 1 # precipitable cm H2O in air column, 0.1 = VERY DRY; 1.0 = MOIST AIR CONDITIONS; 2.0 = HUMID, TROPICAL CONDITIONS (note this is for the whole atmospheric profile, not just near the ground) # Aerosol extinction coefficient profile # the original profile from Elterman, L. 1970. Vertical-attenuation model with eight surface meteorological ranges 2 to 13 kilometers. U. S. Airforce Cambridge Research Laboratory, Bedford, Mass. #TAI <- c(0.0670358341290886,0.0662612704779235,0.065497075238002,0.0647431301168489,0.0639993178022531,0.0632655219571553,0.0625416272145492,0.0611230843885423,0.0597427855962549,0.0583998423063099,0.0570933810229656,0.0558225431259535,0.0545864847111214,0.0533843764318805,0.0522154033414562,0.0499736739981675,0.047855059159556,0.0458535417401334,0.0439633201842001,0.0421788036108921,0.0404946070106968,0.0389055464934382,0.0374066345877315,0.0359930755919066,0.0346602609764008,0.0334037648376212,0.0322193394032758,0.0311029105891739,0.0300505736074963,0.0290585886265337,0.0281233764818952,0.0272415144391857,0.0264097320081524,0.0256249068083005,0.0248840604859789,0.0241843546829336,0.0235230870563317,0.0228976873502544,0.0223057135186581,0.0217448478998064,0.0212128934421699,0.0207077699817964,0.0202275105711489,0.0197702578594144,0.0193342605242809,0.0189178697551836,0.0177713140039894,0.0174187914242432,0.0170790495503944,0.0167509836728154,0.0164335684174899,0.0161258546410128,0.0158269663770596,0.0155360978343254,0.0152525104459325,0.0149755299703076,0.0147045436435285,0.0144389973831391,0.0141783930434343,0.0134220329447663,0.0131772403830191,0.0129356456025128,0.0126970313213065,0.0124612184223418,0.0122280636204822,0.01199745718102,0.0115436048739351,0.0110993711778668,0.0108808815754663,0.0106648652077878,0.0104513876347606,0.0102405315676965,0.00982708969547694,0.00962473896278535,0.00903679230300494,0.00884767454432418,0.0083031278398166,0.00796072474935954,0.00755817587626185,0.00718610751850881,0.00704629977586921,0.00684663903049612,0.00654155580333479,0.00642947339729728,0.00627223096874308,0.00603955966866779,0.00580920937536261,0.00568506186880564,0.00563167068287251,0.00556222005081865,0.00550522989971023,0.00547395763028062,0.0054478983436216,0.00541823364504573,0.00539532163908382,0.00539239864119488,0.00541690124712384,0.00551525885358836,0.00564825853509463,0.00577220185074264,0.00584222986640171,0.00581645238345584,0.00566088137411449,0.00535516862329704,0.00489914757707667,0.00432017939770409,0.0036813032251836,0.00309019064543606,0.00270890436501562,0.00276446109239711,0.00356019862584603) # the values extracted from GADS for Madison TAI <- c(0.269904738,0.266147825,0.262442906,0.258789404,0.255186744,0.251634356,0.248131676,0.2412732,0.234606887,0.228128378,0.221833385,0.215717692,0.20977715,0.204007681,0.198405272,0.187685927,0.177588357,0.168082846,0.159140695,0.150734206,0.142836655,0.135422274,0.128466227,0.12194459,0.115834329,0.110113284,0.104760141,0.099754417,0.09507644,0.090707328,0.086628967,0.082823998,0.07927579,0.075968428,0.072886691,0.070016034,0.067342571,0.064853053,0.062534858,0.060375964,0.058364941,0.056490925,0.054743609,0.053113222,0.051590514,0.050166738,0.046408775,0.045302803,0.044259051,0.043271471,0.042334415,0.041442618,0.040591184,0.039775572,0.038991583,0.038235345,0.037503301,0.036792197,0.036099067,0.034101935,0.033456388,0.032817888,0.032184949,0.031556287,0.030930816,0.030307633,0.029065372,0.027825562,0.027205981,0.026586556,0.025967391,0.025348692,0.024114005,0.023498886,0.021669152,0.021066668,0.019292088,0.018144698,0.016762709,0.015451481,0.014949794,0.014224263,0.013093462,0.012670686,0.012070223,0.011164062,0.010241734,0.009731103,0.009507687,0.009212683,0.008965785,0.008827751,0.008710756,0.008574128,0.008462605,0.008446967,0.008539475,0.009015237,0.009748444,0.010586023,0.011359647,0.011901268,0.012062153,0.011735443,0.010882215,0.009561062,0.007961182,0.006438984,0.005558204,0.006133532,0.009277754)
ALTT <- 226 # elevation (m) slope <- 0 # slope (degrees, range 0-90) azmuth <- 180 # aspect (degrees, 0 = North, range 0-360) hori <- rep(0, 24) # enter the horizon angles (degrees) so that they go from 0 degrees azimuth (north) clockwise in 15 degree intervals VIEWF <- 1-sum(sin(hori*pi/180))/length(hori) # convert horizon angles to radians and calc view factor(s) minshade <- 0 # minimum available shade (%) maxshade <- 90 # maximum available shade (%) PCTWET <- 0 # percentage of surface area acting as a free water surface (%)
DEP <- c(0, 2.5, 5, 10, 15, 20, 30, 50, 100, 200) # Soil nodes (cm) - keep spacing close near the surface, last value is where it is assumed that the soil temperature is at the annual mean air temperature ERR <- 2.0 # Integrator error for soil temperature calculations
TIMINS <- c(0, 0, 1, 1) # time of minima for air temp, wind, humidity and cloud cover (h), air & wind mins relative to sunrise, humidity and cloud cover mins relative to solar noon TIMAXS <- c(1, 1, 0, 0) # time of maxima for air temp, wind, humidity and cloud cover (h), air temp & wind maxs relative to solar noon, humidity and cloud cover maxs relative to sunrise TMINN <- c(-14.3,-12.1,-5.1,1.2,6.9,12.3,15.2,13.6,8.9,3,-3.2,-10.6) # minimum air temperatures (deg C) TMAXX <- c(-3.2,0.1,6.8,14.6,21.3,26.4,29,27.7,23.3,16.6,7.8,-0.4) # maximum air temperatures (deg C) RHMINN <- c(50.2,48.4,48.7,40.8,40,42.1,45.5,47.3,47.6,45,51.3,52.8) # min relative humidity (%) RHMAXX <- c(100,100,100,100,100,100,100,100,100,100,100,100) # max relative humidity (%) WNMINN <- c(4.9,4.8,5.2,5.3,4.6,4.3,3.8,3.7,4,4.6,4.9,4.8) # min wind speed (m/s) WNMAXX <- c(4.9,4.8,5.2,5.3,4.6,4.3,3.8,3.7,4,4.6,4.9,4.8) # max wind speed (m/s) CCMINN <- c(50.3,47,48.2,47.5,40.9,35.7,34.1,36.6,42.6,48.4,61.1,60.1) # min cloud cover (%) CCMAXX <- c(50.3,47,48.2,47.5,40.9,35.7,34.1,36.6,42.6,48.4,61.1,60.1) # max cloud cover (%) RAINFALL <- c(28,28.2,54.6,79.7,81.3,100.1,101.3,102.5,89.7,62.4,54.9,41.2) # monthly mean rainfall (mm) TAIRhr <- rep(0, 24 * ndays) # hourly air temperature input (degrees C) (not used if hourly = 0) RHhr <- rep(0, 24 * ndays) # hourly relative humidity input (%) (not used if hourly = 0) WNhr <- rep(0, 24 * ndays) # hourly wind speed input (m/s) (not used if hourly = 0) CLDhr <- rep(0, 24 * ndays) # hourly cloud cover input (%) (not used if hourly = 0) SOLRhr <- rep(0, 24 * ndays) # hourly solar radiation input (W/m2) (not used if hourly = 0) RAINhr <- rep(0, 24 * ndays) # hourly rainfall input (mm) (not used if hourly = 0) ZENhr <- rep(-1, 24 * ndays) # hourly zenith angle input (degrees) (not used if hourly = 0 or if negative) IRDhr <- rep(-1, 24 * ndays) # hourly longwave radiation input (W/m2) (not used if hourly = 0 or if negative) tannul <- mean(c(TMAXX,TMINN)) # annual mean temperature for getting monthly deep soil temperature (deg C) tannulrun <- rep(tannul,ndays) # monthly deep soil temperature (2m) (deg C) SoilMoist <- c(0.42,0.42,0.42,0.43,0.44,0.44,0.43,0.42,0.41,0.42,0.42,0.43) # soil moisture (decimal %, 1 means saturated) # creating the arrays of environmental variables that are assumed not to change with month for this simulation MAXSHADES <- rep(maxshade,ndays) # daily max shade (%) MINSHADES <- rep(minshade,ndays) # daily min shade (%) SLES <- rep(SLE,ndays) # set up vector of ground emissivities for each day REFLS <- rep(REFL,ndays) # set up vector of soil reflectances for each day PCTWET <- rep(PCTWET,ndays) # set up vector of soil wetness for each day
# set up a profile of soil properites with depth for each day to be run Numtyps <- 2 # number of soil types Nodes <- matrix(data = 0, nrow = 10, ncol = ndays) # array of all possible soil nodes Nodes[1,1:ndays] <- 3 # deepest node for first substrate type Nodes[2,1:ndays] <- 9 # deepest node for second substrate type # soil thermal parameters Thcond <- 1.25 # soil minerals thermal conductivity (W/mC) Density <- 2560 # soil minerals density (kg/m3) SpecHeat <- 870 # soil minerals specific heat (J/kg-K) BulkDensity <- 2560 # soil bulk density (kg/m3) SatWater <- 0.26 # volumetric water content at saturation (0.1 bar matric potential) (m3/m3) Clay <- 20 # clay content for matric potential calculations (%) #SoilMoist <- rep(SoilMoist,timeinterval) # soil moisture Density <- Density/1000 # density of minerals - convert to Mg/m3 BulkDensity <- BulkDensity/1000 # density of minerals - convert to Mg/m3 # now make the depth-specific soil properties matrix # columns are: #1) bulk density (Mg/m3) #2) volumetric water content at saturation (0.1 bar matric potential) (m3/m3) #3) thermal conductivity (W/mK) #4) specific heat capacity (J/kg-K) #5) mineral density (Mg/m3) soilprops <- matrix(data = 0, nrow = 10, ncol = 6) # create an empty soil properties matrix soilprops[1,1] <- BulkDensity # insert soil bulk density to profile 1 soilprops[2,1] <- BulkDensity # insert soil bulk density to profile 2 soilprops[1,2] <- SatWater # insert saturated water content to profile 1 soilprops[2,2] <- SatWater # insert saturated water content to profile 2 soilprops[1,3] <- Thcond # insert thermal conductivity to profile 1 soilprops[2,3] <- Thcond # insert thermal conductivity to profile 2 soilprops[1,4] <- SpecHeat # insert specific heat to profile 1 soilprops[2,4] <- SpecHeat # insert specific heat to profile 2 soilprops[1,5] <- Density # insert mineral density to profile 1 soilprops[2,5] <- Density # insert mineral density to profile 2 soilinit <- rep(tannul,20) # make inital soil temps equal to mean annual
# note that these are set for sand (Table 9.1 in Campbell and Norman, 1995) PE <- rep(0.7,19) #air entry potential J/kg KS <- rep(0.0058,19) #saturated conductivity, kg s/m3 BB <- rep(1.7,19) #soil 'b' parameter BD <- rep(1.3,19) # soil bulk density, Mg/m3 DD <- rep(2.56,19) # soil density (Mg/m3) (19 values descending through soil for specified soil nodes in parameter DEP and points half way between) L <- c(0,0,8.18990859,7.991299442,7.796891252,7.420411664,7.059944542,6.385001059,5.768074989,4.816673431,4.0121088,1.833554792,0.946862989,0.635260544,0.804575,0.43525621,0.366052856,0,0)*10000 # root density at each node, mm/m3 (from Campell 1985 Soil Physics with Basic, p. 131) R1 <- 0.001 # root radius, m RW <- 2.5e+10 # resistance per unit length of root, m3 kg-1 s-1 RL <- 2e+6 # resistance per unit length of leaf, m3 kg-1 s-1 PC <- -1500 # critical leaf water potential for stomatal closure, J kg-1 SP <- 10 # stability parameter for stomatal closure equation, - IM <- 1e-06 # maximum allowable mass balance error, kg MAXCOUNT <- 500 # maximum iterations for mass balance, - LAI <- 0.1 # leaf area index, used to partition traspiration/evaporation from PET rainmult <- 1 # rainfall multiplier to impose catchment maxpool <- 10 # max depth for water pooling on the surface, mm (to account for runoff) rainhourly <- 0 # Is hourly rain input being supplied (1 = yes, 0 = no)? evenrain <- 1 # spread daily rainfall evenly across 24hrs (1) or one event at midnight (2) SoilMoist_Init <- rep(0.2,10) # initial soil water content for each node, m3/m3 moists <- matrix(nrow=10, ncol = ndays, data=0) # set up an empty vector for soil moisture values through time moists[1:10,] <- SoilMoist_Init # insert inital soil moisture
snowtemp <- 1.5 # temperature at which precipitation falls as snow (used for snow model) snowdens <- 0.375 # snow density (mg/m3) densfun <- c(0.5979, 0.2178, 0.001, 0.0038) # slope and intercept of model of snow density as a linear function of snowpack age if first two values are nonzero, and following the exponential function of Sturm et al. 2010 J. of Hydromet. 11:1380-1394 if all values are non-zero; if it is c(0,0,0,0) then fixed density used snowmelt <- 0.9 # proportion of calculated snowmelt that doesn't refreeze undercatch <- 1. # undercatch multipier for converting rainfall to snow rainmelt <- 0.0125 # parameter in equation from Anderson's SNOW-17 model that melts snow with rainfall as a function of air temp snowcond <- 0 # effective snow thermal conductivity W/mC (if zero, uses inbuilt function of density) intercept <- maxshade / 100 * 0.3 # snow interception fraction for when there's shade (0-1) grasshade <- 0 # if 1, means shade is removed when snow is present, because shade is cast by grass/low shrubs
# intertidal simulation input vector (col 1 = tide in(1)/out(0), col 2 = sea water temperature in deg C, col 3 = % wet from wave splash) tides <- matrix(data = 0, nrow = 24*ndays, ncol = 3) # matrix for tides
# input parameter vector microinput <- c(ndays, RUF, ERR, Usrhyt, Refhyt, Numtyps, Z01, Z02, ZH1, ZH2, idayst, ida, HEMIS, ALAT, AMINUT, ALONG, ALMINT, ALREF, slope, azmuth, ALTT, CMH2O, microdaily, tannul, EC, VIEWF, snowtemp, snowdens, snowmelt, undercatch, rainmult, runshade, runmoist, maxpool, evenrain, snowmodel, rainmelt, writecsv, densfun, hourly, rainhourly, lamb, IUV, RW, PC, RL, SP, R1, IM, MAXCOUNT, IR, message, fail, snowcond, intercept, grasshade, solonly, ZH, D0) # Final input list - all these variables are expected by the input argument of the Fortran microclimate subroutine micro <- list(microinput = microinput, tides = tides, doy = doy, SLES = SLES, DEP = DEP, Nodes = Nodes, MAXSHADES = MAXSHADES, MINSHADES = MINSHADES, TIMAXS = TIMAXS, TIMINS = TIMINS, TMAXX = TMAXX, TMINN = TMINN, RHMAXX = RHMAXX, RHMINN = RHMINN, CCMAXX = CCMAXX, CCMINN = CCMINN, WNMAXX = WNMAXX, WNMINN = WNMINN, TAIRhr = TAIRhr, RHhr = RHhr, WNhr = WNhr, CLDhr = CLDhr, SOLRhr = SOLRhr, RAINhr = RAINhr, ZENhr = ZENhr, IRDhr = IRDhr, REFLS = REFLS, PCTWET = PCTWET, soilinit = soilinit, hori = hori, TAI = TAI, soilprops = soilprops, moists = moists, RAINFALL = RAINFALL, tannulrun = tannulrun, PE = PE, KS = KS, BB = BB, BD = BD, DD = DD, L = L, LAI = LAI)
microut <- microclimate(micro) # run the model in Fortran
metout <- as.data.frame(microut$metout[1:(ndays * 24), ]) # retrieve above ground microclimatic conditions, min shade shadmet <- as.data.frame(microut$shadmet[1:(ndays * 24), ]) # retrieve above ground microclimatic conditions, max shade soil <- as.data.frame(microut$soil[1:(ndays * 24), ]) # retrieve soil temperatures, minimum shade shadsoil <- as.data.frame(microut$shadsoil[1:(ndays * 24), ]) # retrieve soil temperatures, maximum shade soilmoist <- as.data.frame(microut$soilmoist[1:(ndays * 24), ]) # retrieve soil moisture, minimum shade shadmoist <- as.data.frame(microut$shadmoist[1:(ndays * 24), ]) # retrieve soil moisture, maximum shade humid <- as.data.frame(microut$humid[1:(ndays * 24), ]) # retrieve soil humidity, minimum shade shadhumid <- as.data.frame(microut$shadhumid[1:(ndays * 24), ]) # retrieve soil humidity, maximum shade soilpot <- as.data.frame(microut$soilpot[1:(ndays * 24), ]) # retrieve soil water potential, minimum shade shadpot <- as.data.frame(microut$shadpot[1:(ndays * 24), ]) # retrieve soil water potential, maximum shade
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.