LandscapeDNDC 1.37.0
wc_canoak.h
1#ifndef M1D_MOD_WATERCYCLE_CANOAK_H
2#define M1D_MOD_WATERCYCLE_CANOAK_H
3
4#define WCC_FLEX 2
5#define sze LConf::flPot + WCC_FLEX // canopy layers is 40 with flex
6#define sze3 3 * LConf::flPot + WCC_FLEX // 3 times canopy layers is 120 with flex
7#define szeang 18 + 1 // number of sky angle classes is 18 with flex ( [sk] flex suddenly 1?)
8#define soilsze 12 // number of soil layers is 10 (100?) with flex
9
10#include "watercyclemodule.h"
11
12using std::string;
13
14
15// structure for input variables
16struct input_variables{
17 int dayy; // day
18 // rg int hhrr; // hour
19 double hhrr; // hour
20 double ta; // air temperature, C
21 double rglobal; // global radiation, W m-2
22 double parin; // photosynthetically active radiation, micromole m-2 s-1
23 double pardif; // diffuse PAR, micromol m-2 s-1
24 double ea; // vapor pressure, kPa
25 double wnd; // wind speed, m s-1
26 double ppt; // precipitation, mm per hour
27 double co2air; // CO2 concentration, ppm
28 double press_mb; // air pressure, mb
29 double tsoil; // soil temperature in 50 cm, degree C
30 double soilmoisture;// soil moisture in 15 cm, %
31 //rg long int flag; // input coding
32};
33
34// structure for time variables
35struct time_variables{
36 double local_time;
37 long int daytime; // day+hour coding, e.g. 1230830
38 int year; // year
39 int days; // day
40 //rg int jdold; // previous day (only needed for calculating averages in the original code)
41 //rg int fileyear; // year for filenames
42 //rg int filestart; // time for filestart
43 //rg int count; // number of iterations
44 //rg int leafout; // day of leaf out
45 //rg int fulleaf; // date of full leaf
46 double lai; // lai as a function of time
47};
48
49// water, carbon and energy exchange
50struct flux_variables {
51 double photosyn; // photosynthesis, um m-2 s-1
52 double nee; // net ecosystem exchange, umol m-2 s-1
53 double potevap; // potential evaporation, mm m-2 s-1 // rg
54 double transp; // canopy transpiration, mm m-2 s-1
55 double soilevap; // soil evaporation, mm m-2 s-1
56 double canopyevap; // canopy evaporation, mm m-2 s-1
57 double through; // throughfall = water flux into soil, mm m-2 s-1
58 double surfrun; // surface runoff = water flux lateral to surface, mm m-2 s-1
59
60 flux_variables ()
61 {
62 photosyn = 0.0;
63 nee = 0.0;
64 potevap = 0.0;
65 transp = 0.0;
66 soilevap = 0.0;
67 canopyevap = 0.0;
68 through = 0.0;
69 surfrun = 0.0;
70 }
71};
72
73// structure for meteorological variables
74struct meteorology{
75 double ustar; // friction velocity, m s-1
76 double ustarnew; // updated friction velocity with new H, m s-1
77 double rhova_g; // absolute humidity, g m-3
78 double rhova_kg; // absolute humidity, kg m-3
79 double sensible_heat_flux; // sensible heat flux, W M-2
80 //rg double H_old; // old sensible heat flux, W m-2 (replaced with new static variable)
81 double air_density; // air density, kg m-3
82 double T_Kelvin; // absolute air temperature, K
83 double dispersion[sze3][sze]; // Lagrangian dispersion matrix, s m-1
84 double zl; // z over L, height normalized by the Obukhov length
85 double press_kpa; // station pressure, kPa
86 double press_bars; // station pressure, bars
87 double press_Pa; // pressure, Pa
88 double pstat273; // gas constant computations
89 double air_density_mole; // air density, mole m-3
90 double relative_humidity; // relative humidity
91
92 meteorology ()
93 {
94 ustar = 0.0;
95 ustarnew = 0.0;
96 rhova_g = 0.0;
97 rhova_kg = 0.0;
98 sensible_heat_flux = 0.0;
99 air_density = 0.0;
100 T_Kelvin = 0.0;
101 zl = 0.0;
102 press_kpa = 0.0;
103 press_bars = 0.0;
104 press_Pa = 0.0;
105 pstat273 = 0.0;
106 air_density_mole = 0.0;
107 relative_humidity = 0.0;
108 }
109};
110
111// structure for surface resistances and conductances
112struct surface_resistances{
113 double gcut; // cuticle conductances, mol m-2 s-1
114 double kballstr; // drought factor to Ball Berry Coefficient
115 double rcuticle; // cuticle resistance, m2 s mol-1
116};
117
118// structure for plant and physical factors
119struct factors{
120 double latent; // latent heat of vaporization, J kg-1
121 double latent18; // latent heat of vaporization times molecular mass of vapor, 18 g mol-1
122 double heatcoef; // factor for sensible heat flux density
123 double a_filt; // filter coefficients
124 double b_filt; // filter coefficients
125 double co2; // CO2 factor, ma/mc * rhoa (mole m-3)
126};
127
128// structure for bole respiration and structure
129struct bole_respiration_structure{
130 double factor; // base respiration, micromoles m-2 s-1, data from Edwards
131 double respiration_mole; // bole respiration, micromol m-2 s-1
132 double respiration_mg; // bole respiration, mg CO2 m-2 s-1
133 double calc; // calculation factor
134 double layer[sze]; // bole pai per layer
135};
136
137// structure for canopy architecture
138struct canopy_architecture{
139 //rg?? double bdens[soilsze]; // probability density of leaf angle
140 double bdens[10]; // probability density of leaf angle
141};
142
143// structure for non dimensional variables
144struct non_dimensional_variables{
145 // Prandtl Number
146 double pr;
147 double pr33;
148 // Schmidt number for vapor
149 double sc;
150 double sc33;
151 // Schmidt number for CO2
152 double scc;
153 double scc33;
154 // Schmidt number for ozone
155 double sco3;
156 double sco333;
157 // Grasshof number
158 double grasshof;
159 // multiplication factors with leaf length and diffusivity
160 double lfddv;
161 double lfddh;
162};
163
164// boundary layer resistances
165struct boundary_layer_resistances{
166 double vapor; // resistance for water vapor, s/m
167 double heat; // resistance for heat, s/m
168 double co2; // resistance for CO2, s/m
169
170 boundary_layer_resistances ()
171 {
172 vapor = 0.0;
173 heat = 0.0;
174 co2 = 0.0;
175 }
176};
177
178struct other_globals{
179 // Parameter are here set on default values (plus for gloabl definition).
180 // Later new parameters are read in from canisotope_param.txt in read_paramter()
181 //rg double parameter[250];// declare parameter array that will be read in from canisotope_param.txt in read_paramter()
182 //rg char parameter_string[5][256];
183 double paraopt[31]; //rg
184
185 // canopy structure variables
186 //rg double ht = 33.0;//changed: Canopy height
187 //rg double lai = 6.4;//changed: original lai=6.0 Leaf area index
188 int jtot;//=40; // number of canopy layers
189 int jktot;//=41; // jtot + 1
190 int jtot3;//=120; // number of layers in the domain, three times canopy height
191 int izref;//=52; // changed: array value of reference ht at 43 m, jtot/ht
192
193 double hm; // canopy height
194 double laiT; // potential leaf area index during the year
195 double pai;// = 1.4; // changed: plant area index
196 double vcopt;// = 66.3; // changed: original: 70, carboxylation rate at optimal temperature, umol m-2 s-1
197 double jmopt;// = 127.5; // changed: original: 170, electron transport rate at optimal temperature, umol m-2 s-1
198 double rd25;// = 0.34; // changed: new: 1.91 original: 0.34, dark respiration at 25 C, rd25= 0.34 umol m-2 s-1
199 double ustar_ref;//=1.00; // reference ustar for Dij
200
201 double delz;// = 0.825; // changed: height of each layer, ht/jtot
202 double zh65;//=0.019697; // changed: 0.65/ht
203
204 // Consts for Photosynthesis model and kinetic equations for Vcmax and Jmax.
205 // Taken from Harley and Baldocchi (1995, PCE)
206 double hkin;// = 200000.0; // enthalpy term
207 double skin;// = 710.0; // entropy term
208 double ejm;// = 55000.0; // activation energy for electron transport
209 double evc;// = 55000.0; // activation energy for carboxylation
210
211 // Enzyme constants & partial pressure of O2 and CO2 Michaelis-Menten K values.
212 // From survey of literature.
213 double kc25;// = 274.6; // kinetic coef for CO2 at 25 C, microbars
214 double ko25;// = 419.8; // kinetic coef for O2 at 25C, millibars
215 double o2;// = 210.0; // oxygen concentration umol mol-1
216
217 /* tau is computed on the basis of the Specificity factor (102.33)
218 times Kco2/Kh2o (28.38) to convert for value in solutiont to that based in air
219 The old value was 2321.1. */
220 double tau25;// = 2904.12; // New value for Quercus robor from Balaguer et al. 1996
221
222 // Arrhenius constants
223 // Eact for Michaelis-Menten const. for KC, KO and dark respiration
224 // These values are from Harley
225 double ekc;// = 80500.0; // Activation energy for K of CO2; J mol-1
226 double eko;// = 14500.0; // Activation energy for K of O2
227 double erd;// = 38000.0; // activation energy for dark respiration, eg Q10=2
228 double ektau;// = -29000.0;
229 double toptvc;// = 311.0; // optimum temperature for maximum carboxylation
230 double toptjm;// = 311.0; // optimum temperature for maximum electron transport
231 double eabole;// = 45162; // activation energy for bole respiration for Q10 = 2.02
232
233 // leaf energy balance
234 double ep;// = 0.98; // emissivity of leaves
235 double epm1;// = 0.02; // 1- ep
236 double epsoil;// = 0.98; // Emissivity of soil
237 double epsigma;// = 5.5566e-8; // ep*sigma
238 double epsigma2;// = 11.1132e-8; // 2.0 * ep * sigma
239 double epsigma4;// = 22.2264e-8; // 4.0 * ep * sigma
240 double epsigma6;// = 33.3396e-8; // 6.0 * ep * sigma
241 double epsigma8;// = 44.448e-8; // 8.0 * ep * sigma
242 double epsigma12;//= 66.6792e-8; // 12.0 * ep * sigma
243 double betfact;// = 1.5; // multiplication factor for aerodynamic sheltering, based on work by Grace and Wilson
244
245 // Ball-Berry stomatal coefficient for stomatal conductance
246 double kball;// = 9.5;
247
248 // intercept of Ball-Berry model, mol m-2 s-1
249 double bprime;// = 0.0175; // intercept for H2O
250 double bprime16;// = bprime/1.577; // 0.0109375; ntercept for CO2, bprime16 = bprime / 1.577;
251
252 // Minimum stomatal resistance, s m-1.
253 double rsm;// = 145.0;
254 double brs;// = 60.0; // curvature coeffient for light response
255
256 // leaf quantum yield, electrons
257 double qalpha;// = 0.22;
258 double qalpha2;// = 0.0484; // qalpha squared, qalpha2 = pow(qalpha, 2.0);
259
260 // leaf clumping factor
261 double markov;// = 0.84;
262
263 // Leaf dimension. geometric mean of length and width (m)
264 double lleaf;// = 0.1; // leaf length, m
265
266 // geografical parameters
267 double latitude;//51.07; // changed 51 degree 04'45,14'' N, latitude
268 double longitude;//10.45; // changed 10 degree 27'07,83'' E, longitude
269 double zone;//-1.0; // changed, minus one hour delay from GMT
270
271 // soil parameters
272 int n_soil;// = 10; // number of soil layers
273 double z_litter; // depth of litter layer
274 double density_litter; // bulk density of litter
275 double density_soil;// = 1.33; // bulk density of soil
276
277 // general Parameters
278
279 // physical parameters
280 double PI; // pi
281 double PI180; // pi divided by 180, radians per degree
282 double PI2; // 2 times pi
283 double PI9; // 0.9 times pi
284
285 // Universal constants
286 double rugc; // universal gas constant (J mole-1 K-1)
287 double rgc1000; // gas constant times 1000.
288 double dtk; // degree Kelvin at 0 degree centigrade
289
290 // Constants for leaf energy balance
291 double sigma; // Stefan-Boltzmann constant W M-2 K-4
292 double cp; // specific heat of air, J KG-1 K-1
293 double mass_air; // molecular mass of dry air, g mole-1
294 double mass_CO2; // molecular mass of CO2, g mole-1
295 double mass_13CO2; // molecular mass of 13CO2, g mole-1
296 double fumol; // conversion factor from Watt in umol PAR (Cox et al. 1998: 4.57) // rg
297 double fpar; // conversion factor for global into PA-radiation (Monteith 1965, Meek et al. 1984)
298
299 // Diffusivity values for 273 K and 1013 mb (STP) using values from Massman (1998) Atmos Environment
300 // These values are for diffusion in air. When used these values must be adjusted for
301 // temperature and pressure
302 // nu, Molecular viscosity
303 double nuvisc; // mm2 s-1
304 double nnu; // m2 s-1
305
306 // Diffusivity of CO2
307 double dc; // mm2 s-1
308 double ddc; // m2 s-1
309
310 // Diffusivity of heat
311 double dh; // mm2 s-1
312 double ddh; // m2 s-1
313
314 // Diffusivity of water vapor
315 double dv; // mm2 s-1
316 double ddv; // m2 s-1
317
318 // Diffusivity of ozone
319 double do3; // mm2 s-1
320 double ddo3; // m2 s-1
321
322 // Diffusivity of 13CO2
323 // Isotope ratio of PeeDee Belimdite standard (PDB)
324 // redefined as 13C/(12C+13C) for PDB, Tans et al 93
325 double Rpdb_CO2;
326
327 // Isotope ratio of PeeDee Belimdite standard (PDB)
328 // defined as 13C/12C, from Farquhar
329 double Rpdb_12C;
330
331 // Start and End of model run in year, format dddhhmm
332 /* rg: not needed
333 long start_run = 0000000;//start on doy 000 at 0000 hours
334 long end_run = 3660000;//end on doy 366 at 0000 hours
335 long start_profiles = 1400000;//start model profiles on doy 140 at 0000 hours
336 long end_profiles = 1450000;//end model profiles on doy 145 at 0000 hours
337 */
338
339 other_globals ()
340 {
341 PI = 3.14159265; // pi
342 PI180 = 0.017453292; // pi divided by 180, radians per degree
343 PI2 = 6.283185307; // 2 times pi
344 PI9 = 2.864788976; // 0.9 times pi
345
346 // Universal constants
347 rugc = 8.3143; // universal gas constant (J mole-1 K-1)
348 rgc1000 = 8314.3; // gas constant times 1000.
349 dtk = 273.16; // degree Kelvin at 0 degree centigrade
350
351 // Constants for leaf energy balance
352 sigma = 5.67e-08; // Stefan-Boltzmann constant W M-2 K-4
353 cp = 1005.; // specific heat of air, J KG-1 K-1
354 mass_air = 28.97; // molecular mass of dry air, g mole-1
355 mass_CO2 = 44.; // molecular mass of CO2, g mole-1
356 mass_13CO2 = 45.; // molecular mass of 13CO2, g mole-1
357 fumol = 4.57; // conversion factor from Watt in umol PAR (Cox et al. 1998: 4.57) // rg
358 fpar = 0.45; // conversion factor for global into PA-radiation (Monteith 1965, Meek et al. 1984)
359
360 // Diffusivity values for 273 K and 1013 mb (STP) using values from Massman (1998) Atmos Environment
361 // These values are for diffusion in air. When used these values must be adjusted for
362 // temperature and pressure
363 // nu, Molecular viscosity
364 nuvisc = 13.27; // mm2 s-1
365 nnu = 0.00001327; // m2 s-1
366
367 // Diffusivity of CO2
368 dc = 13.81; // mm2 s-1
369 ddc = 0.00001381; // m2 s-1
370
371 // Diffusivity of heat
372 dh = 18.69; // mm2 s-1
373 ddh = 0.00001869; // m2 s-1
374
375 // Diffusivity of water vapor
376 dv = 21.78; // mm2 s-1
377 ddv = 0.00002178; // m2 s-1
378
379 // Diffusivity of ozone
380 do3 = 14.44; // mm2 s-1
381 ddo3 = 0.00001444; // m2 s-1
382
383 // Diffusivity of 13CO2
384 // Isotope ratio of PeeDee Belimdite standard (PDB)
385 // redefined as 13C/(12C+13C) for PDB, Tans et al 93
386 Rpdb_CO2 = 0.01115;
387
388 // Isotope ratio of PeeDee Belimdite standard (PDB)
389 // defined as 13C/12C, from Farquhar
390 Rpdb_12C = 0.01124;
391 }
392};
393
394
395// radiation variables, visible, near infrared and infrared
396struct solar_radiation_variables{
397 // profiles of the probabilities of sun and shade leaves
398 double prob_beam[sze]; // probability of beam or sunlit fraction
399 double prob_sh[sze]; // probability of shade
400 double ir_dn[sze]; // downward directed infrared radiation, W m-2
401 double ir_up[sze]; // upward directed infrared radiation. W m-2
402 // inputs of near infrared radiation and components, W m-2
403 double nir_beam; // incoming beam component near infrared radiation, W m-2
404 double nir_diffuse; // incoming diffuse component near infrared radiation, W m-2
405 double nir_total; // incoming total near infrared radiaion, W m-2
406 // computed profiles of near infrared radiation, W m-2
407 double nir_dn[sze]; // downward scattered near infrared radiation
408 double nir_up[sze]; // upward scattered near infrared radiation
409 double nir_sun[sze]; // near infrared radiation on sunlit fraction of layer
410 double nir_sh[sze]; // near infrared radiation on shaded fraction of layer
411 double beam_flux_nir[sze]; // flux density of direct near infrared radiation
412 // leaf and soil optical properities of near infrared radiation
413 double nir_reflect; // leaf reflectance in the near infrared
414 double nir_trans; // leaf transmittance in the near infrared
415 double nir_soil_refl; // soil reflectance in the near infrared
416 double nir_absorbed; // leaf absorptance in the near infrared
417 // inputs of visible light, PAR, W m-2
418 double par_diffuse; // diffuse component of incoming PAR, parin
419 double par_beam; // beam component of incoming PAR, parin
420 // computed profiles of visible radiation, PAR, W m-2
421 double par_shade[sze]; // PAR on shaded fraction of layer
422 double par_sun[sze]; // PAR on sunlit fraction of layer, beam and diffuse
423 double beam_flux_par[sze]; // PAR in the beam component
424 double par_down[sze]; // downward scattered PAR
425 double par_up[sze]; // upward scattered PAR
426 // flux densities of visible quanta on sun and shade leaves for photosynthesis calculations, micromoles m-2 s-1
427 double quantum_sun[sze]; // par on sunlit leaves
428 double quantum_sh[sze]; // par on shaded leaves
429 // optical properties of leaves and soil for PAR
430 double par_absorbed; // PAR leaf absorptance
431 double par_reflect; // PAR leaf reflectance
432 double par_trans; // PAR leaf transmittance
433 double par_soil_refl; // PAR soil reflectance
434 // Net radiation profiles, W m-2
435 double rnet_sun[sze]; // net radiation flux density on sunlit fraction of layer
436 double rnet_sh[sze]; // net radiation flux density on shade fraction of layer
437
438 double exxpdir[sze]; // exponential transmittance of diffuse radiation through a layer
439 double beta_rad; // solar elevation angle, radians
440 double sine_beta; // sine of beta
441 double beta_deg; // solar elevation angle, degrees
442 double ratrad; // radiation ratio to detect cloud amount
443 double ratradnoon; // radiation ratio at noon for guestimating cloud amount at night
444
445 solar_radiation_variables ()
446 {
447 nir_beam = 0.0;
448 nir_diffuse = 0.0;
449 nir_total = 0.0;
450 nir_reflect = 0.0;
451 nir_trans = 0.0;
452 nir_soil_refl = 0.0;
453 nir_absorbed = 0.0;
454 par_diffuse = 0.0;
455 par_beam = 0.0;
456 par_absorbed = 0.0;
457 par_reflect = 0.0;
458 par_trans = 0.0;
459 par_soil_refl = 0.0;
460 beta_rad = 0.0;
461 sine_beta = 0.0;
462 beta_deg = 0.0;
463 ratrad = 0.0;
464 ratradnoon = 0.0;
465
466 for (int n1 = 0; n1 <= sze; n1++){
467 prob_beam[n1] = 0.0;
468 prob_sh[n1] = 0.0 ;
469 ir_dn[n1] = 0.0;
470 ir_up[n1] = 0.0;
471 nir_dn[n1] = 0.0;
472 nir_up[n1] = 0.0;
473 nir_sun[n1] = 0.0;
474 nir_sh[n1] = 0.0;
475 beam_flux_nir[n1] = 0.0;
476 par_shade[n1] = 0.0;
477 par_sun[n1] = 0.0;
478 beam_flux_par[n1] = 0.0;
479 par_down[n1] = 0.0;
480 par_up[n1] = 0.0;
481 rnet_sun[n1] = 0.0;
482 rnet_sh[n1] = 0.0;
483 exxpdir[n1] = 0.0;
484 }
485 }
486};
487
488// physical properties of soil abd soil energy balance variables
489struct soil_variables{
490 // soil properties
491 double z_soil[soilsze]; // depth increments of soil layers
492 double bulk_density[soilsze]; // bulk density of soil
493 double T_soilIni[soilsze]; // initial soil temperature
494 double T_soil[soilsze]; // soil temperature
495 double k_conductivity_soil[soilsze]; // thermal conductivity of soil *dz
496 double k_conductivity_soil_norm[soilsze]; // thermal conductivity of soil
497 double cp_soil[soilsze]; // specific heat of soil layer (J m-2 K-1)
498 double cp_soil_norm[soilsze]; // specific heat of soil (J m-3 K-1)
499 double T_Kelvin; // soil surface temperature in Kelvin
500 double T_air; // air temperature above soil, C
501 double sfc_temperature; // soil surface temperature in C
502 double Temp_ref; // reference soil temperature, annual mean, C
503 double Temp_amp; // amplitude of soil temperature, C
504 double resistance_h2o; // soil resistance for water vapor
505 //rg double water_content_sfc; // volumetric water content of soil surface
506 double water_content_15cm; // vol water content of 15 cm soil layer
507 double water_content_litter; // vol water content of litter
508 double water_content[soilsze];// vol water content of litter
509 double T_base; // base soil temperature
510 double T_15cm; // soil temperature at 15 cm
511 double amplitude; // amplitude of soil temperature cycle
512 double clay_fraction; // clay fraction
513 // soil energy flux densities, W m-2
514 double lout; // longwave efflux from soil
515 double evap; // soil evaporation
516 double heat; // soil sensible heat flux density
517 double rnet; // net radiation budget of the soil
518 double gsoil; // soil heat flux density
519 // soil CO2 respiratory efflux
520 double respiration_mole; // soil respiration, micromol m-2 s-1
521 double respiration_mg; // soil respiration, mg CO2 m-2 s-1
522 double base_respiration; // base rate of soil respiration, micromol m-2 s-1
523 double resp_13; // respiration of 13C micromole m-2 s-1
524 double SR_ref_temp; // refernce temperature for soil respiration, degC
525 double dtsoi; // time step for water balance calculations, s
526 double dtbal; // time step for energy balance calculations, s
527 long int mtime; // number of time steps for energy balance calculatoin per outer time step
528 // from LSM
529 double swp[soilsze]; //soil water potential (kPa)
530 double swp_mm[soilsze]; //soil water potential (mm)
531 double r[soilsze];
532 double a[soilsze];
533 double b[soilsze];
534 double c[soilsze];
535 double dwat[soilsze];
536 double watmin[soilsze]; //minimum volumetric soil water content // rg
537 double watsat[soilsze]; //saturated volumetric soil water content (porosity)
538 double smpsat[soilsze]; //soil matrix potential at saturation (mm)
539 double bch[soilsze]; //clapp and hornberger "b"
540 double hksat[soilsze]; //hydraulic conductivity at saturation (mm h2o/s)
541 double root[soilsze]; //relative root abundance (0 to 1)
542 double h2osoi[soilsze]; //volumetric soil water content
543 double qdrai; //sub-surface runoff (mm h2o s-1)
544 double qinfl; //infiltration rate, mm h2o s-1
545 double clay[soilsze]; //clay fraction in each layer
546 double sand[soilsze]; //sand fraction in each layer
547 double soil_mm; //water content in the soil, mm m-2
548 double qtran; //plant transpiration, mm s-1
549 double qin[soilsze]; //water input into each layer, mm s-1
550 double qout[soilsze]; //water output from each layer, mm s-1
551
552 soil_variables ()
553 {
554 T_Kelvin = 0.0;
555 T_air = 0.0;
556 sfc_temperature = 0.0;
557 Temp_ref = 0.0;
558 Temp_amp = 0.0;
559 resistance_h2o = 0.0;
560 water_content_15cm = 0.0;
561 water_content_litter = 0.0;
562 T_base = 0.0;
563 T_15cm = 0.0;
564 amplitude = 0.0;
565 clay_fraction = 0.0;
566 lout = 0.0;
567 evap = 0.0;
568 heat = 0.0;
569 rnet = 0.0;
570 gsoil = 0.0;
571 respiration_mole = 0.0;
572 respiration_mg = 0.0;
573 base_respiration = 0.0;
574 resp_13 = 0.0;
575 SR_ref_temp = 0.0;
576 dtsoi = 0.0;
577 dtbal = 0.0;
578 qdrai = 0.0;
579 qinfl = 0.0;
580 soil_mm = 0.0;
581 qtran = 0.0;
582 }
583};
584
585// Structure for Profile information,fluxes and concentrations
586struct profile{
587 // microclimate profiles
588 double tair[sze3]; // air temp (C)
589 double tair_filter[sze3]; // numerical filter of Tair
590 double u[sze3]; // wind speed (m/s)
591 double rhov_air[sze3]; // water vapor density
592 double rhov_filter[sze3]; // numerical filter of rhov_air
593 double co2_air[sze3]; // co2 concentration (ppm)
594 // canopy structure profiles
595 double dLAIdz[sze]; // leaf area index of layer (m2/m2)
596 double ht[sze]; // layer height (m)
597 double dPAIdz[sze]; // plant area index of layer
598 double Gfunc_solar[sze]; // leaf-sun direction cosine function
599 double Gfunc_sky[sze][szeang]; // leaf-sky sector direction cosine function
600 double tleaf[sze]; // leaf temperature per layer (sun and shade)
601 double isopreneflux[sze]; // isoprene flux per layer (sun and shade)
602 double vcmax[sze]; // vcmax in per layer
603 double jmax[sze]; // jmax per layer
604 double rd[sze]; // rd per layer
605 // variables for 13C isotopes
606 double c13cnc[sze3]; // concentration of 13C
607 double c12cnc[sze3]; // concentration of 12C
608 double sour13co2[sze]; // source/sink strength 13C
609 double d13C[sze3]; // del 13C
610 double D13C[sze]; // photosynthetic weighted discrimination, 13D
611 double D13C_long[sze]; // photosynthetic weighted discrimination, 13D
612 double d13Cair[sze3]; // del 13C of the air
613 double d13Cplant[sze]; // del 13C of the plant
614 double R13_12_air[sze3]; // 13C/12C ratio of the air
615 double Rplant_sun[sze]; // ratio of discriminated 13C in sunlit leaves
616 double Rplant_shd[sze]; // ratio of discriminated 13C in shaded leaves
617 double Rresp[sze]; // discriminated 13C for later respiration, Ps weight
618 double Rresp_sum[sze]; // summing of Ps weighted values for daily ave
619 double Rresp_ave[sze]; // previous days discriminated 13C for later respiration
620 int cnt_Rresp[sze];
621 double recycle[sze3]; // fraction of recycled CO2, after Yakir and Sternberg
622 // source/sink strengths
623 double source_co2[sze]; // source/sink strength of CO2
624 // fluxes for total layer, with sun/shade fractions
625 double dPsdz_mg[sze]; // layer photosynthesis (mg layer-1 s-1)
626 double dPsdz[sze]; // layer photosynthesis (umol layer-1 s-1)
627 double dHdz[sze]; // layer sensible heat flux
628 double dLEdz[sze]; // layer latent heat flux
629 double dRNdz[sze]; // layer net radiation flux
630 double dRESPdz[sze]; // layer respiration
631 double dStomCondz[sze]; // layer stomatal conductance
632 // sun leaf variables
633 double sun_frac[sze]; // sun leaf fraction
634 double sun_tleaf[sze]; // leaf temp (C)
635 double sun_A[sze]; // layer A flux for sun only (umol m-2 s-1)
636 double sun_gs[sze]; // stomatal conductance (m s-1)
637 double sun_rs[sze]; // stomatal resistance to H2O (s/m)
638 double sun_rbh[sze]; // boundary layer resistance to heat (s/m)
639 double sun_rbv[sze]; // boundary layer resistance to H2O (s/m)
640 double sun_rbco2[sze]; // boundary layer resistance to CO2 (s/m)
641 double sun_cs[sze]; // Cs on sun leaves (CO2 mixing ratio on leaf surface
642 double sun_ci[sze]; // Ci on sun leaves
643 double sun_cc[sze]; // Cc (CO2 mixing ratio at site of carboxylation) on sun leaves
644 double sun_D13[sze]; // discrimination 13C
645 double sun_D13_long[sze]; // discrimination 13C
646 double sun_ccca[sze]; // Cc/Ca on sunlit leaves
647 double sun_cica[sze]; // Ci/Ca on sunlit leaves
648 double sun_lai[sze]; // sunlit lai of layer
649 double sun_T_filter[sze]; // filtered sunlit temperature
650 double sun_wj[sze]; // electron transport rate of Ps for sun leaves
651 double sun_wc[sze]; // carboxylatio velocity for sun leaves
652 double sun_resp[sze]; // respiration
653 double sun_isopreneflux[sze];// isoprene flux per layer for sunleaves
654 double iso_sun[sze]; // isoprene flux per leaf area in the sun
655 // shade leaf variables
656 double shd_frac[sze]; // shade leaf fraction
657 double shd_tleaf[sze]; // temperature of shaded leaves
658 double shd_A[sze]; // photosynthesis of shaded leaves
659 double shd_gs[sze]; // stomatal conductance of shade leaves
660 double shd_rs[sze]; // stomatal resistance of shaded leaves
661 double shd_rbh[sze]; // boundary layer resistance for heat on shade leaves
662 double shd_rbv[sze]; // boundary layer resistance for vapor on shade leaves
663 double shd_rbco2[sze]; // boundary layer resistance for CO2 on shade leaves
664 double shd_cs[sze]; // Cs on shade leaves (CO2 mixing ratio on leaf surface
665 double shd_ci[sze]; // Ci on shade leaves
666 double shd_cc[sze]; // Cc (CO2 mixing ratio at site of carboxylation) on sun leaves
667 double shd_D13[sze]; // del 13C on shaded leaves
668 double shd_D13_long[sze]; // discrimination 13C
669 double shd_ccca[sze]; // Cc/Ca ratio on shaded leaves
670 double shd_cica[sze]; // Ci/Ca ratio on shaded leaves
671 double shd_lai[sze]; // shaded lai of layer
672 double shd_T_filter[sze]; // previous temperature
673 double shd_wc[sze]; // carboxylation rate for shade leaves
674 double shd_wj[sze]; // electron transport rate for shade leaves
675 double shd_resp[sze]; // respiration
676 double shd_isopreneflux[sze];// isoprene flux per layer for shade leaves
677 double iso_shd[sze]; // isoprene flux per leaf area in the shade
678
679 profile ()
680 {
681 for (int n0 = 0; n0 <= sze3; n0++){
682 tair[n0] = 0.0;
683 tair_filter[n0] = 0.0 ;
684 u[n0] = 0.0;
685 rhov_air[n0] = 0.0;
686 rhov_filter[n0] = 0.0;
687 co2_air[n0] = 0.0 ;
688 c13cnc[n0] = 0.0;
689 c12cnc[n0] = 0.0;
690 d13C[n0] = 0.0;
691 d13Cair[n0] = 0.0 ;
692 R13_12_air[n0] = 0.0;
693 recycle[n0] = 0.0;
694 }
695
696 for (int n1 = 0; n1 <= sze; n1++){
697 dLAIdz[n1] = 0.0;
698 ht[n1] = 0.0 ;
699 dPAIdz[n1] = 0.0;
700 Gfunc_solar[n1] = 0.0;
701 tleaf[n1] = 0.0;
702 isopreneflux[n1] = 0.0;
703 vcmax[n1] = 0.0;
704 jmax[n1] = 0.0;
705 rd[n1] = 0.0 ;
706 D13C[n1] = 0.0;
707 D13C_long[n1] = 0.0;
708 d13Cplant[n1] = 0.0;
709 Rplant_shd[n1] = 0.0;
710 Rresp[n1] = 0.0;
711 Rresp_ave[n1] = 0.0;
712 source_co2[n1] = 0.0;
713 cnt_Rresp[n1] = 0;
714 dPsdz_mg[n1] = 0.0;
715 dPsdz[n1] = 0.0 ;
716 dHdz[n1] = 0.0;
717 dLEdz[n1] = 0.0;
718 dRNdz[n1] = 0.0;
719 dRESPdz[n1] = 0.0;
720 dStomCondz[n1] = 0.0;
721 sun_frac[n1] = 0.0;
722 sun_tleaf[n1] = 0.0 ;
723 sun_A[n1] = 0.0;
724 sun_gs[n1] = 0.0;
725 sun_rs[n1] = 0.0;
726 sun_rbh[n1] = 0.0;
727 sun_rbco2[n1] = 0.0;
728 sun_cs[n1] = 0.0;
729 sun_ci[n1] = 0.0;
730 sun_cc[n1] = 0.0;
731 sun_D13[n1] = 0.0;
732 sun_D13_long[n1] = 0.0;
733 sun_ccca[n1] = 0.0;
734 sun_cica[n1] = 0.0;
735 sun_lai[n1] = 0.0;
736 sun_T_filter[n1] = 0.0;
737 sun_wc[n1] = 0.0;
738 sun_wj[n1] = 0.0;
739 sun_resp[n1] = 0.0;
740 sun_isopreneflux[n1] = 0.0;
741 iso_sun[n1] = 0.0;
742 shd_frac[n1] = 0.0;
743 shd_tleaf[n1] = 0.0 ;
744 shd_A[n1] = 0.0;
745 shd_gs[n1] = 0.0;
746 shd_rs[n1] = 0.0;
747 shd_rbh[n1] = 0.0;
748 shd_rbco2[n1] = 0.0;
749 shd_cs[n1] = 0.0;
750 shd_ci[n1] = 0.0;
751 shd_cc[n1] = 0.0;
752 shd_D13[n1] = 0.0;
753 shd_D13_long[n1] = 0.0;
754 shd_ccca[n1] = 0.0;
755 shd_cica[n1] = 0.0;
756 shd_lai[n1] = 0.0;
757 shd_T_filter[n1] = 0.0;
758 shd_wc[n1] = 0.0;
759 shd_wj[n1] = 0.0;
760 shd_resp[n1] = 0.0;
761 shd_isopreneflux[n1] = 0.0;
762 iso_shd[n1] = 0.0;
763 }
764 }
765};
766
767/*rg
768struct switch_variable{ //rg: changed from double to int
769int a;// sets soil respiration reference temperature to (1) input soil temperature or (0) modeled soil temperature at 5 cm
770int b;//
771int c;//
772int d;//
773int e;//
774int f;//
775int g;//
776int h;//
777int i;//
778int j;//
779} set_switch;
780*/
781
782struct isotope_variable{
783 double delta_soil; // delta13C of soil respiration
784 double da_m_day; // slope of daytime regression deltaCa*Ca=m*Ca+b
785 double da_b_day; // intercept of daytime regression deltaCa*Ca=m*Ca+b
786 double da_m_night; // slope of nighttime regression deltaCa*Ca=m*Ca+b
787 double da_b_night; // intercept of nighttime regression deltaCa*Ca=m*Ca+b
788};
789
790
791
792/* WaterCycleCanoak
793 * The class WaterCycleCanoak contains all Methods and data from the WaterCycleCanoak Module.
794 *
795 *
796*/
797
798class WaterCycleCANOAK : public WaterCycleModule {
799 public:
800 struct input_variables input;
801 struct flux_variables flux;
802 struct time_variables time_var;
803 struct meteorology met;
804 struct surface_resistances sfc_res;
805 struct factors fact;
806 struct bole_respiration_structure bole;
807 struct canopy_architecture canopy;
808 struct non_dimensional_variables non_dim;
809 struct boundary_layer_resistances bound_layer_res;
810 struct solar_radiation_variables solar;
811 struct soil_variables soil;
812 struct profile prof;
813 struct isotope_variable Cisotope;
814 struct other_globals og;
815 // VARIABLES and PARAMETERS
816 int ts,nd,ndStart,ndEnd,tsStart,tsEnd,tsMax,ny,doy;
817 bool leaf_out,leaf_full,leaf_fall,leaf_fall_complete;
818 double ts_co2_concentration,ts_airtemperature,ts_shortwaveradiance,ts_airpressure,ts_watervaporsaturationdeficit,ts_precipitation,ts_windspeed;
819 double ZONE,LAT,LONG,TAMP,temp_ny;
820 double TAU,RCMIN,VCMAX25,JMAX25,LLEAF,BFACT1,BFACT2,FSTO,
821 RCPAR1,TCPAR1,SRPAR1,RCNIR1,TCNIR1,SRNIR1,
822 RCPAR2,TCPAR2,SRPAR2,RCNIR2,TCNIR2,SRNIR2,
823 RCPAR3,TCPAR3,SRPAR3,RCNIR3,TCNIR3,SRNIR3,
824 RCPAR4,TCPAR4,SRPAR4,RCNIR4,TCNIR4,SRNIR4;
825 double ht,lai;
826 double h_litter,clayf,bdens;
827 double fH2O,laiPot,wc_sfc,wc_15cm,wc_litter,t_firstmin;
828
829 WaterCycleCANOAK(ModelCore & modelCore, unsigned int timeInterval) : WaterCycleModule(modelCore, timeInterval),
830 modelCore_(modelCore),
831 time_(modelCore.getLandscape().getTime()),
832 si_(modelCore.getInput().getSite()),
833 sipar_(modelCore.getInput().getParameter().getSitePar()),
834 sppar_(modelCore.getInput().getParameter().getSpeciesPar()),
835 ac_(modelCore.getSubstates().getAc()),
836 mc_(const_cast<MicroClimate&>(modelCore.getSubstates().getMc())),
837 sc_(modelCore.getSubstates().getSc()),
838 wc_(const_cast<WaterCycle &>(modelCore.getSubstates().getWc())),
839 ph_(modelCore.getSubstates().getPh()),
840 vs_(modelCore.getSubstates().getVs())
841 {
842 this->setName("WaterCycleCANOAK");
843 }
844
846 ~WaterCycleCANOAK() {
847
848 }
849
850 void solve() {
851 this->run_(
852 sipar_, /*si_.runWC, si_. runPH2,*/ time_.ny(), time_.nd() , time_.ndEnd(), time_.ts(), time_.tsMax(), time_.nd(),
853 /*std::string SITEname,*/
854 time_.ndStart(), time_.tsStart(), time_.tsMax(),
855 si_.vtMax, si_.flMax, si_.slMax, si_.slFloor,
856 si_.name_vt, si_.ZONE, si_.LONG,si_.LAT, si_.TAMP,
857 si_.wcMax_sl, si_.wcMin_sl, si_.h_sl, si_.h_fl,
858 si_.depth_vt, si_.hMax_vt, si_.hMin_vt, si_.clay_sl,
859 si_.dens_sl, si_.poro_sl, si_.stonef_sl, si_.fcOrg_sl,
860 mc_.temp_ny, ac_.ts_co2_concentration, mc_.ts_airpressure, mc_.ts_airtemperature,
861 mc_.ts_shortwaveradiance, mc_.ts_watervaporsaturationdeficit, mc_.ts_windspeed, mc_.ts_precipitation,
862 ph_.lai_vt, ph_.lai_fl, ph_.dvsFlush_vt, ph_.dvsMort_vt,
863 vs_.f_vt, vs_.laiMax_vt, vs_.fFrt_vtsl,
864 mc_.temp_a, mc_.temp_l,mc_.temp_sl, ac_.ts_co2_concentration_fl,
865 wc_.wc_fl, wc_.wc_sl, wc_.ice_sl, wc_.fH2O_vt,
866 wc_.fH2O_sl, wc_.flux_sl, wc_.evapot, wc_.transp,
867 wc_.runoff, wc_.percol, wc_.evacep, wc_.evasoil,
868 wc_.throughf,mc_.dispersion);
869 }
870
871 void CANOAK(double *vpd_layer,double *disper);
872
873 void ROUGH
874 (double HEIGHT,double LAI,double SAI,
875 double &DISP);
876
877 void DISPER
878 (std::string SITEname,double HH,double DD,
879 double *dispersion);
880
881 void INPUT_DATA
882 (int doy,int ts,int tsMax,
883 double ts_airtemperature,double ts_shortwaveradiance,double ts_watervaporsaturationdeficit,
884 double ts_windspeed,double ts_precipitation,double ts_co2_concentration,double ts_airpressure,
885 double t_firstmin,double wc_15cm,double wc_litter);
886
887 void SOIL_RESPIRATION();
888
889 //void BOLE_RESPIRATION();
890 void BOLE_RESPIRATION
891 (double BFACT1,double BFACT2);
892
893 void ENERGY_AND_CARBON_FLUXES(double FSTO,double *vpd_layer);
894
895 void DIFFUSE_DIRECT_RADIATION();
896
897 void RNET();
898
899 void PAR();
900
901 void NIR();
902
903 double SKY_IR(double a);
904
905 void IRFLUX();
906
907 void G_FUNC_DIFFUSE();
908
909 void GFUNC();
910
911 void ANGLE();
912
913 double GAMMAF(double a);
914
915 //void LAI_TIME();
916 void LAI_TIME(double lai);
917
918 void FREQ(double a);
919
920 void STOMATA();
921
922 double TEMP_FUNC(double a,double b,double c,double d, double e);
923
924 double TBOLTZ(double a,double b, double c, double d);
925
926 void BOUNDARY_RESISTANCE(double a, double b);
927
928 void FRICTION_VELOCITY(double &metH_old);
929
930 void CONC(double *ap1, double *ap2, double c, double d, double e);
931
932 void ENERGY_BALANCE
933 (double b,double d,double e,
934 double f,double g,
935 double &apc,double &h,double &i,
936 double &j,double &vpd_leaf);
937
938 double SFC_VPD(double a, double b, double c);
939
940 //rg deleted: void SET_SOIL();
941 //rg deleted: void SET_SOIL_TEMP();
942
943 void SOIL_ENERGY_BALANCE();
944
945 double SOIL_SFC_RESISTANCE(double a);
946
947 void ISOTOPE();
948
949 //rg void INPUT_DATA();
950 //rg deleted: void FILENAMES();
951 //rg deleted void ZERO();
952 //rg deleted: void READ_PARAMETER();
953 //rg void SET_PARAMETER();
954 //rg deleted: int FileCopy(const char *src,const char *dst );
955 //rg deleted: void SOIL_MOISTURE_PROFIL();
956
957 void SET_PARAMETER
958 (double TAU,double RCMIN,double VCMAX25,double JMAX25,double LLEAF,
959 double RCPAR1,double TCPAR1,double SRPAR1,
960 double RCNIR1,double TCNIR1,double SRNIR1,
961 double RCPAR2,double TCPAR2,double SRPAR2,
962 double RCNIR2,double TCNIR2,double SRNIR2,
963 double RCPAR3,double TCPAR3,double SRPAR3,
964 double RCNIR3,double TCNIR3,double SRNIR3,
965 double RCPAR4,double TCPAR4,double SRPAR4,
966 double RCNIR4,double TCNIR4,double SRNIR4,
967 int ndStart,int ndEnd,int tsStart,int tsEnd,int tsMax,int doy,
968 double LAT,double LONG,double ZONE,double TAMP,
969 double temp_ny,double ht,double laiPot,
970 double h_litter,double bdens,double clayf,
971 double wc_15cm,double wc_litter);
972
973 //int SOILH2O();
974 int SOILH2O(int tsMax,double lai);
975
976 int TRIDIA(double &checkChange);
977
978
979 private:
980 ModelCore & modelCore_;
981 Timer const & time_;
982 Site const & si_;
983 SitePar const & sipar_;
984 SpeciesPar const & sppar_;
985
986 AirChemistry const & ac_;
987 MicroClimate & mc_;
988 SoilChemistry const & sc_;
989 WaterCycle & wc_;
990 Physiology const & ph_;
991 VegStructure const & vs_;
992
993 void run_(
994 const SitePar &psi,
995 /*unsigned int runWC, unsigned int runPH2,*/
996 unsigned int ny, unsigned int nd, unsigned int ndEnd, unsigned int ts, unsigned int tsEnd, unsigned int doy,
997 /*std::string SITEname,*/ unsigned int ndStart, unsigned int tsStart, unsigned int tsMax,
998 size_t const & vtMax, size_t const & flMax, size_t const & slMax, size_t const & slFloor,
999 std::string * const& name_vt, const int &ZONE, const double &LONG, const double &LAT,const double &TAMP,
1000 double *const& wcMax_sl,double *const& wcMin_sl,double *const& h_sl,double *const& h_fl,
1001 double *const& depth_vt,double *const& hMax_vt,double *const& hMin_vt,double *const& clay_sl,
1002 double *const& dens_sl,double *const& poro_sl,double *const& stonef_sl,double *const& fcOrg_sl,
1003 const double &temp_ny, const double &ts_co2_concentration, const double &ts_airpressure, const double &ts_airtemperature,
1004 const double &ts_shortwaveradiance, const double &ts_watervaporsaturationdeficit, const double &ts_windspeed, const double &ts_precipitation,
1005 double *const&lai_vt,double *const&lai_fl,double *const&dvsFlush_vt,double *const&dvsMort_vt,
1006 double *const&f_vt,double *const&laiPot_vt,double **const&fFrt_vtsl,
1007 double &temp_a, double & temp_l, double * temp_sl, const double *ts_co2_concentration_fl,
1008 double *wc_fl,double *wc_sl,double *ice_sl,double *fH2O_vt,
1009 double *fH2O_sl,double *flux_sl,double &evapot,double &transp,
1010 double &runoff, double &percol, double &evacep, double &evasoil,
1011 double &throughf, double * const & disper);
1012
1013 int Time_sec();
1014
1015 void PHOTOSYNTHESIS
1016 (int JJ,double zzz,double FSTO,double Iphoton,
1017 double cca,double tlk,double leleaf,
1018 double &rstom,double &A_mg,double &rd,double &ci,
1019 double &cs,double &wj,double &wc);
1020
1021 double DW2DZ(double Z,double HH,double ustar, double sigma_zo);
1022
1023 double RAN0(long *a);
1024
1025 double SIGMA(double Z,double HH,double sigma_zo,double sigma_h,double ustar);
1026
1027 double TL(double Z,double HH,double DD,double sigma_zo,double sigma_h,double ustar);
1028
1029 double RANDGEN(double low,double delta,long &seed);
1030
1031 double UZ(double a);
1032};
1033
1034#endif // M1D_MOD_WATERCYCLE_CANOAK_H
1035