diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index e21e1c56e7..f8ad384e72 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -777,7 +777,7 @@ end subroutine canopy_spread ! ===================================================================================== - subroutine canopy_summarization( nsites, sites, bc_in ) + subroutine canopy_summarization( nsites, sites ) ! ---------------------------------------------------------------------------------- ! Much of this routine was once ed_clm_link minus all the IO and history stuff @@ -793,7 +793,6 @@ subroutine canopy_summarization( nsites, sites, bc_in ) ! !ARGUMENTS integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) - type(bc_in_type) , intent(in) :: bc_in(nsites) ! ! !LOCAL VARIABLES: type (fates_patch_type) , pointer :: currentPatch diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index f509aec000..78dbc94205 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -463,7 +463,7 @@ subroutine PreDisturbanceLitterFluxes( currentSite, currentPatch, bc_in ) !------------------------------------------------------------------------------------ ! Calculate the fragmentation rates - call fragmentation_scaler(currentPatch, bc_in) + call fragmentation_scaler(currentPatch, currentSite%bc_in(currentPatch%patchno)) do el = 1, num_elements diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 2f358cda63..e673110118 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -69,7 +69,7 @@ module FatesPatchMod type (fates_patch_type), pointer :: older => null() ! pointer to next older patch type (fates_patch_type), pointer :: younger => null() ! pointer to next younger patch type (fates_cohort_vec_type), pointer :: co_scr(:) ! Scratch vector of cohort properties - + !--------------------------------------------------------------------------- ! INDICES @@ -296,7 +296,7 @@ subroutine Init(this, num_swb, num_levsoil) allocate(this%sabs_dif(num_swb)) allocate(this%fragmentation_scaler(num_levsoil)) allocate(this%co_scr(max_cohort_per_patch)) - + ! initialize all values to nan call this%NanValues() @@ -534,7 +534,7 @@ subroutine NanValues(this) this%scorch_ht(:) = nan this%tfc_ros = nan this%frac_burnt = nan - + end subroutine NanValues !=========================================================================== @@ -732,11 +732,11 @@ subroutine Create(this, age, area, land_use_label, nocomp_pft, num_swb, num_pft, integer, intent(in) :: num_levsoil ! number of soil layers integer, intent(in) :: current_tod ! time of day [seconds past 0Z] integer, intent(in) :: regeneration_model ! regeneration model version - + ! initialize patch ! sets all values to nan, then some values to zero call this%Init(num_swb, num_levsoil) - + ! initialize running means for patch call this%InitRunningMeans(current_tod, regeneration_model, num_pft) @@ -1341,6 +1341,6 @@ subroutine CheckVars(this, var_aliases, return_code) end subroutine CheckVars - !=========================================================================== +! ====================================================================================== end module FatesPatchMod diff --git a/biogeochem/FatesSoilBGCFluxMod.F90 b/biogeochem/FatesSoilBGCFluxMod.F90 index e0af2a9e0f..4da6180516 100644 --- a/biogeochem/FatesSoilBGCFluxMod.F90 +++ b/biogeochem/FatesSoilBGCFluxMod.F90 @@ -586,7 +586,7 @@ end subroutine EffluxIntoLitterPools ! ===================================================================================== - subroutine FluxIntoLitterPools(csite, bc_in, bc_out) + subroutine FluxIntoLitterPools(csite) ! ----------------------------------------------------------------------------------- ! Created by Charlie Koven and Rosie Fisher, 2014-2015 @@ -619,15 +619,11 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) use FatesConstantsMod, only : itrue use FatesGlobals, only : endrun => fates_endrun use EDParamsMod , only : ED_val_cwd_flig, ED_val_cwd_fcel - - implicit none ! !ARGUMENTS - type(ed_site_type) , intent(inout) :: csite - type(bc_in_type) , intent(in) :: bc_in - type(bc_out_type) , intent(inout),target :: bc_out + type(ed_site_type), target, intent(inout) :: csite ! !LOCAL VARIABLES: type (fates_patch_type), pointer :: currentPatch @@ -635,9 +631,10 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) real(r8), pointer :: flux_cel_si(:) real(r8), pointer :: flux_lab_si(:) real(r8), pointer :: flux_lig_si(:) + real(r8), pointer :: flux_all_si type(litter_type), pointer :: litt - - real(r8) :: surface_prof(bc_in%nlevsoil) ! this array is used to distribute + + real(r8), allocatable :: surface_prof(:) ! this array is used to distribute ! fragmented litter on the surface ! into the soil/decomposition ! layers. It exponentially decays @@ -652,6 +649,7 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) integer :: ic ! CWD type index integer :: ipft ! PFT index + ! The following are used for the MIMICS ligC/N boundary condition real(r8) :: leaf_c, sapw_c ! leaf and sapwood carbon, per plant [kg] real(r8) :: fnrt_c, struct_c ! fineroot and struct carbon, per plant [kg] @@ -671,75 +669,90 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) ! how steep profile is for surface components (1/ e_folding depth) (1/m) real(r8), parameter :: surfprof_exp = 10. - ! This is the number of effective soil layers to transfer from - nlev_eff_soil = max(bc_in%max_rooting_depth_index_col, 1) + ! Loop over patches + currentPatch => csite%oldest_patch + flux_patch_loop: do while (associated(currentPatch)) + + if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then + + associate( & + bc_out => csite%bc_out(currentPatch%patchno), & + bc_in => csite%bc_in(currentPatch%patchno) & + ) + + ! This is the number of effective soil layers to transfer from + nlev_eff_soil = max(bc_in%max_rooting_depth_index_col, 1) - ! The decomposition layers are most likely the exact same layers - ! as the soil layers (same depths also), unless it is a simplified - ! single layer case, where nlevdecomp = 1 + ! The decomposition layers are most likely the exact same layers + ! as the soil layers (same depths also), unless it is a simplified + ! single layer case, where nlevdecomp = 1 - nlev_eff_decomp = min(bc_in%nlevdecomp,nlev_eff_soil) + nlev_eff_decomp = min(bc_in%nlevdecomp,nlev_eff_soil) - ! define a single shallow surface profile for surface additions - ! (leaves, stems, and N deposition). This sends the above ground - ! mass into the soil pools using an exponential depth decay function. - ! Since it is sending an absolute mass [kg] into variable layer - ! widths, we multiply the profile by the layer width, so that - ! wider layers get proportionally more. After the masses - ! are sent, each layer will normalize by depth. + ! define a single shallow surface profile for surface additions + ! (leaves, stems, and N deposition). This sends the above ground + ! mass into the soil pools using an exponential depth decay function. + ! Since it is sending an absolute mass [kg] into variable layer + ! widths, we multiply the profile by the layer width, so that + ! wider layers get proportionally more. After the masses + ! are sent, each layer will normalize by depth. - surface_prof(:) = 0._r8 - z_decomp = 0._r8 - do id = 1,nlev_eff_decomp - z_decomp = z_decomp+0.5*bc_in%dz_decomp_sisl(id) - surface_prof(id) = exp(-surfprof_exp * z_decomp) * bc_in%dz_decomp_sisl(id) - z_decomp = z_decomp+0.5*bc_in%dz_decomp_sisl(id) - end do - surface_prof_tot = sum(surface_prof) - do id = 1,nlev_eff_decomp - surface_prof(id) = surface_prof(id)/surface_prof_tot - end do - - ! Loop over the different elements. - do el = 1, num_elements - - ! Zero out the boundary flux arrays - ! Make a pointer to the cellulose, labile and lignin - ! flux partitions. - - select case (element_list(el)) - case (carbon12_element) - bc_out%litt_flux_cel_c_si(:) = 0.0_r8 - bc_out%litt_flux_lig_c_si(:) = 0.0_r8 - bc_out%litt_flux_lab_c_si(:) = 0.0_r8 - flux_cel_si => bc_out%litt_flux_cel_c_si(:) - flux_lab_si => bc_out%litt_flux_lab_c_si(:) - flux_lig_si => bc_out%litt_flux_lig_c_si(:) - case (nitrogen_element) - bc_out%litt_flux_cel_n_si(:) = 0._r8 - bc_out%litt_flux_lig_n_si(:) = 0._r8 - bc_out%litt_flux_lab_n_si(:) = 0._r8 - flux_cel_si => bc_out%litt_flux_cel_n_si(:) - flux_lab_si => bc_out%litt_flux_lab_n_si(:) - flux_lig_si => bc_out%litt_flux_lig_n_si(:) - case (phosphorus_element) - bc_out%litt_flux_cel_p_si(:) = 0._r8 - bc_out%litt_flux_lig_p_si(:) = 0._r8 - bc_out%litt_flux_lab_p_si(:) = 0._r8 - flux_cel_si => bc_out%litt_flux_cel_p_si(:) - flux_lab_si => bc_out%litt_flux_lab_p_si(:) - flux_lig_si => bc_out%litt_flux_lig_p_si(:) - end select - - currentPatch => csite%oldest_patch - do while (associated(currentPatch)) + allocate(surface_prof(bc_in%nlevsoil)) + surface_prof(:) = 0._r8 + z_decomp = 0._r8 + do id = 1,nlev_eff_decomp + z_decomp = z_decomp+0.5*bc_in%dz_decomp_sisl(id) + surface_prof(id) = exp(-surfprof_exp * z_decomp) * bc_in%dz_decomp_sisl(id) + z_decomp = z_decomp+0.5*bc_in%dz_decomp_sisl(id) + end do + surface_prof_tot = sum(surface_prof) + do id = 1,nlev_eff_decomp + surface_prof(id) = surface_prof(id)/surface_prof_tot + end do + + ! Loop over the different elements. + flux_elem_loop: do el = 1, num_elements + + ! Zero out the boundary flux arrays + ! Make a pointer to the cellulose, labile and lignin + ! flux partitions. + + select case (element_list(el)) + case (carbon12_element) + bc_out%litt_flux_cel_c_si(:) = 0.0_r8 + bc_out%litt_flux_lig_c_si(:) = 0.0_r8 + bc_out%litt_flux_lab_c_si(:) = 0.0_r8 + bc_out%litt_flux_all_c = 0.0_r8 + flux_cel_si => bc_out%litt_flux_cel_c_si(:) + flux_lab_si => bc_out%litt_flux_lab_c_si(:) + flux_lig_si => bc_out%litt_flux_lig_c_si(:) + flux_all_si => bc_out%litt_flux_all_c + case (nitrogen_element) + bc_out%litt_flux_cel_n_si(:) = 0._r8 + bc_out%litt_flux_lig_n_si(:) = 0._r8 + bc_out%litt_flux_lab_n_si(:) = 0._r8 + bc_out%litt_flux_all_n = 0.0_r8 + flux_cel_si => bc_out%litt_flux_cel_n_si(:) + flux_lab_si => bc_out%litt_flux_lab_n_si(:) + flux_lig_si => bc_out%litt_flux_lig_n_si(:) + flux_all_si => bc_out%litt_flux_all_n + case (phosphorus_element) + bc_out%litt_flux_cel_p_si(:) = 0._r8 + bc_out%litt_flux_lig_p_si(:) = 0._r8 + bc_out%litt_flux_lab_p_si(:) = 0._r8 + bc_out%litt_flux_all_p = 0.0_r8 + flux_cel_si => bc_out%litt_flux_cel_p_si(:) + flux_lab_si => bc_out%litt_flux_lab_p_si(:) + flux_lig_si => bc_out%litt_flux_lig_p_si(:) + flux_all_si => bc_out%litt_flux_all_p + end select ! Set a pointer to the litter object ! for the current element on the current ! patch litt => currentPatch%litter(el) area_frac = currentPatch%area/area - + do ic = 1, ncwd do id = 1,nlev_eff_decomp @@ -748,7 +761,6 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) flux_lig_si(id) = flux_lig_si(id) + & litt%ag_cwd_frag(ic) * ED_val_cwd_flig * area_frac * surface_prof(id) - end do do j = 1, nlev_eff_soil @@ -772,13 +784,13 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) do id = 1,nlev_eff_decomp flux_lab_si(id) = flux_lab_si(id) + & - litt%leaf_fines_frag(ilabile) * area_frac* surface_prof(id) + litt%leaf_fines_frag(ilabile) * area_frac * surface_prof(id) flux_cel_si(id) = flux_cel_si(id) + & - litt%leaf_fines_frag(icellulose) * area_frac* surface_prof(id) + litt%leaf_fines_frag(icellulose) * area_frac * surface_prof(id) flux_lig_si(id) = flux_lig_si(id) + & - litt%leaf_fines_frag(ilignin) * area_frac* surface_prof(id) + litt%leaf_fines_frag(ilignin) * area_frac * surface_prof(id) end do @@ -788,15 +800,14 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) do id = 1,nlev_eff_decomp flux_lab_si(id) = flux_lab_si(id) + & (litt%seed_decay(ipft) + litt%seed_germ_decay(ipft)) * & - EDPftvarcon_inst%lf_flab(ipft) * area_frac* surface_prof(id) + EDPftvarcon_inst%lf_flab(ipft) * area_frac * surface_prof(id) flux_cel_si(id) = flux_cel_si(id) + & (litt%seed_decay(ipft) + litt%seed_germ_decay(ipft)) * & - EDPftvarcon_inst%lf_fcel(ipft) * area_frac* surface_prof(id) - + EDPftvarcon_inst%lf_fcel(ipft) * area_frac * surface_prof(id) flux_lig_si(id) = flux_lig_si(id) + & (litt%seed_decay(ipft) + litt%seed_germ_decay(ipft)) * & - EDPftvarcon_inst%lf_flig(ipft) * area_frac* surface_prof(id) + EDPftvarcon_inst%lf_flig(ipft) * area_frac * surface_prof(id) end do end do @@ -810,119 +821,108 @@ subroutine FluxIntoLitterPools(csite, bc_in, bc_out) litt%root_fines_frag(ilignin,j) * area_frac enddo - currentPatch => currentPatch%younger - end do - - ! Normalize all masses over the decomposition layer's depth - ! Convert from kg/m2/day -> g/m3/s + ! Calculate the total flux of C into the litter pools + flux_all_si = sum(flux_cel_si(:)) + sum(flux_lig_si(:)) + sum(flux_lab_si(:)) - do id = 1,nlev_eff_decomp - flux_cel_si(id) = days_per_sec * g_per_kg * & - flux_cel_si(id) / bc_in%dz_decomp_sisl(id) - flux_lig_si(id) = days_per_sec * g_per_kg * & - flux_lig_si(id) / bc_in%dz_decomp_sisl(id) - flux_lab_si(id) = days_per_sec * g_per_kg * & - flux_lab_si(id) / bc_in%dz_decomp_sisl(id) - end do + end do flux_elem_loop + + ! Deallocate temporary array + deallocate(surface_prof) + ! Move to the next patch - end do ! do elements - - ! If we are coupled with MIMICS, then we need some assessment of litter quality - ! ie ligC/totalN. If we are not tracking N in the litter flux (ie C-only model) - ! then we need to approximate this by estimating the mean C:N ratios of each - ! plant organ, and mulitplying that by the different C Fluxes to get a total - ! approximate N flux. Note, in C-only, we will not capture any re-absorption. + ! If we are coupled with MIMICS, then we need some assessment of litter quality + ! ie ligC/totalN. If we are not tracking N in the litter flux (ie C-only model) + ! then we need to approximate this by estimating the mean C:N ratios of each + ! plant organ, and mulitplying that by the different C Fluxes to get a total + ! approximate N flux. Note, in C-only, we will not capture any re-absorption. - if(trim(hlm_decomp).eq.'MIMICS') then - - ! If we track nitrogen (ie cnp or other) then - ! we diagnose the c-lig/n ratio directly from the pools - if(element_pos(nitrogen_element)>0) then - - ! Sum totalN fluxes over depth [g/m2] - sum_N = sum((bc_out%litt_flux_cel_n_si(1:nlev_eff_soil) + & - bc_out%litt_flux_lig_n_si(1:nlev_eff_soil) + & - bc_out%litt_flux_lab_n_si(1:nlev_eff_soil)) * & - bc_in%dz_sisl(1:nlev_eff_soil)) - - else - - ! In this case (Carbon Only), we use the stoichiometry parameters to estimate - ! the C:N of live vegetation and the seedbank, and use that - ! as a proxy for the C:N of the litter flux - - sum_N = 0._r8 - - currentPatch => csite%oldest_patch - do while (associated(currentPatch)) - - litt => currentPatch%litter(element_pos(carbon12_element)) - area_frac = currentPatch%area*area_inv - - tot_leaf_c = 0._r8 - tot_leaf_n = 0._r8 - tot_fnrt_c = 0._r8 - tot_fnrt_n = 0._r8 - tot_wood_c = 0._r8 - tot_wood_n = 0._r8 - - ccohort => currentPatch%tallest - do while (associated(ccohort)) - ipft = ccohort%pft - leaf_c = ccohort%n * area_inv * ccohort%prt%GetState(leaf_organ, carbon12_element) - sapw_c = ccohort%n * area_inv * ccohort%prt%GetState(sapw_organ, carbon12_element) - fnrt_c = ccohort%n * area_inv * ccohort%prt%GetState(fnrt_organ, carbon12_element) - struct_c = ccohort%n * area_inv * ccohort%prt%GetState(struct_organ, carbon12_element) - leaf_n = leaf_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(leaf_organ)) - sapw_n = sapw_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(sapw_organ)) - fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(fnrt_organ)) - struct_n = struct_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(struct_organ)) - tot_leaf_c = tot_leaf_c + leaf_c - tot_leaf_n = tot_leaf_n + leaf_n - tot_fnrt_c = tot_fnrt_c + fnrt_c - tot_fnrt_n = tot_fnrt_n + fnrt_n - tot_wood_c = tot_wood_c + sapw_c + struct_c - tot_wood_n = tot_wood_n + sapw_n + struct_n - ccohort => ccohort%shorter - end do - - if(tot_wood_c>nearzero) then - sum_N = sum_N + area_frac*sum(litt%ag_cwd_frag)*(tot_wood_n/tot_wood_c) - sum_N = sum_N + area_frac*sum(litt%bg_cwd_frag)*(tot_wood_n/tot_wood_c) - end if - if(tot_leaf_c>nearzero)then - sum_N = sum_N + area_frac*sum(litt%leaf_fines_frag)*(tot_leaf_n / tot_leaf_c) - end if - if(tot_fnrt_c>nearzero)then - sum_N = sum_N + area_frac*sum(litt%root_fines_frag)*(tot_fnrt_n / tot_fnrt_c) - end if - do ipft = 1,numpft - sum_N = sum_N + area_frac * currentPatch%nitr_repro_stoich(ipft) * & - (litt%seed_decay(ipft) + litt%seed_germ_decay(ipft)) - end do - - currentPatch => currentPatch%younger - end do - - ! Convert from kg/m2/day -> g/m2/s - sum_N = sum_N * days_per_sec * g_per_kg - - end if - - ! Sum over layers and multiply by depth g/m3/s * m -> g/m2/s - sum_ligC = sum(bc_out%litt_flux_lig_c_si(1:nlev_eff_soil) * bc_in%dz_sisl(1:nlev_eff_soil)) - - if(sum_N>nearzero)then - bc_out%litt_flux_ligc_per_n = sum_ligC / sum_N - else - bc_out%litt_flux_ligc_per_n = 0._r8 - end if - - end if + if(trim(hlm_decomp).eq.'MIMICS') then + + ! If we track nitrogen (ie cnp or other) then + ! we diagnose the c-lig/n ratio directly from the pools + if(element_pos(nitrogen_element)>0) then + + ! Sum totalN fluxes over depth [g/m2] + sum_N = sum((bc_out%litt_flux_cel_n_si(1:nlev_eff_soil) + & + bc_out%litt_flux_lig_n_si(1:nlev_eff_soil) + & + bc_out%litt_flux_lab_n_si(1:nlev_eff_soil)) * & + bc_in%dz_sisl(1:nlev_eff_soil)) + + else + + ! In this case (Carbon Only), we use the stoichiometry parameters to estimate + ! the C:N of live vegetation and the seedbank, and use that + ! as a proxy for the C:N of the litter flux + + sum_N = 0._r8 + + litt => currentPatch%litter(element_pos(carbon12_element)) + + tot_leaf_c = 0._r8 + tot_leaf_n = 0._r8 + tot_fnrt_c = 0._r8 + tot_fnrt_n = 0._r8 + tot_wood_c = 0._r8 + tot_wood_n = 0._r8 + + ccohort => currentPatch%tallest + do while (associated(ccohort)) + ipft = ccohort%pft + leaf_c = ccohort%n * area_inv * ccohort%prt%GetState(leaf_organ, carbon12_element) + sapw_c = ccohort%n * area_inv * ccohort%prt%GetState(sapw_organ, carbon12_element) + fnrt_c = ccohort%n * area_inv * ccohort%prt%GetState(fnrt_organ, carbon12_element) + struct_c = ccohort%n * area_inv * ccohort%prt%GetState(struct_organ, carbon12_element) + leaf_n = leaf_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(leaf_organ)) + sapw_n = sapw_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(sapw_organ)) + fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(fnrt_organ)) + struct_n = struct_c * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(struct_organ)) + tot_leaf_c = tot_leaf_c + leaf_c + tot_leaf_n = tot_leaf_n + leaf_n + tot_fnrt_c = tot_fnrt_c + fnrt_c + tot_fnrt_n = tot_fnrt_n + fnrt_n + tot_wood_c = tot_wood_c + sapw_c + struct_c + tot_wood_n = tot_wood_n + sapw_n + struct_n + ccohort => ccohort%shorter + end do + + if(tot_wood_c>nearzero) then + sum_N = sum_N + area_frac*sum(litt%ag_cwd_frag)*(tot_wood_n/tot_wood_c) + sum_N = sum_N + area_frac*sum(litt%bg_cwd_frag)*(tot_wood_n/tot_wood_c) + end if + if(tot_leaf_c>nearzero)then + sum_N = sum_N + area_frac*sum(litt%leaf_fines_frag)*(tot_leaf_n / tot_leaf_c) + end if + if(tot_fnrt_c>nearzero)then + sum_N = sum_N + area_frac*sum(litt%root_fines_frag)*(tot_fnrt_n / tot_fnrt_c) + end if + do ipft = 1,numpft + sum_N = sum_N + area_frac*currentPatch%nitr_repro_stoich(ipft) * & + (litt%seed_decay(ipft) + litt%seed_germ_decay(ipft)) + end do + + ! Convert from kg/m2/day -> g/m2/s + sum_N = sum_N * days_per_sec * g_per_kg + + end if + + ! Sum over layers and multiply by depth g/m3/s * m -> g/m2/s + sum_ligC = sum(bc_out%litt_flux_lig_c_si(1:nlev_eff_soil) * bc_in%dz_sisl(1:nlev_eff_soil)) + + if(sum_N>nearzero)then + bc_out%litt_flux_ligc_per_n = sum_ligC / sum_N + else + bc_out%litt_flux_ligc_per_n = 0._r8 + end if + + end if + + end associate + end if + + + currentPatch => currentPatch%younger + end do flux_patch_loop - - - return end subroutine FluxIntoLitterPools diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index b4f93ba5c9..22370eaeaf 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -29,7 +29,7 @@ module EDAccumulateFluxesMod !------------------------------------------------------------------------------ - subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + subroutine AccumulateFluxes_ED(nsites, sites, bc_in, dt_time) ! ! !DESCRIPTION: @@ -40,14 +40,13 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) use EDTypesMod , only : ed_site_type, AREA use FatesPatchMod, only : fates_patch_type use FatesCohortMod, only : fates_cohort_type - use FatesInterfaceTypesMod , only : bc_in_type,bc_out_type + use FatesInterfaceTypesMod , only : bc_in_type ! ! !ARGUMENTS integer, intent(in) :: nsites type(ed_site_type), intent(inout), target :: sites(nsites) type(bc_in_type), intent(in) :: bc_in(nsites) - type(bc_out_type), intent(inout) :: bc_out(nsites) real(r8), intent(in) :: dt_time ! timestep interval ! ! !LOCAL VARIABLES: diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 09db71632c..be39037f00 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -53,7 +53,7 @@ module EDInitMod use EDTypesMod , only : phen_dstat_moistoff use EDTypesMod , only : phen_cstat_notcold use EDTypesMod , only : phen_dstat_moiston - use FatesInterfaceTypesMod , only : bc_in_type,bc_out_type + use FatesInterfaceTypesMod , only : bc_in_type use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesInterfaceTypesMod , only : hlm_use_inventory_init use FatesInterfaceTypesMod , only : hlm_use_fixed_biogeog @@ -128,7 +128,7 @@ module EDInitMod ! ============================================================================ - subroutine init_site_vars( site_in, bc_in, bc_out ) + subroutine init_site_vars( site_in, bc_in ) ! ! !DESCRIPTION: ! @@ -136,7 +136,6 @@ subroutine init_site_vars( site_in, bc_in, bc_out ) ! !ARGUMENTS type(ed_site_type), intent(inout) :: site_in type(bc_in_type),intent(in) :: bc_in - type(bc_out_type),intent(in) :: bc_out ! ! !LOCAL VARIABLES: !---------------------------------------------------------------------- diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 7a5bd21840..423fd3e586 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -779,6 +779,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) currentPatch => currentSite%youngest_patch do while(associated(currentPatch)) + + if (currentPatch%nocomp_pft_label .ne. nocomp_bareground) then call GenerateDamageAndLitterFluxes( currentSite, currentPatch) @@ -786,6 +788,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) call PreDisturbanceIntegrateLitter(currentPatch ) + end if + currentPatch => currentPatch%older enddo @@ -794,7 +798,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) ! can remove it completely if/when this call is added in ELM to ! subroutine UpdateLitterFluxes(this,bounds_clump) in elmfates_interfaceMod.F90 - call FluxIntoLitterPools(currentsite, bc_in, bc_out) + call FluxIntoLitterPools(currentsite) ! Update cohort number. diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 90be4df5ec..544defe5e0 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -9,6 +9,7 @@ module EDTypesMod use FatesConstantsMod, only : secondaryland use FatesConstantsMod, only : secondary_age_threshold use FatesConstantsMod, only : nearzero + use FatesConstantsMod , only : n_landuse_cats use FatesGlobals, only : fates_log use FatesHydraulicsMemMod, only : ed_cohort_hydr_type use FatesHydraulicsMemMod, only : ed_site_hydr_type @@ -27,15 +28,17 @@ module EDTypesMod use FatesConstantsMod, only : days_per_year use FatesRunningMeanMod, only : rmean_type,rmean_arr_type use FatesConstantsMod, only : fates_unset_r8 + use FatesConstantsMod, only : fates_unset_int use FatesInterfaceTypesMod,only : bc_in_type use FatesInterfaceTypesMod,only : bc_out_type - use FatesConstantsMod , only : n_landuse_cats use FatesInterfaceTypesMod,only : hlm_parteh_mode + use FatesInterfaceTypesMod,only : fates_interface_registry_type use FatesCohortMod, only : fates_cohort_type use FatesPatchMod, only : fates_patch_type use EDParamsMod, only : nclmax, nlevleaf, maxpft use FatesConstantsMod, only : n_dbh_bins, n_dist_types use shr_log_mod, only : errMsg => shr_log_errMsg + use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use SFFireWeatherMod, only : fire_weather implicit none @@ -327,10 +330,19 @@ module EDTypesMod type, public :: ed_site_type - ! POINTERS + !! POINTERS + + ! patch pointers type (fates_patch_type), pointer :: oldest_patch => null() ! pointer to oldest patch at the site type (fates_patch_type), pointer :: youngest_patch => null() ! pointer to yngest patch at the site + ! Boundary conditions + type(bc_in_type), allocatable :: bc_in(:) + type(bc_out_type), allocatable :: bc_out(:) + + ! Registry index array + integer, allocatable, private :: ridx(:) + ! Resource management type (ed_resources_management_type) :: resources_management ! resources_management at the site @@ -603,8 +615,11 @@ module EDTypesMod contains + procedure, public :: AllocateRegistryIndexArray procedure, public :: get_current_landuse_statevector procedure, public :: get_secondary_young_fraction + procedure, public :: GetRegistryIndex + procedure, public :: SetRegistryIndex end type ed_site_type @@ -615,6 +630,54 @@ module EDTypesMod contains + ! ============================================================================ + + subroutine AllocateRegistryIndexArray(this, patches_per_site) + + ! Allocates an array that maps the local patch number on a site to its + ! associated registry index + + ! Arguments + class(ed_site_type), intent(inout) :: this + integer, intent(in) :: patches_per_site + + ! Allocate the registry index array + allocate(this%ridx(patches_per_site)) + this%ridx = fates_unset_int + + end subroutine AllocateRegistryIndexArray + + ! ============================================================================ + + integer function GetRegistryIndex(this, ifp) result(ridx) + + ! Gets the registry index for a given local patch number from the + ! site object registry index mapping array + + ! Arguments + class(ed_site_type), intent(in) :: this + integer, intent(in) :: ifp + + ridx = this%ridx(ifp) + + end function GetRegistryIndex + + ! ============================================================================ + + subroutine SetRegistryIndex(this, ifp, ridx) + + ! Sets the registry index for a given local patch number in the + ! site object registry index mapping array + + ! Arguments + class(ed_site_type), intent(inout) :: this + integer, intent(in) :: ifp + integer, intent(in) :: ridx + + this%ridx(ifp) = ridx + + end subroutine SetRegistryIndex + ! ============================================================================ subroutine set_patchno( currentSite, check , call_id) diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 4d535f9de4..743a6c22df 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -108,7 +108,7 @@ module FatesConstantsMod integer, parameter, public :: ican_ustory = 2 ! nominal index for diagnostics that refer to understory layers ! (all layers that are not the top canopy layer) - ! Flags specifying how phosphorous uptake and turnover interacts + ! Flags specifying how phosphorus uptake and turnover interacts ! with the host model. integer, public, parameter :: prescribed_p_uptake = 1 integer, public, parameter :: coupled_p_uptake = 2 diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 6982dd700b..6cdde56391 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -46,7 +46,8 @@ module FatesInterfaceMod use FatesGlobals , only : fates_global_verbose use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun - use FatesConstantsMod , only : fates_unset_r8 + use FatesConstantsMod , only : fates_unset_r8 + use FatesConstantsMod , only : fates_unset_int use FatesLitterMod , only : ncwd use FatesLitterMod , only : ndcmpy use EDPftvarcon , only : FatesReportPFTParams @@ -114,6 +115,7 @@ module FatesInterfaceMod use FatesTwoStreamUtilsMod, only : TransferRadParams use LeafBiophysicsMod , only : lb_params use LeafBiophysicsMod , only : FvCB1980 + ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -134,6 +136,7 @@ module FatesInterfaceMod ! grid-cell, this is intended to be migrated to columns integer :: nsites + integer :: npatches type(ed_site_type), pointer :: sites(:) @@ -160,7 +163,28 @@ module FatesInterfaceMod type(bc_pconst_type) :: bc_pconst + ! This is the interface registry which associates variables with a common keyword + type(fates_interface_registry_type), allocatable :: registry(:) + + ! Index filter array for registries that are active (i.e. has a FATES patch) + integer, allocatable :: filter_registry_active(:) + + ! Active vegetated patches + integer :: num_active_patches + + + + contains + procedure :: CheckInterfaceVariables + procedure, public :: InitializeInterfaceRegistry + procedure, public :: InitializeFatesSites + procedure, public :: InitializeBoundaryConditions + procedure :: SetRegistryActiveState + procedure :: SetRegistryLastState + procedure, public :: UpdateInterfaceVariables + procedure, public :: UpdateLitterFluxes + end type fates_interface_type character(len=*), parameter :: sourcefile = & @@ -186,7 +210,7 @@ module FatesInterfaceMod public :: DetermineGridCellNeighbors logical :: debug = .false. ! for debugging this module - + contains ! ==================================================================================== @@ -302,8 +326,6 @@ subroutine zero_bcs(fates,s) fates%bc_in(s)%tot_litc = 0.0_r8 fates%bc_in(s)%snow_depth_si = 0.0_r8 fates%bc_in(s)%frac_sno_eff_si = 0.0_r8 - fates%bc_in(s)%w_scalar_sisl(:) = 0.0_r8 - fates%bc_in(s)%t_scalar_sisl(:) = 0.0_r8 if(do_fates_salinity)then fates%bc_in(s)%salinity_sl(:) = 0.0_r8 @@ -500,8 +522,6 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats, allocate(bc_in%z_sisl(nlevsoil_in)) allocate(bc_in%decomp_id(nlevsoil_in)) allocate(bc_in%dz_decomp_sisl(nlevdecomp_in)) - allocate(bc_in%w_scalar_sisl(nlevsoil_in)) - allocate(bc_in%t_scalar_sisl(nlevsoil_in)) ! Lightning (or successful ignitions) and population density ! Fire related variables @@ -2234,6 +2254,115 @@ end subroutine set_fates_ctrlparms ! ==================================================================================== + subroutine SetRegistryActiveState(this) + + ! This subroutine determines if a registry is active, i.e. whether it is associated + ! with a FATES patch. Registries are allocated during initialization to the maximum + ! number of possible patches, but not all FATES patches will necessarily be created + ! particularly during a cold start. As such, we need a procedure to determine with + ! each API update, which registries are active to avoid passing invalid data. + + ! Argument + class(fates_interface_type), intent(inout) :: this + + ! Locals + type(fates_patch_type), pointer :: currentPatch + + integer :: s ! site index + integer :: r ! registry index + integer :: i ! filter index + + ! Set all registries to inactive by default + do r = 1, this%npatches + call this%registry(r)%SetActiveState(active_state=.false.) + end do + + ! Set the active registry index filter to unset + this%filter_registry_active = fates_unset_int + + ! Set the active registry counter and filter index iterator to zero + this%num_active_patches = 0 + i = 0 + + ! Loop over sites and patches to set active registries and iterate the counter + do s = 1, this%nsites + currentPatch => this%sites(s)%oldest_patch + do while (associated(currentPatch)) + + if (currentPatch%nocomp_pft_label .ne. nocomp_bareground) then + + ! Get the registry index for the current site + patch combo and set it to active + r = this%sites(s)%GetRegistryIndex(currentPatch%patchno) + call this%registry(r)%SetActiveState(active_state=.true.) + + ! Increment the active patch counter and update the active index filter + this%num_active_patches = this%num_active_patches + 1 + i = i + 1 + this%filter_registry_active(i) = r + + end if + + ! Move to the next patch + currentPatch => currentPatch%younger + end do + end do + + end subroutine SetRegistryActiveState + + ! ==================================================================================== + + subroutine SetRegistryLastState(this) + + ! This procedure determines when a registry is the last in a series to be associated + ! with a specific subgrid object index (e.g. column index). The API design currently + ! implements any unit conversion or scaling as a last step after all data associated + ! with the specific subgrid object has been passed to the HLM. As such, we need to + ! identify when the next registry index in the series is associated with a different + ! subgrid object index. + + use FatesInterfaceParametersMod, only : registry_var_intid_column + + ! Argument + class(fates_interface_type), intent(inout) :: this + + ! Locals + integer :: n ! loop index + integer :: r, r_next ! registry index + integer :: c, c_next ! column index + + ! Loop over all active registries and set the last state to the current state + do n = 1, this%num_active_patches + + ! Set the registry indices for the current and previous active registry + r = this%filter_registry_active(n) + + ! Set last in subgrid state to false by default + call this%registry(r)%SetLastState() + + ! Get the column index for the current and previous active registry + c = this%registry(r)%GetColumnIndex() + + ! If the next index is greater than the current, set the previous last state to true + if (n .lt. this%num_active_patches) then + r_next = this%filter_registry_active(n+1) + c_next = this%registry(r_next)%GetColumnIndex() + if (c_next .ne. c) call this%registry(r)%SetLastState(registry_var_intid_column) + ! if (c_next .ne. c) write(fates_log(),*) 'SRLS: n, r, c: ', n, r, c + else + call this%registry(r)%SetLastState(registry_var_intid_column) + ! write(fates_log(),*) 'SRLS: n, r, c: ', n, r, c + end if + + ! ! If this is the last active registry, set its last state to true for all subgrid cases + ! if (n .eq. this%num_active_patches) call this%registry(r)%SetLastState(registry_var_intid_column) + ! if (n .eq. this%num_active_patches) write(fates_log(),*) 'SRLS: n, rp, cp: ', n, r_prev, c_prev + + end do + + end subroutine SetRegistryLastState + + ! ==================================================================================== + subroutine FatesReportParameters(masterproc) ! ----------------------------------------------------- @@ -2702,7 +2831,372 @@ subroutine FatesReadParameters(param_reader) call fates_params%Destroy() deallocate(fates_params) + +end subroutine FatesReadParameters + +! ====================================================================================== + +subroutine InitializeInterfaceRegistry(this, num_veg_patches, patchlist) + + ! This procedure intializes an interface registry for each patch index on the clump - end subroutine FatesReadParameters + ! Arguments + class(fates_interface_type), intent(inout) :: this ! fates interface + integer, intent(in) :: num_veg_patches ! number of veg patches in this clump + integer, intent(in) :: patchlist(:) ! list of hlm patches for registry index + + ! Locals + integer :: r ! registry index + + ! Allocate interface registries for each vegetated patch on the clump + allocate(this%registry(num_veg_patches)) + + ! Allocate the active registry filter array to the maximum number of possible active patches + allocate(this%filter_registry_active(num_veg_patches)) + + ! Set the number of vegetated patches to the interface type level + this%npatches = num_veg_patches + + ! Initialize and set patch index for each registry which is associated with a vegetated patch + do r = 1, num_veg_patches + + ! Initialize each registry with a dictionary of keys to register fates and hlm variables against + ! The keys are defined in the registry type-bound procedures + call this%registry(r)%InitializeInterfaceRegistry() + + ! Set the HLM patch index with the current registry + call this%registry(r)%SetSubgridIndices(hlmpatch=patchlist(r)) + + end do + +end subroutine InitializeInterfaceRegistry + +! ====================================================================================== + +subroutine InitializeFatesSites(this, patches_per_site) + + ! This procedure initializes the FATES sites by iterating through the number of + ! vegetated patches and incrementing the number of sites when a new gridcell index + ! is encountered. + + class(fates_interface_type), intent(inout) :: this ! fates interface + integer, intent(in) :: patches_per_site ! number of patches per site + + ! Local + integer :: r ! interface registry index + integer :: g ! gridcell index + integer :: gc ! current gridcell index + integer :: s ! site counter/index + integer :: ifp ! fates patch counter/index + + ! Initialize the current gridcell index and the fates site counter + gc = fates_unset_int + s = 0 + ifp = 0 + + ! Iterate over the number of vegetated patches and determine + do r = 1, this%npatches + + ! Get the gridcell index + g = this%registry(r)%GetGridcellIndex() + + ! Update the fates counter + ifp = ifp + 1 + + ! Iterate the fates site and reset the fates patch counter for each new gridcell index + if (gc /= g) then + gc = g + s = s + 1 + ifp = 1 + end if + + ! Set the site and fates patch index for the current registry + call this%registry(r)%SetSubgridIndices(fatespatch=ifp, site=s) + + end do + + ! Set the number of fates sites for the interface + this%nsites = s + + ! Allocate the sites + allocate(this%sites(this%nsites)) + + ! Allocate the registry index array for the the sites + do s = 1, this%nsites + call this%sites(s)%AllocateRegistryIndexArray(patches_per_site) + end do + + ! Iterate through the registries again and store the registry indices for each site + do r = 1, this%npatches + + ! Get the site index for the current registry + s = this%registry(r)%GetSiteIndex() + ifp = this%registry(r)%GetFatesPatchIndex() + + ! Store the registry index for the current site and fates patch index + call this%sites(s)%SetRegistryIndex(ifp, r) + + end do + + +end subroutine InitializeFatesSites + +! ====================================================================================== + +subroutine InitializeBoundaryConditions(this, patches_per_site) + + ! This procedure initializes the boundary conditions arrays with the maximum number + ! of possible patches for each FATES site. It then registers the boundary condition + ! variables for each patch. Conversion factors are applied as necessary during + ! registration, which will be utilized during the update API calls. + + use FatesInterfaceParametersMod + use FatesConstantsMod , only : days_per_sec + + ! Arguments + class(fates_interface_type), intent(inout) :: this ! fates interface type + integer, intent(in) :: patches_per_site ! number of patches per site + + ! Locals + integer :: r ! registery iterator + integer :: s ! site iterator + integer :: ifp ! boundary condition index + integer :: nlevdecomp + type(bc_in_type), pointer :: bc_in + type(bc_out_type), pointer :: bc_out + + ! Register the input boundary conditions use for BC allocations + do s = 1, this%nsites + + ! Allocate boundary conditions for all sites with the maximum number of patches per site + allocate(this%sites(s)%bc_in(patches_per_site)) + allocate(this%sites(s)%bc_out(patches_per_site)) + + ! Iterate over the maximum number of patches for the current site + do ifp = 1, patches_per_site + + ! Create convenience pointers to the current boundary conditions + bc_in => this%sites(s)%bc_in(ifp) + bc_out => this%sites(s)%bc_out(ifp) + + ! Get the site associated with this registry + r = this%sites(s)%GetRegistryIndex(ifp) + + ! Register the boundary conditions that are necessary for allocating other boundary conditions first + call this%registry(r)%Register(key=hlm_fates_decomp_max, data=bc_in%nlevdecomp_full, hlm_flag=.false.) + call this%registry(r)%Register(key=hlm_fates_decomp, data=bc_in%nlevdecomp, hlm_flag=.false.) + call this%registry(r)%Register(key=hlm_fates_soil_level, data=bc_in%nlevsoil, hlm_flag=.false.) + + ! Initialize the interface variables necessary for allocating boundary conditions dimensions + call this%registry(r)%InitializeInterfaceVariablesDimensions() + + ! Initialize the currently registered boundary conditions + call bc_in%Initialize() + call bc_out%Initialize(bc_in) + + ! Register the remaining boundary conditions + ! bc_in + call this%registry(r)%Register(key=hlm_fates_thaw_max_depth_index, & + data=bc_in%max_thaw_depth_index, hlm_flag=.false.) + call this%registry(r)%Register(key=hlm_fates_decomp_thickness, & + data=bc_in%dz_decomp_sisl, hlm_flag=.false.) + + call this%registry(r)%Register(key=hlm_fates_decomp_frac_moisture, & + data=bc_in%w_scalar_sisl, hlm_flag=.false.) + call this%registry(r)%Register(key=hlm_fates_decomp_frac_temperature, & + data=bc_in%t_scalar_sisl, hlm_flag=.false.) + + ! bc_out + nlevdecomp = bc_in%nlevdecomp + call this%registry(r)%Register(key=hlm_fates_litter_carbon_cellulose, & + data=bc_out%litt_flux_cel_c_si(1:nlevdecomp), hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + + call this%registry(r)%Register(key=hlm_fates_litter_carbon_lignin, & + data=bc_out%litt_flux_lig_c_si(1:nlevdecomp), hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + + call this%registry(r)%Register(key=hlm_fates_litter_carbon_labile, & + data=bc_out%litt_flux_lab_c_si(1:nlevdecomp), hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + + call this%registry(r)%Register(key=hlm_fates_litter_carbon_total, & + data=bc_out%litt_flux_all_c, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + + if (hlm_parteh_mode == prt_cnp_flex_allom_hyp) then + call this%registry(r)%Register(key=hlm_fates_litter_phosphorus_cellulose, & + data=bc_out%litt_flux_cel_p_si, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + call this%registry(r)%Register(key=hlm_fates_litter_phosphorus_lignin, & + data=bc_out%litt_flux_lig_p_si, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + call this%registry(r)%Register(key=hlm_fates_litter_phosphorus_labile, & + data=bc_out%litt_flux_lab_p_si, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + call this%registry(r)%Register(key=hlm_fates_litter_phosphorus_total, & + data=bc_out%litt_flux_all_p, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + + call this%registry(r)%Register(key=hlm_fates_litter_nitrogen_cellulose, & + data=bc_out%litt_flux_cel_n_si, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + call this%registry(r)%Register(key=hlm_fates_litter_nitrogen_lignin, & + data=bc_out%litt_flux_lig_n_si, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + call this%registry(r)%Register(key=hlm_fates_litter_nitrogen_labile, & + data=bc_out%litt_flux_lab_n_si, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + call this%registry(r)%Register(key=hlm_fates_litter_nitrogen_total, & + data=bc_out%litt_flux_all_n, hlm_flag=.false., & + conversion_factor=days_per_sec*g_per_kg) + end if + end do + end do + +end subroutine InitializeBoundaryConditions + +! ====================================================================================== + +subroutine UpdateInterfaceVariables(this, initialize, restarting) + + ! This procedure is responsible for updating the bulk of the interface variables by + ! iterating through each registry entry. It also handles special updates during + ! initialization and restart relating to the maximum rooting depth index, nlevdecomp + ! and the decomp_id. + + ! Arguments + class(fates_interface_type), intent(inout) :: this + logical, intent(in), optional :: initialize + logical, intent(in), optional :: restarting + + ! Locals + type(bc_in_type), pointer :: bc_in + type(bc_out_type), pointer :: bc_out + + logical :: initialize_local + logical :: restarting_local + + integer :: r ! registry interface index + integer :: s ! site index + integer :: ifp ! fates patch index + integer :: i ! layer index + + ! Set the default initialize flag to false + initialize_local = .false. + if (present(initialize)) then + initialize_local = initialize + end if + + ! Set the default restart flag to false + ! If we are restarting we need to initialize as well + restarting_local = .false. + if (present(restarting)) then + restarting_local = restarting + initialize_local = .true. + end if + + do r = 1, this%npatches + + ! Update the interface variables for the current registry + call this%registry(r)%Update() + + ! Get the site associated with this registry + s = this%registry(r)%GetSiteIndex() + ifp = this%registry(r)%GetFatesPatchIndex() + + bc_in => this%sites(s)%bc_in(ifp) + + ! Set the max rooting depth index to either to the soil depth or the max thaw depth last year, whichever is shallower + bc_in%max_rooting_depth_index_col = min(bc_in%nlevsoil, bc_in%max_thaw_depth_index) + + ! Calculate various bc_in variables that are based on other variables or namelist states + if (initialize_local) then + + ! Check vertical soil carbon decomposition usage + if (hlm_use_vertsoilc == itrue) then + if(bc_in%nlevdecomp .ne. bc_in%nlevsoil) then + write(fates_log(), *) 'The host has signaled a vertically resolved' + write(fates_log(), *) 'soil decomposition model. Therefore, the ' + write(fates_log(), *) 'total number of soil layers should equal the' + write(fates_log(), *) 'total number of decomposition layers.' + write(fates_log(), *) 'nlevdecomp: ',bc_in%nlevdecomp + write(fates_log(), *) 'nlevsoil: ',bc_in%nlevsoil + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Set all decomposition layer ids to their respective soil layer index + do i = 1, bc_in%nlevsoil + bc_in%decomp_id(i) = i + end do + + else ! No vertical soil carbon decomposition usage + if(bc_in%nlevdecomp .ne. 1)then + write(fates_log(), *) 'The host has signaled a non-vertically resolved' + write(fates_log(), *) 'soil decomposition model. Therefore, the ' + write(fates_log(), *) 'total number of decomposition layers should be 1.' + write(fates_log(), *) 'nlevdecomp: ',bc_in%nlevdecomp + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Set all decomposition layer ids to 1 + bc_in%decomp_id(:) = 1 + + end if + + ! On initialization, set the max rooting depth index to the maximum decomposition level + ! unless we are restarting + if (.not. restarting_local) bc_in%max_rooting_depth_index_col = bc_in%nlevdecomp + + end if + end do + +end subroutine UpdateInterfaceVariables + +! ====================================================================================== + +subroutine UpdateLitterFluxes(this, dtime) + + ! This procedure handles the updating of litter fluxes from FATES to the HLM. + + class(fates_interface_type), intent(inout) :: this + real(r8), intent(in) :: dtime ! HLM timestep + + ! Locals + integer :: n ! active registry index iterator + integer :: r ! registry index + integer :: next_r + + ! Set the registry active state + call this%SetRegistryActiveState() + + ! Set the registry last state + call this%SetRegistryLastState() + + ! Loop through the active registries and update the litter fluxes + do n = 1, this%num_active_patches + r = this%filter_registry_active(n) + + call this%registry(r)%UpdateLitterFluxes() + + end do + +end subroutine UpdateLitterFluxes + +! ====================================================================================== + +subroutine CheckInterfaceVariables(this) + class(fates_interface_type), intent(inout) :: this + + ! Locals + integer :: r + + do r = 1, this%npatches + call this%registry(r)%CheckInterfaceVariables() + end do + +end subroutine CheckInterfaceVariables + +! ====================================================================================== end module FatesInterfaceMod diff --git a/main/FatesInterfaceParametersMod.F90 b/main/FatesInterfaceParametersMod.F90 new file mode 100644 index 0000000000..31e4d8c0ae --- /dev/null +++ b/main/FatesInterfaceParametersMod.F90 @@ -0,0 +1,43 @@ +module FatesInterfaceParametersMod + + implicit none + private + + ! Registry key parameters + character(len=*), parameter, public :: hlm_fates_decomp= 'decomp_layers' + character(len=*), parameter, public :: hlm_fates_decomp_max = 'max_decomp_layers' + character(len=*), parameter, public :: hlm_fates_decomp_thickness= 'decomp_thickness' + character(len=*), parameter, public :: hlm_fates_thaw_max_depth_index = 'prioryear_thaw_max_depth_index' + character(len=*), parameter, public :: hlm_fates_soil_level = 'soil_level_number' + character(len=*), parameter, public :: hlm_fates_decomp_frac_moisture = 'decomp_frac_moisture' + character(len=*), parameter, public :: hlm_fates_decomp_frac_temperature = 'decomp_frac_temperature' + character(len=*), parameter, public :: hlm_fates_litter_carbon_cellulose = 'litter_carbon_cellulose' + character(len=*), parameter, public :: hlm_fates_litter_carbon_lignin = 'litter_carbon_lignin' + character(len=*), parameter, public :: hlm_fates_litter_carbon_labile = 'litter_carbon_labile' + character(len=*), parameter, public :: hlm_fates_litter_carbon_total= 'litter_carbon_total' + character(len=*), parameter, public :: hlm_fates_litter_phosphorus_cellulose = 'litter_phosphorus_cellulose' + character(len=*), parameter, public :: hlm_fates_litter_phosphorus_lignin = 'litter_phosphorus_lignin' + character(len=*), parameter, public :: hlm_fates_litter_phosphorus_labile = 'litter_phosphorus_labile' + character(len=*), parameter, public :: hlm_fates_litter_phosphorus_total= 'litter_phosphorus_total' + character(len=*), parameter, public :: hlm_fates_litter_nitrogen_cellulose = 'litter_nitrogen_cellulose' + character(len=*), parameter, public :: hlm_fates_litter_nitrogen_lignin = 'litter_nitrogen_lignin' + character(len=*), parameter, public :: hlm_fates_litter_nitrogen_labile = 'litter_nitrogen_labile' + character(len=*), parameter, public :: hlm_fates_litter_nitrogen_total= 'litter_nitrogen_total' + + ! Registry update frequency parameters + integer, parameter, public :: registry_update_init_dims = 0 ! variable dimension that needs to be updated during initialization + integer, parameter, public :: registry_update_init = 1 ! variable only needs to be updated during initialization + integer, parameter, public :: registry_update_daily = 2 ! variable needs to be updated daily + integer, parameter, public :: registry_update_timestep = 3 ! variable needs to be updated at each timestep + + ! Registry boundary condition parameters + integer, parameter, public :: registry_bc_in = 0 + integer, parameter, public :: registry_bc_out = 1 + + ! Registry subgrid object parameters + integer, parameter, public :: registry_var_intid_gridcell = 0 + integer, parameter, public :: registry_var_intid_topounit = 1 + integer, parameter, public :: registry_var_intid_landunit = 2 + integer, parameter, public :: registry_var_intid_column = 3 + +end module FatesInterfaceParametersMod \ No newline at end of file diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 22f032728c..aaed3585fe 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -2,11 +2,15 @@ module FatesInterfaceTypesMod use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse + use FatesConstantsMod , only : fates_unset_int use FatesGlobals , only : fates_global_verbose use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + use PRTGenericMod , only : prt_cnp_flex_allom_hyp + use FatesInterfaceVariableTypeMod, only : fates_interface_variable_type + use FatesInterfaceParametersMod implicit none @@ -398,7 +402,8 @@ module FatesInterfaceTypesMod ! Soil layer structure integer :: nlevsoil ! the number of soil layers in this column - integer :: nlevdecomp ! the number of soil layers in the column + integer :: nlevdecomp ! the number of active soil layers in the column + integer :: nlevdecomp_full ! the maximum possible soil layers for any column ! that are biogeochemically active real(r8),allocatable :: zi_sisl(:) ! interface level below a "z" level (m) ! this contains a zero index for surface. @@ -524,6 +529,10 @@ module FatesInterfaceTypesMod ! due to permafrost or bedrock constraints integer :: max_rooting_depth_index_col + ! The prior year maximum thaw depth index + ! Used to determine max_rooting_depth_index_col + integer :: max_thaw_depth_index + ! BGC Accounting real(r8) :: tot_het_resp ! total heterotrophic respiration (gC/m2/s) @@ -605,6 +614,10 @@ module FatesInterfaceTypesMod real(r8),allocatable :: hlm_sp_tsai(:) ! Interpolated sailt total SAI (stem area index) input from HLM per patch/pft real(r8),allocatable :: hlm_sp_htop(:) ! Interpolated daily canopy vegetation height input from HLM per patch/pft + contains + + procedure, public :: Initialize => InitializeBCIn + end type bc_in_type @@ -684,8 +697,10 @@ module FatesInterfaceTypesMod real(r8), allocatable :: litt_flux_cel_p_si(:) ! cellulose phosphorus litter, fates->BGC g/m3/s real(r8), allocatable :: litt_flux_lig_p_si(:) ! lignin phosphorus litter, fates->BGC g/m3/s real(r8), allocatable :: litt_flux_lab_p_si(:) ! labile phosphorus litter, fates->BGC g/m3/s + real(r8) :: litt_flux_all_c ! total litterfall carbon + real(r8) :: litt_flux_all_n ! total litterfall nitrogen + real(r8) :: litt_flux_all_p ! total litterfall phosphorus - ! MIMICS Boundary Conditions ! ----------------------------------------------------------------------------------- real(r8) :: litt_flux_ligc_per_n ! lignin carbon per total nitrogen @@ -806,6 +821,10 @@ module FatesInterfaceTypesMod real(r8) :: litter_cwd_c_si ! Total litter plus CWD carbon [Site-Level, gC m-2] real(r8) :: seed_c_si ! Total seed carbon [Site-Level, gC m-2] + contains + + procedure, public :: Initialize => InitializeBCOut + end type bc_out_type @@ -843,13 +862,106 @@ module FatesInterfaceTypesMod end type bc_pconst_type + ! Base type to be extended for the API registry + type, public :: fates_interface_registry_type + + ! Is registry have a FATES patch that exists? + logical, private :: active + + ! Container array of interface variables indexed by key + type(fates_interface_variable_type), allocatable :: hlm_vars(:) + type(fates_interface_variable_type), allocatable :: fates_vars(:) + + ! Array of keys associated with the interface variables + character(len=48), allocatable, private :: key(:) + + ! Variable regsitry metadata + integer :: num_api_vars ! number of variables in the registry + integer :: num_api_vars_update_init_dims ! number of variables dimensions needed during initialization + integer :: num_api_vars_update_init ! number of variables that update only at initialization + integer :: num_api_vars_update_daily ! number of variables that update daily + integer :: num_api_vars_update_timestep ! number of variables that update on the model timestep + integer :: num_api_vars_bc_in ! number of variables that are bc_in associated + integer :: num_api_vars_bc_out ! number of variables that are bc_ associated + integer :: num_api_vars_litter_flux ! number of variables that deal with all litter fluxes + + ! Array of update frequency values for each regsitry index + integer, allocatable :: update_frequency(:) + + ! Array of boundary condition directions for each registry index + integer, allocatable :: bc_dir(:) + + ! Arrays that hold the registry indices of variables based on update frequency + integer, allocatable :: filter_init_dims(:) ! registry index of variables dimensions that update at initialization + integer, allocatable :: filter_init(:) ! registry index of variables that update only at initialization + integer, allocatable :: filter_daily(:) ! registry index of variables that update daily + integer, allocatable :: filter_timestep(:) ! registry index of variables that update at each timestep + + ! Filter arrays that hold the registry indices for litter fluxes + integer, allocatable :: filter_litter_flux(:) + + ! Index arrays that map to the boundary condition types + integer, allocatable :: filter_bc_in(:) + integer, allocatable :: filter_bc_out(:) + + ! Subgrid index data + integer, private :: gidx + integer, private :: tidx + integer, private :: lidx + integer, private :: cidx + integer, private :: sidx + integer, private :: hpidx + integer, private :: fpidx + logical, private :: bareground + + ! Subgrid last state data + integer, private :: is_last_in_gridcell + integer, private :: is_last_in_topounit + integer, private :: is_last_in_landunit + integer, private :: is_last_in_column + + contains + + procedure :: CheckInterfaceVariables + procedure :: GetActivateState + procedure :: GetGridcellIndex + procedure :: GetLandunitIndex + procedure :: GetColumnIndex + procedure :: GetSiteIndex + procedure :: GetHLMPatchIndex + procedure :: GetFatesPatchIndex + procedure :: InitializeInterfaceRegistry + procedure :: InitializeInterfaceVariablesDimensions + procedure :: InitializeInterfaceVariables + procedure :: IsBareground => HLMPatchIsBareground + procedure :: SetSubgridIndices + procedure :: SetActiveState + procedure :: SetLastState + procedure :: UpdateLitterFluxes + procedure :: Update => UpdateInterfaceVariables + + generic :: Register => RegisterInterfaceVariables_0d, & + RegisterInterfaceVariables_1d, & + RegisterInterfaceVariables_2d + procedure, private :: RegisterInterfaceVariables_0d + procedure, private :: RegisterInterfaceVariables_1d + procedure, private :: RegisterInterfaceVariables_2d + + procedure, private :: DefineInterfaceRegistry + procedure, private :: DefineInterfaceVariable + procedure, private :: SetFilterMapArrays + procedure, private :: GetRegistryVariableIndex + procedure, private :: GetRegistryVariableKey + + end type fates_interface_registry_type + public :: ZeroBCOutCarbonFluxes contains - - ! ====================================================================================== + + ! ====================================================================================== - subroutine ZeroBCOutCarbonFluxes(bc_out) + subroutine ZeroBCOutCarbonFluxes(bc_out) ! !ARGUMENTS type(bc_out_type), intent(inout) :: bc_out @@ -861,4 +973,1046 @@ subroutine ZeroBCOutCarbonFluxes(bc_out) end subroutine ZeroBCOutCarbonFluxes + ! ====================================================================================== + + subroutine InitializeBCIn(this) + + ! This procedure initializes the input boundary condition structure by allocating + ! the associated variables arrays and unsetting their values. Note that some of the + ! input boundary conditions are used to initialize the output boundary conditions + ! and as such, this initialization must occur prior to the output boundary condition + ! initialization. + + ! Arguments + class(bc_in_type), intent(inout) :: this + + ! Allocate boundary condition arrays + allocate(this%decomp_id(this%nlevsoil)) + allocate(this%dz_decomp_sisl(this%nlevdecomp_full)) + allocate(this%w_scalar_sisl(this%nlevdecomp_full)) + allocate(this%t_scalar_sisl(this%nlevdecomp_full)) + + ! Unset variables + this%decomp_id = fates_unset_int + this%dz_decomp_sisl = nan + this%w_scalar_sisl = nan + this%t_scalar_sisl = nan + this%max_thaw_depth_index = fates_unset_int + + end subroutine InitializeBCIn + + ! ====================================================================================== + + subroutine InitializeBCOut(this, bc_in) + + ! This procedure initializes the output boundary condition structure by allocating + ! the associated variables arrays and unsetting their values. Note that the input + ! boundary conditions are needed to determine the set the size of some of the output + ! boundary condition arrays and as such, this procedure should be called after + ! the input boundary condition initialization. + + ! Arguments + class(bc_out_type), intent(inout) :: this + type(bc_in_type), intent(in) :: bc_in + + ! Allocate boundary condition arrays + allocate(this%litt_flux_cel_c_si(bc_in%nlevdecomp_full)) + allocate(this%litt_flux_lig_c_si(bc_in%nlevdecomp_full)) + allocate(this%litt_flux_lab_c_si(bc_in%nlevdecomp_full)) + + ! Unset the arrays + this%litt_flux_cel_c_si = nan + this%litt_flux_lig_c_si = nan + this%litt_flux_lab_c_si = nan + this%litt_flux_all_c = nan + + if (hlm_parteh_mode == prt_cnp_flex_allom_hyp) then + allocate(this%litt_flux_cel_n_si(bc_in%nlevdecomp_full)) + allocate(this%litt_flux_lig_n_si(bc_in%nlevdecomp_full)) + allocate(this%litt_flux_lab_n_si(bc_in%nlevdecomp_full)) + allocate(this%litt_flux_cel_p_si(bc_in%nlevdecomp_full)) + allocate(this%litt_flux_lig_p_si(bc_in%nlevdecomp_full)) + allocate(this%litt_flux_lab_p_si(bc_in%nlevdecomp_full)) + + this%litt_flux_cel_n_si = nan + this%litt_flux_lig_n_si = nan + this%litt_flux_lab_n_si = nan + this%litt_flux_all_n = nan + this%litt_flux_cel_p_si = nan + this%litt_flux_lig_p_si = nan + this%litt_flux_lab_p_si = nan + this%litt_flux_all_p = nan + end if + + end subroutine InitializeBCOut + + ! ====================================================================================== + + subroutine InitializeInterfaceRegistry(this) + + ! This initializes the interface registry with the associated variable counters, filters, + ! and metadata. This procedure calls the interface registry definition procedure which + ! is necessary to allow for the registration of the HLM and FATES variables. + + class(fates_interface_registry_type), intent(inout) :: this + + logical :: initialize + + ! Initialize registry counters + this%num_api_vars = 0 + this%num_api_vars_update_init_dims = 0 + this%num_api_vars_update_init = 0 + this%num_api_vars_update_daily = 0 + this%num_api_vars_update_timestep = 0 + this%num_api_vars_bc_in = 0 + this%num_api_vars_bc_out = 0 + this%num_api_vars_litter_flux = 0 + + ! First count up the keys defined in the registry and the registry counters + call this%DefineInterfaceRegistry(initialize=.false.) + + ! Allocate the registry variables arrays + allocate(this%fates_vars(this%num_api_vars)) + allocate(this%hlm_vars(this%num_api_vars)) + allocate(this%key(this%num_api_vars)) + allocate(this%update_frequency(this%num_api_vars)) + allocate(this%bc_dir(this%num_api_vars)) + + ! Allocate the index filter maps + allocate(this%filter_init_dims(this%num_api_vars_update_init_dims)) + allocate(this%filter_init(this%num_api_vars_update_init)) + allocate(this%filter_daily(this%num_api_vars_update_daily)) + allocate(this%filter_timestep(this%num_api_vars_update_timestep)) + + ! Allocate the boundary condition filter maps + allocate(this%filter_bc_in(this%num_api_vars_bc_in)) + allocate(this%filter_bc_out(this%num_api_vars_bc_out)) + + ! Allocate the litter flux filter + allocate(this%filter_litter_flux(this%num_api_vars_litter_flux)) + + ! Unset the allocatables not including the interface variables + this%update_frequency(:) = fates_unset_int + this%bc_dir(:) = fates_unset_int + this%filter_init_dims(:) = fates_unset_int + this%filter_init(:) = fates_unset_int + this%filter_daily(:) = fates_unset_int + this%filter_timestep(:) = fates_unset_int + this%filter_litter_flux(:) = fates_unset_int + + ! Now initialize the registry keys + call this%DefineInterfaceRegistry(initialize=.true.) + + ! Set filter map arrays + call this%SetFilterMapArrays() + + end subroutine InitializeInterfaceRegistry + + ! ====================================================================================== + + subroutine CheckInterfaceVariables(this) + + ! This procedure checks the registered HLM and FATES interface variables for consistency + + class(fates_interface_registry_type), intent(inout) :: this + + integer :: i + + do i = 1, this%num_api_vars + + + ! Check that variable keys match + if (this%hlm_vars(i)%key /= this%fates_vars(i)%key) then + write(*,*) "Key mismatch for variable: ", this%key(i), this%hlm_vars(i)%key, this%fates_vars(i)%key + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + ! Check that the rank matches + if (this%hlm_vars(i)%data_rank /= this%fates_vars(i)%data_rank) then + write(*,*) "Rank mismatch for variable: ", this%key(i) + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + ! Check that the bounds match + call this%hlm_vars(i)%CheckBounds(this%fates_vars(i)) + + ! Check that the size of the interface variables match + if (this%hlm_vars(i)%data_rank > 0) then + if (any(this%hlm_vars(i)%data_size(:) /= this%fates_vars(i)%data_size(:))) then + write(*,*) "Size mismatch: key, hlm size: ", this%key(i), this%hlm_vars(i)%data_size + write(*,*) "Size mismatch: key, fates size: ", this%key(i), this%fates_vars(i)%data_size + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + end if + + end do + + end subroutine CheckInterfaceVariables + + ! ====================================================================================== + + subroutine DefineInterfaceRegistry(this, initialize) + + ! This procedure defines the list of common names (i.e. "keys") to be associated with + ! FATES and HLM variables. All new keys, regardless of whether or not they are used + ! in a particular run mode or by a specific host land model, must be defined here. + ! For a future update, this procedure could be split in a manner similar to the fates + ! history interface module, where related keys are grouped together in separate procedure + ! calls to provide clarity to the developer and/or definion based on run mode . + + class(fates_interface_registry_type), intent(inout) :: this + + logical, intent(in) :: initialize ! false = count up the keys in the registry + + integer :: index ! Index to be incremented for each call to DefineInterfaceVariable() + + ! Initialize the index + index = 0 + + associate(bc_in => registry_bc_in, & + bc_out => registry_bc_out) + + ! Define the interface registry names and indices + ! Variables that need to be updated during initialization and are necessary for other boundary conditions + ! such as dimensions + call this%DefineInterfaceVariable(key=hlm_fates_decomp_max, initialize=initialize, index=index, & + update_frequency=registry_update_init_dims) + call this%DefineInterfaceVariable(key=hlm_fates_decomp, initialize=initialize, index=index, & + update_frequency=registry_update_init_dims) + call this%DefineInterfaceVariable(key=hlm_fates_soil_level, initialize=initialize, index=index, & + update_frequency=registry_update_init_dims) + + ! Variables that only need to be updated at initialization + call this%DefineInterfaceVariable(key=hlm_fates_decomp_thickness, initialize=initialize, index=index, & + update_frequency=registry_update_init) + + ! Variables that need to be updated daily + call this%DefineInterfaceVariable(key=hlm_fates_decomp_frac_moisture, initialize=initialize, index=index) + call this%DefineInterfaceVariable(key=hlm_fates_decomp_frac_temperature, initialize=initialize, index=index) + call this%DefineInterfaceVariable(key=hlm_fates_thaw_max_depth_index, initialize=initialize, index=index) + + ! Variables that need to be updated with each timestep + call this%DefineInterfaceVariable(key=hlm_fates_litter_carbon_cellulose, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_carbon_lignin, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_carbon_labile, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_carbon_total, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + + ! Define the N and P litter fluxes if in CNP mode + ! We could define the interface variables always, even if not registered, but this helps reduce the memory needs + if (hlm_parteh_mode == prt_cnp_flex_allom_hyp) then + call this%DefineInterfaceVariable(key=hlm_fates_litter_phosphorus_cellulose, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_phosphorus_lignin, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_phosphorus_labile, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_phosphorus_total, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + + call this%DefineInterfaceVariable(key=hlm_fates_litter_nitrogen_cellulose, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_nitrogen_lignin, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_nitrogen_labile, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + call this%DefineInterfaceVariable(key=hlm_fates_litter_nitrogen_total, initialize=initialize, index=index, & + update_frequency=registry_update_timestep, bc_dir=bc_out) + end if + + end associate + + end subroutine DefineInterfaceRegistry + + ! ====================================================================================== + + subroutine DefineInterfaceVariable(this, key, initialize, index, update_frequency, bc_dir) + + ! This procedure initializes the interface variables for this registry passing the + ! key associated with the variables along with other relavant data. The procedure + ! also increments various registry counters depending on the input arguments. + + class(fates_interface_registry_type), intent(inout) :: this + + character(len=*), intent(in) :: key + logical, intent(in) :: initialize + integer, intent(inout) :: index + integer, intent(in), optional :: update_frequency + integer, intent(in), optional :: bc_dir ! 0 = bc_in, 1 = bc_out, + + ! Local variables + integer :: index_type + integer :: update_frequency_local + integer :: bc_dir_local + + ! Increment the index + index = index + 1 + + ! If not initializing, increment the registry count variables, otherwise initialize the variable at the correct index + if (initialize) then + + ! Set the key for each index + this%key(index) = key + + ! Set the update frequency, default to daily update frequency + update_frequency_local = registry_update_daily + if (present(update_frequency)) then + select case (update_frequency) + case (registry_update_init_dims) + update_frequency_local = registry_update_init_dims + case (registry_update_init) + update_frequency_local = registry_update_init + case (registry_update_daily) + update_frequency_local = registry_update_daily + case (registry_update_timestep) + update_frequency_local = registry_update_timestep + case default + write(fates_log(),*) 'ERROR: Unrecognized update frequency in DefineInterfaceVariable(): ', update_frequency + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + else + end if + this%update_frequency(index) = update_frequency_local + + ! Set the boundary condition directionality for the variable index defaulting to bc_in + bc_dir_local = registry_bc_in + if (present(bc_dir)) then + bc_dir_local = bc_dir + end if + this%bc_dir(index) = bc_dir_local + + ! Initialize the interface variables and pass the key and update frequency to each for metadata + call this%hlm_vars(index)%Initialize(key, update_frequency_local, bc_dir_local) + call this%fates_vars(index)%Initialize(key, update_frequency_local, bc_dir_local) + + ! Not initializing, just counting the variables + else + + ! Increment the total API count + this%num_api_vars = this%num_api_vars + 1 + + ! Increment the count for the update frequency counts, defaulting to daily if not specified + if (present(update_frequency)) then + select case (update_frequency) + case (registry_update_init_dims) + this%num_api_vars_update_init_dims = this%num_api_vars_update_init_dims + 1 + case (registry_update_init) + this%num_api_vars_update_init = this%num_api_vars_update_init + 1 + case (registry_update_daily) + this%num_api_vars_update_daily = this%num_api_vars_update_daily + 1 + case (registry_update_timestep) + this%num_api_vars_update_timestep = this%num_api_vars_update_timestep + 1 + case default + write(fates_log(),*) 'ERROR: Unrecognized update frequency in DefineInterfaceVariable(): ', update_frequency + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + else + ! Default to daily update frequency + this%num_api_vars_update_daily = this%num_api_vars_update_daily + 1 + end if + + ! Increment the count for the boundary condition counters + if (present(bc_dir)) then + select case (bc_dir) + case (registry_bc_in) + this%num_api_vars_bc_in = this%num_api_vars_bc_in + 1 + case (registry_bc_out) + this%num_api_vars_bc_out = this%num_api_vars_bc_out + 1 + case default + write(fates_log(), *) 'ERROR: Unrecognized bc_dir in DefineInterfaceVariable(): ', bc_dir + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + else + ! defaults to bc_in + this%num_api_vars_bc_in = this%num_api_vars_bc_in + 1 + end if + + + ! Update the litter flux counters not including the total flux counter + if (key == hlm_fates_litter_carbon_cellulose .or. & + key == hlm_fates_litter_carbon_labile .or. & + key == hlm_fates_litter_carbon_lignin .or. & + key == hlm_fates_litter_nitrogen_cellulose .or. & + key == hlm_fates_litter_nitrogen_labile .or. & + key == hlm_fates_litter_nitrogen_lignin .or. & + key == hlm_fates_litter_phosphorus_cellulose .or. & + key == hlm_fates_litter_phosphorus_labile .or. & + key == hlm_fates_litter_phosphorus_lignin) then + this%num_api_vars_litter_flux= this%num_api_vars_litter_flux + 1 + end if + + end if + + end subroutine DefineInterfaceVariable + + ! ====================================================================================== + + subroutine SetSubgridIndices(this, gridcell, topounit, landunit, column, hlmpatch, fatespatch, site, bareground) + + ! This procedure sets the associated HLM subgrid index values for the calling registry. + ! This is provided so that the FATES registry and interface variable subroutines can + ! have access to the HLM subgrid indices as needed internal to the FATES code. + + class(fates_interface_registry_type), intent(inout) :: this + integer, intent(in), optional :: gridcell + integer, intent(in), optional :: topounit + integer, intent(in), optional :: landunit + integer, intent(in), optional :: column + integer, intent(in), optional :: hlmpatch + integer, intent(in), optional :: fatespatch + integer, intent(in), optional :: site + logical, intent(in), optional :: bareground + + if (present(gridcell)) this%gidx = gridcell + if (present(topounit)) this%tidx = topounit + if (present(landunit)) this%lidx = landunit + if (present(column)) this%cidx = column + if (present(hlmpatch)) this%hpidx = hlmpatch + if (present(fatespatch)) this%fpidx = fatespatch + if (present(site)) this%sidx = site + if (present(bareground)) this%bareground = bareground + + end subroutine SetSubgridIndices + + ! ====================================================================================== + + subroutine SetActiveState(this, active_state) + + ! This procedure sets the registry active state flag for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + logical, intent(in) :: active_state + + this%active = active_state + + end subroutine SetActiveState + + ! ====================================================================================== + + subroutine SetLastState(this, subgrid_type) + + ! Set the registry last state flag based on the subgrid type being processed + ! If a subgrid_type is not provided, then reset all last state flags to false + + class(fates_interface_registry_type), intent(inout) :: this + integer, intent(in), optional :: subgrid_type + + ! Local variables + integer :: n + logical :: last_state_local + integer :: subgrid_type_local + + ! If a subgrid type is provided, we set the last state flag for that type at the register level + ! and update the associated local subgrid type variable. + if (present(subgrid_type)) then + select case (subgrid_type) + case(registry_var_intid_gridcell) + this%is_last_in_gridcell = .true. + subgrid_type_local = registry_var_intid_gridcell + case(registry_var_intid_topounit) + this%is_last_in_topounit = .true. + subgrid_type_local = registry_var_intid_topounit + case(registry_var_intid_landunit) + this%is_last_in_landunit = .true. + subgrid_type_local = registry_var_intid_landunit + case(registry_var_intid_column) + this%is_last_in_column = .true. + subgrid_type_local = registry_var_intid_column + case default + write(fates_log(),*) 'ERROR: unrecognised subgrid_type provided to SetLastState()' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + last_state_local = .true. + else + this%is_last_in_gridcell = .false. + this%is_last_in_topounit = .false. + this%is_last_in_landunit = .false. + this%is_last_in_column = .false. + last_state_local = .false. + subgrid_type_local = fates_unset_int + end if + + ! Set the last state flags in the interface variables associated with the subgrid type + ! or if no subgrid type is provided, set all last state flags to false + do n = 1, this%num_api_vars + if (this%hlm_vars(n)%IsSubgridType(subgrid_type_local) .or. subgrid_type_local == fates_unset_int) then + ! write(fates_log(),*) 'SLS: n, lsl, stl: ' , n, last_state_local, subgrid_type_local + call this%hlm_vars(n)%SetLastState(last_state_local) + end if + end do + + end subroutine SetLastState + + ! ====================================================================================== + + logical function GetActivateState(this) result(active_state) + + ! This procedure gets the registry active state flag for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + active_state = this%active + + end function GetActivateState + + ! ====================================================================================== + + integer function GetGridcellIndex(this) result(gidx) + + ! This procedure gets the HLM gridcell index for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + gidx = this%gidx + + end function GetGridcellIndex + + ! ====================================================================================== + + integer function GetLandunitIndex(this) result(lidx) + + ! This procedure gets the HLM landunit index for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + lidx = this%lidx + + end function GetLandunitIndex + + ! ====================================================================================== + + integer function GetColumnIndex(this) result(cidx) + + ! This procedure gets the HLM column index for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + cidx = this%cidx + + end function GetColumnIndex + + ! ====================================================================================== + + integer function GetHLMPatchIndex(this) result(hpidx) + + ! This procedure gets the HLM patch index for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + hpidx = this%hpidx + + end function GetHLMPatchIndex + + ! ====================================================================================== + + integer function GetSiteIndex(this) result(sidx) + + ! This procedure gets the FATES site index for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + sidx = this%sidx + + end function GetSiteIndex + + ! ====================================================================================== + + integer function GetFatesPatchIndex(this) result(fpidx) + + ! This procedure gets the FATES patch index for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + fpidx = this%fpidx + + end function GetFatesPatchIndex + + ! ====================================================================================== + + logical function HLMPatchIsBareGround(this) result(bareground) + + ! This procedure gets the HLM bareground flag for the calling registry. + + class(fates_interface_registry_type), intent(inout) :: this + + bareground = this%bareground + + end function HLMPatchIsBareGround + + ! ====================================================================================== + + subroutine SetFilterMapArrays(this) + + ! This procedure sets the filter mapping arrays for the interface registry. These + ! filters are provided as a mean to quickly access the registry indices for specific + ! variable types (e.g. those that update on the model timestep) + + class(fates_interface_registry_type), intent(inout) :: this + + integer :: index + integer :: count_init_dims + integer :: count_init + integer :: count_daily + integer :: count_timestep + integer :: count_bc_in + integer :: count_bc_out + integer :: count_litter_flux + + ! Initialize counters + count_init_dims = 0 + count_init = 0 + count_daily = 0 + count_timestep = 0 + count_bc_in = 0 + count_bc_out = 0 + count_litter_flux= 0 + + ! Iterate over all registered variables and populate the filter maps accordingly + do index = 1, this%num_api_vars + + ! Frequency update + if (this%update_frequency(index) == registry_update_init_dims) then + count_init_dims = count_init_dims + 1 + this%filter_init_dims(count_init_dims) = index + else if (this%update_frequency(index) == registry_update_init) then + count_init = count_init + 1 + this%filter_init(count_init) = index + else if (this%update_frequency(index) == registry_update_daily) then + count_daily = count_daily + 1 + this%filter_daily(count_daily) = index + else if (this%update_frequency(index) == registry_update_timestep) then + count_timestep = count_timestep + 1 + this%filter_timestep(count_timestep) = index + else + write(fates_log(),*) 'ERROR: Unrecognized update frequency in SetFilterMapArrays(): ', this%update_frequency(index) + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + ! Boundary condition filter update + if (this%bc_dir(index) == registry_bc_in) then + count_bc_in = count_bc_in + 1 + this%filter_bc_in(count_bc_in) = index + else if (this%bc_dir(index) == registry_bc_out) then + count_bc_out = count_bc_out + 1 + this%filter_bc_out(count_bc_out) = index + end if + + ! Litter flux update + if (this%key(index) == hlm_fates_litter_carbon_cellulose .or. & + this%key(index) == hlm_fates_litter_nitrogen_cellulose .or. & + this%key(index) == hlm_fates_litter_phosphorus_cellulose .or. & + this%key(index) == hlm_fates_litter_carbon_labile .or. & + this%key(index) == hlm_fates_litter_nitrogen_labile .or. & + this%key(index) == hlm_fates_litter_phosphorus_labile .or. & + this%key(index) == hlm_fates_litter_carbon_lignin .or. & + this%key(index) == hlm_fates_litter_nitrogen_lignin .or. & + this%key(index) == hlm_fates_litter_phosphorus_lignin) then + count_litter_flux = count_litter_flux + 1 + this%filter_litter_flux(count_litter_flux) = index + end if + + + end do + + ! Check that the counts match the expected sizes + if (count_init_dims /= this%num_api_vars_update_init_dims .or. & + count_init /= this%num_api_vars_update_init .or. & + count_daily /= this%num_api_vars_update_daily .or. & + count_timestep /= this%num_api_vars_update_timestep .or. & + count_litter_flux /= this%num_api_vars_litter_flux) then + + write(fates_log(),*) 'ERROR: Mismatch in initialization counts in SetFilterMapArrays(): ' + write(fates_log(),*) ' count_init = ', count_init, ' expected = ', this%num_api_vars_update_init + write(fates_log(),*) ' count_daily = ', count_daily, ' expected = ', this%num_api_vars_update_daily + write(fates_log(),*) ' count_timestep = ', count_timestep, ' expected = ', this%num_api_vars_update_timestep + call endrun(msg=errMsg(__FILE__, __LINE__)) + + end if + + end subroutine SetFilterMapArrays + + ! ====================================================================================== + + subroutine RegisterInterfaceVariables_0d(this, key, data, hlm_flag, accumulate, is_first, subgrid_type, conversion_factor) + + ! This procedure associates a scalar data variable with a particular registry key + + ! Arguments + class(fates_interface_registry_type), intent(inout) :: this + class(*), target, intent(in) :: data ! data to be associated with key + character(len=*), intent(in) :: key ! variable registry key + logical, intent(in) :: hlm_flag ! Is the variable being register from the HLM? + logical, intent(in), optional :: accumulate ! Should the variable accumulate during the update? + logical, intent(in), optional :: is_first ! Should the variable be zeroed first? + integer, intent(in), optional :: subgrid_type ! The subgrid integer id associated with this HLM variable + real(r8), intent(in), optional :: conversion_factor ! Conversion factor for the variable + + ! Local + logical :: accumulate_local + logical :: is_first_local + integer :: subgrid_type_local + real(r8) :: conversion_factor_local + + ! Default accumulate to false + if (present(accumulate)) then + accumulate_local = accumulate + else + accumulate_local = .false. + end if + + ! Default is_first to false + if (present(is_first)) then + is_first_local = is_first + else + is_first_local = .false. + end if + + ! Default is_last to false + if (present(subgrid_type)) then + subgrid_type_local = subgrid_type + else + subgrid_type_local = fates_unset_int + end if + + ! Default conversion factor to 1.0 + if (present(conversion_factor)) then + conversion_factor_local = conversion_factor + else + conversion_factor_local = 1.0_r8 + end if + + ! Get index from registry key and associate the given data pointer + if (hlm_flag) then + call this%hlm_vars(this%GetRegistryVariableIndex(key))%Register(data, active=.true., & + accumulate=accumulate_local, & + is_first=is_first_local, & + subgrid_type=subgrid_type_local, & + conversion_factor=conversion_factor_local) + else + call this%fates_vars(this%GetRegistryVariableIndex(key))%Register(data, active=.true., & + accumulate=accumulate_local, & + is_first=is_first_local, & + subgrid_type=subgrid_type_local, & + conversion_factor=conversion_factor_local) + end if + + + end subroutine RegisterInterfaceVariables_0d + + ! ====================================================================================== + + subroutine RegisterInterfaceVariables_1d(this, key, data, hlm_flag, accumulate, is_first, subgrid_type, conversion_factor) + + ! This procedure associates 1D array data variable with a particular registry key + + class(fates_interface_registry_type), intent(inout) :: this + class(*), target, intent(in) :: data(:) ! data to be associated with key + character(len=*), intent(in) :: key ! variable registry key + logical, intent(in) :: hlm_flag ! Is the variable being register from the HLM? + logical, intent(in), optional :: accumulate ! Should the variable accumulate during the update? + logical, intent(in), optional :: is_first ! Should the variable be zeroed first? + integer, intent(in), optional :: subgrid_type ! The subgrid integer id associated with this HLM variable + real(r8), intent(in), optional :: conversion_factor ! Conversion factor for the variable + + ! Local + logical :: accumulate_local + logical :: is_first_local + integer :: subgrid_type_local + real(r8) :: conversion_factor_local + + ! Default accumulate to false + if (present(accumulate)) then + accumulate_local = accumulate + else + accumulate_local = .false. + end if + + ! Default is_first to false + if (present(is_first)) then + is_first_local = is_first + else + is_first_local = .false. + end if + + ! Default subgrid_type to fates_unset_int + if (present(subgrid_type)) then + subgrid_type_local = subgrid_type + else + subgrid_type_local = fates_unset_int + end if + + ! Default conversion factor to 1.0 + if (present(conversion_factor)) then + conversion_factor_local = conversion_factor + else + conversion_factor_local = 1.0_r8 + end if + + ! Default conversion factor to 1.0 + if (present(conversion_factor)) then + conversion_factor_local = conversion_factor + else + conversion_factor_local = 1.0_r8 + end if + + ! Get index from registry key and associate the given data pointer + if (hlm_flag) then + call this%hlm_vars(this%GetRegistryVariableIndex(key))%Register(data(:), active=.true., & + accumulate=accumulate_local, & + is_first=is_first_local, & + subgrid_type=subgrid_type_local, & + conversion_factor=conversion_factor_local) + else + call this%fates_vars(this%GetRegistryVariableIndex(key))%Register(data(:), active=.true., & + accumulate=accumulate_local, & + is_first=is_first_local, & + subgrid_type=subgrid_type_local, & + conversion_factor=conversion_factor_local) + end if + + end subroutine RegisterInterfaceVariables_1d + + ! ====================================================================================== + + subroutine RegisterInterfaceVariables_2d(this, key, data, hlm_flag, accumulate, is_first, subgrid_type, conversion_factor) + + ! This procedure associates 2D array data variable with a particular registry key + + class(fates_interface_registry_type), intent(inout) :: this + class(*), target, intent(in) :: data(:,:) ! data to be associated with key + character(len=*), intent(in) :: key ! variable registry key + logical, intent(in) :: hlm_flag ! Is the variable being register from the HLM? + logical, intent(in), optional :: accumulate ! Should the variable accumulate during the update? + logical, intent(in), optional :: is_first ! Should the variable be zeroed first? + integer, intent(in), optional :: subgrid_type ! The subgrid integer id associated with this HLM variable + real(r8), intent(in), optional :: conversion_factor ! Conversion factor for the variable + + ! Local + logical :: accumulate_local + logical :: is_first_local + integer :: subgrid_type_local + real(r8) :: conversion_factor_local + + ! Default accumulate to false + if (present(accumulate)) then + accumulate_local = accumulate + else + accumulate_local = .false. + end if + + ! Default is_first to false + if (present(is_first)) then + is_first_local = is_first + else + is_first_local = .false. + end if + + ! Default subgrid_type to fates_unset_int + if (present(subgrid_type)) then + subgrid_type_local = subgrid_type + else + subgrid_type_local = fates_unset_int + end if + + ! Default conversion factor to 1.0 + if (present(conversion_factor)) then + conversion_factor_local = conversion_factor + else + conversion_factor_local = 1.0_r8 + end if + + ! Get index from registry key and associate the given data pointer + if (hlm_flag) then + call this%hlm_vars(this%GetRegistryVariableIndex(key))%Register(data(:,:), active=.true., & + accumulate=accumulate_local, & + is_first=is_first_local, & + subgrid_type=subgrid_type_local, & + conversion_factor=conversion_factor_local) + else + call this%fates_vars(this%GetRegistryVariableIndex(key))%Register(data(:,:), active=.true., & + accumulate=accumulate_local, & + is_first=is_first_local, & + subgrid_type=subgrid_type_local, & + conversion_factor=conversion_factor_local) + end if + + end subroutine RegisterInterfaceVariables_2d + + ! ====================================================================================== + + subroutine InitializeInterfaceVariablesDimensions(this) + + ! This procedure initializes the interface variables that are used as dimensions + ! for initializing other interface variables. + + ! Arguments + class(fates_interface_registry_type), intent(inout) :: this ! registry being initialized + + ! Locals + integer :: i ! initialization iterator + integer :: j ! variable index + + ! Update the boundary conditions necessary during initialization only + do i = 1, this%num_api_vars_update_init_dims + + ! Get the variable index from the init filter + j = this%filter_init_dims(i) + + ! Update the variables + call this%fates_vars(j)%Update(this%hlm_vars(j)) + + end do + + end subroutine InitializeInterfaceVariablesDimensions + + ! ====================================================================================== + + subroutine InitializeInterfaceVariables(this) + + ! This procedure updates the interface variables that are used during initialization. + + ! Arguments + class(fates_interface_registry_type), intent(inout) :: this ! registry being initialized + + ! Locals + integer :: i ! initialization iterator + integer :: j ! variable index + + ! Update the boundary conditions necessary during initialization only + do i = 1, this%num_api_vars_update_init + + ! Get the variable index from the init filter + j = this%filter_init(i) + + ! Update the variables + call this%fates_vars(j)%Update(this%hlm_vars(j)) + + end do + + end subroutine InitializeInterfaceVariables + + ! ====================================================================================== + + subroutine UpdateInterfaceVariables(this) + + ! This procedure updates BC input interface variables + + class(fates_interface_registry_type), intent(inout) :: this + + integer :: n + integer :: ibc + + ! Iterate over all registered variables + do n = 1, this%num_api_vars_bc_in + ibc = this%filter_bc_in(n) + call this%fates_vars(ibc)%Update(this%hlm_vars(ibc)) + end do + + end subroutine UpdateInterfaceVariables + + ! ====================================================================================== + + ! subroutine UpdateLitterFluxes(this, dtime) + subroutine UpdateLitterFluxes(this) + + ! This procedure updates the litter flux output boundary conditions + ! This procedure is separated from other update calls as it happens + ! on the model time-step and has specific update calls for the total + ! litterfall per nutrient type based on run mode + + use FatesConstantsMod, only : days_per_sec + use FatesConstantsMod, only : g_per_kg + + ! Arguments + class(fates_interface_registry_type), intent(inout) :: this + ! real(r8), intent(in) :: dtime + + ! Locals + integer :: i + integer :: j + integer :: k + logical :: normalization_flag + + ! Iterate over the litter flux filter to update the individual litter types + do i = 1, this%num_api_vars_litter_flux + j = this%filter_litter_flux(i) + + ! Update the hlm variables with the fates variables + call this%hlm_vars(j)%Update(this%fates_vars(j), normalization_flag) + + ! If the normalization flag is set, scale the updated HLM litter flux interface variable + ! by the decomposition layer thickness. + ! The normalization factor could be set during the registration call in a manner similar + ! to the unit conversation factor, but given that not all interface variables may need normalization + ! and the normalization factor is not necessarily a simple standard constant, we handle it here for now. + if (normalization_flag) then + + ! Get the index for the decomposition thickness key + k = this%GetRegistryVariableIndex(hlm_fates_decomp_thickness) + + ! Normalize the litter fluxes against the decomposition layer thicknesses + call this%hlm_vars(j)%Normalize(this%fates_vars(k)) + + end if + + end do + + ! Update the HLM variable with the total litterfall + ! These are not currently included in the litter flux filter as these are not normalized + j = this%GetRegistryVariableIndex(hlm_fates_litter_carbon_total) + call this%hlm_vars(j)%Update(this%fates_vars(j)) + + if (hlm_parteh_mode == prt_cnp_flex_allom_hyp) then + j = this%GetRegistryVariableIndex(hlm_fates_litter_phosphorus_total) + call this%hlm_vars(j)%Update(this%fates_vars(j)) + + j = this%GetRegistryVariableIndex(hlm_fates_litter_nitrogen_total) + call this%hlm_vars(j)%Update(this%fates_vars(j)) + end if + + + end subroutine UpdateLitterFluxes + + ! ====================================================================================== + + integer function GetRegistryVariableIndex(this, key) result(index) + + ! This procedure returns the index associated with the key provided + + class(fates_interface_registry_type), intent(in) :: this + + character(len=*), intent(in) :: key ! variable registry key to search + + integer :: ivar ! Iterator + + ! Iterate over the registry variables until the associated key is found + do ivar = 1, this%num_api_vars + if (this%key(ivar) == key) then + index = ivar + return + end if + end do + + end function GetRegistryVariableIndex + + ! ====================================================================================== + + function GetRegistryVariableKey(this, index) result(key) + + ! This procedure returns the index associated with the key provided + + class(fates_interface_registry_type), intent(in) :: this + + integer, intent(in) :: index ! variable registry index + character(len=:), allocatable :: key + + key = this%key(index) + + end function GetRegistryVariableKey + + ! ====================================================================================== + end module FatesInterfaceTypesMod diff --git a/main/FatesInterfaceVarTypeMod.F90 b/main/FatesInterfaceVarTypeMod.F90 new file mode 100644 index 0000000000..fbb965e033 --- /dev/null +++ b/main/FatesInterfaceVarTypeMod.F90 @@ -0,0 +1,605 @@ +module FatesInterfaceVariableTypeMod + + ! This module contains the type definition and associated type-bound procedures + ! used to create an indexed list of associated HLM and FATES variables that are + ! related across the HLM-FATES application programming interface (API). + ! This method is largely inspired by the FATES history infrastructure + + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) + + use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : fates_long_string_length + use FatesConstantsMod, only : fates_unset_int + + implicit none + private + + ! Interface registry variable type + ! This defined type holds the pointer to HLM or FATES variable that is + ! associated with a common "key" name (defined in the interface parameter types). + ! It also includes a set of data characterizing the type of data not defined by + ! type or rank (e.g. HLM subgrid association, update frequency, etc). + + type, public :: fates_interface_variable_type + + character(len=48) :: key ! common registry key + class(*), pointer :: data0d ! scalar polymorphic data pointer + class(*), pointer :: data1d(:) ! 1D polymorphic data pointer + class(*), pointer :: data2d(:,:) ! 2D polymorphic data pointer + class(*), pointer :: data3d(:,:,:) ! 3D polymorphic data pointer + logical :: active ! true if the variable is used by the host land model + logical :: accumulate ! If true, this variable should add the source to the target + logical :: zero_first ! If true, zero the target variable before accumulation + logical :: is_last ! True if the variable is associated with the last patch for the associated subgrid unit + integer :: bc_dir ! 0 if bc_in, 1 if bc_out + integer :: data_rank ! rank of the variable (0, 1, 2, or 3) + integer :: update_frequency ! frequency of updates + integer :: subgrid_type ! subgrid integer ID associated with this variable + real(r8) :: conversion_factor ! conversion factor to adjust units as necessary + integer, allocatable :: data_size(:) ! size of the first dimension of the variable + + contains + procedure :: CheckBounds + procedure :: Convert => ConvertScaleInterfaceVariable + procedure :: Initialize => InitializeInterfaceVariable + procedure :: IsSubgridType + procedure :: Normalize => NormalizeLitterVariable + procedure :: Update => UpdateInterfaceVariable + procedure :: Dump + procedure :: SetLastState + + generic :: Register => RegisterInterfaceVariable_0d, & + RegisterInterfaceVariable_1d, & + RegisterInterfaceVariable_2d + procedure, private :: RegisterInterfaceVariable_0d + procedure, private :: RegisterInterfaceVariable_1d + procedure, private :: RegisterInterfaceVariable_2d + + procedure, private :: CompareRegistryVariableSizes + + end type fates_interface_variable_type + + contains + + ! ==================================================================================== + + subroutine CheckBounds(this, var) + + ! This procedure checks that the data bounds of the calling variable are consistent + ! with the data bounds of the input argument variable + + class(fates_interface_variable_type), intent(in) :: this + class(fates_interface_variable_type), intent(in) :: var + + ! Locals + integer :: ub1, ub2 + integer :: lb1, lb2 + logical :: bounds_mismatch + + bounds_mismatch = .false. + + if (this%data_rank >= 1) then + ub1 = ubound(this%data1d, dim=1) + lb1 = lbound(this%data1d, dim=1) + ub2 = ubound(var%data1d, dim=1) + lb2 = lbound(var%data1d, dim=1) + if (ub1 /= ub2 .or. lb1 /= lb2) then + write(fates_log(),*) 'Dimension 1 bounds mismatch for variable: ', this%key + write(fates_log(),*) 'Upper bounds: ', ub1, ', ', ub2 + write(fates_log(),*) 'Lower bounds: ', lb1, ', ', lb2 + bounds_mismatch = .true. + end if + else if (this%data_rank >= 2) then + ub1 = ubound(this%data2d, dim=2) + lb1 = lbound(this%data2d, dim=2) + ub2 = ubound(var%data2d, dim=2) + lb2 = lbound(var%data2d, dim=2) + if (ub1 /= ub2 .or. lb1 /= lb2) then + write(fates_log(),*) 'Dimension 2 bounds mismatch for variable: ', this%key + write(fates_log(),*) 'Upper bounds: ', ub1, ', ', ub2 + write(fates_log(),*) 'Lower bounds: ', lb1, ', ', lb2 + bounds_mismatch = .true. + end if + end if + + if (bounds_mismatch) then + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + end subroutine CheckBounds + + ! ==================================================================================== + + subroutine ConvertScaleInterfaceVariable(this, scalar) + + ! This subroutine scales a given interface variable data. This + ! routine is provided primarily to allow for post accumulation conversion + ! of units and scaling of data as necessary. Note that this only provides + ! scaling for real data types currently. + + class(fates_interface_variable_type), intent(inout) :: this + real(r8), intent(in) :: scalar + + ! Locals + integer :: i + character(len=fates_long_string_length) :: err_msg = 'FATES Error: invalid interface type being scaled' + + select case (this%data_rank) + case(0) + select type(var => this%data0d) + type is (real(r8)) + var = var * scalar + class default + write(fates_log(),*), err_msg + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case(1) + select type(var => this%data1d) + type is (real(r8)) + var = var * scalar + class default + write(fates_log(),*), err_msg + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case(2) + select type(var => this%data2d) + type is (real(r8)) + var = var * scalar + class default + write(fates_log(),*), err_msg + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case default + write(fates_log(),*) 'FATES ERROR: Unsupported interface variable rank in ConvertScaleInterfaceVariable' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + + end subroutine ConvertScaleInterfaceVariable + + ! ==================================================================================== + + subroutine Dump(this) + + ! This procedure will print output to the log file. Developers should use this + ! procedure for diagnostic purposes when attempt to review interface variable + ! data as it includes select type statements to resolve the polymorphic data types. + + class(fates_interface_variable_type), intent(in) :: this + + write(fates_log(),*) 'FATES Interface Variable Dump:' + write(fates_log(),*) ' Key: ', this%key + write(fates_log(),*) ' Active: ', this%active + write(fates_log(),*) ' Accumulate: ', this%accumulate + write(fates_log(),*) ' Data Rank: ', this%data_rank + write(fates_log(),*) ' Data Size: ', this%data_size + + select case (this%data_rank) + case(0) + select type(var => this%data0d) + type is (real(r8)) + write(fates_log(),*) ' Data (real): ', var + type is (integer) + write(fates_log(),*) ' Data (integer): ', var + class default + write(fates_log(),*), 'FATES ERROR: Unsupported interface variable type' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case(1) + select type(var => this%data1d) + type is (real(r8)) + write(fates_log(),*) ' Data (real): ', var + type is (integer) + write(fates_log(),*) ' Data (integer): ', var + class default + write(fates_log(),*), 'FATES ERROR: Unsupported interface variable type' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case(2) + select type(var => this%data2d) + type is (real(r8)) + write(fates_log(),*) ' Data (real): ', var + type is (integer) + write(fates_log(),*) ' Data (integer): ', var + class default + write(fates_log(),*), 'FATES ERROR: Unsupported interface variable type' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case default + write(fates_log(),*) 'FATES ERROR: Unsupported interface variable rank in Dump' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + + end subroutine Dump + + ! ==================================================================================== + + subroutine InitializeInterfaceVariable(this, key, update_frequency, bc_dir) + + ! This procedure initializes an interface variable assigning its key + ! update frequency and boundary condition direction. It also unsets + ! all other variables including the data pointers. + + class(fates_interface_variable_type), intent(inout) :: this + character(len=*), intent(in) :: key + integer, intent(in) :: update_frequency + integer, intent(in) :: bc_dir + + + allocate(this%data_size(3)) + + ! Initialize components that are set later + this%data_size = fates_unset_int + this%data_rank = fates_unset_int + this%data0d => null() + this%data1d => null() + this%data2d => null() + this%data3d => null() + this%active = .false. + this%accumulate = .false. + this%zero_first = .false. + this%is_last = .false. + this%conversion_factor = nan + this%subgrid_type = fates_unset_int + + ! Initialize registry variable components that are updated at variable definition + this%key = key + this%update_frequency = update_frequency + this%bc_dir = bc_dir + + end subroutine InitializeInterfaceVariable + + ! ==================================================================================== + + logical function IsSubgridType(this, subgrid_type) + + ! This procedure will check if the subgrid type of the interface variable + ! matches the input argument subgrid type + + class(fates_interface_variable_type), intent(in) :: this + integer, intent(in) :: subgrid_type + + IsSubgridType = (this%subgrid_type == subgrid_type) + + end function IsSubgridType + + ! ==================================================================================== + + subroutine NormalizeLitterVariable(this, norm_var) + + ! This subroutine is specifically set up to normalize the litter flux interface + ! variables by the decomposition thickness variable passed in as norm_factor. + ! It could be extended later to handle other normalization needs as necessary. + + class(fates_interface_variable_type), intent(inout) :: this + class(fates_interface_variable_type), intent(in) :: norm_var ! normalization interface variable + + ! Locals + integer :: i + character(len=fates_long_string_length) :: err_msg = 'FATES Error: invalid interface type being normalized' + + select case (this%data_rank) + case(1) + select type(var => this%data1d) + type is (real(r8)) + select type(norm => norm_var%data1d) + type is (real(r8)) + do i = lbound(var, dim=1), ubound(var, dim=1) + var(i) = var(i) / norm(i) + end do + class default + write(fates_log(),*), err_msg + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + class default + write(fates_log(),*), err_msg + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case default + write(fates_log(),*) 'FATES ERROR: Unsupported interface variable rank in NormalizeLitterVariable' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + + end subroutine NormalizeLitterVariable + + ! ==================================================================================== + + subroutine RegisterInterfaceVariable_0d(this, data, active, accumulate, is_first, subgrid_type, conversion_factor) + + ! This procedure registers the interface variable by associating the data pointer with + ! the scalar input data argument. It also sets a number of required and optional arguments + ! for metadata associated with the variable. + + class(fates_interface_variable_type), intent(inout) :: this + class(*), target, intent(in) :: data + logical, intent(in) :: active + logical, intent(in) :: accumulate + logical, intent(in) :: is_first + integer, intent(in) :: subgrid_type + real(r8), intent(in) :: conversion_factor + + this%data0d => data + this%active = active + this%accumulate = accumulate + this%zero_first = is_first + this%data_rank = rank(data) + this%conversion_factor = conversion_factor + this%subgrid_type = subgrid_type + + end subroutine RegisterInterfaceVariable_0d + + ! ==================================================================================== + + subroutine RegisterInterfaceVariable_1d(this, data, active, accumulate, is_first, subgrid_type, conversion_factor) + + ! This procedure registers the interface variable by associating the data pointer with + ! the 1D array input data argument. It also sets a number of required and optional arguments + ! for metadata associated with the variable. + + class(fates_interface_variable_type), intent(inout) :: this + + class(*), target, intent(in) :: data(:) + logical, intent(in) :: active + logical, intent(in) :: accumulate + logical, intent(in) :: is_first + integer, intent(in) :: subgrid_type + real(r8), intent(in) :: conversion_factor + + this%data1d => data(:) + this%active = active + this%accumulate = accumulate + this%zero_first = is_first + this%data_rank = rank(data) + this%data_size(1) = size(data, dim=1) + this%conversion_factor = conversion_factor + this%subgrid_type = subgrid_type + + end subroutine RegisterInterfaceVariable_1d + + ! ==================================================================================== + + subroutine RegisterInterfaceVariable_2d(this, data, active, accumulate, is_first, subgrid_type, conversion_factor) + + ! This procedure registers the interface variable by associating the data pointer with + ! the 2D array input data argument. It also sets a number of required and optional arguments + ! for metadata associated with the variable. + + class(fates_interface_variable_type), intent(inout) :: this + class(*), target, intent(in) :: data(:,:) + logical, intent(in) :: active + logical, intent(in) :: accumulate + logical, intent(in) :: is_first + integer, intent(in) :: subgrid_type + real(r8), intent(in) :: conversion_factor + + this%data2d => data(:,:) + this%active = active + this%accumulate = accumulate + this%zero_first = is_first + this%data_rank = rank(data) + this%data_size(1) = size(data, dim=1) + this%data_size(2) = size(data, dim=2) + this%subgrid_type = subgrid_type + this%conversion_factor = conversion_factor + + end subroutine RegisterInterfaceVariable_2d + + ! ==================================================================================== + + subroutine SetLastState(this, last_state) + + ! This procedure sets the last state logical metadata for the calling interface variable + + class(fates_interface_variable_type), intent(inout) :: this + logical, intent(in) :: last_state + + this%is_last = last_state + + end subroutine SetLastState + + ! ==================================================================================== + + subroutine UpdateInterfaceVariable(this, var, is_last) + + ! This is the main procedure that drives the updates between the HLM and FATES. + ! It updates the calling interface variable data with data from the argument + ! interface variable. Directionality of the update is from the argument to the + ! caller. The bulk of the code here is type selects to handle the polymorphic + ! data pointers in the interface variables types. + ! This procedure utilizes metadata stored in the interface variables to determine + ! how updates should be handled, for example whether or not to zero the calling + ! interface variable data prior to adding the argument variable data to it. + + ! Arguments + class(fates_interface_variable_type), intent(inout) :: this ! variable being updated + class(fates_interface_variable_type), intent(in) :: var ! variable update source + logical, intent(out), optional :: is_last ! true if last patch for subgrid unit + + ! Locals + character(len=fates_long_string_length) :: msg_mismatch = 'FATES ERROR: Mismatched interface variable types' + + ! Check that the dimensions of the source and target match + call this%CompareRegistryVariableSizes(var) + + ! If the is_last argument is present, output the is_last flag state for this variable + if (present(is_last)) is_last = this%is_last + + ! Update the data of the target variable using the source variable data pointer + ! Make sure the types match for the polymorphic data to allow for copying from the + ! source to the target. + ! Note that due to the use of polymorphic pointers, we must use select type constructs + ! to determine the actual type of the data being pointed to allowing for type-safe assignment. + ! This currently only supports real and integer types and no conversion between types + ! should be performed + select case (this%data_rank) + case(0) + select type(dest => this%data0d) + type is (real(r8)) + select type(source => var%data0d) + type is (real(r8)) + if (this%accumulate) then + if (this%zero_first) then + dest = 0.0_r8 + end if + dest = dest + source + else + dest = source + end if + ! Apply conversion factor if this is the last patch associated with the subgrid unit + if (this%is_last) then + dest = dest * var%conversion_factor + end if + class default + write(fates_log(),*), msg_mismatch + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + type is (integer) + select type(source => var%data0d) + type is (integer) + if (this%accumulate) then + if (this%zero_first) then + dest = 0 + end if + dest = dest + source + else + dest = source + end if + class default + write(fates_log(),*), msg_mismatch + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + class default + write(fates_log(),*), 'FATES ERROR: Unsupported interface variable type' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + + case(1) + + select type(dest => this%data1d) + type is (real(r8)) + select type(source => var%data1d) + type is (real(r8)) + if (this%accumulate) then + if (this%zero_first) then + dest = 0.0_r8 + end if + dest = dest + source + else + dest = source + end if + if (this%is_last) then + dest = dest * var%conversion_factor + end if + class default + write(fates_log(),*), msg_mismatch + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + type is (integer) + select type(source => var%data1d) + type is (integer) + if (this%accumulate) then + if (this%zero_first) then + dest = 0 + end if + dest = dest + source + else + dest = source + end if + class default + write(fates_log(),*), msg_mismatch + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + class default + write(fates_log(),*), 'FATES ERROR: Unsupported interface variable type' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + + case(2) + + select type(dest => this%data2d) + type is (real(r8)) + select type(source => var%data2d) + type is (real(r8)) + if (this%accumulate) then + if (this%zero_first) then + dest = 0.0_r8 + end if + dest = dest + source + else + dest = source + end if + if (this%is_last) then + dest = dest * var%conversion_factor + end if + class default + write(fates_log(),*), msg_mismatch + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + type is (integer) + select type(source => var%data2d) + type is (integer) + if (this%accumulate) then + if (this%zero_first) then + dest = 0 + end if + dest = dest + source + else + dest = source + end if + class default + write(fates_log(),*), msg_mismatch + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + class default + write(fates_log(),*), 'FATES ERROR: Unsupported interface variable type' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + case default + write(fates_log(),*) 'FATES ERROR: Unsupported interface variable target rank in UpdateInterfaceVariable' + call endrun(msg=errMsg(__FILE__, __LINE__)) + end select + + end subroutine UpdateInterfaceVariable + + ! ==================================================================================== + + subroutine CompareRegistryVariableSizes(this, var) + + ! This procedure checks to make sure that the ranks and size of the + ! data associated with the interface variables is the same. + + class(fates_interface_variable_type), intent(in) :: this ! variable being updated + class(fates_interface_variable_type), intent(in) :: var ! variable update source + + if (this%data_size(1) /= var%data_size(1) .or. & + this%data_size(2) /= var%data_size(2) .or. & + this%data_size(3) /= var%data_size(3)) then + + write(fates_log(),*) 'FATES ERROR: Mismatched interface variable sizes in UpdateInterfaceVariable' + + if (this%data_rank == 1) then + write(fates_log(),*) ' Target, size: ', this%key, this%data_size(1) + write(fates_log(),*) ' Source, size: ', var%key, var%data_size(1) + else if (this%data_rank == 2) then + write(fates_log(),*) ' Target, size: ', this%key, this%data_size(1), this%data_size(2) + write(fates_log(),*) ' Source, size: ', var%key, var%data_size(1), var%data_size(2) + else if (this%data_rank == 3) then + write(fates_log(),*) ' Target, size: ', this%key, this%data_size(1), this%data_size(2), this%data_size(3) + write(fates_log(),*) ' Source, size: ', var%key, var%data_size(1), var%data_size(2), var%data_size(3) + else + write(fates_log(),*) ' Unsupported interface variable rank in UpdateErrorMessage' + write(fates_log(),*) ' Target key, rank: ', this%key, this%data_rank + write(fates_log(),*) ' Source key, rank: ', var%key, var%data_rank + end if + + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + end subroutine CompareRegistryVariableSizes + + ! ==================================================================================== + +end module FatesInterfaceVariableTypeMod \ No newline at end of file diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index d1d689fa62..c9e98c49ad 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -3085,7 +3085,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in, bc_out) io_idx_si = this%restart_map(nc)%site_index(s) io_idx_co_1st = this%restart_map(nc)%cohort1_index(s) - call init_site_vars( sites(s), bc_in(s), bc_out(s) ) + call init_site_vars( sites(s), bc_in(s) ) call zero_site( sites(s) ) if ( rio_npatch_si(io_idx_si)<0 .or. rio_npatch_si(io_idx_si) > 10000 ) then