diff --git a/cime_config/tests.py b/cime_config/tests.py index df2c581938e5..14a384000615 100644 --- a/cime_config/tests.py +++ b/cime_config/tests.py @@ -49,6 +49,9 @@ "ERS.f19_g16.I20TRGSWCNPECACNTBC.elm-eca_f19_g16_I20TRGSWCNPECACNTBC", "ERS.f19_g16.I20TRGSWCNPRDCTCBC.elm-ctc_f19_g16_I20TRGSWCNPRDCTCBC", "ERS.r05_r05.ICNPRDCTCBC.elm-cbudget", + "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-snowveg_arctic", + "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_default_I1850CNPRDCTCBC", + "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_codetest_I1850CNPRDCTCBC", ) }, @@ -94,8 +97,6 @@ "SMS.r05_r05.IELM.elm-topounit", "ERS.ELM_USRDAT.I1850ELM.elm-usrdat", "ERS.r05_r05.IELM.elm-lnd_rof_2way", - "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_default_I1850CNPRDCTCBC", - "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_codetest_I1850CNPRDCTCBC", "ERS.r05_r05.IELM.elm-V2_ELM_MOSART_features", "ERS.ELM_USRDAT.IELM.elm-surface_water_dynamics" ) diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/shell_commands b/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/shell_commands new file mode 100644 index 000000000000..396d8cea8d66 --- /dev/null +++ b/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/shell_commands @@ -0,0 +1,10 @@ +./xmlchange LND_DOMAIN_FILE=domain_42_FLUXNETSITES_simyr1850_c170912.nc +./xmlchange ATM_DOMAIN_FILE=domain_42_FLUXNETSITES_simyr1850_c170912.nc +./xmlchange LND_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm" +./xmlchange ATM_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm" +./xmlchange DATM_MODE=CLM1PT +./xmlchange DATM_CLMNCEP_YR_START=2000 +./xmlchange DATM_CLMNCEP_YR_END=2000 +./xmlchange ELM_USRDAT_NAME=42_FLUXNETSITES +./xmlchange NTASKS=1 +./xmlchange NTHRDS=1 diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/user_nl_elm new file mode 100644 index 000000000000..bd07678a1345 --- /dev/null +++ b/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/user_nl_elm @@ -0,0 +1,8 @@ +fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_42_FLUXNETSITES_simyr1850_c170912.nc' +! this test has a parameter file with four new parameters defined for NGEE Arctic IM3 +! meant to improve snow-vegetation interactions. +! bendresist -> varies between 0 and 1, represents how much vegetation bends under snow loading +! vegshape -> suggested value 1 (parabolic), but 2 also used previously (Liston and Heimstra, 2011) +! taper -> deadstem height/radius ratio, moved from hardcoded in VegStructUpdateMod to pftvarcon +! stocking -> individual density on landscape (plants/m2), moved from hardcoded in VegStructUpdateMod to pftvarcon +paramfile = '$DIN_LOC_ROOT/lnd/clm2/paramdata/clm_params_ngeea-im3_c240822.nc' diff --git a/components/elm/src/biogeochem/AllocationMod.F90 b/components/elm/src/biogeochem/AllocationMod.F90 index 141194397fb0..5c315fc1769b 100644 --- a/components/elm/src/biogeochem/AllocationMod.F90 +++ b/components/elm/src/biogeochem/AllocationMod.F90 @@ -36,11 +36,11 @@ module AllocationMod use elm_varctl , only : NFIX_PTASE_plant use ELMFatesInterfaceMod , only : hlm_fates_interface_type use elm_varctl , only: iulog - use elm_varctl , only : carbon_only - use elm_varctl , only : carbonnitrogen_only + use elm_varctl , only : carbon_only + use elm_varctl , only : carbonnitrogen_only use elm_varctl , only : carbonphosphorus_only use shr_infnan_mod , only: nan => shr_infnan_nan, assignment(=) - + ! implicit none save @@ -114,23 +114,23 @@ module AllocationMod ! to toggle and update which processes are active. ! This will get set to false ! after ad_carbon_only is complete. - - + + logical :: crop_supln = .false. !Prognostic crop receives supplemental Nitrogen - + real(r8), allocatable,target :: veg_rootc_bigleaf(:,:) ! column-level fine-root biomas kgc/m3 integer, pointer :: ft_index_bigleaf(:) ! array holding the pft index of each competitor ! ECA parameters - ! scaling factor for plant fine root biomass to calculate nutrient carrier enzyme abundance - real(r8), parameter :: e_plant_scalar = 0.0000125_r8 - - ! scaling factor for plant fine root biomass to calculate nutrient carrier enzyme abundance - real(r8), parameter :: e_decomp_scalar = 0.05_r8 + ! scaling factor for plant fine root biomass to calculate nutrient carrier enzyme abundance + real(r8), parameter :: e_plant_scalar = 0.0000125_r8 + + ! scaling factor for plant fine root biomass to calculate nutrient carrier enzyme abundance + real(r8), parameter :: e_decomp_scalar = 0.05_r8 !$acc declare create(e_decomp_scalar) !$acc declare create(e_plant_scalar) - + !$acc declare copyin(crop_supln) !----------------------------------------------------------------------- @@ -213,8 +213,8 @@ subroutine AllocationInit ( bounds, elm_fates) use elm_time_manager, only: get_step_size use elm_varpar , only: crop_prog use elm_varctl , only: iulog - use elm_varctl , only : carbon_only - use elm_varctl , only : carbonnitrogen_only + use elm_varctl , only : carbon_only + use elm_varctl , only : carbonnitrogen_only use elm_varctl , only : carbonphosphorus_only @@ -234,9 +234,9 @@ subroutine AllocationInit ( bounds, elm_fates) integer :: max_comps ! maximum number of possible plant competitors ! elm big-leaf: number of pfts/patches ! fates: number of cohorts in the column - - + + !----------------------------------------------------------------------- if ( crop_prog )then @@ -253,7 +253,7 @@ subroutine AllocationInit ( bounds, elm_fates) end if end if - + ! set time steps dt = real( get_step_size(), r8 ) @@ -374,7 +374,7 @@ subroutine EvaluateSupplStatus() end subroutine EvaluateSupplStatus !------------------------------------------------------------------------------------------------- - + subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & photosyns_vars, crop_vars, canopystate_vars, cnstate_vars, dt, yr) ! PHASE-1 of Allocation: loop over patches to assess the total plant N demand and P demand @@ -430,7 +430,7 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) froot_leaf => veg_vp%froot_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new fine root C per new leaf C (gC/gC) croot_stem => veg_vp%croot_stem , & ! Input: [real(r8) (:) ] allocation parameter: new coarse root C per new stem C (gC/gC) stem_leaf => veg_vp%stem_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new stem c per new leaf C (gC/gC) @@ -971,7 +971,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & ! Fractional uptake profiles, that are proportional to root density real(r8):: nuptake_prof(bounds%begc:bounds%endc,1:nlevdecomp) real(r8):: puptake_prof(bounds%begc:bounds%endc,1:nlevdecomp) - integer, allocatable :: filter_pcomp(:) ! this is a plant competitor map for FATES/ELM-BL w/ ECA + integer, allocatable :: filter_pcomp(:) ! this is a plant competitor map for FATES/ELM-BL w/ ECA real(r8), allocatable,target :: plant_nh4demand_vr_fates(:,:) ! nh4 demand per competitor per soil layer real(r8), allocatable,target :: plant_no3demand_vr_fates(:,:) ! no3 demand per competitor per soil layer real(r8), allocatable,target :: plant_pdemand_vr_fates(:,:) ! p demand per competitor per soil layer @@ -990,7 +990,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & real(r8):: cp_stoich_var=0.4 ! variability of CP ratio - + !----------------------------------------------------------------------- associate( & @@ -1152,7 +1152,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & n_pcomp = elm_fates%fates(ci)%bc_out(s)%num_plant_comps pci = 1 pcf = n_pcomp - + if( nu_com.eq.'RD') then ! Overwrite the column level demands, since fates plants are all sharing @@ -1160,13 +1160,13 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & ! to scale up to column plant_ndemand_col(c) = 0._r8 plant_pdemand_col(c) = 0._r8 - + ! We fill the vertically resolved array to simplify some jointly used code do j = 1, nlevdecomp col_plant_ndemand_vr(c,j) = 0._r8 col_plant_pdemand_vr(c,j) = 0._r8 - + do f = 1,n_pcomp ft = elm_fates%fates(ci)%bc_out(s)%ft_index(f) @@ -1184,10 +1184,10 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & ! [gN/m2/s] plant_ndemand_col(c) = plant_ndemand_col(c) + col_plant_ndemand_vr(c,j)*dzsoi_decomp(j) plant_pdemand_col(c) = plant_pdemand_col(c) + col_plant_pdemand_vr(c,j)*dzsoi_decomp(j) - + end do - + else !(ECA) do f = 1,n_pcomp @@ -1195,7 +1195,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & end do veg_rootc_ptr => elm_fates%fates(ci)%bc_out(s)%veg_rootc - ft_index_ptr => elm_fates%fates(ci)%bc_out(s)%ft_index ! Should be + ft_index_ptr => elm_fates%fates(ci)%bc_out(s)%ft_index ! Should be decompmicc(:) = elm_fates%fates(ci)%bc_out(s)%decompmicc(:) ! Should be (nlevdecomp) cn_scalar_runmean_ptr => elm_fates%fates(ci)%bc_out(s)%cn_scalar ! This is 1.0 @@ -1216,7 +1216,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & pci = col_pp%pfti(c) ! Initial plant competitor index pcf = col_pp%pftf(c) ! Final plant competitor index - + if (nu_com .eq. 'RD') then do j = 1, nlevdecomp @@ -1225,7 +1225,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & end do else - + f = 0 decompmicc(:) = 0._r8 do p = col_pp%pfti(c), col_pp%pftf(c) @@ -1246,7 +1246,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & end if end do n_pcomp = f - + ft_index_ptr => ft_index_bigleaf veg_rootc_ptr => veg_rootc_bigleaf @@ -1263,7 +1263,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & (leafcn(ivt(p)) - leafcn(ivt(p))*(1- cn_stoich_var)),0.0_r8),1.0_r8) end do end if - + km_nh4_ptr => km_plant_nh4 vmax_nh4_ptr => vmax_plant_nh4 cn_scalar_runmean_ptr => cn_scalar_runmean @@ -1298,14 +1298,14 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & ! (1) add nitrogen and phosphorus competition ! (2) nitrogen and phosphorus uptake is based on root kinetics ! (3) no second pass nutrient uptake for plants - ! ============================================================= - + ! ============================================================= + if (nu_com .eq. 'RD') then ! Estimate actual allocation rates via Relative Demand ! approach (RD) - + call NAllocationRD(col_plant_ndemand_vr(c,:), & ! IN potential_immob_vr(c,:), & ! IN AllocParamsInst%compet_plant_nh4, & ! IN @@ -1332,19 +1332,19 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & ! Estimate actual allocation rates via Capacitance Aquisition ! approach (ECA/CA) - + call NAllocationECAMIC(pci,dt, & ! IN bd(c,:), & ! IN h2osoi_vol(c,:), & ! IN t_scalar(c,:), & ! IN - n_pcomp, & ! IN + n_pcomp, & ! IN filter_pcomp(1:n_pcomp), & ! IN veg_rootc_ptr(pci:pcf,:), & ! IN ft_index_ptr(pci:pcf), & ! IN cn_scalar_runmean_ptr(pci:pcf), & ! IN decompmicc, & ! IN smin_nh4_vr(c,:), & ! IN - nu_com, & ! IN + nu_com, & ! IN km_nh4_ptr, & ! IN vmax_nh4_ptr, & ! IN km_decomp_nh4, & ! IN @@ -1389,7 +1389,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & ! NO3 flux demands. supplement_to_sminn_vr(c,j) = 0._r8 if (carbon_only .or. carbonphosphorus_only) then - + if ( fpi_no3_vr(j) + fpi_nh4_vr(j) < 1._r8 ) then fpi_vr(c,j) = 1._r8 fpi_nh4_vr(j) = 1.0_r8 - fpi_no3_vr(j) @@ -1409,9 +1409,9 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & supplement_to_sminn_vr(c,j) = supplement_to_sminn_vr(c,j) + col_plant_ndemand_vr(c,j) smin_nh4_to_plant_vr(c,j) = col_plant_ndemand_vr(c,j) - smin_no3_to_plant_vr(c,j) end if - - + + end if ! sum up nitrogen limitation to decomposition @@ -1422,14 +1422,14 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & actual_immob_vr(c,j) = actual_immob_no3_vr(c,j) + actual_immob_nh4_vr(c,j) end do - + ! Starting resolving P limitation !!! ! ============================================================= if (nu_com .eq. 'RD') then ! Relative Demand (RD) - + call PAllocationRD(col_plant_pdemand_vr(c,:), & ! IN potential_immob_p_vr(c,:), & ! IN solutionp_vr(c,:), & ! IN @@ -1438,38 +1438,38 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & actual_immob_p_vr(c,:), & ! OUT sminp_to_plant_vr(c,:), & ! OUT supplement_to_sminp_vr(c,:)) ! OUT - + else call PAllocationECAMIC(pci,dt, & ! IN h2osoi_vol(c,:), & ! IN - t_scalar(c,:), & ! IN - gross_pmin_vr(c,:), & ! IN - potential_immob_p_vr(c,:), & ! IN - biochem_pmin_vr_col(c,:), & ! IN - primp_to_labilep_vr_col(c,:), & ! IN - pdep_to_sminp(c), & ! IN - pdep_prof(c,:), & ! IN - vmax_minsurf_p_vr(isoilorder(c),:), & ! IN - km_minsurf_p_vr(isoilorder(c),:), & ! IN - solutionp_vr(c,:), & ! IN - nu_com, & ! IN - n_pcomp, & ! IN + t_scalar(c,:), & ! IN + gross_pmin_vr(c,:), & ! IN + potential_immob_p_vr(c,:), & ! IN + biochem_pmin_vr_col(c,:), & ! IN + primp_to_labilep_vr_col(c,:), & ! IN + pdep_to_sminp(c), & ! IN + pdep_prof(c,:), & ! IN + vmax_minsurf_p_vr(isoilorder(c),:), & ! IN + km_minsurf_p_vr(isoilorder(c),:), & ! IN + solutionp_vr(c,:), & ! IN + nu_com, & ! IN + n_pcomp, & ! IN filter_pcomp(1:n_pcomp), & ! IN - veg_rootc_ptr(pci:pcf,:), & ! IN - ft_index_ptr(pci:pcf), & ! IN - decompmicc, & ! IN - cp_scalar_runmean_ptr(pci:pcf), & ! IN - km_p_ptr(:), & ! IN - vmax_p_ptr(:), & ! IN - km_decomp_p, & ! IN - labilep_vr(c,:), & ! IN + veg_rootc_ptr(pci:pcf,:), & ! IN + ft_index_ptr(pci:pcf), & ! IN + decompmicc, & ! IN + cp_scalar_runmean_ptr(pci:pcf), & ! IN + km_p_ptr(:), & ! IN + vmax_p_ptr(:), & ! IN + km_decomp_p, & ! IN + labilep_vr(c,:), & ! IN plant_pdemand_vr_ptr(pci:pcf,:), & ! INOUT - col_plant_pdemand_vr(c,:), & ! OUT + col_plant_pdemand_vr(c,:), & ! OUT adsorb_to_labilep_vr(c,:), & ! OUT fpi_p_vr(c,:), & ! OUT actual_immob_p_vr(c,:), & ! OUT - sminp_to_plant_vr(c,:), & ! OUT + sminp_to_plant_vr(c,:), & ! OUT desorb_to_solutionp_vr(c,:), & ! OUT supplement_to_sminp_vr(c,:)) ! OUT @@ -1482,13 +1482,13 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & supplement_to_sminp_vr(c,j) = col_plant_pdemand_vr(c,j) end do end if - + end if ! end of P competition ! resolving N limitation vs. P limitation for decomposition ! update (1) actual immobilization for N and P (2) sminn_to_plant and sminp_to_plant ! We only resolve co-limitations when are supplementing neither element - + np_bothactive: if ( .not.carbon_only .and. & .not.carbonphosphorus_only .and. & .not.carbonnitrogen_only ) then @@ -1593,7 +1593,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & actual_immob_vr(c,j) = potential_immob_vr(c,j) * fpi_p_vr(c,j) end do end if - + ! sum up plant N/P uptake at column level and patch level ! sum up N fluxes to plant after initial competition sminn_to_plant(c) = 0._r8 @@ -1602,7 +1602,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & sminn_to_plant(c) = sminn_to_plant(c) + sminn_to_plant_vr(c,j) * dzsoi_decomp(j) sminp_to_plant(c) = sminp_to_plant(c) + sminp_to_plant_vr(c,j) * dzsoi_decomp(j) end do - + ! update column plant N/P demand, pft level plant NP uptake for ECA and MIC mode eca_filter: if (nu_com .eq. 'ECA' .or. nu_com .eq. 'MIC') then @@ -1648,14 +1648,14 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & end if end do end if - + end do end if eca_filter end do col_loop - + if ((nu_com .eq. 'ECA' .or. nu_com .eq. 'MIC')) then deallocate(filter_pcomp) if(.not.use_fates)then @@ -1750,7 +1750,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & ! Set the FATES N and P uptake fluxes - + if(use_fates)then do fc=1,num_soilc c = filter_soilc(fc) @@ -1771,15 +1771,15 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & (elm_fates%fates(ci)%bc_pconst%vmax_nh4(ft)+elm_fates%fates(ci)%bc_pconst%vmax_no3(ft)) * & dzsoi_decomp(j) end do - + do j = 1,nlevdecomp - elm_fates%fates(ci)%bc_in(s)%plant_nh4_uptake_flux(f,1) = & + elm_fates%fates(ci)%bc_in(s)%plant_nh4_uptake_flux(f,1) = & elm_fates%fates(ci)%bc_in(s)%plant_nh4_uptake_flux(f,1) + & smin_nh4_to_plant_vr(c,j)*dt*dzsoi_decomp(j) * & (ndemand/plant_ndemand_col(c)) - elm_fates%fates(ci)%bc_in(s)%plant_no3_uptake_flux(f,1) = & + elm_fates%fates(ci)%bc_in(s)%plant_no3_uptake_flux(f,1) = & elm_fates%fates(ci)%bc_in(s)%plant_no3_uptake_flux(f,1) + & smin_no3_to_plant_vr(c,j)*dt*dzsoi_decomp(j) * & (ndemand/plant_ndemand_col(c)) @@ -1787,12 +1787,12 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & end do end do end if - + if( plant_pdemand_col(c)>tiny(plant_pdemand_col(c)) ) then do f = 1,n_pcomp ft = elm_fates%fates(ci)%bc_out(s)%ft_index(f) - + pdemand=0._r8 do j = 1,nlevdecomp ! [gP/m2/s] @@ -1800,14 +1800,14 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & elm_fates%fates(ci)%bc_pconst%vmax_p(ft) * & dzsoi_decomp(j) end do - + do j = 1,nlevdecomp ! [gP/m2/step] - elm_fates%fates(ci)%bc_in(s)%plant_p_uptake_flux(f,1) = & + elm_fates%fates(ci)%bc_in(s)%plant_p_uptake_flux(f,1) = & elm_fates%fates(ci)%bc_in(s)%plant_p_uptake_flux(f,1) + & sminp_to_plant_vr(c,j)*dt*dzsoi_decomp(j) * & (pdemand/plant_pdemand_col(c)) - + end do end do end if @@ -1817,21 +1817,21 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & do f = 1,n_pcomp do j = 1,nlevdecomp - elm_fates%fates(ci)%bc_in(s)%plant_nh4_uptake_flux(f,1) = & - elm_fates%fates(ci)%bc_in(s)%plant_nh4_uptake_flux(f,1) + & + elm_fates%fates(ci)%bc_in(s)%plant_nh4_uptake_flux(f,1) = & + elm_fates%fates(ci)%bc_in(s)%plant_nh4_uptake_flux(f,1) + & plant_nh4demand_vr_fates(f,j) * fpg_nh4_vr(c,j) * dzsoi_decomp(j) * dt - - elm_fates%fates(ci)%bc_in(s)%plant_no3_uptake_flux(f,1) = & - elm_fates%fates(ci)%bc_in(s)%plant_no3_uptake_flux(f,1) + & + + elm_fates%fates(ci)%bc_in(s)%plant_no3_uptake_flux(f,1) = & + elm_fates%fates(ci)%bc_in(s)%plant_no3_uptake_flux(f,1) + & plant_no3demand_vr_fates(f,j) * fpg_no3_vr(c,j) * dzsoi_decomp(j) * dt - - elm_fates%fates(ci)%bc_in(s)%plant_p_uptake_flux(f,1) = & - elm_fates%fates(ci)%bc_in(s)%plant_p_uptake_flux(f,1) + & + + elm_fates%fates(ci)%bc_in(s)%plant_p_uptake_flux(f,1) = & + elm_fates%fates(ci)%bc_in(s)%plant_p_uptake_flux(f,1) + & (plant_pdemand_vr_fates(f,j) * fpg_p_vr(c,j)) * dzsoi_decomp(j) * dt - + end do end do - + end if end do @@ -1840,7 +1840,7 @@ subroutine Allocation2_ResolveNPLimit (bounds, num_soilc, filter_soilc , & deallocate(plant_no3demand_vr_fates) deallocate(plant_pdemand_vr_fates) end if - + end if ! if(use_fates) end associate @@ -1907,7 +1907,7 @@ subroutine Allocation3_PlantCNPAlloc (bounds , & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) froot_leaf => veg_vp%froot_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new fine root C per new leaf C (gC/gC) croot_stem => veg_vp%croot_stem , & ! Input: [real(r8) (:) ] allocation parameter: new coarse root C per new stem C (gC/gC) stem_leaf => veg_vp%stem_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new stem c per new leaf C (gC/gC) @@ -2959,7 +2959,7 @@ subroutine Allocation3_PlantCNPAlloc (bounds , & end associate end subroutine Allocation3_PlantCNPAlloc - + ! ====================================================================================== subroutine NAllocationECAMIC(pci,dt, & ! IN @@ -2976,14 +2976,14 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN nu_com, & ! IN km_nh4_plant, & ! IN (pft) vmax_nh4_plant, & ! IN (pft) - km_decomp_nh4, & ! IN + km_decomp_nh4, & ! IN potential_immob_vr, & ! IN (j) plant_nh4demand_vr, & ! INOUT (i,j) col_plant_nh4demand_vr, & ! OUT (j) fpi_nh4_vr, & ! OUT (j) actual_immob_nh4_vr, & ! OUT (j) smin_nh4_to_plant_vr, & ! OUT (j) - smin_no3_vr, & ! IN (j) + smin_no3_vr, & ! IN (j) km_no3_plant, & ! IN (pft) vmax_no3_plant, & ! IN (pft) km_decomp_no3, & ! IN (j) @@ -3004,12 +3004,12 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN ! kinetics following Zhu et al., 2016 DOI: 10.1002/2016JG003554 ! ------------------------------------------------------------------------------------ use elm_varpar , only: nlevdecomp - + integer, intent(in) :: pci ! First index of plant comp arrays real(r8), intent(in) :: dt ! Time step duration [s] real(r8), intent(in) :: bd(:) ! Bulk density of dry soil material [kg m-3] real(r8), intent(in) :: h2osoi_vol(:) ! Vol. Soil Water in each layer [m3] - real(r8), intent(in) :: t_scalar(:) ! fraction by which decomposition is limited by temperature + real(r8), intent(in) :: t_scalar(:) ! fraction by which decomposition is limited by temperature integer, intent(in) :: n_pcomp ! number of plant competitors integer, intent(in) :: filter_pcomp(:) ! plant competition filter real(r8), intent(in) :: veg_rootc(pci:,:) ! total fine-root biomass of each competitor [gC/m3] @@ -3065,7 +3065,7 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN integer :: i,ip ! loop index for competitors integer :: ft ! loop index for pfts - ! 2.76 consider soil adsorption effect on [NH4+] availability, + ! 2.76 consider soil adsorption effect on [NH4+] availability, ! based on Zhu et al., 2016 DOI: 10.1002/2016JG003554 real(r8), parameter :: adsorp_nh4_eff = 2.76_r8 @@ -3073,14 +3073,14 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN do j = 1, nlevdecomp - ! Plant, microbial decomposers compete for NH4. Thus loop over each + ! Plant, microbial decomposers compete for NH4. Thus loop over each ! plant competitor in this competitive space (column). - ! Calculate competition coefficients for N/P, first need to convert - ! concentration to per soil water based + ! Calculate competition coefficients for N/P, first need to convert + ! concentration to per soil water based ! concentration of mineralized nutrient, per soil water solution_conc = smin_nh4_vr(j) / (bd(j)*adsorp_nh4_eff*m3_per_liter + h2osoi_vol(j)) - + e_km = 0._r8 do i = 1, n_pcomp ip = filter_pcomp(i) @@ -3089,14 +3089,14 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN end do e_km = e_km + e_decomp_scalar*decompmicc(j)*(1._r8/km_decomp_nh4 + 1._r8/km_nit) - + do i = 1, n_pcomp ip = filter_pcomp(i) ft = ft_index(ip) - compet_plant(i) = solution_conc / & + compet_plant(i) = solution_conc / & ( km_nh4_plant(ft) * (1._r8 + solution_conc/km_nh4_plant(ft) + e_km)) end do - + compet_decomp = solution_conc / (km_decomp_nh4 * (1._r8 + solution_conc/km_decomp_nh4 + e_km)) compet_nit = solution_conc / (km_nit * (1._r8 + solution_conc/km_nit + e_km)) @@ -3110,7 +3110,7 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN ip = filter_pcomp(i) ft = ft_index(ip) - ! This is the demand per m3 of the column (not patch) + ! This is the demand per m3 of the column (not patch) ! (for native ELM divide through by the patch weight to get per m3 of patch) plant_nh4demand_vr(ip,j) = max(0._r8,vmax_nh4_plant(ft) * veg_rootc(ip,j) * & cn_scalar_runmean(ip) * t_scalar(j) * compet_plant(i)) @@ -3130,7 +3130,7 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN else ! 'MIC' mode - sum_nh4_demand_scaled = potential_immob_vr(j)*compet_decomp + & + sum_nh4_demand_scaled = potential_immob_vr(j)*compet_decomp + & pot_f_nit_vr(j)*compet_nit end if @@ -3199,10 +3199,10 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN do i = 1, n_pcomp ip = filter_pcomp(i) ft = ft_index(ip) - compet_plant(i) = solution_conc / & + compet_plant(i) = solution_conc / & ( km_no3_plant(ft) * (1._r8 + solution_conc/km_no3_plant(ft) + e_km)) end do - + compet_decomp = solution_conc / (km_decomp_no3 * (1._r8 + solution_conc/km_decomp_no3 + e_km)) compet_denit = solution_conc / (km_den * (1._r8 + solution_conc/km_den + e_km)) @@ -3215,7 +3215,7 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN ip = filter_pcomp(i) ft = ft_index(ip) - ! This is the demand per m3 of the column (not patch) + ! This is the demand per m3 of the column (not patch) ! (for native ELM divide through by the patch weight to get per m3 of patch) plant_no3demand_vr(ip,j) = max(0._r8,vmax_no3_plant(ft) * veg_rootc(ip,j) * & cn_scalar_runmean(ip) * t_scalar(j) * compet_plant(i)) @@ -3245,7 +3245,7 @@ subroutine NAllocationECAMIC(pci,dt, & ! IN smin_no3_to_plant_vr(j) = col_plant_no3demand_vr(j) f_denit_vr(j) = pot_f_denit_vr(j) - else + else ! NO3 availability can not satisfy the sum of immobilization, denitrification, and ! plant growth demands, so these three demands compete for available @@ -3284,34 +3284,34 @@ end subroutine NAllocationECAMIC subroutine PAllocationECAMIC(pci, & dt, & - h2osoi_vol, & - t_scalar, & - gross_pmin_vr, & - potential_immob_p_vr, & - biochem_pmin_vr_col, & - primp_to_labilep_vr_col, & - pdep_to_sminp, & - pdep_prof, & + h2osoi_vol, & + t_scalar, & + gross_pmin_vr, & + potential_immob_p_vr, & + biochem_pmin_vr_col, & + primp_to_labilep_vr_col, & + pdep_to_sminp, & + pdep_prof, & vmax_minsurf_p_vr, & km_minsurf_p_vr, & solutionp_vr, & nu_com, & n_pcomp, & - filter_pcomp, & - veg_rootc, & + filter_pcomp, & + veg_rootc, & ft_index, & decompmicc, & - cp_scalar_runmean, & + cp_scalar_runmean, & km_plant_p, & vmax_plant_p, & - km_decomp_p, & + km_decomp_p, & labilep_vr, & - plant_pdemand_vr_patch, & - col_plant_pdemand_vr, & + plant_pdemand_vr_patch, & + col_plant_pdemand_vr, & adsorb_to_labilep_vr, & fpi_p_vr, & actual_immob_p_vr, & - sminp_to_plant_vr, & + sminp_to_plant_vr, & desorb_to_solutionp_vr, & supplement_to_sminp_vr) @@ -3370,7 +3370,7 @@ subroutine PAllocationECAMIC(pci, & ! plant P uptake, microbial P uptake/release ! secondary P desorption is assumed to go into solution P pool - do j = 1, nlevdecomp + do j = 1, nlevdecomp ! plant, microbial decomposer, mineral surface compete for P ! loop over each pft within the same column @@ -3390,27 +3390,27 @@ subroutine PAllocationECAMIC(pci, & do i = 1,n_pcomp ip = filter_pcomp(i) ft = ft_index(ip) - compet_plant(i) = solution_pconc / & + compet_plant(i) = solution_pconc / & (km_plant_p(ft)*(1._r8 + solution_pconc/km_plant_p(ft) + e_km_p)) end do - + compet_decomp_p = solution_pconc / & (km_decomp_p * (1._r8 + solution_pconc/km_decomp_p + e_km_p)) - compet_minsurf_p = solution_pconc/ & + compet_minsurf_p = solution_pconc/ & (km_minsurf_p_vr(j) * (1._r8 + solution_pconc/km_minsurf_p_vr(j) + e_km_p)) col_plant_pdemand_vr(j) = 0._r8 do i = 1,n_pcomp ip = filter_pcomp(i) ft = ft_index(ip) - plant_pdemand_vr_patch(ip,j) = max(0._r8,vmax_plant_p(ft) * veg_rootc(ip,j) * & + plant_pdemand_vr_patch(ip,j) = max(0._r8,vmax_plant_p(ft) * veg_rootc(ip,j) * & cp_scalar_runmean(ip) * t_scalar(j) * compet_plant(i)) col_plant_pdemand_vr(j) = col_plant_pdemand_vr(j) + plant_pdemand_vr_patch(ip,j) end do ! potential adsorption rate without plant and microbial interaction - ! including weathering, deposition, phosphatase, mineralization, + ! including weathering, deposition, phosphatase, mineralization, ! immobilization, plant uptake dsolutionp_dt = gross_pmin_vr(j) -potential_immob_p_vr(j) - & col_plant_pdemand_vr(j) + biochem_pmin_vr_col(j) + & @@ -3477,7 +3477,7 @@ subroutine PAllocationECAMIC(pci, & if (nu_com .eq. 'MIC') sminp_to_plant_vr(j) = min(max( 0._r8, & (solutionp_vr(j)/dt) - actual_immob_p_vr(j) - adsorb_to_labilep_vr(j) ), & - col_plant_pdemand_vr(j)) + col_plant_pdemand_vr(j)) end if end do @@ -3488,7 +3488,7 @@ end subroutine PAllocationECAMIC subroutine NAllocationRD(col_plant_ndemand_vr, &! IN (j) potential_immob_vr, & ! IN (j) - compet_plants_nh4, & ! IN + compet_plants_nh4, & ! IN compet_decomp_nh4, & ! IN dt, & ! IN smin_nh4_vr, & ! IN (j) @@ -3641,7 +3641,7 @@ end subroutine NAllocationRD ! ====================================================================================== - subroutine PAllocationRD(col_plant_pdemand_vr, & ! IN + subroutine PAllocationRD(col_plant_pdemand_vr, & ! IN potential_immob_p_vr, & ! IN (j) solutionp_vr, & ! IN (j) dt, & ! IN @@ -3684,7 +3684,7 @@ subroutine PAllocationRD(col_plant_pdemand_vr, & ! IN actual_immob_p_vr(j) = potential_immob_p_vr(j) sminp_to_plant_vr(j) = col_plant_pdemand_vr(j) supplement_to_sminp_vr(j) = sum_pdemand - (solutionp_vr(j)/dt) - + else ! P availability can not satisfy the sum of immobilization and ! plant growth demands, so these two demands compete for @@ -3702,7 +3702,7 @@ subroutine PAllocationRD(col_plant_pdemand_vr, & ! IN fpi_p_vr(j) = 0.0_r8 end if - sminp_to_plant_vr(j) = max( 0._r8,(solutionp_vr(j)/dt) - actual_immob_p_vr(j) ) + sminp_to_plant_vr(j) = max( 0._r8,(solutionp_vr(j)/dt) - actual_immob_p_vr(j) ) end if end do @@ -3749,7 +3749,7 @@ subroutine calc_nuptake_prof(bounds, num_soilc, filter_soilc, cnstate_vars, nupt c = filter_soilc(fc) sminn_vr_loc(c,j) = smin_no3_vr(c,j) + smin_nh4_vr(c,j) - + if(use_pflotran .and. pf_cmode) then sminn_tot(c) = sminn_tot(c) + sminn_vr_loc(c,j) * dzsoi_decomp(j) & *(nfixation_prof(c,j)*dzsoi_decomp(j)) ! weighted by froot fractions in annual max. active layers diff --git a/components/elm/src/biogeochem/CNAllocationBetrMod.F90 b/components/elm/src/biogeochem/CNAllocationBetrMod.F90 index 8cfc342bb404..32f9fc6eb23b 100644 --- a/components/elm/src/biogeochem/CNAllocationBetrMod.F90 +++ b/components/elm/src/biogeochem/CNAllocationBetrMod.F90 @@ -26,13 +26,13 @@ module CNAllocationBeTRMod use PhotosynthesisType , only : photosyns_type use CropType , only : crop_type use VegetationPropertiesType, only : veg_vp - use LandunitType , only : lun_pp + use LandunitType , only : lun_pp use ColumnType , only : col_pp - use ColumnDataType , only : col_cf, col_ns, col_nf, col_ps, col_pf + use ColumnDataType , only : col_cf, col_ns, col_nf, col_ps, col_pf use VegetationType , only : veg_pp use VegetationDataType , only : veg_cs, veg_ns, veg_nf, veg_ps, veg_pf - use VegetationDataType , only : veg_cf, c13_veg_cf, c14_veg_cf - + use VegetationDataType , only : veg_cf, c13_veg_cf, c14_veg_cf + ! bgc interface & pflotran module switches use elm_varctl , only : nu_com use SoilStatetype , only : soilstate_type @@ -359,7 +359,7 @@ subroutine Allocation1_PlantNPDemand (bounds, num_soilc, filter_soilc, num_soilp associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) froot_leaf => veg_vp%froot_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new fine root C per new leaf C (gC/gC) croot_stem => veg_vp%croot_stem , & ! Input: [real(r8) (:) ] allocation parameter: new coarse root C per new stem C (gC/gC) stem_leaf => veg_vp%stem_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new stem c per new leaf C (gC/gC) @@ -1167,7 +1167,7 @@ subroutine Allocation3_PlantCNPAlloc (bounds , & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type ! - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) froot_leaf => veg_vp%froot_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new fine root C per new leaf C (gC/gC) croot_stem => veg_vp%croot_stem , & ! Input: [real(r8) (:) ] allocation parameter: new coarse root C per new stem C (gC/gC) stem_leaf => veg_vp%stem_leaf , & ! Input: [real(r8) (:) ] allocation parameter: new stem c per new leaf C (gC/gC) diff --git a/components/elm/src/biogeochem/CNGapMortalityBeTRMod.F90 b/components/elm/src/biogeochem/CNGapMortalityBeTRMod.F90 index 3bbbf1db1f60..ed8811f61f89 100644 --- a/components/elm/src/biogeochem/CNGapMortalityBeTRMod.F90 +++ b/components/elm/src/biogeochem/CNGapMortalityBeTRMod.F90 @@ -18,9 +18,9 @@ module CNGapMortalityBeTRMod use ColumnType , only : col_pp use ColumnDataType , only : col_cf, col_nf, col_pf use VegetationPropertiesType , only : veg_vp - use VegetationType , only : veg_pp - use VegetationDataType , only : veg_cs, veg_cf, veg_ns, veg_nf - use VegetationDataType , only : veg_ps, veg_pf + use VegetationType , only : veg_pp + use VegetationDataType , only : veg_cs, veg_cf, veg_ns, veg_nf + use VegetationDataType , only : veg_ps, veg_pf use PhosphorusFluxType , only : phosphorusflux_type use PhosphorusStateType , only : phosphorusstate_type @@ -121,7 +121,7 @@ subroutine CNGapMortality (& associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody & ! Input: [real(r8) (:) ] binary flag for woody lifeform + woody => veg_vp%woody & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) ) ! set the mortality rate based on annual rate diff --git a/components/elm/src/biogeochem/CNNStateUpdate1BeTRMod.F90 b/components/elm/src/biogeochem/CNNStateUpdate1BeTRMod.F90 index 465c92f82b4c..b785f52549e7 100644 --- a/components/elm/src/biogeochem/CNNStateUpdate1BeTRMod.F90 +++ b/components/elm/src/biogeochem/CNNStateUpdate1BeTRMod.F90 @@ -59,7 +59,7 @@ subroutine NStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp, & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Input: [integer (:) ] which pool is C taken from for a given decomposition step cascade_receiver_pool => decomp_cascade_con%cascade_receiver_pool , & ! Input: [integer (:) ] which pool is C added to for a given decomposition step diff --git a/components/elm/src/biogeochem/CNPhenologyBeTRMod.F90 b/components/elm/src/biogeochem/CNPhenologyBeTRMod.F90 index ac3b58538e8f..066dcfb62c82 100644 --- a/components/elm/src/biogeochem/CNPhenologyBeTRMod.F90 +++ b/components/elm/src/biogeochem/CNPhenologyBeTRMod.F90 @@ -525,7 +525,7 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp , & prev_dayl => grc_pp%prev_dayl , & ! Input: [real(r8) (:) ] daylength from previous time step (s) season_decid => veg_vp%season_decid , & ! Input: [real(r8) (:) ] binary flag for seasonal-deciduous leaf habit (0 or 1) - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd) @@ -875,7 +875,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , & dayl => grc_pp%dayl , & ! Input: [real(r8) (:) ] daylength (s) leaf_long => veg_vp%leaf_long , & ! Input: [real(r8) (:) ] leaf longevity (yrs) - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) stress_decid => veg_vp%stress_decid , & ! Input: [real(r8) (:) ] binary flag for stress-deciduous leaf habit (0 or 1) soilpsi => soilstate_vars%soilpsi_col , & ! Input: [real(r8) (:,:) ] soil water potential in each soil layer (MPa) @@ -2295,7 +2295,7 @@ subroutine CNOnsetGrowth (num_soilp, filter_soilp, & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) onset_flag => cnstate_vars%onset_flag_patch , & ! Input: [real(r8) (:) ] onset flag onset_counter => cnstate_vars%onset_counter_patch , & ! Input: [real(r8) (:) ] onset days counter @@ -2880,7 +2880,7 @@ subroutine CNLivewoodTurnover (num_soilp, filter_soilp, & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) livewdcn => veg_vp%livewdcn , & ! Input: [real(r8) (:) ] live wood (phloem and ray parenchyma) C:N (gC/gN) deadwdcn => veg_vp%deadwdcn , & ! Input: [real(r8) (:) ] dead wood (xylem and heartwood) C:N (gC/gN) diff --git a/components/elm/src/biogeochem/CarbonStateUpdate1Mod.F90 b/components/elm/src/biogeochem/CarbonStateUpdate1Mod.F90 index 1eb7f64d025a..5133575abd30 100644 --- a/components/elm/src/biogeochem/CarbonStateUpdate1Mod.F90 +++ b/components/elm/src/biogeochem/CarbonStateUpdate1Mod.F90 @@ -200,7 +200,7 @@ subroutine CarbonStateUpdate1(bounds, & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Input: [integer (:) ] which pool is C taken from for a given decomposition step cascade_receiver_pool => decomp_cascade_con%cascade_receiver_pool , & ! Input: [integer (:) ] which pool is C added to for a given decomposition step harvdate => crop_vars%harvdate_patch & ! Input: [integer (:) ] harvest date diff --git a/components/elm/src/biogeochem/GrowthRespMod.F90 b/components/elm/src/biogeochem/GrowthRespMod.F90 index cb7ad522265f..cc6707d723b0 100644 --- a/components/elm/src/biogeochem/GrowthRespMod.F90 +++ b/components/elm/src/biogeochem/GrowthRespMod.F90 @@ -51,7 +51,7 @@ subroutine GrowthResp(num_soilp, filter_soilp) associate( & ivt => veg_pp%itype , & ! Input: [integer (:)] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:)] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:)] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) cpool_to_leafc => veg_cf%cpool_to_leafc , & ! Input: [real(r8) (:)] cpool_to_leafc_storage => veg_cf%cpool_to_leafc_storage , & ! Input: [real(r8) (:)] @@ -98,7 +98,7 @@ subroutine GrowthResp(num_soilp, filter_soilp) ) ! Loop through patches - ! start pft loop + ! start pft loop do fp = 1,num_soilp p = filter_soilp(fp) diff --git a/components/elm/src/biogeochem/MaintenanceRespMod.F90 b/components/elm/src/biogeochem/MaintenanceRespMod.F90 index 54005275b296..829e5c742a06 100644 --- a/components/elm/src/biogeochem/MaintenanceRespMod.F90 +++ b/components/elm/src/biogeochem/MaintenanceRespMod.F90 @@ -70,7 +70,7 @@ subroutine readMaintenanceRespParams ( ncid ) call ncd_io(varname=trim(tString),data=tempr, flag='read', ncid=ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(__FILE__, __LINE__)) br_mr_Inst = tempr - + end subroutine readMaintenanceRespParams !----------------------------------------------------------------------- @@ -107,7 +107,7 @@ subroutine MaintenanceResp(bounds, & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] patch vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) br_xr => veg_vp%br_xr , & ! Input: [real(r8) (:) ] base rate for excess respiration frac_veg_nosno => canopystate_vars%frac_veg_nosno_patch , & ! Input: [integer (:) ] fraction of vegetation not covered by snow (0 OR 1) [-] laisun => canopystate_vars%laisun_patch , & ! Input: [real(r8) (:) ] sunlit projected leaf area index @@ -160,7 +160,7 @@ subroutine MaintenanceResp(bounds, & ! calculate temperature corrections for each soil layer, for use in ! estimating fine root maintenance respiration with depth tcsoi(c,j) = Q10**((t_soisno(c,j)-SHR_CONST_TKFRZ - 20.0_r8)/10.0_r8) - + end do end do diff --git a/components/elm/src/biogeochem/NitrogenStateUpdate1Mod.F90 b/components/elm/src/biogeochem/NitrogenStateUpdate1Mod.F90 index c1b2bfc5d65a..f95b9a258cb1 100644 --- a/components/elm/src/biogeochem/NitrogenStateUpdate1Mod.F90 +++ b/components/elm/src/biogeochem/NitrogenStateUpdate1Mod.F90 @@ -130,7 +130,7 @@ subroutine NitrogenStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Input: [integer (:) ] which pool is C taken from for a given decomposition step cascade_receiver_pool => decomp_cascade_con%cascade_receiver_pool , & ! Input: [integer (:) ] which pool is C added to for a given decomposition step @@ -147,7 +147,7 @@ subroutine NitrogenStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp do j = 1, nlevdecomp do fc = 1,num_soilc c = filter_soilc(fc) - + ! N deposition and fixation (put all into NH4 pool) col_ns%smin_nh4_vr(c,j) = col_ns%smin_nh4_vr(c,j) + col_nf%ndep_to_sminn(c)*dt * ndep_prof(c,j) col_ns%smin_nh4_vr(c,j) = col_ns%smin_nh4_vr(c,j) + col_nf%nfix_to_sminn(c)*dt * nfixation_prof(c,j) @@ -219,7 +219,7 @@ subroutine NitrogenStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp end do end if end do - + do j = 1, nlevdecomp ! column loop do fc = 1,num_soilc @@ -254,8 +254,8 @@ subroutine NitrogenStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp end do ! end of column loop end do - - endif !end if is_active_betr_bgc + + endif !end if is_active_betr_bgc ! forest fertilization call get_curr_date(kyr, kmo, kda, mcsec) diff --git a/components/elm/src/biogeochem/PhenologyFluxLimitMod.F90 b/components/elm/src/biogeochem/PhenologyFluxLimitMod.F90 index 87b9637d5512..d80b844e6009 100644 --- a/components/elm/src/biogeochem/PhenologyFluxLimitMod.F90 +++ b/components/elm/src/biogeochem/PhenologyFluxLimitMod.F90 @@ -397,7 +397,7 @@ subroutine InitPhenoFluxLimiter() call spm_list_insert(spm_list, -1._r8, f_gresp_storage_to_xfer , s_gresp_storage,nelms) !turn the list into sparse matrix form call spm_list_to_mat(spm_list, spm_carbon_d, nelms, f_gresp_storage_to_xfer) - + !initialize stoichiometric relationship between carbon production flux and corresponding state varaibles call spm_list_init(spm_list, 1._r8, f_cpool_to_leafc , s_leafc, nelms) call spm_list_insert(spm_list, 1._r8, f_leafc_xfer_to_leafc , s_leafc, nelms) @@ -432,7 +432,7 @@ subroutine InitPhenoFluxLimiter() call spm_list_insert(spm_list, 1._r8, f_gresp_storage_to_xfer , s_gresp_xfer,nelms) call spm_list_insert(spm_list, 1._r8, f_cpool_to_gresp_storage , s_gresp_storage, nelms) - + !turn the list into sparse matrix form call spm_list_to_mat(spm_list, spm_carbon_p, nelms, f_gresp_storage_to_xfer) !initialize stoichiometry relationship between nutrient consumption and corresponding state variables @@ -474,7 +474,7 @@ subroutine InitPhenoFluxLimiter() call spm_list_insert(spm_list, -1._r8, f_grainn_to_food , s_grainn, nelms) call spm_list_insert(spm_list, -1._r8, f_grainn_xfer_to_grainn , s_grainn_xfer,nelms) call spm_list_insert(spm_list, -1._r8, f_retransn_to_npool , s_retransn, nelms) - !turn the list into sparse matrix form + !turn the list into sparse matrix form call spm_list_to_mat(spm_list, spm_nutrient_d, nelms, f_supplement_to_plantn) !initialize stoichiometry relationship between nutrient production and corresponding state variables call spm_list_init(spm_list, 1._r8, f_retransn_to_npool , s_npool, nelms) @@ -512,7 +512,7 @@ subroutine InitPhenoFluxLimiter() call spm_list_insert(spm_list, 1._r8, f_frootn_to_retransn , s_retransn, nelms) call spm_list_insert(spm_list, 1._r8, f_livestemn_to_retransn , s_retransn, nelms) call spm_list_insert(spm_list, 1._r8, f_livecrootn_to_retransn , s_retransn, nelms) - !turn the list into sparse matrix form + !turn the list into sparse matrix form call spm_list_to_mat(spm_list, spm_nutrient_p, nelms, f_supplement_to_plantn) end subroutine InitPhenoFluxLimiter !--------------------------------------------------------------------------- @@ -613,7 +613,7 @@ subroutine carbon_flux_limiter(bounds, num_soilc, filter_soilc,& real(r8) :: ar_p associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) harvdate => crop_vars%harvdate_patch & ! Input: [integer (:) ] harvest date ) ! set time steps @@ -847,7 +847,7 @@ subroutine nitrogen_flux_limiter(bounds, num_soilc, filter_soilc,& real(r8) :: dt associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) nf => veg_nf , & ns => veg_ns & ) @@ -1031,7 +1031,7 @@ subroutine phosphorus_flux_limiter(bounds, num_soilc, filter_soilc,& associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) pf => veg_pf , & ps => veg_ps & ) diff --git a/components/elm/src/biogeochem/PhenologyMod.F90 b/components/elm/src/biogeochem/PhenologyMod.F90 index 08e9b94da524..c46e251b60cd 100644 --- a/components/elm/src/biogeochem/PhenologyMod.F90 +++ b/components/elm/src/biogeochem/PhenologyMod.F90 @@ -163,7 +163,7 @@ subroutine readPhenolParams ( ncid ) allocate(PhenolParamsInst%lwtop ) ! ! read in parameters - ! + ! tString='crit_dayl' call ncd_io(varname=trim(tString),data=tempr, flag='read', ncid=ncid, readvar=readv) if ( .not. readv ) call endrun( msg=trim(errCode)//trim(tString)//errMsg(__FILE__, __LINE__)) @@ -175,7 +175,7 @@ subroutine readPhenolParams ( ncid ) else tString='crit_dayl_stress' call ncd_io(varname=trim(tString),data=tempr, flag='read', ncid=ncid, readvar=readv) - if ( .not. readv ) then + if ( .not. readv ) then PhenolParamsInst%crit_dayl_stress = secspqtrday else PhenolParamsInst%crit_dayl_stress = tempr @@ -238,7 +238,7 @@ subroutine readPhenolParams ( ncid ) tString='lwtop_ann' call ncd_io(varname=trim(tString),data=tempr, flag='read', ncid=ncid, readvar=readv) if ( .not. readv ) call endrun( msg=trim(errCode)//trim(tString)//errMsg(__FILE__, __LINE__)) - PhenolParamsInst%lwtop=tempr + PhenolParamsInst%lwtop=tempr !!!!========== Update to device ========= !!! !$acc update device(PhenolParamsInst%crit_dayl, & @@ -601,7 +601,7 @@ subroutine CNSeasonDecidPhenology (num_soilp, filter_soilp, cnstate_vars) prev_dayl => grc_pp%prev_dayl , & ! Input: [real(r8) (:) ] daylength from previous time step (s) season_decid => veg_vp%season_decid , & ! Input: [real(r8) (:) ] binary flag for seasonal-deciduous leaf habit (0 or 1) - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd) @@ -946,7 +946,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , & leaf_long => veg_vp%leaf_long , & ! Input: [real(r8) (:) ] leaf longevity (yrs) froot_long => veg_vp%froot_long , & ! Input: [real(r8) (:) ] fine root longevity (yrs) - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) stress_decid => veg_vp%stress_decid , & ! Input: [real(r8) (:) ] binary flag for stress-deciduous leaf habit (0 or 1) soilpsi => soilstate_vars%soilpsi_col , & ! Input: [real(r8) (:,:) ] soil water potential in each soil layer (MPa) @@ -1445,7 +1445,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp, & froot_long => veg_vp%froot_long , & ! Input: [real(r8) (:) ] fine root longevity (yrs) leafcn => veg_vp%leafcn , & ! Input: [real(r8) (:) ] leaf C:N (gC/gN) - manunitro => veg_vp%manunitro , & ! Input: max manure to apply (kgN/m2) + manunitro => veg_vp%manunitro , & ! Input: max manure to apply (kgN/m2) t_ref2m_min => veg_es%t_ref2m_min , & ! Input: [real(r8) (:) ] daily minimum of average 2 m height surface air temperature (K) t10 => veg_es%t_a10 , & ! Input: [real(r8) (:) ] 10-day running mean of the 2 m temperature (K) a5tmin => veg_es%t_a5min , & ! Input: [real(r8) (:) ] 5-day running mean of min 2-m temperature @@ -1988,7 +1988,7 @@ subroutine PerennialCropPhenology(num_ppercropp, filter_ppercropp, & froot_long => veg_vp%froot_long , & ! Input: [real(r8) (:) ] fine root longevity (yrs) leafcn => veg_vp%leafcn , & ! Input: [real(r8) (:) ] leaf C:N (gC/gN) leafcp => veg_vp%leafcp , & ! Input: [real(r8) (:) ] leaf C:P (gC/gP) - manunitro => veg_vp%manunitro , & ! Input: max manure to apply (kgN/m2) + manunitro => veg_vp%manunitro , & ! Input: max manure to apply (kgN/m2) t10 => veg_es%t_a10 , & ! Input: [real(r8) (:) ] 10-day running mean of the 2 m temperature (K) a10tmin => veg_es%t_a10min , & ! Input: [real(r8) (:) ] 10-day running mean of min 2-m temperature fertnitro => crop_vars%fertnitro_patch , & ! Input: [real(r8) (:) ] max fertilizer to be applied in total (kgN/m2) @@ -2016,8 +2016,8 @@ subroutine PerennialCropPhenology(num_ppercropp, filter_ppercropp, & crop_seedc_to_leaf => veg_cf%crop_seedc_to_leaf , & ! Output: [real(r8) (:) ] (gC/m2/s) seed source to PFT-level - synthfert => veg_nf%synthfert , & ! Output: [real(r8) (:) ] (gN/m2/s) fertilizer applied each timestep - manure => veg_nf%manure , & ! Output: [real(r8) (:) ] (gN/m2/s) manure applied each timestep + synthfert => veg_nf%synthfert , & ! Output: [real(r8) (:) ] (gN/m2/s) fertilizer applied each timestep + manure => veg_nf%manure , & ! Output: [real(r8) (:) ] (gN/m2/s) manure applied each timestep fert_p => veg_pf%fert_p , & ! Output: [real(r8) (:) ] (gP/m2/s) phosphorus fertilizer applied each timestep fert_counter => veg_nf%fert_counter , & ! Output: [real(r8) (:) ] >0 fertilize; <=0 not (seconds) @@ -2550,22 +2550,22 @@ subroutine CropPlantDate (num_soilp, filter_soilp, num_pcropp, filter_pcropp, & xt(p,kmo) = xt(p,kmo) + t_ref2m(p) * fracday/ndaypm(kmo) ! monthly average temperature xp(p,kmo) = xp(p,kmo) + (forc_rain(t)+forc_snow(t))*dt ! monthly average precipitation - ! calculate the potential evapotranspiration + ! calculate the potential evapotranspiration netrad = fsa(p) + eflx_lwrad_net(p) ! moved this here because it is calculated too late call calculate_eto(t_ref2m(p), netrad, eflx_soil_grnd(p), forc_pbot(t), forc_rh(t), forc_wind(t), dt, ETout) ! monthly ETo ETo(p,kmo) = ETo(p,kmo) + ETout - + ! calculate the P:PET for each month - if ( abs(ETo(p,kmo)) > 0._r8) then + if ( abs(ETo(p,kmo)) > 0._r8) then p2ETo(p,kmo) = xp(p,kmo)/ETo(p,kmo) else ! P:PET is undefined. ! Setting to a fill value ( 'spval' ) would - ! require nested if statements due to + ! require nested if statements due to ! the weighting of previous years (i.e., p2ETo and prev_p2ETo_bar ) ! So, set to zero for simplicity. p2ETo(p,kmo) = 0._r8 - end if + end if if (nyrs_crop_active(p) == 0) then ! for the first year, use last years values prev_xt_bar(p,kmo) = xt(p,kmo) @@ -2617,7 +2617,7 @@ subroutine CNOnsetGrowth (num_soilp, filter_soilp, & associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) onset_flag => cnstate_vars%onset_flag_patch , & ! Input: [real(r8) (:) ] onset flag onset_counter => cnstate_vars%onset_counter_patch , & ! Input: [real(r8) (:) ] onset days counter @@ -3277,7 +3277,7 @@ subroutine CNLivewoodTurnover (num_soilp, filter_soilp) associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) livewdcn => veg_vp%livewdcn , & ! Input: [real(r8) (:) ] live wood (phloem and ray parenchyma) C:N (gC/gN) deadwdcn => veg_vp%deadwdcn , & ! Input: [real(r8) (:) ] dead wood (xylem and heartwood) C:N (gC/gN) diff --git a/components/elm/src/biogeochem/PhosphorusStateUpdate1Mod.F90 b/components/elm/src/biogeochem/PhosphorusStateUpdate1Mod.F90 index 6e1111169221..9297ab41e1bf 100644 --- a/components/elm/src/biogeochem/PhosphorusStateUpdate1Mod.F90 +++ b/components/elm/src/biogeochem/PhosphorusStateUpdate1Mod.F90 @@ -75,7 +75,7 @@ subroutine PhosphorusStateUpdateDynPatch(bounds, num_soilc_with_inactive,& c = filter_soilc_with_inactive(fc) col_ps%prod10p(c) = col_ps%prod10p(c) + col_pf%dwt_prod10p_gain(c)*dt col_ps%prod100p(c) = col_ps%prod100p(c) + col_pf%dwt_prod100p_gain(c)*dt - col_ps%prod1p(c) = col_ps%prod1p(c) + col_pf%dwt_crop_productp_gain(c)*dt + col_ps%prod1p(c) = col_ps%prod1p(c) + col_pf%dwt_crop_productp_gain(c)*dt do j = 1,nlevdecomp @@ -125,7 +125,7 @@ subroutine PhosphorusStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soi associate( & ivt => veg_pp%itype , & ! Input: [integer (:) ] pft vegetation type - woody => veg_vp%woody , & ! Input: [real(r8) (:) ] binary flag for woody lifeform (1=woody, 0=not woody) + woody => veg_vp%woody , & ! Input: [real(r8) (:) ] woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Input: [integer (:) ] which pool is C taken from for a given decomposition step cascade_receiver_pool => decomp_cascade_con%cascade_receiver_pool , & ! Input: [integer (:) ] which pool is C added to for a given decomposition step @@ -147,18 +147,18 @@ subroutine PhosphorusStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soi do j = 1, nlevdecomp do fc = 1,num_soilc c = filter_soilc(fc) - + ! plant to litter fluxes ! phenology and dynamic landcover fluxes col_pf%decomp_ppools_sourcesink(c,j,i_met_lit) = & col_pf%phenology_p_to_litr_met_p(c,j) * dt - + col_pf%decomp_ppools_sourcesink(c,j,i_cel_lit) = & col_pf%phenology_p_to_litr_cel_p(c,j) * dt - + col_pf%decomp_ppools_sourcesink(c,j,i_lig_lit) = & col_pf%phenology_p_to_litr_lig_p(c,j) * dt - + end do end do end if diff --git a/components/elm/src/biogeochem/VegStructUpdateMod.F90 b/components/elm/src/biogeochem/VegStructUpdateMod.F90 index 5943da355a04..74e135f6f0ea 100644 --- a/components/elm/src/biogeochem/VegStructUpdateMod.F90 +++ b/components/elm/src/biogeochem/VegStructUpdateMod.F90 @@ -40,6 +40,7 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & use pftvarcon , only : noveg, woody, iscft, crop use pftvarcon , only : ncorn, ncornirrig, ztopmx, laimx use pftvarcon , only : nmiscanthus, nmiscanthusirrig, nswitchgrass, nswitchgrassirrig + use pftvarcon , only : bendresist, vegshape, stocking, taper use elm_time_manager , only : get_rad_step_size use elm_varctl , only : spinup_state, spinup_mortality_factor ! @@ -58,8 +59,6 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & ! !LOCAL VARIABLES: integer :: p,c,g ! indices integer :: fp ! lake filter indices - real(r8) :: taper ! ratio of height:radius_breast_height (tree allometry) - real(r8) :: stocking ! #stems / ha (stocking density) real(r8) :: ol ! thickness of canopy layer covered by snow (m) real(r8) :: fb ! fraction of canopy layer covered by snow real(r8) :: tlai_old ! for use in Zeng tsai formula @@ -88,6 +87,9 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & z0mr => veg_vp%z0mr , & ! Input: [real(r8) (:) ] ratio of momentum roughness length to canopy top height (-) displar => veg_vp%displar , & ! Input: [real(r8) (:) ] ratio of displacement height to canopy top height (-) dwood => veg_vp%dwood , & ! Input: [real(r8) (:) ] density of wood (gC/m^3) + bendresist => veg_vp%bendresist , & ! Input: [real(r8) (:) ] resistance to bending under snow loading Sturm et al. 2005 (-) [0,1] + vegshape => veg_vp%vegshape , & ! Input: [real(r8) (:) ] vegetation shape for calculating snow buried fraction only (-) (Liston and Heimstra, 2011) + stocking => veg_vp%stocking , & ! Input: [real(r8) (:) ] Stocking density [stems / hectare] snow_depth => col_ws%snow_depth , & ! Input: [real(r8) (:) ] snow height (m) @@ -113,13 +115,6 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & dt = real( get_rad_step_size(), r8 ) - ! constant allometric parameters - taper = 200._r8 - stocking = 1000._r8 - - ! convert from stems/ha -> stems/m^2 - stocking = stocking / 10000._r8 - ! patch loop do fp = 1,num_soilp p = filter_soilp(fp) @@ -160,24 +155,16 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & if (woody(ivt(p)) >= 1.0_r8) then - ! trees and shrubs - - ! if shrubs have a squat taper - if (woody(ivt(p)) == 2.0_r8) then - taper = 10._r8 - ! otherwise have a tall taper - else - taper = 200._r8 - end if - - ! trees and shrubs for now have a very simple allometry, with hard-wired - ! stem taper (height:radius) and hard-wired stocking density (#individuals/area) + ! trees and shrubs currently have a very simple allometry, based on + ! stem taper (height:radius) and stocking density (#individuals/area) + ! taper and stocking density can be set as input variables now to + ! change from default values set in pftvarcon.F90 if (spinup_state >= 1) then - htop(p) = ((3._r8 * deadstemc(p) * spinup_mortality_factor * taper * taper)/ & - (SHR_CONST_PI * stocking * dwood(ivt(p))))**(1._r8/3._r8) + htop(p) = ((3._r8 * deadstemc(p) * spinup_mortality_factor * taper(ivt(p)) * taper(ivt(p)))/ & + (SHR_CONST_PI * stocking(ivt(p)) * dwood(ivt(p))))**(1._r8/3._r8) else - htop(p) = ((3._r8 * deadstemc(p) * taper * taper)/ & - (SHR_CONST_PI * stocking * dwood(ivt(p))))**(1._r8/3._r8) + htop(p) = ((3._r8 * deadstemc(p) * taper(ivt(p)) * taper(ivt(p)))/ & + (SHR_CONST_PI * stocking(ivt(p)) * dwood(ivt(p))))**(1._r8/3._r8) end if ! Peter Thornton, 5/3/2004 @@ -246,10 +233,12 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, & ! adjust lai and sai for burying by snow. ! snow burial fraction for short vegetation (e.g. grasses) as in - ! Wang and Zeng, 2007. - if (woody(ivt(p)) >= 1.0_r8 ) then + ! Wang and Zeng et al 2007. + ! Taller vegetation (trees and shrubs) have been updated to use formulation similar to + ! Sturm et al. 2005; Liston and Hiemstra, 2011; and Belke-Brea et al. 2020 + if ( woody(ivt(p)) >= 1.0_r8 ) then ol = min( max(snow_depth(c)-hbot(p), 0._r8), htop(p)-hbot(p)) - fb = 1._r8 - ol / max(1.e-06_r8, htop(p)-hbot(p)) + fb = 1._r8 - (ol / max(1.e-06_r8, bendresist(ivt(p)) * (htop(p)-hbot(p)))) ** vegshape(ivt(p)) else fb = 1._r8 - max(min(snow_depth(c),0.2_r8),0._r8)/0.2_r8 ! 0.2m is assumed !depth of snow required for complete burial of grasses diff --git a/components/elm/src/cpl/lnd_downscale_atm_forcing.F90 b/components/elm/src/cpl/lnd_downscale_atm_forcing.F90 index 99b03a2e25ae..3b69158f49eb 100644 --- a/components/elm/src/cpl/lnd_downscale_atm_forcing.F90 +++ b/components/elm/src/cpl/lnd_downscale_atm_forcing.F90 @@ -35,10 +35,13 @@ module lnd_downscale_atm_forcing ! ! !PUBLIC MEMBER FUNCTIONS: public :: downscale_atm_forcing_to_topounit ! Calls downscaling subroutines of forcing fields from gridcell to topounit + public :: downscale_atm_forcing_to_topounit_cpl_bypass ! Calls downscaling subroutines of forcing fields from gridcell to topounit, for use within the CPL_BYPASS code ! ! !PRIVATE MEMBER FUNCTIONS: private :: downscale_atm_state_to_topounit ! Downscale atmosperic state fields from gridcell to topounit + private :: downscale_atm_state_to_topounit_cpl_bypass ! Downscale atmosperic state fields from gridcell to topounit, for cpl_bypass code private :: downscale_longwave_to_topounit ! Downscale longwave radiation field from gridcell to topounit + private :: downscale_longwave_to_topounit_cpl_bypass ! Downscale longwave radiation field from gridcell to topounit, for cpl_bypass code private :: downscale_precip_to_topounit_FNM ! Downscale precipitation field from gridcell to topounit using Froude number method (FNM) private :: downscale_precip_to_topounit_ERMM ! Downscale precipitation field from gridcell to topounit using elevation ration with maximum elevation method (ERMM) private :: build_normalization ! Compute normalization factors so that downscaled fields are conservative @@ -787,4 +790,526 @@ subroutine build_normalization(orig_field, sum_field, sum_wts, norms) end subroutine build_normalization + subroutine downscale_atm_forcing_to_topounit_cpl_bypass(g, atm2lnd_vars, lnd2atm_vars) + ! + ! !DESCRIPTION: + ! Downscale fields from gridcell to topounit + ! + ! Downscaling is done over topounits if the number of topounits > 1. + ! + ! !USES: + use elm_time_manager, only : get_nstep + use elm_varcon , only : rair, cpair, grav, lapse_glcmec + use elm_varcon , only : glcmec_rain_snow_threshold, o2_molar_const + use shr_const_mod , only : SHR_CONST_TKFRZ + use landunit_varcon , only : istice_mec + use elm_varctl , only : glcmec_downscale_rain_snow_convert + use domainMod , only : ldomain + use QsatMod , only : Qsat + use FrictionVelocityMod, only: atm_gustiness + ! + ! !ARGUMENTS: + integer , intent(in) :: g + type(atm2lnd_type) , intent(in) :: atm2lnd_vars + type(lnd2atm_type) , intent(in) :: lnd2atm_vars + + ! + ! !LOCAL VARIABLES: + integer :: t, l, c, fc, t2 ! indices + integer :: clo, cc + integer :: numt_pg ! Number of topounits per grid + + ! temporaries for topo downscaling + real(r8) :: rain_g, snow_g + real(r8) :: mxElv ! Maximum elevation value per grid + real(r8) :: uovern_t ! Froude Number + real(r8) :: grdElv ! Grid elevation + real(r8) :: topoElv ! Topounit elevation + real(r8) :: max_tpuElv ! Maximum topounit elevation for calculating elevation range + real(r8) :: min_tpuElv ! Minimum topounit elevation for calculating elevation range + real(r8) :: crnt_temp_t ! Current downscaled topounit temperature + real(r8) :: temp_r ! Temporary topounit rainfall + real(r8) :: temp_s ! Temporary topounit snowfall + real(r8) :: t_th ! Temperature threshold for snowfall + real(r8) :: Ta_th1 ! Temperature at 0.5 Celsius in K 273.65 + real(r8) :: Ta_th2 ! Temperature at 2.0 Celsius in K 275.15 + real(r8) :: Ta_th3 ! Temperature at 2.5 Celsius in K 275.65 + real(r8) :: tmp_Snow_frc ! Current snow fraction + + real(r8) :: e + real(r8) :: qvsat + + real(r8) :: sum_qbot_g ! weighted sum of column-level lwrad + real(r8) :: sum_wtsq_g ! sum of weights that contribute to sum_lwrad_g + real(r8) :: qbot_norm_g ! normalization factors + real(r8) :: sum_lwrad_g ! weighted sum of column-level lwrad + real(r8) :: sum_wtslw_g ! sum of weights that contribute to sum_lwrad_g + real(r8) :: lwrad_norm_g ! normalization factors + real(r8) :: esatw ! saturation vapor pressure over water (Pa) + real(r8) :: esati ! saturation vapor pressure over ice (Pa) + real(r8) :: a0,a1,a2,a3,a4,a5,a6 ! coefficients for esat over water + real(r8) :: b0,b1,b2,b3,b4,b5,b6 ! coefficients for esat over ice + real(r8) :: tdc, temp ! Kelvins to Celcius function and its input + real(r8) :: vp ! water vapor pressure (Pa) + + real(r8), allocatable :: deltaRain(:) ! Deviation of subgrid rain from grid rain + real(r8), allocatable :: deltaSnow(:) ! Deviation of subgrid snow from grid snow + real(r8) :: deltaR ! Temporary deltaRain + real(r8) :: deltaS ! Temporary deltaSnow + real(r8) :: sum_of_hrise ! Sum of height rise of air parcel of all subgrids of a grid + real(r8) :: hrise ! Temporary height rise + real(r8) :: elvrnge ! Elevation range between lowest and highest tpu + real(r8) :: ave_elv + integer :: elv_flag ! Elevation flag to trac grids with +ve grid elevation and -ve tpu elevation + + integer :: uaflag = 0 + integer :: precip_dwn = 0 ! Used to turn on/off the downscaling of precipitation 0 = on; 1 = off + integer :: other_forcing_dwn = 0 ! Used to turn on/off the downscaling of other forcing 0 = on; 1 = off + + character(len=*), parameter :: subname = 'downscale_atm_forcing_to_topounit' + !---------------------------------------------------------------------------------------- + + ! Constants to compute vapor pressure + parameter (a0=6.107799961_r8 , a1=4.436518521e-01_r8, & + a2=1.428945805e-02_r8, a3=2.650648471e-04_r8, & + a4=3.031240396e-06_r8, a5=2.034080948e-08_r8, & + a6=6.136820929e-11_r8) + + parameter (b0=6.109177956_r8 , b1=5.034698970e-01_r8, & + b2=1.886013408e-02_r8, b3=4.176223716e-04_r8, & + b4=5.824720280e-06_r8, b5=4.838803174e-08_r8, & + b6=1.838826904e-10_r8) + + ! + ! function declarations + ! + tdc(temp) = min( 50._r8, max(-50._r8,(temp-SHR_CONST_TKFRZ)) ) ! Taken from lnd_import_export.F90 + esatw(temp) = 100._r8*(a0+temp*(a1+temp*(a2+temp*(a3+temp*(a4+temp*(a5+temp*a6)))))) ! Taken from lnd_import_export.F90 + esati(temp) = 100._r8*(b0+temp*(b1+temp*(b2+temp*(b3+temp*(b4+temp*(b5+temp*b6)))))) ! Taken from lnd_import_export.F90 + !----------------------------------------------------------------------- + ! Get required inputs + numt_pg = grc_pp%ntopounits(g) ! Number of topounits per grid + grdElv = grc_pp%elevation(g) ! Grid level sfc elevation + mxElv = grc_pp%MaxElevation(g) ! Maximum src elevation per grid obtained from the highest elevation topounit + uovern_t = atm2lnd_vars%forc_uovern(g) ! Froude Number + + snow_g = atm2lnd_vars%forc_snow_not_downscaled_grc(g) + rain_g = atm2lnd_vars%forc_rain_not_downscaled_grc(g) + + sum_qbot_g = 0._r8 + sum_wtsq_g = 0._r8 + sum_lwrad_g = 0._r8 + sum_wtslw_g = 0._r8 + + sum_of_hrise = 0._r8 + t_th = 273.15_r8 ! Freezing temperature in K + Ta_th1 = 273.65_r8 ! Lowest threshold for snow calculation + Ta_th2 = 275.15_r8 ! Middle threshold for rain/snow partitioning + Ta_th3 = 275.65_r8 ! Highest threshold for rain/snow partitioning + tmp_Snow_frc = 0._r8 ! Snow fraction + elv_flag = 0 + elvrnge = 0._r8 + ave_elv = 0._r8 + max_tpuElv = 0._r8 + min_tpuElv = 0._r8 + if (numt_pg > 1) then !downscaling is done only if a grid has more than 1 topounits + if (precip_downscaling_method == 'FNM') then + allocate(deltaRain(numt_pg)) + deltaRain(:) = 0._r8 + allocate(deltaSnow(numt_pg)) + deltaSnow(:) = 0._r8 + hrise = 0._r8 + end if + + ! calculate elevation range and track grids with +ve elevation but have -ve tpu elevation + max_tpuElv = -100000._r8 + min_tpuElv = 100000._r8 + do t = grc_pp%topi(g), grc_pp%topf(g) ! Check occurence of grid elevation is +ve while tpu elevation is -ve + topoElv = top_pp%elevation(t) + if (topoElv > max_tpuElv) then + max_tpuElv = topoElv + end if + if (topoElv < min_tpuElv) then + min_tpuElv = topoElv + end if + end do + elvrnge = max_tpuElv - min_tpuElv + + do t = grc_pp%topi(g), grc_pp%topf(g) + t2 = t - grc_pp%topi(g) + 1 + topoElv = top_pp%elevation(t) ! Topounit sfc elevation + + ! Downscale precipitation + if (mxElv == 0._r8 .or. precip_dwn == 1) then ! avoid dividing by 0 + top_af%rain(t) = rain_g + top_af%snow(t) = snow_g + else + if (precip_downscaling_method == 'FNM') then + call downscale_precip_to_topounit_FNM(mxElv,uovern_t,grdElv,topoElv,rain_g,snow_g,deltaR,deltaS,hrise) !Use FNM method + deltaRain(t2) = deltaR + deltaSnow(t2) = deltaS + sum_of_hrise = sum_of_hrise + hrise + else + call downscale_precip_to_topounit_ERMM(t,mxElv,grdElv,topoElv,rain_g, snow_g,elv_flag,elvrnge) ! Use ERMM method + end if + end if + + ! Downscale other fluxes + if (other_forcing_dwn == 1) then ! flag to turn on or off downscaling of other forcing + top_af%rain(t) = rain_g + top_af%snow(t) = snow_g + top_af%lwrad(t) = atm2lnd_vars%forc_lwrad_not_downscaled_grc(g) + ! Update top_as + top_as%tbot(t) = atm2lnd_vars%forc_t_not_downscaled_grc(g) ! forc_txy Atm state K + top_as%thbot(t) = atm2lnd_vars%forc_th_not_downscaled_grc(g) ! forc_thxy Atm state K + top_as%pbot(t) = atm2lnd_vars%forc_pbot_not_downscaled_grc(g) ! ptcmxy Atm state Pa + top_as%qbot(t) = atm2lnd_vars%forc_q_not_downscaled_grc(g) ! forc_qxy Atm state kg/kg + top_as%ubot(t) = atm2lnd_vars%forc_u_grc(g) ! forc_uxy Atm state m/s + top_as%vbot(t) = atm2lnd_vars%forc_v_grc(g) ! forc_vxy Atm state m/s + top_as%zbot(t) = atm2lnd_vars%forc_hgt_grc(g) ! zgcmxy Atm state m + + ! assign the state forcing fields derived from other inputs + ! Horizontal windspeed (m/s) + top_as%windbot(t) = sqrt(top_as%ubot(t)**2 + top_as%vbot(t)**2) + ! Relative humidity (percent) + if (top_as%tbot(t) > SHR_CONST_TKFRZ) then + e = esatw(tdc(top_as%tbot(t))) + else + e = esati(tdc(top_as%tbot(t))) + end if + qvsat = 0.622_r8*e / (top_as%pbot(t) - 0.378_r8*e) + top_as%rhbot(t) = 100.0_r8*(top_as%qbot(t) / qvsat) + ! partial pressure of oxygen (Pa) + top_as%po2bot(t) = o2_molar_const * top_as%pbot(t) + ! air density (kg/m**3) - uses a temporary calculation + ! of water vapor pressure (Pa) + vp = top_as%qbot(t) * top_as%pbot(t) / (0.622_r8 + 0.378_r8 * top_as%qbot(t)) + top_as%rhobot(t) = (top_as%pbot(t) - 0.378_r8 * vp) / (rair * top_as%tbot(t)) + + top_af%solad(t,2) = atm2lnd_vars%forc_solad_grc(g,2) + top_af%solad(t,1) = atm2lnd_vars%forc_solad_grc(g,1) + top_af%solai(t,2) = atm2lnd_vars%forc_solai_grc(g,2) + top_af%solai(t,1) = atm2lnd_vars%forc_solai_grc(g,1) + ! derived flux forcings + top_af%solar(t) = top_af%solad(t,2) + top_af%solad(t,1) + & + top_af%solai(t,2) + top_af%solai(t,1) + else + call downscale_atm_state_to_topounit_cpl_bypass(g, t, atm2lnd_vars) + call downscale_longwave_to_topounit_cpl_bypass(g, t, atm2lnd_vars) + top_as%ubot(t) = atm2lnd_vars%forc_u_grc(g) ! forc_uxy Atm state m/s + top_as%vbot(t) = atm2lnd_vars%forc_v_grc(g) ! forc_vxy Atm state m/s + top_as%zbot(t) = atm2lnd_vars%forc_hgt_grc(g) ! zgcmxy Atm state m + + sum_qbot_g = sum_qbot_g + top_pp%wtgcell(t)*top_as%qbot(t) + sum_wtsq_g = sum_wtsq_g + top_pp%wtgcell(t) + + ! assign the state forcing fields derived from other inputs + ! Horizontal windspeed (m/s) + top_as%windbot(t) = sqrt(top_as%ubot(t)**2 + top_as%vbot(t)**2) + ! partial pressure of oxygen (Pa) + top_as%po2bot(t) = o2_molar_const * top_as%pbot(t) + ! air density (kg/m**3) - uses a temporary calculation + ! of water vapor pressure (Pa) + vp = top_as%qbot(t) * top_as%pbot(t) / (0.622_r8 + 0.378_r8 * top_as%qbot(t)) + top_as%rhobot(t) = (top_as%pbot(t) - 0.378_r8 * vp) / (rair * top_as%tbot(t)) + top_af%solad(t,2) = atm2lnd_vars%forc_solad_grc(g,2) + top_af%solad(t,1) = atm2lnd_vars%forc_solad_grc(g,1) + top_af%solai(t,2) = atm2lnd_vars%forc_solai_grc(g,2) + top_af%solai(t,1) = atm2lnd_vars%forc_solai_grc(g,1) + ! derived flux forcings + top_af%solar(t) = top_af%solad(t,2) + top_af%solad(t,1) + & + top_af%solai(t,2) + top_af%solai(t,1) + + ! Keep track of the gridcell-level weighted sum for later normalization. + ! + ! This gridcell-level weighted sum just includes points for which we do the + ! downscaling (e.g., glc_mec points). Thus the contributing weights + ! generally do not add to 1. So to do the normalization properly, we also + ! need to keep track of the weights that have contributed to this sum. + sum_lwrad_g = sum_lwrad_g + top_pp%wtgcell(t)*top_af%lwrad(t) + sum_wtslw_g = sum_wtslw_g + top_pp%wtgcell(t) + end if + end do + if (precip_downscaling_method == 'FNM') then + do t = grc_pp%topi(g), grc_pp%topf(g) + t2 = t - grc_pp%topi(g) + 1 + if (mxElv == 0.) then ! avoid dividing by 0 + top_af%rain(t) = rain_g + top_af%snow(t) = snow_g + else + top_af%rain(t) = rain_g + (deltaRain(t2) - (rain_g/mxElv)*(sum_of_hrise/numt_pg)) + top_af%snow(t) = snow_g + (deltaSnow(t2) - (snow_g/mxElv)*(sum_of_hrise/numt_pg)) + end if + + end do + deallocate(deltaRain) + deallocate(deltaSnow) + end if + + ! Precipitation partitioning using simple method following Jordan (1991) + do t = grc_pp%topi(g), grc_pp%topf(g) + crnt_temp_t = top_as%tbot(t) + if (crnt_temp_t > Ta_th3) then ! No snow or all rain + temp_r = top_af%rain(t) + top_af%snow(t) + temp_s = 0._r8 + else if (crnt_temp_t >= Ta_th2 .and. crnt_temp_t <= Ta_th3) then + temp_s = (top_af%snow(t) + top_af%rain(t))*0.6_r8 ! 0.6 fraction of the total precip is snow + temp_r = (top_af%snow(t) + top_af%rain(t)) - temp_s + else if (crnt_temp_t < Ta_th2 .and. crnt_temp_t > Ta_th1) then + tmp_Snow_frc = (crnt_temp_t-Ta_th2)*(0.4_r8/(Ta_th1-Ta_th2))+0.6_r8 ! Snow fraction value 1-0.6 = 0.4 snowfrc is 1 at t = Ta_th1 + temp_s = (top_af%snow(t) + top_af%rain(t))*tmp_Snow_frc ! 0.6 is snow fraction at t = Ta_th2 + temp_r = (top_af%snow(t) + top_af%rain(t)) - temp_s + else ! crnt_temp_t <= Ta_th1 ==> all snow + temp_s = top_af%snow(t) + top_af%rain(t) + temp_r = 0._r8 + + end if + + top_af%rain(t) = temp_r + top_af%snow(t) = temp_s + + end do + + else !grid has a single topounit + ! update top_af using grid level values + t = grc_pp%topi(g) + top_af%rain(t) = rain_g + top_af%snow(t) = snow_g + top_af%lwrad(t) = atm2lnd_vars%forc_lwrad_not_downscaled_grc(g) + + ! Update top_as + top_as%tbot(t) = atm2lnd_vars%forc_t_not_downscaled_grc(g) ! forc_txy Atm state K + top_as%thbot(t) = atm2lnd_vars%forc_th_not_downscaled_grc(g) ! forc_thxy Atm state K + top_as%pbot(t) = atm2lnd_vars%forc_pbot_not_downscaled_grc(g) ! ptcmxy Atm state Pa + top_as%qbot(t) = atm2lnd_vars%forc_q_not_downscaled_grc(g) ! forc_qxy Atm state kg/kg + top_as%ubot(t) = atm2lnd_vars%forc_u_grc(g) ! forc_uxy Atm state m/s + top_as%vbot(t) = atm2lnd_vars%forc_v_grc(g) ! forc_vxy Atm state m/s + top_as%zbot(t) = atm2lnd_vars%forc_hgt_grc(g) ! zgcmxy Atm state m + + ! assign the state forcing fields derived from other inputs + ! Horizontal windspeed (m/s) + top_as%windbot(t) = sqrt(top_as%ubot(t)**2 + top_as%vbot(t)**2) + if (atm_gustiness) then + top_as%windbot(t) = sqrt(top_as%windbot(t)**2 + top_as%ugust(t)**2) + end if + ! Relative humidity (percent) + if (top_as%tbot(t) > SHR_CONST_TKFRZ) then + e = esatw(tdc(top_as%tbot(t))) + else + e = esati(tdc(top_as%tbot(t))) + end if + qvsat = 0.622_r8*e / (top_as%pbot(t) - 0.378_r8*e) + top_as%rhbot(t) = 100.0_r8*(top_as%qbot(t) / qvsat) + ! partial pressure of oxygen (Pa) + top_as%po2bot(t) = o2_molar_const * top_as%pbot(t) + ! air density (kg/m**3) - uses a temporary calculation + ! of water vapor pressure (Pa) + vp = top_as%qbot(t) * top_as%pbot(t) / (0.622_r8 + 0.378_r8 * top_as%qbot(t)) + top_as%rhobot(t) = (top_as%pbot(t) - 0.378_r8 * vp) / (rair * top_as%tbot(t)) + + top_af%solad(t,2) = atm2lnd_vars%forc_solad_grc(g,2) + top_af%solad(t,1) = atm2lnd_vars%forc_solad_grc(g,1) + top_af%solai(t,2) = atm2lnd_vars%forc_solai_grc(g,2) + top_af%solai(t,1) = atm2lnd_vars%forc_solai_grc(g,1) + ! derived flux forcings + top_af%solar(t) = top_af%solad(t,2) + top_af%solad(t,1) + & + top_af%solai(t,2) + top_af%solai(t,1) + + end if + + if (numt_pg > 1) then + + ! Normalize forc_qbot to conserve energy + + call build_normalization(orig_field=atm2lnd_vars%forc_q_not_downscaled_grc(g), & + sum_field=sum_qbot_g, sum_wts=sum_wtsq_g, norms=qbot_norm_g) + + do t = grc_pp%topi(g), grc_pp%topf(g) + top_as%qbot(t) = top_as%qbot(t) * qbot_norm_g + + ! Relative humidity (percent) + if (top_as%tbot(t) > SHR_CONST_TKFRZ) then + e = esatw(tdc(top_as%tbot(t))) + else + e = esati(tdc(top_as%tbot(t))) + end if + qvsat = 0.622_r8*e / (top_as%pbot(t) - 0.378_r8*e) + top_as%rhbot(t) = 100.0_r8*(top_as%qbot(t) / qvsat) + ! partial pressure of oxygen (Pa) + top_as%po2bot(t) = o2_molar_const * top_as%pbot(t) + ! air density (kg/m**3) - uses a temporary calculation + ! of water vapor pressure (Pa) + vp = top_as%qbot(t) * top_as%pbot(t) / (0.622_r8 + 0.378_r8 * top_as%qbot(t)) + top_as%rhobot(t) = (top_as%pbot(t) - 0.378_r8 * vp) / (rair * top_as%tbot(t)) + + end do + + + ! Normalize forc_lwrad_c(c) to conserve energy + + call build_normalization(orig_field=atm2lnd_vars%forc_lwrad_not_downscaled_grc(g), & + sum_field=sum_lwrad_g, sum_wts=sum_wtslw_g, norms=lwrad_norm_g) + + do t = grc_pp%topi(g), grc_pp%topf(g) + top_af%lwrad(t) = top_af%lwrad(t) * lwrad_norm_g + end do + + end if + + end subroutine downscale_atm_forcing_to_topounit_cpl_bypass + + !----------------------------------------------------------------------- + ! Downscale other atmospheric state variables + !----------------------------------------------------------------------- + subroutine downscale_atm_state_to_topounit_cpl_bypass(g, t, atm2lnd_vars) + ! + ! !DESCRIPTION: + ! Downscale atmospheric forcing fields from gridcell to topounit + ! + ! Downscaling is done over topounits. + ! + ! !USES: + use elm_time_manager, only : get_nstep + use elm_varcon , only : rair, cpair, grav, lapse_glcmec + use elm_varcon , only : glcmec_rain_snow_threshold + use landunit_varcon , only : istice_mec + use elm_varctl , only : glcmec_downscale_rain_snow_convert + use domainMod , only : ldomain + use QsatMod , only : Qsat + ! + ! !ARGUMENTS: + integer , intent(in) :: g + integer , intent(in) :: t + type(atm2lnd_type) , intent(in) :: atm2lnd_vars + ! + ! !LOCAL VARIABLES: + integer :: l, c, fc ! indices + integer :: clo, cc + integer :: nstep + + ! temporaries for topo downscaling + real(r8) :: hsurf_g,hsurf_t,Hbot + real(r8) :: zbot_g, tbot_g, pbot_g, thbot_g, qbot_g, qs_g, es_g + real(r8) :: zbot_t, tbot_t, pbot_t, thbot_t, qbot_t, qs_t, es_t + real(r8) :: egcm_t, rhos_t + real(r8) :: dum1, dum2 + + character(len=*), parameter :: subname = 'downscale_atm_state_to_topounit' + !----------------------------------------------------------------------- + + nstep = get_nstep() + + ! Downscale forc_t, forc_th, forc_q, forc_pbot, and forc_rho to columns. + ! For glacier_mec columns the downscaling is based on surface elevation. + ! For other columns the downscaling is a simple copy (above). + + ! This is a simple downscaling procedure + ! Note that forc_hgt, forc_u, and forc_v are not downscaled. + + hsurf_g = grc_pp%elevation(g) ! gridcell sfc elevation + hsurf_t = top_pp%elevation(t) ! topounit sfc elevation + tbot_g = atm2lnd_vars%forc_t_not_downscaled_grc(g) ! atm sfc temp + thbot_g = atm2lnd_vars%forc_th_not_downscaled_grc(g) ! atm sfc pot temp + qbot_g = atm2lnd_vars%forc_q_not_downscaled_grc(g) ! atm sfc spec humid + pbot_g = atm2lnd_vars%forc_pbot_not_downscaled_grc(g) ! atm sfc pressure + zbot_g = atm2lnd_vars%forc_hgt_grc(g) ! atm ref height + + zbot_t = zbot_g + tbot_t = tbot_g-lapse_glcmec*(hsurf_t-hsurf_g) ! sfc temp for column + + Hbot = rair*0.5_r8*(tbot_g+tbot_t)/grav ! scale ht at avg temp + pbot_t = pbot_g*exp(-(hsurf_t-hsurf_g)/Hbot) ! column sfc press + + ! Derivation of potential temperature calculation: + ! + ! The textbook definition would be: + ! thbot_c = tbot_c * (p0/pbot_c)^(rair/cpair) + ! + ! Note that pressure is related to scale height as: + ! pbot_c = p0 * exp(-zbot_c/H) + ! + ! Using Hbot in place of H, we get: + ! pbot_c = p0 * exp(-zbot_c/Hbot) + ! + ! Plugging this in to the textbook definition, then manipulating, we get: + ! thbot_c = tbot_c * (p0/(p0*exp(-zbot_c/Hbot)))^(rair/cpair) + ! = tbot_c * (1/exp(-zbot_c/Hbot))^(rair/cpair) + ! = tbot_c * (exp(zbot_c/Hbot))^(rair/cpair) + ! = tbot_c * exp((zbot_c/Hbot) * (rair/cpair)) + + thbot_t= tbot_t*exp((zbot_t/Hbot)*(rair/cpair)) ! pot temp calc + + call Qsat(tbot_g,pbot_g,es_g,dum1,qs_g,dum2) + call Qsat(tbot_t,pbot_t,es_t,dum1,qs_t,dum2) + + qbot_t = qbot_g*(qs_t/qs_g) + egcm_t = qbot_t*pbot_t/(0.622_r8+0.378_r8*qbot_t) + rhos_t = (pbot_t-0.378_r8*egcm_t) / (rair*tbot_t) + + top_as%tbot(t) = tbot_t + top_as%thbot(t) = thbot_t + top_as%qbot(t) = qbot_t + top_as%pbot(t) = pbot_t + +! call check_downscale_consistency(bounds, atm2lnd_vars) + + end subroutine downscale_atm_state_to_topounit_cpl_bypass + + !------------------------------------------------------- + ! Downscale longwave radiation place holder + subroutine downscale_longwave_to_topounit_cpl_bypass(g, t, atm2lnd_vars) + + ! !DESCRIPTION: + ! Downscale longwave radiation from gridcell to column + ! Must be done AFTER temperature downscaling + + ! !USES: + use elm_time_manager, only : get_nstep + use domainMod , only : ldomain + use landunit_varcon , only : istice_mec + use elm_varcon , only : lapse_glcmec + use elm_varctl , only : glcmec_downscale_longwave + + ! !ARGUMENTS: + integer , intent(in) :: g + integer , intent(in) :: t + type(atm2lnd_type) , intent(in) :: atm2lnd_vars + + ! !LOCAL VARIABLES: + integer :: c,l,fc ! indices + integer :: nstep + real(r8) :: hsurf_t ! column-level elevation (m) + real(r8) :: hsurf_g ! gridcell-level elevation (m) + + real(r8) :: tair_g ! original gridcell mean air temperature + real(r8) :: tair_t ! downscaled topounit air temperature + real(r8) :: tsfc_g ! original gridcell surface temperature + real(r8) :: tsfc_t ! downscaled topounit surface temperature + real(r8) :: lwrad_g ! original gridcell mean LW radiation + real(r8) :: lwrad_t ! downscaled topounit LW radiation + real(r8) :: newsum_lwrad_g ! weighted sum of column-level lwrad after normalization + + character(len=*), parameter :: subname = 'downscale_longwave_to_topounit' + !----------------------------------------------------------------------- + + nstep = get_nstep() + + ! Do the downscaling + hsurf_g = grc_pp%elevation(g) + hsurf_t = top_pp%elevation(t) + + ! Here we assume that deltaLW = (dLW/dT)*(dT/dz)*deltaz + ! We get dLW/dT = 4*eps*sigma*T^3 = 4*LW/T from the Stefan-Boltzmann law, + ! evaluated at the mean temp. + ! We assume the same temperature lapse rate as above. + + tair_g = atm2lnd_vars%forc_t_not_downscaled_grc(g) + tair_t = top_as%tbot(t) + lwrad_g = atm2lnd_vars%forc_lwrad_not_downscaled_grc(g) + top_af%lwrad(t) = lwrad_g - & + 4.0_r8 * lwrad_g/(0.5_r8*(tair_t+tair_g)) * & + lapse_glcmec * (hsurf_t - hsurf_g) + + end subroutine downscale_longwave_to_topounit_cpl_bypass + end module lnd_downscale_atm_forcing diff --git a/components/elm/src/cpl/lnd_import_export.F90 b/components/elm/src/cpl/lnd_import_export.F90 index ad81846196d0..27c01e668184 100644 --- a/components/elm/src/cpl/lnd_import_export.F90 +++ b/components/elm/src/cpl/lnd_import_export.F90 @@ -1030,52 +1030,68 @@ subroutine lnd_import( bounds, x2l, atm2lnd_vars, glc2lnd_vars, lnd2atm_vars) end if end if - !set the topounit-level atmospheric state and flux forcings (bypass mode) + ! Adding topographic downscaling capability within CPL_BYPASS code block + ! PET, 7/12/2024 + if (use_atm_downscaling_to_topunit) then + atm2lnd_vars%forc_uovern = x2l(index_x2l_Sa_uovern,i) + atm2lnd_vars%forc_rain_not_downscaled_grc = forc_rainc + forc_rainl + atm2lnd_vars%forc_snow_not_downscaled_grc = forc_snowc + forc_snowl + + if(atm_gustiness) then + call endrun("Error: atm_gustiness not yet supported with multiple topounits (in CPL_BYPASS)") + end if + do topo = grc_pp%topi(g) , grc_pp%topf(g) + top_as%ugust(topo) = 0._r8 + end do + + call downscale_atm_forcing_to_topounit_cpl_bypass(g, atm2lnd_vars, lnd2atm_vars) + else + do topo = grc_pp%topi(g), grc_pp%topf(g) + top_as%tbot(topo) = atm2lnd_vars%forc_t_not_downscaled_grc(g) ! forc_txy Atm state K + top_as%thbot(topo) = atm2lnd_vars%forc_th_not_downscaled_grc(g) ! forc_thxy Atm state K + top_as%pbot(topo) = atm2lnd_vars%forc_pbot_not_downscaled_grc(g) ! ptcmxy Atm state Pa + top_as%qbot(topo) = atm2lnd_vars%forc_q_not_downscaled_grc(g) ! forc_qxy Atm state kg/kg + top_as%ubot(topo) = atm2lnd_vars%forc_u_grc(g) ! forc_uxy Atm state m/s + top_as%vbot(topo) = atm2lnd_vars%forc_v_grc(g) ! forc_vxy Atm state m/s + top_as%zbot(topo) = atm2lnd_vars%forc_hgt_grc(g) ! zgcmxy Atm state m + top_as%windbot(topo) = sqrt(top_as%ubot(topo)**2 + top_as%vbot(topo)**2) + ! Relative humidity (percent) + if (top_as%tbot(topo) > SHR_CONST_TKFRZ) then + e = esatw(tdc(top_as%tbot(topo))) + else + e = esati(tdc(top_as%tbot(topo))) + end if + qsat = 0.622_r8*e / (top_as%pbot(topo) - 0.378_r8*e) + top_as%rhbot(topo) = 100.0_r8*(top_as%qbot(topo) / qsat) + ! partial pressure of oxygen (Pa) + top_as%po2bot(topo) = o2_molar_const * top_as%pbot(topo) + ! air density (kg/m**3) - uses a temporary calculation of water vapor pressure (Pa) + vp = top_as%qbot(topo) * top_as%pbot(topo) / (0.622_r8 + 0.378_r8 * top_as%qbot(topo)) + top_as%rhobot(topo) = (top_as%pbot(topo) - 0.378_r8 * vp) / (rair * top_as%tbot(topo)) + top_af%rain(topo) = forc_rainc + forc_rainl ! sum of convective and large-scale rain + top_af%snow(topo) = forc_snowc + forc_snowl ! sum of convective and large-scale snow + top_af%solad(topo,2) = atm2lnd_vars%forc_solad_grc(g,2) ! forc_sollxy Atm flux W/m^2 + top_af%solad(topo,1) = atm2lnd_vars%forc_solad_grc(g,1) ! forc_solsxy Atm flux W/m^2 + top_af%solai(topo,2) = atm2lnd_vars%forc_solai_grc(g,2) ! forc_solldxy Atm flux W/m^2 + top_af%solai(topo,1) = atm2lnd_vars%forc_solai_grc(g,1) ! forc_solsdxy Atm flux W/m^2 + top_af%lwrad(topo) = atm2lnd_vars%forc_lwrad_not_downscaled_grc(g) ! flwdsxy Atm flux W/m^2 + ! derived flux forcings + top_af%solar(topo) = top_af%solad(topo,2) + top_af%solad(topo,1) + & + top_af%solai(topo,2) + top_af%solai(topo,1) + end do + end if + + !set the topounit-level atmospheric variables that are not handled in downscaling code do topo = grc_pp%topi(g), grc_pp%topf(g) ! first, all the state forcings - top_as%tbot(topo) = atm2lnd_vars%forc_t_not_downscaled_grc(g) ! forc_txy Atm state K - top_as%thbot(topo) = atm2lnd_vars%forc_th_not_downscaled_grc(g) ! forc_thxy Atm state K - top_as%pbot(topo) = atm2lnd_vars%forc_pbot_not_downscaled_grc(g) ! ptcmxy Atm state Pa - top_as%qbot(topo) = atm2lnd_vars%forc_q_not_downscaled_grc(g) ! forc_qxy Atm state kg/kg - top_as%ubot(topo) = atm2lnd_vars%forc_u_grc(g) ! forc_uxy Atm state m/s - top_as%vbot(topo) = atm2lnd_vars%forc_v_grc(g) ! forc_vxy Atm state m/s - if (implicit_stress) then + if (implicit_stress) then top_as%wsresp(topo) = 0._r8 ! Atm state m/s/Pa top_as%tau_est(topo) = 0._r8 ! Atm state Pa end if top_as%ugust(topo) = 0._r8 ! Atm state m/s - top_as%zbot(topo) = atm2lnd_vars%forc_hgt_grc(g) ! zgcmxy Atm state m - ! assign the state forcing fields derived from other inputs - ! Horizontal windspeed (m/s) - top_as%windbot(topo) = sqrt(top_as%ubot(topo)**2 + top_as%vbot(topo)**2) if (atm_gustiness) then top_as%windbot(topo) = sqrt(top_as%windbot(topo)**2 + top_as%ugust(topo)**2) end if - ! Relative humidity (percent) - if (top_as%tbot(topo) > SHR_CONST_TKFRZ) then - e = esatw(tdc(top_as%tbot(topo))) - else - e = esati(tdc(top_as%tbot(topo))) - end if - qsat = 0.622_r8*e / (top_as%pbot(topo) - 0.378_r8*e) - top_as%rhbot(topo) = 100.0_r8*(top_as%qbot(topo) / qsat) - ! partial pressure of oxygen (Pa) - top_as%po2bot(topo) = o2_molar_const * top_as%pbot(topo) - ! air density (kg/m**3) - uses a temporary calculation of water vapor pressure (Pa) - vp = top_as%qbot(topo) * top_as%pbot(topo) / (0.622_r8 + 0.378_r8 * top_as%qbot(topo)) - top_as%rhobot(topo) = (top_as%pbot(topo) - 0.378_r8 * vp) / (rair * top_as%tbot(topo)) - - ! second, all the flux forcings - top_af%rain(topo) = forc_rainc + forc_rainl ! sum of convective and large-scale rain - top_af%snow(topo) = forc_snowc + forc_snowl ! sum of convective and large-scale snow - top_af%solad(topo,2) = atm2lnd_vars%forc_solad_grc(g,2) ! forc_sollxy Atm flux W/m^2 - top_af%solad(topo,1) = atm2lnd_vars%forc_solad_grc(g,1) ! forc_solsxy Atm flux W/m^2 - top_af%solai(topo,2) = atm2lnd_vars%forc_solai_grc(g,2) ! forc_solldxy Atm flux W/m^2 - top_af%solai(topo,1) = atm2lnd_vars%forc_solai_grc(g,1) ! forc_solsdxy Atm flux W/m^2 - top_af%lwrad(topo) = atm2lnd_vars%forc_lwrad_not_downscaled_grc(g) ! flwdsxy Atm flux W/m^2 - ! derived flux forcings - top_af%solar(topo) = top_af%solad(topo,2) + top_af%solad(topo,1) + & - top_af%solai(topo,2) + top_af%solai(topo,1) end do !----------------------------------------------------------------------------------------------------- diff --git a/components/elm/src/data_types/VegetationPropertiesType.F90 b/components/elm/src/data_types/VegetationPropertiesType.F90 index 12f9a7b22ffa..d62167301cd4 100644 --- a/components/elm/src/data_types/VegetationPropertiesType.F90 +++ b/components/elm/src/data_types/VegetationPropertiesType.F90 @@ -43,7 +43,7 @@ module VegetationPropertiesType real(r8), pointer :: dsladlai (:) => null() ! dSLA/dLAI, projected area basis [m^2/gC] real(r8), pointer :: leafcn (:) => null() ! leaf C:N (gC/gN) real(r8), pointer :: flnr (:) => null() ! fraction of leaf N in the Rubisco enzyme (gN Rubisco / gN leaf) - real(r8), pointer :: woody (:) => null() ! binary flag for woody lifeform (1=woody, 0=not woody) + real(r8), pointer :: woody (:) => null() ! woody lifeform flag (0 = non-woody, 1 = tree, 2 = shrub) real(r8), pointer :: lflitcn (:) => null() ! leaf litter C:N (gC/gN) real(r8), pointer :: frootcn (:) => null() ! fine root C:N (gC/gN) real(r8), pointer :: livewdcn (:) => null() ! live wood (phloem and ray parenchyma) C:N (gC/gN) @@ -115,7 +115,7 @@ module VegetationPropertiesType real(r8), pointer :: lamda_ptase => null()! critical value that incur biochemical production real(r8), pointer :: i_vc(:) => null() ! intercept of photosynthesis vcmax ~ leaf n content regression model real(r8), pointer :: s_vc(:) => null() ! slope of photosynthesis vcmax ~ leaf n content regression model - real(r8), pointer :: nsc_rtime(:) => null() ! non-structural carbon residence time + real(r8), pointer :: nsc_rtime(:) => null() ! non-structural carbon residence time real(r8), pointer :: pinit_beta1(:) => null() ! shaping parameter for P initialization real(r8), pointer :: pinit_beta2(:) => null() ! shaping parameter for P initialization real(r8), pointer :: alpha_nfix(:) => null() ! fraction of fixed N goes directly to plant @@ -151,6 +151,11 @@ module VegetationPropertiesType real(r8), pointer :: needleleaf(:) => null() !needleleaf or broadleaf real(r8), pointer :: nfixer(:) => null() !cablity of nitrogen fixation from atm. N2 + ! NGEE Arctic snow-vegetation interactions + real(r8), pointer :: bendresist(:) ! vegetation resistance to bending under snow loading, 0 to 1 (e.g., Liston and Hiemstra 2011) + real(r8), pointer :: vegshape(:) ! shape parameter to modify shrub burial by snow (1 = parabolic, 2 = hemispheric) + real(r8), pointer :: stocking(:) ! stocking density for pft (stems / hectare) + real(r8), pointer :: taper(:) ! ratio of height:radius_breast_height (woody vegetation allometry) contains procedure, public :: Init => veg_vp_init @@ -188,8 +193,10 @@ subroutine veg_vp_init(this) use pftvarcon , only : fnr, act25, kcha, koha, cpha, vcmaxha, jmaxha, tpuha use pftvarcon , only : lmrha, vcmaxhd, jmaxhd, tpuhd, lmrse, qe, theta_cj use pftvarcon , only : bbbopt, mbbopt, nstor, br_xr, tc_stress, lmrhd - ! new properties for flexible PFT + ! new properties for flexible PFT (NGEE Arctic IM4) use pftvarcon , only : climatezone, nonvascular, graminoid, iscft,needleleaf, nfixer + ! snow/vegetation interactions (NGEE Arctic IM3) + use pftvarcon , only : bendresist, stocking, vegshape, taper ! class (vegetation_properties_type) :: this @@ -322,6 +329,11 @@ subroutine veg_vp_init(this) allocate( this%needleleaf(0:numpft)) ; this%needleleaf(:) =spval allocate( this%nfixer(0:numpft)) ; this%nfixer(:) =spval ! ----------------------------------------------------------------------------------------------------------- + ! NGEE Arctic snow-vegetation interactions + allocate(this%bendresist(0:numpft)) ; this%bendresist(:) =spval + allocate(this%vegshape(0:numpft)) ; this%vegshape(:) =spval + allocate(this%stocking(0:numpft)) ; this%stocking(:) =spval + allocate(this%taper(0:numpft)) ; this%taper(:) =spval do m = 0,numpft @@ -472,6 +484,13 @@ subroutine veg_vp_init(this) this%lamda_ptase = lamda_ptase this%tc_stress = tc_stress + ! NGEE Arctic - snow/vegetation interactions + do m = 0, numpft ! RPF - move up to earlier pft loops? + this%bendresist(m) = bendresist(m) + this%vegshape(m) = vegshape(m) + this%stocking(m) = stocking(m) + this%taper(m) = taper(m) + end do end subroutine veg_vp_init end module VegetationPropertiesType diff --git a/components/elm/src/main/atm2lndType.F90 b/components/elm/src/main/atm2lndType.F90 index a00140ab2dd6..f2facc776224 100644 --- a/components/elm/src/main/atm2lndType.F90 +++ b/components/elm/src/main/atm2lndType.F90 @@ -146,6 +146,10 @@ module atm2lndType real(r8) , pointer :: t_mo_patch (:) => null() ! patch 30-day average temperature (Kelvin) real(r8) , pointer :: t_mo_min_patch (:) => null() ! patch annual min of t_mo (Kelvin) + ! Needed for FNM precip downscaling, when used within CPL_BYPASS + real(r8), pointer :: forc_uovern (:) => null() ! Froude number (dimensionless) + + contains procedure, public :: Init @@ -310,6 +314,7 @@ subroutine InitAllocate(this, bounds) allocate(this%forc_ndep_nitr_grc (begg:endg)) ; this%forc_ndep_nitr_grc (:) = ival allocate(this%forc_soilph_grc (begg:endg)) ; this%forc_soilph_grc (:) = ival end if + allocate(this%forc_uovern (begg:endg)) ; this%forc_uovern (:) = ival end subroutine InitAllocate diff --git a/components/elm/src/main/pftvarcon.F90 b/components/elm/src/main/pftvarcon.F90 index fc5bfd54e7d3..66c2adda90f8 100644 --- a/components/elm/src/main/pftvarcon.F90 +++ b/components/elm/src/main/pftvarcon.F90 @@ -29,7 +29,7 @@ module pftvarcon ! ! Vegetation type constants ! - integer :: noveg !value for not vegetated + integer :: noveg !value for not vegetated integer :: ndllf_evr_tmp_tree !value for Needleleaf evergreen temperate tree integer :: ndllf_evr_brl_tree !value for Needleleaf evergreen boreal tree integer :: ndllf_dcd_brl_tree !value for Needleleaf deciduous boreal tree @@ -76,7 +76,7 @@ module pftvarcon integer :: nrtubers !value for root tubers, rain fed (rf) integer :: nrtubersirrig !value for root tubers, irrigated (ir) integer :: nsugarcane !value for sugarcane, rain fed (rf) - integer :: nsugarcaneirrig !value for sugarcane, irrigated (ir) + integer :: nsugarcaneirrig !value for sugarcane, irrigated (ir) integer :: nmiscanthus !value for miscanthus, rain fed (rf) integer :: nmiscanthusirrig !value for miscanthus, irrigated (ir) integer :: nswitchgrass !value for switchgrass, rain fed (rf) @@ -84,8 +84,8 @@ module pftvarcon integer :: npoplar !value for poplar, rain fed (rf) integer :: npoplarirrig !value for poplar, irrigated (ir) integer :: nwillow !value for willow, rain fed (rf) - integer :: nwillowirrig !value for willow, irrigated (ir) - + integer :: nwillowirrig !value for willow, irrigated (ir) + ! Number of crop functional types actually used in the model. This includes each CFT for ! which is_pft_known_to_model is true. Note that this includes irrigated crops even if ! irrigation is turned off in this run: it just excludes crop types that aren't handled @@ -171,7 +171,7 @@ module pftvarcon real(r8), allocatable :: minplanttemp(:) !mininum planting temperature used in Phenology (K) real(r8), allocatable :: senestemp(:) !senescence temperature for perennial crops used in Phenology (K) real(r8), allocatable :: min_days_senes(:) !minimum leaf age to allow for leaf senescence - real(r8), allocatable :: froot_leaf(:) !allocation parameter: new fine root C per new leaf C (gC/gC) + real(r8), allocatable :: froot_leaf(:) !allocation parameter: new fine root C per new leaf C (gC/gC) real(r8), allocatable :: stem_leaf(:) !allocation parameter: new stem c per new leaf C (gC/gC) real(r8), allocatable :: croot_stem(:) !allocation parameter: new coarse root C per new stem C (gC/gC) real(r8), allocatable :: flivewd(:) !allocation parameter: fraction of new wood that is live (phloem and ray parenchyma) (no units) @@ -215,7 +215,7 @@ module pftvarcon real(r8), allocatable :: fyield(:) !fraction of grain that is actually harvested real(r8), allocatable :: root_dmx(:) !maximum root depth - integer, parameter :: pftname_len = 40 ! max length of pftname + integer, parameter :: pftname_len = 40 ! max length of pftname character(len=:), allocatable :: pftname(:) !PFT description real(r8), parameter :: reinickerp = 1.6_r8 !parameter in allometric equation @@ -276,7 +276,7 @@ module pftvarcon real(r8), allocatable :: deadwdcp_obs_flex(:,:) !upper and lower range of dead wood (xylem and heartwood) C:P (gC/gP) ! Photosynthesis parameters real(r8), allocatable :: fnr(:) !fraction of nitrogen in RuBisCO - real(r8), allocatable :: act25(:) + real(r8), allocatable :: act25(:) real(r8), allocatable :: kcha(:) !Activation energy for kc real(r8), allocatable :: koha(:) !Activation energy for ko real(r8), allocatable :: cpha(:) !Activation energy for cp @@ -293,7 +293,7 @@ module pftvarcon real(r8), allocatable :: theta_cj(:) ! real(r8), allocatable :: bbbopt(:) !Ball-Berry stomatal conductance intercept real(r8), allocatable :: mbbopt(:) !Ball-Berry stomatal conductance slope - real(r8), allocatable :: nstor(:) !Nitrogen storage pool timescale + real(r8), allocatable :: nstor(:) !Nitrogen storage pool timescale real(r8), allocatable :: br_xr(:) !Base rate for excess respiration real(r8) :: tc_stress !Critial temperature for moisture stress real(r8), allocatable :: vcmax_np1(:) !vcmax~np relationship coefficient @@ -312,6 +312,13 @@ module pftvarcon real(r8), allocatable :: gcbr_p(:) !effectiveness of roots in reducing rainfall-driven erosion real(r8), allocatable :: gcbr_q(:) !effectiveness of roots in reducing runoff-driven erosion + ! NGEE Arctic snow-vegetation interactions + real(r8), allocatable :: bendresist(:) ! vegetation resistance to bending under snow loading, 0 to 1 (e.g., Liston and Hiemstra 2011; Sturm et al. 2005) + real(r8), allocatable :: vegshape(:) ! shape parameter to modify shrub burial by snow (1 = parabolic, 2 = hemispheric) + real(r8), allocatable :: stocking(:) ! stocking density for pft (stems / hectare) + real(r8), allocatable :: taper(:) ! ratio of height:radius_breast_height (woody vegetation allometry) + logical :: taper_defaults ! set flag to use taper defaults if not on params file (necessary as import and set values are in different places) + ! new pft properties, together with woody, crop, percrop, evergreen, stress_decid, season_decid, defined above, ! are introduced to define vegetation properties. This will be well defineing a pft so that no indices needed for codes. real(r8), allocatable :: climatezone(:) ! distributed climate zone (0 = any zone, 1 = tropical, 2 = temperate, 3 = boreal, 4 = arctic) @@ -379,7 +386,7 @@ subroutine pftconrd ! and finally crops, ending with soybean ! DO NOT CHANGE THE ORDER -- WITHOUT MODIFYING OTHER PARTS OF THE CODE WHERE THE ORDER MATTERS! ! - character(len=pftname_len) :: expected_pftnames(0:mxpft) + character(len=pftname_len) :: expected_pftnames(0:mxpft) !----------------------------------------------------------------------- expected_pftnames( 0) = 'not_vegetated ' @@ -452,95 +459,95 @@ subroutine pftconrd ! but cannot be less, i.e. there cannot be more pfts in the surface file than there are set of pft parameters, if running the non-crop version of model. if (npft mxpft_nc) EXIT ! exit the do loop - + ! (FATES-INTERF) Later, depending on how the team plans to structure the crop model ! or other modules that co-exist while FATES is on, we may want to preserve these pft definitions ! on non-fates columns. For now, they are incompatible, and this check is warranted (rgk 04-2017) @@ -1507,6 +1542,17 @@ subroutine pftconrd end do end if + ! reassign taper values for shrubs - RPF + if (taper_defaults) then + do i = 0, npft-1 + if (woody(i) == 2._r8) then + taper(i) = 10._r8 ! shrubs + else if (woody(i) == 1._r8) then + taper(i) = 200._r8 + end if + end do + end if + if (masterproc) then write(iulog,*) 'Successfully read PFT physiological data' write(iulog,*)