This repository has been archived by the owner on Sep 14, 2018. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 12
/
Copy pathEDPhysiologyMod.F90
2419 lines (2003 loc) · 111 KB
/
EDPhysiologyMod.F90
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
module EDPhysiologyMod
#include "shr_assert.h"
! ============================================================================
! Miscellaneous physiology routines from ED.
! ============================================================================
use FatesGlobals, only : fates_log
use FatesInterfaceMod, only : hlm_days_per_year
use FatesInterfaceMod, only : hlm_model_day
use FatesInterfaceMod, only : hlm_freq_day
use FatesInterfaceMod, only : hlm_day_of_year
use FatesInterfaceMod, only : numpft
use FatesInterfaceMod, only : hlm_use_planthydro
use FatesConstantsMod, only : r8 => fates_r8
use EDPftvarcon , only : EDPftvarcon_inst
use FatesInterfaceMod, only : bc_in_type
use EDCohortDynamicsMod , only : zero_cohort
use EDCohortDynamicsMod , only : create_cohort, sort_cohorts
use FatesAllometryMod , only : tree_lai
use FatesAllometryMod , only : tree_sai
use EDTypesMod , only : numWaterMem
use EDTypesMod , only : dl_sf, dinc_ed
use EDTypesMod , only : external_recruitment
use EDTypesMod , only : ncwd
use EDTypesMod , only : nlevleaf
use EDTypesMod , only : senes
use EDTypesMod , only : maxpft
use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type
use EDTypesMod , only : dump_cohort
use shr_log_mod , only : errMsg => shr_log_errMsg
use FatesGlobals , only : fates_log
use FatesGlobals , only : endrun => fates_endrun
use EDParamsMod , only : fates_mortality_disturbance_fraction
use FatesPlantHydraulicsMod , only : AccumulateMortalityWaterStorage
use FatesConstantsMod , only : itrue,ifalse
use FatesConstantsMod , only : calloc_abs_error
use FatesAllometryMod , only : h_allom
use FatesAllometryMod , only : h2d_allom
use FatesAllometryMod , only : bagw_allom
use FatesAllometryMod , only : bsap_allom
use FatesAllometryMod , only : bleaf
use FatesAllometryMod , only : bfineroot
use FatesAllometryMod , only : bdead_allom
use FatesAllometryMod , only : bstore_allom
use FatesAllometryMod , only : bbgw_allom
use FatesAllometryMod , only : carea_allom
use FatesAllometryMod , only : CheckIntegratedAllometries
use FatesAllometryMod , only : StructureResetOfDH
use FatesIntegratorsMod, only : RKF45
use FatesIntegratorsMod, only : Euler
implicit none
private
public :: non_canopy_derivs
public :: trim_canopy
public :: phenology
private :: phenology_leafonoff
public :: PlantGrowth
public :: recruitment
private :: cwd_input
private :: cwd_out
private :: fragmentation_scaler
private :: seeds_in
private :: seed_decay
private :: seed_germination
public :: flux_into_litter_pools
logical, parameter :: DEBUG = .false. ! local debug flag
character(len=*), parameter, private :: sourcefile = &
__FILE__
integer, parameter :: i_dbh = 1 ! Array index associated with dbh
integer, parameter :: i_cleaf = 2 ! Array index associated with leaf carbon
integer, parameter :: i_cfroot = 3 ! Array index associated with fine-root carbon
integer, parameter :: i_csap = 4 ! Array index associated with sapwood carbon
integer, parameter :: i_cstore = 5 ! Array index associated with storage carbon
integer, parameter :: i_cdead = 6 ! Array index associated with structural carbon
integer, parameter :: i_crepro = 7 ! Array index associated with reproductive carbon
integer, parameter :: n_cplantpools = 7 ! Size of the carbon only integration framework
! ============================================================================
contains
! ============================================================================
subroutine non_canopy_derivs( currentSite, currentPatch, bc_in )
!
! !DESCRIPTION:
! Returns time differentials of the state vector
!
! !USES:
use EDTypesMod, only : AREA
!
! !ARGUMENTS
type(ed_site_type), intent(inout), target :: currentSite
type(ed_patch_type), intent(inout) :: currentPatch
type(bc_in_type), intent(in) :: bc_in
!
! !LOCAL VARIABLES:
integer c,p
!----------------------------------------------------------------------
currentPatch%leaf_litter_in(:) = 0.0_r8
currentPatch%root_litter_in(:) = 0.0_r8
currentPatch%dleaf_litter_dt(:) = 0.0_r8
currentPatch%droot_litter_dt(:) = 0.0_r8
currentPatch%leaf_litter_out(:) = 0.0_r8
currentPatch%root_litter_out(:) = 0.0_r8
currentPatch%cwd_AG_in(:) = 0.0_r8
currentPatch%cwd_BG_in(:) = 0.0_r8
currentPatch%cwd_AG_out(:) = 0.0_r8
currentPatch%cwd_BG_out(:) = 0.0_r8
currentPatch%seeds_in(:) = 0.0_r8
currentPatch%seed_decay(:) = 0.0_r8
currentPatch%seed_germination(:) = 0.0_r8
! update seed fluxes
call seeds_in(currentSite, currentPatch)
call seed_decay(currentSite, currentPatch)
call seed_germination(currentSite, currentPatch)
! update fragmenting pool fluxes
call cwd_input( currentSite, currentPatch)
call cwd_out( currentSite, currentPatch, bc_in)
do p = 1,numpft
currentSite%dseed_dt(p) = currentSite%dseed_dt(p) + &
(currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - &
currentPatch%seed_germination(p)) * currentPatch%area/AREA
enddo
do c = 1,ncwd
currentPatch%dcwd_AG_dt(c) = currentPatch%cwd_AG_in(c) - currentPatch%cwd_AG_out(c)
currentPatch%dcwd_BG_dt(c) = currentPatch%cwd_BG_in(c) - currentPatch%cwd_BG_out(c)
enddo
do p = 1,numpft
currentPatch%dleaf_litter_dt(p) = currentPatch%leaf_litter_in(p) - &
currentPatch%leaf_litter_out(p)
currentPatch%droot_litter_dt(p) = currentPatch%root_litter_in(p) - &
currentPatch%root_litter_out(p)
enddo
end subroutine non_canopy_derivs
! ============================================================================
subroutine trim_canopy( currentSite )
!
! !DESCRIPTION:
! Canopy trimming / leaf optimisation. Removes leaves in negative annual carbon balance.
!
! !USES:
!
!
! !ARGUMENTS
type (ed_site_type),intent(inout), target :: currentSite
!
! !LOCAL VARIABLES:
type (ed_cohort_type) , pointer :: currentCohort
type (ed_patch_type) , pointer :: currentPatch
integer :: z ! leaf layer
integer :: ipft ! pft index
integer :: trimmed ! was this layer trimmed in this year? If not expand the canopy.
real(r8) :: tar_bl ! target leaf biomass (leaves flushed, trimmed)
real(r8) :: tar_bfr ! target fine-root biomass (leaves flushed, trimmed)
real(r8) :: bfr_per_bleaf ! ratio of fine root per leaf biomass
!----------------------------------------------------------------------
currentPatch => currentSite%youngest_patch
do while(associated(currentPatch))
currentCohort => currentPatch%tallest
do while (associated(currentCohort))
trimmed = 0
ipft = currentCohort%pft
call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area)
currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%status_coh, currentCohort%pft, &
currentCohort%c_area, currentCohort%n )
currentCohort%treesai = tree_sai(currentCohort%dbh, currentCohort%pft, currentCohort%canopy_trim, &
currentCohort%c_area, currentCohort%n)
currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed)
if (currentCohort%nv > nlevleaf)then
write(fates_log(),*) 'nv > nlevleaf',currentCohort%nv,currentCohort%treelai,currentCohort%treesai, &
currentCohort%c_area,currentCohort%n,currentCohort%bl
endif
call bleaf(currentcohort%dbh,ipft,currentcohort%canopy_trim,tar_bl)
if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then
! only query fine root biomass if using a fine root allometric model that takes leaf trim into account
call bfineroot(currentcohort%dbh,ipft,currentcohort%canopy_trim,tar_bfr)
bfr_per_bleaf = tar_bfr/tar_bl
endif
!Leaf cost vs netuptake for each leaf layer.
do z = 1,nlevleaf
if (currentCohort%year_net_uptake(z) /= 999._r8)then !there was activity this year in this leaf layer.
!Leaf Cost kgC/m2/year-1
!decidous costs.
if (EDPftvarcon_inst%season_decid(ipft) == 1.or. &
EDPftvarcon_inst%stress_decid(ipft) == 1)then
currentCohort%leaf_cost = 1._r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8)
if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then
! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment
! to the leaf increment; otherwise do not.
currentCohort%leaf_cost = currentCohort%leaf_cost + &
1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
endif
currentCohort%leaf_cost = currentCohort%leaf_cost * &
(EDPftvarcon_inst%grperc(ipft) + 1._r8)
else !evergreen costs
currentCohort%leaf_cost = 1.0_r8/(EDPftvarcon_inst%slatop(ipft)* &
EDPftvarcon_inst%leaf_long(ipft)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1
if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then
! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment
! to the leaf increment; otherwise do not.
currentCohort%leaf_cost = currentCohort%leaf_cost + &
1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
endif
currentCohort%leaf_cost = currentCohort%leaf_cost * &
(EDPftvarcon_inst%grperc(ipft) + 1._r8)
endif
if (currentCohort%year_net_uptake(z) < currentCohort%leaf_cost)then
if (currentCohort%canopy_trim > EDPftvarcon_inst%trim_limit(ipft))then
if ( DEBUG ) then
write(fates_log(),*) 'trimming leaves', &
currentCohort%canopy_trim,currentCohort%leaf_cost
endif
! keep trimming until none of the canopy is in negative carbon balance.
if (currentCohort%hite > EDPftvarcon_inst%hgt_min(ipft))then
currentCohort%canopy_trim = currentCohort%canopy_trim - &
EDPftvarcon_inst%trim_inc(ipft)
if (EDPftvarcon_inst%evergreen(ipft) /= 1)then
currentCohort%laimemory = currentCohort%laimemory * &
(1.0_r8 - EDPftvarcon_inst%trim_inc(ipft))
endif
trimmed = 1
endif
endif
endif
endif !leaf activity?
enddo !z
currentCohort%year_net_uptake(:) = 999.0_r8
if (trimmed == 0.and.currentCohort%canopy_trim < 1.0_r8)then
currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(ipft)
endif
if ( DEBUG ) then
write(fates_log(),*) 'trimming',currentCohort%canopy_trim
endif
! currentCohort%canopy_trim = 1.0_r8 !FIX(RF,032414) this turns off ctrim for now.
currentCohort => currentCohort%shorter
enddo
currentPatch => currentPatch%older
enddo
end subroutine trim_canopy
! ============================================================================
subroutine phenology( currentSite, bc_in )
!
! !DESCRIPTION:
! Phenology.
!
! !USES:
use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm
use EDParamsMod, only : ED_val_phen_drought_threshold, ED_val_phen_doff_time
use EDParamsMod, only : ED_val_phen_a, ED_val_phen_b, ED_val_phen_c, ED_val_phen_chiltemp
use EDParamsMod, only : ED_val_phen_mindayson, ED_val_phen_ncolddayslim, ED_val_phen_coldtemp
!
! !ARGUMENTS:
type(ed_site_type), intent(inout), target :: currentSite
type(bc_in_type), intent(in) :: bc_in
!
! !LOCAL VARIABLES:
integer :: t ! day of year
integer :: ncolddays ! no days underneath the threshold for leaf drop
integer :: i
integer :: timesincedleafon,timesincedleafoff,timesinceleafon,timesinceleafoff
integer :: refdate
integer :: curdate
integer :: yr ! year (0, ...)
integer :: mon ! month (1, ..., 12)
integer :: day ! day of month (1, ..., 31)
integer :: sec ! seconds of the day
real(r8) :: gdd_threshold
integer :: ncdstart ! beginning of counting period for chilling degree days.
integer :: gddstart ! beginning of counting period for growing degree days.
real(r8) :: temp_in_C ! daily averaged temperature in celcius
real(r8), parameter :: canopy_leaf_lifespan = 365.0_r8 ! Mean lifespan canopy leaves
! FIX(RGK 07/10/17)
! This is a band-aid on unusual code
! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414)
! - this is arbitrary and poorly understood. Needs work. ED_
!Parameters: defaults from Botta et al. 2000 GCB,6 709-725
!Parameters, default from from SDGVM model of senesence
t = hlm_day_of_year
temp_in_C = bc_in%t_veg24_si - tfrz
!-----------------Cold Phenology--------------------!
!Zero growing degree and chilling day counters
if (currentSite%lat > 0)then
ncdstart = 270 !Northern Hemisphere begining November
gddstart = 1 !Northern Hemisphere begining January
else
ncdstart = 120 !Southern Hemisphere beginning May
gddstart = 181 !Northern Hemisphere begining July
endif
! FIX(SPM,032414) - this will only work for the first year, no?
if (t == ncdstart)then
currentSite%ncd = 0._r8
endif
!Accumulate growing/chilling days after start of counting period
if (temp_in_C < ED_val_phen_chiltemp)then
currentSite%ncd = currentSite%ncd + 1.0_r8
endif
!GDD accumulation function, which also depends on chilling days.
gdd_threshold = ED_val_phen_a + ED_val_phen_b*exp(ED_val_phen_c*currentSite%ncd)
!Accumulate temperature of last 10 days.
currentSite%last_n_days(2:senes) = currentSite%last_n_days(1:senes-1)
currentSite%last_n_days(1) = temp_in_C
!count number of days for leaves off
ncolddays = 0
do i = 1,senes
if (currentSite%last_n_days(i) < ED_val_phen_coldtemp)then
ncolddays = ncolddays + 1
endif
enddo
! Here is where we do the GDD accumulation calculation
!
! reset GDD on set dates
if (t == gddstart)then
currentSite%ED_GDD_site = 0._r8
endif
!
! accumulate the GDD using daily mean temperatures
if (bc_in%t_veg24_si .gt. tfrz) then
currentSite%ED_GDD_site = currentSite%ED_GDD_site + bc_in%t_veg24_si - tfrz
endif
timesinceleafoff = hlm_model_day - currentSite%leafoffdate
!LEAF ON: COLD DECIDUOUS. Needs to
!1) have exceeded the growing degree day threshold
!2) The leaves should not be on already
!3) There should have been at least on chilling day in the counting period.
if (currentSite%ED_GDD_site > gdd_threshold)then
if (currentSite%status == 1) then
if (currentSite%ncd >= 1) then
currentSite%status = 2 !alter status of site to 'leaves on'
! NOTE(bja, 2015-01) should leafondate = model_day to be consistent with leaf off?
currentSite%leafondate = t !record leaf on date
if ( DEBUG ) write(fates_log(),*) 'leaves on'
endif !ncd
endif !status
endif !GDD
timesinceleafon = hlm_model_day - currentSite%leafondate
!LEAF OFF: COLD THRESHOLD
!Needs to:
!1) have exceeded the number of cold days threshold
!2) have exceeded the minimum leafon time.
!3) The leaves should not be off already
!4) The day of the year should be larger than the counting period. (not sure if we need this/if it will break the restarting)
if (ncolddays > ED_val_phen_ncolddayslim)then
if (timesinceleafon > ED_val_phen_mindayson)then
if (currentSite%status == 2)then
currentSite%status = 1 !alter status of site to 'leaves on'
currentSite%leafoffdate = hlm_model_day !record leaf off date
if ( DEBUG ) write(fates_log(),*) 'leaves off'
endif
endif
endif
!LEAF OFF: COLD LIFESPAN THRESHOLD
if(timesinceleafoff > 400)then !remove leaves after a whole year when there is no 'off' period.
if(currentSite%status == 2)then
currentSite%status = 1 !alter status of site to 'leaves on'
currentSite%leafoffdate = hlm_model_day !record leaf off date
if ( DEBUG ) write(fates_log(),*) 'leaves off'
endif
endif
!-----------------Drought Phenology--------------------!
! Principles of drought-deciduos phenology model...
! The 'dstatus' flag is 2 when leaves are on, and 1 when leaves area off.
! The following sets those site-level flags, which are acted on in phenology_deciduos.
! A* The leaves live for either the length of time the soil moisture is over the threshold
! or the lifetime of the leaves, whichever is shorter.
! B*: If the soil is only wet for a very short time, then the leaves stay on for 100 days
! C*: The leaves are only permitted to come ON for a 60 day window around when they last came on,
! to prevent 'flickering' on in response to wet season storms
! D*: We don't allow anything to happen in the first ten days to allow the water memory window to come into equlibirum.
! E*: If the soil is always wet, the leaves come on at the beginning of the window, and then last for their lifespan.
! ISSUES
! 1. It's not clear what water content we should track. Here we are tracking the top layer,
! but we probably should track something like BTRAN,
! but BTRAN is defined for each PFT, and there could potentially be more than one stress-dec PFT.... ?
! 2. In the beginning, the window is set at an arbitrary time of the year, so the leaves might come on
! in the dry season, using up stored reserves
! for the stress-dec plants, and potentially killing them. To get around this, we need to read in the
! 'leaf on' date from some kind of start-up file
! but we would need that to happen for every resolution, etc.
! 3. Will this methodology properly kill off the stress-dec trees where there is no water stress?
! What about where the wet period coincides with the
! warm period? We would just get them overlapping with the cold-dec trees, even though that isn't appropriate....
! Why don't the drought deciduous trees grow
! in the North? Is cold decidousness maybe even the same as drought deciduosness there (and so does this
! distinction actually matter??)....
!Accumulate surface water memory of last 10 days.
do i = 1,numWaterMem-1 !shift memory along one
currentSite%water_memory(numWaterMem+1-i) = currentSite%water_memory(numWaterMem-i)
enddo
currentSite%water_memory(1) = bc_in%h2o_liqvol_sl(1) !waterstate_inst%h2osoi_vol_col(coli,1)
!In drought phenology, we often need to force the leaves to stay on or off as moisture fluctuates...
timesincedleafoff = 0
if (currentSite%dstatus == 1)then !the leaves are off. How long have they been off?
!leaves have come on, but last year, so at a later date than now.
if (currentSite%dleafoffdate > 0.and.currentSite%dleafoffdate > t)then
timesincedleafoff = t + (360 - currentSite%dleafoffdate)
else
timesincedleafoff = t - currentSite%dleafoffdate
endif
endif
timesincedleafon = 0
!the leaves are on. How long have they been on?
if (currentSite%dstatus == 2)then
!leaves have come on, but last year, so at a later date than now.
if (currentSite%dleafondate > 0.and.currentSite%dleafondate > t)then
timesincedleafon = t + (360 - currentSite%dleafondate)
else
timesincedleafon = t - currentSite%dleafondate
endif
endif
!LEAF ON: DROUGHT DECIDUOUS WETNESS
!Here, we used a window of oppurtunity to determine if we are close to the time when then leaves came on last year
if ((t >= currentSite%dleafondate - 30.and.t <= currentSite%dleafondate + 30).or.(t > 360 - 15.and. &
currentSite%dleafondate < 15))then ! are we in the window?
! TODO: CHANGE THIS MATH, MOVE THE DENOMENATOR OUTSIDE OF THE SUM (rgk 01-2017)
if (sum(currentSite%water_memory(1:numWaterMem)/dble(numWaterMem)) &
>= ED_val_phen_drought_threshold.and.currentSite%dstatus == 1.and.t >= 10)then
! leave some minimum time between leaf off and leaf on to prevent 'flickering'.
if (timesincedleafoff > ED_val_phen_doff_time)then
currentSite%dstatus = 2 !alter status of site to 'leaves on'
currentSite%dleafondate = t !record leaf on date
endif
endif
endif
!we still haven't done budburst by end of window
if (t == currentSite%dleafondate+30.and.currentSite%dstatus == 1)then
currentSite%dstatus = 2 ! force budburst!
currentSite%dleafondate = t ! record leaf on date
endif
!LEAF OFF: DROUGHT DECIDUOUS LIFESPAN - if the leaf gets to the end of its useful life. A*, E*
if (currentSite%dstatus == 2.and.t >= 10)then !D*
!Are the leaves at the end of their lives?
!FIX(RF,0401014)- this is hardwiring....
!FIX(RGK:changed from hard-coded pft 7 leaf lifespan to labeled constant (1 year)
if ( timesincedleafon > canopy_leaf_lifespan )then
currentSite%dstatus = 1 !alter status of site to 'leaves on'
currentSite%dleafoffdate = t !record leaf on date
endif
endif
!LEAF OFF: DROUGHT DECIDUOUS DRYNESS - if the soil gets too dry, and the leaves have already been on a while...
if (currentSite%dstatus == 2.and.t >= 10)then !D*
if (sum(currentSite%water_memory(1:10)/10._r8) <= ED_val_phen_drought_threshold)then
if (timesincedleafon > 100)then !B* Have the leaves been on for some reasonable length of time? To prevent flickering.
currentSite%dstatus = 1 !alter status of site to 'leaves on'
currentSite%dleafoffdate = t !record leaf on date
endif
endif
endif
call phenology_leafonoff(currentSite)
end subroutine phenology
! ============================================================================
subroutine phenology_leafonoff(currentSite)
!
! !DESCRIPTION:
! Controls the leaf on and off economics
!
! !USES:
!
! !ARGUMENTS:
type(ed_site_type), intent(inout), target :: currentSite
!
! !LOCAL VARIABLES:
type(ed_patch_type) , pointer :: currentPatch
type(ed_cohort_type), pointer :: currentCohort
real(r8) :: store_output ! the amount of the store to put into leaves -
! is a barrier against negative storage and C starvation.
!------------------------------------------------------------------------
currentPatch => CurrentSite%oldest_patch
store_output = 0.5_r8
do while(associated(currentPatch))
currentCohort => currentPatch%tallest
do while(associated(currentCohort))
currentCohort%leaf_litter = 0.0_r8 !zero leaf litter for today.
!COLD LEAF ON
if (EDPftvarcon_inst%season_decid(currentCohort%pft) == 1)then
if (currentSite%status == 2)then !we have just moved to leaves being on .
if (currentCohort%status_coh == 1)then !Are the leaves currently off?
currentCohort%status_coh = 2 !Leaves are on, so change status to stop flow of carbon out of bstore.
if (currentCohort%laimemory <= currentCohort%bstore)then
currentCohort%bl = currentCohort%laimemory !extract stored carbon to make new leaves.
else
! we can only put on as much carbon as there is in the store...
! nb. Putting all of bstore into leaves is C-starvation suicidal.
! The tendency for this could be parameterized
currentCohort%bl = currentCohort%bstore * store_output
endif
if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 1 ',currentCohort%bstore
currentCohort%bstore = currentCohort%bstore - currentCohort%bl ! Drain store
if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 2 ',currentCohort%bstore
currentCohort%laimemory = 0.0_r8
endif !pft phenology
endif ! growing season
!COLD LEAF OFF
! currentCohort%leaf_litter = 0.0_r8 !zero leaf litter for today.
if (currentSite%status == 1)then !past leaf drop day? Leaves still on tree?
if (currentCohort%status_coh == 2)then ! leaves have not dropped
currentCohort%status_coh = 1
!remember what the lai was this year to put the same amount back on in the spring...
currentCohort%laimemory = currentCohort%bl
! add lost carbon to litter
currentCohort%leaf_litter = currentCohort%bl
currentCohort%bl = 0.0_r8
endif !leaf status
endif !currentSite status
endif !season_decid
!DROUGHT LEAF ON
if (EDPftvarcon_inst%stress_decid(currentCohort%pft) == 1)then
if (currentSite%dstatus == 2)then !we have just moved to leaves being on .
if (currentCohort%status_coh == 1)then !is it the leaf-on day? Are the leaves currently off?
currentCohort%status_coh = 2 !Leaves are on, so change status to stop flow of carbon out of bstore.
if (currentCohort%laimemory <= currentCohort%bstore)then
currentCohort%bl = currentCohort%laimemory !extract stored carbon to make new leaves.
else
!we can only put on as much carbon as there is in the store.
currentCohort%bl = currentCohort%bstore * store_output
endif
if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 3 ',currentCohort%bstore
currentCohort%bstore = currentCohort%bstore - currentCohort%bl ! empty store
if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 4 ',currentCohort%bstore
currentCohort%laimemory = 0.0_r8
endif !currentCohort status again?
endif !currentSite status
!DROUGHT LEAF OFF
if (currentSite%dstatus == 1)then
if (currentCohort%status_coh == 2)then ! leaves have not dropped
currentCohort%status_coh = 1
currentCohort%laimemory = currentCohort%bl
! add retranslocated carbon (very small) to store.
currentCohort%bstore = currentCohort%bstore
! add falling leaves to litter pools . convert to KgC/m2
currentCohort%leaf_litter = currentCohort%bl
currentCohort%bl = 0.0_r8
endif
endif !status
endif !drought dec.
currentCohort => currentCohort%shorter
enddo !currentCohort
currentPatch => currentPatch%younger
enddo !currentPatch
end subroutine phenology_leafonoff
! ============================================================================
subroutine seeds_in( currentSite, cp_pnt )
!
! !DESCRIPTION:
! Flux from plants into seed pool.
!
! !USES:
use EDTypesMod, only : AREA
use EDTypesMod, only : homogenize_seed_pfts
!
! !ARGUMENTS
type(ed_site_type), intent(inout), target :: currentSite
type(ed_patch_type), intent(inout), target :: cp_pnt ! seeds go to these patches.
!
! !LOCAL VARIABLES:
type(ed_patch_type), pointer :: currentPatch
type(ed_cohort_type), pointer :: currentCohort
integer :: p
logical :: pft_present(maxpft)
real(r8) :: npfts_present
!----------------------------------------------------------------------
currentPatch => cp_pnt
currentPatch%seeds_in(:) = 0.0_r8
if ( homogenize_seed_pfts ) then
! special mode to remove intergenerational filters on PFT existence: each PFT seeds all PFTs
! first loop over all patches and cohorts to see what and how many PFTs are present on this site
pft_present(:) = .false.
npfts_present = 0._r8
currentPatch => currentSite%oldest_patch
do while(associated(currentPatch))
currentCohort => currentPatch%tallest
do while (associated(currentCohort))
p = currentCohort%pft
if (.not. pft_present(p)) then
pft_present(p) = .true.
npfts_present = npfts_present + 1._r8
endif
currentCohort => currentCohort%shorter
enddo !cohort loop
currentPatch => currentPatch%younger
enddo ! patch loop
! now calculate the homogenized seed flux into each PFT pool
currentPatch => cp_pnt
currentCohort => currentPatch%tallest
do while (associated(currentCohort))
do p = 1, numpft
if (pft_present(p)) then
currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + currentCohort%seed_prod * currentCohort%n / &
(currentPatch%area * npfts_present)
endif
end do
currentCohort => currentCohort%shorter
enddo !cohort loop
else
! normal case: each PFT seeds its own type
currentCohort => currentPatch%tallest
do while (associated(currentCohort))
p = currentCohort%pft
currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + &
currentCohort%seed_prod * currentCohort%n/currentPatch%area
currentCohort => currentCohort%shorter
enddo !cohort loop
endif
currentPatch => currentSite%oldest_patch
do while(associated(currentPatch))
if (external_recruitment == 1) then !external seed rain - needed to prevent extinction
do p = 1,numpft
currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + &
EDPftvarcon_inst%seed_rain(p) !KgC/m2/year
currentSite%seed_rain_flux(p) = currentSite%seed_rain_flux(p) + &
EDPftvarcon_inst%seed_rain(p) * currentPatch%area/AREA !KgC/m2/year
enddo
endif
currentPatch => currentPatch%younger
enddo
end subroutine seeds_in
! ============================================================================
subroutine seed_decay( currentSite, currentPatch )
!
! !DESCRIPTION:
! Flux from seed pool into leaf litter pool
!
! !USES:
use EDPftvarcon , only : EDPftvarcon_inst
!
! !ARGUMENTS
type(ed_site_type), intent(inout), target :: currentSite
type(ed_patch_type),intent(inout) :: currentPatch ! seeds go to these patches.
!
! !LOCAL VARIABLES:
integer :: p
!----------------------------------------------------------------------
! default value from Liscke and Loffler 2006 ; making this a PFT-specific parameter
! decays the seed pool according to exponential model
! seed_decay_turnover is in yr-1
do p = 1,numpft
currentPatch%seed_decay(p) = currentSite%seed_bank(p) * EDPftvarcon_inst%seed_decay_turnover(p)
enddo
end subroutine seed_decay
! ============================================================================
subroutine seed_germination( currentSite, currentPatch )
!
! !DESCRIPTION:
! Flux from seed pool into sapling pool
!
! !USES:
use EDPftvarcon , only : EDPftvarcon_inst
!
! !ARGUMENTS
type(ed_site_type), intent(inout), target :: currentSite
type(ed_patch_type),intent(inout) :: currentPatch ! seeds go to these patches.
!
! !LOCAL VARIABLES:
integer :: p
real(r8) max_germination !cap on germination rates. KgC/m2/yr Lishcke et al. 2009
!----------------------------------------------------------------------
max_germination = 1.0_r8 !this is arbitrary
! germination_timescale is being pulled to PFT parameter; units are 1/yr
! thus the mortality rate of seed -> recruit (in units of carbon)
! is seed_decay_turnover(p)/germination_timescale(p)
! and thus the mortlaity rate (in units of individuals) is the product of
! that times the ratio of (hypothetical) seed mass to recruit biomass
do p = 1,numpft
currentPatch%seed_germination(p) = min(currentSite%seed_bank(p) * &
EDPftvarcon_inst%germination_timescale(p),max_germination)
enddo
end subroutine seed_germination
! ============================================================================
subroutine PlantGrowth( currentSite, currentCohort, bc_in )
!
! !DESCRIPTION:
! Main subroutine for plant allocation and growth
!
! !USES:
! Original: Rosie Fisher
! Updated: Ryan Knox
use FatesInterfaceMod, only : hlm_use_ed_prescribed_phys
!
! !ARGUMENTS
type(ed_site_type), intent(inout), target :: currentSite
type(ed_cohort_type),intent(inout), target :: currentCohort
type(bc_in_type), intent(in) :: bc_in
!
! !LOCAL VARIABLES:
integer :: ipft ! PFT index
real(r8) :: carbon_balance ! daily carbon balance for this cohort
! Per plant allocation targets
real(r8) :: bt_leaf ! leaf biomass (kgC)
real(r8) :: dbt_leaf_dd ! change in leaf biomass wrt diameter (kgC/cm)
real(r8) :: bt_fineroot ! fine root biomass (kgC)
real(r8) :: dbt_fineroot_dd ! change in fine root biomass wrt diameter (kgC/cm)
real(r8) :: bt_sap ! sapwood biomass (kgC)
real(r8) :: dbt_sap_dd ! change in sapwood biomass wrt diameter (kgC/cm)
real(r8) :: bt_agw ! above ground biomass (kgC/cm)
real(r8) :: dbt_agw_dd ! change in above ground biomass wrt diameter (kgC/cm)
real(r8) :: bt_bgw ! coarse root biomass (kgC)
real(r8) :: dbt_bgw_dd ! change in coarse root biomass (kgC/cm)
real(r8) :: bt_dead ! dead (structural) biomass (kgC)
real(r8) :: dbt_dead_dd ! change in dead biomass wrt diameter (kgC/cm)
real(r8) :: bt_store ! target storage biomass (kgC)
real(r8) :: dbt_store_dd ! target rate of change in storage (kgC/cm)
real(r8) :: dbt_total_dd ! total target biomass rate of change (kgC/cm)
real(r8) :: leaf_below_target ! leaf biomass below target amount [kgC]
real(r8) :: froot_below_target ! fineroot biomass below target amount [kgC]
real(r8) :: sap_below_target ! sapwood biomass below target amount [kgC]
real(r8) :: store_below_target ! storage biomass below target amount [kgC]
real(r8) :: dead_below_target ! dead (structural) biomass below target amount [kgC]
real(r8) :: total_below_target ! total biomass below the allometric target [kgC]
real(r8) :: bstore_flux ! carbon fluxing into storage [kgC]
real(r8) :: bl_flux ! carbon fluxing into leaves [kgC]
real(r8) :: br_flux ! carbon fluxing into fineroots [kgC]
real(r8) :: bsw_flux ! carbon fluxing into sapwood [kgC]
real(r8) :: bdead_flux ! carbon fluxing into structure [kgC]
real(r8) :: brepro_flux ! carbon fluxing into reproductive tissues [kgC]
real(r8) :: flux_adj ! adjustment made to growth flux term to minimize error [kgC]
real(r8) :: store_target_fraction ! ratio between storage and leaf biomass when on allometry [kgC]
real(r8) :: repro_fraction ! fraction of carbon gain sent to reproduction when on-allometry
real(r8) :: leaf_turnover_demand ! leaf carbon that is demanded to replace maintenance turnover [kgC]
real(r8) :: root_turnover_demand ! fineroot carbon that is demanded to replace
! maintenance turnover [kgC]
real(r8) :: total_turnover_demand ! total carbon that is demanded to replace maintenance turnover [kgC]
real(r8),dimension(n_cplantpools) :: c_pool ! Vector of carbon pools passed to integrator
real(r8),dimension(n_cplantpools) :: c_pool_out ! Vector of carbon pools passed back from integrator
logical,dimension(n_cplantpools) :: c_mask ! Mask of active pools during integration
logical :: step_pass ! Did the integration step pass?
logical :: grow_leaf ! Are leaves at allometric target and should be grown?
logical :: grow_froot ! Are fine-roots at allometric target and should be grown?
logical :: grow_sap ! Is sapwood at allometric target and should be grown?
logical :: grow_store ! Is storage at allometric target and should be grown?
! integrator variables
real(r8) :: deltaC ! trial value for substep
integer :: ierr ! error flag for allometric growth step
integer :: nsteps ! number of sub-steps
integer :: istep ! current substep index
real(r8) :: totalC ! total carbon allocated over alometric growth step
real(r8) :: dbh_sub ! substep dbh
real(r8) :: h_sub ! substep h
real(r8) :: bl_sub ! substep leaf biomass
real(r8) :: br_sub ! substep root biomass
real(r8) :: bsw_sub ! substep sapwood biomass
real(r8) :: bstore_sub ! substep storage biomass
real(r8) :: bdead_sub ! substep structural biomass
real(r8) :: brepro_sub ! substep reproductive biomass
! Woody turnover timescale [years]
real(r8), parameter :: cbal_prec = 1.0e-15_r8 ! Desired precision in carbon balance
! non-integrator part
integer , parameter :: max_substeps = 300 ! Number of step attempts before
! giving up
real(r8), parameter :: max_trunc_error = 1.0_r8 ! allowable numerical truncation error
integer, parameter :: ODESolve = 2 ! 1=RKF45, 2=Euler
ipft = currentCohort%pft
! Initialize seed production
currentCohort%seed_prod = 0.0_r8
! Initialize NPP flux diagnostics
currentCohort%npp_stor = 0.0_r8
currentCohort%npp_leaf = 0.0_r8
currentCohort%npp_fnrt = 0.0_r8
currentCohort%npp_dead = 0.0_r8
currentCohort%npp_seed = 0.0_r8
currentCohort%npp_sapw = 0.0_r8
! Initialize rates of change
currentCohort%dhdt = 0.0_r8
currentCohort%dbdeaddt = 0.0_r8
currentCohort%dbstoredt = 0.0_r8
currentCohort%ddbhdt = 0.0_r8
! If the cohort has grown, it is not new
currentCohort%isnew=.false.
! -----------------------------------------------------------------------------------
! I. Identify the net carbon gain for this dynamics interval
! Set the available carbon pool, identify allocation portions, and decrement
! the available carbon pool to zero.
! -----------------------------------------------------------------------------------
! convert from kgC/indiv/day into kgC/indiv/year
! <x>_acc_hold is remembered until the next dynamics step (used for I/O)
! <x>_acc will be reset soon and will be accumulated on the next leaf photosynthesis
! step
if (hlm_use_ed_prescribed_phys .eq. itrue) then
if (currentCohort%canopy_layer .eq. 1) then
currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_canopy(ipft) &
* currentCohort%c_area / currentCohort%n
! add these for balance checking purposes
currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year
else
currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_understory(ipft) &
* currentCohort%c_area / currentCohort%n
! add these for balance checking purposes
currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year
endif
else
currentCohort%npp_acc_hold = currentCohort%npp_acc * dble(hlm_days_per_year)
currentCohort%gpp_acc_hold = currentCohort%gpp_acc * dble(hlm_days_per_year)
currentCohort%resp_acc_hold = currentCohort%resp_acc * dble(hlm_days_per_year)
endif
currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n
! Available carbon for growth [kgC]
carbon_balance = currentCohort%npp_acc
! -----------------------------------------------------------------------------------
! II. Calculate target size of living biomass compartment for a given dbh.
! -----------------------------------------------------------------------------------
! Target sapwood biomass and deriv. according to allometry and trimming [kgC, kgC/cm]
call bsap_allom(currentCohort%dbh,ipft,currentCohort%canopy_trim,bt_sap,dbt_sap_dd)
! Target total above ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm]
call bagw_allom(currentCohort%dbh,ipft,bt_agw,dbt_agw_dd)
! Target total below ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm]
call bbgw_allom(currentCohort%dbh,ipft,bt_bgw,dbt_bgw_dd)
! Target total dead (structrual) biomass and deriv. [kgC, kgC/cm]
call bdead_allom( bt_agw, bt_bgw, bt_sap, ipft, bt_dead, &
dbt_agw_dd, dbt_bgw_dd, dbt_sap_dd, dbt_dead_dd )
! ------------------------------------------------------------------------------------
! If structure is larger than target, then we need to correct some integration errors
! by slightly increasing dbh to match it.
! -----------------------------------------------------------------------------------
if( ((currentCohort%bdead-bt_dead) > calloc_abs_error) .and. &
(EDPftvarcon_inst%woody(ipft) == itrue) ) then
call StructureResetOfDH( currentCohort%bdead, ipft, &
currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite )
! Re-calculate the sapwood and structural wood targets based on the new dbh
! ------------------------------------------------------------------------------------------
! Target sapwood biomass and deriv. according to allometry and trimming [kgC, kgC/cm]
call bsap_allom(currentCohort%dbh,ipft,currentCohort%canopy_trim,bt_sap,dbt_sap_dd)
! Target total above ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm]
call bagw_allom(currentCohort%dbh,ipft,bt_agw,dbt_agw_dd)
! Target total below ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm]
call bbgw_allom(currentCohort%dbh,ipft,bt_bgw,dbt_bgw_dd)
! Target total dead (structrual) biomass and deriv. [kgC, kgC/cm]
call bdead_allom( bt_agw, bt_bgw, bt_sap, ipft, bt_dead, &
dbt_agw_dd, dbt_bgw_dd, dbt_sap_dd, dbt_dead_dd )
end if
! Target leaf biomass according to allometry and trimming
call bleaf(currentCohort%dbh,ipft,currentCohort%canopy_trim,bt_leaf,dbt_leaf_dd)
! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm]