Skip to content

Commit 0090297

Browse files
rgknoxrgknox
authored
Merge pull request #357 from rgknox/rgknox-profile-cleaning
refactors and some numerical protections on lai profile
2 parents 1b93458 + f5d5a7c commit 0090297

File tree

6 files changed

+500
-402
lines changed

6 files changed

+500
-402
lines changed

‎biogeochem/EDCanopyStructureMod.F90‎

Lines changed: 451 additions & 359 deletions
Large diffs are not rendered by default.

‎biogeochem/EDPatchDynamicsMod.F90‎

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1219,9 +1219,8 @@ subroutine zero_patch(cp_p)
12191219
currentPatch%age = nan
12201220
currentPatch%age_class = 1
12211221
currentPatch%area = nan
1222-
currentPatch%canopy_layer_lai(:) = nan
1222+
currentPatch%canopy_layer_tai(:) = nan
12231223
currentPatch%total_canopy_area = nan
1224-
currentPatch%canopy_area = nan
12251224
currentPatch%bare_frac_area = nan
12261225

12271226
currentPatch%tlai_profile(:,:,:) = nan
@@ -1248,10 +1247,9 @@ subroutine zero_patch(cp_p)
12481247
currentPatch%fabd(:) = nan ! fraction of incoming direct radiation that is absorbed by the canopy
12491248
currentPatch%fabi(:) = nan ! fraction of incoming diffuse radiation that is absorbed by the canopy
12501249

1251-
currentPatch%present(:,:) = 999 ! is there any of this pft in this layer?
1250+
currentPatch%canopy_mask(:,:) = 999 ! is there any of this pft in this layer?
12521251
currentPatch%nrad(:,:) = 999 ! number of exposed leaf layers for each canopy layer and pft
12531252
currentPatch%ncan(:,:) = 999 ! number of total leaf layers for each canopy layer and pft
1254-
currentPatch%lai = nan ! leaf area index of patch
12551253
currentPatch%pft_agb_profile(:,:) = nan
12561254

12571255
! DISTURBANCE
@@ -1302,7 +1300,7 @@ subroutine zero_patch(cp_p)
13021300
currentPatch%burnt_frac_litter(:) = 0.0_r8
13031301
currentPatch%btran_ft(:) = 0.0_r8
13041302

1305-
currentPatch%canopy_layer_lai(:) = 0.0_r8
1303+
currentPatch%canopy_layer_tai(:) = 0.0_r8
13061304

13071305
currentPatch%seeds_in(:) = 0.0_r8
13081306
currentPatch%seed_decay(:) = 0.0_r8
@@ -1714,13 +1712,13 @@ subroutine terminate_patches(cs_pnt)
17141712
! This is only really meant for very old patches.
17151713
if(associated(currentPatch%older) )then
17161714
write(fates_log(),*) 'fusing to older patch because this one is too small',&
1717-
currentPatch%area, currentPatch%lai, &
1718-
currentPatch%older%area,currentPatch%older%lai
1715+
currentPatch%area, &
1716+
currentPatch%older%area
17191717
call fuse_2_patches(currentPatch%older, currentPatch)
17201718
write(fates_log(),*) 'after fusion to older patch',currentPatch%area
17211719
else
17221720
write(fates_log(),*) 'fusing to younger patch because oldest one is too small',&
1723-
currentPatch%area, currentPatch%lai
1721+
currentPatch%area
17241722
tmpptr => currentPatch%younger
17251723
call fuse_2_patches(currentPatch, currentPatch%younger)
17261724
write(fates_log(),*) 'after fusion to younger patch'

‎biogeophys/EDSurfaceAlbedoMod.F90‎

Lines changed: 14 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -197,10 +197,10 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
197197
! Is this pft/canopy layer combination present in this patch?
198198
do L = 1,nclmax
199199
do ft = 1,numpft
200-
currentPatch%present(L,ft) = 0
200+
currentPatch%canopy_mask(L,ft) = 0
201201
do iv = 1, currentPatch%nrad(L,ft)
202202
if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then
203-
currentPatch%present(L,ft) = 1
203+
currentPatch%canopy_mask(L,ft) = 1
204204
!I think 'present' is only used here...
205205
endif
206206
end do !iv
@@ -260,7 +260,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
260260
weighted_dif_ratio(L,1:hlm_numSWb) = 0._r8
261261
!Each canopy layer (canopy, understorey) has multiple 'parallel' pft's
262262
do ft =1,numpft
263-
if (currentPatch%present(L,ft) == 1)then !only do calculation if there are the appropriate leaves.
263+
if (currentPatch%canopy_mask(L,ft) == 1)then !only do calculation if there are the appropriate leaves.
264264
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
265265
! Diffuse transmittance, tr_dif, do each layer with thickness elai_z.
266266
! Estimated do nine sky angles in increments of 10 degrees
@@ -392,7 +392,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
392392

393393
do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up.
394394
do ft = 1,numpft
395-
if (currentPatch%present(L,ft) == 1)then
395+
if (currentPatch%canopy_mask(L,ft) == 1)then
396396
!==============================================================================!
397397
! Iterative solution do scattering
398398
!==============================================================================!
@@ -440,7 +440,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
440440
dif_ratio(L,ft,1,ib) * ftweight(L,ft,1)
441441
!instance where the first layer ftweight is used a proxy for the whole column. FTWA
442442
end do!hlm_numSWb
443-
endif ! currentPatch%present
443+
endif ! currentPatch%canopy_mask
444444
end do!ft
445445
end do!L
446446

@@ -450,7 +450,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
450450
do L = 1, currentPatch%NCL_p !work down from the top of the canopy.
451451
weighted_dif_down(L) = 0._r8
452452
do ft = 1, numpft
453-
if (currentPatch%present(L,ft) == 1)then
453+
if (currentPatch%canopy_mask(L,ft) == 1)then
454454
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
455455
! First estimates do downward and upward diffuse flux
456456
!
@@ -506,7 +506,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
506506
do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom.
507507
weighted_dif_up(L) = 0._r8
508508
do ft = 1, numpft
509-
if (currentPatch%present(L,ft) == 1)then
509+
if (currentPatch%canopy_mask(L,ft) == 1)then
510510
!Bounce diffuse radiation off soil surface.
511511
iv = currentPatch%nrad(L,ft) + 1
512512
if (L==currentPatch%NCL_p)then !is this the bottom layer ?
@@ -562,7 +562,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
562562
do L = 1,currentPatch%NCL_p !working from the top down
563563
weighted_dif_down(L) = 0._r8
564564
do ft =1,numpft
565-
if (currentPatch%present(L,ft) == 1)then
565+
if (currentPatch%canopy_mask(L,ft) == 1)then
566566
! forward diffuse flux within the canopy and at soil, working forward through canopy
567567
! with Dif_up -from previous iteration-. Dif_dn(1) is the forward diffuse flux onto the canopy.
568568
! Note: down = forward flux onto next layer
@@ -618,7 +618,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
618618
do L = 1, currentPatch%NCL_p ! working from the top down.
619619
weighted_dif_up(L) = 0._r8
620620
do ft =1,numpft
621-
if (currentPatch%present(L,ft) == 1)then
621+
if (currentPatch%canopy_mask(L,ft) == 1)then
622622
! Upward diffuse flux at soil or from lower canopy (forward diffuse and unscattered direct beam)
623623
iv = currentPatch%nrad(L,ft) + 1
624624
if (L==currentPatch%NCL_p)then !In the bottom canopy layer, reflect off the soil
@@ -670,7 +670,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
670670
abs_dir_z(:,:) = 0._r8
671671
abs_dif_z(:,:) = 0._r8
672672
do ft =1,numpft
673-
if (currentPatch%present(L,ft) == 1)then
673+
if (currentPatch%canopy_mask(L,ft) == 1)then
674674
!==============================================================================!
675675
! Compute absorbed flux densities
676676
!==============================================================================!
@@ -799,7 +799,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
799799
currentPatch%tr_soil_dir(ib)* &
800800
(1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1))
801801
write(fates_log(),*) 'albedos',currentPatch%sabs_dir(ib) ,currentPatch%tr_soil_dir(ib), &
802-
(1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%lai
802+
(1.0_r8-bc_in(s)%albgr_dir_rb(ib))
803803

804804
do ft =1,3
805805
iv = currentPatch%nrad(1,ft) + 1
@@ -826,7 +826,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
826826
lai_reduction(:) = 0.0_r8
827827
do L = 1, currentPatch%NCL_p
828828
do ft =1,numpft
829-
if (currentPatch%present(L,ft) == 1)then
829+
if (currentPatch%canopy_mask(L,ft) == 1)then
830830
do iv = 1, currentPatch%nrad(L,ft)
831831
if (lai_change(L,ft,iv) > 0.0_r8)then
832832
lai_reduction(L) = max(lai_reduction(L),lai_change(L,ft,iv))
@@ -879,7 +879,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out )
879879
write(fates_log(),*) 'bc_in(s)%albgr_dif_rb(ib)',bc_in(s)%albgr_dif_rb(ib)
880880
write(fates_log(),*) 'rhol',rhol(1:numpft,:)
881881
write(fates_log(),*) 'ftw',sum(ftweight(1,1:numpft,1)),ftweight(1,1:numpft,1)
882-
write(fates_log(),*) 'present',currentPatch%present(1,1:numpft)
882+
write(fates_log(),*) 'present',currentPatch%canopy_mask(1,1:numpft)
883883
write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft,1)
884884

885885
bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error
@@ -1003,7 +1003,7 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out)
10031003

10041004
if(bc_out(s)%fsun_pa(ifp) > 1._r8)then
10051005
write(fates_log(),*) 'too much leaf area in profile', bc_out(s)%fsun_pa(ifp), &
1006-
cpatch%lai,sunlai,shalai
1006+
sunlai,shalai
10071007
endif
10081008

10091009
elai = calc_areaindex(cpatch,'elai')

‎biogeophys/FatesPlantRespPhotosynthMod.F90‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -242,7 +242,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
242242
! And then identify which layer/pft combinations have things in them.
243243
! Output:
244244
! currentPatch%ncan(:,:)
245-
! currentPatch%present(:,:)
245+
! currentPatch%canopy_mask(:,:)
246246
call UpdateCanopyNCanNRadPresent(currentPatch)
247247

248248

@@ -322,12 +322,12 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
322322

323323

324324
! are there any leaves of this pft in this layer?
325-
if(currentPatch%present(cl,ft) == 1)then
325+
if(currentPatch%canopy_mask(cl,ft) == 1)then
326326

327327
if(cl==NCL_p)then !are we in the top canopy layer or a shaded layer?
328328
laican = 0._r8
329329
else
330-
laican = sum(currentPatch%canopy_layer_lai(cl+1:NCL_p))
330+
laican = sum(currentPatch%canopy_layer_tai(cl+1:NCL_p))
331331
end if
332332

333333
! Loop over leaf-layers
@@ -484,7 +484,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
484484
currentCohort%gscan = 0.0_r8
485485
currentCohort%ts_net_uptake(:) = 0.0_r8
486486

487-
end if ! if(currentPatch%present(cl,ft) == 1)then
487+
end if ! if(currentPatch%canopy_mask(cl,ft) == 1)then
488488

489489

490490
! ------------------------------------------------------------------
@@ -1281,7 +1281,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch)
12811281
! ---------------------------------------------------------------------------------
12821282
! This subroutine calculates two patch level quanities:
12831283
! currentPatch%ncan and
1284-
! currentPatch%present
1284+
! currentPatch%canopy_mask
12851285
!
12861286
! currentPatch%ncan(:,:) is a two dimensional array that indicates
12871287
! the total number of leaf layers (including those that are not exposed to light)
@@ -1291,7 +1291,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch)
12911291
! the total number of EXPOSED leaf layers, but for all intents and purposes
12921292
! in the photosynthesis routine, this appears to be the same as %ncan...
12931293
!
1294-
! currentPatch%present(:,:) has the same dimensions, is binary, and
1294+
! currentPatch%canopy_mask(:,:) has the same dimensions, is binary, and
12951295
! indicates whether or not leaf layers are present (by evaluating the canopy area
12961296
! profile).
12971297
! ---------------------------------------------------------------------------------
@@ -1334,10 +1334,10 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch)
13341334
! Now loop through and identify which layer and pft combo has scattering elements
13351335
do cl = 1,nclmax
13361336
do ft = 1,numpft
1337-
currentPatch%present(cl,ft) = 0
1337+
currentPatch%canopy_mask(cl,ft) = 0
13381338
do iv = 1, currentPatch%nrad(cl,ft);
13391339
if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then
1340-
currentPatch%present(cl,ft) = 1
1340+
currentPatch%canopy_mask(cl,ft) = 1
13411341
end if
13421342
end do !iv
13431343
enddo !ft

‎main/EDTypesMod.F90‎

Lines changed: 18 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -296,22 +296,31 @@ module EDTypesMod
296296

297297
! LEAF ORGANIZATION
298298
real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2
299-
real(r8) :: canopy_layer_lai(nclmax) ! lai that is shading this canopy layer: m2/m2
299+
real(r8) :: canopy_layer_tai(nclmax) ! total area index of each canopy layer
300+
! used to determine attenuation of parameters during
301+
! photosynthesis m2 veg / m2 of canopy area (patch without bare ground)
300302
real(r8) :: total_canopy_area ! area that is covered by vegetation : m2
301303
real(r8) :: total_tree_area ! area that is covered by woody vegetation : m2
302-
real(r8) :: canopy_area ! area that is covered by vegetation : m2 (is this different to total_canopy_area?
303304
real(r8) :: bare_frac_area ! bare soil in this patch expressed as a fraction of the total soil surface.
304-
real(r8) :: lai ! leaf area index of patch
305305
real(r8) :: zstar ! height of smallest canopy tree -- only meaningful in "strict PPA" mode
306306

307-
real(r8) :: tlai_profile(nclmax,maxpft,nlevleaf) ! total leaf area in each canopy layer, pft, and leaf layer. m2/m2
308-
real(r8) :: elai_profile(nclmax,maxpft,nlevleaf) ! exposed leaf area in each canopy layer, pft, and leaf layer. m2/m2
309-
real(r8) :: tsai_profile(nclmax,maxpft,nlevleaf) ! total stem area in each canopy layer, pft, and leaf layer. m2/m2
310-
real(r8) :: esai_profile(nclmax,maxpft,nlevleaf) ! exposed stem area in each canopy layer, pft, and leaf layer. m2/m2
307+
308+
! UNITS for the ai profiles
309+
! [ m2 leaf / m2 contributing crown footprints]
310+
real(r8) :: tlai_profile(nclmax,maxpft,nlevleaf) ! total leaf area in each canopy layer, pft, and leaf layer.
311+
real(r8) :: elai_profile(nclmax,maxpft,nlevleaf) ! exposed leaf area in each canopy layer, pft, and leaf layer
312+
real(r8) :: tsai_profile(nclmax,maxpft,nlevleaf) ! total stem area in each canopy layer, pft, and leaf layer
313+
real(r8) :: esai_profile(nclmax,maxpft,nlevleaf) ! exposed stem area in each canopy layer, pft, and leaf layer
314+
311315
real(r8) :: layer_height_profile(nclmax,maxpft,nlevleaf)
312-
real(r8) :: canopy_area_profile(nclmax,maxpft,nlevleaf) ! fraction of canopy in each canopy
316+
real(r8) :: canopy_area_profile(nclmax,maxpft,nlevleaf) ! fraction of crown area per canopy area in each layer
317+
! they will sum to 1.0 in the fully closed canopy layers
318+
! but only in leaf-layers that contain contributions
319+
! from all cohorts that donate to canopy_area
320+
321+
313322
! layer, pft, and leaf layer:-
314-
integer :: present(nclmax,maxpft) ! is there any of this pft in this canopy layer?
323+
integer :: canopy_mask(nclmax,maxpft) ! is there any of this pft in this canopy layer?
315324
integer :: nrad(nclmax,maxpft) ! number of exposed leaf layers for each canopy layer and pft
316325
integer :: ncan(nclmax,maxpft) ! number of total leaf layers for each canopy layer and pft
317326

@@ -682,9 +691,7 @@ subroutine dump_patch(cpatch)
682691
write(fates_log(),*) 'pa%ncl_p = ',cpatch%ncl_p
683692
write(fates_log(),*) 'pa%total_canopy_area = ',cpatch%total_canopy_area
684693
write(fates_log(),*) 'pa%total_tree_area = ',cpatch%total_tree_area
685-
write(fates_log(),*) 'pa%canopy_area = ',cpatch%canopy_area
686694
write(fates_log(),*) 'pa%bare_frac_area = ',cpatch%bare_frac_area
687-
write(fates_log(),*) 'pa%lai = ',cpatch%lai
688695
write(fates_log(),*) 'pa%zstar = ',cpatch%zstar
689696
write(fates_log(),*) 'pa%disturbance_rate = ',cpatch%disturbance_rate
690697
write(fates_log(),*) '----------------------------------------'

0 commit comments

Comments
 (0)