forked from ESCOMP/CTSM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCNDVLightMod.F90
232 lines (189 loc) · 9.98 KB
/
CNDVLightMod.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
module CNDVLightMod
!-----------------------------------------------------------------------
! !DESCRIPTION:
! Calculate light competition
! Update fpc for establishment routine
! Called once per year
!
! !USES:
use shr_kind_mod , only: r8 => shr_kind_r8
use shr_const_mod , only : SHR_CONST_PI
use decompMod , only : bounds_type
use pftconMod , only : pftcon
use CNDVType , only : dgv_ecophyscon, dgvs_type
use CNVegCarbonStateType , only : cnveg_carbonstate_type
use PatchType , only : patch
!
! !PUBLIC TYPES:
implicit none
!
! !PUBLIC MEMBER FUNCTIONS:
public :: Light
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
subroutine Light(bounds, num_natvegp, filter_natvegp, &
cnveg_carbonstate_inst, dgvs_inst)
!
! !DESCRIPTION:
! Calculate light competition and update fpc for establishment routine
! Called once per year
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
integer , intent(in) :: num_natvegp ! number of naturally-vegetated patches in filter
integer , intent(in) :: filter_natvegp(:) ! patch filter for naturally-vegetated points
type(cnveg_carbonstate_type) , intent(in) :: cnveg_carbonstate_inst
type(dgvs_type) , intent(inout) :: dgvs_inst
!
! !LOCAL VARIABLES:
real(r8), parameter :: fpc_tree_max = 0.95_r8 !maximum total tree FPC
integer :: p,fp, g ! indices
real(r8) :: fpc_tree_total(bounds%begg:bounds%endg)
real(r8) :: fpc_inc_tree(bounds%begg:bounds%endg)
real(r8) :: fpc_inc(bounds%begp:bounds%endp) ! foliar projective cover increment (fraction)
real(r8) :: fpc_grass_total(bounds%begg:bounds%endg)
real(r8) :: fpc_shrub_total(bounds%begg:bounds%endg)
real(r8) :: fpc_grass_max(bounds%begg:bounds%endg)
real(r8) :: fpc_shrub_max(bounds%begg:bounds%endg)
integer :: numtrees(bounds%begg:bounds%endg)
real(r8) :: excess
real(r8) :: nind_kill
real(r8) :: lai_ind
real(r8) :: fpc_ind
real(r8) :: fpcgrid_old
real(r8) :: lm_ind ! leaf carbon (gC/individual)
real(r8) :: stemdiam ! stem diameter
real(r8) :: stocking ! #stems / ha (stocking density)
real(r8) :: taper ! ratio of height:radius_breast_height (tree allometry)
!-----------------------------------------------------------------------
associate( &
ivt => patch%itype , & ! Input: [integer (:) ] patch vegetation type
crownarea_max => dgv_ecophyscon%crownarea_max , & ! Input: [real(r8) (:) ] ecophys const - tree maximum crown a
reinickerp => dgv_ecophyscon%reinickerp , & ! Input: [real(r8) (:) ] ecophys const - parameter in allomet
allom1 => dgv_ecophyscon%allom1 , & ! Input: [real(r8) (:) ] ecophys const - parameter in allomet
dwood => pftcon%dwood , & ! Input: wood density (gC/m3)
slatop => pftcon%slatop , & ! Input: specific leaf area at top of canopy, projected area basis (m2/gC)
dsladlai => pftcon%dsladlai , & ! Input: dSLA/dLAI, projected area basis (m2/gC)
woody => pftcon%woody , & ! Input: woody patch or not
is_tree => pftcon%is_tree , & ! Input: tree patch or not
is_shrub => pftcon%is_shrub , & ! Input: shrub patch or not
deadstemc => cnveg_carbonstate_inst%deadstemc_patch , & ! Input: [real(r8) (:) ] (gC/m2) dead stem C
leafcmax => cnveg_carbonstate_inst%leafcmax_patch , & ! Input: [real(r8) (:) ] (gC/m2) leaf C storage
crownarea => dgvs_inst%crownarea_patch , & ! Output: [real(r8) (:) ] area that each individual tree takes up (m^2)
nind => dgvs_inst%nind_patch , & ! Output: [real(r8) (:) ] number of individuals
fpcgrid => dgvs_inst%fpcgrid_patch & ! Output: [real(r8) (:) ] foliar projective cover on gridcell (fraction)
)
taper = 200._r8 ! make a global constant; used in Establishment + ?
! Initialize gridcell-level metrics
do g = bounds%begg, bounds%endg
fpc_tree_total(g) = 0._r8
fpc_inc_tree(g) = 0._r8
fpc_grass_total(g) = 0._r8
fpc_shrub_total(g) = 0._r8
numtrees(g) = 0
end do
do fp = 1,num_natvegp
p = filter_natvegp(fp)
g = patch%gridcell(p)
! Update LAI and FPC as in the last lines of DGVMAllocation
if (woody(ivt(p))==1._r8) then
if (fpcgrid(p) > 0._r8 .and. nind(p) > 0._r8) then
stocking = nind(p)/fpcgrid(p) !#ind/m2 nat veg area -> #ind/m2 patch area
! stemdiam derived here from cn's formula for htop found in
! CNVegStructUpdate and cn's assumption stemdiam=2*htop/taper
! this derivation neglects upper htop limit enforced elsewhere
stemdiam = (24._r8 * deadstemc(p) / (SHR_CONST_PI * stocking * dwood(ivt(p)) * taper))**(1._r8/3._r8)
else
stemdiam = 0._r8
end if
crownarea(p) = min(crownarea_max(ivt(p)), allom1(ivt(p))*stemdiam**reinickerp(ivt(p))) ! Eqn D (from Establishment)
!else ! crownarea is 1 and does not need updating
end if
if (crownarea(p) > 0._r8 .and. nind(p) > 0._r8) then
lm_ind = leafcmax(p) * fpcgrid(p) / nind(p)
if (dsladlai(ivt(p)) > 0._r8) then
lai_ind = max(0.001_r8,((exp(lm_ind*dsladlai(ivt(p)) + log(slatop(ivt(p)))) - &
slatop(ivt(p)))/dsladlai(ivt(p))) / crownarea(p))
else
lai_ind = lm_ind * slatop(ivt(p)) / crownarea(p)
end if
else
lai_ind = 0._r8
end if
fpc_ind = 1._r8 - exp(-0.5_r8*lai_ind)
fpcgrid_old = fpcgrid(p)
fpcgrid(p) = crownarea(p) * nind(p) * fpc_ind
fpc_inc(p) = max(0._r8, fpcgrid(p) - fpcgrid_old)
if (woody(ivt(p)) == 1._r8) then
if (is_tree(ivt(p))) then
numtrees(g) = numtrees(g) + 1
fpc_tree_total(g) = fpc_tree_total(g) + fpcgrid(p)
fpc_inc_tree(g) = fpc_inc_tree(g) + fpc_inc(p)
else if (is_shrub(ivt(p))) then
fpc_shrub_total(g) = fpc_shrub_total(g) + fpcgrid(p)
end if
else ! if grass
fpc_grass_total(g) = fpc_grass_total(g) + fpcgrid(p)
end if
end do
do g = bounds%begg, bounds%endg
fpc_grass_max(g) = 1._r8 - min(fpc_tree_total(g), fpc_tree_max)
fpc_shrub_max(g) = max(0._r8, fpc_grass_max(g) - fpc_grass_total(g))
end do
! The gridcell level metrics are now in place; continue...
! slevis replaced the previous code that updated pfpcgrid
! with a simpler way of doing so:
! fpcgrid(p) = fpcgrid(p) - excess
! Later we may wish to update this subroutine
! according to Strassmann's recommendations (see relevant pdf)
do fp = 1,num_natvegp
p = filter_natvegp(fp)
g = patch%gridcell(p)
! light competition
if (woody(ivt(p))==1._r8 .and. is_tree(ivt(p))) then
if (fpc_tree_total(g) > fpc_tree_max) then
if (fpc_inc_tree(g) > 0._r8) then
excess = (fpc_tree_total(g) - fpc_tree_max) * &
fpc_inc(p) / fpc_inc_tree(g)
else
excess = (fpc_tree_total(g) - fpc_tree_max) / &
real(numtrees(g))
end if
! Reduce individual density (and thereby gridcell-level biomass)
! so that total tree FPC reduced to 'fpc_tree_max'
if (fpcgrid(p) > 0._r8) then
nind_kill = nind(p) * excess / fpcgrid(p)
nind(p) = max(0._r8, nind(p) - nind_kill)
fpcgrid(p) = max(0._r8, fpcgrid(p) - excess)
else
nind(p) = 0._r8
fpcgrid(p) = 0._r8
end if
! Transfer lost biomass to litter
end if ! if tree cover exceeds max allowed
else if (woody(ivt(p))==0._r8) then ! grass
if (fpc_grass_total(g) > fpc_grass_max(g)) then
! grass competes with itself if total fpc exceeds 1
excess = (fpc_grass_total(g) - fpc_grass_max(g)) * fpcgrid(p) / fpc_grass_total(g)
fpcgrid(p) = max(0._r8, fpcgrid(p) - excess)
end if
else if (woody(ivt(p))==1._r8 .and. is_shrub(ivt(p))) then ! shrub
if (fpc_shrub_total(g) > fpc_shrub_max(g)) then
excess = 1._r8 - fpc_shrub_max(g) / fpc_shrub_total(g)
! Reduce individual density (and thereby gridcell-level biomass)
! so that total shrub FPC reduced to fpc_shrub_max(g)
if (fpcgrid(p) > 0._r8) then
nind_kill = nind(p) * excess / fpcgrid(p)
nind(p) = max(0._r8, nind(p) - nind_kill)
fpcgrid(p) = max(0._r8, fpcgrid(p) - excess)
else
nind(p) = 0._r8
fpcgrid(p) = 0._r8
end if
end if
end if ! end of if-tree
end do
end associate
end subroutine Light
end module CNDVLightMod