diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index e859a7339f..d41a4b98b4 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1977,6 +1977,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) bc_out(s)%z0m_pa(ifp) = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp) bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp) bc_out(s)%dleaf_pa(ifp) = EDPftvarcon_inst%dleaf(1) + bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp) = 1 endif ! ----------------------------------------------------------------------------- @@ -2002,6 +2003,12 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label + if(currentPatch%nocomp_pft_label.gt.0)then + bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp) = EDPftvarcon_inst%voc_pftindex(currentPatch%nocomp_pft_label) + else + bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp) = 1 ! dummy for bare ground. + endif + ! Calculate area indices for output boundary to HLM ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index b37e3630a0..5f4e4a4eff 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -23,6 +23,7 @@ module FATESPlantRespPhotosynthMod use FatesGlobals, only : endrun => fates_endrun use FatesGlobals, only : fates_log use FatesGlobals, only : FatesWarn,N2S,A2S,I2S + use FatesInterfaceTypesMod , only : hlm_use_nocomp use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : itrue use FatesConstantsMod, only : nearzero @@ -193,6 +194,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! net leaf photosynthesis averaged over sun and shade leaves. [umol CO2/m**2/s] real(r8) :: anet_av_z(nlevleaf,maxpft,nclmax) + ! internal leaf Co2 memory variable (for passing into MEGAN) Pa + real(r8) :: internal_co2_z(nlevleaf,maxpft,nclmax) + ! Photosynthesis [umol /m2 /s] real(r8) :: psn_z(nlevleaf,maxpft,nclmax) @@ -355,9 +359,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! --------------------------------------------------------------------------- bc_out(s)%rssun_pa(ifp) = 0._r8 bc_out(s)%rssha_pa(ifp) = 0._r8 - - psn_z(:,:,:) = 0._r8 + bc_out(s)%ci_pa(ifp) = 0._r8 + psn_z(:,:,:) = 0._r8 + !internal_co2_z(:,:,:) = 0._r8 g_sb_leaves = 0._r8 patch_la = 0._r8 @@ -731,6 +736,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) psn_z(iv,ft,cl), & ! out rs_z(iv,ft,cl), & ! out anet_av_z(iv,ft,cl), & ! out + internal_co2_z(iv,ft,cl), & ! out c13disc_z(cl,ft,iv)) ! out rate_mask_z(iv,ft,cl) = .true. @@ -1080,7 +1086,13 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! when it comes time to calculate a flux rate per unit ground bc_out(s)%rssun_pa(ifp) = r_stomata bc_out(s)%rssha_pa(ifp) = r_stomata - + + if(hlm_use_nocomp)then + bc_out(s)%ci_pa(ifp) = internal_co2_z(1,currentpatch%nocomp_pft_label,1) + else + bc_out(s)%ci_pa(ifp) = -999 + endif + ! This value is used for diagnostics, the molar form of conductance ! is what is used in the field usually, so we track that form currentPatch%c_stomata = cf / r_stomata @@ -1204,6 +1216,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in psn_out, & ! out rstoma_out, & ! out anet_av_out, & ! out + internal_co2_out, & ! out c13disc_z) ! out @@ -1263,6 +1276,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in real(r8), intent(out) :: rstoma_out ! stomatal resistance (1/gs_lsl) (s/m) real(r8), intent(out) :: anet_av_out ! net leaf photosynthesis (umol CO2/m**2/s) ! averaged over sun and shade leaves. + real(r8), intent(out) :: internal_co2_out ! internal co2 diagnostic (Pa) real(r8), intent(out) :: c13disc_z ! carbon 13 in newly assimilated carbon @@ -1354,6 +1368,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in anet_av_out = -lmr psn_out = 0._r8 + internal_co2_out = init_co2_inter_c ! The cuticular conductance already factored in maximum resistance as a bound ! no need to re-bound it @@ -1372,6 +1387,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. anet_av_out = 0._r8 + internal_co2_out = 0._r8 gstoma = 0._r8 do sunsha = 1,2 @@ -1599,10 +1615,12 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in if(sunsha == 1)then !sunlit psn_out = psn_out + agross * f_sun_lsl anet_av_out = anet_av_out + anet * f_sun_lsl + internal_co2_out = internal_co2_out + co2_inter_c * f_sun_lsl gstoma = gstoma + 1._r8/(min(1._r8/gs, rsmax0)) * f_sun_lsl else psn_out = psn_out + agross * (1.0_r8-f_sun_lsl) anet_av_out = anet_av_out + anet * (1.0_r8-f_sun_lsl) + internal_co2_out = internal_co2_out + co2_inter_c * (1.0_r8 -f_sun_lsl) gstoma = gstoma + & 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) end if @@ -1648,7 +1666,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in psn_out = 0._r8 anet_av_out = 0._r8 - + internal_co2_out = 0._r8 rstoma_out = min(rsmax0,cf/(stem_cuticle_loss_frac*stomatal_intercept(ft))) c13disc_z = 0.0_r8 diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 3d91f2f1b9..07c7922751 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -147,7 +147,8 @@ module EDPftvarcon real(r8), allocatable :: rhos(:, :) ! Stem reflectance; second dim: 1 = vis, 2 = nir real(r8), allocatable :: taul(:, :) ! Leaf transmittance; second dim: 1 = vis, 2 = nir real(r8), allocatable :: taus(:, :) ! Stem transmittance; second dim: 1 = vis, 2 = nir - + real(r8), allocatable :: voc_pftindex(:) ! Index for MEGAN parameters + ! Fire Parameters (No PFT vector capabilities in their own routines) ! See fire/SFParamsMod.F90 for bulk of fire parameters ! ------------------------------------------------------------------------------------------- @@ -454,6 +455,10 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_voc_pftindex' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_nonhydro_smpso' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -921,6 +926,10 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=this%c3psn) + name = 'fates_voc_pftindex' + call fates_params%RetrieveParameterAllocate(name=name, & + data=this%voc_pftindex) + name = 'fates_nonhydro_smpso' call fates_params%RetrieveParameterAllocate(name=name, & data=this%smpso) @@ -1728,6 +1737,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'xl = ',EDPftvarcon_inst%xl write(fates_log(),fmt0) 'clumping_index = ',EDPftvarcon_inst%clumping_index write(fates_log(),fmt0) 'c3psn = ',EDPftvarcon_inst%c3psn + write(fates_log(),fmt0) 'voc_pftindex = ',EDPftvarcon_inst%voc_pftindex write(fates_log(),fmt0) 'vcmax25top = ',EDPftvarcon_inst%vcmax25top write(fates_log(),fmt0) 'smpso = ',EDPftvarcon_inst%smpso write(fates_log(),fmt0) 'smpsc = ',EDPftvarcon_inst%smpsc @@ -2229,7 +2239,12 @@ subroutine FatesCheckParams(is_master) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + write(*,*) EDPftvarcon_inst%voc_pftindex(:) + if(EDPftvarcon_inst%voc_pftindex(ipft) .le. 0 .or. EDPftvarcon_inst%voc_pftindex(ipft) .gt. 16 ) then + write(fates_log(),*) 'MEGAN indices must be between 1 and 16',ipft,EDPftvarcon_inst%voc_pftindex(ipft) + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + if( hlm_use_fixed_biogeog .eq. itrue ) then ! check that the host-fates PFT map adds to one along HLM dimension so that all the HLM area ! goes to a FATES PFT. Each FATES PFT can get < or > 1 of an HLM PFT. diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 12eb83bcc0..7e2e017d87 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -325,6 +325,7 @@ subroutine zero_bcs(fates,s) ! Output boundaries fates%bc_out(s)%active_suction_sl(:) = .false. fates%bc_out(s)%fsun_pa(:) = 0.0_r8 + fates%bc_out(s)%ci_pa(:) = 0.0_r8 fates%bc_out(s)%laisun_pa(:) = 0.0_r8 fates%bc_out(s)%laisha_pa(:) = 0.0_r8 fates%bc_out(s)%rootr_pasl(:,:) = 0.0_r8 @@ -391,6 +392,8 @@ subroutine zero_bcs(fates,s) fates%bc_out(s)%z0m_pa(:) = 0.0_r8 fates%bc_out(s)%dleaf_pa(:) = 0.0_r8 fates%bc_out(s)%nocomp_pft_label_pa(:) = 0 + fates%bc_out(s)%nocomp_MEGAN_pft_label_pa(:) = 0 + fates%bc_out(s)%canopy_fraction_pa(:) = 0.0_r8 fates%bc_out(s)%frac_veg_nosno_alb_pa(:) = 0.0_r8 @@ -605,6 +608,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in) ! Radiation allocate(bc_out%fsun_pa(maxpatch_total)) + allocate(bc_out%ci_pa(maxpatch_total)) allocate(bc_out%laisun_pa(maxpatch_total)) allocate(bc_out%laisha_pa(maxpatch_total)) @@ -718,6 +722,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in) allocate(bc_out%frac_veg_nosno_alb_pa(maxpatch_total)) allocate(bc_out%nocomp_pft_label_pa(maxpatch_total)) + allocate(bc_out%nocomp_MEGAN_pft_label_pa(maxpatch_total)) ! Plant-Hydro BC's if (hlm_use_planthydro.eq.itrue) then diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 75e64307a5..0b84d89605 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -600,6 +600,10 @@ module FatesInterfaceTypesMod ! Shaded canopy LAI real(r8),allocatable :: laisha_pa(:) + + ! Average internal CO2 concentration + real(r8), allocatable :: ci_pa(:) + ! Logical stating whether a soil layer can have water uptake by plants ! The only condition right now is that liquid water exists @@ -750,6 +754,8 @@ module FatesInterfaceTypesMod integer, allocatable :: nocomp_pft_label_pa(:) ! in nocomp and SP mode, each patch has a PFT identity. + integer, allocatable :: nocomp_MEGAN_pft_label_pa(:) ! Index to map from FATES NOCOMP PFT identity into MEGAN PFT space. + ! FATES Hydraulics diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 1f77e53e11..90a2e3f3cb 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -672,6 +672,9 @@ variables: double fates_woody(fates_pft) ; fates_woody:units = "logical flag" ; fates_woody:long_name = "Binary woody lifeform flag" ; + double fates_voc_pftindex(fates_pft) ; + fates_voc_pftindex: units = "Integer PFT index" ; + fates_voc_pftindex: long_name = "integer index for MEGAN PFT definitions" ; double fates_hlm_pft_map(fates_hlm_pftno, fates_pft) ; fates_hlm_pft_map:units = "area fraction" ; fates_hlm_pft_map:long_name = "In fixed biogeog mode, fraction of HLM area associated with each FATES PFT" ; @@ -1669,7 +1672,9 @@ data: fates_wood_density = 0.548327, 0.44235, 0.454845, 0.754336, 0.548327, 0.566452, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7 ; - fates_woody = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0 ; + fates_woody = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0 ; + + fates_voc_pftindex = 4, 1, 3, 5, 6, 7, 9, 10, 10, 9, 11, 12, 13, 14 ; fates_hlm_pft_map = 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,