From 1296c8618882f5fa2409aaccf87d58462ada0c07 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 30 Jan 2026 15:17:53 -0500 Subject: [PATCH] applying direct to parameter datastructure usage in allometry code --- biogeochem/FatesAllometryMod.F90 | 275 +++++++++++++----------------- biogeophys/LeafBiophysicsMod.F90 | 3 +- main/FatesInterfaceMod.F90 | 11 +- main/FatesParameterDerivedMod.F90 | 32 ++-- main/FatesParametersInterface.F90 | 195 +++++++++++++++++++++ main/JSONParameterUtilsMod.F90 | 43 +++-- radiation/TwoStreamMLPEMod.F90 | 1 - 7 files changed, 367 insertions(+), 193 deletions(-) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 94b58a881a..9a7feb018e 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -25,7 +25,6 @@ ! bbgw: biomass below ground woody tissues [kgC] ! similar to AGBW, this essentially encompasses ! all non-fineroot belowground tissues -! bfrmax: biomass in fine roots when "on allometry" [kgC] ! bsap: biomass in sapwood (above and below) [kgC] ! bdead: biomass (above and below) in the structural pool [kgC] ! @@ -83,7 +82,8 @@ module FatesAllometryMod ! If this is a unit-test, these globals will be provided by a wrapper - use PRTParametersMod, only : prt_params + ! This contains the pstruct data structure, and anything pid_ + use FatesParametersInterface use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : i4 => fates_int use FatesConstantsMod, only : g_per_kg @@ -199,7 +199,6 @@ subroutine CheckIntegratedAllometries(dbh,ipft,crowndamage, & logical,intent(out) :: l_pass ! Error flag (pass=true,no-pass=false) - real(r8) :: height ! diagnosed height [m] real(r8) :: bl_diag ! diagnosed leaf biomass [kgC] real(r8) :: bfr_diag ! diagnosed fine-root biomass [kgC] real(r8) :: bsap_diag ! diagnosed sapwood biomass [kgC] @@ -209,8 +208,6 @@ subroutine CheckIntegratedAllometries(dbh,ipft,crowndamage, & real(r8) :: bagw_diag ! diagnosed agbw [kgC] real(r8) :: bbgw_diag ! diagnosed below ground wood [kgC] - - l_pass = .true. ! Default assumption is that step passed if (grow_leaf) then @@ -303,10 +300,10 @@ subroutine h2d_allom(h,ipft,d,dddh) real(r8),intent(out) :: d ! plant diameter [cm] real(r8),intent(out),optional :: dddh ! change in diameter per height [cm/m] - associate( p1 => prt_params%allom_d2h1(ipft), & - p2 => prt_params%allom_d2h2(ipft), & - p3 => prt_params%allom_d2h3(ipft), & - allom_hmode => prt_params%allom_hmode(ipft)) + associate(p1 => pstruct%parameters(pid_allom_d2h1)%r_data_1d(ipft), & + p2 => pstruct%parameters(pid_allom_d2h2)%r_data_1d(ipft), & + p3 => pstruct%parameters(pid_allom_d2h3)%r_data_1d(ipft), & + allom_hmode => pstruct%parameters(pid_allom_hmode)%i_data_1d(ipft)) select case(allom_hmode) case (1) ! O'Brien et al 1995, BCI @@ -340,12 +337,12 @@ subroutine h_allom(d,ipft,h,dhdd) integer(i4),intent(in) :: ipft ! PFT index real(r8),intent(out) :: h ! plant height [m] real(r8),intent(out),optional :: dhdd ! change in height per diameter [m/cm] - - associate( dbh_maxh => prt_params%allom_dbh_maxheight(ipft), & - p1 => prt_params%allom_d2h1(ipft), & - p2 => prt_params%allom_d2h2(ipft), & - p3 => prt_params%allom_d2h3(ipft), & - allom_hmode => prt_params%allom_hmode(ipft)) + + associate( dbh_maxh => pstruct%parameters(pid_allom_dbh_maxh)%r_data_1d(ipft), & + p1 => pstruct%parameters(pid_allom_d2h1)%r_data_1d(ipft), & + p2 => pstruct%parameters(pid_allom_d2h2)%r_data_1d(ipft), & + p3 => pstruct%parameters(pid_allom_d2h3)%r_data_1d(ipft), & + allom_hmode => pstruct%parameters(pid_allom_hmode)%i_data_1d(ipft)) select case(allom_hmode) case (1) ! "obrien" @@ -388,22 +385,24 @@ subroutine bagw_allom(d,ipft,crowndamage, elongf_stem, bagw,dbagwdd) real(r8) :: dhdd ! change in height wrt d real(r8) :: crown_reduction ! crown reduction from damage real(r8) :: branch_frac ! fraction of aboveground woody biomass in branches - - associate( p1 => prt_params%allom_agb1(ipft), & - p2 => prt_params%allom_agb2(ipft), & - p3 => prt_params%allom_agb3(ipft), & - p4 => prt_params%allom_agb4(ipft), & - wood_density => prt_params%wood_density(ipft), & - c2b => prt_params%c2b(ipft), & - agb_frac => prt_params%allom_agb_frac(ipft), & - allom_amode => prt_params%allom_amode(ipft)) + + !pstruct%parameters(pid_allom_hmode)%i_data_1d + + associate( p1 => pstruct%parameters(pid_allom_agb1)%r_data_1d(ipft), & + p2 => pstruct%parameters(pid_allom_agb2)%r_data_1d(ipft), & + p3 => pstruct%parameters(pid_allom_agb3)%r_data_1d(ipft), & + p4 => pstruct%parameters(pid_allom_agb4)%r_data_1d(ipft), & + wood_density => pstruct%parameters(pid_wood_density)%r_data_1d(ipft), & + c2b => pstruct%parameters(pid_c2b)%r_data_1d(ipft), & + agb_frac => pstruct%parameters(pid_allom_agbfrac)%r_data_1d(ipft), & + allom_amode => pstruct%parameters(pid_allom_amode)%i_data_1d(ipft)) branch_frac = param_derived%branch_frac(ipft) select case(allom_amode) case (1) !"salda") call h_allom(d,ipft,h,dhdd) - call dh2bagw_salda(d,h,dhdd,p1,p2,p3,p4,wood_density,c2b,agb_frac,bagw,dbagwdd) + call dh2bagw_salda(d,h,dhdd,p1,p2,p3,p4,wood_density,agb_frac,bagw,dbagwdd) case (2) !"2par_pwr") ! Switch for woodland dbh->drc call d2bagw_2pwr(d,p1,p2,c2b,bagw,dbagwdd) @@ -456,18 +455,18 @@ subroutine blmax_allom(d,ipft,blmax,dblmaxdd) real(r8) :: h ! height real(r8) :: dhdd ! change in height wrt d - associate( dbh_maxh => prt_params%allom_dbh_maxheight(ipft), & - rho => prt_params%wood_density(ipft), & - slatop => prt_params%slatop(ipft), & - c2b => prt_params%c2b(ipft), & - allom_lmode => prt_params%allom_lmode(ipft), & - p1 => prt_params%allom_d2bl1(ipft), & - p2 => prt_params%allom_d2bl2(ipft), & - p3 => prt_params%allom_d2bl3(ipft)) - + associate(dbh_maxh => pstruct%parameters(pid_allom_dbh_maxh)%r_data_1d(ipft), & + rho => pstruct%parameters(pid_wood_density)%r_data_1d(ipft), & + slatop => pstruct%parameters(pid_leaf_slatop)%r_data_1d(ipft), & + c2b => pstruct%parameters(pid_c2b)%r_data_1d(ipft), & + allom_lmode => pstruct%parameters(pid_allom_lmode)%i_data_1d(ipft), & + p1 => pstruct%parameters(pid_allom_d2l1)%r_data_1d(ipft), & + p2 => pstruct%parameters(pid_allom_d2l2)%r_data_1d(ipft), & + p3 => pstruct%parameters(pid_allom_d2l3)%r_data_1d(ipft)) + select case(allom_lmode) case(1) !"salda") - call d2blmax_salda(d,p1,p2,p3,rho,dbh_maxh,c2b,blmax,dblmaxdd) + call d2blmax_salda(d,p1,p2,p3,rho,dbh_maxh,blmax,dblmaxdd) case(2) !"2par_pwr") call d2blmax_2pwr(d,p1,p2,c2b,blmax,dblmaxdd) case(3) ! dh2blmax_2pwr @@ -511,12 +510,12 @@ subroutine carea_allom(dbh,nplant,site_spread,ipft,crowndamage,c_area,inverse) ! crown area at height, we need to make ! special considerations - associate( dbh_maxh => prt_params%allom_dbh_maxheight(ipft), & - allom_lmode => prt_params%allom_lmode(ipft), & - d2bl_p2 => prt_params%allom_d2bl2(ipft), & - d2bl_ediff => prt_params%allom_blca_expnt_diff(ipft), & - d2ca_min => prt_params%allom_d2ca_coefficient_min(ipft), & - d2ca_max => prt_params%allom_d2ca_coefficient_max(ipft)) + associate( dbh_maxh => pstruct%parameters(pid_allom_dbh_maxh)%r_data_1d(ipft), & + allom_lmode => pstruct%parameters(pid_allom_lmode)%i_data_1d(ipft), & + d2bl_p2 => pstruct%parameters(pid_allom_d2l2)%r_data_1d(ipft), & + d2bl_ediff => pstruct%parameters(pid_allom_blca_ediff)%r_data_1d(ipft), & + d2ca_min => pstruct%parameters(pid_allom_d2ca_min)%r_data_1d(ipft), & + d2ca_max => pstruct%parameters(pid_allom_d2ca_max)%r_data_1d(ipft)) if( .not. present(inverse) ) then do_inverse = .false. @@ -686,8 +685,6 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25 real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground. real(r8) :: slat ! the sla of the top leaf layer. m2/kgC real(r8) :: canopy_lai_above ! total LAI of canopy layer overlying this tree - real(r8) :: vai_per_lai ! ratio of vegetation area index (ie. sai+lai) - ! to lai for individual tree real(r8) :: kn ! coefficient for exponential decay of 1/sla and ! vcmax with canopy depth real(r8) :: sla_max ! Observational constraint on how large sla @@ -704,7 +701,7 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25 call endrun(msg=errMsg(sourcefile, __LINE__)) endif - slat = g_per_kg * prt_params%slatop(pft) ! m2/g to m2/kg + slat = g_per_kg * pstruct%parameters(pid_leaf_slatop)%r_data_1d(pft) ! m2/g to m2/kg leafc_per_unitarea = leaf_c/(c_area/nplant) !KgC/m2 if(leafc_per_unitarea > 0.0_r8)then @@ -718,12 +715,12 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25 ! Coefficient for exponential decay of 1/sla with canopy depth: kn = DecayCoeffVcmax(vcmax25top, & - prt_params%leafn_vert_scaler_coeff1(pft), & - prt_params%leafn_vert_scaler_coeff2(pft)) + pstruct%parameters(pid_leafn_vert_scaler1)%r_data_1d(pft), & + pstruct%parameters(pid_leafn_vert_scaler2)%r_data_1d(pft)) ! take PFT-level maximum SLA value, even if under a thick canopy (which has units of m2/gC), ! and put into units of m2/kgC - sla_max = g_per_kg*prt_params%slamax(pft) + sla_max = g_per_kg*pstruct%parameters(pid_leaf_slamax)%r_data_1d(pft) ! Leafc_per_unitarea at which sla_max is reached due to exponential sla profile in canopy: leafc_slamax = (slat - sla_max * exp(-1.0_r8 * kn * canopy_lai_above)) / & (-1.0_r8 * kn * slat * sla_max) @@ -798,7 +795,7 @@ end function tree_lai ! ============================================================================ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, elongf_stem, c_area, nplant, & - cl, canopy_lai, treelai, vcmax25top, call_id ) + cl, canopy_lai, treelai, vcmax25top) ! ============================================================================ ! SAI of individual trees is a function of the LAI of individual trees @@ -816,9 +813,6 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, elongf_stem, c_ar ! each canopy layer real(r8), intent(in) :: treelai ! tree LAI for checking purposes only real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown - integer,intent(in) :: call_id ! flag specifying where this is called - ! from - real(r8) :: h real(r8) :: target_lai real(r8) :: target_bleaf @@ -829,7 +823,7 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, elongf_stem, c_ar target_lai = tree_lai(target_bleaf, pft, c_area, nplant, cl,& canopy_lai, vcmax25top) - tree_sai = elongf_stem * prt_params%allom_sai_scaler(pft) * target_lai + tree_sai = elongf_stem * pstruct%parameters(pid_allom_sai_scaler)%r_data_1d(pft) * target_lai return end function tree_sai @@ -837,7 +831,7 @@ end function tree_sai ! ============================================================================ subroutine tree_lai_sai(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top, & - dbh, crowndamage, canopy_trim, elongf_stem, call_id, & + dbh, crowndamage, canopy_trim, elongf_stem, & treelai, treesai) ! This is the public wrapper for calling plant lai and sai @@ -857,9 +851,6 @@ subroutine tree_lai_sai(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top, integer, intent(in) :: crowndamage real(r8), intent(in) :: canopy_trim ! trimming function (0-1) real(r8), intent(in) :: elongf_stem ! Elongation factor for stems. - integer,intent(in) :: call_id ! flag specifying where this is called from - ! from - real(r8), intent(out) :: treelai ! plant LAI [m2 leaf area/m2 crown area] real(r8), intent(out) :: treesai ! plant SAI [m2 stem area/m2 crown area] @@ -871,13 +862,13 @@ subroutine tree_lai_sai(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top, treelai = tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top) treesai = tree_sai( pft, dbh, crowndamage, canopy_trim, elongf_stem, c_area, nplant, & - cl, canopy_lai, treelai, vcmax25top, call_id ) + cl, canopy_lai, treelai, vcmax25top) ! Don't allow lai+sai to exceed the vertical discretization bounds if( do_vai_capping ) then if( (treelai + treesai) > (sum(dinc_vai)) )then - treelai = sum(dinc_vai) * (1._r8 - prt_params%allom_sai_scaler(pft)) - nearzero - treesai = sum(dinc_vai) * prt_params%allom_sai_scaler(pft) - nearzero + treelai = sum(dinc_vai) * (1._r8 - pstruct%parameters(pid_allom_sai_scaler)%r_data_1d(pft)) - nearzero + treesai = sum(dinc_vai) * pstruct%parameters(pid_allom_sai_scaler)%r_data_1d(pft) - nearzero end if end if @@ -905,18 +896,13 @@ real(r8) function leafc_from_treelai( treelai, treesai, pft, c_area, nplant, cl, ! top, ref 25C ! !LOCAL VARIABLES: - real(r8) :: leaf_c ! plant leaf carbon [kg] real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground. real(r8) :: slat ! the sla of the top leaf layer. m2/kgC - real(r8) :: vai_per_lai ! ratio of vegetation area index (ie. sai+lai) - ! to lai for individual tree real(r8) :: kn ! coefficient for exponential decay of 1/sla and ! vcmax with canopy depth real(r8) :: sla_max ! Observational constraint on how large sla ! (m2/gC) can become real(r8) :: leafc_slamax ! Leafc_per_unitarea at which sla_max is reached - real(r8) :: clim ! Upper limit for leafc_per_unitarea in exponential - ! tree_lai function real(r8) :: tree_lai_at_slamax ! lai at which we reach the maximum sla value. real(r8) :: leafc_linear_phase ! amount of leaf carbon needed to get to the target treelai ! when the slamax value has been reached (i.e. deep layers with unchanging sla) @@ -943,13 +929,13 @@ real(r8) function leafc_from_treelai( treelai, treesai, pft, c_area, nplant, cl, end if ! convert PFT-level canopy top and maximum SLA values and convert from m2/gC to m2/kgC - slat = g_per_kg * prt_params%slatop(pft) - sla_max = g_per_kg * prt_params%slamax(pft) + slat = g_per_kg * pstruct%parameters(pid_leaf_slatop)%r_data_1d(pft) + sla_max = g_per_kg * pstruct%parameters(pid_leaf_slamax)%r_data_1d(pft) ! Coefficient for exponential decay of 1/sla with canopy depth: kn = DecayCoeffVcmax(vcmax25top, & - prt_params%leafn_vert_scaler_coeff1(pft), & - prt_params%leafn_vert_scaler_coeff2(pft)) + pstruct%parameters(pid_leafn_vert_scaler1)%r_data_1d(pft), & + pstruct%parameters(pid_leafn_vert_scaler2)%r_data_1d(pft)) if(treelai > 0.0_r8)then ! Leafc_per_unitarea at which sla_max is reached due to exponential sla profile in canopy: @@ -1023,11 +1009,11 @@ subroutine bsap_allom(d,ipft,crowndamage,canopy_trim,elongf_stem, sapw_area,bsap ! X% of total woody/fibrous (ie non leaf/fineroot) tissues real(r8),parameter :: max_frac = 0.95_r8 - agb_frac = prt_params%allom_agb_frac(ipft) + agb_frac = pstruct%parameters(pid_allom_agbfrac)%r_data_1d(ipft) branch_frac = param_derived%branch_frac(ipft) - select case(prt_params%allom_smode(ipft)) + select case(pstruct%parameters(pid_allom_smode)%i_data_1d(ipft)) ! --------------------------------------------------------------------- ! Currently only one sapwood allometry model. the slope ! of the la:sa to diameter line is zero. @@ -1039,7 +1025,7 @@ subroutine bsap_allom(d,ipft,crowndamage,canopy_trim,elongf_stem, sapw_area,bsap ! (but could be modulated by stem phenology). call h_allom(d,ipft,h,dhdd) call bleaf(d,ipft,1,canopy_trim,1.0_r8,bl,dbldd) - call bsap_ltarg_slatop(d,h,dhdd,bl,dbldd,ipft,sapw_area,bsap,dbsapdd) + call bsap_ltarg_slatop(h,dhdd,bl,dbldd,ipft,sapw_area,bsap,dbsapdd) ! if trees are damaged reduce bsap by percent crown loss * ! fraction of biomass that would be in branches (pft specific) @@ -1098,7 +1084,7 @@ subroutine bsap_allom(d,ipft,crowndamage,canopy_trim,elongf_stem, sapw_area,bsap case DEFAULT write(fates_log(),*) 'An undefined sapwood allometry was specified: ', & - prt_params%allom_smode(ipft) + pstruct%parameters(pid_allom_smode)%i_data_1d(ipft) write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end select @@ -1122,17 +1108,17 @@ subroutine bbgw_allom(d,ipft,elongf_stem,bbgw,dbbgwdd) real(r8) :: bagw ! above ground biomass [kgC] real(r8) :: dbagwdd ! change in agb per diameter [kgC/cm] - select case(prt_params%allom_cmode(ipft)) + select case(pstruct%parameters(pid_allom_cmode)%i_data_1d(ipft)) case(1) !"constant") ! bbgw not affected by damage so use target allometry no damage. But note that bbgw ! is affected by stem phenology (typically applied only to grasses). We do not need ! to account for stem phenology in bbgw_const because bbgw will be proportional to ! bagw, and bagw is downscaled due to stem phenology. call bagw_allom(d,ipft,1, elongf_stem, bagw,dbagwdd) - call bbgw_const(d,bagw,dbagwdd,ipft,bbgw,dbbgwdd) + call bbgw_const(bagw,dbagwdd,ipft,bbgw,dbbgwdd) case DEFAULT write(fates_log(),*) 'An undefined coarse root allometry was specified: ', & - prt_params%allom_cmode(ipft) + pstruct%parameters(pid_allom_cmode)%i_data_1d(ipft) write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end select @@ -1163,11 +1149,8 @@ subroutine bfineroot(d,ipft,canopy_trim,l2fr,elongf_fnrt,bfr,dbfrdd) real(r8) :: blmax ! maximum leaf biomss per allometry real(r8) :: dblmaxdd - real(r8) :: bfrmax - real(r8) :: dbfrmaxdd - real(r8) :: slascaler - - select case(prt_params%allom_fmode(ipft)) + + select case( pstruct%parameters(pid_allom_fmode)%i_data_1d(ipft) ) case(1) ! "constant proportionality with TRIMMED target bleaf" call blmax_allom(d,ipft,blmax,dblmaxdd) @@ -1189,7 +1172,7 @@ subroutine bfineroot(d,ipft,canopy_trim,l2fr,elongf_fnrt,bfr,dbfrdd) case DEFAULT write(fates_log(),*) 'An undefined fine root allometry was specified: ', & - prt_params%allom_fmode(ipft) + pstruct%parameters(pid_allom_fmode)%i_data_1d(ipft) write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end select @@ -1201,11 +1184,9 @@ subroutine bfineroot(d,ipft,canopy_trim,l2fr,elongf_fnrt,bfr,dbfrdd) dbfrdd = elongf_fnrt * dbfrdd end if - return end subroutine bfineroot - ! ============================================================================ ! Storage biomass interface ! ============================================================================ @@ -1225,19 +1206,19 @@ subroutine bstore_allom(d,ipft,crowndamage, canopy_trim,bstore,dbstoredd) real(r8) :: dblmaxdd ! Allometric target change in leaf biomass per cm (UNTRIMMED) - associate( allom_stmode => prt_params%allom_stmode(ipft), & - cushion => prt_params%cushion(ipft) ) + associate( allom_stmode => pstruct%parameters(pid_allom_stmode)%i_data_1d(ipft), & + cushion => pstruct%parameters(pid_storage_cushion)%r_data_1d(ipft) ) select case(allom_stmode) case(1) ! Storage is constant proportionality of trimmed maximum leaf ! biomass (ie cushion * bleaf), and thus leaf phenology is ignored. call bleaf(d,ipft, crowndamage, canopy_trim, 1.0_r8, bl, dbldd) - call bstore_blcushion(d,bl,dbldd,cushion,ipft,bstore,dbstoredd) + call bstore_blcushion(bl,dbldd,cushion,bstore,dbstoredd) case(2) ! Storage is constant proportionality of untrimmed maximum leaf ! biomass (ie cushion * bleaf_max) call blmax_allom(d,ipft,blmax,dblmaxdd) - call bstore_blcushion(d,blmax,dblmaxdd,cushion,ipft,bstore,dbstoredd) + call bstore_blcushion(blmax,dblmaxdd,cushion,bstore,dbstoredd) case DEFAULT write(fates_log(),*) 'An undefined fine storage allometry was specified: ', & @@ -1255,7 +1236,7 @@ end subroutine bstore_allom ! ============================================================================ ! Dead biomass interface ! ============================================================================ - + subroutine bdead_allom(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeaddd) real(r8),intent(in) :: bagw ! biomass above ground wood (agb) [kgC] @@ -1277,9 +1258,9 @@ subroutine bdead_allom(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeadd ! bbgw. Therefore, it is not removed from AGB and BBGW in the calculation of dead mass. - associate( agb_fraction => prt_params%allom_agb_frac(ipft)) + associate(agb_fraction => pstruct%parameters(pid_allom_agbfrac)%r_data_1d(ipft)) - select case(prt_params%allom_amode(ipft)) + select case(pstruct%parameters(pid_allom_amode)%i_data_1d(ipft)) case(1) ! Saldariagga mass allometry originally calculated bdead directly. ! we assume proportionality between bdead and bagw @@ -1299,7 +1280,7 @@ subroutine bdead_allom(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeadd case DEFAULT write(fates_log(),*) 'An undefined AGB allometry was specified: ',& - prt_params%allom_amode(ipft) + pstruct%parameters(pid_allom_amode)%i_data_1d(ipft) write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -1313,16 +1294,15 @@ end subroutine bdead_allom ! Specific bbgw relationships ! ============================================================================ - subroutine bbgw_const(d,bagw,dbagwdd,ipft,bbgw,dbbgwdd) + subroutine bbgw_const(bagw,dbagwdd,ipft,bbgw,dbbgwdd) - real(r8),intent(in) :: d ! plant diameter [cm] real(r8),intent(in) :: bagw ! above ground biomass [kg] real(r8),intent(in) :: dbagwdd ! change in agb per diameter [kg/cm] integer(i4),intent(in) :: ipft ! PFT index real(r8),intent(out) :: bbgw ! coarse root biomass [kg] real(r8),intent(out),optional :: dbbgwdd ! change croot bio per diam [kg/cm] - associate( agb_fraction => prt_params%allom_agb_frac(ipft) ) + associate( agb_fraction => pstruct%parameters(pid_allom_agbfrac)%r_data_1d(ipft) ) bbgw = (1.0_r8/agb_fraction-1.0_r8)*bagw @@ -1340,7 +1320,7 @@ end subroutine bbgw_const ! Specific d2bsap relationships ! ============================================================================ - subroutine bsap_ltarg_slatop(d,h,dhdd,bleaf,dbleafdd,ipft, & + subroutine bsap_ltarg_slatop(h,dhdd,bleaf,dbleafdd,ipft, & sapw_area,bsap,dbsapdd) ! ------------------------------------------------------------------------- @@ -1358,7 +1338,6 @@ subroutine bsap_ltarg_slatop(d,h,dhdd,bleaf,dbleafdd,ipft, & ! ! ------------------------------------------------------------------------- - real(r8),intent(in) :: d ! plant diameter [cm] real(r8),intent(in) :: h ! plant height [m] real(r8),intent(in) :: dhdd ! change in height per diameter [m/cm] real(r8),intent(in) :: bleaf ! plant leaf target biomass [kgC] @@ -1377,12 +1356,12 @@ subroutine bsap_ltarg_slatop(d,h,dhdd,bleaf,dbleafdd,ipft, & real(r8) :: hbl2bsap ! sapwood biomass per lineal height - associate ( la_per_sa_int => prt_params%allom_la_per_sa_int(ipft), & - la_per_sa_slp => prt_params%allom_la_per_sa_slp(ipft), & - slatop => prt_params%slatop(ipft), & - wood_density => prt_params%wood_density(ipft), & - c2b => prt_params%c2b(ipft), & - agb_fraction => prt_params%allom_agb_frac(ipft) ) + associate ( la_per_sa_int => pstruct%parameters(pid_allom_lapersa_int)%r_data_1d(ipft), & + la_per_sa_slp => pstruct%parameters(pid_allom_lapersa_slp)%r_data_1d(ipft), & + slatop => pstruct%parameters(pid_leaf_slatop)%r_data_1d(ipft), & + wood_density => pstruct%parameters(pid_wood_density)%r_data_1d(ipft), & + c2b => pstruct%parameters(pid_c2b)%r_data_1d(ipft), & + agb_fraction => pstruct%parameters(pid_allom_agbfrac)%r_data_1d(ipft)) ! Calculate sapwood biomass per linear height and kgC of leaf [m-1] @@ -1490,17 +1469,15 @@ end subroutine SapwoodAreaGrass ! Specific storage relationships ! ============================================================================ - subroutine bstore_blcushion(d,bl,dbldd,cushion,ipft,bstore,dbstoredd) + subroutine bstore_blcushion(bl,dbldd,cushion,bstore,dbstoredd) ! This discracefully simple subroutine calculates allometric target ! storage biomass based on a constant-specified ratio (cushion) ! of storage to target allometricc leaf biomass - real(r8),intent(in) :: d ! plant diameter [cm] real(r8),intent(in) :: bl ! plant leaf biomass [kgC] real(r8),intent(in) :: dbldd ! change in blmax per diam [kgC/cm] real(r8),intent(in) :: cushion ! simple constant ration bstore/bleaf - integer(i4),intent(in) :: ipft ! PFT index real(r8),intent(out) :: bstore ! plant leaf biomass [kgC] real(r8),intent(out),optional :: dbstoredd ! change leaf bio per diameter [kgC/cm] @@ -1519,7 +1496,7 @@ end subroutine bstore_blcushion ! Specific d2blmax relationships ! ============================================================================ - subroutine d2blmax_salda(d,p1,p2,p3,rho,dbh_maxh,c2b,blmax,dblmaxdd) + subroutine d2blmax_salda(d,p1,p2,p3,rho,dbh_maxh,blmax,dblmaxdd) real(r8),intent(in) :: d ! plant diameter [cm] @@ -1528,7 +1505,6 @@ subroutine d2blmax_salda(d,p1,p2,p3,rho,dbh_maxh,c2b,blmax,dblmaxdd) real(r8),intent(in) :: p3 real(r8),intent(in) :: rho ! plant's wood specific gravity real(r8),intent(in) :: dbh_maxh ! dbh at which max height occurs - real(r8),intent(in) :: c2b ! c to biomass multiplier (~2.0) real(r8),intent(out) :: blmax ! plant leaf biomass [kg] real(r8),intent(out),optional :: dblmaxdd ! change leaf bio per diam [kgC/cm] @@ -2305,7 +2281,7 @@ end subroutine d2bagw_2pwr subroutine dh2bagw_salda(d,h,dhdd,p1,p2,p3,p4, & - wood_density,c2b,allom_agb_frac,bagw,dbagwdd) + wood_density,allom_agb_frac,bagw,dbagwdd) ! -------------------------------------------------------------------- ! Calculate stem biomass from height(m) dbh(cm) and wood density(g/cm3) @@ -2323,7 +2299,6 @@ subroutine dh2bagw_salda(d,h,dhdd,p1,p2,p3,p4, & real(r8),intent(in) :: p2 ! = 0.572_r8 real(r8),intent(in) :: p3 ! = 1.94_r8 real(r8),intent(in) :: p4 ! = 0.931_r8 - real(r8),intent(in) :: c2b ! carbon 2 biomass ratio real(r8),intent(in) :: wood_density real(r8),intent(in) :: allom_agb_frac real(r8),intent(out) :: bagw ! plant biomass [kgC/indiv] @@ -2376,7 +2351,7 @@ subroutine h2d_chave2014(h,p1,p2,p3,de,ddedh) real(r8),intent(out) :: de ! effective plant diameter [cm] real(r8),intent(out),optional :: ddedh ! effective change in d per height [cm/m] - real(r8) :: p1e, eroot, dbh1,dhpdd + real(r8) :: p1e, eroot, dhpdd p1e = p1 !-eclim (assumed that p1 already has eclim removed) eroot = (-p2 + sqrt(p2**2.0_r8 + 4.0_r8*log(h*exp(-p1e))*p3)) & @@ -2567,9 +2542,9 @@ subroutine CrownDepth(height,ipft,crown_depth) integer ,intent(in) :: ipft ! functional type index real(r8),intent(out) :: crown_depth ! The depth of the crown [m] - associate( p1 => prt_params%allom_h2cd1(ipft), & - p2 => prt_params%allom_h2cd2(ipft), & - allom_dmode => prt_params%allom_dmode(ipft)) + associate( p1 => pstruct%parameters(pid_allom_h2cd1)%r_data_1d(ipft), & + p2 => pstruct%parameters(pid_allom_h2cd2)%r_data_1d(ipft), & + allom_dmode => pstruct%parameters(pid_allom_dmode)%i_data_1d(ipft)) select case (allom_dmode) case (1) ! Default, linear relationship with height @@ -2777,17 +2752,8 @@ subroutine set_root_fraction(root_fraction, ft, zi, max_nlevroot) ! like permafrost. If so, compress profile over the maximum depth integer,optional, intent(in) :: max_nlevroot - - ! locals - real(r8) :: a_par ! local temporary for "a" parameter - real(r8) :: b_par ! "" "b" parameter - ! Parameters ! - ! TO-DO: NEXT TIME WE ROLL OUT A NEW PARAMETER INTERFACE, ADD - ! PROFILE SWAPPING FLAGS. OR IF THERE IS NO DEMAND< LEAVE AS IS. - ! - ! ! Two context exist 'hydraulic' and 'biomass'. This allows us to ! allow different profiles for how water is drawn from the soil ! and different profiles to define the biomass for litter flux. @@ -2801,7 +2767,6 @@ subroutine set_root_fraction(root_fraction, ft, zi, max_nlevroot) integer, parameter :: exponential_1p_profile_type = 2 integer, parameter :: exponential_2p_profile_type = 3 - integer :: root_profile_type integer :: corr_id(1) ! This is the bin with largest fraction ! add/subtract any corrections there integer :: nlevroot @@ -2830,14 +2795,17 @@ subroutine set_root_fraction(root_fraction, ft, zi, max_nlevroot) nlevroot = min(max_nlevroot,nlevroot) end if - select case(nint(prt_params%fnrt_prof_mode(ft))) + select case(pstruct%parameters(pid_allom_fnrt_prof_mode)%i_data_1d(ft)) case ( exponential_1p_profile_type ) - call exponential_1p_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), prt_params%fnrt_prof_a(ft)) + call exponential_1p_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), & + pstruct%parameters(pid_allom_fnrt_prof_a)%r_data_1d(ft)) case ( jackson_beta_profile_type ) - call jackson_beta_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), prt_params%fnrt_prof_a(ft)) + call jackson_beta_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), & + pstruct%parameters(pid_allom_fnrt_prof_b)%r_data_1d(ft)) case ( exponential_2p_profile_type ) call exponential_2p_root_profile(root_fraction(1:nlevroot), zi(0:nlevroot), & - prt_params%fnrt_prof_a(ft),prt_params%fnrt_prof_b(ft)) + pstruct%parameters(pid_allom_fnrt_prof_a)%r_data_1d(ft), & + pstruct%parameters(pid_allom_fnrt_prof_b)%r_data_1d(ft)) case default write(fates_log(),*) 'An undefined root profile type was specified' @@ -2924,7 +2892,6 @@ subroutine exponential_1p_root_profile(root_fraction, zi, a) ! LOCAL VARIABLES: integer :: lev ! soil depth layer index integer :: nlevsoil ! number of soil layers - real(r8) :: depth ! Depth to middle of layer [m] real(r8) :: sum_rootfr ! sum of rooting profile for normalization ! Typical default parameter is a = 3. @@ -3031,7 +2998,7 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d ! Do reduce "if" calls, we break this call into two parts - if ( prt_params%woody(ipft) == itrue ) then + if ( pstruct%parameters(pid_woody)%i_data_1d(ipft) == itrue ) then if(.not.present(bdead)) then write(fates_log(),*) 'woody plants must use structure for dbh reset' @@ -3121,9 +3088,9 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d call h_allom(d,ipft,h) if(counter>20)then write(fates_log(),*) 'dbh counter: ',counter,' is woody: ',& - (prt_params%woody(ipft) == itrue) + (pstruct%parameters(pid_woody)%i_data_1d(ipft) == itrue) - if(prt_params%woody(ipft)==itrue)then + if(pstruct%parameters(pid_woody)%i_data_1d(ipft)==itrue)then warn_msg = 'dbh counter: '//trim(I2S(counter))//' is woody' else warn_msg = 'dbh counter: '//trim(I2S(counter))//' is not woody' @@ -3170,8 +3137,6 @@ subroutine VegAreaLayer(tree_lai,tree_sai,tree_height,iv,nv,pft,snow_depth, & 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 @@ -3248,7 +3213,7 @@ end subroutine VegAreaLayer ! ========================================================================= - subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx) + !!subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx) ! ============================================================================ ! This subroutine performs a cubic spline interpolation between known @@ -3262,29 +3227,29 @@ subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx) ! Arguments - real(r8),intent(in) :: x1 ! Lower endpoint independent - real(r8),intent(in) :: x2 ! Upper endpoint independent - real(r8),intent(in) :: y1 ! Lower endpoint dependent - real(r8),intent(in) :: y2 ! Upper endpoint dependent - real(r8),intent(in) :: dydx1 ! Lower endpoint slope - real(r8),intent(in) :: dydx2 ! Upper endpoint slope - real(r8),intent(in) :: x ! Independent - real(r8),intent(out) :: y ! Dependent - real(r8),intent(out) :: dydx ! Slope + !! real(r8),intent(in) :: x1 ! Lower endpoint independent + !! real(r8),intent(in) :: x2 ! Upper endpoint independent + !! real(r8),intent(in) :: y1 ! Lower endpoint dependent + !! real(r8),intent(in) :: y2 ! Upper endpoint dependent + !! real(r8),intent(in) :: dydx1 ! Lower endpoint slope + !! real(r8),intent(in) :: dydx2 ! Upper endpoint slope + !! real(r8),intent(in) :: x ! Independent + !! real(r8),intent(out) :: y ! Dependent + !! real(r8),intent(out) :: dydx ! Slope ! Temps - real(r8) :: t - real(r8) :: a - real(r8) :: b - - t = (x-x1)/(x2-x1) - a = dydx1*(x2-x1) - (y2-y1) - b = -dydx2*(x2-x1) + (y2-y1) - - y = (1.0_r8-t)*y1 + t*y2 + t*(1.0_r8-t)*(a*(1.0_r8-t) + b*t) - dydx = (y2-y1)/(x2-x1) + (1.0_r8-2.0_r8*t)*(a*(1.0_r8-t)+b*t)/(x2-x1) + t*(1.0_r8-t)*(b-a)/(x2-x1) - return - end subroutine cspline + !! real(r8) :: t + !! real(r8) :: a + !! real(r8) :: b + + !! t = (x-x1)/(x2-x1) + !! a = dydx1*(x2-x1) - (y2-y1) + !! b = -dydx2*(x2-x1) + (y2-y1) + + !! y = (1.0_r8-t)*y1 + t*y2 + t*(1.0_r8-t)*(a*(1.0_r8-t) + b*t) + !! dydx = (y2-y1)/(x2-x1) + (1.0_r8-2.0_r8*t)*(a*(1.0_r8-t)+b*t)/(x2-x1) + t*(1.0_r8-t)*(b-a)/(x2-x1) + !! return + !!end subroutine cspline ! ============================================================================ diff --git a/biogeophys/LeafBiophysicsMod.F90 b/biogeophys/LeafBiophysicsMod.F90 index c238eebf16..fa0edd77cf 100644 --- a/biogeophys/LeafBiophysicsMod.F90 +++ b/biogeophys/LeafBiophysicsMod.F90 @@ -950,7 +950,6 @@ subroutine CiFunc(ci, & real(r8), intent(out) :: fval ! ci_input - ci_updated (Pa) !real(r8) :: veg_esat ! Saturation vapor pressure at leaf-surface [Pa] - real(r8) :: veg_qs ! DUMMY, specific humidity at leaf-surface [kg/kg] real(r8) :: a_gs ! The assimilation (a) for calculating conductance (gs) ! is either = to anet or agross real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) @@ -1478,7 +1477,7 @@ function LeafHumidityStomaResis(leaf_psi, k_lwp, veg_tempk, can_vpress, can_pres lwp_star = 1._r8 end if - ! call QSat(veg_tempk, can_press, qsat_alt, veg_esat) + call QSat(veg_tempk, can_press, qsat_alt, veg_esat) qsat_alt = qsat_alt * g_per_kg diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3ca8220ddb..9b6ca44b91 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -821,11 +821,15 @@ subroutine SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft,paramfile) end do write(fates_log(),*) '============ End FATES Parameter Info =========' end if + + ! This assigns named constants to the indices of each + ! parameter datastructure, allowing data to be accessed + ! directly from the structure + call SetParameterIndices() ! This call transfers parameters from the pstruct data-structure ! into the specific datastructures where parameters have there ! own primitive arrays - call FatesTransferParameters() fates_numpft = size(prt_params%wood_density,dim=1) @@ -2338,8 +2342,9 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) do while (associated(ccohort)) ! call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) if(.not.ccohort%isnew)then - ! [kgC/plant/yr] -> [gC/m2/yr] - site_npp = site_npp + ccohort%npp_acc_hold * ccohort%n*area_inv * g_per_kg + ! [kgC/plant/yr] -> [gC/m2/s] + site_npp = site_npp + ccohort%npp_acc_hold * ccohort%n*area_inv * & + g_per_kg * hlm_days_per_year / sec_per_day end if ccohort => ccohort%shorter end do diff --git a/main/FatesParameterDerivedMod.F90 b/main/FatesParameterDerivedMod.F90 index a044c33721..4e4204550f 100644 --- a/main/FatesParameterDerivedMod.F90 +++ b/main/FatesParameterDerivedMod.F90 @@ -15,8 +15,9 @@ module FatesParameterDerivedMod use FatesInterfaceTypesMod, only : nleafage use FatesInterfaceTypesMod, only : nlevdamage use FatesGlobals , only : fates_log - use EDParamsMod , only : ED_val_history_damage_bin_edges - + use FatesParametersInterface + + implicit none private @@ -81,21 +82,21 @@ end subroutine InitAllocateDamageTransitions ! ===================================================================================== - subroutine Init(this,numpft) + subroutine Init(this) - use EDPftvarcon, only: EDPftvarcon_inst - use SFParamsMod, only: SF_val_CWD_frac use FatesLitterMod, only : ncwd class(param_derived_type), intent(inout) :: this - integer, intent(in) :: numpft ! local variables integer :: ft ! pft index integer :: iage ! leaf age class index + integer :: numpft - associate( vcmax25top => EDPftvarcon_inst%vcmax25top ) - + associate( vcmax25top => pstruct%parameters(pid_vcmax25top)%r_data_2d ) + + numpft = size(vcmax25top,dim=2) + call this%InitAllocate(numpft) call this%InitDamageTransitions(numpft) @@ -114,14 +115,14 @@ subroutine Init(this,numpft) ! jmax25top(ft) = & ! (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrzc),11._r8),35._r8)) * vcmax25top(ft) - this%jmax25top(ft,iage) = 1.67_r8 * vcmax25top(ft,iage) - this%tpu25top(ft,iage) = 0.167_r8 * vcmax25top(ft,iage) - this%kp25top(ft,iage) = 20000._r8 * vcmax25top(ft,iage) + this%jmax25top(ft,iage) = 1.67_r8 * vcmax25top(iage,ft) + this%tpu25top(ft,iage) = 0.167_r8 * vcmax25top(iage,ft) + this%kp25top(ft,iage) = 20000._r8 * vcmax25top(iage,ft) end do ! Allocate fraction of aboveground woody biomass in branches - this%branch_frac(ft) = sum(SF_val_CWD_frac(1:3)) + this%branch_frac(ft) = sum(pstruct%parameters(pid_cwd_frac)%r_data_1d(1:3),dim=1) end do !ft @@ -133,9 +134,6 @@ end subroutine Init subroutine InitDamageTransitions(this, numpft) - use EDPftvarcon, only: EDPftvarcon_inst - - class(param_derived_type), intent(inout) :: this integer, intent(in) :: numpft @@ -155,7 +153,7 @@ subroutine InitDamageTransitions(this, numpft) ! class widths ! append 100 to ED_val_history_damage_bin_edges do j = 1,nlevdamage - damage_bin_edges_ex(j) = ED_val_history_damage_bin_edges(j) + damage_bin_edges_ex(j) = pstruct%parameters(pid_damage_bins)%r_data_1d(j) !ED_val_history_damage_bin_edges(j) end do damage_bin_edges_ex(j) = 100.0_r8 @@ -165,7 +163,7 @@ subroutine InitDamageTransitions(this, numpft) do ft = 1, numpft - damage_frac = EDPftvarcon_inst%damage_frac(ft) + damage_frac = pstruct%parameters(pid_damage_frac)%r_data_1d(ft) do i = 1, nlevdamage diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index ed3fa26591..8832e37c1b 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -8,6 +8,7 @@ module FatesParametersInterface ! only uses shr libraries use JSONParameterUtilsMod, only: params_type + use JSONParameterUtilsMod, only: param_type use FatesConstantsMod, only : r8 => fates_r8 implicit none @@ -15,13 +16,207 @@ module FatesParametersInterface type(params_type) :: pstruct + ! Parameter indexes + + integer, public :: pid_vcmax25top + integer, public :: pid_cwd_frac + integer, public :: pid_damage_frac + integer, public :: pid_damage_bins + integer, public :: pid_alloc_organ_id + integer, public :: pid_phen_leaf_habit + integer, public :: pid_phen_stem_drop_frac + integer, public :: pid_phen_fnrt_drop_frac + integer, public :: pid_phen_mindaysoff + integer, public :: pid_phen_drought_thresh + integer, public :: pid_phen_moist_thresh + integer, public :: pid_leaf_slamax + integer, public :: pid_leaf_slatop + integer, public :: pid_allom_sai_scaler + integer, public :: pid_allom_fnrt_prof_a + integer, public :: pid_allom_fnrt_prof_b + integer, public :: pid_allom_fnrt_prof_mode + integer, public :: pid_woody + integer, public :: pid_wood_density + integer, public :: pid_recruit_seed_dbh_thresh + integer, public :: pid_storage_cushion + integer, public :: pid_store_priority_frac + integer, public :: pid_turnover_senleaf_fdrought + integer, public :: pid_turnover_fnrt + integer, public :: pid_leafn_vert_scaler1 + integer, public :: pid_leafn_vert_scaler2 + integer, public :: pid_recruit_seed_alloc_mature + integer, public :: pid_recruit_seed_alloc + integer, public :: pid_trs_repro_alloc_a + integer, public :: pid_trs_repro_alloc_b + integer, public :: pid_c2b + integer, public :: pid_grperc + integer, public :: pid_allom_dbh_maxh + integer, public :: pid_allom_hmode + integer, public :: pid_allom_lmode + integer, public :: pid_allom_fmode + integer, public :: pid_allom_amode + integer, public :: pid_allom_stmode + integer, public :: pid_allom_cmode + integer, public :: pid_allom_smode + integer, public :: pid_allom_dmode + integer, public :: pid_allom_lapersa_int + integer, public :: pid_allom_lapersa_slp + integer, public :: pid_allom_l2fr + integer, public :: pid_cnp_kp + integer, public :: pid_cnp_ki + integer, public :: pid_cnp_kd + integer, public :: pid_cnp_store_ovrflw_frac + integer, public :: pid_cnp_nfix1 + integer, public :: pid_allom_d2h1 + integer, public :: pid_allom_d2h2 + integer, public :: pid_allom_d2h3 + integer, public :: pid_allom_agbfrac + integer, public :: pid_allom_d2l1 + integer, public :: pid_allom_d2l2 + integer, public :: pid_allom_d2l3 + integer, public :: pid_allom_blca_ediff + integer, public :: pid_allom_d2ca_max + integer, public :: pid_allom_d2ca_min + integer, public :: pid_allom_agb1 + integer, public :: pid_allom_agb2 + integer, public :: pid_allom_agb3 + integer, public :: pid_allom_agb4 + integer, public :: pid_allom_h2cd1 + integer, public :: pid_allom_h2cd2 + integer, public :: pid_allom_zroot_maxd + integer, public :: pid_allom_zroot_maxz + integer, public :: pid_allom_zroot_mind + integer, public :: pid_allom_zroot_minz + integer, public :: pid_allom_zroot_k + integer, public :: pid_turnover_branch + integer, public :: pid_cnp_nitr_store_ratio + integer, public :: pid_cnp_phos_store_ratio + integer, public :: pid_turnover_leaf_canopy + integer, public :: pid_turnover_leaf_ustory + integer, public :: pid_stoich_nitr + integer, public :: pid_stoich_phos + integer, public :: pid_alloc_organ_priority + integer, public :: pid_cnp_nitr_retrans + integer, public :: pid_cnp_phos_retrans + public :: pstruct + public :: GetParameterIndices public :: Transp2dInt public :: Transp2dReal contains + subroutine GetParameterIndices() + + ! Query the parameter data structure for the names + ! of known parameters, and identify the data-structure + ! index of each, into a named integer constant. These + ! integer constants will be used for retrieval. + ! This assignment happens once, here. + + ! Scalar Parameters + pid_damage_bins = pstruct%GetIndexFromName('fates_history_damage_bin_edges') + + pid_vcmax25top = pstruct%GetIndexFromName('fates_leaf_vcmax25top') + pid_cwd_frac = pstruct%GetIndexFromName('fates_frag_cwd_frac') + pid_damage_frac = pstruct%GetIndexFromName('fates_damage_frac') + + pid_alloc_organ_id = pstruct%GetIndexFromName('fates_alloc_organ_id') + pid_phen_leaf_habit = pstruct%GetIndexFromName('fates_phen_leaf_habit') + pid_phen_stem_drop_frac = pstruct%GetIndexFromName('fates_phen_stem_drop_fraction') + pid_phen_fnrt_drop_frac = pstruct%GetIndexFromName('fates_phen_fnrt_drop_fraction') + pid_phen_mindaysoff = pstruct%GetIndexFromName('fates_phen_mindaysoff') + pid_phen_drought_thresh = pstruct%GetIndexFromName('fates_phen_drought_threshold') + pid_phen_moist_thresh = pstruct%GetIndexFromName('fates_phen_moist_threshold') + pid_leaf_slamax = pstruct%GetIndexFromName('fates_leaf_slamax') + pid_leaf_slatop = pstruct%GetIndexFromName('fates_leaf_slatop') + pid_allom_sai_scaler = pstruct%GetIndexFromName('fates_allom_sai_scaler') + pid_allom_fnrt_prof_a = pstruct%GetIndexFromName('fates_allom_fnrt_prof_a') + pid_allom_fnrt_prof_b = pstruct%GetIndexFromName('fates_allom_fnrt_prof_b') + pid_allom_fnrt_prof_mode = pstruct%GetIndexFromName('fates_allom_fnrt_prof_mode') + pid_woody = pstruct%GetIndexFromName('fates_woody') + pid_wood_density = pstruct%GetIndexFromName('fates_wood_density') + pid_recruit_seed_dbh_thresh = pstruct%GetIndexFromName('fates_recruit_seed_dbh_repro_threshold') + pid_storage_cushion = pstruct%GetIndexFromName('fates_alloc_storage_cushion') + pid_store_priority_frac = pstruct%GetIndexFromName('fates_alloc_store_priority_frac') + pid_turnover_senleaf_fdrought = pstruct%GetIndexFromName('fates_turnover_senleaf_fdrought') + pid_turnover_fnrt = pstruct%GetIndexFromName('fates_turnover_fnrt') + pid_leafn_vert_scaler1 = pstruct%GetIndexFromName('fates_leafn_vert_scaler_coeff1') + pid_leafn_vert_scaler2 = pstruct%GetIndexFromName('fates_leafn_vert_scaler_coeff2') + pid_recruit_seed_alloc_mature = pstruct%GetIndexFromName('fates_recruit_seed_alloc_mature') + pid_recruit_seed_alloc = pstruct%GetIndexFromName('fates_recruit_seed_alloc') + pid_trs_repro_alloc_a = pstruct%GetIndexFromName('fates_trs_repro_alloc_a') + pid_trs_repro_alloc_b = pstruct%GetIndexFromName('fates_trs_repro_alloc_b') + pid_c2b = pstruct%GetIndexFromName('fates_c2b') + pid_grperc = pstruct%GetIndexFromName('fates_grperc') + pid_allom_dbh_maxh = pstruct%GetIndexFromName('fates_allom_dbh_maxheight') + pid_allom_hmode = pstruct%GetIndexFromName('fates_allom_hmode') + pid_allom_lmode = pstruct%GetIndexFromName('fates_allom_lmode') + pid_allom_fmode = pstruct%GetIndexFromName('fates_allom_fmode') + pid_allom_amode = pstruct%GetIndexFromName('fates_allom_amode') + pid_allom_stmode = pstruct%GetIndexFromName('fates_allom_stmode') + pid_allom_cmode = pstruct%GetIndexFromName('fates_allom_cmode') + pid_allom_smode = pstruct%GetIndexFromName('fates_allom_smode') + pid_allom_dmode = pstruct%GetIndexFromName('fates_allom_dmode') + pid_allom_lapersa_int = pstruct%GetIndexFromName('fates_allom_la_per_sa_int') + pid_allom_lapersa_slp = pstruct%GetIndexFromName('fates_allom_la_per_sa_slp') + pid_allom_l2fr = pstruct%GetIndexFromName('fates_allom_l2fr') + pid_cnp_kp = pstruct%GetIndexFromName('fates_cnp_pid_kp') + pid_cnp_ki = pstruct%GetIndexFromName('fates_cnp_pid_ki') + pid_cnp_kd = pstruct%GetIndexFromName('fates_cnp_pid_kd') + pid_cnp_store_ovrflw_frac = pstruct%GetIndexFromName('fates_cnp_store_ovrflw_frac') + pid_cnp_nfix1 = pstruct%GetIndexFromName('fates_cnp_nfix1') + pid_allom_agbfrac = pstruct%GetIndexFromName('fates_allom_agb_frac') + pid_allom_d2h1 = pstruct%GetIndexFromName('fates_allom_d2h1') + pid_allom_d2h2 = pstruct%GetIndexFromName('fates_allom_d2h2') + pid_allom_d2h3 = pstruct%GetIndexFromName('fates_allom_d2h3') + pid_allom_d2l1 = pstruct%GetIndexFromName('fates_allom_d2bl1') + pid_allom_d2l2 = pstruct%GetIndexFromName('fates_allom_d2bl2') + pid_allom_d2l3 = pstruct%GetIndexFromName('fates_allom_d2bl3') + pid_allom_blca_ediff = pstruct%GetIndexFromName('fates_allom_blca_expnt_diff') + pid_allom_d2ca_max = pstruct%GetIndexFromName('fates_allom_d2ca_coefficient_max') + pid_allom_d2ca_min = pstruct%GetIndexFromName('fates_allom_d2ca_coefficient_min') + pid_allom_agb1 = pstruct%GetIndexFromName('fates_allom_agb1') + pid_allom_agb2 = pstruct%GetIndexFromName('fates_allom_agb2') + pid_allom_agb3 = pstruct%GetIndexFromName('fates_allom_agb3') + pid_allom_agb4 = pstruct%GetIndexFromName('fates_allom_agb4') + pid_allom_h2cd1 = pstruct%GetIndexFromName('fates_allom_h2cd1') + pid_allom_h2cd2 = pstruct%GetIndexFromName('fates_allom_h2cd2') + pid_allom_zroot_maxd = pstruct%GetIndexFromName('fates_allom_zroot_max_dbh') + pid_allom_zroot_maxz = pstruct%GetIndexFromName('fates_allom_zroot_max_z') + pid_allom_zroot_mind = pstruct%GetIndexFromName('fates_allom_zroot_min_dbh') + pid_allom_zroot_minz = pstruct%GetIndexFromName('fates_allom_zroot_min_z') + pid_allom_zroot_k = pstruct%GetIndexFromName('fates_allom_zroot_k') + pid_turnover_branch = pstruct%GetIndexFromName('fates_turnover_branch') + pid_cnp_nitr_store_ratio = pstruct%GetIndexFromName('fates_cnp_nitr_store_ratio') + pid_cnp_phos_store_ratio = pstruct%GetIndexFromName('fates_cnp_phos_store_ratio') + pid_turnover_leaf_canopy = pstruct%GetIndexFromName('fates_turnover_leaf_canopy') + pid_turnover_leaf_ustory = pstruct%GetIndexFromName('fates_turnover_leaf_ustory') + pid_stoich_nitr = pstruct%GetIndexFromName('fates_stoich_nitr') + pid_stoich_phos = pstruct%GetIndexFromName('fates_stoich_phos') + pid_alloc_organ_priority = pstruct%GetIndexFromName('fates_alloc_organ_priority') + pid_cnp_nitr_retrans = pstruct%GetIndexFromName('fates_cnp_turnover_nitr_retrans') + pid_cnp_phos_retrans = pstruct%GetIndexFromName('fates_cnp_turnover_phos_retrans') + + end subroutine GetParameterIndices + subroutine CheckParameters + + type(param_type), pointer :: param(:) + integer :: corr_id + real(r8) :: correction + + param => pstruct%parameters + + ! Apply correction to cwd fractions + correction = 1._r8 - sum(param(pid_cwd_frac)%r_data_1d(:),dim=1) + corr_id = maxloc(param(pid_cwd_frac)%r_data_1d(:),dim=1) + param(pid_cwd_frac)%r_data_1d(corr_id) = param(pid_cwd_frac)%r_data_1d(corr_id) + correction + + + end subroutine CheckParameters + + subroutine Transp2dInt(i_2d_in,i_2d_out) ! The FATES JSON parameter files have a legacy from the netcdf diff --git a/main/JSONParameterUtilsMod.F90 b/main/JSONParameterUtilsMod.F90 index 0e43cf4f36..bfc5aec240 100644 --- a/main/JSONParameterUtilsMod.F90 +++ b/main/JSONParameterUtilsMod.F90 @@ -56,15 +56,15 @@ module JSONParameterUtilsMod integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real ! Data types - integer, parameter :: r_scalar_type = 1 - integer, parameter :: r_1d_type = 2 - integer, parameter :: r_2d_type = 3 - integer, parameter :: i_scalar_type = 4 - integer, parameter :: i_1d_type = 5 - integer, parameter :: i_2d_type = 6 - integer, parameter :: c_solo_type = 7 - integer, parameter :: c_1d_type = 8 - integer, parameter :: c_2d_type = 9 + integer, public, parameter :: r_scalar_type = 1 + integer, public, parameter :: r_1d_type = 2 + integer, public, parameter :: r_2d_type = 3 + integer, public, parameter :: i_scalar_type = 4 + integer, public, parameter :: i_1d_type = 5 + integer, public, parameter :: i_2d_type = 6 + integer, public, parameter :: c_solo_type = 7 + integer, public, parameter :: c_1d_type = 8 + integer, public, parameter :: c_2d_type = 9 logical, parameter :: debug = .false. @@ -161,6 +161,7 @@ module JSONParameterUtilsMod contains procedure :: GetDimSizeFromName procedure :: GetParamFromName + procedure :: GetIndexFromName procedure :: ReportAccessCounts procedure :: Destroy end type params_type @@ -1133,26 +1134,38 @@ end function GetDimSizeFromName ! ===================================================================================== - function GetParamFromName(this,param_name) result(param_ptr) - + function GetIndexFromName(this,param_name) result(i) + class(params_type) :: this character(len=*) :: param_name - type(param_type),pointer :: param_ptr integer :: i - nullify(param_ptr) loop_params: do i = 1,size(this%parameters) if(trim(param_name)==this%parameters(i)%name)then - param_ptr=>this%parameters(i) this%parameters(i)%access_count = this%parameters(i)%access_count + 1 return end if end do loop_params - + write(log_unit,*)'Error finding parameter by name,scanned ',size(this%parameters),' parameters' write(log_unit,*)'Cant find: ',trim(param_name) call shr_sys_abort() + end function GetIndexFromName + + ! ============================================================================= + + function GetParamFromName(this,param_name) result(param_ptr) + + class(params_type) :: this + character(len=*) :: param_name + type(param_type),pointer :: param_ptr + integer :: i + + nullify(param_ptr) + i = this%GetIndexFromName(param_name) + param_ptr=>this%parameters(i) + end function GetParamFromName ! ===================================================================================== diff --git a/radiation/TwoStreamMLPEMod.F90 b/radiation/TwoStreamMLPEMod.F90 index 22bf303649..58c7df8f2d 100644 --- a/radiation/TwoStreamMLPEMod.F90 +++ b/radiation/TwoStreamMLPEMod.F90 @@ -910,7 +910,6 @@ subroutine ZenithPrep(this,cosz_in) integer :: ican ! scattering element canopy layer index (top down) integer :: icol ! scattering element column integer :: ib ! band index, matches indexing of rad_params - integer :: ib2 ! band inner loop index while testing for singularity real(r8) :: asu ! single scattering albedo real(r8) :: gdir real(r8) :: tmp0,tmp1,tmp2