forked from perabrahamsen/daisy-model
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbioclimate_std.C
1662 lines (1499 loc) · 64.4 KB
/
bioclimate_std.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// bioclimate_std.C --- The default biclimate model
//
// Copyright 1996-2001 Per Abrahamsen and Søren Hansen
// Copyright 2000-2001 KVL.
//
// This file is part of Daisy.
//
// Daisy is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser P //ublic License as published by
// the Free Software Foundation; either version 2.1 of the License, or
// (at your option) any later version.
//
// Daisy is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser Public License for more details.
//
// You should have received a copy of the GNU Lesser Public License
// along with Daisy; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#define BUILD_DLL
#include "bioclimate.h"
#include "block_model.h"
#include "surface.h"
#include "weather.h"
#include "plf.h"
#include "frame.h"
#include "geometry.h"
#include "soil.h"
#include "soil_heat.h"
#include "snow.h"
#include "log.h"
#include "mathlib.h"
#include "net_radiation.h"
#include "pet.h"
#include "difrad.h"
#include "raddist.h"
#include "deposition.h"
#include "svat.h"
#include "vegetation.h"
#include "litter.h"
#include "time.h"
#include "check.h"
#include "fao.h"
#include "librarian.h"
#include "treelog_store.h"
#include "resistance.h"
#include "im.h"
#include "soil_water.h"
#include <sstream>
struct BioclimateStandard : public Bioclimate
{
const Metalib& metalib;
static const double reference_height; // [m]
// Canopy State.
const long No; // No of intervals in canopy discretization.
double LAI_; // Total LAI of all crops on this column. [0-]
double cover_; // Fraction of soil covered by vegetation.
double sun_LAI_fraction_total_;// Total sun LAI fraction of all crops on this column.
std::vector<double> Height; // Height in cm of each endpoint in c.d.
std::vector<double> total_PAR_;// Total PAR of each interval of c.d. [W/m2]
std::vector<double> sun_PAR_; //Sunlit fraction of PAR of each interval of c.d.[W/m2]
std::vector<double> total_NIR_;// Total NIR of each interval of c.d. [W/m2]
std::vector<double> sun_NIR_; //Sunlit fraction of NIR of each interval of c.d[W/m2].
std::vector<double> sun_LAI_fraction_;//Fraction of sunlit LAI of each interval of cd
double shared_light_fraction_; // Fraction of field with light competition.
// Canopy Tick.
void CanopyStructure (const Vegetation&);
// External water sinks and sources.
const std::unique_ptr<NetRadiation> net_radiation;
std::unique_ptr<Pet> pet; // Potential Evapotranspiration model.
double total_ep_; // Potential evapotranspiration [mm/h]
double total_ea_; // Actual evapotranspiration [mm/h]
double direct_rain_; // Rain hitting soil directly [mm/h]
double irrigation_overhead; // Irrigation above canopy [mm/h]
double irrigation_overhead_temperature; // Water temperature [dg C]
double irrigation_surface; // Irrigation below canopy [mm/h]
double irrigation_surface_temperature; // Water temperature [dg C]
double irrigation_subsoil; // Irrigation incorporated in soil.
double irrigation_subsoil_permanent; // Irrigation incorporated in soil.
double tillage_water; // Water added to surface due to tillage [mm/h]
// Water in snowpack.
Snow snow;
double snow_ep; // Potential snow evaporation [mm/h]
double snow_ea_; // Actual snow evaporation [mm/h]
double snow_water_in; // Water entering snow pack [mm/h]
double snow_water_in_temperature; // Incoming water temperature [dg C]
double snow_water_out; // Water leaving snow pack [mm/h]
double snow_water_out_temperature; // Temperature of water leaving [dg C]
// Water intercepted on canopy.
double canopy_ep; // Potential canopy evaporation [mm/h]
double canopy_ea_; // Actual canopy evaporation [mm/h]
double canopy_water_capacity; // Pot. intercepted water on canopy [mm]
double canopy_water_storage; // Intercepted water on canopy [mm]
double canopy_water_temperature; // Temperature of incoming water [dg C]
double canopy_water_in; // Water entering canopy [mm/h]
double canopy_water_out; // Canopy drip throughfall [mm/h]
double canopy_water_bypass; // Water from above bypassing the canopy [mm/h]
double canopy_water_below; // Total water input below canopy [mm/h]
// Water intercepted on litter.
double litter_ep; // Potential litter evaporation [mm/h]
double litter_ea; // Actual litter evaporation [mm/h]
double litter_water_capacity; // Pot. intercepted water in litter [mm]
double litter_water_storage; // Intercepted water in litter [mm]
double litter_water_temperature; // Temperature of water in litter [dg C]
double litter_water_in; // Water entering litter [mm/h]
double litter_water_out; // Water leaving litter [mm/h]
// Water in pond.
double pond_ep; // Potential evaporation from pond [mm/h]
double pond_ea_; // Actual evaporation from pond [mm/h]
// Water going through soil surface.
double soil_ep; // Potential exfiltration. [mm/h]
double soil_ea_; // Actual exfiltration. [mm/h]
// Water transpirated through plant roots.
const int max_svat_iterations; // Max number of iterations with SVAT model.
const double max_svat_absolute_difference; // Max difference. [mm/h]
const double maxTdiff; // Max temperature diff. [dg C]
const double maxEdiff; // Max pressure diff [Pa]
const std::unique_ptr<SVAT> svat; // Soil Vegetation Atmosphere model.
size_t svat_fail; // Number of times the svat loop failed.
size_t svat_total; // Total number of times the svat loop entered.
bool svat_failed; // True iff the SVAT loop failed.
double crop_ep_; // Potential transpiration. [mm/h]
double crop_ea_soil; // Crop limited transpiration. [mm/h]
double crop_ea_svat; // SVAT limited transpiration. [mm/h]
double crop_ea_; // Actual transpiration. [mm/h]
double production_stress; // Stress calculated by SVAT module.
// Bioclimate canopy
double wind_speed_field_; // wind speed at screen height above the canopy [m/s]
double wind_speed_weather; // measured wind speed [m/s]
void WaterDistribution (const Time&,
Surface& surface, const Weather& weather,
Vegetation& vegetation, const Litter& litter,
const Movement&,
const Geometry&, const Soil& soil,
const SoilWater& soil_water, const SoilHeat&,
const double T_bottom,
double dt, Treelog&);
// Radiation.
static double find_albedo (const Vegetation& crops, const Litter& litter,
const Surface& surface,
const Geometry&, const Soil&, const SoilWater&);
double albedo; // Reflection factor []
const std::unique_ptr<Raddist> raddist;// Radiation distribution model.
const double min_sin_beta_; // Sinus to lowest sun angle for some models.
void RadiationDistribution (const Vegetation&, double pF, double sin_beta, Treelog&);
std::unique_ptr<Difrad> difrad; // Diffuse radiation model.
double difrad0; // Diffuse radiation above canopy [W/m2]
double absorbed_total_PAR_canopy; // Canopy absorbed PAR (sun+shade) [W/m2]
double absorbed_total_NIR_canopy; // Canopy absorbed NIR (sun+shade) [W/m2]
double absorbed_total_Long_canopy;// Canopy absorbed Long wave radiation (sun+shade)[W/m2]
double absorbed_total_PAR_soil; // Soil absorbed PAR [W/m2]
double absorbed_total_NIR_soil; // Soil absorbed NIR [W/m2]
double absorbed_total_Long_soil; // Soil absorbed Long wave radiation [W/m2]
double absorbed_sun_PAR_canopy; // Canopy absorbed PAR on sunlit leaves [W/m2]
double absorbed_sun_NIR_canopy; // Canopy absorbed NIR on sunlit leaves [W/m2]
double absorbed_sun_Long_canopy; // Canopy absorbed Long wave radiatio on sunlit leaves [W/m2]
double absorbed_shadow_PAR_canopy;// Canopy absorbed PAR on shadow leaves [W/m2]
double absorbed_shadow_NIR_canopy;// Canopy absorbed NIR on shadow leaves [W/m2]
double absorbed_shadow_Long_canopy;// Canopy absorbed Long wave radiation on shadow leaves [W/m2]
double rad_abs_soil_; // Soil absorbed radiation [W/m2]
double rad_abs_sun_canopy_; // Canopy absorbed radiatio on sunlit leaves [W/m2]
double rad_abs_shadow_canopy_; // Canopy absorbed radiation on shadow leaves [W/m2]
double incoming_Long_radiation; // Incoming longwave radiation [W/m2]
double incoming_PAR_radiation; // Incoming PAR radiation [W/m2]
double incoming_NIR_radiation; // Incoming NIR radiation [W/m2]
double incoming_Total_radiation; // Incoming radiation sum of shortwave and
// longwave [W/m2]
double sin_beta_;
// Weather.
double daily_air_temperature_; // Air temperature in canopy [dg C]
double daily_precipitation_; // From weather [mm]
double day_length_; // From weather (does not belong here) [h]
double daily_global_radiation_; // From weather [W/m2]
double global_radiation_; // From weather [W/m2]
double atmospheric_CO2_; // From weather [Pa]
double atmospheric_O2_; // From weather [Pa]
double air_pressure_; // From weather [Pa]
// Deposition.
const std::unique_ptr<Deposition> deposition; // Deposition model.
// Initialization.
const bool fixed_pet;
const bool fixed_difrad;
bool has_reference_evapotranspiration;
bool has_vapor_pressure;
bool has_wind;
bool has_min_max_temperature;
bool has_diffuse_radiation;
Weatherdata::surface_t old_surface;
// Simulation
void tick (const Time&, Surface&, const Weather&,
Vegetation&, const Litter& litter,
const Movement&, const Geometry&,
const Soil&, SoilWater&, const SoilHeat&, const double T_bottom,
double dt, Treelog&);
void clear ();
void output (Log&) const;
// Canopy.
const std::vector<double>& height () const
{ return Height; }
const std::vector<double>& PAR () const
{ return total_PAR_; }
const std::vector<double>& sun_PAR () const
{ return sun_PAR_; }
const std::vector<double>& NIR () const
{ return total_NIR_; }
const std::vector<double>& sun_NIR () const
{ return sun_NIR_; }
const std::vector<double>& sun_LAI_fraction () const
{ return sun_LAI_fraction_; }
double LAI () const
{ return LAI_; }
double cover () const
{ return cover_; }
double sun_LAI_fraction_total () const
{ return sun_LAI_fraction_total_; }
double wind_speed_field () const
{ return wind_speed_field_; }
double rad_abs_soil() const
{ return rad_abs_soil_; }
double rad_abs_sun_canopy() const
{ return rad_abs_sun_canopy_; }
double rad_abs_shadow_canopy() const
{ return rad_abs_shadow_canopy_; }
double sin_beta() const
{ return sin_beta_; }
double shared_light_fraction () const
{ return shared_light_fraction_; }
// Utils.
double LAI_sun () const // [area LEAF/FIELD]
{ return LAI () * sun_LAI_fraction_total (); }
double LAI_shadow () const // [area LEAF/FIELD]
{ return LAI () - LAI_sun (); }
double Collatz_gbw (const double tleaf) const // [m/s LEAF]
{
const double Ptot = air_pressure (); // [Pa]
const double gbw_molly = 2.0; // [mol/m^2/s LEAF]
return Resistance::molly2ms (tleaf, Ptot, gbw_molly); // [m/s LEAF]
}
// Weather.
double daily_rain_temperature () const
{ return std::max (daily_air_temperature_, 0.1); }
double daily_air_temperature () const
{ return daily_air_temperature_; }
double canopy_temperature () const
{ return svat->CanopyTemperature (); }
double canopy_vapour_pressure () const
{ return svat->CanopyVapourPressure (); }
double sun_leaf_temperature () const
{ return svat->SunLeafTemperature (); }
double shadow_leaf_temperature () const
{ return svat->ShadowLeafTemperature (); }
double sun_boundary_layer_water_conductivity () const
{
double gbw = svat->SunBoundaryLayerWaterConductivity ();
if (gbw > 0.0)
return gbw;
return Collatz_gbw (sun_leaf_temperature ()) * LAI_sun ();
}
double shadow_boundary_layer_water_conductivity () const
{
double gbw = svat->ShadowBoundaryLayerWaterConductivity ();
if (gbw > 0.0)
return gbw;
return Collatz_gbw (shadow_leaf_temperature ()) * LAI_shadow ();
}
double daily_precipitation () const
{ return daily_precipitation_; }
double day_length () const
{ return day_length_; }
double daily_global_radiation () const
{ return daily_global_radiation_; }
double global_radiation () const
{ return global_radiation_; }
double direct_rain () const
{ return direct_rain_; }
double atmospheric_CO2 () const
{ return atmospheric_CO2_; }
double atmospheric_O2 () const
{ return atmospheric_O2_; }
double air_pressure () const
{ return air_pressure_; }
// Manager.
void irrigate_overhead (double flux, double temp);
void irrigate_surface (double flux, double temp);
void irrigate_overhead (double flux);
void irrigate_surface (double flux);
void irrigate_subsoil (double flux);
void set_subsoil_irrigation (double flux);
void add_tillage_water (double amount);
// Communication with svat and external model.
double total_ep () const // [mm/h]
{ return total_ep_; }
double total_ea () const // [mm/h]
{ return total_ea_; }
double snow_ea () const // [mm/h]
{ return snow_ea_; }
double pond_ea () const // [mm/h]
{ return pond_ea_; }
double soil_ea () const // [mm/h]
{ return soil_ea_; }
double soil_surface_ea () const // [mm/h]
{ return soil_ea_ + litter_ea + pond_ea_ + snow_ea_; }
double crop_ep () const // [mm/h]
{ return crop_ep_; }
double crop_ea () const // [mm/h]
{ return crop_ea_; }
double canopy_ea () const // [mm/h]
{ return canopy_ea_; }
double min_sin_beta () const // []
{ return min_sin_beta_; }
double get_intercepted_water () const // [mm]
{ return canopy_water_storage; }
double get_snow_storage () const // [mm]
{ return snow.storage (); }
double snow_leak_rate (const double dt) const
{
const double snow_new = snow.storage ();
if (snow_new < 1e-5)
return 1.0 / dt;
if (snow_water_out < 1e-8)
return 0.0;
const double snow_old = snow_new + snow_water_out * dt;
return snow_water_out / snow_old;
}
double canopy_leak_rate (const double dt) const
{
const double epsilon_storage = 1e-8; // [mm]
const double epsilon_out = 1e-8 / dt; // [mm/h]
if (canopy_water_out < epsilon_out)
return 0.0;
if (canopy_water_storage < epsilon_storage)
return -1.0;
// Should be storage at "end of timestep" because of the way the
// rate is used.
const double clr = canopy_water_out / canopy_water_storage;
daisy_assert (clr > 0.0);
return clr;
}
double canopy_leak () const // [mm/h]
{ return canopy_water_out; }
double litter_leak_rate (const double dt) const
{
if (litter_water_out < 1e-8)
return 0.0;
if (litter_water_storage < 0.01 * litter_water_out * dt)
return 1.0 / dt;
const double litter_water_old
= litter_water_storage + litter_water_out * dt;
return litter_water_out / litter_water_old;
}
const IM& deposit () const
{ return deposition->deposit (); }
// Create.
void initialize (const Weather&, Treelog&);
void reset_weather (const Weather&, Treelog&);
bool check (const Weather&, Treelog&) const;
BioclimateStandard (const BlockModel&);
void summarize (Treelog&) const;
~BioclimateStandard ();
};
const double
BioclimateStandard::reference_height = 0.12; // [m]
void
BioclimateStandard::initialize (const Weather& weather, Treelog& msg)
{
TREELOG_MODEL (msg);
reset_weather (weather, msg);
}
void
BioclimateStandard::reset_weather (const Weather& weather, Treelog& msg)
{
// Old weather.
has_reference_evapotranspiration
= weather.has_reference_evapotranspiration ();
has_vapor_pressure = weather.has_vapor_pressure ();
has_wind = weather.has_wind ();
has_min_max_temperature = weather.has_min_max_temperature ();
has_diffuse_radiation = weather.has_diffuse_radiation ();
old_surface = weather.surface ();
const double timestep = weather.timestep ();
// Potential evapotranspiration model.
if (!fixed_pet) // Explicit.
{
symbol type;
if (has_reference_evapotranspiration)
type = "weather";
else if (has_vapor_pressure && has_wind)
{
if (old_surface == Weatherdata::field)
type = symbol ("PM");
else if (timestep < 4.0)
type = symbol ("FAO_PM_hourly");
else
type = symbol ("FAO_PM");
}
else if (has_min_max_temperature)
type = symbol ("Hargreaves");
else
type = symbol ("deBruin87");
msg.message ("Pet choosen: " + type);
pet.reset (Librarian::build_stock<Pet> (metalib, msg, type, "pet"));
}
// Diffuse radiation model.
if (!fixed_difrad) // Explicit.
{
symbol type;
if (has_diffuse_radiation)
type = symbol ("weather");
else
type = symbol ("DPF");
msg.debug ("Difrad choosen: " + type);
difrad.reset (Librarian::build_stock<Difrad> (metalib, msg, type,
"difrad"));
}
}
bool
BioclimateStandard::check (const Weather& weather, Treelog& msg) const
{
TREELOG_MODEL (msg);
bool ok = true;
if (!pet->check (weather, msg))
ok = false;
if (!svat->check (weather, msg))
ok = false;
if (weather.surface () != Weatherdata::field
&& weather.screen_height () < reference_height)
{
std::ostringstream tmp;
tmp << "Weather ScreenHeight is " << weather.screen_height ()
<< " m, should be > " << reference_height << " m";
msg.error (tmp.str ());
ok = false;
}
return ok;
}
BioclimateStandard::BioclimateStandard (const BlockModel& al)
: Bioclimate (al),
metalib (al.metalib ()),
No (al.integer ("NoOfIntervals")),
LAI_ (0.0),
cover_ (0.0),
sun_LAI_fraction_total_ (0.0),
Height (No + 1),
total_PAR_ (No + 1),
sun_PAR_ (No + 1),
total_NIR_ (No + 1),
sun_NIR_ (No + 1),
sun_LAI_fraction_ (No),
shared_light_fraction_ (1.0),
net_radiation (Librarian::build_item<NetRadiation> (al, "net_radiation")),
pet (al.check ("pet")
? Librarian::build_item<Pet> (al, "pet")
: NULL),
total_ep_ (0.0),
total_ea_ (0.0),
direct_rain_ (0.0),
irrigation_overhead (0.0),
irrigation_overhead_temperature (0.0),
irrigation_surface (0.0),
irrigation_surface_temperature (0.0),
irrigation_subsoil (0.0),
irrigation_subsoil_permanent (al.number ("irrigation_subsoil_permanent")),
tillage_water (0.0),
snow (al.submodel ("Snow")),
snow_ep (0.0),
snow_ea_ (0.0),
snow_water_in (0.0),
snow_water_in_temperature (0.0),
snow_water_out (0.0),
snow_water_out_temperature (0.0),
canopy_ep (0.0),
canopy_ea_ (0.0),
canopy_water_capacity (0.0),
canopy_water_storage (al.number ("canopy_water_storage")),
canopy_water_temperature (0.0),
canopy_water_in (0.0),
canopy_water_out (0.0),
canopy_water_bypass (0.0),
canopy_water_below (0.0),
litter_ep (0.0),
litter_ea (0.0),
litter_water_capacity (0.0),
litter_water_storage (al.number ("litter_water_storage")),
litter_water_temperature (0.0),
litter_water_in (0.0),
litter_water_out (0.0),
pond_ep (0.0),
pond_ea_ (0.0),
soil_ep (0.0),
soil_ea_ (0.0),
max_svat_iterations (al.integer ("max_svat_iterations")),
max_svat_absolute_difference (al.number ("max_svat_absolute_difference")),
maxTdiff (al.number ("maxTdiff")),
maxEdiff (al.number ("maxEdiff")),
svat (Librarian::build_item<SVAT> (al, "svat")),
svat_fail (0),
svat_total (0),
svat_failed (true),
crop_ep_ (0.0),
crop_ea_soil (0.0),
crop_ea_svat (0.0),
crop_ea_ (0.0),
production_stress (-1.0),
albedo (NAN),
raddist (Librarian::build_item<Raddist> (al, "raddist")),
min_sin_beta_ (std::sin (al.number ("min_sun_angle"))),
difrad (al.check ("difrad")
? Librarian::build_item<Difrad> (al, "difrad")
: NULL),
difrad0 (0.0),
rad_abs_soil_ (0.0),
rad_abs_sun_canopy_ (0.0),
rad_abs_shadow_canopy_ (0.0),
sin_beta_ (0.0),
// BUG: These should really be part of the state, for checkpoints.
daily_air_temperature_ (0.0),
daily_precipitation_ (0.0),
day_length_ (0.0),
daily_global_radiation_ (0.0),
global_radiation_ (0.0),
atmospheric_CO2_ (-42.42e42),
atmospheric_O2_ (-42.42e42),
air_pressure_ (-42.42e42),
// Deposition.
deposition (Librarian::build_item<Deposition> (al, "deposition")),
// For initialization and weather data shifts.
fixed_pet (pet.get ()),
fixed_difrad (difrad.get ()),
has_reference_evapotranspiration (false),
has_vapor_pressure (false),
has_wind (false),
has_min_max_temperature (false),
has_diffuse_radiation (false),
old_surface (Weatherdata::field)
{ }
void
BioclimateStandard::summarize (Treelog& msg) const
{
TREELOG_MODEL (msg);
if (svat_fail > 0)
{
TREELOG_MODEL (msg);
daisy_assert (svat_total > 0);
std::ostringstream tmp;
tmp << "Convergence of svat loop failed " << svat_fail
<< " times out of " << svat_total << "; or "
<< (100.0 * svat_fail / (svat_total + 0.0)) << "%";
msg.warning (tmp.str ());
msg.message ("See 'daisy.log' for details");
}
svat->summarize (msg);
}
BioclimateStandard::~BioclimateStandard ()
{ }
void
BioclimateStandard::CanopyStructure (const Vegetation& vegetation)
// Calculate values for the total crop canopy.
{
shared_light_fraction_ = vegetation.shared_light_fraction ();
// Update vegetation state.
const PLF& HvsLAI = vegetation.HvsLAI ();
LAI_ = vegetation.LAI ();
cover_ = vegetation.cover ();
// Reset PAR intervals.
std::fill (Height.begin (), Height.end (), 0.0);
if (LAI_ > 0.0)
{
// Count height of each interval. Interval 0 is the top of the crop
// and interval "No" is group zero (no bomb intended).
double dLAI = LAI_ / No;
for (int i = 0; i <= No; i++)
Height[i] = HvsLAI ((No - i) * dLAI);
daisy_assert (iszero (Height[No]));
// daisy_assert (approximate (Height[0], MxH));
Height[0] = vegetation.height ();
}
}
double
BioclimateStandard::find_albedo (const Vegetation& crops, const Litter& litter,
const Surface& surface,
const Geometry& geo,
const Soil& soil, const SoilWater& soil_water)
{
const double surface_albedo = surface.albedo (geo, soil, soil_water);
const double litter_albedo = litter.albedo ();
const double litter_cover = litter.cover ();
// Find albedo below crops.
const double below_albedo = (litter_albedo < 0.0)
? surface_albedo
: litter_albedo * litter_cover + surface_albedo * (1.0 - litter_cover);
const double crop_cover = crops.cover ();
return crops.albedo () * crop_cover
+ below_albedo * (1.0 - crop_cover);
}
void
BioclimateStandard::RadiationDistribution (const Vegetation& vegetation, const double pF,
const double sin_beta_, Treelog& msg)
{
TREELOG_MODEL (msg);
raddist->tick(sun_LAI_fraction_, sun_PAR_, total_PAR_, sun_NIR_, total_NIR_,
global_radiation (), difrad0, min_sin_beta_, sin_beta_,
vegetation, pF, msg);
//Absorbed PAR in the canopy and the soil:
incoming_PAR_radiation = total_PAR_[0]; // [W/m2]
absorbed_total_PAR_canopy = incoming_PAR_radiation - total_PAR_[No]; // [W/m2]
absorbed_total_PAR_soil = incoming_PAR_radiation - absorbed_total_PAR_canopy;//[W/m2]
absorbed_sun_PAR_canopy = sun_PAR_[0] - sun_PAR_[No]; // [W/m2]
absorbed_shadow_PAR_canopy = absorbed_total_PAR_canopy - absorbed_sun_PAR_canopy; // [W/m2]
//Absorbed NIR in the canopy and the soil:
incoming_NIR_radiation = total_NIR_[0]; // [W/m2]
absorbed_total_NIR_canopy = incoming_NIR_radiation - total_NIR_[No];
absorbed_total_NIR_soil = incoming_NIR_radiation - absorbed_total_NIR_canopy;
absorbed_sun_NIR_canopy = sun_NIR_[0] - sun_NIR_[No];
absorbed_shadow_NIR_canopy = absorbed_total_NIR_canopy - absorbed_sun_NIR_canopy;
//Absorbed Long wave radiation in the canopy and the soil:
incoming_Long_radiation = net_radiation->incoming_longwave_radiation();
sun_LAI_fraction_total_ = 0.0;
for (int i = 0; i < No; i++)
sun_LAI_fraction_total_ += sun_LAI_fraction_[i] / No;
daisy_assert (sun_LAI_fraction_total_ <= 1.0);
absorbed_total_Long_soil = incoming_Long_radiation * (1-cover ());
absorbed_total_Long_canopy = incoming_Long_radiation * (cover ());
absorbed_sun_Long_canopy = absorbed_total_Long_canopy * sun_LAI_fraction_total_;
absorbed_shadow_Long_canopy = absorbed_total_Long_canopy - absorbed_sun_Long_canopy;
// Variables for SVAT
rad_abs_soil_ = absorbed_total_Long_soil + absorbed_total_NIR_soil
+ absorbed_total_PAR_soil;
rad_abs_sun_canopy_ = absorbed_sun_Long_canopy + absorbed_sun_NIR_canopy
+ absorbed_sun_PAR_canopy;
rad_abs_shadow_canopy_ = absorbed_shadow_Long_canopy
+ absorbed_shadow_NIR_canopy + absorbed_shadow_PAR_canopy;
// Variable for log
incoming_Total_radiation = incoming_Long_radiation + incoming_PAR_radiation
+ incoming_NIR_radiation; // [W/m2]
}
void
BioclimateStandard::WaterDistribution (const Time& time, Surface& surface,
const Weather& weather,
Vegetation& vegetation,
const Litter& litter,
const Movement& movement,
const Geometry& geo,
const Soil& soil,
const SoilWater& soil_water,
const SoilHeat& soil_heat,
const double T_bottom,
const double dt, Treelog& msg)
{
TREELOG_MODEL (msg);
// Overview.
//
// First we calculate the external water sources (precipitation,
// irrigation) and sinks (potential evapotranspiration). Then we
// update the sources, sinks and storage for each of the following
// four "containers": The snow pack, water intercepted on the canopy,
// pond, and soil surface. It is assumed that each the flows for
// each container only depends on the container above, so we
// calculate the changes from the top and downwards.
//
// A final and fifth "container" is the crop transpiration.
// 1 External water sinks and sources.
// Net Radiation.
const double Cloudiness = weather.cloudiness ();
const double air_temperature = weather.air_temperature ();//[dg C]
const double VaporPressure = weather.vapor_pressure ();
const double Si = weather.global_radiation ();
albedo = find_albedo (vegetation, litter, surface, geo, soil, soil_water);
net_radiation->tick (Cloudiness, air_temperature, VaporPressure, Si, albedo,
msg);
const double Rn = net_radiation->net_radiation ();
// 1.1 Fluxify management operations.
tillage_water /= dt;
// 1.2 Weather.
const double rain = weather.rain ();
// 1.3 Evapotranspiration
const double total_input
= rain + weather.snow () + irrigation_overhead + irrigation_surface;
const double free_water = total_input * dt
+ snow.storage () + surface.ponding_average ()
+ canopy_water_storage + litter_water_storage;
daisy_assert (pet.get () != NULL);
pet->tick (time,
weather, Rn, vegetation, surface,
geo, soil, soil_heat, soil_water, msg);
if (free_water > 0.01)
total_ep_ = pet->wet ();
else
total_ep_ = pet->dry ();
const double total_ep_dry = pet->dry ();
daisy_assert (std::isfinite (total_ep_));
daisy_assert (total_ep_ >= 0.0);
total_ea_ = 0.0; // To be calculated.
// 2 Snow Pack
const double rain_temperature = std::max (air_temperature, 0.1);
snow_ep = total_ep_ - total_ea_;
daisy_assert (snow_ep >= 0.0);
snow_water_in = rain + irrigation_overhead;
daisy_assert (snow_water_in >= 0.0);
if (irrigation_overhead > 0.01)
snow_water_in_temperature
= (irrigation_overhead * irrigation_overhead_temperature
+ rain * rain_temperature) / snow_water_in;
else if (rain > 0.01)
snow_water_in_temperature = rain_temperature;
else
snow_water_in_temperature = air_temperature;
snow.tick (msg, movement, soil, soil_water, soil_heat,
weather.global_radiation (), 0.0,
snow_water_in, weather.snow (),
surface.ponding_average (),
snow_water_in_temperature, snow_ep, dt);
snow_ea_ = snow.evaporation ();
daisy_assert (snow_ea_ >= 0.0);
total_ea_ += snow_ea_;
daisy_assert (total_ea_ >= 0.0);
snow_water_out = snow.percolation ();
if (snow_water_out < 0.0)
{
double adjusted_pond = snow_water_out * dt + surface.ponding_average ();
if (adjusted_pond < 0.0)
{
if (approximate (-snow_water_out * dt, surface.ponding_average ()))
adjusted_pond = 0.0;
else
{
std::ostringstream tmp;
tmp << "snow_water_out (" << snow_water_out
<< ") * dt (" << dt
<< ") + pond (" << surface.ponding_average ()
<< ") = " << adjusted_pond << ") < 0";
msg.warning (tmp.str ());
}
}
surface.put_ponding (adjusted_pond);
snow_water_out = 0.0;
}
snow_water_out_temperature = snow.temperature ();
const double below_snow_ep = total_ep_ - snow_ea_;
// 3 Water intercepted on canopy
const double canopy_water_capacity = vegetation.interception_capacity ();
daisy_assert (canopy_water_capacity >= 0.0);
const double canopy_cover = cover ();
daisy_assert (snow_ea_ <= total_ep_);
canopy_ep = below_snow_ep * canopy_cover;
const double below_canopy_ep = below_snow_ep - canopy_ep;
const double canopy_ep_dry = total_ep_dry * canopy_cover;
const double below_canopy_ep_dry = total_ep_dry - canopy_ep_dry;
if (below_canopy_ep_dry < 0.0 || !std::isfinite (below_canopy_ep_dry))
{
std::ostringstream tmp;
tmp << "below_canopy_ep_dry = " << below_canopy_ep_dry
<< ", total_ep_dry = " << total_ep_dry
<< ", canopy_cover = " << canopy_cover;
daisy_warning (tmp.str ());
}
daisy_assert (canopy_ep >= 0.0);
if (snow_water_out < 0.0)
{
canopy_water_in = 0.0;
canopy_water_bypass = snow_water_out;
}
else
{
canopy_water_in = snow_water_out * canopy_cover;
canopy_water_bypass = snow_water_out - canopy_water_in;
}
daisy_assert (canopy_water_in >= 0.0);
if (canopy_water_in > 0.01)
canopy_water_temperature
= (canopy_water_storage * rain_temperature
+ canopy_water_in * snow_water_out_temperature * dt)
/ (canopy_water_storage + canopy_water_in * dt);
else
canopy_water_temperature = air_temperature;
canopy_ea_ = std::min (canopy_ep,
canopy_water_storage / dt + canopy_water_in);
daisy_assert (canopy_ea_ >= 0.0);
total_ea_ += canopy_ea_;
daisy_assert (total_ea_ >= 0.0);
canopy_water_storage += (canopy_water_in - canopy_ea_) * dt;
if (canopy_water_storage < 0.0)
{
daisy_assert (canopy_water_storage > -1e-8);
canopy_water_storage = 0.0;
}
if (canopy_water_storage > canopy_water_capacity + 1e-8)
{
canopy_water_out = (canopy_water_storage - canopy_water_capacity) / dt;
canopy_water_storage = canopy_water_capacity;
}
else
{
canopy_water_out = 0.0;
}
daisy_assert (canopy_water_out >= 0.0);
// 4 Water intercepted by litter
const double litter_cover = litter.cover ();
litter_ep = below_canopy_ep * litter_cover;
daisy_assert (snow_ea_ <= total_ep_);
if (litter_ep < 0.0)
{
std::ostringstream tmp;
tmp << "BUG:\nlitter_ep = " << litter_ep << "\n"
<< "total_ep = " << total_ep_ << "\n"
<< "snow_ea_ = " << snow_ea_ << "\n"
<< "cover = " << canopy_cover;
msg.error (tmp.str ());
litter_ep = 0.0;
}
canopy_water_below = canopy_water_out + canopy_water_bypass
+ irrigation_surface + tillage_water;
const double canopy_water_below_temperature
= canopy_water_below > 0.01
? ((canopy_water_bypass * snow_water_out_temperature
+ canopy_water_out * canopy_water_temperature
+ irrigation_surface * irrigation_surface_temperature
+ tillage_water * rain_temperature)
/ canopy_water_below)
: air_temperature;
litter_water_in = canopy_water_below * litter_cover;
const double litter_water_bypass = canopy_water_below - litter_water_in;
if (litter_water_in > 0.01)
litter_water_temperature
= (litter_water_storage * rain_temperature
+ litter_water_in * dt * canopy_water_below_temperature)
/ (litter_water_in * dt + litter_water_storage);
else
litter_water_temperature = air_temperature;
daisy_assert (litter_water_storage >= 0.0);
litter_ea = std::max (std::min (litter_ep,
litter_water_storage / dt + litter_water_in),
0.0);
total_ea_ += litter_ea;
daisy_assert (total_ea_ >= 0.0);
litter_water_capacity = litter.water_capacity ();
litter_water_storage += (litter_water_in - litter_ea) * dt;
if (litter_water_storage < 0.0)
{
if (litter_water_storage < -1e-8)
{
std::ostringstream tmp;
tmp << "BUG:\n"
<< "litter_water_storage = " << litter_water_storage << "\n"
<< "litter_water_capacity = " << litter_water_capacity << "\n"
<< "litter_water_in = " << litter_water_in << "\n"
<< "litter_ea = " << litter_ea << "\n"
<< "litter_ep = " << litter_ep;
msg.error (tmp.str ());
}
litter_water_out = litter_water_storage / dt;
litter_water_storage = 0.0;
}
else if (litter_water_storage > litter_water_capacity + 1e-8)
{
litter_water_out = (litter_water_storage - litter_water_capacity) / dt;
litter_water_storage = litter_water_capacity;
}
else
litter_water_out = 0.0;
// 5 Ponding
pond_ep = (litter_ep - litter_ea) * litter.vapor_flux_factor()
+ below_canopy_ep * (1.0 - litter_cover);
daisy_assert (litter_ea <= litter_ep);
if (pond_ep < 0.0)
{
std::ostringstream tmp;
tmp << "BUG:\npond_ep = " << pond_ep << "\n"
<< "litter_ep = " << litter_ep << "\n"
<< "litter_ea = " << litter_ea << "\n"
<< "litter_factor = " << litter.vapor_flux_factor ();
msg.error (tmp.str ());
pond_ep = 0.0;
}
const double pond_in = litter_water_out + litter_water_bypass;
const double pond_in_temperature = (pond_in > 0.01)
? ((litter_water_out * litter_water_temperature
+ litter_water_bypass * canopy_water_below_temperature) / pond_in)
: air_temperature;
const double soil_T
= geo.content_hood (soil_heat, &SoilHeat::T, Geometry::cell_above);
surface.tick (msg, pond_ep, below_canopy_ep_dry,
pond_in, pond_in_temperature,
geo, soil, soil_water, soil_T, dt);
pond_ea_ = surface.evap_pond (dt, msg);
#ifdef NO_NEGATIVE_POND
daisy_assert (pond_ea_ >= 0.0);
#endif
total_ea_ += pond_ea_;
#ifdef NO_NEGATIVE_POND
daisy_assert (total_ea_ >= 0.0);
#endif
// 6 Soil
soil_ep = bound (0.0, pond_ep - pond_ea_,
std::max (0.0, below_canopy_ep_dry));
soil_ea_ = surface.exfiltration (dt);
daisy_assert (soil_ea_ >= 0.0);