diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 86d3bbeddb..65752ddca1 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -9,6 +9,7 @@ module FatesHistoryInterfaceMod use FatesConstantsMod , only : mg_per_kg use FatesConstantsMod , only : pi_const use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : rsnbl_math_prec use FatesConstantsMod , only : t_water_freeze_k_1atm use FatesConstantsMod , only : n_term_mort_types use FatesConstantsMod , only : i_term_mort_type_cstarv @@ -138,6 +139,7 @@ module FatesHistoryInterfaceMod use FatesSizeAgeTypeIndicesMod, only : get_cdamagesize_class_index use FatesSizeAgeTypeIndicesMod, only : get_cdamagesizepft_class_index use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index + use FatesInterfaceTypesMod , only : hlm_use_luh implicit none private ! By default everything is private @@ -362,6 +364,12 @@ module FatesHistoryInterfaceMod integer :: ih_burnedarea_si_landuse integer :: ih_gpp_si_landuse integer :: ih_npp_si_landuse + integer :: ih_tveg_si_landuse + integer :: ih_tsa_si_landuse + integer :: ih_sw_abs_si_landuse + integer :: ih_lw_net_si_landuse + integer :: ih_shflux_si_landuse + integer :: ih_lhflux_si_landuse ! land use by land use variables integer :: ih_disturbance_rate_si_lulu @@ -838,6 +846,7 @@ module FatesHistoryInterfaceMod procedure :: update_history_hifrq_sitelevel procedure :: update_history_hifrq_subsite procedure :: update_history_hifrq_subsite_ageclass + procedure :: update_history_hifrq_landuse procedure :: update_history_hydraulics procedure :: update_history_nutrflux @@ -5007,6 +5016,9 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) if(hlm_hist_level_hifrq>0) then call update_history_hifrq_sitelevel(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) + if (hlm_use_luh .eq. itrue) then + call update_history_hifrq_landuse(this,nc,nsites,sites,bc_in,dt_tstep) + end if if(hlm_hist_level_hifrq>1) then call update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) call update_history_hifrq_subsite_ageclass(this,nsites,sites,dt_tstep) @@ -5240,6 +5252,161 @@ end subroutine update_history_hifrq_sitelevel ! =============================================================================================== + subroutine update_history_hifrq_landuse(this,nc,nsites,sites,bc_in,dt_tstep) + + ! + ! Arguments + class(fates_history_interface_type) :: this + integer , intent(in) :: nc ! clump index + integer , intent(in) :: nsites + type(ed_site_type) , intent(inout), target :: sites(nsites) + type(bc_in_type) , intent(in) :: bc_in(nsites) + real(r8) , intent(in) :: dt_tstep + + ! Locals + integer :: s ! The local site index + integer :: io_si ! The site index of the IO array + + real(r8) :: landuse_statevector(n_landuse_cats) + real(r8) :: canopy_area_bylanduse(n_landuse_cats) + integer :: i_lu + logical :: foundbaregroundpatch + + type(fates_patch_type),pointer :: cpatch + type(fates_cohort_type),pointer :: ccohort + real(r8) :: dt_tstep_inv ! Time step in frequency units (/s) + real(r8) :: vegarea_per_patcharea ! temporary weighting variable (unitless) + + associate( hio_tveg_si_landuse => this%hvars(ih_tveg_si_landuse)%r82d,& + hio_gpp_si_landuse => this%hvars(ih_gpp_si_landuse)%r82d, & + hio_tsa_si_landuse => this%hvars(ih_tsa_si_landuse)%r82d,& + hio_sw_abs_si_landuse => this%hvars(ih_sw_abs_si_landuse)%r82d,& + hio_lw_net_si_landuse => this%hvars(ih_lw_net_si_landuse)%r82d,& + hio_shflux_si_landuse => this%hvars(ih_shflux_si_landuse)%r82d,& + hio_lhflux_si_landuse => this%hvars(ih_lhflux_si_landuse)%r82d) + + do_sites: do s = 1,nsites + + io_si = sites(s)%h_gid + + dt_tstep_inv = 1.0_r8/dt_tstep + + ! biophysical properties that are indexed by land use + landuse_statevector(:) = sites(s)%get_current_landuse_statevector() * AREA + + ! get the total canopy area for each land use type + canopy_area_bylanduse(:) = 0._r8 + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + if ( cpatch%land_use_label .ne. nocomp_bareground_land ) then + canopy_area_bylanduse(cpatch%land_use_label) = canopy_area_bylanduse(cpatch%land_use_label) + & + cpatch%total_canopy_area + endif + cpatch => cpatch%younger + end do + + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + if (cpatch%total_canopy_area .gt. rsnbl_math_prec) then + ! for TVEG, since it is only defined on vegetated area of vegetated patches, normalize by the total vegetated area + hio_tveg_si_landuse(io_si,cpatch%land_use_label) = hio_tveg_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%t_veg_pa(cpatch%patchno) * cpatch%total_canopy_area/canopy_area_bylanduse(cpatch%land_use_label) + + ! for the rest of these, first weight by the vegetated area of each patch over the total patch area for each land use type + vegarea_per_patcharea = cpatch%total_canopy_area/landuse_statevector(cpatch%land_use_label) + + hio_tsa_si_landuse(io_si,cpatch%land_use_label) = hio_tsa_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%t2m_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_sw_abs_si_landuse(io_si,cpatch%land_use_label) = hio_sw_abs_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%swabs_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_lw_net_si_landuse(io_si,cpatch%land_use_label) = hio_lw_net_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%netlw_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_shflux_si_landuse(io_si,cpatch%land_use_label) = hio_shflux_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%shflux_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_lhflux_si_landuse(io_si,cpatch%land_use_label) = hio_lhflux_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%lhflux_pa(cpatch%patchno) * vegarea_per_patcharea + endif + cpatch => cpatch%younger + end do + + ! for all the land-use indexed variables, except for TVEG, also add in the component for the unvegetated area of each land use + foundbaregroundpatch = .false. + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + if (cpatch%land_use_label .eq. nocomp_bareground_land .and. .not. foundbaregroundpatch) then + foundbaregroundpatch = .true. + do i_lu = 1, n_landuse_cats + if ( landuse_statevector(i_lu) .gt. rsnbl_math_prec ) then + hio_tsa_si_landuse(io_si,i_lu) = hio_tsa_si_landuse(io_si,i_lu) + & + bc_in(s)%t2m_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_sw_abs_si_landuse(io_si,i_lu) = hio_sw_abs_si_landuse(io_si,i_lu) + & + bc_in(s)%swabs_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_lw_net_si_landuse(io_si,i_lu) = hio_lw_net_si_landuse(io_si,i_lu) + & + bc_in(s)%netlw_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_shflux_si_landuse(io_si,i_lu) = hio_shflux_si_landuse(io_si,i_lu) + & + bc_in(s)%shflux_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_lhflux_si_landuse(io_si,i_lu) = hio_lhflux_si_landuse(io_si,i_lu) + & + bc_in(s)%lhflux_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + end if + end do + end if + cpatch => cpatch%younger + end do + + ! instead of leaving the values for unoccupied areas as zero, set as missing values + do i_lu = 1, n_landuse_cats + + ! if a given land use type is not present, set the value as missing + if ( landuse_statevector(i_lu) .le. rsnbl_math_prec ) then + hio_tsa_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_sw_abs_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_lw_net_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_shflux_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_lhflux_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + end if + + ! for tveg, ignore if there is no vegetation present on any patches of a given land use type + if ( canopy_area_bylanduse(i_lu) .le. rsnbl_math_prec ) then + hio_tveg_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + end if + + end do + + ! for GPP by land use, we need to loop over both patches and cohorts + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + ccohort => cpatch%shortest + do while(associated(ccohort)) + if ( (.not. ccohort%isnew) .and. (cpatch%land_use_label .gt. nocomp_bareground_land) ) then + hio_gpp_si_landuse(io_si,cpatch%land_use_label) = hio_gpp_si_landuse(io_si,cpatch%land_use_label) & + + ccohort%gpp_tstep * ccohort%n * dt_tstep_inv + end if + ccohort => ccohort%taller + end do + cpatch => cpatch%younger + end do + + end do do_sites + + end associate + return + end subroutine update_history_hifrq_landuse + + ! =============================================================================================== + subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) ! --------------------------------------------------------------------------------- @@ -5298,7 +5465,6 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tst hio_froot_mr_understory_si_scls => this%hvars(ih_froot_mr_understory_si_scls)%r82d, & hio_resp_g_understory_si_scls => this%hvars(ih_resp_g_understory_si_scls)%r82d, & hio_resp_m_understory_si_scls => this%hvars(ih_resp_m_understory_si_scls)%r82d, & - hio_gpp_si_landuse => this%hvars(ih_gpp_si_landuse)%r82d, & hio_parsun_z_si_cnlf => this%hvars(ih_parsun_z_si_cnlf)%r82d, & hio_parsha_z_si_cnlf => this%hvars(ih_parsha_z_si_cnlf)%r82d, & hio_ts_net_uptake_si_cnlf => this%hvars(ih_ts_net_uptake_si_cnlf)%r82d, & @@ -5387,11 +5553,6 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,bc_out,dt_tst hio_ar_frootm_si_scpf(io_si,scpf) = hio_ar_frootm_si_scpf(io_si,scpf) + & ccohort%froot_mr * n_perm2 - if (cpatch%land_use_label .gt. nocomp_bareground_land) then - hio_gpp_si_landuse(io_si,cpatch%land_use_label) = hio_gpp_si_landuse(io_si,cpatch%land_use_label) & - + ccohort%gpp_tstep * ccohort%n * dt_tstep_inv - end if - ! accumulate fluxes on canopy- and understory- separated fluxes if (ccohort%canopy_layer .eq. 1) then @@ -8743,6 +8904,51 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & ivar=ivar, initialize=initialize_variables, index = ih_tveg_si ) + if (hlm_use_luh .eq. itrue) then + ! biophysics variables that are indexed by land use type + call this%set_history_var(vname='FATES_TVEG_LU', units='degrees Kelvin', & + long='fates instantaneous mean vegetation temperature by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_tveg_si_landuse ) + + call this%set_history_var(vname='FATES_TSA_LU', units='degrees Kelvin', & + long='fates instantaneous mean near-surface (2m) air temperature by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_tsa_si_landuse ) + + call this%set_history_var(vname='FATES_SWABS_LU', units='W m-2', & + long='fates absorbed shortwave radiation by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_sw_abs_si_landuse ) + + call this%set_history_var(vname='FATES_NETLW_LU', units='W m-2', & + long='fates net longwave flux by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_lw_net_si_landuse ) + + call this%set_history_var(vname='FATES_SHFLUX_LU', units='W m-2', & + long='fates sensible heat flux by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_shflux_si_landuse ) + + call this%set_history_var(vname='FATES_LHFLUX_LU', units='W m-2', & + long='fates latent heat flux by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_lhflux_si_landuse ) + + call this%set_history_var(vname='FATES_GPP_LU', units='kg m-2 s-1', & + long='gross primary productivity by land use type in kg carbon per m2 per second', & + use_default='inactive', avgflag='A', vtype=site_landuse_r8, & + hlms='CLM:ALM', upfreq=group_hifr_simple, ivar=ivar, initialize=initialize_variables, & + index = ih_gpp_si_landuse) + endif + call this%set_history_var(vname='FATES_VIS_RAD_ERROR', units='-', & long='mean two-stream solver error for VIS', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & @@ -8902,12 +9108,6 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=group_hifr_complx, ivar=ivar, initialize=initialize_variables, & index = ih_gpp_si_age) - call this%set_history_var(vname='FATES_GPP_LU', units='kg m-2 s-1', & - long='gross primary productivity by land use type in kg carbon per m2 per second', & - use_default='inactive', avgflag='A', vtype=site_landuse_r8, & - hlms='CLM:ALM', upfreq=group_hifr_complx, ivar=ivar, initialize=initialize_variables, & - index = ih_gpp_si_landuse) - call this%set_history_var(vname='FATES_RDARK_USTORY_SZ', & units = 'kg m-2 s-1', & long='dark respiration for understory plants in kg carbon per m2 per second by size', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 555fa526c1..17557d6964 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -321,6 +321,12 @@ subroutine zero_bcs(fates,s) fates%bc_in(s)%hksat_sisl(:) = 0.0_r8 end if + ! Diagnostic quantities for outputting FATES patch-resolved + fates%bc_in(s)%lhflux_pa(:) = 0._r8 + fates%bc_in(s)%shflux_pa(:) = 0._r8 + fates%bc_in(s)%swabs_pa(:) = 0._r8 + fates%bc_in(s)%netlw_pa(:) = 0._r8 + fates%bc_in(s)%t2m_pa(:) = 0._r8 ! Output boundaries fates%bc_out(s)%active_suction_sl(:) = .false. @@ -552,6 +558,12 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats, end if ! Land use + ! Diagnostic quantities for outputting FATES patch-resolved + allocate(bc_in%lhflux_pa(0:maxpatch_total)) + allocate(bc_in%shflux_pa(0:maxpatch_total)) + allocate(bc_in%swabs_pa(0:maxpatch_total)) + allocate(bc_in%netlw_pa(0:maxpatch_total)) + allocate(bc_in%t2m_pa(0:maxpatch_total)) ! harvest flag denote data from hlm, ! while the logging flag signifies only that logging is occurring (which could just be FATES logging) diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 416af2728a..a49cf824bd 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -593,6 +593,13 @@ module FatesInterfaceTypesMod real(r8),allocatable :: h2o_liq_sisl(:) ! Liquid water mass in each layer (kg/m2) real(r8) :: smpmin_si ! restriction for min of soil potential (mm) + ! Diagnostic quantities for outputting FATES patch-resolved + real(r8),allocatable :: lhflux_pa(:) ! latent heat flux + real(r8),allocatable :: shflux_pa(:) ! sensible heat flux + real(r8),allocatable :: swabs_pa(:) ! shortwave absorbed radiation + real(r8),allocatable :: netlw_pa(:) ! longwave net radiation + real(r8),allocatable :: t2m_pa(:) ! 2m air temeprature + ! Land use ! --------------------------------------------------------------------------------- real(r8),allocatable :: hlm_harvest_rates(:) ! annual harvest rate per cat from hlm for a site