-
Notifications
You must be signed in to change notification settings - Fork 92
/
FatesPlantRespPhotosynthMod.F90
2513 lines (2015 loc) · 126 KB
/
FatesPlantRespPhotosynthMod.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 FATESPlantRespPhotosynthMod
!-------------------------------------------------------------------------------------
! !DESCRIPTION:
! Calculates the plant respiration and photosynthetic fluxes for the FATES model
! This code is similar to and was originally based off of the 'photosynthesis'
! subroutine in the CLM model.
!
! Parameter for activation and deactivation energies were taken from:
! Activation energy, from:
! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259
! Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430
! except TPU from: Harley et al (1992) Plant, Cell and Environment 15:271-282
! High temperature deactivation, from:
! Leuning (2002) Plant, Cell and Environment 25:1205-1210
! The factor "c" scales the deactivation to a value of 1.0 at 25C
! Photosynthesis and stomatal conductance parameters, from:
! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593
! ------------------------------------------------------------------------------------
! !USES:
use FatesGlobals, only : endrun => fates_endrun
use FatesGlobals, only : fates_log
use FatesGlobals, only : FatesWarn,N2S,A2S,I2S
use FatesConstantsMod, only : r8 => fates_r8
use FatesConstantsMod, only : itrue
use FatesConstantsMod, only : nearzero
use FatesConstantsMod, only : molar_mass_ratio_vapdry
use FatesConstantsMod, only : molar_mass_water
use FatesConstantsMod, only : rgas_J_K_mol
use FatesConstantsMod, only : fates_unset_r8
use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm
use FatesConstantsMod, only : nocomp_bareground
use FatesConstantsMod, only : photosynth_acclim_model_none
use FatesConstantsMod, only : photosynth_acclim_model_kumarathunge_etal_2019
use FatesInterfaceTypesMod, only : hlm_use_planthydro
use FatesInterfaceTypesMod, only : hlm_parteh_mode
use FatesInterfaceTypesMod, only : numpft
use FatesInterfaceTypesMod, only : nleafage
use FatesUtilsMod, only : QuadraticRoots => QuadraticRootsSridharachary
use EDParamsMod, only : maxpft
use EDParamsMod, only : nlevleaf
use EDParamsMod, only : nclmax
use PRTGenericMod, only : max_nleafage
use EDTypesMod, only : do_fates_salinity
use EDParamsMod, only : q10_mr
use FatesPatchMod, only : fates_patch_type
use FatesCohortMod, only : fates_cohort_type
use EDParamsMod, only : maintresp_leaf_model
use FatesConstantsMod, only : lmrmodel_ryan_1991
use FatesConstantsMod, only : lmrmodel_atkin_etal_2017
use PRTGenericMod, only : prt_carbon_allom_hyp
use PRTGenericMod, only : prt_cnp_flex_allom_hyp
use PRTGenericMod, only : carbon12_element
use PRTGenericMod, only : nitrogen_element
use PRTGenericMod, only : leaf_organ
use PRTGenericMod, only : fnrt_organ
use PRTGenericMod, only : sapw_organ
use PRTGenericMod, only : store_organ
use PRTGenericMod, only : repro_organ
use PRTGenericMod, only : struct_organ
use EDParamsMod, only : maintresp_nonleaf_baserate
use EDParamsMod, only : stomatal_model
use EDParamsMod, only : stomatal_assim_model
use EDParamsMod, only : dayl_switch
use EDParamsMod, only : photo_tempsens_model
use PRTParametersMod, only : prt_params
use EDPftvarcon , only : EDPftvarcon_inst
use TemperatureType, only : temperature_type
use FatesRadiationMemMod, only : norman_solver,twostr_solver
use EDParamsMod, only : radiation_model
use FatesRadiationMemMod, only : ipar
use FatesTwoStreamUtilsMod, only : FatesGetCohortAbsRad
use FatesAllometryMod , only : VegAreaLayer
use FatesAllometryMod, only : decay_coeff_vcmax
! CIME Globals
use shr_log_mod , only : errMsg => shr_log_errMsg
implicit none
private
public :: FatesPlantRespPhotosynthDrive ! Called by the HLM-Fates interface
character(len=*), parameter, private :: sourcefile = &
__FILE__
character(len=1024) :: warn_msg ! for defining a warning message
!-------------------------------------------------------------------------------------
! maximum stomatal resistance [s/m] (used across several procedures)
real(r8),parameter :: rsmax0 = 2.e8_r8
logical :: debug = .false.
!-------------------------------------------------------------------------------------
! Ratio of H2O/CO2 gas diffusion in stomatal airspace (approximate)
real(r8),parameter :: h2o_co2_stoma_diffuse_ratio = 1.6_r8
! Ratio of H2O/CO2 gass diffusion in the leaf boundary layer (approximate)
real(r8),parameter :: h2o_co2_bl_diffuse_ratio = 1.4_r8
! Constants used to define C3 versus C4 photosynth pathways
integer, parameter :: c3_path_index = 1
integer, parameter :: c4_path_index = 0
! Constants used to define conductance models
integer, parameter :: medlyn_model = 2
integer, parameter :: ballberry_model = 1
! Alternatively, Gross Assimilation can be used to estimate
! leaf co2 partial pressure and therefore conductance. The default
! is to use anet
integer, parameter :: net_assim_model = 1
integer, parameter :: gross_assim_model = 2
logical, parameter :: preserve_b4b = .true.
contains
!--------------------------------------------------------------------------------------
subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
! -----------------------------------------------------------------------------------
! !DESCRIPTION:
! Leaf photosynthesis and stomatal conductance calculation as described by
! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 and extended to
! a multi-layer canopy
! -----------------------------------------------------------------------------------
! !USES:
use EDTypesMod , only : ed_site_type
use EDParamsMod , only : dinc_vai
use EDParamsMod , only : dlower_vai
use FatesInterfaceTypesMod , only : bc_in_type
use FatesInterfaceTypesMod , only : bc_out_type
use EDCanopyStructureMod, only : calc_areaindex
use FatesConstantsMod, only : umolC_to_kgC
use FatesConstantsMod, only : umol_per_mmol
use FatesConstantsMod, only : rgas => rgas_J_K_kmol
use FatesParameterDerivedMod, only : param_derived
use FatesAllometryMod, only : bleaf, bstore_allom
use FatesAllometryMod, only : storage_fraction_of_target
use FatesAllometryMod, only : set_root_fraction
use DamageMainMod, only : GetCrownReduction
use FatesInterfaceTypesMod, only : hlm_use_tree_damage
! ARGUMENTS:
! -----------------------------------------------------------------------------------
integer,intent(in) :: nsites
type(ed_site_type),intent(inout),target :: sites(nsites)
type(bc_in_type),intent(in) :: bc_in(nsites)
type(bc_out_type),intent(inout) :: bc_out(nsites)
real(r8),intent(in) :: dtime
! LOCAL VARIABLES:
! -----------------------------------------------------------------------------------
type (fates_patch_type) , pointer :: currentPatch
type (fates_cohort_type), pointer :: currentCohort
! -----------------------------------------------------------------------------------
! These three arrays hold leaf-level biophysical rates that are calculated
! in one loop and then sent to the cohorts in another loop. If hydraulics are
! on, we calculate a unique solution for each level-cohort-layer combination.
! If we are not using hydraulics, we calculate a unique solution for each
! level-pft-layer combination. Thus the following three arrays are statically
! allocated for the maximum space of the two cases (numCohortsPerPatch)
! The "_z" suffix indicates these variables are discretized at the "leaf_layer"
! scale.
! Note: For these temporary arrays, we have the leaf layer dimension first
! and the canopy layer last. This order is chosen for efficiency. The arrays
! such as leaf area that are bound to the patch structure DO NOT follow this order
! as they are used in many other parts of the code with different looping, we
! are not modifying its order now.
! -----------------------------------------------------------------------------------
! leaf maintenance (dark) respiration [umol CO2/m**2/s]
real(r8) :: lmr_z(nlevleaf,maxpft,nclmax)
! stomatal resistance [s/m]
real(r8) :: rs_z(nlevleaf,maxpft,nclmax)
! net leaf photosynthesis averaged over sun and shade leaves. [umol CO2/m**2/s]
real(r8) :: anet_av_z(nlevleaf,maxpft,nclmax)
! Mask used to determine which leaf-layer biophysical rates have been
! used already
logical :: rate_mask_z(nlevleaf,maxpft,nclmax)
real(r8) :: vcmax_z ! leaf layer maximum rate of carboxylation
! (umol co2/m**2/s)
real(r8) :: jmax_z ! leaf layer maximum electron transport rate
! (umol electrons/m**2/s)
real(r8) :: kp_z ! leaf layer initial slope of CO2 response
! curve (C4 plants)
real(r8) :: c13disc_z(nclmax,maxpft,nlevleaf) ! carbon 13 in newly assimilated carbon at leaf level
real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa)
real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa)
real(r8) :: co2_cpoint ! CO2 compensation point (Pa)
real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1)
real(r8) :: stomatal_intercept_btran ! water-stressed minimum stomatal conductance (umol H2O/m**2/s)
real(r8) :: kn ! leaf nitrogen decay coefficient
real(r8) :: cf ! s m**2/umol -> s/m (ideal gas conversion) [umol/m3]
real(r8) :: gb_mol ! leaf boundary layer conductance (molar form: [umol /m**2/s])
real(r8) :: ceair ! vapor pressure of air, constrained (Pa)
real(r8) :: nscaler ! leaf nitrogen scaling coefficient
real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass
real(r8) :: tcsoi ! Temperature response function for root respiration.
real(r8) :: tcwood ! Temperature response function for wood
real(r8) :: patch_la ! exposed leaf area (patch scale)
real(r8) :: live_stem_n ! Live stem (above-ground sapwood)
! nitrogen content (kgN/plant)
real(r8) :: live_croot_n ! Live coarse root (below-ground sapwood)
! nitrogen content (kgN/plant)
real(r8) :: sapw_c ! Sapwood carbon (kgC/plant)
real(r8) :: store_c_target ! Target storage carbon (kgC/plant)
real(r8) :: fnrt_c ! Fine root carbon (kgC/plant)
real(r8) :: fnrt_n ! Fine root nitrogen content (kgN/plant)
real(r8) :: leaf_c ! Leaf carbon (kgC/plant)
real(r8) :: leaf_n ! leaf nitrogen content (kgN/plant)
real(r8) :: g_sb_leaves ! Mean combined (stomata+boundary layer) leaf conductance [m/s]
! over all of the patch's leaves. The "sb" refers to the combined
! "s"tomatal and "b"oundary layer.
! This quantity is relevant on leaf surfaces. It does not
! have units of /m2 leaf per say, but is implicitly on leaf surfaces
real(r8) :: r_sb_leaves ! Mean leaf resistance over all the patch's leaves [s/m]
! This is the direct reciprocal of g_sb_leaves
real(r8) :: r_stomata ! Mean stomatal resistance across all leaves in the patch [s/m]
real(r8) :: maintresp_reduction_factor ! factor by which to reduce maintenance
! respiration when storage pools are low
real(r8) :: b_leaf ! leaf biomass kgC
real(r8) :: frac ! storage pool as a fraction of target leaf biomass
! over each cohort x layer.
real(r8) :: cohort_eleaf_area ! This is the effective leaf area [m2] reported by each cohort
real(r8) :: lnc_top ! Leaf nitrogen content per unit area at canopy top [gN/m2]
real(r8) :: lmr25top ! canopy top leaf maint resp rate at 25C
! for this plant or pft (umol CO2/m**2/s)
real(r8) :: leaf_inc ! LAI-only portion of the vegetation increment of dinc_vai
real(r8) :: lai_canopy_above ! the LAI in the canopy layers above the layer of interest
real(r8) :: lai_layers_above ! the LAI in the leaf layers, within the current canopy,
! above the leaf layer of interest
real(r8) :: lai_current ! the LAI in the current leaf layer
real(r8) :: cumulative_lai ! the cumulative LAI, top down, to the leaf layer of interest
real(r8) :: leaf_psi ! leaf xylem matric potential [MPa] (only meaningful/used w/ hydro)
real(r8) :: fnrt_mr_layer ! fine root maintenance respiation per layer [kgC/plant/s]
real(r8) :: fnrt_mr_nfix_layer ! fineroot maintenance respiration specifically for symbiotic fixation [kgC/plant/layer/s]
real(r8) :: nfix_layer ! Nitrogen fixed in each layer this timestep [kgN/plant/layer/timestep]
real(r8), allocatable :: rootfr_ft(:,:) ! Root fractions per depth and PFT
real(r8) :: agb_frac ! fraction of biomass aboveground
real(r8) :: branch_frac ! fraction of aboveground woody biomass in branches
real(r8) :: crown_reduction ! reduction in crown biomass from damage
real(r8) :: sapw_c_bgw ! belowground sapwood
real(r8) :: sapw_c_agw ! aboveground sapwood
real(r8) :: sapw_c_undamaged ! the target sapwood of an undamaged tree
real(r8) :: sapw_n ! sapwood nitrogen
real(r8) :: sapw_n_bgw ! nitrogen in belowground portion of sapwood
real(r8) :: sapw_n_agw ! nitrogen in aboveground portion of sapwood
real(r8) :: sapw_n_undamaged ! nitrogen in sapwood of undamaged tree
real(r8) :: rd_abs_leaf, rb_abs_leaf, r_abs_stem, r_abs_snow, rb_abs, rd_abs
real(r8) :: fsun
real(r8) :: par_per_sunla, par_per_shala ! PAR per sunlit and shaded leaf area [W/m2 leaf]
real(r8),dimension(75) :: cohort_vaitop
real(r8),dimension(75) :: cohort_vaibot
real(r8),dimension(75) :: cohort_layer_elai
real(r8),dimension(75) :: cohort_layer_esai
real(r8),dimension(75) :: cohort_layer_tlai
real(r8),dimension(75) :: cohort_layer_tsai
real(r8) :: cohort_elai
real(r8) :: cohort_esai
real(r8) :: laisun,laisha
real(r8) :: canopy_area
real(r8) :: elai
! -----------------------------------------------------------------------------------
! Keeping these two definitions in case they need to be added later
!
! -----------------------------------------------------------------------------------
!real(r8) :: psncanopy_pa ! patch sunlit leaf photosynthesis (umol CO2 /m**2/ s)
!real(r8) :: lmrcanopy_pa ! patch sunlit leaf maintenance respiration rate (umol CO2/m**2/s)
integer :: cl,s,iv,j,ps,ft,ifp ! indices
integer :: nv ! number of leaf layers
integer :: NCL_p ! number of canopy layers in patch
integer :: iage ! loop counter for leaf age classes
! Parameters
!
! Base rate is at 20C. Adjust to 25C using the CN Q10 = 1.5
! (gC/gN/s)
! ------------------------------------------------------------------------
! -----------------------------------------------------------------------------------
! Photosynthesis and stomatal conductance parameters, from:
! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593
! -----------------------------------------------------------------------------------
associate( &
c3psn => EDPftvarcon_inst%c3psn , &
slatop => prt_params%slatop , & ! specific leaf area at top of canopy,
! projected area basis [m^2/gC]
woody => prt_params%woody, & ! Is vegetation woody or not?
stomatal_intercept => EDPftvarcon_inst%stomatal_intercept ) !Unstressed minimum stomatal conductance
do s = 1,nsites
! Multi-layer parameters scaled by leaf nitrogen profile.
! Loop through each canopy layer to calculate nitrogen profile using
! cumulative lai at the midpoint of the layer
! Pre-process some variables that are PFT dependent
! but not environmentally dependent
! ------------------------------------------------------------------------
allocate(rootfr_ft(numpft, bc_in(s)%nlevsoil))
do ft = 1,numpft
call set_root_fraction(rootfr_ft(ft,:), ft, &
bc_in(s)%zi_sisl, &
bc_in(s)%max_rooting_depth_index_col)
end do
ifp = 0
currentpatch => sites(s)%oldest_patch
do while (associated(currentpatch))
if_notbare: if(currentpatch%nocomp_pft_label.ne.nocomp_bareground)then
ifp = ifp+1
NCL_p = currentPatch%NCL_p
! Part I. Zero output boundary conditions
! ---------------------------------------------------------------------------
bc_out(s)%rssun_pa(ifp) = 0._r8
bc_out(s)%rssha_pa(ifp) = 0._r8
g_sb_leaves = 0._r8
patch_la = 0._r8
! Part II. Filter out patches
! Patch level filter flag for photosynthesis calculations
! has a short memory, flags:
! 1 = patch has not been called
! 2 = patch is currently marked for photosynthesis
! 3 = patch has been called for photosynthesis already
! ---------------------------------------------------------------------------
if_filter2: if(bc_in(s)%filter_photo_pa(ifp)==2)then
! Part III. Calculate the number of sublayers for each pft and layer.
! And then identify which layer/pft combinations have things in them.
! Output:
! currentPatch%ncan(:,:)
! currentPatch%canopy_mask(:,:)
call UpdateCanopyNCanNRadPresent(currentPatch)
! Part IV. Identify some environmentally derived parameters:
! These quantities are biologically irrelevant
! Michaelis-Menten constant for CO2 (Pa)
! Michaelis-Menten constant for O2 (Pa)
! CO2 compensation point (Pa)
! leaf boundary layer conductance of h20
! constrained vapor pressure
call GetCanopyGasParameters(bc_in(s)%forc_pbot, & ! in
bc_in(s)%oair_pa(ifp), & ! in
bc_in(s)%t_veg_pa(ifp), & ! in
bc_in(s)%tgcm_pa(ifp), & ! in
bc_in(s)%eair_pa(ifp), & ! in
bc_in(s)%esat_tv_pa(ifp), & ! in
bc_in(s)%rb_pa(ifp), & ! in
mm_kco2, & ! out
mm_ko2, & ! out
co2_cpoint, & ! out
cf, & ! out
gb_mol, & ! out
ceair) ! out
! ------------------------------------------------------------------------
! Part VI: Loop over all leaf layers.
! The concept of leaf layers is a result of the radiative transfer scheme.
! A leaf layer has uniform radiation environment. Leaf layers are a group
! of vegetation surfaces (stems and leaves) which inhabit the same
! canopy-layer "CL", have the same functional type "ft" and within those
! two partitions are further partitioned into vertical layers where
! downwelling radiation attenuates in order.
! In this phase we loop over the leaf layers and calculate the
! photosynthesis and respiration of the layer (since all biophysical
! properties are homogeneous). After this step, we can loop through
! our cohort list, associate each cohort with its list of leaf-layers
! and transfer these quantities to the cohort.
! With plant hydraulics, we must realize that photosynthesis and
! respiration will be different for leaves of each cohort in the leaf
! layers, as they will have there own hydraulic limitations.
! NOTE: Only need to flush mask on the number of used pfts, not the whole
! scratch space.
! ------------------------------------------------------------------------
rate_mask_z(:,1:numpft,:) = .false.
if_any_cohorts: if(currentPatch%countcohorts > 0.0)then
currentCohort => currentPatch%tallest
do_cohort_drive: do while (associated(currentCohort)) ! Cohort loop
! Identify the canopy layer (cl), functional type (ft)
! and the leaf layer (IV) for this cohort
ft = currentCohort%pft
cl = currentCohort%canopy_layer
! Calculate the cohort specific elai profile
! And the top and bottom edges of the veg area index
! of each layer bin are. Note, if the layers
! sink below the ground snow line, then the effective
! LAI and SAI start to shrink to zero, as well as
! the difference between vaitop and vaibot.
if(currentCohort%treesai>0._r8)then
do iv = 1,currentCohort%nv
call VegAreaLayer(currentCohort%treelai, &
currentCohort%treesai, &
currentCohort%height, &
iv, &
currentCohort%nv, &
currentCohort%pft, &
sites(s)%snow_depth, &
cohort_vaitop(iv), &
cohort_vaibot(iv), &
cohort_layer_elai(iv), &
cohort_layer_esai(iv))
end do
cohort_elai = sum(cohort_layer_elai(1:currentCohort%nv))
cohort_esai = sum(cohort_layer_esai(1:currentCohort%nv))
else
cohort_layer_elai(:) = 0._r8
cohort_layer_esai(:) = 0._r8
cohort_vaitop(:) = 0._r8
cohort_vaibot(:) = 0._r8
cohort_elai = 0._r8
cohort_esai = 0._r8
end if
! MLO. Assuming target to be related to leaf biomass when leaves are fully
! flushed. But unsure whether this call is correct or not, shouldn't we get
! the target value directly from the bstore_allom?
call bleaf(currentCohort%dbh,currentCohort%pft,&
currentCohort%crowndamage,currentCohort%canopy_trim,1.0_r8,store_c_target)
! call bstore_allom(currentCohort%dbh,currentCohort%pft, &
! currentCohort%canopy_trim,store_c_target)
call storage_fraction_of_target(store_c_target, &
currentCohort%prt%GetState(store_organ, carbon12_element), &
frac)
call lowstorage_maintresp_reduction(frac,currentCohort%pft, &
maintresp_reduction_factor)
! are there any leaves of this pft in this layer?
canopy_mask_if: if(currentPatch%canopy_mask(cl,ft) == 1)then
! Loop over leaf-layers
leaf_layer_loop : do iv = 1,currentCohort%nv
! ------------------------------------------------------------
! If we are doing plant hydro-dynamics (or any run-type
! where cohorts may generate different photosynthetic rates
! of other cohorts in the same canopy-pft-layer combo),
! we re-calculate the leaf biophysical rates for the
! cohort-layer combo of interest.
! but in the vanilla case, we only re-calculate if it has
! not been done yet.
! Other cases where we need to solve for every cohort
! in every leaf layer: nutrient dynamic mode, multiple leaf
! age classes
! ------------------------------------------------------------
rate_mask_if: if ( .not.rate_mask_z(iv,ft,cl) .or. &
(hlm_use_planthydro.eq.itrue) .or. &
(radiation_model .eq. twostr_solver ) .or. &
(nleafage > 1) .or. &
(hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then
if (hlm_use_planthydro.eq.itrue ) then
stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentCohort%co_hydr%btran )
btran_eff = currentCohort%co_hydr%btran
! dinc_vai(:) is the total vegetation area index of each "leaf" layer
! we convert to the leaf only portion of the increment
! ------------------------------------------------------
leaf_inc = dinc_vai(iv) * &
currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
! Now calculate the cumulative top-down lai of the current layer's midpoint
lai_canopy_above = sum(currentPatch%canopy_layer_tlai(1:cl-1))
lai_layers_above = (dlower_vai(iv) - dinc_vai(iv)) * &
currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
lai_current = min(leaf_inc, currentCohort%treelai - lai_layers_above)
cumulative_lai = lai_canopy_above + lai_layers_above + 0.5*lai_current
leaf_psi = currentCohort%co_hydr%psi_ag(1)
else
stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentPatch%btran_ft(ft) )
btran_eff = currentPatch%btran_ft(ft)
! For consistency sake, we use total LAI here, and not exposed
! if the plant is under-snow, it will be effectively dormant for
! the purposes of nscaler
cumulative_lai = sum(currentPatch%canopy_layer_tlai(1:cl-1)) + &
sum(currentPatch%tlai_profile(cl,ft,1:iv-1)) + &
0.5*currentPatch%tlai_profile(cl,ft,iv)
leaf_psi = fates_unset_r8
end if
if(do_fates_salinity)then
btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft)
endif
! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used
! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al
! (2010) Biogeosciences, 7, 1833-1859
kn = decay_coeff_vcmax(currentCohort%vcmax25top, &
prt_params%leafn_vert_scaler_coeff1(ft), &
prt_params%leafn_vert_scaler_coeff2(ft))
! Scale for leaf nitrogen profile
nscaler = exp(-kn * cumulative_lai)
! Leaf maintenance respiration to match the base rate used in CN
! but with the new temperature functions for C3 and C4 plants.
! CN respiration has units: g C / g N [leaf] / s. This needs to be
! converted from g C / g N [leaf] / s to umol CO2 / m**2 [leaf] / s
! Then scale this value at the top of the canopy for canopy depth
! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf)
select case(hlm_parteh_mode)
case (prt_carbon_allom_hyp)
lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft)
case (prt_cnp_flex_allom_hyp)
leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element)
if( (leaf_c*slatop(ft)) > nearzero) then
leaf_n = currentCohort%prt%GetState(leaf_organ, nitrogen_element)
lnc_top = leaf_n / (slatop(ft) * leaf_c )
else
lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft)
end if
! If one wants to break coupling with dynamic N conentrations,
! use the stoichiometry parameter
! lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft)
end select
! Part VII: Calculate dark respiration (leaf maintenance) for this layer
select case (maintresp_leaf_model)
case (lmrmodel_ryan_1991)
call LeafLayerMaintenanceRespiration_Ryan_1991( lnc_top, & ! in
nscaler, & ! in
ft, & ! in
bc_in(s)%t_veg_pa(ifp), & ! in
lmr_z(iv,ft,cl)) ! out
case (lmrmodel_atkin_etal_2017)
call LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & ! in
cumulative_lai, & ! in
currentCohort%vcmax25top, & ! in
ft, & ! in
bc_in(s)%t_veg_pa(ifp), & ! in
currentPatch%tveg_lpa%GetMean(), & ! in
lmr_z(iv,ft,cl)) ! out
case default
write (fates_log(),*)'error, incorrect leaf respiration model specified'
call endrun(msg=errMsg(sourcefile, __LINE__))
end select
! Pre-process PAR absorbed per unit leaf area for different schemes
! par_per_sunla = [W absorbed beam+diffuse radiation / m2 of sunlit leaves]
! par_per_shala = [W absorbed diffuse radiation / m2 of shaded leaves]
! fsun = [m2 of sunlit leaves / m2 of total leaves]
! laisun: m2 of exposed leaf, per m2 of crown. If this is the lowest layer
! for the pft/canopy group, than the m2 per crown is probably not
! as large as the layer above.
! ------------------------------------------------------------------
if_radsolver: if(radiation_model.eq.norman_solver) then
laisun = currentPatch%ed_laisun_z(cl,ft,iv)
laisha = currentPatch%ed_laisha_z(cl,ft,iv)
par_per_sunla = currentPatch%ed_parsun_z(cl,ft,iv)
par_per_shala = currentPatch%ed_parsha_z(cl,ft,iv)
canopy_area = currentPatch%canopy_area_profile(cl,ft,iv)
fsun = currentPatch%f_sun(cl,ft,iv)
else ! Two-stream
if(cohort_layer_elai(iv) > nearzero .and. currentPatch%solar_zenith_flag) then
call FatesGetCohortAbsRad(currentPatch, currentCohort, ipar, &
cohort_vaitop(iv), cohort_vaibot(iv), cohort_elai, cohort_esai, &
rb_abs, rd_abs, rb_abs_leaf, rd_abs_leaf, fsun)
! rd_abs_leaf: Watts of diffuse light absorbed by leaves over this
! depth interval and ground footprint (m2)
! rd_abs_leaf*fsun Watts of diffuse light absorbed by sunlit leaves
! over this depth interval and ground footprint (m2)
! rb_abs_leaf Watts of beam absorbed by sunlit leaves over this
! depth interval and ground footprint (m2)
! cohort_layer_elai*fsun Leaf area in sunlight within this interval and ground footprint
! cohort_layer_elai*(1-fsun) Leaf area in shade within this interval and ground footprint
laisun = (fsun*cohort_layer_elai(iv))
laisha = ((1._r8 - fsun)*cohort_layer_elai(iv))
if(fsun>nearzero) then
par_per_sunla = (rd_abs_leaf*fsun + rb_abs_leaf)! / laisun
else
par_per_sunla = 0._r8
end if
par_per_shala = rd_abs_leaf*(1._r8-fsun) !/ laisha
canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv)
else
par_per_sunla = 0._r8
par_per_shala = 0._r8
laisun = 0.5_r8*cohort_layer_elai(iv)
laisha = 0.5_r8*cohort_layer_elai(iv)
canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv)
fsun = 0.5_r8 !avoid div0, should have no impact
end if
end if if_radsolver
! Part VII: Calculate (1) maximum rate of carboxylation (vcmax),
! (2) maximum electron transport rate, (3) triose phosphate
! utilization rate and (4) the initial slope of CO2 response curve
! (C4 plants). Earlier we calculated their base rates as dictated
! by their plant functional type and some simple scaling rules for
! nitrogen limitation baesd on canopy position (not prognostic).
! These rates are the specific rates used in the actual photosynthesis
! calculations that take localized environmental effects (temperature)
! into consideration.
call LeafLayerBiophysicalRates(par_per_sunla, & ! in
ft, & ! in
currentCohort%vcmax25top, & ! in
currentCohort%jmax25top, & ! in
currentCohort%kp25top, & ! in
nscaler, & ! in
bc_in(s)%t_veg_pa(ifp), & ! in
bc_in(s)%dayl_factor_pa(ifp), & ! in
currentPatch%tveg_lpa%GetMean(), & ! in
currentPatch%tveg_longterm%GetMean(),& ! in
btran_eff, & ! in
vcmax_z, & ! out
jmax_z, & ! out
kp_z ) ! out
! Part IX: This call calculates the actual photosynthesis for the
! leaf layer, as well as the stomatal resistance and the net assimilated carbon.
call LeafLayerPhotosynthesis(fsun, & ! in
par_per_sunla, & ! in
par_per_shala, & ! in
laisun, & ! in
laisha, & ! in
canopy_area, & ! in
ft, & ! in
vcmax_z, & ! in
jmax_z, & ! in
kp_z, & ! in
bc_in(s)%t_veg_pa(ifp), & ! in
bc_in(s)%esat_tv_pa(ifp), & ! in
bc_in(s)%forc_pbot, & ! in
bc_in(s)%cair_pa(ifp), & ! in
bc_in(s)%oair_pa(ifp), & ! in
btran_eff, & ! in
stomatal_intercept_btran, & ! in
cf, & ! in
gb_mol, & ! in
ceair, & ! in
mm_kco2, & ! in
mm_ko2, & ! in
co2_cpoint, & ! in
lmr_z(iv,ft,cl), & ! in
leaf_psi, & ! in
bc_in(s)%rb_pa(ifp), & ! in
currentPatch%psn_z(cl,ft,iv), & ! out
rs_z(iv,ft,cl), & ! out
anet_av_z(iv,ft,cl), & ! out
c13disc_z(cl,ft,iv)) ! out
rate_mask_z(iv,ft,cl) = .true.
end if rate_mask_if
end do leaf_layer_loop
! Zero cohort flux accumulators.
currentCohort%npp_tstep = 0.0_r8
currentCohort%resp_tstep = 0.0_r8
currentCohort%gpp_tstep = 0.0_r8
currentCohort%rdark = 0.0_r8
currentCohort%resp_m = 0.0_r8
currentCohort%ts_net_uptake = 0.0_r8
currentCohort%c13disc_clm = 0.0_r8
! ---------------------------------------------------------------
! Part VII: Transfer leaf flux rates (like maintenance respiration,
! carbon assimilation and conductance) that are defined by the
! leaf layer (which is area independent, ie /m2) onto each cohort
! (where the rates become per cohort, ie /individual). Most likely
! a sum over layers.
! ---------------------------------------------------------------
nv = currentCohort%nv
! Temporary bypass to preserve B4B behavior
if(radiation_model.eq.norman_solver) then
call ScaleLeafLayerFluxToCohort(nv, & !in
currentPatch%psn_z(cl,ft,1:nv), & !in
lmr_z(1:nv,ft,cl), & !in
rs_z(1:nv,ft,cl), & !in
currentPatch%elai_profile(cl,ft,1:nv), & !in
c13disc_z(cl, ft, 1:nv), & !in
currentCohort%c_area, & !in
currentCohort%n, & !in
bc_in(s)%rb_pa(ifp), & !in
maintresp_reduction_factor, & !in
currentCohort%g_sb_laweight, & !out
currentCohort%gpp_tstep, & !out
currentCohort%rdark, & !out
currentCohort%c13disc_clm, & !out
cohort_eleaf_area) !out
else
call ScaleLeafLayerFluxToCohort(nv, & !in
currentPatch%psn_z(cl,ft,1:nv), & !in
lmr_z(1:nv,ft,cl), & !in
rs_z(1:nv,ft,cl), & !in
cohort_layer_elai(1:nv), & !in
c13disc_z(cl, ft, 1:nv), & !in
currentCohort%c_area, & !in
currentCohort%n, & !in
bc_in(s)%rb_pa(ifp), & !in
maintresp_reduction_factor, & !in
currentCohort%g_sb_laweight, & !out
currentCohort%gpp_tstep, & !out
currentCohort%rdark, & !out
currentCohort%c13disc_clm, & !out
cohort_eleaf_area) !out
end if
! Net Uptake does not need to be scaled, just transfer directly
currentCohort%ts_net_uptake(1:nv) = anet_av_z(1:nv,ft,cl) * umolC_to_kgC
else
! In this case, the cohort had no leaves,
! so no productivity,conductance, transpiration uptake
! or dark respiration
cohort_eleaf_area = 0.0_r8
currentCohort%gpp_tstep = 0.0_r8
currentCohort%rdark = 0.0_r8
currentCohort%g_sb_laweight = 0.0_r8
currentCohort%ts_net_uptake(:) = 0.0_r8
end if canopy_mask_if
! ------------------------------------------------------------------
! Part VIII: Calculate maintenance respiration in the sapwood and
! fine root pools.
! ------------------------------------------------------------------
! Calculate the amount of nitrogen in the above and below ground
! stem and root pools, used for maint resp
! We are using the fine-root C:N ratio as an approximation for
! the sapwood pools.
! Units are in (kgN/plant)
! ------------------------------------------------------------------
sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element)
fnrt_c = currentCohort%prt%GetState(fnrt_organ, carbon12_element)
if (hlm_use_tree_damage .eq. itrue) then
! Crown damage currenly only reduces the aboveground portion of
! sapwood. Therefore we calculate the aboveground and the belowground portion
! sapwood for use in stem respiration.
call GetCrownReduction(currentCohort%crowndamage, crown_reduction)
else
crown_reduction = 0.0_r8
end if
! If crown reduction is zero, undamaged sapwood target will equal sapwood carbon
agb_frac = prt_params%allom_agb_frac(currentCohort%pft)
branch_frac = param_derived%branch_frac(currentCohort%pft)
sapw_c_undamaged = sapw_c / (1.0_r8 - (agb_frac * branch_frac * crown_reduction))
! Undamaged below ground portion
sapw_c_bgw = sapw_c_undamaged * (1.0_r8 - agb_frac)
! Damaged aboveground portion
sapw_c_agw = sapw_c - sapw_c_bgw
select case(hlm_parteh_mode)
case (prt_carbon_allom_hyp)
live_stem_n = sapw_c_agw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ))
live_croot_n = sapw_c_bgw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ))
fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ))
case(prt_cnp_flex_allom_hyp)
live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * &
currentCohort%prt%GetState(sapw_organ, nitrogen_element)
live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * &
currentCohort%prt%GetState(sapw_organ, nitrogen_element)
fnrt_n = currentCohort%prt%GetState(fnrt_organ, nitrogen_element)
if (hlm_use_tree_damage .eq. itrue) then
sapw_n = currentCohort%prt%GetState(sapw_organ, nitrogen_element)
sapw_n_undamaged = sapw_n / &
(1.0_r8 - (agb_frac * branch_frac * crown_reduction))
sapw_n_bgw = sapw_n_undamaged * (1.0_r8 - agb_frac)
sapw_n_agw = sapw_n - sapw_n_bgw
live_croot_n = sapw_n_bgw
live_stem_n = sapw_n_agw
end if
! If one wants to break coupling with dynamic N conentrations,
! use the stoichiometry parameter
!
! live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * &
! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ))
! live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * &
! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ))
! fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ))
case default
end select
!------------------------------------------------------------------------------
! Calculate Whole Plant Respiration
! (this doesn't really need to be in this iteration at all, surely?)
! Response: (RGK 12-2016): I think the positioning of these calls is
! appropriate as of now. Maintenance calculations in sapwood and roots
! vary by cohort and with changing temperature at the minimum, and there are
! no sub-pools chopping up those pools any finer that need to be dealt with.
!------------------------------------------------------------------------------
! Live stem MR (kgC/plant/s) (above ground sapwood)
! ------------------------------------------------------------------
if ( int(woody(ft)) == itrue) then
tcwood = q10_mr**((bc_in(s)%t_veg_pa(ifp)-tfrz - 20.0_r8)/10.0_r8)
! kgC/s = kgN * kgC/kgN/s
currentCohort%livestem_mr = live_stem_n * maintresp_nonleaf_baserate * tcwood * maintresp_reduction_factor
else
currentCohort%livestem_mr = 0._r8
end if
! Fine Root MR (kgC/plant/s)
! and calculate the N fixation rate as a function of the fixation-specific root respiration
! for now use dev_arbitrary_pft as scaling term between 0 and 1 as additional increment of root respiration used for N fixation
! ------------------------------------------------------------------
currentCohort%froot_mr = 0._r8
currentCohort%sym_nfix_tstep = 0._r8
! n_fixation is integrated over the course of the day
! this variable is zeroed at the end of the FATES dynamics sequence
do j = 1,bc_in(s)%nlevsoil
tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8)
fnrt_mr_layer = fnrt_n * maintresp_nonleaf_baserate * tcsoi * rootfr_ft(ft,j) * maintresp_reduction_factor
! calculate the cost of carbon for N fixation in each soil layer and calculate N fixation rate based on that [kgC / kgN]
call RootLayerNFixation(bc_in(s)%t_soisno_sl(j),ft,dtime,fnrt_mr_layer,fnrt_mr_nfix_layer,nfix_layer)
currentCohort%froot_mr = currentCohort%froot_mr + fnrt_mr_nfix_layer + fnrt_mr_layer
currentCohort%sym_nfix_tstep = currentCohort%sym_nfix_tstep + nfix_layer
enddo
! Coarse Root MR (kgC/plant/s) (below ground sapwood)
! ------------------------------------------------------------------
if ( int(woody(ft)) == itrue) then
currentCohort%livecroot_mr = 0._r8
do j = 1,bc_in(s)%nlevsoil
! Soil temperature used to adjust base rate of MR
tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8)
currentCohort%livecroot_mr = currentCohort%livecroot_mr + &
live_croot_n * maintresp_nonleaf_baserate * tcsoi * &
rootfr_ft(ft,j) * maintresp_reduction_factor
enddo
else
currentCohort%livecroot_mr = 0._r8
end if
! ------------------------------------------------------------------
! Part IX: Perform some unit conversions (rate to integrated) and
! calcualate some fluxes that are sums and nets of the base fluxes
! ------------------------------------------------------------------
if ( debug ) write(fates_log(),*) 'EDPhoto 904 ', currentCohort%resp_m
if ( debug ) write(fates_log(),*) 'EDPhoto 905 ', currentCohort%rdark
if ( debug ) write(fates_log(),*) 'EDPhoto 906 ', currentCohort%livestem_mr
if ( debug ) write(fates_log(),*) 'EDPhoto 907 ', currentCohort%livecroot_mr
if ( debug ) write(fates_log(),*) 'EDPhoto 908 ', currentCohort%froot_mr
! add on whole plant respiration values in kgC/indiv/s-1
currentCohort%resp_m = currentCohort%livestem_mr + &
currentCohort%livecroot_mr + &
currentCohort%froot_mr
! no drought response right now.. something like:
! resp_m = resp_m * (1.0_r8 - currentPatch%btran_ft(currentCohort%pft) * &
! EDPftvarcon_inst%resp_drought_response(ft))
currentCohort%resp_m = currentCohort%resp_m + currentCohort%rdark
! save as a diagnostic the un-throttled maintenance respiration to be able to know how strong this is
currentCohort%resp_m_unreduced = currentCohort%resp_m / maintresp_reduction_factor
! convert from kgC/indiv/s to kgC/indiv/timestep
currentCohort%resp_m = currentCohort%resp_m * dtime
currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime
currentCohort%ts_net_uptake = currentCohort%ts_net_uptake * dtime
if ( debug ) write(fates_log(),*) 'EDPhoto 911 ', currentCohort%gpp_tstep
if ( debug ) write(fates_log(),*) 'EDPhoto 912 ', currentCohort%resp_tstep
if ( debug ) write(fates_log(),*) 'EDPhoto 913 ', currentCohort%resp_m
currentCohort%resp_g_tstep = prt_params%grperc(ft) * &
(max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m))