Skip to content

Commit

Permalink
Merge pull request #1141 from rgknox/two-stream-clean-b4b-nohist
Browse files Browse the repository at this point in the history 
two-stream - b4b Norman and no history changes
  • Loading branch information
@rgknox
rgknox authored Jan 16, 2024
2 parents 698a8df + 7b2cdbe commit 7d518ac
Show file tree
Hide file tree
Showing 34 changed files with 6,529 additions and 2,241 deletions.
291 changes: 169 additions & 122 deletions biogeochem/EDCanopyStructureMod.F90
View file

Large diffs are not rendered by default.

9 changes: 5 additions & 4 deletions biogeochem/EDPatchDynamicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ module EDPatchDynamicsMod
use FatesConstantsMod , only : fates_tiny
use FatesConstantsMod , only : nocomp_bareground
use FatesInterfaceTypesMod , only : hlm_use_planthydro
use FatesInterfaceTypesMod , only : hlm_numSWb
use FatesInterfaceTypesMod , only : bc_in_type
use FatesInterfaceTypesMod , only : numpft
use FatesInterfaceTypesMod , only : hlm_stepsize
Expand Down Expand Up @@ -106,6 +105,7 @@ module EDPatchDynamicsMod
use FatesRunningMeanMod, only : ema_sdlng_mdd
use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par, ema_sdlng2sap_par
use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm
use FatesRadiationMemMod, only : num_swb

! CIME globals
use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
Expand Down Expand Up @@ -603,7 +603,7 @@ subroutine spawn_patches( currentSite, bc_in)
allocate(newPatch)

call newPatch%Create(age, site_areadis, i_landusechange_receiverpatchlabel, i_nocomp_pft, &
hlm_numSWb, numpft, currentSite%nlevsoil, hlm_current_tod, &
num_swb, numpft, currentSite%nlevsoil, hlm_current_tod, &
regeneration_model)

! Initialize the litter pools to zero, these
Expand Down Expand Up @@ -2773,8 +2773,9 @@ subroutine fuse_2_patches(csite, dp, rp)
rp%zstar = (dp%zstar*dp%area + rp%zstar*rp%area) * inv_sum_area
rp%c_stomata = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area
rp%c_lblayer = (dp%c_lblayer*dp%area + rp%c_lblayer*rp%area) * inv_sum_area
rp%radiation_error = (dp%radiation_error*dp%area + rp%radiation_error*rp%area) * inv_sum_area

rp%rad_error(1) = (dp%rad_error(1)*dp%area + rp%rad_error(1)*rp%area) * inv_sum_area
rp%rad_error(2) = (dp%rad_error(2)*dp%area + rp%rad_error(2)*rp%area) * inv_sum_area

rp%area = rp%area + dp%area !THIS MUST COME AT THE END!

!insert donor cohorts into recipient patch
Expand Down
116 changes: 112 additions & 4 deletions biogeochem/FatesAllometryMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,12 +92,13 @@ module FatesAllometryMod
use FatesConstantsMod, only : calloc_abs_error
use FatesConstantsMod, only : fates_unset_r8
use FatesConstantsMod, only : itrue
use FatesConstantsMod, only : nearzero
use shr_log_mod , only : errMsg => shr_log_errMsg
use FatesGlobals , only : fates_log
use FatesGlobals , only : endrun => fates_endrun
use FatesGlobals , only : FatesWarn,N2S,A2S,I2S
use EDParamsMod , only : nlevleaf, dinc_vai
use EDParamsMod , only : nclmax
use EDParamsMod , only : nlevleaf,dinc_vai,dlower_vai
use EDParamsMod , only : nclmax
use DamageMainMod , only : GetCrownReduction

implicit none
Expand Down Expand Up @@ -125,7 +126,8 @@ module FatesAllometryMod
public :: set_root_fraction ! Generic wrapper to calculate normalized
! root profiles
public :: leafc_from_treelai ! Calculate target leaf carbon for a given treelai for SP mode

public :: VegAreaLayer

logical , parameter :: verbose_logging = .false.
character(len=*), parameter :: sourcefile = __FILE__

Expand Down Expand Up @@ -2436,7 +2438,8 @@ real(r8) function decay_coeff_kn(pft,vcmax25top)
end function decay_coeff_kn

! =====================================================================================
subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl )

subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl )

! =========================================================================
! This subroutine estimates the diameter based on either the structural biomass
Expand Down Expand Up @@ -2586,6 +2589,111 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d
return
end subroutine ForceDBH

! =========================================================================

subroutine VegAreaLayer(tree_lai,tree_sai,tree_height,iv,nv,pft,snow_depth, &
vai_top,vai_bot, elai_layer,esai_layer,tlai_layer,tsai_layer)

! -----------------------------------------------------------------------------------
! This routine returns the exposed leaf and stem areas (m2 of leaf and stem) per m2 of
! ground inside the crown, for the leaf-layer specified.
! -----------------------------------------------------------------------------------

real(r8),intent(in) :: tree_lai ! the in-crown leaf area index for the plant
! [m2 leaf/m2 crown footprint]
real(r8),intent(in) :: tree_sai ! the in-crown stem area index for the plant
! [m2 stem/m2 crown footprint]
real(r8),intent(in) :: tree_height ! the height of the plant [m]
integer,intent(in) :: iv ! vegetation layer index
integer,intent(in) :: nv ! this plants total number of veg layers
integer,intent(in) :: pft ! plant functional type index
real(r8),intent(in) :: snow_depth ! the depth of snow on the ground [m]
real(r8),intent(out) :: vai_top
real(r8),intent(out) :: vai_bot ! the VAI of the bin top and bottom
real(r8),intent(out) :: elai_layer ! exposed leaf area index of the layer
real(r8),intent(out) :: esai_layer ! exposed stem area index of the layer
real(r8),optional,intent(out) :: tlai_layer ! total leaf area index of the layer
real(r8),optional,intent(out) :: tsai_layer ! total stem area index of the layer

! [m2 of leaf in bin / m2 crown footprint]
real(r8) :: tree_vai ! the in-crown veg area index for the plant
real(r8) :: fraction_exposed ! fraction of the veg media that is above snow
real(r8) :: layer_top_height ! Physical height of the layer top relative to ground [m]
real(r8) :: layer_bot_height ! Physical height of the layer bottom relative to ground [m]
real(r8) :: tlai,tsai ! temporary total area indices [m2/m2]
real(r8) :: fleaf ! fraction of biomass in layer that is leaf
real(r8) :: remainder ! old-method: remainder of biomass in last bin
integer, parameter :: layer_height_const_depth = 1 ! constant physical depth assumption
integer, parameter :: layer_height_const_lad = 2 ! constant leaf area depth assumption
integer, parameter :: layer_height_method = layer_height_const_depth

tree_vai = tree_lai + tree_sai

if_any_vai: if(tree_vai>0._r8)then

if(iv==0)then
vai_top = 0.0
vai_bot = tree_vai
else

if(iv>1)then
vai_top = dlower_vai(iv) - dinc_vai(iv)
else
vai_top = 0._r8
end if

if(iv<nv) then
vai_bot = dlower_vai(iv)
else
vai_bot = tree_vai
end if
end if

if(layer_height_method .eq. layer_height_const_depth)then
if(iv==0)then
layer_top_height = tree_height
layer_bot_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft))
else
layer_top_height = tree_height*(1._r8 - real(iv-1,r8)/real(nv,r8)*prt_params%crown_depth_frac(pft))
layer_bot_height = tree_height*(1._r8 - real(iv,r8)/real(nv,r8)*prt_params%crown_depth_frac(pft))
end if
else
layer_top_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft)*vai_top/tree_vai)
layer_bot_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft)*vai_bot/tree_vai)
end if

fraction_exposed = 1._r8 - max(0._r8,(min(1._r8, (snow_depth-layer_bot_height)/(layer_top_height-layer_bot_height))))

tlai = (vai_bot-vai_top) * tree_lai / tree_vai
tsai = (vai_bot-vai_top) * tree_sai / tree_vai

if(present(tlai_layer)) tlai_layer = tlai
if(present(tsai_layer)) tsai_layer = tsai

elai_layer = fraction_exposed * tlai
esai_layer = fraction_exposed * tsai

! Update the vai at the bottom to be removed/decreased if there is no exposure
! set the vai top and bottoms to the snow layer if below
!vai_top = min(vai_top,fraction_exposed*tree_vai)

vai_bot = vai_top + fraction_exposed*(vai_bot-vai_top)

else

if(present(tlai_layer)) tlai_layer = 0._r8
if(present(tsai_layer)) tsai_layer = 0._r8
elai_layer = 0._r8
esai_layer = 0._r8
vai_bot = 0._r8
vai_top = 0._r8

end if if_any_vai


return
end subroutine VegAreaLayer

! =========================================================================

subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx)
Expand Down
2 changes: 2 additions & 0 deletions biogeochem/FatesCohortMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,6 +200,8 @@ module FatesCohortMod

!---------------------------------------------------------------------------

integer :: twostr_col ! The column index in the two-stream solution that this cohort is part of

! RESPIRATION COMPONENTS
real(r8) :: rdark ! dark respiration [kgC/indiv/s]
real(r8) :: resp_g_tstep ! growth respiration [kgC/indiv/timestep]
Expand Down
66 changes: 35 additions & 31 deletions biogeochem/FatesPatchMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,16 @@ module FatesPatchMod
use FatesLitterMod, only : litter_type
use PRTGenericMod, only : num_elements
use PRTGenericMod, only : element_list
use EDParamsMod, only : maxSWb, nlevleaf, nclmax, maxpft
use EDParamsMod, only : nlevleaf, nclmax, maxpft
use FatesConstantsMod, only : n_dbh_bins, n_dist_types
use FatesConstantsMod, only : n_rad_stream_types
use FatesConstantsMod, only : t_water_freeze_k_1atm
use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm
use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par
use FatesRunningMeanMod, only : ema_sdlng2sap_par, ema_sdlng_mdd

use TwoStreamMLPEMod, only : twostream_type
use FatesRadiationMemMod,only : num_swb
use FatesRadiationMemMod,only : num_rad_stream_types
use FatesInterfaceTypesMod,only : hlm_hio_ignore_val
use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=)
use shr_log_mod, only : errMsg => shr_log_errMsg

Expand Down Expand Up @@ -111,40 +113,36 @@ module FatesPatchMod
real(r8) :: c_stomata ! mean stomatal conductance of all leaves in the patch [umol/m2/s]
real(r8) :: c_lblayer ! mean boundary layer conductance of all leaves in the patch [umol/m2/s]

!TODO - can we delete these?
real(r8) :: layer_height_profile(nclmax,maxpft,nlevleaf)
real(r8) :: psn_z(nclmax,maxpft,nlevleaf)
real(r8) :: nrmlzd_parprof_pft_dir_z(n_rad_stream_types,nclmax,maxpft,nlevleaf)
real(r8) :: nrmlzd_parprof_pft_dif_z(n_rad_stream_types,nclmax,maxpft,nlevleaf)
real(r8) :: nrmlzd_parprof_dir_z(n_rad_stream_types,nclmax,nlevleaf)
real(r8) :: nrmlzd_parprof_dif_z(n_rad_stream_types,nclmax,nlevleaf)
real(r8) :: nrmlzd_parprof_pft_dir_z(num_rad_stream_types,nclmax,maxpft,nlevleaf)
real(r8) :: nrmlzd_parprof_pft_dif_z(num_rad_stream_types,nclmax,maxpft,nlevleaf)

!---------------------------------------------------------------------------

! RADIATION
real(r8) :: radiation_error ! radiation error [W/m2]
real(r8) :: rad_error(num_swb) ! radiation consv error by band [W/m2]
real(r8) :: fcansno ! fraction of canopy covered in snow [0-1]
logical :: solar_zenith_flag ! integer flag specifying daylight (based on zenith angle)
real(r8) :: solar_zenith_angle ! solar zenith angle [radians]
real(r8) :: gnd_alb_dif(maxSWb) ! ground albedo for diffuse rad, both bands [0-1]
real(r8) :: gnd_alb_dir(maxSWb) ! ground albedo for direct rad, both bands [0-1]

real(r8) :: gnd_alb_dif(num_swb) ! ground albedo for diffuse rad, both bands [0-1]
real(r8) :: gnd_alb_dir(num_swb) ! ground albedo for direct rad, both bands [0-1]


! organized by canopy layer, pft, and leaf layer
real(r8) :: fabd_sun_z(nclmax,maxpft,nlevleaf) ! sun fraction of direct light absorbed [0-1]
real(r8) :: fabd_sha_z(nclmax,maxpft,nlevleaf) ! shade fraction of direct light absorbed [0-1]
real(r8) :: fabi_sun_z(nclmax,maxpft,nlevleaf) ! sun fraction of indirect light absorbed [0-1]
real(r8) :: fabi_sha_z(nclmax,maxpft,nlevleaf) ! shade fraction of indirect light absorbed [0-1]
real(r8) :: ed_parsun_z(nclmax,maxpft,nlevleaf) ! PAR absorbed in the sun [W/m2]
real(r8) :: ed_parsha_z(nclmax,maxpft,nlevleaf) ! PAR absorbed in the shade [W/m2]
real(r8) :: f_sun(nclmax,maxpft,nlevleaf) ! fraction of leaves in the sun [0-1]
real(r8) :: ed_laisun_z(nclmax,maxpft,nlevleaf)
real(r8) :: ed_laisha_z(nclmax,maxpft,nlevleaf)
real(r8) :: f_sun(nclmax,maxpft,nlevleaf) ! fraction of leaves in the sun [0-1]


! radiation profiles for comparison against observations
real(r8) :: parprof_pft_dir_z(nclmax,maxpft,nlevleaf) ! direct-beam PAR profile through canopy, by canopy, PFT, leaf level [W/m2]
real(r8) :: parprof_pft_dif_z(nclmax,maxpft,nlevleaf) ! diffuse PAR profile through canopy, by canopy, PFT, leaf level [W/m2]
real(r8) :: parprof_dir_z(nclmax,nlevleaf) ! direct-beam PAR profile through canopy, by canopy, leaf level [W/m2]
real(r8) :: parprof_dif_z(nclmax,nlevleaf) ! diffuse PAR profile through canopy, by canopy, leaf level [W/m2]

real(r8), allocatable :: tr_soil_dir(:) ! fraction of incoming direct radiation transmitted to the soil as direct, by numSWB [0-1]
real(r8), allocatable :: tr_soil_dif(:) ! fraction of incoming diffuse radiation that is transmitted to the soil as diffuse [0-1]
Expand All @@ -155,6 +153,10 @@ module FatesPatchMod
real(r8), allocatable :: sabs_dir(:) ! fraction of incoming direct radiation that is absorbed by the canopy
real(r8), allocatable :: sabs_dif(:) ! fraction of incoming diffuse radiation that is absorbed by the canopy

! Twostream data structures
type(twostream_type) :: twostr ! This holds all two-stream data and procedures


!---------------------------------------------------------------------------

! ROOTS
Expand Down Expand Up @@ -318,16 +320,13 @@ subroutine NanValues(this)
this%ncan(:,:) = fates_unset_int
this%c_stomata = nan
this%c_lblayer = nan
this%layer_height_profile(:,:,:) = nan

this%psn_z(:,:,:) = nan
this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = nan
this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = nan
this%nrmlzd_parprof_dir_z(:,:,:) = nan
this%nrmlzd_parprof_dir_z(:,:,:) = nan

! RADIATION
this%radiation_error = nan
this%rad_error(:) = nan
this%fcansno = nan
this%solar_zenith_flag = .false.
this%solar_zenith_angle = nan
Expand All @@ -336,16 +335,14 @@ subroutine NanValues(this)
this%fabd_sun_z(:,:,:) = nan
this%fabd_sha_z(:,:,:) = nan
this%fabi_sun_z(:,:,:) = nan
this%fabi_sha_z(:,:,:) = nan
this%fabi_sha_z(:,:,:) = nan
this%ed_laisun_z(:,:,:) = nan
this%ed_laisha_z(:,:,:) = nan
this%ed_parsun_z(:,:,:) = nan
this%ed_parsha_z(:,:,:) = nan
this%f_sun(:,:,:) = nan
this%parprof_pft_dir_z(:,:,:) = nan
this%parprof_pft_dif_z(:,:,:) = nan
this%parprof_dir_z(:,:) = nan
this%parprof_dif_z(:,:) = nan
this%tr_soil_dir(:) = nan
this%tr_soil_dif(:) = nan
this%tr_soil_dir_dif(:) = nan
Expand Down Expand Up @@ -418,19 +415,17 @@ subroutine ZeroValues(this)
this%psn_z(:,:,:) = 0.0_r8
this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0.0_r8
this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0.0_r8
this%nrmlzd_parprof_dir_z(:,:,:) = 0.0_r8
this%nrmlzd_parprof_dif_z(:,:,:) = 0.0_r8

! RADIATION
this%radiation_error = 0.0_r8
this%rad_error(:) = 0.0_r8
this%fabd_sun_z(:,:,:) = 0.0_r8
this%fabd_sha_z(:,:,:) = 0.0_r8
this%fabi_sun_z(:,:,:) = 0.0_r8
this%fabi_sha_z(:,:,:) = 0.0_r8
this%ed_parsun_z(:,:,:) = 0.0_r8
this%ed_parsha_z(:,:,:) = 0.0_r8
this%ed_laisun_z(:,:,:) = 0.0_r8
this%ed_laisha_z(:,:,:) = 0.0_r8
this%ed_parsha_z(:,:,:) = 0.0_r8
this%ed_laisun_z(:,:,:) = 0._r8
this%ed_laisha_z(:,:,:) = 0._r8
this%f_sun = 0.0_r8
this%tr_soil_dir_dif(:) = 0.0_r8
this%fab(:) = 0.0_r8
Expand Down Expand Up @@ -591,7 +586,11 @@ subroutine Create(this, age, area, label, nocomp_pft, num_swb, num_pft, &

! initialize litter
call this%InitLitter(num_pft, num_levsoil)


this%twostr%scelg => null() ! The radiation module will check if this
! is associated, since it is not, it will then
! initialize and allocate

! assign known patch attributes
this%age = age
this%age_class = 1
Expand Down Expand Up @@ -646,7 +645,12 @@ subroutine FreeMemory(this, regeneration_model, numpft)
endif
ccohort => ncohort
end do


! Deallocate Radiation scattering elements
if(associated(this%twostr%scelg)) then
call this%twostr%DeallocTwoStream()
end if

! deallocate all litter objects
do el=1,num_elements
call this%litter(el)%DeallocateLitt()
Expand Down
Loading

0 comments on commit 7d518ac

Please sign in to comment.