-
Notifications
You must be signed in to change notification settings - Fork 1
/
micro_aust2.R
1597 lines (1517 loc) · 76.4 KB
/
micro_aust2.R
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
#' Australian implementation of the microclimate model.
#'
#' An implementation of the Niche Mapper microclimate model that uses the AWAP daily weather database
#' @param loc Either a longitude and latitude (decimal degrees) or a place name to search for on Google Earth
#' @param timeinterval The number of time intervals to generate predictions for over a year (must be 12 <= x <=365)
#' @param ystart First year to run
#' @param yfinish Last year to run
#' @param REFL Soil solar reflectance, decimal \%
#' @param slope Slope in degrees
#' @param aspect Aspect in degrees (0 = north)
#' @param DEP Soil depths at which calculations are to be made (cm), must be 10 values starting from 0, and more closely spaced near the surface
#' @param soiltype Soil type: Rock = 0, sand = 1, loamy sand = 2, sandy loam = 3, loam = 4, silt loam = 5, sandy clay loam = 6, clay loam = 7, silt clay loam = 8, sandy clay = 9, silty clay = 10, clay = 11, user-defined = 12, based on Campbell and Norman 1990 Table 9.1.
#' @param minshade Minimum shade level to use (\%)
#' @param maxshade Maximum shade level to us (\%)
#' @param Usrhyt Local height (m) at which air temperature, wind speed and humidity are to be computed for organism of interest
#' @param ... Additional arguments, see Details
#' @return metout The above ground micrometeorological conditions under the minimum specified shade
#' @return shadmet The above ground micrometeorological conditions under the maximum specified shade
#' @return soil Hourly predictions of the soil temperatures under the minimum specified shade
#' @return shadsoil Hourly predictions of the soil temperatures under the maximum specified shade
#' @return soilmoist Hourly predictions of the soil moisture under the minimum specified shade
#' @return shadmoist Hourly predictions of the soil moisture under the maximum specified shade
#' @return soilpot Hourly predictions of the soil water potential under the minimum specified shade
#' @return shadpot Hourly predictions of the soil water potential under the maximum specified shade
#' @return humid Hourly predictions of the soil humidity under the minimum specified shade
#' @return shadhumid Hourly predictions of the soil humidity under the maximum specified shade
#' @usage micro_aust(loc = "Melbourne, Australia", timeinterval = 365, ystart = 1990, yfinish = 1990, soiltype = 4,
#' REFL = 0.15, slope = 0, aspect = 0, DEP = c(0., 2.5, 5., 10., 15, 20, 30, 50, 100, 200), minshade = 0, maxshade = 90,
#' Usrhyt = 0.01, ...)
#' @export
#' @details
#'
#' \strong{ Parameters controling how the model runs:}
#'
#' \code{runshade}{ = 1, Run the microclimate model twice, once for each shade level (1) or just once for the minimum shade (0)?}\cr\cr
#' \code{rungads}{ = 1, Use the Global Aerosol Database? 1=yes, 0=no}\cr\cr
#' \code{write_input}{ = 0, Write csv files of final input to folder 'csv input' in working directory? 1=yes, 0=no}\cr\cr
#' \code{writecsv}{ = 0, Make Fortran code write output as csv files? 1=yes, 0=no}\cr\cr
#' \code{manualshade}{ = 1, Use CSIRO Soil and Landscape Grid of Australia? 1=yes, 0=no}\cr\cr
#' \code{soildata}{ = 1, Use CSIRO Soil and Landscape Grid of Australia? 1=yes, 0=no}\cr\cr
#' \code{terrain}{ = 0, Use 250m resolution terrain data? 1=yes, 0=no}\cr\cr
#' \code{dailywind}{ = 1, Make Fortran code write output as csv files? 1=yes, 0=no}\cr\cr
#' \code{adiab_cor}{ = 1, use adiabatic lapse rate correction? 1=yes, 0=no}\cr\cr
#' \code{warm}{ = 0, uniform warming, deg C}\cr\cr
#' \code{spatial}{ = "c:/Australian Environment/", choose location of terrain data}\cr\cr
#' \code{vlsci}{ = 0, running on the VLSCI system? 1=yes, 0=no}\cr\cr
#' \code{loop}{ = 0, if doing multiple years, this shifts the starting year by the integer value}\cr\cr
#'
#' \strong{ General additional parameters:}\cr\cr
#' \code{ERR}{ = 1.5, Integrator error tolerance for soil temperature calculations}\cr\cr
#' \code{Refhyt}{ = 1.2, Reference height (m), reference height at which air temperature, wind speed and relative humidity input data are measured}\cr\cr
#' \code{RUF}{ = 0.004, Roughness height (m), e.g. smooth desert is 0.0003, closely mowed grass may be 0.001, bare tilled soil 0.002-0.006, current allowed range: 0.00001 (snow) - 0.02 m.}\cr\cr
#' \code{Z01}{ = 0, Top (1st) segment roughness height(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used)}\cr\cr
#' \code{Z02}{ = 0, 2nd segment roughness height(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used).}\cr\cr
#' \code{ZH1}{ = 0, Top of (1st) segment, height above surface(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used).}\cr\cr
#' \code{ZH2}{ = 0, 2nd segment, height above surface(m) - IF NO EXPERIMENTAL WIND PROFILE DATA SET THIS TO ZERO! (then RUF and Refhyt used).}\cr\cr
#' \code{EC}{ = 0.0167238, Eccenricity of the earth's orbit (current value 0.0167238, ranges between 0.0034 to 0.058)}\cr\cr
#' \code{SLE}{ = 0.95, Substrate longwave IR emissivity (decimal \%), typically close to 1}\cr\cr
#' \code{Thcond}{ = 2.5, Soil minerals thermal conductivity (W/mK)}\cr\cr
#' \code{Density}{ = 2560, Soil minerals density (kg/m3)}\cr\cr
#' \code{SpecHeat}{ = 870, Soil minerals specific heat (J/kg-K)}\cr\cr
#' \code{BulkDensity}{ = 1300, Soil bulk density (kg/m3)}\cr\cr
#' \code{PCTWET}{ = 0, \% of ground surface area acting as a free water surface}\cr\cr
#' \code{rainwet}{ = 1.5, mm of rainfall causing the ground to be 90% wet for the day}\cr\cr
#' \code{cap}{ = 1, organic cap present on soil surface? (cap has lower conductivity - 0.2 W/mC - and higher specific heat 1920 J/kg-K)}\cr\cr
#' \code{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)}\cr\cr
#' \code{hori}{ = rep(0,24), Horizon angles (degrees), from 0 degrees azimuth (north) clockwise in 15 degree intervals}\cr\cr
#' \code{TIMAXS}{ = c(1.0, 1.0, 0.0, 0.0), Time of Maximums for Air Wind RelHum Cloud (h), air & Wind max's relative to solar noon, humidity and cloud cover max's relative to sunrise}\cr\cr
#' \code{TIMINS}{ = c(0, 0, 1, 1), Time of Minimums for Air Wind RelHum Cloud (h), air & Wind min's relative to sunrise, humidity and cloud cover min's relative to solar noon}\cr\cr
#' \code{timezone}{ = 0, Use GNtimezone function in package geonames to correct to local time zone (excluding daylight saving correction)? 1=yes, 0=no}\cr\cr
#'
#' \strong{ Soil moisture mode parameters:}
#'
#' \code{runmoist}{ = 0, Run soil moisture model? 1=yes, 0=no 1=yes, 0=no (note that this may cause slower runs)}\cr\cr
#' \code{PE}{ = rep(1.1,19), Air entry potential (J/kg) (19 values descending through soil for specified soil nodes in parameter}
#' \code{DEP}
#' { and points half way between)}\cr\cr
#' \code{KS}{ = rep(0.0037,19), Saturated conductivity, (kg s/m3) (19 values descending through soil for specified soil nodes in parameter}
#' \code{DEP}
#' { and points half way between)}\cr\cr
#' \code{BB}{ = rep(4.5,19), Campbell's soil 'b' parameter (-) (19 values descending through soil for specified soil nodes in parameter}
#' \code{DEP}
#' { and points half way between)}\cr\cr
#' \code{BD}{ = rep(1.3,19), Soil bulk density (kg/m3) (19 values descending through soil for specified soil nodes in parameter}
#' \code{DEP}
#' { and points half way between)}\cr\cr
#' \code{Clay}{ = 20, Clay content for matric potential calculations (\%)}\cr\cr
#' \code{SatWater}{ = rep(0.26,10), # volumetric water content at saturation (0.1 bar matric potential) (m3/m3)}\cr\cr
#' \code{maxpool}{ = 10000, Max depth for water pooling on the surface (mm), to account for runoff}\cr\cr
#' \code{rainmult}{ = 1, Rain multiplier for surface soil moisture (-), used to induce runon}\cr\cr
#' \code{evenrain}{ = 0, Spread daily rainfall evenly across 24hrs (1) or one event at midnight (0)}\cr\cr
#' \code{SoilMoist_Init}{ = c(0.1,0.12,0.15,0.2,0.25,0.3,0.3,0.3,0.3,0.3), initial soil water content at each soil node, m3/m3}\cr\cr
#' \code{L}{ = c(0,0,8.2,8.0,7.8,7.4,7.1,6.4,5.8,4.8,4.0,1.8,0.9,0.6,0.8,0.4,0.4,0,0)*10000, root density (m/m3), (19 values descending through soil for specified soil nodes in parameter}
#' \code{DEP}
#' { and points half way between)}\cr\cr
#' \code{LAI}{ = 0.1, leaf area index, used to partition traspiration/evaporation from PET}\cr\cr
#'
#' \strong{ Snow mode parameters:}
#'
#' \code{snowmodel}{ = 0, run the snow model 1=yes, 0=no (note that this may cause slower runs)}\cr\cr
#' \code{snowtemp}{ = 1.5, Temperature (deg C) at which precipitation falls as snow}\cr\cr
#' \code{snowdens}{ = 0.375, snow density (mg/m3), overridden by }
#' \code{densfun}\cr\cr
#' \code{densfun}{ = c(0,0), slope and intercept of linear model of snow density as a function of day of year - if it is c(0,0) then fixed density used}\cr\cr
#' \code{snowmelt}{ = 0.9, proportion of calculated snowmelt that doesn't refreeze}\cr\cr
#' \code{undercatch}{ = 1, undercatch multipier for converting rainfall to snow}\cr\cr
#' \code{rainmelt}{ = 0.0125, paramter in equation that melts snow with rainfall as a function of air temp}\cr\cr
#' \code{rainfrac}{ = 0.5, fraction of rain that falls on the first day of the month (decimal \% with 0 meaning rain falls evenly) - this parameter allows something other than an even intensity of rainfall when interpolating the montly rainfall data)}\cr\cr
#'
#' \strong{ Intertidal mode parameters:}
#'
#' \code{shore}{ Include tide effects? If 1, the matrix}
#' \code{tides}
#' { is used to specify tide presence, sea water temperature and presence of wavesplash}\cr\cr
#' \code{tides}{ = matrix(data = 0., nrow = 24*timeinterval*nyears, ncol = 3), matrix for each how of the simulation of 1. tide state (0=out, 1=in), 2. Water temperature (deg C) and 3. Wave splash (0=yes, 1=no)}\cr\cr
#'
#' \strong{Outputs:}
#' metout/shadmet variables:
#' \itemize{
#' \item 1 JULDAY - day of year
#' \item 2 TIME - time of day (mins)
#' \item 3 TALOC - air temperature (deg C) at local height (specified by 'Usrhyt' variable)
#' \item 4 TAREF - air temperature (deg C) at reference height (1.2m)
#' \item 5 RHLOC - relative humidity (\%) at local height (specified by 'Usrhyt' variable)
#' \item 6 RH - relative humidity (\%) at reference height (1.2m)
#' \item 7 VLOC - wind speed (m/s) at local height (specified by 'Usrhyt' variable)
#' \item 8 VREF - wind speed (m/s) at reference height (1.2m)
#' \item 9 SNOWMELT - snowmelt (mm)
#' \item 10 POOLDEP - water pooling on surface (mm)
#' \item 11 PCTWET - soil surface wetness (\%)
#' \item 12 ZEN - zenith angle of sun (degrees - 90 = below the horizon)
#' \item 13 SOLR - solar radiation (W/m2)
#' \item 14 TSKYC - sky radiant temperature (deg C)
#' \item 15 DEW - dew presence (0 or 1)
#' \item 16 FROST - frost presence (0 or 1)
#' \item 17 SNOWFALL - snow predicted to have fallen (mm)
#' \item 18 SNOWDEP - predicted snow depth (cm)
#'}
#' soil and shadsoil variables:
#' \itemize{
#' \item 1 JULDAY - day of year
#' \item 2 TIME - time of day (mins)
#' \item 3-12 D0cm ... - soil temperatures at each of the 10 specified depths
#'
#' if soil moisture model is run i.e. parameter runmoist = 1\cr
#'
#' soilmoist and shadmoist variables:
#' \itemize{
#' \item 1 JULDAY - day of year
#' \item 2 TIME - time of day (mins)
#' \item 3-12 WC0cm ... - soil moisuture (m3/m3) at each of the 10 specified depths
#'}
#' soilpot and shadpot variables:
#' \itemize{
#' \item 1 JULDAY - day of year
#' \item 2 TIME - time of day (mins)
#' \item 3-12 PT0cm ... - soil water potential (J/kg = kpa = bar/100) at each of the 10 specified depths
#' }
#'
#' humid and shadhumid variables:
#' \itemize{
#' \item 1 JULDAY - day of year
#' \item 2 TIME - time of day (mins)
#' \item 3-12 RH0cm ... - soil relative humidity (decimal \%), at each of the 10 specified depths
#' }
#' }
#' @examples
#'micro<-micro_global() # run the model with default location and settings
#'
#'metout<-as.data.frame(micro$metout) # above ground microclimatic conditions, min shade
#'shadmet<-as.data.frame(micro$shadmet) # above ground microclimatic conditions, max shade
#'soil<-as.data.frame(micro$soil) # soil temperatures, minimum shade
#'shadsoil<-as.data.frame(micro$shadsoil) # soil temperatures, maximum shade
#'
#'# append dates
#'days<-rep(seq(1,12),24)
#'days<-days[order(days)]
#'dates<-days+metout$TIME/60/24-1 # dates for hourly output
#'dates2<-seq(1,12,1) # dates for daily output
#'
#'plotmetout<-cbind(dates,metout)
#'plotsoil<-cbind(dates,soil)
#'plotshadmet<-cbind(dates,shadmet)
#'plotshadsoil<-cbind(dates,shadsoil)
#'
#'minshade<-micro$minshade
#'maxshade<-micro$maxshade
#'
#'# plotting above-ground conditions in minimum shade
#'with(plotmetout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Air Temperature (deg C)"
#', type = "l",main=paste("air temperature, ",minshade,"% shade",sep=""))})
#'with(plotmetout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Air Temperature (deg C)"
#', type = "l",lty=2,col='blue')})
#'with(plotmetout,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (%)"
#', type = "l",ylim=c(0,100),main=paste("humidity, ",minshade,"% shade",sep=""))})
#'with(plotmetout,{points(RH ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (%)"
#', type = "l",col='blue',lty=2,ylim=c(0,100))})
#'with(plotmetout,{plot(TSKYC ~ dates,xlab = "Date and Time", ylab = "Sky Temperature (deg C)"
#', type = "l",main=paste("sky temperature, ",minshade,"% shade",sep=""))})
#'with(plotmetout,{plot(VREF ~ dates,xlab = "Date and Time", ylab = "Wind Speed (m/s)"
#', type = "l",main="wind speed")})
#'with(plotmetout,{points(VLOC ~ dates,xlab = "Date and Time", ylab = "Wind Speed (m/s)"
#', type = "l",lty=2,col='blue')})
#'with(plotmetout,{plot(ZEN ~ dates,xlab = "Date and Time", ylab = "Zenith Angle of Sun (deg)"
#', type = "l",main="solar angle, sun")})
#'with(plotmetout,{plot(SOLR ~ dates,xlab = "Date and Time", ylab = "Solar Radiation (W/m2)"
#', type = "l",main="solar radiation")})
#'
#'# plotting soil temperature for minimum shade
#'for(i in 1:10){
#' if(i==1){
#' plot(plotsoil[,i+3]~plotsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature (deg C)"
#' ,col=i,type = "l",main=paste("soil temperature ",minshade,"% shade",sep=""))
#' }else{
#' points(plotsoil[,i+3]~plotsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature
#' (deg C)",col=i,type = "l")
#' }
#'}
#'
#'# plotting above-ground conditions in maximum shade
#'with(plotshadmet,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Air Temperature (deg C)"
#', type = "l",main="air temperature, sun")})
#'with(plotshadmet,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Air Temperature (deg C)"
#', type = "l",lty=2,col='blue')})
#'with(plotshadmet,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (%)"
#', type = "l",ylim=c(0,100),main="humidity, shade")})
#'with(plotshadmet,{points(RH ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (%)"
#', type = "l",col='blue',lty=2,ylim=c(0,100))})
#'with(plotshadmet,{plot(TSKYC ~ dates,xlab = "Date and Time", ylab = "Sky Temperature (deg C)",
#' type = "l",main="sky temperature, shade")})
#'
#'# plotting soil temperature for maximum shade
#'for(i in 1:10){
#' if(i==1){
#' plot(plotshadsoil[,i+3]~plotshadsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature
#' (deg C)",col=i,type = "l",main=paste("soil temperature ",maxshade,"% shade",sep=""))
#' }else{
#' points(plotshadsoil[,i+3]~plotshadsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature
#' (deg C)",col=i,type = "l")
#' }
#'}
micro_aust2 <- function(loc="Nyrripi, Northern Territory",timeinterval=365,ystart=1990,yfinish=1990,
nyears=1,soiltype=4,REFL=0.15,slope=0,aspect=0,
DEP=c(0., 2.5, 5., 10., 15, 20, 30, 50, 100, 200),
minshade=0,maxshade=90,Refhyt=1.2,Usrhyt=0.01,Z01=0,Z02=0,ZH1=0,ZH2=0,
runshade=1,rungads=1,write_input=0,writecsv=0,manualshade=1,
soildata=1,terrain=0,dailywind=1,adiab_cor=1,warm=0,spatial="c:/Australian Environment/",vlsci=0,
ERR=1.5,RUF=0.004,EC=0.0167238,SLE=0.95,Thcond=2.5,Density=2560,SpecHeat=870,BulkDensity=1300,
PCTWET=0,rainwet=1.5,cap=1,CMH2O=1.,hori=rep(0,24),
TIMAXS=c(1.0, 1.0, 0.0, 0.0),TIMINS=c(0, 0, 1, 1),timezone=0,
runmoist=1,PE=rep(1.1,19),KS=rep(0.0037,19),BB=rep(4.5,19),BD=rep(1.3,19),Clay=20,
SatWater=rep(0.26,10),maxpool=10000,rainmult=1,evenrain=0,
SoilMoist_Init=c(0.1,0.12,0.15,0.3,0.4,0.4,0.4,0.4,0.4,0.4),
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,
LAI=0.1,
snowmodel=0,snowtemp=1.5,snowdens=0.375,densfun=c(0,0),snowmelt=0.9,undercatch=1,rainmelt=0.0125,
rainfrac=0.5,
shore=0,tides=matrix(data = 0., nrow = 24*timeinterval*nyears, ncol = 3),loop=0, scenario="",barcoo="",quadrangle=1) {
#
# loc="Nyrripi, Northern Territory"
# timeinterval=365
# ystart=1990
# yfinish=1990
# soiltype=4
# REFL=0.15
# slope=0
# aspect=0
# DEP=c(0., 2.5, 5., 10., 15., 20., 30., 50., 100., 200.)
# minshade=0
# maxshade=90
# Usrhyt=.01
# runshade=1
# rungads=1
# write_input=0
# writecsv=0
# manualshade=1
# soildata=1
# terrain=0
# dailywind=1
# adiab_cor=1
# warm=0
# spatial="c:/Australian Environment/"
# vlsci=0
# loop=0
# ERR=1.5
# RUF=0.004
# EC=0.0167238
# SLE=0.95
# Thcond=2.5
# Density=2560
# SpecHeat=870
# BulkDensity=1300
# PCTWET=0
# rainwet=1.5
# cap=1
# CMH2O=1
# hori=rep(0,24)
# TIMAXS=c(1.0, 1.0, 0.0, 0.0)
# TIMINS=c(0, 0, 1, 1)
# timezone=0
# runmoist=1
# PE=rep(1.1,19)
# KS=rep(0.0037,19)
# BB=rep(4.5,19)
# BD=rep(1.3,19)
# Clay=20
# SatWater=rep(0.26,10)
# maxpool=10000
# rainmult=1
# evenrain=0
# SoilMoist_Init=c(0.1,0.12,0.15,0.2,0.25,0.3,0.3,0.3,0.3,0.3)
# 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
# LAI=0.1
# snowmodel=0
# snowtemp=1.5
# snowdens=0.375
# densfun=c(0,0)
# snowmelt=0.9
# undercatch=1
# rainmelt=0.0125
# rainfrac=0.5
# shore=0
# tides=matrix(data = 0., nrow = 24*timeinterval*nyears, ncol = 3)
if(vlsci==0){
library(RODBC)
}
errors<-0
# error trapping - originally inside the Fortran code, but now checking before executing Fortran
if(DEP[2]-DEP[1]>3 | DEP[3]-DEP[2]>3){
cat("warning, nodes might be too far apart near the surface, try a different spacing if the program is crashing \n")
}
if(timeinterval<12 | timeinterval > 365){
cat("ERROR: the variable 'timeinterval' is out of bounds.
Please enter a correct value (12 - 365).", '\n')
errors<-1
}
if(is.numeric(loc[1])){
if(loc[1]>180 | loc[2] > 90){
cat("ERROR: Latitude or longitude (longlat) is out of bounds.
Please enter a correct value.", '\n')
errors<-1
}
}
if(timezone%in%c(0,1)==FALSE){
cat("ERROR: the variable 'timezone' be either 0 or 1.
Please correct.", '\n')
errors<-1
}
if(rungads%in%c(0,1)==FALSE){
cat("ERROR: the variable 'rungads' be either 0 or 1.
Please correct.", '\n')
errors<-1
}
if(Clay<0){
cat("ERROR: Clay density value (Clay) is negative.
Please input a positive value.", '\n')
errors<-1
}
if(write_input%in%c(0,1)==FALSE){
cat("ERROR: the variable 'write_input' be either 0 or 1.
Please correct.", '\n')
errors<-1
}
if(EC<0.0034 | EC > 0.058){
cat("ERROR: the eccentricity variable (EC) is out of bounds.
Please enter a correct value (0.0034 - 0.058).", '\n')
errors<-1
}
if(RUF<0.0001){
cat("ERROR: The roughness height (RUF) is too small ( < 0.0001).
Please enter a larger value.", '\n')
errors<-1
}
if(RUF>2){
cat("ERROR: The roughness height (RUF) is too large ( > 2).
Please enter a smaller value.", '\n')
errors<-1
}
if(DEP[1]!=0){
cat("ERROR: First soil node (DEP[1]) must = 0 cm.
Please correct", '\n')
errors<-1
}
if(length(DEP)!=10){
cat("ERROR: You must enter 10 different soil depths.", '\n')
errors<-1
}
for(i in 1:9){
if(DEP[i+1]<=DEP[i]){
cat("ERROR: Soil depth (DEP array) is not in ascending size", '\n')
errors<-1
}
}
if(DEP[10]>500){
cat("ERROR: Deepest soil depth (DEP array) is too large (<=500 cm)", '\n')
errors<-1
}
if(Thcond<0){
cat("ERROR: Thermal variable conductivity (THCOND) is negative.
Please input a positive value.", '\n')
errors<-1
}
if(Density<0){
cat("ERROR: Density variable (Density) is negative.
Please input a positive value.", '\n')
errors<-1
}
if(SpecHeat<0){
cat("ERROR: Specific heat variable (SpecHeat) is negative.
Please input a positive value.", '\n')
errors<-1
}
if(BulkDensity<0){
cat("ERROR: Bulk density value (BulkDensity) is negative.
Please input a positive value.", '\n')
errors<-1
}
if(REFL<0 | REFL>1){
cat("ERROR: Soil reflectivity value (REFL) is out of bounds.
Please input a value between 0 and 1.", '\n')
errors<-1
}
if(slope<0 | slope>90){
cat("ERROR: Slope value (slope) is out of bounds.
Please input a value between 0 and 90.", '\n')
errors<-1
}
if(aspect<0 | aspect>365){
cat("ERROR: Aspect value (aspect) is out of bounds.
Please input a value between 0 and 365.", '\n')
errors<-1
}
if(max(hori)>90 | min(hori)<0){
cat("ERROR: At least one of your horizon angles (hori) is out of bounds.
Please input a value between 0 and 90", '\n')
errors<-1
}
if(length(hori)!=24){
cat("ERROR: You must enter 24 horizon angle values.", '\n')
errors<-1
}
if(SLE<0.5 | SLE > 1){
cat("ERROR: Emissivity (SLE) is out of bounds.
Please enter a correct value (0.05 - 1.00).", '\n')
errors<-1
}
if(ERR<0){
cat("ERROR: Error bound (ERR) is too small.
Please enter a correct value (> 0.00).", '\n')
errors<-1
}
if(Usrhyt<RUF){
cat("ERROR: Reference height (Usrhyt) smaller than roughness height (RUF).
Please use a larger height above the surface.", '\n')
errors<-1
}
if(Usrhyt<0.005 | Usrhyt>Refhyt){
cat("ERROR: Local height (Usrhyt) is out of bounds.
Please enter a correct value (0.005 - Refhyt).", '\n')
errors<-1
}
if(CMH2O<0.5 | CMH2O>2){
cat("ERROR: Preciptable water in air column (CMH2O) is out of bounds.
Please enter a correct value (0.1 - 2cm).", '\n')
errors<-1
}
if(max(TIMAXS)>24 | min(TIMAXS)<0){
cat("ERROR: At least one of your times of weather maxima (TIMAXS) is out of bounds.
Please input a value between 0 and 24", '\n')
errors<-1
}
if(max(TIMINS)>24 | min(TIMINS)<0){
cat("ERROR: At least one of your times of weather minima (TIMINS) is out of bounds.
Please input a value between 0 and 24", '\n')
errors<-1
}
if(minshade>maxshade | minshade==maxshade){
cat("ERROR: Your value for minimum shade (minshade) is greater than or equal to the maximum shade (maxshade).
Please correct this.", '\n')
errors<-1
}
if(minshade>100 | minshade<0){
cat("ERROR: Your value for minimum shade (minshade) is out of bounds.
Please input a value between 0 and 100.", '\n')
errors<-1
}
if(maxshade>100 | maxshade<0){
cat("ERROR: Your value for maximum shade (maxshade) is out of bounds.
Please input a value between 0 and 100.", '\n')
errors<-1
}
if(soiltype<0 | soiltype>11){
cat("ERROR: the soil type must range between 1 and 11.
Please correct.", '\n')
errors<-1
}
# end error trapping
if(errors==0){ # continue
################## time related variables #################################
nyears<-yfinish-ystart+1
if(vlsci==1){
ndays<-365*nyears
}
juldays12<-c(15.,46.,74.,105.,135.,166.,196.,227.,258.,288.,319.,349.) # middle day of each month
juldaysn<-juldays12 # variable of juldays for when doing multiple years
if(nyears>1 & timeinterval==365){ # create sequence of days for splining across multiple years
for(i in 1:(nyears-1)){
juldaysn<-c(juldaysn,(juldays12+365*i))
}
}
if(timeinterval<365){
microdaily<-0 # run microclimate model as normal, where each day is iterated 3 times starting with the initial condition of uniform soil temp at mean monthly temperature
}else{
microdaily<-1 # run microclimate model where one iteration of each day occurs and last day gives initial conditions for present day with an initial 3 day burn in
}
# now check if doing something other than middle day of each month, and create appropriate vector of julian days
daystart<-as.integer(ceiling(365/timeinterval/2))
if(timeinterval!=12){
juldays<-seq(daystart,365,as.integer(floor(365/timeinterval)))
}else{
juldays<-juldaysn
}
julnum <- timeinterval*nyears # total days to do
julday <- subset(juldays, juldays!=0) # final vector of julian days
julday<-rep(julday,nyears)
idayst <- 1 # start day
ida<-timeinterval*nyears # end day
dates<-Sys.time()-60*60*24
curyear<-as.numeric(format(dates,"%Y"))
################## location and terrain #################################
if(vlsci==0){
f1 <- paste(spatial,"ausclim_rowids.nc",sep="");
f2 <- paste(spatial,"ausdem_shift1.tif",sep="");
f3 <- paste(spatial,"agg_9secdem.nc",sep="");
f4 <- paste(spatial,"Aust9secDEM.tif",sep="");
}else{
f1 <- "/vlsci/VR0212/shared/Spatial_Data/ausclim_rowids.nc";
f2 <- "/vlsci/VR0212/shared/Spatial_Data/ausdem_shift1.tif";
f3 <- "/vlsci/VR0212/shared/Spatial_Data/agg_9secdem.nc";
f4 <- "/vlsci/VR0212/shared/Spatial_Data/Aust9secDEM.grd";
}
if(is.numeric(loc)==FALSE){ # use geocode to get location from site name via googlemaps
if (!requireNamespace("dismo", quietly = TRUE)) {
stop("dismo needed for the place name geocode function to work. Please install it.",
call. = FALSE)
}
if (!requireNamespace("XML", quietly = TRUE)) {
stop("XML needed for the place name geocode function to work. Please install it.",
call. = FALSE)
}
longlat <- dismo::geocode(loc)[3:4] # assumes first geocode match is correct
if(nrow(longlat>1)){longlat<-longlat[1,]}
x <- t(as.matrix(as.numeric(c(longlat[1,1],longlat[1,2]))))
}else{
longlat <- loc
x <- t(as.matrix(as.numeric(c(loc[1],loc[2]))))
}
# get the local timezone reference longitude
if(timezone==1){ # this now requires registration
if(!require(geonames, quietly = TRUE)){
stop('package "geonames" is required to do a specific time zone (timezone=1). Please install it.')
}
ALREF<-(geonames::GNtimezone(longlat[2],longlat[1])[4])*-15
}else{ # just use local solar noon
ALREF <- abs(trunc(x[1]))
}
HEMIS <- ifelse(x[2]<0,2.,1.) # 1 is northern hemisphere
# break decimal degree lat/lon into deg and min
ALAT <- abs(trunc(x[2]))
AMINUT <- (abs(x[2])-ALAT)*60
ALONG <- abs(trunc(x[1]))
ALMINT <- (abs(x[1])-ALONG)*60
azmuth<-aspect
if(soildata==0){
soilprop<-cbind(0,0)
# creating the shade array
MAXSHADES <- rep(0,(timeinterval*nyears))+maxshade # daily max shade (%)
MINSHADES <- rep(0,(timeinterval*nyears))+minshade # daily min shade (%)
}
if(soildata==1){
cat("extracting soil data", '\n')
if(vlsci==0){
static_soil<-paste(spatial,"static_soil.nc",sep="")
emissivities<-paste(spatial,"aus_emissivities.nc",sep="")
}else{
static_soil<-'/vlsci/VR0212/shared/Spatial_Data/static_soil.nc'
emissivities<-'/vlsci/VR0212/shared/Spatial_Data/aus_emissivities.nc'
}
# read data in from netcdf file
static_soil_data<-raster::brick(static_soil)
static_soil_vars <- raster::extract(static_soil_data,x)
labels<-c('albedo','FAPAR1','FAPAR2','FAPAR3','FAPAR4','FAPAR5','FAPAR6','FAPAR7','FAPAR8','FAPAR9','FAPAR10','FAPAR11','FAPAR12','volwater_Upper','volwater_lower','thick_upper','thick_lower','code')
colnames(static_soil_vars)<-labels
emissivities_data<-raster::brick(emissivities)
SLES2 <- raster::extract(emissivities_data,x)
# read in other soil related files for working out lumped soil type and properties
# such as clay % for getting water potential
if(vlsci==0){
filename<-paste(spatial,"ppfInterpAll.txt",sep="")
ppf<-as.data.frame(read.table(file = filename, sep = ",", header=TRUE))
filename<-paste(spatial,"Lumped soil types.txt",sep="")
lumped.soil<-as.data.frame(read.table(file = filename, sep = ","))
filename<-paste(spatial,"SoilTypeLUT_725_AWAP.csv",sep="")
soiltype<-as.data.frame(read.table(file = filename, sep = ","))
}else{
filename<-'/vlsci/VR0212/shared/Spatial_Data/ppfInterpAll.txt'
ppf<-as.data.frame(read.table(file = filename, sep = ",", header=TRUE))
filename<-'/vlsci/VR0212/shared/Spatial_Data/Lumped soil types.txt'
lumped.soil<-as.data.frame(read.table(file = filename, sep = ","))
filename<-'/vlsci/VR0212/shared/Spatial_Data/SoilTypeLUT_725_AWAP.csv'
soiltype<-as.data.frame(read.table(file = filename, sep = ","))
}
soilcode<-subset(soiltype, soiltype[1]==static_soil_vars[18])
lumped<-subset(lumped.soil, V4==as.character(soilcode[1,2]))
soiltype<-lumped[1,6]
soilprop<-subset(ppf, ppf==soilcode[1,2])
}else{
SLES2 <- rep(SLE,timeinterval*nyears)
if(manualshade==0){
cat("extracting shade data", '\n')
if(vlsci==0){
static_soil<-paste(spatial,"static_soil.nc",sep="")
emissivities<-paste(spatial,"aus_emissivities.nc",sep="")
}else{
static_soil<-'/vlsci/VR0212/shared/Spatial_Data/static_soil.nc'
emissivities<-'/vlsci/VR0212/shared/Spatial_Data/aus_emissivities.nc'
}
# read data in from netcdf file
static_soil_data<-raster::brick(static_soil)
static_soil_vars <- raster::extract(static_soil_data,x)
labels<-c('albedo','FAPAR1','FAPAR2','FAPAR3','FAPAR4','FAPAR5','FAPAR6','FAPAR7','FAPAR8','FAPAR9','FAPAR10','FAPAR11','FAPAR12','volwater_Upper','volwater_lower','thick_upper','thick_lower','code')
colnames(static_soil_vars)<-labels
}
}
if(terrain==1){
cat("extracting terrain data \n")
e<-extent(x[1]-0.05,x[1]+0.05,x[2]-0.05,x[2]+0.05)
for(i in 1:24){
if(vlsci==0){
horifile<-paste(spatial,'horizon',i,'.tif',sep="")
}else{
horifile<-paste('/vlsci/VR0212/shared/Spatial_Data/','horizon',i,'.tif',sep="")
}
horiz<-raster::crop(raster::raster(horifile),e)
if(i==1){
horizons_data<-horiz
}else{
horizons_data<-raster::stack(horizons_data,horiz)
}
}
HORIZONS <- t(raster::extract(horizons_data,x))
if(vlsci==0){
elev1<-crop(raster::raster(paste(spatial,'elev.tif',sep="")),e)
slope1<-crop(raster::raster(paste(spatial,'slope.tif',sep="")),e)
aspect1<-crop(raster::raster(paste(spatial,'aspect.tif',sep="")),e)
elevslpasp<-raster::stack(elev1,slope1,aspect1)
}else{
elev1<-raster::crop(raster(paste('/vlsci/VR0212/shared/Spatial_Data/','elev.tif',sep="")),e)
slope1<-raster::crop(raster(paste('/vlsci/VR0212/shared/Spatial_Data/','slope.tif',sep="")),e)
aspect1<-raster::crop(raster(paste('/vlsci/VR0212/shared/Spatial_Data/','aspect.tif',sep="")),e)
elevslpasp<-raster::stack(elev1,slope1,aspect1)
}
ELEVSLPASP <- raster::extract(elevslpasp,x)
ELEVSLPASP<-as.matrix((ifelse(is.na(ELEVSLPASP),0,ELEVSLPASP)))
ALTITUDES <- ELEVSLPASP[,1]
SLOPES <- ELEVSLPASP[,2]
AZMUTHS <- ELEVSLPASP[,3]
# the horizons have been arranged so that they go from 0 degrees azimuth (north) clockwise - r.horizon starts
# in the east and goes counter clockwise!
HORIZONS <- (ifelse(is.na(HORIZONS),0,HORIZONS))/10 # get rid of na and get back to floating point
HORIZONS <- data.frame(HORIZONS)
VIEWF_all <- 1-rowSums(sin(t(HORIZONS)*pi/180))/length(t(HORIZONS)) # convert horizon angles to radians and calc view factor(s)
r1 <- raster::raster(f1)
r2 <- raster::raster(f2)
r3 <- raster::raster(f3)
dbrow <- raster::extract(r1, x)
AUSDEM <- raster::extract(r2, x)
AGG <- raster::extract(r3, x)
}else{
r1 <- raster::raster(f1)
r2 <- raster::raster(f2)
r3 <- raster::raster(f3)
r4 <- raster::raster(f4)
dbrow <- raster::extract(r1, x)
AUSDEM <- raster::extract(r2, x)
AGG <- raster::extract(r3, x)
ALTITUDES <- raster::extract(r4, x)
#ALTITUDES <- AUSDEM
#cat("using 0.05 res DEM!")
HORIZONS <- hori
HORIZONS <- data.frame(HORIZONS)
VIEWF_all <- rep(1,length(x[,1]))
SLOPES<-rep(slope,length(x[,1]))
AZMUTHS<-rep(aspect,length(x[,1]))
}
hori<-HORIZONS
row.names(hori)<-NULL
hori<-as.numeric(as.matrix(hori))
if(soildata==1){
VIEWF<-VIEWF_all
SLES<-SLES2
}else{
VIEWF<-VIEWF_all
}
# setting up for temperature correction using lapse rate given difference between 9sec DEM value and 0.05 deg value
if(AUSDEM==-9999 | is.na(AUSDEM)=='TRUE'){
delta_elev = AGG - ALTITUDES
}else{
delta_elev = AUSDEM - ALTITUDES
}
adiab_corr = delta_elev * 0.0058 # Adiabatic temperature correction for elevation (C), mean for Australian Alps
adiab_corr_max = delta_elev * 0.0077 # Adiabatic temperature correction for elevation (C), mean for Australian Alps
adiab_corr_min = delta_elev * 0.0039 # Adiabatic temperature correction for elevation (C), mean for Australian Alps
cat("extracting climate data", '\n')
# connect to server
if(vlsci==0){
channel2 <- RODBC::odbcConnect("ausclim_predecol")
channel <- RODBC::odbcConnect("AWAPDaily")
# preliminary test for incomplete year, if simulation includes the present year
yearlist<-seq(ystart,(ystart+(nyears-1)),1)
for(j in 1:nyears){ # start loop through years
yeartodo<-yearlist[j]
lat1<-x[2]-0.024
lat2<-x[2]+0.025
lon1<-x[1]-0.024
lon2<-x[1]+0.025
query<-paste("SELECT a.latitude, a.longitude, b.*
FROM [AWAPDaily].[dbo].[latlon] as a
, [AWAPDaily].[dbo].[",yeartodo,"] as b
where (a.id = b.id) and (a.latitude between ",lat1," and ",lat2,") and (a.longitude between ",lon1," and ",lon2,")
order by b.day",sep="")
output<- sqlQuery(channel,query)
if(yearlist[j]<1971){
output$vpr<-output$tmin/output$tmin-1
}
if(yearlist[j]>1989){
output$sol<-as.numeric(as.character(output$sol))
}else{
output$sol<-output$tmin/output$tmin-1
}
if(nrow(output)>365){
# fix leap years
output<-output[-60,]
}
if(j==1){
results<-output
}else{
results<-rbind(results,output)
}
}
nyears2<-nrow(results)/365
ndays<-nrow(results)
juldaysn2<-juldaysn[juldaysn<=ndays]
juldaysn2<-juldaysn[1:round(nyears2*12)]
} #end vlsci check
dim<-ndays
juldays<-seq(daystart,dim,1)
julday<-rep(seq(1,365),nyears)
juday<-julday[1:dim]
ida<-ndays
idayst <- 1 # start month
# end preliminary test for incomplete year, if simulation includes the present year
if((soildata==1 & nrow(soilprop)>0)|soildata==0){
if(soildata==1){
# get static soil data into arrays
REFL <- static_soil_vars[,1] # albedo/reflectances
maxshades <- static_soil_vars[,2:13] # assuming FAPAR represents shade
shademax<-maxshades
}else{
if(manualshade==0){
maxshades <- static_soil_vars[,2:13] # assuming FAPAR represents shade
}
shademax<-maxshade
}
if(is.na(dbrow)!=TRUE & is.na(ALTITUDES)!=TRUE){
if(rungads==1){
####### get solar attenuation due to aerosols with program GADS #####################
relhum<-1.
optdep.summer<-as.data.frame(rungads(longlat[2],longlat[1],relhum,0))
optdep.winter<-as.data.frame(rungads(longlat[2],longlat[1],relhum,1))
optdep<-cbind(optdep.winter[,1],rowMeans(cbind(optdep.summer[,2],optdep.winter[,2])))
optdep<-as.data.frame(optdep)
colnames(optdep)<-c("LAMBDA","OPTDEPTH")
a<-lm(OPTDEPTH~poly(LAMBDA, 6, raw=TRUE),data=optdep)
LAMBDA<-c(290,295,300,305,310,315,320,330,340,350,360,370,380,390,400,420,440,460,480,500,520,540,560,580,600,620,640,660,680,700,720,740,760,780,800,820,840,860,880,900,920,940,960,980,1000,1020,1080,1100,1120,1140,1160,1180,1200,1220,1240,1260,1280,1300,1320,1380,1400,1420,1440,1460,1480,1500,1540,1580,1600,1620,1640,1660,1700,1720,1780,1800,1860,1900,1950,2000,2020,2050,2100,2120,2150,2200,2260,2300,2320,2350,2380,2400,2420,2450,2490,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000)
TAI<-predict(a,data.frame(LAMBDA))
################ end GADS ##################################################
}else{ # use a suitable one for Australia (same as around Adelaide/Melbourne)
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)
} #end check if running gads
if(vlsci==0){
yearlist<-seq(ystart,(ystart+(nyears-1)),1)
for(j in 1:nyears){ # start loop through years
yeartodo<-yearlist[j]
lat1<-x[2]-0.024
lat2<-x[2]+0.025
lon1<-x[1]-0.024
lon2<-x[1]+0.025
query<-paste("SELECT a.latitude, a.longitude, b.*
FROM [AWAPDaily].[dbo].[latlon] as a
, [AWAPDaily].[dbo].[",yeartodo,"] as b
where (a.id = b.id) and (a.latitude between ",lat1," and ",lat2,") and (a.longitude between ",lon1," and ",lon2,")
order by b.day",sep="")
output<- sqlQuery(channel,query)
if(yearlist[j]<1971){
output$vpr<-output$tmin/output$tmin-1
}
if(yearlist[j]>1989){
output$sol<-as.numeric(as.character(output$sol))
}else{
output$sol<-output$tmin/output$tmin-1
}
if(nrow(output)>365){
# fix leap years
output<-output[-60,]
}
if(j==1){
results<-output
}else{
results<-rbind(results,output)
}
}
if(dailywind==1){
channel <- RODBC::odbcConnect("dailywind")
if(min(yearlist)<1975){
# get mean of 1975-1984
for(j in 1:10){ # start loop through years
yeartodo<-1974+j
lat1<-x[2]-0.024
lat2<-x[2]+0.025
lon1<-x[1]-0.024
lon2<-x[1]+0.025
query<-paste("SELECT a.latitude, a.longitude, b.*
FROM [dailywind].[dbo].[latlon] as a
, [dailywind].[dbo].[",yeartodo,"] as b
where (a.id = b.id) and (a.latitude between ",lat1," and ",lat2,") and (a.longitude between ",lon1," and ",lon2,")
order by b.day",sep="")
output<- sqlQuery(channel,query)
if(nrow(output)>365){
# fix leap years
output<-output[-60,]
}
if(j==1){
dwindmean<-output
}else{
dwindmean<-cbind(dwindmean,output[,5])
}
}
dwindmean<-cbind(dwindmean[,1:4],rowMeans(dwindmean[,5:14]))
colnames(dwindmean)[5]<-'wind'
}
for(j in 1:nyears){ # start loop through years
yeartodo<-yearlist[j]
if(yeartodo<1975){
output<-dwindmean
}else{
lat1<-x[2]-0.024
lat2<-x[2]+0.025
lon1<-x[1]-0.024
lon2<-x[1]+0.025
query<-paste("SELECT a.latitude, a.longitude, b.*
FROM [dailywind].[dbo].[latlon] as a
, [dailywind].[dbo].[",yeartodo,"] as b
where (a.id = b.id) and (a.latitude between ",lat1," and ",lat2,") and (a.longitude between ",lon1," and ",lon2,")
order by b.day",sep="")
output<- sqlQuery(channel,query)
if(nrow(output)>365){
# fix leap years
output<-output[-60,]
}
}
if(j==1){
dwind<-output
}else{
dwind<-rbind(dwind,output)
}
}
dwind<-dwind$wind/15.875
}
}else{ #vlsci==1
load(paste(barcoo,'TMAXX.bin',sep=''))
TMAXX <- as.matrix(data[quadrangle,2:7301])
load(paste(barcoo,'TMINN.bin',sep=''))
TMINN <- as.matrix(data[quadrangle,2:7301])
load(paste(barcoo,'RHMAXX.bin',sep=''))
RHMAXX <- as.numeric(as.matrix(data[quadrangle,2:7301]))
load(paste(barcoo,'RHMINN.bin',sep=''))
RHMINN <- as.numeric(as.matrix(data[quadrangle,2:7301]))
load(paste(barcoo,'CCMAXX.bin',sep=''))
CCMAXX <- as.numeric(as.matrix(data[quadrangle,2:7301]))
load(paste(barcoo,'CCMINN.bin',sep=''))
CCMINN <- as.numeric(as.matrix(data[quadrangle,2:7301]))
load(paste(barcoo,'RAINFALL.bin',sep=''))
RAINFALL <- as.matrix(data[quadrangle,2:7301])
RHMAXX<-t(RHMAXX)
RHMINN<-t(RHMINN)
CCMAXX<-t(CCMAXX)
CCMINN<-t(CCMINN)
TMAXX<-t(as.data.frame(TMAXX[((ystart-1989)*365-364):((yfinish-1989)*365)]))
TMINN<-t(as.data.frame(TMINN[((ystart-1989)*365-364):((yfinish-1989)*365)]))
RHMAXX<-t(as.data.frame(RHMAXX[((ystart-1989)*365-364):((yfinish-1989)*365)]))
RHMINN<-t(as.data.frame(RHMINN[((ystart-1989)*365-364):((yfinish-1989)*365)]))
CCMAXX<-t(as.data.frame(CCMAXX[((ystart-1989)*365-364):((yfinish-1989)*365)]))
CCMINN<-t(as.data.frame(CCMINN[((ystart-1989)*365-364):((yfinish-1989)*365)]))
RAINFALL<-t(as.data.frame(RAINFALL[((ystart-1989)*365-364):((yfinish-1989)*365)]))
if(dailywind==1){
load(paste(barcoo,'dwind.bin',sep=''))
dwind<- as.matrix(data[quadrangle,2:7301])
dwind<-t(as.data.frame(dwind[((ystart-1989)*365-364):((yfinish-1989)*365)]))
#WNMINN<-WNMAXX
}
if(adiab_cor==1){
TMAXX<-TMAXX+adiab_corr_max
TMINN<-TMINN+adiab_corr_min
}
} #end vlsci check
if(vlsci==0){
if(adiab_cor==1){
TMAXX<-as.matrix(results$tmax+adiab_corr_max)
TMINN<-as.matrix(results$tmin+adiab_corr_min)
}else{
TMAXX<-as.matrix(results$tmax)
TMINN<-as.matrix(results$tmin)
}
RAINFALL<-results$rr
dim<-length(RAINFALL)
output_AWAPDaily<-results
}
# cloud cover
if(vlsci==0){
if(ystart>1989 & sum(results[,9],na.rm=TRUE)>0){ # solar radiation data available
query<-paste("SELECT a.*
FROM [ausclim].[dbo].[clearskysol] as a
where (a.latitude between ",lat1," and ",lat2,") and (a.longitude between ",lon1," and ",lon2,")
",sep="")
output_ausclim<- sqlQuery(channel,query)
if(nrow(output_ausclim)==0){ #no satellite coverage, get data from ausclim
clouds<-paste("select cloud1,cloud2,cloud3,cloud4,cloud5,cloud6,cloud7,cloud8,cloud9,cloud10,cloud11,cloud12 FROM cloudcover WHERE i = ",dbrow,sep="")
#CCMAXX <- dbGetQuery(con,statement=clouds)*100
if(vlsci==0){
CCMAXX<- sqlQuery(channel2,clouds)*100
}
CCMINN <- CCMAXX
CCMAXX1 <-suppressWarnings(spline(juldays12,CCMAXX,n=timeinterval,xmin=1,xmax=365,method="periodic"))
CCMAXX <- rep(CCMAXX1$y,nyears)
CCMINN <- CCMAXX
}else{
weekly_sol<-cbind(1:52,t(output_ausclim[3:54]))
daily_sol <-suppressWarnings(spline(seq(3,361,7),weekly_sol[,2],n=365,xmin=1,xmax=365,method="periodic"))
daily_sol<-as.numeric(daily_sol$y)
if(yfinish==curyear){ # if it's the current year, then make sure daily_sol goes however far into the year we are
daily_sol<-c(rep(daily_sol,nyears-1),daily_sol[1:(ndays-(nyears-1)*365)])
}else{
daily_sol<-rep(daily_sol,nyears)
}
if(is.na(output_AWAPDaily[1,9])==TRUE){
output_AWAPDaily[1,9]=mean(output_AWAPDaily[,9],na.rm=TRUE)
}
if(is.na(output_AWAPDaily[dim,9])==TRUE){
output_AWAPDaily[nrow(output_AWAPDaily),9]=mean(output_AWAPDaily[,9],na.rm=TRUE)
}
solar<-zoo::na.approx(output_AWAPDaily[,9])