Bug in two-stream diagnostic arrays for normalized radiation profiles sets seedling-layer PAR to ~0

[r-ward]

Describe the issue

While adding a history variable that tracks seedling-layer light by patch age, I noticed patches > 1yr of age receive essentially 0 PAR. This matters when using the Tree Recruitment Scheme (TRS) because light mediates PFT-specific germination and seedling survival.

I think I've found the bug in radiation/FatesRadiationDriveMod.F90.

The FatesNormalizedCanopyRadiation subroutine fills diagnostic arrays of the fraction of incoming direct/diffuse radiation reaching each canopy layer:

fates/radiation/FatesRadiationDriveMod.F90

Line 61 in 121723f

subroutine FatesNormalizedCanopyRadiation(sites, bc_in, bc_out )

Direct and diffuse top-of-canopy boundary conditions need to be set to 1 before filling arrays (nrmlzd_parprof_pft_dir_z/dif_z) because the boundary is used by the radiation functions (GetRb/GetRdDn/GetRdUp) to scale per-layer fractions.

The debug loop does this:

fates/radiation/FatesRadiationDriveMod.F90

Lines 185 to 191 in 121723f

	if(debug) then
	currentPatch%twostr%band(ib)%Rbeam_atm = 1._r8
	currentPatch%twostr%band(ib)%Rdiff_atm = 1._r8
	call CheckPatchRadiationBalance(currentPatch, sites(s)%snow_depth, &
	ib, bc_out(s)%fabd_parb(ifp,ib),bc_out(s)%fabi_parb(ifp,ib))
	currentPatch%twostr%band(ib)%Rbeam_atm = fates_unset_r8
	currentPatch%twostr%band(ib)%Rdiff_atm = fates_unset_r8

But the main loop does not:

fates/radiation/FatesRadiationDriveMod.F90

Lines 202 to 220 in 121723f

	! Fill in the diagnostic arrays for normalized radiation profiles
	do_cl: do cl = 1,twostr%n_lyr
	do_icol: do icol = 1,twostr%n_col(cl)
	ft = twostr%scelg(cl,icol)%pft
	if_notair: if (ft>0) then
	area_frac = twostr%scelg(cl,icol)%area
	vai = twostr%scelg(cl,icol)%sai+twostr%scelg(cl,icol)%lai
	nv = GetNVegLayers(vai)
	do iv = 1, nv
	vai_top = dlower_vai(iv)
	currentPatch%nrmlzd_parprof_pft_dir_z(cl,ft,iv) = currentPatch%nrmlzd_parprof_pft_dir_z(cl,ft,iv) + &
	area_frac*twostr%GetRb(cl,icol,ivis,vai_top)
	currentPatch%nrmlzd_parprof_pft_dif_z(cl,ft,iv) = currentPatch%nrmlzd_parprof_pft_dif_z(cl,ft,iv) + &
	area_frac*twostr%GetRdDn(cl,icol,ivis,vai_top) + &
	area_frac*twostr%GetRdUp(cl,icol,ivis,vai_top)
	end do
	end if if_notair
	end do do_icol
	end do do_cl

Because GetRb/GetRdDn/GetRdUp scale each layer's normalized value by the top of canopy boundary, if unset the diagnostic arrays are scaled by fates_unset_r8 (=-1e36), and later truncated to 0.

These diagnostic arrays are only used in:

• /FatesInterfaceMod.F90 to calculate PAR at the seedling layer
• /FatesHistoryInterfaceMod.F90 to populate the FATES_PARPROF_DIR/_DIF diagnostic history variables

Carbon/structure are unaffected because they use outputs from Solve():

fates/radiation/FatesRadiationDriveMod.F90

Line 170 in 121723f

call twostr%Solve(ib, & ! in

The fix is to set the top-of-canopy boundary to 1 and reset afterwards, as in the debug loop. I'm working on a PR that will fix this and add FATES-TRS history vars including FATES_SEEDLING_LAYER_PAR_AP.

Relevant log output

FATES tag

sci.1.87.5_api.41.0.0-36-g66cdbd41

Host land model tag

ctsm5.3.077

Machine

derecho

Other supported machine name

No response

Additional context

[rgknox]

Great find @r-ward , yeah, I think you are correct and the atmospheric flux needs to be initialized before calling those GetXX functions:

currentPatch%twostr%band(ib)%Rbeam_atm = 1._r8 
    currentPatch%twostr%band(ib)%Rdiff_atm = 1._r8 

@r-ward and I are planning to talk off-line today about this.