diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 6847a37867..ddab8e74af 100755 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -7,16 +7,16 @@ module EDCanopyStructureMod use FatesConstantsMod , only : r8 => fates_r8 use FatesGlobals , only : fates_log - use EDPftvarcon , only : EDPftvarcon_inst + use EDPftvarcon , only : EDPftvarcon_inst use EDGrowthFunctionsMod , only : c_area use EDCohortDynamicsMod , only : copy_cohort, terminate_cohorts, fuse_cohorts use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, ncwd use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf - use EDTypesMod , only : numpft_ed use EDtypesMod , only : AREA use FatesGlobals , only : endrun => fates_endrun use FatesInterfaceMod , only : hlm_days_per_year + use FatesInterfaceMod , only : numpft ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -825,7 +825,6 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use EDPatchDynamicsMod , only : set_root_fraction use EDTypesMod , only : sizetype_class_index use EDGrowthFunctionsMod , only : tree_lai, c_area - use EDEcophysConType , only : EDecophyscon use EDtypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst @@ -941,7 +940,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) use EDGrowthFunctionsMod , only : tree_lai, tree_sai, c_area use EDtypesMod , only : area, dinc_ed, hitemax, n_hite_bins - use EDEcophysConType , only : EDecophyscon ! ! !ARGUMENTS @@ -1051,19 +1049,19 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort => currentPatch%shortest do while(associated(currentCohort)) ft = currentCohort%pft - min_chite = currentCohort%hite - currentCohort%hite * EDecophyscon%crown(ft) + min_chite = currentCohort%hite - currentCohort%hite * EDPftvarcon_inst%crown(ft) max_chite = currentCohort%hite do iv = 1,N_HITE_BINS frac_canopy(iv) = 0.0_r8 ! this layer is in the middle of the canopy if(max_chite > maxh(iv).and.min_chite < minh(iv))then - frac_canopy(iv)= min(1.0_r8,dh / (currentCohort%hite*EDecophyscon%crown(ft))) + frac_canopy(iv)= min(1.0_r8,dh / (currentCohort%hite*EDPftvarcon_inst%crown(ft))) ! this is the layer with the bottom of the canopy in it. elseif(min_chite < maxh(iv).and.min_chite > minh(iv).and.max_chite > maxh(iv))then - frac_canopy(iv) = (maxh(iv) -min_chite ) / (currentCohort%hite*EDecophyscon%crown(ft)) + frac_canopy(iv) = (maxh(iv) -min_chite ) / (currentCohort%hite*EDPftvarcon_inst%crown(ft)) ! this is the layer with the top of the canopy in it. elseif(max_chite > minh(iv).and.max_chite < maxh(iv).and.min_chite < minh(iv))then - frac_canopy(iv) = (max_chite - minh(iv)) / (currentCohort%hite*EDecophyscon%crown(ft)) + frac_canopy(iv) = (max_chite - minh(iv)) / (currentCohort%hite*EDPftvarcon_inst%crown(ft)) elseif(max_chite < maxh(iv).and.min_chite > minh(iv))then !the whole cohort is within this layer. frac_canopy(iv) = 1.0_r8 endif @@ -1111,7 +1109,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort => currentCohort%taller enddo !currentCohort lai = 0.0_r8 - do ft = 1,numpft_ed + do ft = 1,numpft lai = lai+ sum(currentPatch%tlai_profile(1,ft,:)) enddo @@ -1159,9 +1157,9 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! what is the height of this layer? (for snow burial purposes...) ! EDPftvarcon_inst%vertical_canopy_frac(ft))! fudge - this should be pft specific but i cant get it to compile. layer_top_hite = currentCohort%hite-((iv/currentCohort%NV) * currentCohort%hite * & - EDecophyscon%crown(currentCohort%pft) ) + EDPftvarcon_inst%crown(currentCohort%pft) ) layer_bottom_hite = currentCohort%hite-(((iv+1)/currentCohort%NV) * currentCohort%hite * & - EDecophyscon%crown(currentCohort%pft)) ! EDPftvarcon_inst%vertical_canopy_frac(ft)) + EDPftvarcon_inst%crown(currentCohort%pft)) ! EDPftvarcon_inst%vertical_canopy_frac(ft)) fraction_exposed =1.0_r8 @@ -1192,10 +1190,10 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) iv = currentCohort%NV ! EDPftvarcon_inst%vertical_canopy_frac(ft))! fudge - this should be pft specific but i cant get it to compile. layer_top_hite = currentCohort%hite-((iv/currentCohort%NV) * currentCohort%hite * & - EDecophyscon%crown(currentCohort%pft) ) + EDPftvarcon_inst%crown(currentCohort%pft) ) ! EDPftvarcon_inst%vertical_canopy_frac(ft)) layer_bottom_hite = currentCohort%hite-(((iv+1)/currentCohort%NV) * currentCohort%hite * & - EDecophyscon%crown(currentCohort%pft)) + EDPftvarcon_inst%crown(currentCohort%pft)) fraction_exposed = 1.0_r8 !default. snow_depth_avg = snow_depth_si * frac_sno_eff_si @@ -1251,7 +1249,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo !cohort do L = 1,currentPatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft do iv = 1,currentPatch%nrad(L,ft) !account for total canopy area currentPatch%tlai_profile(L,ft,iv) = currentPatch%tlai_profile(L,ft,iv) / & @@ -1279,7 +1277,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentPatch%nrad = currentPatch%ncan do L = 1,currentPatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft if(currentPatch%nrad(L,ft) > 30)then write(fates_log(), *) 'ED: issue w/ nrad' endif @@ -1291,24 +1289,24 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) end do !iv enddo !ft - if ( L == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft_ed,1))) < 0.99999 & + if ( L == 1 .and. abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999 & .and. currentPatch%NCL_p > 1 ) then - write(fates_log(), *) 'ED: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft_ed,1)) - write(fates_log(), *) 'ED: cohort areas', currentPatch%canopy_area_profile(1,1:numpft_ed,:) + write(fates_log(), *) 'ED: canopy area too small',sum(currentPatch%canopy_area_profile(1,1:numpft,1)) + write(fates_log(), *) 'ED: cohort areas', currentPatch%canopy_area_profile(1,1:numpft,:) endif if (L == 1 .and. currentPatch%NCL_p > 1 .and. & - abs(sum(currentPatch%canopy_area_profile(1,1:numpft_ed,1))) < 0.99999) then + abs(sum(currentPatch%canopy_area_profile(1,1:numpft,1))) < 0.99999) then write(fates_log(), *) 'ED: not enough area in the top canopy', & - sum(currentPatch%canopy_area_profile(L,1:numpft_ed,1)), & - currentPatch%canopy_area_profile(L,1:numpft_ed,1) + sum(currentPatch%canopy_area_profile(L,1:numpft,1)), & + currentPatch%canopy_area_profile(L,1:numpft,1) endif - if(abs(sum(currentPatch%canopy_area_profile(L,1:numpft_ed,1))) > 1.00001)then + if(abs(sum(currentPatch%canopy_area_profile(L,1:numpft,1))) > 1.00001)then write(fates_log(), *) 'ED: canopy-area-profile wrong', & - sum(currentPatch%canopy_area_profile(L,1:numpft_ed,1)), & + sum(currentPatch%canopy_area_profile(L,1:numpft,1)), & currentPatch%patchno, L - write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(L,1:numpft_ed,1),currentPatch%patchno + write(fates_log(), *) 'ED: areas',currentPatch%canopy_area_profile(L,1:numpft,1),currentPatch%patchno currentCohort => currentPatch%shortest @@ -1328,7 +1326,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) enddo ! loop over L do L = 1,currentPatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft if(currentPatch%present(L,FT) > 1)then write(fates_log(), *) 'ED: present issue',L,ft,currentPatch%present(L,FT) currentPatch%present(L,ft) = 1 @@ -1500,28 +1498,28 @@ function calc_areaindex(cpatch,ai_type) result(ai) ai = 0._r8 if (trim(ai_type) == 'elai') then do cl = 1,cpatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'tlai') then do cl = 1,cpatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%tlai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'esai') then do cl = 1,cpatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%esai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'tsai') then do cl = 1,cpatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%tsai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 2fbf2be010..4dccf83b20 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -13,7 +13,6 @@ module EDCohortDynamicsMod use FatesConstantsMod , only : itrue use FatesInterfaceMod , only : hlm_days_per_year use EDPftvarcon , only : EDPftvarcon_inst - use EDEcophysContype , only : EDecophyscon use EDGrowthFunctionsMod , only : c_area, tree_lai use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax @@ -676,8 +675,8 @@ subroutine fuse_cohorts(patchptr, bc_in) type (ed_cohort_type) , pointer :: currentCohort, nextc, nextnextc integer :: i integer :: fusion_took_place - integer :: maxcohorts !maximum total no of cohorts. Needs to be >numpft_edx2 - integer :: iterate !do we need to keep fusing to get below maxcohorts? + integer :: maxcohorts ! maximum total no of cohorts. + integer :: iterate ! do we need to keep fusing to get below maxcohorts? integer :: nocohorts real(r8) :: newn real(r8) :: diff diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 2a19bf51b2..533073b0b8 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -8,17 +8,18 @@ module EDPatchDynamicsMod use EDPftvarcon , only : EDPftvarcon_inst use EDCohortDynamicsMod , only : fuse_cohorts, sort_cohorts, insert_cohort use EDtypesMod , only : ncwd, n_dbh_bins, ntol, area, dbhmax - use EDTypesMod , only : numpft_ed use EDTypesMod , only : maxPatchesPerSite use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : min_patch_area use EDTypesMod , only : nclmax + use EDTypesMod , only : maxpft use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_numlevgrnd use FatesInterfaceMod , only : hlm_numlevsoil use FatesInterfaceMod , only : hlm_numSWb use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_days_per_year + use FatesInterfaceMod , only : numpft use FatesGlobals , only : endrun => fates_endrun use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue @@ -206,8 +207,8 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: age ! notional age of this patch in years integer :: tnull ! is there a tallest cohort? integer :: snull ! is there a shortest cohort? - real(r8) :: root_litter_local(numpft_ed) ! initial value of root litter. KgC/m2 - real(r8) :: leaf_litter_local(numpft_ed) ! initial value of leaf litter. KgC/m2 + real(r8) :: root_litter_local(maxpft) ! initial value of root litter. KgC/m2 + real(r8) :: leaf_litter_local(maxpft) ! initial value of leaf litter. KgC/m2 real(r8) :: cwd_ag_local(ncwd) ! initial value of above ground coarse woody debris. KgC/m2 real(r8) :: cwd_bg_local(ncwd) ! initial value of below ground coarse woody debris. KgC/m2 real(r8) :: spread_local(nclmax) ! initial value of canopy spread parameter.no units @@ -523,7 +524,7 @@ subroutine average_patch_properties( currentPatch, newPatch, patch_site_areadis newPatch%cwd_bg(c) = newPatch%cwd_bg(c) + currentPatch%cwd_bg(c) * patch_site_areadis/newPatch%area enddo - do p = 1,numpft_ed !move litter pool en mass into the new patch + do p = 1,numpft !move litter pool en mass into the new patch newPatch%root_litter(p) = newPatch%root_litter(p) + currentPatch%root_litter(p) * patch_site_areadis/newPatch%area newPatch%leaf_litter(p) = newPatch%leaf_litter(p) + currentPatch%leaf_litter(p) * patch_site_areadis/newPatch%area @@ -589,7 +590,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + burned_litter * new_patch%area !kG/site/day enddo - do p = 1,numpft_ed + do p = 1,numpft burned_litter = new_patch%leaf_litter(p) * patch_site_areadis/new_patch%area * currentPatch%burnt_frac_litter(dl_sf) new_patch%leaf_litter(p) = new_patch%leaf_litter(p) - burned_litter currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day @@ -684,7 +685,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si enddo !burned leaves. - do p = 1,numpft_ed + do p = 1,numpft currentSite%leaf_litter_burned(p) = currentSite%leaf_litter_burned(p) + & dead_tree_density * currentCohort%bl * currentCohort%cfa @@ -763,8 +764,8 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat real(r8) :: np_mult !Fraction of the new patch which came from the current patch (and so needs the same litter) integer :: p,c real(r8) :: canopy_mortality_woody_litter ! flux of wood litter in to litter pool: KgC/m2/day - real(r8) :: canopy_mortality_leaf_litter(numpft_ed) ! flux in to leaf litter from tree death: KgC/m2/day - real(r8) :: canopy_mortality_root_litter(numpft_ed) ! flux in to froot litter from tree death: KgC/m2/day + real(r8) :: canopy_mortality_leaf_litter(maxpft) ! flux in to leaf litter from tree death: KgC/m2/day + real(r8) :: canopy_mortality_root_litter(maxpft) ! flux in to froot litter from tree death: KgC/m2/day real(r8) :: mean_agb_frac ! mean fraction of AGB to total woody biomass (stand mean) !--------------------------------------------------------------------- @@ -830,7 +831,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat ! For the new patch, only some fraction of its land area (patch_areadis/np%area) is derived from the current patch ! so we need to multiply by patch_areadis/np%area - mean_agb_frac = sum(EDPftvarcon_inst%allom_agb_frac(1:numpft_ed))/dble(numpft_ed) + mean_agb_frac = sum(EDPftvarcon_inst%allom_agb_frac(1:numpft))/dble(numpft) do c = 1,ncwd @@ -848,7 +849,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat SF_val_CWD_frac(c) * canopy_mortality_woody_litter * hlm_days_per_year * (1.0_r8 - mean_agb_frac) / AREA enddo - do p = 1,numpft_ed + do p = 1,numpft new_patch%leaf_litter(p) = new_patch%leaf_litter(p) + canopy_mortality_leaf_litter(p) / litter_area * np_mult new_patch%root_litter(p) = new_patch%root_litter(p) + canopy_mortality_root_litter(p) / litter_area * np_mult @@ -897,8 +898,8 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_ allocate(new_patch%fabi(hlm_numSWb)) allocate(new_patch%sabs_dir(hlm_numSWb)) allocate(new_patch%sabs_dif(hlm_numSWb)) - allocate(new_patch%rootfr_ft(numpft_ed,hlm_numlevgrnd)) - allocate(new_patch%rootr_ft(numpft_ed,hlm_numlevgrnd)) + allocate(new_patch%rootfr_ft(numpft,hlm_numlevgrnd)) + allocate(new_patch%rootr_ft(numpft,hlm_numlevgrnd)) call zero_patch(new_patch) !The nan value in here is not working?? @@ -1152,7 +1153,7 @@ subroutine fuse_patches( csite, bc_in ) !---------------------------------------------------------------------! ! Calculate the difference criteria for each pft and dbh class ! !---------------------------------------------------------------------! - do ft = 1,numpft_ed ! loop over pfts + do ft = 1,numpft ! loop over pfts do z = 1,n_dbh_bins ! loop over hgt bins !is there biomass in this category? if(currentPatch%pft_agb_profile(ft,z) > 0.0_r8.or.tpp%pft_agb_profile(ft,z) > 0.0_r8)then @@ -1268,7 +1269,7 @@ subroutine fuse_2_patches(dp, rp) rp%cwd_bg(c) = (dp%cwd_bg(c)*dp%area + rp%cwd_bg(c)*rp%area) * inv_sum_area enddo - do p = 1,numpft_ed + do p = 1,numpft rp%seeds_in(p) = (rp%seeds_in(p)*rp%area + dp%seeds_in(p)*dp%area) * inv_sum_area rp%seed_decay(p) = (rp%seed_decay(p)*rp%area + dp%seed_decay(p)*dp%area) * inv_sum_area rp%seed_germination(p) = (rp%seed_germination(p)*rp%area + dp%seed_germination(p)*dp%area) * inv_sum_area @@ -1531,11 +1532,7 @@ subroutine patch_pft_size_profile(cp_pnt) delta_dbh = (DBHMAX/N_DBH_BINS) - do p = 1,numpft_ed - do j = 1,N_DBH_BINS - currentPatch%pft_agb_profile(p,j) = 0.0_r8 - enddo - enddo + currentPatch%pft_agb_profile(:,:) = 0.0_r8 do j = 1,N_DBH_BINS if (j == 1) then @@ -1617,7 +1614,7 @@ subroutine set_root_fraction( cpatch , zi ) integer :: lev,p,c,ft !---------------------------------------------------------------------- - do ft = 1,numpft_ed + do ft = 1,numpft do lev = 1, hlm_numlevgrnd cpatch%rootfr_ft(ft,lev) = 0._r8 enddo diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 24b37794b3..7bc59ec483 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -11,10 +11,9 @@ module EDPhysiologyMod use FatesInterfaceMod, only : hlm_model_day use FatesInterfaceMod, only : hlm_freq_day use FatesInterfaceMod, only : hlm_day_of_year + use FatesInterfaceMod, only : numpft use FatesConstantsMod, only : r8 => fates_r8 - use EDEcophysContype , only : EDecophyscon use EDPftvarcon , only : EDPftvarcon_inst - use EDEcophysContype , only : EDecophyscon use FatesInterfaceMod, only : bc_in_type use EDCohortDynamicsMod , only : allocate_live_biomass, zero_cohort use EDCohortDynamicsMod , only : create_cohort, sort_cohorts @@ -24,8 +23,8 @@ module EDPhysiologyMod use EDTypesMod , only : external_recruitment use EDTypesMod , only : ncwd use EDTypesMod , only : nlevleaf - use EDTypesMod , only : numpft_ed use EDTypesMod , only : senes + use EDTypesMod , only : maxpft use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use shr_log_mod , only : errMsg => shr_log_errMsg @@ -128,7 +127,7 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in ) call cwd_input(currentPatch) call cwd_out( currentSite, currentPatch, bc_in) - do p = 1,numpft_ed + do p = 1,numpft currentSite%dseed_dt(p) = currentSite%dseed_dt(p) + & (currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - & currentPatch%seed_germination(p)) * currentPatch%area/AREA @@ -139,7 +138,7 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in ) currentPatch%dcwd_BG_dt(c) = currentPatch%cwd_BG_in(c) - currentPatch%cwd_BG_out(c) enddo - do p = 1,numpft_ed + do p = 1,numpft currentPatch%dleaf_litter_dt(p) = currentPatch%leaf_litter_in(p) - & currentPatch%leaf_litter_out(p) currentPatch%droot_litter_dt(p) = currentPatch%root_litter_in(p) - & @@ -193,14 +192,14 @@ subroutine trim_canopy( currentSite ) currentCohort%leaf_cost = 1._r8/(EDPftvarcon_inst%slatop(currentCohort%pft)*1000.0_r8) currentCohort%leaf_cost = currentCohort%leaf_cost + & 1.0_r8/(EDPftvarcon_inst%slatop(currentCohort%pft)*1000.0_r8) * & - EDPftvarcon_inst%allom_l2fr(currentCohort%pft) / EDecophyscon%root_long(currentCohort%pft) + EDPftvarcon_inst%allom_l2fr(currentCohort%pft) / EDPftvarcon_inst%root_long(currentCohort%pft) currentCohort%leaf_cost = currentCohort%leaf_cost * (EDPftvarcon_inst%grperc(currentCohort%pft) + 1._r8) else !evergreen costs currentCohort%leaf_cost = 1.0_r8/(EDPftvarcon_inst%slatop(currentCohort%pft)* & EDPftvarcon_inst%leaf_long(currentCohort%pft)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1 currentCohort%leaf_cost = currentCohort%leaf_cost + & 1.0_r8/(EDPftvarcon_inst%slatop(currentCohort%pft)*1000.0_r8) * & - EDPftvarcon_inst%allom_l2fr(currentCohort%pft) / EDecophyscon%root_long(currentCohort%pft) + EDPftvarcon_inst%allom_l2fr(currentCohort%pft) / EDPftvarcon_inst%root_long(currentCohort%pft) currentCohort%leaf_cost = currentCohort%leaf_cost * (EDPftvarcon_inst%grperc(currentCohort%pft) + 1._r8) endif if (currentCohort%year_net_uptake(z) < currentCohort%leaf_cost)then @@ -211,7 +210,7 @@ subroutine trim_canopy( currentSite ) endif ! keep trimming until none of the canopy is in negative carbon balance. - if (currentCohort%hite > EDecophyscon%hgt_min(currentCohort%pft))then + if (currentCohort%hite > EDPftvarcon_inst%hgt_min(currentCohort%pft))then currentCohort%canopy_trim = currentCohort%canopy_trim - EDPftvarcon_inst%trim_inc(currentCohort%pft) if (EDPftvarcon_inst%evergreen(currentCohort%pft) /= 1)then currentCohort%laimemory = currentCohort%laimemory*(1.0_r8 - EDPftvarcon_inst%trim_inc(currentCohort%pft)) @@ -630,7 +629,7 @@ subroutine seeds_in( currentSite, cp_pnt ) type(ed_patch_type), pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort integer :: p - logical :: pft_present(numpft_ed) + logical :: pft_present(maxpft) real(r8) :: npfts_present !---------------------------------------------------------------------- @@ -661,7 +660,7 @@ subroutine seeds_in( currentSite, cp_pnt ) currentPatch => cp_pnt currentCohort => currentPatch%tallest do while (associated(currentCohort)) - do p = 1, numpft_ed + do p = 1, numpft if (pft_present(p)) then currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + currentCohort%seed_prod * currentCohort%n / & (currentPatch%area * npfts_present) @@ -686,11 +685,11 @@ subroutine seeds_in( currentSite, cp_pnt ) do while(associated(currentPatch)) if (external_recruitment == 1) then !external seed rain - needed to prevent extinction - do p = 1,numpft_ed + do p = 1,numpft currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + & - EDecophyscon%seed_rain(p) !KgC/m2/year + EDPftvarcon_inst%seed_rain(p) !KgC/m2/year currentSite%seed_rain_flux(p) = currentSite%seed_rain_flux(p) + & - EDecophyscon%seed_rain(p) * currentPatch%area/AREA !KgC/m2/year + EDPftvarcon_inst%seed_rain(p) * currentPatch%area/AREA !KgC/m2/year enddo endif currentPatch => currentPatch%younger @@ -718,7 +717,7 @@ subroutine seed_decay( currentSite, currentPatch ) ! default value from Liscke and Loffler 2006 ; making this a PFT-specific parameter ! decays the seed pool according to exponential model ! seed_decay_turnover is in yr-1 - do p = 1,numpft_ed + do p = 1,numpft currentPatch%seed_decay(p) = currentSite%seed_bank(p) * EDPftvarcon_inst%seed_decay_turnover(p) enddo @@ -747,7 +746,7 @@ subroutine seed_germination( currentSite, currentPatch ) ! germination_timescale is being pulled to PFT parameter; units are 1/yr ! thus the mortality rate of seed -> recruit (in units of carbon) is seed_decay_turnover(p)/germination_timescale(p) ! and thus the mortlaity rate (in units of individuals) is the product of that times the ratio of (hypothetical) seed mass to recruit biomass - do p = 1,numpft_ed + do p = 1,numpft currentPatch%seed_germination(p) = min(currentSite%seed_bank(p) * & EDPftvarcon_inst%germination_timescale(p),max_germination) enddo @@ -826,7 +825,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! Maintenance demands if (EDPftvarcon_inst%evergreen(currentCohort%pft) == 1)then !grass and EBT currentCohort%leaf_md = currentCohort%bl / EDPftvarcon_inst%leaf_long(currentCohort%pft) - currentCohort%root_md = currentCohort%br / EDecophyscon%root_long(currentCohort%pft) + currentCohort%root_md = currentCohort%br / EDPftvarcon_inst%root_long(currentCohort%pft) currentCohort%md = currentCohort%root_md + currentCohort%leaf_md endif @@ -836,13 +835,13 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !are still in an expansion phase. if (EDPftvarcon_inst%season_decid(currentCohort%pft) == 1)then - currentCohort%root_md = currentCohort%br /EDecophyscon%root_long(currentCohort%pft) + currentCohort%root_md = currentCohort%br /EDPftvarcon_inst%root_long(currentCohort%pft) currentCohort%leaf_md = 0._r8 currentCohort%md = currentCohort%root_md + currentCohort%leaf_md endif if (EDPftvarcon_inst%stress_decid(currentCohort%pft) == 1)then - currentCohort%root_md = currentCohort%br /EDecophyscon%root_long(currentCohort%pft) + currentCohort%root_md = currentCohort%br /EDPftvarcon_inst%root_long(currentCohort%pft) currentCohort%leaf_md = 0._r8 currentCohort%md = currentCohort%root_md + currentCohort%leaf_md endif @@ -864,17 +863,17 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! this is the fraction of maintenance demand we -have- to do... if ( DEBUG ) write(fates_log(),*) 'EDphys 760 ',currentCohort%npp_acc_hold, currentCohort%md, & - EDecophyscon%leaf_stor_priority(currentCohort%pft) + EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft) currentCohort%carbon_balance = currentCohort%npp_acc_hold - & - currentCohort%md * EDecophyscon%leaf_stor_priority(currentCohort%pft) + currentCohort%md * EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft) ! Allowing only carbon from NPP pool to account for npp flux into the maintenance turnover pools ! ie this does not include any use of storage carbon or balive to make up for missing carbon balance in the transfer currentCohort%npp_leaf = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%leaf_md/currentCohort%md, & - currentCohort%leaf_md*EDecophyscon%leaf_stor_priority(currentCohort%pft))) + currentCohort%leaf_md*EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft))) currentCohort%npp_froot = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%root_md/currentCohort%md, & - currentCohort%root_md*EDecophyscon%leaf_stor_priority(currentCohort%pft))) + currentCohort%root_md*EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft))) if (Bleaf(currentCohort) > 0._r8)then @@ -883,7 +882,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) if (currentCohort%carbon_balance > 0._r8)then !spend C on growing and storing !what fraction of the target storage do we have? - frac = max(0.0_r8,currentCohort%bstore/(Bleaf(currentCohort) * EDecophyscon%cushion(currentCohort%pft))) + frac = max(0.0_r8,currentCohort%bstore/(Bleaf(currentCohort) * EDPftvarcon_inst%cushion(currentCohort%pft))) ! FIX(SPM,080514,fstore never used ) f_store = max(exp(-1.*frac**4._r8) - exp( -1.0_r8 ),0.0_r8) !what fraction of allocation do we divert to storage? @@ -915,14 +914,14 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !Do we have enough carbon left over to make up the rest of the turnover demand? balive_loss = 0._r8 - if (currentCohort%carbon_balance > currentCohort%md*(1.0_r8- EDecophyscon%leaf_stor_priority(currentCohort%pft)))then ! Yes... + if (currentCohort%carbon_balance > currentCohort%md*(1.0_r8- EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)))then ! Yes... currentCohort%carbon_balance = currentCohort%carbon_balance - currentCohort%md * (1.0_r8 - & - EDecophyscon%leaf_stor_priority(currentCohort%pft)) + EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)) currentCohort%npp_leaf = currentCohort%npp_leaf + & - currentCohort%leaf_md * (1.0_r8-EDecophyscon%leaf_stor_priority(currentCohort%pft)) + currentCohort%leaf_md * (1.0_r8-EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)) currentCohort%npp_froot = currentCohort%npp_froot + & - currentCohort%root_md * (1.0_r8-EDecophyscon%leaf_stor_priority(currentCohort%pft)) + currentCohort%root_md * (1.0_r8-EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft)) else ! we can't maintain constant leaf area and root area. Balive is reduced @@ -931,7 +930,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) currentCohort%npp_froot = currentCohort%npp_froot + & max(0.0_r8,currentCohort%carbon_balance*(currentCohort%root_md/currentCohort%md)) - balive_loss = currentCohort%md *(1.0_r8- EDecophyscon%leaf_stor_priority(currentCohort%pft))- currentCohort%carbon_balance + balive_loss = currentCohort%md *(1.0_r8- EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft))- currentCohort%carbon_balance currentCohort%carbon_balance = 0._r8 endif @@ -943,7 +942,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !only if carbon balance is +ve if ((currentCohort%balive >= target_balive).AND.(currentCohort%carbon_balance > 0._r8))then ! fraction of carbon going into active vs structural carbon - if (currentCohort%dbh <= EDecophyscon%max_dbh(currentCohort%pft))then ! cap on leaf biomass + if (currentCohort%dbh <= EDPftvarcon_inst%max_dbh(currentCohort%pft))then ! cap on leaf biomass dbldbd = dDbhdBd(currentCohort)/dDbhdBl(currentCohort) dbrdbd = EDPftvarcon_inst%allom_l2fr(currentCohort%pft) * dbldbd dhdbd_fn = dhdbd(currentCohort) @@ -951,12 +950,12 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) u = 1.0_r8 / (dbldbd + dbrdbd + dbswdbd) va = 1.0_r8 / (1.0_r8 + u) vs = u / (1.0_r8 + u) - gr_fract = 1.0_r8 - EDecophyscon%seed_alloc(currentCohort%pft) + gr_fract = 1.0_r8 - EDPftvarcon_inst%seed_alloc(currentCohort%pft) else dbldbd = 0._r8; dbrdbd = 0._r8 ;dbswdbd = 0._r8 va = 0.0_r8 vs = 1.0_r8 - gr_fract = 1.0_r8 - (EDecophyscon%seed_alloc(currentCohort%pft) + EDecophyscon%clone_alloc(currentCohort%pft)) + gr_fract = 1.0_r8 - (EDPftvarcon_inst%seed_alloc(currentCohort%pft) + EDPftvarcon_inst%clone_alloc(currentCohort%pft)) endif !FIX(RF,032414) - to fix high bl's. needed to prevent numerical errors without the ODEINT. @@ -986,7 +985,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) write(fates_log(),*) 'cohort fluxes',currentCohort%pft,currentCohort%canopy_layer,currentCohort%n, & currentCohort%npp_acc_hold,currentCohort%dbalivedt,balive_loss, & currentCohort%dbdeaddt,currentCohort%dbstoredt,currentCohort%seed_prod,currentCohort%md * & - EDecophyscon%leaf_stor_priority(currentCohort%pft) + EDPftvarcon_inst%leaf_stor_priority(currentCohort%pft) write(fates_log(),*) 'proxies' ,target_balive,currentCohort%balive,currentCohort%dbh,va,vs,gr_fract endif @@ -1037,16 +1036,16 @@ subroutine recruitment( t, currentSite, currentPatch, bc_in ) allocate(temp_cohort) ! create temporary cohort call zero_cohort(temp_cohort) - do ft = 1,numpft_ed + do ft = 1,numpft temp_cohort%canopy_trim = 0.8_r8 !starting with the canopy not fully expanded temp_cohort%pft = ft - temp_cohort%hite = EDecophyscon%hgt_min(ft) + temp_cohort%hite = EDPftvarcon_inst%hgt_min(ft) temp_cohort%dbh = Dbh(temp_cohort) temp_cohort%bdead = Bdead(temp_cohort) temp_cohort%balive = Bleaf(temp_cohort)*(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) & + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite) - temp_cohort%bstore = EDecophyscon%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) & + temp_cohort%bstore = EDPftvarcon_inst%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) & + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite)) temp_cohort%n = currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day & / (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore) @@ -1170,7 +1169,7 @@ subroutine CWD_Input( currentPatch) enddo ! end loop over cohorts - do p = 1,numpft_ed + do p = 1,numpft currentPatch%leaf_litter_in(p) = currentPatch%leaf_litter_in(p) + currentPatch%seed_decay(p) !KgC/m2/yr enddo @@ -1238,7 +1237,7 @@ subroutine fragmentation_scaler( currentPatch, bc_in) !BTRAN APPROACH - is quite simple, but max's out decomp at all unstressed !soil moisture values, which is not realistic. !litter decomp is proportional to water limitation on average... - w_scalar = sum(currentPatch%btran_ft(1:numpft_ed))/numpft_ed + w_scalar = sum(currentPatch%btran_ft(1:numpft))/numpft currentPatch%fragmentation_scaler = min(1.0_r8,max(0.0_r8,t_scalar * w_scalar)) @@ -1274,8 +1273,8 @@ subroutine cwd_out( currentSite, currentPatch, bc_in ) currentPatch%cwd_ag_out(1:ncwd) = 0.0_r8 currentPatch%cwd_bg_out(1:ncwd) = 0.0_r8 - currentPatch%leaf_litter_out(1:numpft_ed) = 0.0_r8 - currentPatch%root_litter_out(1:numpft_ed) = 0.0_r8 + currentPatch%leaf_litter_out(:) = 0.0_r8 + currentPatch%root_litter_out(:) = 0.0_r8 do c = 1,ncwd currentPatch%cwd_ag_out(c) = max(0.0_r8, currentPatch%cwd_ag(c) * & @@ -1289,7 +1288,7 @@ subroutine cwd_out( currentSite, currentPatch, bc_in ) ! thick leaves can dry out before they are decomposed, for example. ! this section needs further scientific input. - do ft = 1,numpft_ed + do ft = 1,numpft currentPatch%leaf_litter_out(ft) = max(0.0_r8,currentPatch%leaf_litter(ft)* SF_val_max_decomp(dl_sf) * & currentPatch%fragmentation_scaler ) currentPatch%root_litter_out(ft) = max(0.0_r8,currentPatch%root_litter(ft)* SF_val_max_decomp(dl_sf) * & @@ -1331,7 +1330,6 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use EDTypesMod, only : AREA - use EDTypesMod, only : numpft_ed use FatesInterfaceMod, only : hlm_numlevdecomp_full use FatesInterfaceMod, only : hlm_numlevdecomp use EDPftvarcon, only : EDPftvarcon_inst @@ -1361,11 +1359,11 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) integer :: begp,endp integer :: begc,endc !bounds !------------------------------------------------------------------------ - real(r8) :: cinput_rootfr(1:numpft_ed, 1:hlm_numlevdecomp_full) ! column by pft root fraction used for calculating inputs + real(r8) :: cinput_rootfr(1:maxpft, 1:hlm_numlevdecomp_full) ! column by pft root fraction used for calculating inputs real(r8) :: croot_prof_perpatch(1:hlm_numlevdecomp_full) real(r8) :: surface_prof(1:hlm_numlevdecomp_full) integer :: ft - real(r8) :: rootfr_tot(1:numpft_ed), biomass_bg_ft(1:numpft_ed) + real(r8) :: rootfr_tot(1:maxpft), biomass_bg_ft(1:maxpft) real(r8) :: surface_prof_tot, leaf_prof_sum, stem_prof_sum, froot_prof_sum, biomass_bg_tot real(r8) :: delta @@ -1392,7 +1390,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) integer, parameter :: rooting_profile_varindex_water = 1 real(r8) :: leaf_prof(1:nsites, 1:hlm_numlevdecomp) - real(r8) :: froot_prof(1:nsites, 1:numpft_ed, 1:hlm_numlevdecomp) + real(r8) :: froot_prof(1:nsites, 1:maxpft, 1:hlm_numlevdecomp) real(r8) :: croot_prof(1:nsites, 1:hlm_numlevdecomp) real(r8) :: stem_prof(1:nsites, 1:hlm_numlevdecomp) @@ -1418,7 +1416,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! initialize profiles to zero leaf_prof(1:nsites, :) = 0._r8 - froot_prof(1:nsites, 1:numpft_ed, :) = 0._r8 + froot_prof(1:nsites, 1:maxpft, :) = 0._r8 stem_prof(1:nsites, :) = 0._r8 do s = 1,nsites @@ -1428,20 +1426,20 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) surface_prof(j) = exp(-surfprof_exp * bc_in(s)%z_sisl(j)) / bc_in(s)%dz_decomp_sisl(j) end do - cinput_rootfr(1:numpft_ed, :) = 0._r8 + cinput_rootfr(:,:) = 0._r8 ! calculate pft-specific rooting profiles in the absence of permafrost or bedrock limitations if ( exponential_rooting_profile ) then if ( .not. pftspecific_rootingprofile ) then ! define rooting profile from exponential parameters - do ft = 1, numpft_ed + do ft = 1, numpft do j = 1, hlm_numlevdecomp cinput_rootfr(ft,j) = exp(-rootprof_exp * bc_in(s)%z_sisl(j)) / bc_in(s)%dz_decomp_sisl(j) end do end do else ! use beta distribution parameter from Jackson et al., 1996 - do ft = 1, numpft_ed + do ft = 1, numpft do j = 1, hlm_numlevdecomp cinput_rootfr(ft,j) = & ( EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & @@ -1453,7 +1451,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end do endif else - do ft = 1,numpft_ed + do ft = 1,numpft do j = 1, hlm_numlevdecomp ! use standard CLM root fraction profiles; cinput_rootfr(ft,j) = ( .5_r8*( & @@ -1469,22 +1467,21 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! ! now add permafrost constraint: integrate rootfr over active layer of soil site, ! truncate below permafrost or bedrock table where present, and rescale so that integral = 1 - do ft = 1,numpft_ed - rootfr_tot(ft) = 0._r8 - end do + rootfr_tot(:) = 0._r8 + surface_prof_tot = 0._r8 ! do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), hlm_numlevdecomp) surface_prof_tot = surface_prof_tot + surface_prof(j) * bc_in(s)%dz_decomp_sisl(j) end do - do ft = 1,numpft_ed + do ft = 1,numpft do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), hlm_numlevdecomp) rootfr_tot(ft) = rootfr_tot(ft) + cinput_rootfr(ft,j) * bc_in(s)%dz_decomp_sisl(j) end do end do ! ! rescale the fine root profile - do ft = 1,numpft_ed + do ft = 1,numpft if ( (bc_in(s)%max_rooting_depth_index_col > 0) .and. (rootfr_tot(ft) > 0._r8) ) then ! where there is not permafrost extending to the surface, integrate the profiles over the active layer ! this is equivalent to integrating over all soil layers outside of permafrost regions @@ -1521,7 +1518,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! for one layer decomposition model, set profiles to unity leaf_prof(1:nsites, :) = 1._r8 - froot_prof(1:nsites, 1:numpft_ed, :) = 1._r8 + froot_prof(1:nsites, 1:numpft, :) = 1._r8 stem_prof(1:nsites, :) = 1._r8 end if @@ -1546,7 +1543,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) call endrun(msg=errMsg(sourcefile, __LINE__)) endif ! now check each fine root profile - do ft = 1,numpft_ed + do ft = 1,numpft froot_prof_sum = 0._r8 do j = 1, hlm_numlevdecomp froot_prof_sum = froot_prof_sum + froot_prof(s,ft,j) * bc_in(s)%dz_decomp_sisl(j) @@ -1583,7 +1580,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! the CWD pools lose information about which PFT they came from; for the stems this doesn't matter as they all have the same profile, ! however for the coarse roots they may have different profiles. to approximately recover this information, loop over all cohorts in patch ! to calculate the total root biomass in that patch of each pft, and then rescale the croot_prof as the weighted average of the froot_prof - biomass_bg_ft(1:numpft_ed) = 0._r8 + biomass_bg_ft(:) = 0._r8 currentCohort => currentPatch%tallest do while(associated(currentCohort)) biomass_bg_ft(currentCohort%pft) = biomass_bg_ft(currentCohort%pft) + & @@ -1592,7 +1589,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) enddo !currentCohort ! biomass_bg_tot = 0._r8 - do ft = 1,numpft_ed + do ft = 1,numpft biomass_bg_tot = biomass_bg_tot + biomass_bg_ft(ft) end do ! @@ -1602,7 +1599,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end do ! if ( biomass_bg_tot .gt. 0._r8) then - do ft = 1,numpft_ed + do ft = 1,numpft do j = 1, hlm_numlevdecomp croot_prof_perpatch(j) = croot_prof_perpatch(j) + froot_prof(s,ft,j) * biomass_bg_ft(ft) / biomass_bg_tot end do @@ -1623,7 +1620,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! write(fates_log(),*)'cdk CWD_AG_out', c, currentpatch%CWD_AG_out(c), ED_val_cwd_fcel, currentpatch%area/AREA ! write(fates_log(),*)'cdk CWD_BG_out', c, currentpatch%CWD_BG_out(c), ED_val_cwd_fcel, currentpatch%area/AREA ! end do - ! do ft = 1,numpft_ed + ! do ft = 1,numpft ! write(fates_log(),*)'cdk leaf_litter_out', ft, currentpatch%leaf_litter_out(ft), ED_val_cwd_fcel, currentpatch%area/AREA ! write(fates_log(),*)'cdk root_litter_out', ft, currentpatch%root_litter_out(ft), ED_val_cwd_fcel, currentpatch%area/AREA ! end do @@ -1644,7 +1641,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end do ! leaf and fine root pools. - do ft = 1,numpft_ed + do ft = 1,numpft do j = 1, hlm_numlevdecomp bc_out(s)%FATES_c_to_litr_lab_c_col(j) = bc_out(s)%FATES_c_to_litr_lab_c_col(j) + & currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(s,j) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 602df3fe3e..47a70b16df 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -11,11 +11,12 @@ module EDBtranMod use EDTypesMod , only : ed_site_type, & ed_patch_type, & ed_cohort_type, & - numpft_ed + maxpft use FatesInterfaceMod , only : hlm_numlevgrnd use shr_kind_mod , only : r8 => shr_kind_r8 use FatesInterfaceMod , only : bc_in_type, & - bc_out_type + bc_out_type, & + numpft use FatesInterfaceMod , only : hlm_use_planthydro use FatesGlobals , only : fates_log @@ -111,9 +112,10 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) real(r8) :: smp_node ! matrix potential real(r8) :: rresis ! suction limitation to transpiration independent ! of root density - real(r8) :: pftgs(numpft_ed) ! pft weighted stomatal conductance s/m + real(r8) :: pftgs(maxpft) ! pft weighted stomatal conductance s/m real(r8) :: temprootr real(r8) :: balive_patch + real(r8) :: sum_pftgs ! sum of weighted conductances (for normalization) !------------------------------------------------------------------------------ associate( & @@ -132,7 +134,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK) - do ft = 1,numpft_ed + do ft = 1,numpft cpatch%btran_ft(ft) = 0.0_r8 do j = 1,hlm_numlevgrnd @@ -184,17 +186,18 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) ! sink term across the different layers in driver/host. Photosynthesis will ! pass the host a total transpiration for the patch. This needs rootr to be ! distributed over the soil layers. - + sum_pftgs = sum(pftgs(1:numpft)) + do j = 1,hlm_numlevgrnd bc_out(s)%rootr_pagl(ifp,j) = 0._r8 - do ft = 1,numpft_ed - if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail + do ft = 1,numpft + if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail !bit-retart test as a result? FIX(RF,032414) bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & - cpatch%rootr_ft(ft,j) * pftgs(ft)/sum(pftgs) + cpatch%rootr_ft(ft,j) * pftgs(ft)/sum_pftgs else bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & - cpatch%rootr_ft(ft,j) * 1./numpft_ed + cpatch%rootr_ft(ft,j) * 1./numpft end if enddo enddo @@ -206,12 +209,12 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) if(hlm_use_planthydro.eq.ifalse) then !weight patch level output BTRAN for the bc_out(s)%btran_pa(ifp) = 0.0_r8 - do ft = 1,numpft_ed - if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail + do ft = 1,numpft + if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * pftgs(ft)/sum(pftgs) + bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * pftgs(ft)/sum_pftgs else - bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft_ed + bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft end if enddo end if @@ -219,7 +222,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) temprootr = sum(bc_out(s)%rootr_pagl(ifp,1:hlm_numlevgrnd)) if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then - write(fates_log(),*) 'error with rootr in canopy fluxes',temprootr,sum(pftgs) + write(fates_log(),*) 'error with rootr in canopy fluxes',temprootr,sum_pftgs do j = 1,hlm_numlevgrnd bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr enddo @@ -238,117 +241,5 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) end subroutine btran_ed - ! ========================================================================================= - - !--------------------------------------------------------------------------------------- - ! SPA based recalculation of BTRAN and water uptake. - !--------------------------------------------------------------------------------------- - -! if (SPA_soil) then ! normal case don't run this. -! rootr(p,:) = 0._r8 -! do FT = 1,numpft_ed - -! ! Soil Physics -! do j = 1,nlevgrnd -! ! CLM water retention curve. Clapp and Hornberger equation. -! s1 = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8) -! s1 = min(1.0_r8,s1) -! smp_node = -sucsat(c,j)*s1**(-bsw(c,j)) -! swp_mpa(j) = smp_node *10.0_r8/1000000.0_r8 !convert from mm to Mpa - -! ! CLM hydraulic conductivity curve. -! ! As opposed to the Richard's equation solution in SoilHydrology.Mod -! ! the conductivity here is defined in the middle of the layer in question, not at the edge... -! xksat = 0.0070556_r8 * (10._r8**(-0.884_r8+0.0153_r8*sand(p)) ) -! hk(j) = xksat*s1**(2._r8*bsw(c,j)+2._r8) !removed the ice from here to avoid 1st ts crashing -! enddo - -! ! Root resistance -! rootxsecarea=3.14159*rootrad**2 -! do j = 1,nlevgrnd -! rootmass(j) = EDecophyscon%soilbeta(FT) * cpatch%rootfr_ft(FT,j) -! rootlength(j) = rootmass(j)/(rootdens*rootxsecarea) !m m-3 soil -! Lsoil(j) = hk(j)/1000/head !converts from mms-1 to ms-1 and then to m2 s-1 MPa-1 -! if(Lsoil(j) < 1e-35_r8.or.cpatch%rootfr_ft(ft,j) <= 0.0_r8)then !prevent floating point error -! soilr_z(j) = 1e35_r8 -! soilr2(j) = 1e35_r8 -! else -! ! Soil-to-root water uptake from Newman (1969). -! rs = sqrt (1._r8 / (rootlength(j) * pi)) -! soilr1(j) = log(rs/rootrad) / (2.0_r8 * pi * rootlength(j) * Lsoil(j) * dz(c,j)) -! ! convert from MPa s m2 m-3 to MPa s m2 mmol-1 -! soilr1(j) = soilr1(j) * 1E-6_r8 * 18_r8 * 0.001_r8 -! ! second component of below ground resistance is related to root hydraulics -! soilr2(j) = EDecophyscon%rootresist(FT)/(rootmass(j)*dz(c,j)) -! soilr_z(j) = soilr1(j)+soilr2(j) -! end if -! enddo - - ! Aggregate soil layers -! totestevap=0._r8 -! weighted_SWP=0._r8 -! estevap=0._r8 -! fraction_uptake=0._r8 -! canopy_soil_resistance=0._r8 !Reset Counters -! totmaxevap = 0._r8 - - ! Estimated max transpiration from LWP gradient / soil resistance -! do j = 1,nlevgrnd -! estevap(j) = (swp_mpa(j) - minlwp)/(soilr_z(j)) -! estevap(j) = max(0._r8,estevap(j)) ! no negative uptake -! maxevap(j) = (0.0_r8 - minlwp)/(soilr2(j)) -! enddo -! totestevap = sum(estevap) -! totmaxevap = sum(maxevap) - - ! Weighted soil water potential -! do j = 1,nlevgrnd -! if(totestevap > 0._r8)then -! fraction_uptake(j) = estevap(j)/totestevap !Fraction of total ET taken from this soil layer -! else -! estevap(j) = 0._r8 -! fraction_uptake(j)=1._r8/nlevgrnd -! end if -! weighted_SWP = weighted_SWP + swp_mpa(j) * estevap(j) -! enddo - -! if(totestevap > 0._r8)then -! weighted_swp = weighted_swp/totestevap -! ! weight SWP for the total evaporation -! else -! write(fates_log(),*) 'empty soil', totestevap -! ! error check -! weighted_swp = minlwp -! end if - - ! Weighted soil-root resistance. Aggregate the conductances (1/soilR) for each soil layer -! do iv = 1,nv !leaf layers -! fleaf = 1.0_r8/nv -! do j = 1,nlevgrnd !root layers -! ! Soil resistance for each canopy layer is related to leaf area -! ! The conductance of the root system to the -! ! whole canopy is reduced by the fraction of leaves in this layer... -! canopy_soil_resistance(iv) = canopy_soil_resistance(iv)+fleaf * 1.0_r8/(soilr_z(j)) -! enddo -! ! Turn aggregated conductance back into resistance. mmol MPa-1 s-1 m-2 to MPa s m2 mmol-1 -! canopy_soil_resistance(iv) = 1./canopy_soil_resistance(iv) -! enddo -! -! cpatch%btran_ft(FT) = totestevap/totmaxevap -! do j = 1,nlevgrnd -! if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail -! !bit-retart test as a result? FIX(RF,032414) -! rootr(p,j) = rootr(p,j) + fraction_uptake(j) * pftgs(ft)/sum(pftgs) -! else -! rootr(p,j) = rootr(p,j) + fraction_uptake(j) * 1./numpft_ed -! end if -! enddo -! enddo !pft loop -! end if ! - !--------------------------------------------------------------------------------------- - ! end of SPA based recalculation of BTRAN and water uptake. - !--------------------------------------------------------------------------------------- - - end module EDBtranMod diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index 6da2a28040..7ee743cba1 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -11,15 +11,15 @@ module EDSurfaceRadiationMod #include "shr_assert.h" use EDTypesMod , only : ed_patch_type, ed_site_type - use EDTypesMod , only : numpft_ed use EDTypesMod , only : maxPatchesPerSite + use EDTypesMod , only : maxpft use FatesConstantsMod , only : r8 => fates_r8 use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : bc_out_type use FatesInterfaceMod , only : hlm_numSWb + use FatesInterfaceMod , only : numpft use EDTypesMod , only : maxSWb use EDTypesMod , only : nclmax - use EDTypesMod , only : numpft_ed use EDTypesMod , only : nlevleaf use EDCanopyStructureMod, only: calc_areaindex use FatesGlobals , only : fates_log @@ -74,10 +74,10 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) real(r8) :: sb real(r8) :: error ! Error check real(r8) :: down_rad, up_rad ! Iterative solution do Dif_dn and Dif_up - real(r8) :: ftweight(nclmax,numpft_ed,nlevleaf) - real(r8) :: k_dir(numpft_ed) ! Direct beam extinction coefficient - real(r8) :: tr_dir_z(nclmax,numpft_ed,nlevleaf) ! Exponential transmittance of direct beam radiation through a single layer - real(r8) :: tr_dif_z(nclmax,numpft_ed,nlevleaf) ! Exponential transmittance of diffuse radiation through a single layer + real(r8) :: ftweight(nclmax,maxpft,nlevleaf) + real(r8) :: k_dir(maxpft) ! Direct beam extinction coefficient + real(r8) :: tr_dir_z(nclmax,maxpft,nlevleaf) ! Exponential transmittance of direct beam radiation through a single layer + real(r8) :: tr_dif_z(nclmax,maxpft,nlevleaf) ! Exponential transmittance of diffuse radiation through a single layer real(r8) :: forc_dir(maxPatchesPerSite,maxSWb) real(r8) :: forc_dif(maxPatchesPerSite,maxSWb) real(r8) :: weighted_dir_tr(nclmax) @@ -85,26 +85,30 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) real(r8) :: weighted_dif_ratio(nclmax,maxSWb) real(r8) :: weighted_dif_down(nclmax) real(r8) :: weighted_dif_up(nclmax) - real(r8) :: refl_dif(nclmax,numpft_ed,nlevleaf,maxSWb) ! Term for diffuse radiation reflected by laye - real(r8) :: tran_dif(nclmax,numpft_ed,nlevleaf,maxSWb) ! Term for diffuse radiation transmitted by layer - real(r8) :: dif_ratio(nclmax,numpft_ed,nlevleaf,maxSWb) ! Ratio of upward to forward diffuse fluxes - real(r8) :: Dif_dn(nclmax,numpft_ed,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) - real(r8) :: Dif_up(nclmax,numpft_ed,nlevleaf) ! Upward diffuse flux above canopy layer J (W/m**2 ground area) - real(r8) :: lai_change(nclmax,numpft_ed,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) - real(r8) :: f_not_abs(numpft_ed,maxSWb) ! Fraction reflected + transmitted. 1-absorbtion. - real(r8) :: Abs_dir_z(numpft_ed,nlevleaf) - real(r8) :: Abs_dif_z(numpft_ed,nlevleaf) + real(r8) :: refl_dif(nclmax,maxpft,nlevleaf,maxSWb) ! Term for diffuse radiation reflected by laye + real(r8) :: tran_dif(nclmax,maxpft,nlevleaf,maxSWb) ! Term for diffuse radiation transmitted by layer + real(r8) :: dif_ratio(nclmax,maxpft,nlevleaf,maxSWb) ! Ratio of upward to forward diffuse fluxes + real(r8) :: Dif_dn(nclmax,maxpft,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) + real(r8) :: Dif_up(nclmax,maxpft,nlevleaf) ! Upward diffuse flux above canopy layer J (W/m**2 ground area) + real(r8) :: lai_change(nclmax,maxpft,nlevleaf) ! Forward diffuse flux onto canopy layer J (W/m**2 ground area) + real(r8) :: f_not_abs(maxpft,maxSWb) ! Fraction reflected + transmitted. 1-absorbtion. + real(r8) :: Abs_dir_z(maxpft,nlevleaf) + real(r8) :: Abs_dif_z(maxpft,nlevleaf) real(r8) :: abs_rad(maxSWb) !radiation absorbed by soil real(r8) :: tr_soili ! Radiation transmitted to the soil surface. real(r8) :: tr_soild ! Radiation transmitted to the soil surface. - real(r8) :: phi1b(maxPatchesPerSite,numpft_ed) ! Radiation transmitted to the soil surface. - real(r8) :: phi2b(maxPatchesPerSite,numpft_ed) + real(r8) :: phi1b(maxPatchesPerSite,maxpft) ! Radiation transmitted to the soil surface. + real(r8) :: phi2b(maxPatchesPerSite,maxpft) real(r8) :: laisum ! cumulative lai+sai for canopy layer (at middle of layer) real(r8) :: angle real(r8),parameter :: tolerance = 0.000000001_r8 real(r8), parameter :: pi = 3.141592654 ! PI + + integer, parameter :: max_diag_nlevleaf = 4 + integer, parameter :: diag_nlevleaf = min(nlevleaf,max_diag_nlevleaf) ! for diagnostics, write a small number of leaf layers + real(r8) :: denom real(r8) :: lai_reduction(2) @@ -192,7 +196,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) ! Is this pft/canopy layer combination present in this patch? do L = 1,nclmax - do ft = 1,numpft_ed + do ft = 1,numpft currentPatch%present(L,ft) = 0 do iv = 1, currentPatch%nrad(L,ft) if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then @@ -219,7 +223,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) !Extract information that needs to be provided by ED into local array. ftweight(:,:,:) = 0._r8 do L = 1,currentPatch%NCL_p - do ft = 1,numpft_ed + do ft = 1,numpft do iv = 1, currentPatch%nrad(L,ft) !this is already corrected for area in CLAP ftweight(L,ft,iv) = currentPatch%canopy_area_profile(L,ft,iv) @@ -237,7 +241,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) ! Direct beam extinction coefficient, k_dir. PFT specific. !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! cosz = max(0.001_r8, bc_in(s)%coszen_pa(ifp)) !copied from previous radiation code... - do ft = 1,numpft_ed + do ft = 1,numpft sb = (90._r8 - (acos(cosz)*180/pi)) * (pi / 180._r8) chil(ifp) = xl(ft) !min(max(xl(ft), -0.4_r8), 0.6_r8 ) if (abs(chil(ifp)) <= 0.01_r8) then @@ -255,7 +259,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) weighted_fsun(L) = 0._r8 weighted_dif_ratio(L,1:hlm_numSWb) = 0._r8 !Each canopy layer (canopy, understorey) has multiple 'parallel' pft's - do ft =1,numpft_ed + do ft =1,numpft if (currentPatch%present(L,ft) == 1)then !only do calculation if there are the appropriate leaves. !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! Diffuse transmittance, tr_dif, do each layer with thickness elai_z. @@ -387,7 +391,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) end do !L do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up. - do ft = 1,numpft_ed + do ft = 1,numpft if (currentPatch%present(L,ft) == 1)then !==============================================================================! ! Iterative solution do scattering @@ -445,7 +449,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) Dif_up(:,:,:) = 0.00_r8 do L = 1, currentPatch%NCL_p !work down from the top of the canopy. weighted_dif_down(L) = 0._r8 - do ft = 1, numpft_ed + do ft = 1, numpft if (currentPatch%present(L,ft) == 1)then !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! First estimates do downward and upward diffuse flux @@ -501,7 +505,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom. weighted_dif_up(L) = 0._r8 - do ft = 1, numpft_ed + do ft = 1, numpft if (currentPatch%present(L,ft) == 1)then !Bounce diffuse radiation off soil surface. iv = currentPatch%nrad(L,ft) + 1 @@ -535,11 +539,11 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? !Add on the radiation coming up through the canopy gaps. !diffuse to diffuse - weighted_dif_up(L) = weighted_dif_up(L) +(1.0-sum(ftweight(L,:,1))) * & + weighted_dif_up(L) = weighted_dif_up(L) +(1.0-sum(ftweight(L,1:numpft,1))) * & weighted_dif_down(L-1) * bc_in(s)%albgr_dif_rb(ib) !direct to diffuse weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(ifp,ib) * & - weighted_dir_tr(L-1) * (1.0-sum(ftweight(L,:,1)))*bc_in(s)%albgr_dir_rb(ib) + weighted_dir_tr(L-1) * (1.0-sum(ftweight(L,1:numpft,1)))*bc_in(s)%albgr_dir_rb(ib) endif end do !L !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! @@ -557,7 +561,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) irep = 0 do L = 1,currentPatch%NCL_p !working from the top down weighted_dif_down(L) = 0._r8 - do ft =1,numpft_ed + do ft =1,numpft if (currentPatch%present(L,ft) == 1)then ! forward diffuse flux within the canopy and at soil, working forward through canopy ! with Dif_up -from previous iteration-. Dif_dn(1) is the forward diffuse flux onto the canopy. @@ -606,13 +610,14 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) endif !present end do!ft if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? - weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1)*(1.0-sum(ftweight(L,:,1))) + weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1) * & + (1.0-sum(ftweight(L,1:numpft,1))) end if end do ! do L loop do L = 1, currentPatch%NCL_p ! working from the top down. weighted_dif_up(L) = 0._r8 - do ft =1,numpft_ed + do ft =1,numpft if (currentPatch%present(L,ft) == 1)then ! Upward diffuse flux at soil or from lower canopy (forward diffuse and unscattered direct beam) iv = currentPatch%nrad(L,ft) + 1 @@ -650,10 +655,10 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? !Add on the radiation coming up through the canopy gaps. - weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,:,1))) * & + weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * & weighted_dif_down(L-1) * bc_in(s)%albgr_dif_rb(ib) weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(ifp,ib) * & - weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,:,1)))*bc_in(s)%albgr_dir_rb(ib) + weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))*bc_in(s)%albgr_dir_rb(ib) end if end do!L end do ! do while over iter @@ -664,7 +669,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) do L = 1, currentPatch%NCL_p !working from the top down. abs_dir_z(:,:) = 0._r8 abs_dif_z(:,:) = 0._r8 - do ft =1,numpft_ed + do ft =1,numpft if (currentPatch%present(L,ft) == 1)then !==============================================================================! ! Compute absorbed flux densities @@ -758,11 +763,11 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) !radiation absorbed from fluxes through unfilled part of lower canopy. if (currentPatch%NCL_p > 1.and.L == currentPatch%NCL_p)then abs_rad(ib) = abs_rad(ib) + weighted_dif_down(L-1) * & - (1.0_r8-sum(ftweight(L,:,1)))*(1.0_r8-bc_in(s)%albgr_dif_rb(ib)) + (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-bc_in(s)%albgr_dif_rb(ib)) abs_rad(ib) = abs_rad(ib) + forc_dir(ifp,ib) * weighted_dir_tr(L-1) * & - (1.0_r8-sum(ftweight(L,:,1)))*(1.0_r8-bc_in(s)%albgr_dir_rb(ib)) - tr_soili = tr_soili + weighted_dif_down(L-1) * (1.0_r8-sum(ftweight(L,:,1))) - tr_soild = tr_soild + forc_dir(ifp,ib) * weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,:,1))) + (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-bc_in(s)%albgr_dir_rb(ib)) + tr_soili = tr_soili + weighted_dif_down(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) + tr_soild = tr_soild + forc_dir(ifp,ib) * weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) endif if (radtype == 1)then @@ -792,7 +797,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) if ( abs(error) > 0.0001)then write(fates_log(),*)'dir ground absorption error',ifp,s,error,currentPatch%sabs_dir(ib), & currentPatch%tr_soil_dir(ib)* & - (1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%NCL_p,ib,sum(ftweight(1,:,1)) + (1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) write(fates_log(),*) 'albedos',currentPatch%sabs_dir(ib) ,currentPatch%tr_soil_dir(ib), & (1.0_r8-bc_in(s)%albgr_dir_rb(ib)),currentPatch%lai @@ -807,7 +812,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) (1.0_r8-bc_in(s)%albgr_dif_rb(ib)))) > 0.0001)then write(fates_log(),*)'dif ground absorption error',ifp,s,currentPatch%sabs_dif(ib) , & (currentPatch%tr_soil_dif(ib)* & - (1.0_r8-bc_in(s)%albgr_dif_rb(ib))),currentPatch%NCL_p,ib,sum(ftweight(1,:,1)) + (1.0_r8-bc_in(s)%albgr_dif_rb(ib))),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) endif endif @@ -820,7 +825,7 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) endif lai_reduction(:) = 0.0_r8 do L = 1, currentPatch%NCL_p - do ft =1,numpft_ed + do ft =1,numpft if (currentPatch%present(L,ft) == 1)then do iv = 1, currentPatch%nrad(L,ft) if (lai_change(L,ft,iv) > 0.0_r8)then @@ -831,25 +836,6 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) enddo enddo - if (lai_change(1,2,1).gt.0.0.and.lai_change(1,2,2).gt.0.0)then - ! write(fates_log(),*) 'lai_change(1,2,12)',lai_change(1,2,1:4) - endif - if (lai_change(1,2,2).gt.0.0.and.lai_change(1,2,3).gt.0.0)then - ! write(fates_log(),*) ' lai_change (1,2,23)',lai_change(1,2,1:4) - endif - if (lai_change(1,1,3).gt.0.0.and.lai_change(1,1,2).gt.0.0)then - ! NO-OP - ! write(fates_log(),*) 'first layer of lai_change 2 3',lai_change(1,1,1:3) - endif - if (lai_change(1,1,3).gt.0.0.and.lai_change(1,1,4).gt.0.0)then - ! NO-OP - ! write(fates_log(),*) 'first layer of lai_change 3 4',lai_change(1,1,1:4) - endif - if (lai_change(1,1,4).gt.0.0.and.lai_change(1,1,5).gt.0.0)then - ! NO-OP - ! write(fates_log(),*) 'first layer of lai_change 4 5',lai_change(1,1,1:5) - endif - if (radtype == 1)then !here we are adding a within-ED radiation scheme tolerance, and then adding the diffrence onto the albedo !it is important that the lower boundary for this is ~1000 times smaller than the tolerance in surface albedo. @@ -866,10 +852,10 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) write(fates_log(),*) 'Large Dir Radn consvn error',error ,ifp,ib write(fates_log(),*) 'diags', bc_out(s)%albd_parb(ifp,ib), bc_out(s)%ftdd_parb(ifp,ib), & bc_out(s)%ftid_parb(ifp,ib), bc_out(s)%fabd_parb(ifp,ib) - write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft_ed,1:4) - write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) - write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) - write(fates_log(),*) 'ftweight',ftweight(1,1:numpft_ed,1:4) + write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'ftweight',ftweight(1,1:numpft,1:diag_nlevleaf) write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno write(fates_log(),*) 'bc_in(s)%albgr_dir_rb(ib)',bc_in(s)%albgr_dir_rb(ib) @@ -885,16 +871,16 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) write(fates_log(),*) '>5% Dif Radn consvn error',error ,ifp,ib write(fates_log(),*) 'diags', bc_out(s)%albi_parb(ifp,ib), bc_out(s)%ftii_parb(ifp,ib), & bc_out(s)%fabi_parb(ifp,ib) - write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft_ed,1:4) - write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) - write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft_ed,1:4) - write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:numpft_ed,1:4) + write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno write(fates_log(),*) 'bc_in(s)%albgr_dif_rb(ib)',bc_in(s)%albgr_dif_rb(ib) - write(fates_log(),*) 'rhol',rhol(1:numpft_ed,:) - write(fates_log(),*) 'ftw',sum(ftweight(1,:,1)),ftweight(1,1:numpft_ed,1) - write(fates_log(),*) 'present',currentPatch%present(1,1:numpft_ed) - write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft_ed,1) + write(fates_log(),*) 'rhol',rhol(1:numpft,:) + write(fates_log(),*) 'ftw',sum(ftweight(1,1:numpft,1)),ftweight(1,1:numpft,1) + write(fates_log(),*) 'present',currentPatch%present(1,1:numpft) + write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft,1) bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error end if @@ -961,7 +947,7 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out) ifp=ifp+1 - if( DEBUG ) write(fates_log(),*) 'edsurfRad_5600',ifp,s,cpatch%NCL_p,numpft_ed + if( DEBUG ) write(fates_log(),*) 'edsurfRad_5600',ifp,s,cpatch%NCL_p,numpft ! zero out various datas cpatch%ed_parsun_z(:,:,:) = 0._r8 @@ -982,7 +968,7 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out) ! cpatch%f_sun is calculated in the surface_albedo routine... do CL = 1, cpatch%NCL_p - do FT = 1,numpft_ed + do FT = 1,numpft if( DEBUG ) write(fates_log(),*) 'edsurfRad_5601',CL,FT,cpatch%nrad(CL,ft) @@ -1029,10 +1015,10 @@ subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out) ! If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo ! are canopy integrated so that layer values equal big leaf values. - if ( DEBUG ) write(fates_log(),*) 'edsurfRad 645 ',cpatch%NCL_p,numpft_ed + if ( DEBUG ) write(fates_log(),*) 'edsurfRad 645 ',cpatch%NCL_p,numpft do CL = 1, cpatch%NCL_p - do FT = 1,numpft_ed + do FT = 1,numpft if ( DEBUG ) write(fates_log(),*) 'edsurfRad 649 ',cpatch%nrad(CL,ft) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index dbc863ced8..e68291057a 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -24,8 +24,6 @@ module FatesPlantHydraulicsMod use FatesInterfaceMod , only : bc_out_type use FatesInterfaceMod , only : hlm_numlevsoil - use EDEcophysconType, only : EDecophyscon - use FatesHydraulicsMemMod, only: ed_site_hydr_type use FatesHydraulicsMemMod, only: ed_patch_hydr_type use FatesHydraulicsMemMod, only: ed_cohort_hydr_type @@ -42,6 +40,8 @@ module FatesPlantHydraulicsMod use FatesHydraulicsMemMod, only: porous_media use FatesHydraulicsMemMod, only: nlevsoi_hyd + use EDPftvarcon, only : EDPftvarcon_inst + ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : isnan => shr_infnan_isnan @@ -125,7 +125,6 @@ subroutine initTreeHydStates(cc_p, bc_in) ! !DESCRIPTION: ! ! !USES: - use EDEcophysConType , only : EDecophyscon ! !ARGUMENTS: type(ed_cohort_type), intent(inout), target :: cc_p ! current cohort pointer @@ -146,7 +145,6 @@ subroutine updateSizeDepTreeHydProps(cc_p,bc_in) ! ! !USES: use FatesConstantsMod , only : pi_const - use EDEcophysConType , only : EDecophyscon use shr_sys_mod , only : shr_sys_abort ! ! !ARGUMENTS: diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 4901ac30a4..0cfd778323 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -25,7 +25,8 @@ module FATESPlantRespPhotosynthMod use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : itrue use FatesInterfaceMod, only : hlm_use_planthydro - use EDTypesMod, only : numpft_ed + use FatesInterfaceMod, only : numpft + use EDTypesMod, only : maxpft use EDTypesMod, only : nlevleaf use EDTypesMod, only : nclmax @@ -272,7 +273,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! but not environmentally dependent ! ------------------------------------------------------------------------ - do ft = 1,numpft_ed + do ft = 1,numpft ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al @@ -308,7 +309,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! NOTE: Only need to flush mask on the number of used pfts, not the whole ! scratch space. ! ------------------------------------------------------------------------ - rate_mask_z(:,1:numpft_ed,:) = .false. + rate_mask_z(:,1:numpft,:) = .false. if(currentPatch%countcohorts > 0.0)then ! Ignore empty patches @@ -1343,7 +1344,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! Now loop through and identify which layer and pft combo has scattering elements do cl = 1,nclmax - do ft = 1,numpft_ed + do ft = 1,numpft currentPatch%present(cl,ft) = 0 do iv = 1, currentPatch%nrad(cl,ft); if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index fe51009cd1..b47740d1f1 100755 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -14,7 +14,6 @@ module SFMainMod use FatesInterfaceMod , only : bc_in_type use EDPftvarcon , only : EDPftvarcon_inst - use EDEcophysconType , only : EDecophyscon use EDtypesMod , only : ed_site_type use EDtypesMod , only : ed_patch_type @@ -882,18 +881,18 @@ subroutine crown_damage ( currentSite ) if (EDPftvarcon_inst%woody(currentCohort%pft) == 1) then !trees only ! Flames lower than bottom of canopy. ! c%hite is height of cohort - if (currentPatch%SH < (currentCohort%hite-currentCohort%hite*EDecophyscon%crown(currentCohort%pft))) then + if (currentPatch%SH < (currentCohort%hite-currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft))) then currentCohort%cfa = 0.0_r8 else ! Flames part of way up canopy. ! Equation 17 in Thonicke et al. 2010. ! flames over bottom of canopy but not over top. if ((currentCohort%hite > 0.0_r8).and.(currentPatch%SH >= & - (currentCohort%hite-currentCohort%hite*EDecophyscon%crown(currentCohort%pft)))) then + (currentCohort%hite-currentCohort%hite*EDPftvarcon_inst%crown(currentCohort%pft)))) then currentCohort%cfa = (currentPatch%SH-currentCohort%hite*(1- & - EDecophyscon%crown(currentCohort%pft)))/(currentCohort%hite* & - EDecophyscon%crown(currentCohort%pft)) + EDPftvarcon_inst%crown(currentCohort%pft)))/(currentCohort%hite* & + EDPftvarcon_inst%crown(currentCohort%pft)) else ! Flames over top of canopy. @@ -942,7 +941,7 @@ subroutine cambial_damage_kill ( currentSite ) do while(associated(currentCohort)) if (EDPftvarcon_inst%woody(currentCohort%pft) == 1) then !trees only ! Equation 21 in Thonicke et al 2010 - bt = EDecophyscon%bark_scaler(currentCohort%pft)*currentCohort%dbh ! bark thickness. + bt = EDPftvarcon_inst%bark_scaler(currentCohort%pft)*currentCohort%dbh ! bark thickness. ! Equation 20 in Thonicke et al. 2010. tau_c = 2.9_r8*bt**2.0_r8 !calculate time it takes to kill cambium (min) ! Equation 19 in Thonicke et al. 2010 @@ -994,7 +993,7 @@ subroutine post_fire_mortality ( currentSite ) currentCohort%crownfire_mort = 0.0_r8 if (EDPftvarcon_inst%woody(currentCohort%pft) == 1) then ! Equation 22 in Thonicke et al. 2010. - currentCohort%crownfire_mort = EDecophyscon%crown_kill(currentCohort%pft)*currentCohort%cfa**3.0_r8 + currentCohort%crownfire_mort = EDPftvarcon_inst%crown_kill(currentCohort%pft)*currentCohort%cfa**3.0_r8 ! Equation 18 in Thonicke et al. 2010. currentCohort%fire_mort = currentCohort%crownfire_mort+currentCohort%cambial_mort- & (currentCohort%crownfire_mort*currentCohort%cambial_mort) !joint prob. diff --git a/main/EDEcophysConType.F90 b/main/EDEcophysConType.F90 deleted file mode 100644 index 8c405ed987..0000000000 --- a/main/EDEcophysConType.F90 +++ /dev/null @@ -1,236 +0,0 @@ -module EDEcophysConType - - !---------------------------------------------------- - ! ED ecophysiological constants - !---------------------------------------------------- - ! - ! !USES: - use shr_kind_mod , only : r8 => shr_kind_r8 - use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - use shr_log_mod , only : errMsg => shr_log_errMsg - use FatesGlobals, only : endrun => fates_endrun - use FatesGlobals, only : fates_log - - use FatesHydraulicsMemMod , only : n_porous_media - use FatesHydraulicsMemMod , only : porous_media - use FatesHydraulicsMemMod , only : npool_tot - use FatesHydraulicsMemMod , only : npool_leaf - use FatesHydraulicsMemMod , only : npool_stem - use FatesHydraulicsMemMod , only : npool_aroot - use FatesHydraulicsMemMod , only : npool_troot - use FatesConstantsMod, only : itrue,ifalse - use FatesInterfaceMod, only : hlm_use_planthydro - - ! - implicit none - save - private - - character(len=*), parameter, private :: sourcefile = & - __FILE__ - - ! - ! !PUBLIC MEMBER FUNCTIONS: - public :: EDecophysconInit - ! - ! !PUBLIC TYPES: - type, public :: EDecophyscon_type - real(r8), pointer :: max_dbh (:) ! maximum dbh at which height growth ceases... - real(r8), pointer :: freezetol (:) ! minimum temperature tolerance... - real(r8), pointer :: wood_density (:) ! wood density g cm^-3 ... - real(r8), pointer :: hgt_min (:) ! sapling height m - real(r8), pointer :: cushion (:) ! labile carbon storage target as multiple of leaf pool. - real(r8), pointer :: leaf_stor_priority (:) ! leaf turnover vs labile carbon use prioritisation. ! (1=lose leaves, 0=use store). - real(r8), pointer :: crown (:) ! fraction of the height of the plant that is occupied by crown. For fire model. - real(r8), pointer :: bark_scaler (:) ! scaler from dbh to bark thickness. For fire model. - real(r8), pointer :: crown_kill (:) ! scaler on fire death. For fire model. - real(r8), pointer :: initd (:) ! initial seedling density - real(r8), pointer :: seed_rain (:) ! seeds that come from outside the gridbox. - real(r8), pointer :: BB_slope (:) ! ball berry slope parameter - real(r8), pointer :: root_long (:) ! root longevity (yrs) - real(r8), pointer :: clone_alloc (:) ! fraction of carbon balance allocated to clonal reproduction. - real(r8), pointer :: seed_alloc (:) ! fraction of carbon balance allocated to seeds. - - ! pft parameters for plant hydraulics (PFT) - real(r8), pointer :: wd (:) ! wood density (distinct from wood_density for testing) [g m-3] - real(r8), pointer :: lma (:) ! leaf mass per area [g m-2] - ! ~ 90 for tropical angiosperms, cf Patino et al. 2012 - ! (existing param 'slatop' is biased high) - real(r8), pointer :: n (:) ! leaf nitrogen [mg g-1] - real(r8), pointer :: p (:) ! leaf phosphorus [mg g-1] - real(r8), pointer :: ldmc (:) ! leaf dry matter content [g g-1] - real(r8), pointer :: lmv (:) ! leaf mass per volume [g m-3] - real(r8), pointer :: psi0 (:) ! sapwood water potential at saturation [MPa] - real(r8), pointer :: psicap (:) ! sapwood water potential at rwcft [MPa] - ! BOC...rhoc, rint_petiole, rint_jansenchoat, ccontent and (maybe) - ! rs2 and rfrac_stem should really be global constants, not pft parameters - real(r8), pointer :: rhoc (:) ! dry matter (or cell wall) density of wood [g cm-3] Siau 1984 - real(r8), pointer :: rint_petiole (:) ! radius of xylem conduits in petioles [um] - real(r8), pointer :: rint_jansenchoat (:) ! average radius of xylem conduits where ks mmts were made [um] - ! taken from choat & jansen XFT database for tropical angiosperms only - real(r8), pointer :: Amaxh (:) ! light-saturated photosynthesis rate [umol m-2 s-1] - real(r8), pointer :: rs2 (:) ! mean absorbing fine root radius [m] ~ 0.001 m? - real(r8), pointer :: srl (:) ! specific root length [m kg-1] - ! ~ 15000 for tropical angiosperms, cf Metcalfe et al. 2008 Plant Soil Fig. 2b; - ! rootdens = 500 kg m-3 is biased high by an order of magnitude - ! (cf Comas et al. 2002 Oecologia); SPA rootdens implies a SRL of only 637 m kg-1. - real(r8), pointer :: ccontent (:) ! carbon content (fraction of dry mass) [-] - ! ~ 0.47 for tropical angiosperms, cf Thomas & Martin (2012) Forests - real(r8), pointer :: latosa (:) ! leaf to sapwood area ratio [m2 m-2] - ! ~ 8000 for tropical angiosperms, cf Patino et al. 2012 - real(r8), pointer :: rfrac_stem (:) ! fraction of total tree resistance (under well-watered conditions) - ! from troot to canopy (i.e., aboveground) [-] ~ 0.625 for tropical angiosperms, - ! cf BC re-analysis of Fisher et al. 2006 - real(r8), pointer :: rootshoot (:) ! root:shoot ratio (belowground-to-aboveground biomass) [-] - ! ~ 0.20 for tropical forests (see Houghton et al. 2001 Table 3, - ! Cairns et al. 1997 Table 2, Jackson et al. 1996 Table 3) - real(r8), pointer :: avuln_gs (:) ! stomata PLC: vulnerability curve shape parameter [-] - real(r8), pointer :: p50_gs (:) ! stomata PLC: water potential at 50% loss of conductivity [Pa] - - ! pft parameters for plant hydraulics (PFT x tissue type (leaf, stem, troot, aroot)) - real(r8), pointer :: kmax_node (:,:) ! xylem PLC: maximum xylem hydraulic conductivity [kg m-1 s-1 Pa-1] - real(r8), pointer :: avuln_node (:,:) ! xylem PLC: vulnerability curve shape parameter [-] - real(r8), pointer :: p50_node (:,:) ! xylem PLC: water potential at 50% loss of conductivity [Pa] - real(r8), pointer :: thetas_node (:,:) ! P-V curve: saturated volumetric water content for node [m3 m-3] - real(r8), pointer :: epsil_node (:,:) ! P-V curve: bulk elastic modulus [MPa] - real(r8), pointer :: pinot_node (:,:) ! P-V curve: osmotic potential at full turgor [MPa] - real(r8), pointer :: pitlp_node (:,:) ! P-V curve: osmotic potential at turgor loss [MPa] - real(r8), pointer :: resid_node (:,:) ! P-V curve: residual fraction [-] - real(r8), pointer :: rwctlp_node (:,:) ! P-V curve: total relative water content at turgor loss [g or m3 H2O / g or m3 H2O, sat] - real(r8), pointer :: fcap_node (:,:) ! P-V curve: fraction of (1-resid_node) that is capillary in source [-] - real(r8), pointer :: rwcft_node (:,:) ! P-V curve: total RWC @ which elastic drainage begins [-] - real(r8), pointer :: rwccap_node (:,:) ! P-V curve: total RWC @ which capillary reserves exhausted - real(r8), pointer :: slp_node (:,:) ! P-V curve: slope of capillary region of curve (sapwood only) - real(r8), pointer :: intercept_node (:,:) ! P-V curve: intercept of capillary region of curve (sapwood only) - real(r8), pointer :: corrInt_node (:,:) ! P-V curve: correction for nonzero psi0 - - - - end type EDecophyscon_type - - type(EDecophyscon_type), public :: EDecophyscon ! ED ecophysiological constants structure - !------------------------------------------------------------------------ - - - - - -contains - - !------------------------------------------------------------------------ - subroutine EDecophysconInit(EDpftvarcon_inst, numpft) - ! - ! !USES: - use EDPftvarcon, only : EDPftvarcon_type - ! - ! !ARGUMENTS: - type(EDpftVarCon_type) , intent(in) :: EDpftvarcon_inst - integer , intent(in) :: numpft - ! - ! !LOCAL VARIABLES: - integer :: m, ib, n, k - integer :: lb - !------------------------------------------------------------------------ - - - if( lbound(EDPftvarcon_inst%max_dbh,dim=1) .eq. 0) then - lb = 0 - else - lb = 1 - end if - - allocate( EDecophyscon%max_dbh (0:numpft)); EDecophyscon%max_dbh (:) = nan - allocate( EDecophyscon%freezetol (0:numpft)); EDecophyscon%freezetol (:) = nan - allocate( EDecophyscon%wood_density (0:numpft)); EDecophyscon%wood_density (:) = nan - allocate( EDecophyscon%hgt_min (0:numpft)); EDecophyscon%hgt_min (:) = nan - allocate( EDecophyscon%cushion (0:numpft)); EDecophyscon%cushion (:) = nan - allocate( EDecophyscon%leaf_stor_priority (0:numpft)); EDecophyscon%leaf_stor_priority (:) = nan - allocate( EDecophyscon%crown (0:numpft)); EDecophyscon%crown (:) = nan - allocate( EDecophyscon%bark_scaler (0:numpft)); EDecophyscon%bark_scaler (:) = nan - allocate( EDecophyscon%crown_kill (0:numpft)); EDecophyscon%crown_kill (:) = nan - allocate( EDecophyscon%initd (0:numpft)); EDecophyscon%initd (:) = nan - allocate( EDecophyscon%seed_rain (0:numpft)); EDecophyscon%seed_rain (:) = nan - allocate( EDecophyscon%BB_slope (0:numpft)); EDecophyscon%BB_slope (:) = nan - allocate( EDecophyscon%root_long (0:numpft)); EDecophyscon%root_long (:) = nan - allocate( EDecophyscon%seed_alloc (0:numpft)); EDecophyscon%seed_alloc (:) = nan - allocate( EDecophyscon%clone_alloc (0:numpft)); EDecophyscon%clone_alloc (:) = nan - - do m = lb,numpft - EDecophyscon%max_dbh(m) = EDPftvarcon_inst%max_dbh(m) - EDecophyscon%freezetol(m) = EDPftvarcon_inst%freezetol(m) - EDecophyscon%wood_density(m) = EDPftvarcon_inst%wood_density(m) - EDecophyscon%hgt_min(m) = EDPftvarcon_inst%hgt_min(m) - EDecophyscon%cushion(m) = EDPftvarcon_inst%cushion(m) - EDecophyscon%leaf_stor_priority(m) = EDPftvarcon_inst%leaf_stor_priority(m) - EDecophyscon%crown(m) = EDPftvarcon_inst%crown(m) - EDecophyscon%bark_scaler(m) = EDPftvarcon_inst%bark_scaler(m) - EDecophyscon%crown_kill(m) = EDPftvarcon_inst%crown_kill(m) - EDecophyscon%initd(m) = EDPftvarcon_inst%initd(m) - EDecophyscon%seed_rain(m) = EDPftvarcon_inst%seed_rain(m) - EDecophyscon%bb_slope(m) = EDPftvarcon_inst%bb_slope(m) - EDecophyscon%root_long(m) = EDPftvarcon_inst%root_long(m) - EDecophyscon%seed_alloc(m) = EDPftvarcon_inst%seed_alloc(m) - EDecophyscon%clone_alloc(m) = EDPftvarcon_inst%clone_alloc(m) - end do - - - if (hlm_use_planthydro.eq.itrue) then - allocate( EDecophyscon%wd (0:numpft) ); EDecophyscon%wd (:) = nan - allocate( EDecophyscon%lma (0:numpft) ); EDecophyscon%lma (:) = nan - allocate( EDecophyscon%n (0:numpft) ); EDecophyscon%n (:) = nan - allocate( EDecophyscon%p (0:numpft) ); EDecophyscon%p (:) = nan - allocate( EDecophyscon%ldmc (0:numpft) ); EDecophyscon%ldmc (:) = nan - allocate( EDecophyscon%lmv (0:numpft) ); EDecophyscon%lmv (:) = nan - allocate( EDecophyscon%psi0 (0:numpft) ); EDecophyscon%psi0 (:) = nan - allocate( EDecophyscon%psicap (0:numpft) ); EDecophyscon%psicap (:) = nan - allocate( EDecophyscon%rhoc (0:numpft) ); EDecophyscon%rhoc (:) = nan - allocate( EDecophyscon%rint_petiole (0:numpft) ); EDecophyscon%rint_petiole (:) = nan - allocate( EDecophyscon%rint_jansenchoat (0:numpft) ); EDecophyscon%rint_jansenchoat (:) = nan - allocate( EDecophyscon%Amaxh (0:numpft) ); EDecophyscon%Amaxh (:) = nan - allocate( EDecophyscon%rs2 (0:numpft) ); EDecophyscon%rs2 (:) = nan - allocate( EDecophyscon%srl (0:numpft) ); EDecophyscon%srl (:) = nan - allocate( EDecophyscon%ccontent (0:numpft) ); EDecophyscon%ccontent (:) = nan - allocate( EDecophyscon%latosa (0:numpft) ); EDecophyscon%latosa (:) = nan - allocate( EDecophyscon%rfrac_stem (0:numpft) ); EDecophyscon%rfrac_stem (:) = nan - allocate( EDecophyscon%rootshoot (0:numpft) ); EDecophyscon%rootshoot (:) = nan - allocate( EDecophyscon%avuln_gs (0:numpft) ); EDecophyscon%avuln_gs (:) = nan - allocate( EDecophyscon%p50_gs (0:numpft) ); EDecophyscon%p50_gs (:) = nan - - allocate( EDecophyscon%kmax_node (0:numpft,1:n_porous_media) ); EDecophyscon%kmax_node (:,:) = nan - allocate( EDecophyscon%avuln_node (0:numpft,1:n_porous_media) ); EDecophyscon%avuln_node (:,:) = nan - allocate( EDecophyscon%p50_node (0:numpft,1:n_porous_media) ); EDecophyscon%p50_node (:,:) = nan - allocate( EDecophyscon%thetas_node (0:numpft,1:n_porous_media) ); EDecophyscon%thetas_node (:,:) = nan - allocate( EDecophyscon%epsil_node (0:numpft,1:n_porous_media) ); EDecophyscon%epsil_node (:,:) = nan - allocate( EDecophyscon%pinot_node (0:numpft,1:n_porous_media) ); EDecophyscon%pinot_node (:,:) = nan - allocate( EDecophyscon%pitlp_node (0:numpft,1:n_porous_media) ); EDecophyscon%pitlp_node (:,:) = nan - allocate( EDecophyscon%resid_node (0:numpft,1:n_porous_media) ); EDecophyscon%resid_node (:,:) = nan - allocate( EDecophyscon%rwctlp_node (0:numpft,1:n_porous_media) ); EDecophyscon%rwctlp_node (:,:) = nan - allocate( EDecophyscon%fcap_node (0:numpft,1:n_porous_media) ); EDecophyscon%fcap_node (:,:) = nan - allocate( EDecophyscon%rwcft_node (0:numpft,1:n_porous_media) ); EDecophyscon%rwcft_node (:,:) = nan - allocate( EDecophyscon%rwccap_node (0:numpft,1:n_porous_media) ); EDecophyscon%rwccap_node (:,:) = nan - allocate( EDecophyscon%slp_node (0:numpft,1:n_porous_media) ); EDecophyscon%slp_node (:,:) = nan - allocate( EDecophyscon%intercept_node (0:numpft,1:n_porous_media) ); EDecophyscon%intercept_node (:,:) = nan - allocate( EDecophyscon%corrInt_node (0:numpft,1:n_porous_media) ); EDecophyscon%corrInt_node (:,:) = nan - - ! ------------------------------------------------------------------------------------------------ - ! Until the hydraulics parameter are added to the parameter file, they need a location to be set. - ! This happens here until further notice. - ! ------------------------------------------------------------------------------------------------ - call SetHydraulicsTestingParams(EDecophyscon) - - end if - - end subroutine EDecophysconInit - - subroutine SetHydraulicsTestingParams(EDEcophyscon) - - ! Arguments - type(EDecophyscon_type), intent(inout) :: EDEcophyscon - - write(fates_log(),*) 'FATES Plant Hydraulics is still under development, ending run.' - call endrun(msg=errMsg(sourcefile, __LINE__)) - - - end subroutine SetHydraulicsTestingParams - -end module EDEcophysConType diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 811151023f..73bda67f5d 100755 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -12,17 +12,17 @@ module EDInitMod use FatesGlobals , only : fates_log use FatesInterfaceMod , only : hlm_is_restart use EDPftvarcon , only : EDPftvarcon_inst - use EDEcophysConType , only : EDecophyscon use EDGrowthFunctionsMod , only : bdead, bleaf, dbh use EDCohortDynamicsMod , only : create_cohort, fuse_cohorts, sort_cohorts use EDPatchDynamicsMod , only : create_patch use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, area use EDTypesMod , only : ncwd use EDTypesMod , only : nuMWaterMem - use EDTypesMod , only : numpft_ed + use EDTypesMod , only : maxpft use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_use_inventory_init + use FatesInterfaceMod , only : numpft ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -216,8 +216,8 @@ subroutine init_patches( nsites, sites, bc_in) real(r8) :: cwd_ag_local(ncwd) real(r8) :: cwd_bg_local(ncwd) real(r8) :: spread_local(nclmax) - real(r8) :: leaf_litter_local(numpft_ed) - real(r8) :: root_litter_local(numpft_ed) + real(r8) :: leaf_litter_local(maxpft) + real(r8) :: root_litter_local(maxpft) real(r8) :: age !notional age of this patch type(ed_patch_type), pointer :: newp @@ -302,17 +302,17 @@ subroutine init_cohorts( patch_in, bc_in) patch_in%tallest => null() patch_in%shortest => null() - do pft = 1,numpft_ed !FIX(RF,032414) - turning off veg dynamics + do pft = 1,numpft - if(EDecophyscon%initd(pft)>1.0E-7) then + if(EDPftvarcon_inst%initd(pft)>1.0E-7) then allocate(temp_cohort) ! temporary cohort temp_cohort%pft = pft - temp_cohort%n = EDecophyscon%initd(pft) * patch_in%area - temp_cohort%hite = EDecophyscon%hgt_min(pft) - !temp_cohort%n = 0.5_r8 * 0.0028_r8 * patch_in%area ! BOC for fixed size runs EDecophyscon%initd(pft) * patch_in%area - !temp_cohort%hite = 28.65_r8 ! BOC translates to DBH of 50cm. EDecophyscon%hgt_min(pft) + temp_cohort%n = EDPftvarcon_inst%initd(pft) * patch_in%area + temp_cohort%hite = EDPftvarcon_inst%hgt_min(pft) + !temp_cohort%n = 0.5_r8 * 0.0028_r8 * patch_in%area ! BOC for fixed size runs EDPftvarcon_inst%initd(pft) * patch_in%area + !temp_cohort%hite = 28.65_r8 ! BOC translates to DBH of 50cm. EDPftvarcon_inst%hgt_min(pft) temp_cohort%dbh = Dbh(temp_cohort) ! FIX(RF, 090314) - comment out addition of ' + 0.0001_r8*pft ' - seperate out PFTs a little bit... temp_cohort%canopy_trim = 1.0_r8 temp_cohort%bdead = Bdead(temp_cohort) @@ -321,13 +321,13 @@ subroutine init_cohorts( patch_in, bc_in) temp_cohort%b = temp_cohort%balive + temp_cohort%bdead if( EDPftvarcon_inst%evergreen(pft) == 1) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(pft) + temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(pft) temp_cohort%laimemory = 0._r8 cstatus = 2 endif if( EDPftvarcon_inst%season_decid(pft) == 1 ) then !for dorment places - temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(pft) !stored carbon in new seedlings. + temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(pft) !stored carbon in new seedlings. if(patch_in%siteptr%status == 2)then temp_cohort%laimemory = 0.0_r8 else @@ -339,7 +339,7 @@ subroutine init_cohorts( patch_in, bc_in) endif if ( EDPftvarcon_inst%stress_decid(pft) == 1 ) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(pft) + temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(pft) temp_cohort%laimemory = Bleaf(temp_cohort) temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory cstatus = patch_in%siteptr%dstatus diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 1db955f715..a2a204454f 100755 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -17,6 +17,7 @@ module EDMainMod use FatesInterfaceMod , only : hlm_use_ed_st3 use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_masterproc + use FatesInterfaceMod , only : numpft use EDCohortDynamicsMod , only : allocate_live_biomass use EDCohortDynamicsMod , only : terminate_cohorts use EDCohortDynamicsMod , only : fuse_cohorts @@ -34,7 +35,6 @@ module EDMainMod use SFMainMod , only : fire_model use EDTypesMod , only : get_age_class_index use EDtypesMod , only : ncwd - use EDTypesMod , only : numpft_ed use EDtypesMod , only : ed_site_type use EDtypesMod , only : ed_patch_type use EDtypesMod , only : ed_cohort_type @@ -228,9 +228,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) small_no = 0.0000000000_r8 ! Obviously, this is arbitrary. RF - changed to zero - do ft = 1,numpft_ed - currentSite%dseed_dt(ft) = 0._r8 ! zero the dseed_dt at the site level before looping through patches and adding the fluxes from each patch - end do + currentSite%dseed_dt(:) = 0._r8 currentSite%seed_rain_flux(:) = 0._r8 currentPatch => currentSite%youngest_patch @@ -309,7 +307,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentPatch%cwd_bg(c) = currentPatch%cwd_bg(c) + currentPatch%dcwd_bg_dt(c)* hlm_freq_day enddo - do ft = 1,numpft_ed + do ft = 1,numpft currentPatch%leaf_litter(ft) = currentPatch%leaf_litter(ft) + currentPatch%dleaf_litter_dt(ft)* hlm_freq_day currentPatch%root_litter(ft) = currentPatch%root_litter(ft) + currentPatch%droot_litter_dt(ft)* hlm_freq_day enddo @@ -325,7 +323,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) endif enddo - do ft = 1,numpft_ed + do ft = 1,numpft if(currentPatch%leaf_litter(ft) fates_endrun diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index cfadc1ee0c..d4b9fdf18a 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -36,12 +36,12 @@ module EDPftvarcon real(r8), allocatable :: cushion (:) ! labile carbon storage target as multiple of leaf pool. real(r8), allocatable :: leaf_stor_priority (:) ! leaf turnover vs labile carbon use prioritisation ! (1 = lose leaves, 0 = use store). - real(r8), allocatable :: crown (:) - real(r8), allocatable :: bark_scaler (:) - real(r8), allocatable :: crown_kill (:) - real(r8), allocatable :: initd (:) - real(r8), allocatable :: seed_rain (:) - real(r8), allocatable :: BB_slope (:) + real(r8), allocatable :: crown (:) ! fraction of the height of the plant that is occupied by crown. For fire model. + real(r8), allocatable :: bark_scaler (:) ! scaler from dbh to bark thickness. For fire model. + real(r8), allocatable :: crown_kill (:) ! scaler on fire death. For fire model. + real(r8), allocatable :: initd (:) ! initial seedling density + real(r8), allocatable :: seed_rain (:) ! seeds that come from outside the gridbox. + real(r8), allocatable :: BB_slope (:) ! ball berry slope parameter real(r8), allocatable :: root_long (:) ! root longevity (yrs) real(r8), allocatable :: clone_alloc (:) ! fraction of carbon balance allocated to clonal reproduction. real(r8), allocatable :: seed_alloc (:) ! fraction of carbon balance allocated to seeds. diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index b60d31108c..e0b7f2b44c 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -11,6 +11,7 @@ module EDTypesMod save integer, parameter :: maxPatchesPerSite = 10 ! maximum number of patches to live on a site + integer, parameter :: maxCohortsPerPatch = 160 ! maximum number of cohorts per patch integer, parameter :: nclmax = 2 ! Maximum number of canopy layers integer, parameter :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter :: ican_ustory = 2 ! Nominal index for understory in two-canopy system @@ -21,9 +22,6 @@ module EDTypesMod ! are used, but this helps allocate scratch ! space and output arrays. - integer, parameter :: numpft_ed = 2 ! number of PFTs used in ED. - - integer, parameter :: maxCohortsPerPatch = nclmax * numpft_ed * nlevleaf ! maximum number of cohorts to live on a patch ! TODO: we use this cp_maxSWb only because we have a static array q(size=2) of ! land-ice abledo for vis and nir. This should be a parameter, which would @@ -125,30 +123,7 @@ module EDTypesMod (/"background","hydraulic ","carbon ","impact ","fire "/) - ! ------------------------------------------------------------------------------------- - ! These vectors are used for history output mapping - ! CLM/ALM have limited support for multi-dimensional history output arrays. - ! FATES structure and composition is multi-dimensional, so we end up "multi-plexing" - ! multiple dimensions into one dimension. These new dimensions need definitions, - ! mapping to component dimensions, and definitions for those component dimensions as - ! well. - ! ------------------------------------------------------------------------------------- - - real(r8) ,allocatable :: fates_hdim_levsclass(:) ! plant size class lower bound dimension - integer , allocatable :: fates_hdim_pfmap_levscpf(:) ! map of pfts into size-class x pft dimension - integer , allocatable :: fates_hdim_scmap_levscpf(:) ! map of size-class into size-class x pft dimension - real(r8), allocatable :: fates_hdim_levage(:) ! patch age lower bound dimension - integer , allocatable :: fates_hdim_levpft(:) ! plant pft dimension - integer , allocatable :: fates_hdim_levfuel(:) ! fire fuel class dimension - integer , allocatable :: fates_hdim_levcwdsc(:) ! cwd class dimension - integer , allocatable :: fates_hdim_levcan(:) ! canopy-layer dimension - integer , allocatable :: fates_hdim_canmap_levcnlf(:) ! canopy-layer map into the canopy-layer x leaf-layer dimension - integer , allocatable :: fates_hdim_lfmap_levcnlf(:) ! leaf-layer map into the canopy-layer x leaf-layer dimension - integer , allocatable :: fates_hdim_canmap_levcnlfpf(:) ! canopy-layer map into the canopy-layer x pft x leaf-layer dimension - integer , allocatable :: fates_hdim_lfmap_levcnlfpf(:) ! leaf-layer map into the canopy-layer x pft x leaf-layer dimension - integer , allocatable :: fates_hdim_pftmap_levcnlfpf(:) ! pft map into the canopy-layer x pft x leaf-layer dimension - integer , allocatable :: fates_hdim_scmap_levscag(:) ! map of size-class into size-class x patch age dimension - integer , allocatable :: fates_hdim_agmap_levscag(:) ! map of patch-age into size-class x patch age dimension + !************************************ !** COHORT type structure ** @@ -324,7 +299,7 @@ module EDTypesMod ! LEAF ORGANIZATION real(r8) :: spread(nclmax) ! dynamic ratio of dbh to canopy area: cm/m2 - real(r8) :: pft_agb_profile(numpft_ed,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 + real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 real(r8) :: canopy_layer_lai(nclmax) ! lai that is shading this canopy layer: m2/m2 real(r8) :: total_canopy_area ! area that is covered by vegetation : m2 real(r8) :: total_tree_area ! area that is covered by woody vegetation : m2 @@ -333,33 +308,33 @@ module EDTypesMod real(r8) :: lai ! leaf area index of patch real(r8) :: zstar ! height of smallest canopy tree -- only meaningful in "strict PPA" mode - real(r8) :: tlai_profile(nclmax,numpft_ed,nlevleaf) ! total leaf area in each canopy layer, pft, and leaf layer. m2/m2 - real(r8) :: elai_profile(nclmax,numpft_ed,nlevleaf) ! exposed leaf area in each canopy layer, pft, and leaf layer. m2/m2 - real(r8) :: tsai_profile(nclmax,numpft_ed,nlevleaf) ! total stem area in each canopy layer, pft, and leaf layer. m2/m2 - real(r8) :: esai_profile(nclmax,numpft_ed,nlevleaf) ! exposed stem area in each canopy layer, pft, and leaf layer. m2/m2 - real(r8) :: layer_height_profile(nclmax,numpft_ed,nlevleaf) - real(r8) :: canopy_area_profile(nclmax,numpft_ed,nlevleaf) ! fraction of canopy in each canopy + real(r8) :: tlai_profile(nclmax,maxpft,nlevleaf) ! total leaf area in each canopy layer, pft, and leaf layer. m2/m2 + real(r8) :: elai_profile(nclmax,maxpft,nlevleaf) ! exposed leaf area in each canopy layer, pft, and leaf layer. m2/m2 + real(r8) :: tsai_profile(nclmax,maxpft,nlevleaf) ! total stem area in each canopy layer, pft, and leaf layer. m2/m2 + real(r8) :: esai_profile(nclmax,maxpft,nlevleaf) ! exposed stem area in each canopy layer, pft, and leaf layer. m2/m2 + real(r8) :: layer_height_profile(nclmax,maxpft,nlevleaf) + real(r8) :: canopy_area_profile(nclmax,maxpft,nlevleaf) ! fraction of canopy in each canopy ! layer, pft, and leaf layer:- - integer :: present(nclmax,numpft_ed) ! is there any of this pft in this canopy layer? - integer :: nrad(nclmax,numpft_ed) ! number of exposed leaf layers for each canopy layer and pft - integer :: ncan(nclmax,numpft_ed) ! number of total leaf layers for each canopy layer and pft + integer :: present(nclmax,maxpft) ! is there any of this pft in this canopy layer? + integer :: nrad(nclmax,maxpft) ! number of exposed leaf layers for each canopy layer and pft + integer :: ncan(nclmax,maxpft) ! number of total leaf layers for each canopy layer and pft !RADIATION FLUXES - real(r8) :: fabd_sun_z(nclmax,numpft_ed,nlevleaf) ! sun fraction of direct light absorbed by each canopy + real(r8) :: fabd_sun_z(nclmax,maxpft,nlevleaf) ! sun fraction of direct light absorbed by each canopy ! layer, pft, and leaf layer:- - real(r8) :: fabd_sha_z(nclmax,numpft_ed,nlevleaf) ! shade fraction of direct light absorbed by each canopy + real(r8) :: fabd_sha_z(nclmax,maxpft,nlevleaf) ! shade fraction of direct light absorbed by each canopy ! layer, pft, and leaf layer:- - real(r8) :: fabi_sun_z(nclmax,numpft_ed,nlevleaf) ! sun fraction of indirect light absorbed by each canopy + real(r8) :: fabi_sun_z(nclmax,maxpft,nlevleaf) ! sun fraction of indirect light absorbed by each canopy ! layer, pft, and leaf layer:- - real(r8) :: fabi_sha_z(nclmax,numpft_ed,nlevleaf) ! shade fraction of indirect light absorbed by each canopy + real(r8) :: fabi_sha_z(nclmax,maxpft,nlevleaf) ! shade fraction of indirect light absorbed by each canopy ! layer, pft, and leaf layer:- - real(r8) :: ed_laisun_z(nclmax,numpft_ed,nlevleaf) ! amount of LAI in the sun in each canopy layer, + real(r8) :: ed_laisun_z(nclmax,maxpft,nlevleaf) ! amount of LAI in the sun in each canopy layer, ! pft, and leaf layer. m2/m2 - real(r8) :: ed_laisha_z(nclmax,numpft_ed,nlevleaf) ! amount of LAI in the shade in each canopy layer, - real(r8) :: ed_parsun_z(nclmax,numpft_ed,nlevleaf) ! PAR absorbed in the sun in each canopy layer, - real(r8) :: ed_parsha_z(nclmax,numpft_ed,nlevleaf) ! PAR absorbed in the shade in each canopy layer, - real(r8) :: f_sun(nclmax,numpft_ed,nlevleaf) ! fraction of leaves in the sun in each canopy layer, pft, + real(r8) :: ed_laisha_z(nclmax,maxpft,nlevleaf) ! amount of LAI in the shade in each canopy layer, + real(r8) :: ed_parsun_z(nclmax,maxpft,nlevleaf) ! PAR absorbed in the sun in each canopy layer, + real(r8) :: ed_parsha_z(nclmax,maxpft,nlevleaf) ! PAR absorbed in the shade in each canopy layer, + real(r8) :: f_sun(nclmax,maxpft,nlevleaf) ! fraction of leaves in the sun in each canopy layer, pft, ! and leaf layer. m2/m2 real(r8),allocatable :: tr_soil_dir(:) ! fraction of incoming direct radiation that (cm_numSWb) @@ -376,20 +351,20 @@ module EDTypesMod !SEED BANK - real(r8) :: seeds_in(numpft_ed) ! seed production KgC/m2/year - real(r8) :: seed_decay(numpft_ed) ! seed decay in KgC/m2/year - real(r8) :: seed_germination(numpft_ed) ! germination rate of seed pool in KgC/m2/year + real(r8) :: seeds_in(maxpft) ! seed production KgC/m2/year + real(r8) :: seed_decay(maxpft) ! seed decay in KgC/m2/year + real(r8) :: seed_germination(maxpft) ! germination rate of seed pool in KgC/m2/year ! PHOTOSYNTHESIS - real(r8) :: psn_z(nclmax,numpft_ed,nlevleaf) ! carbon assimilation in each canopy layer, pft, and leaf layer. umolC/m2/s + real(r8) :: psn_z(nclmax,maxpft,nlevleaf) ! carbon assimilation in each canopy layer, pft, and leaf layer. umolC/m2/s ! real(r8) :: gpp ! total patch gpp: KgC/m2/year ! real(r8) :: npp ! total patch npp: KgC/m2/year ! ROOTS real(r8), allocatable :: rootfr_ft(:,:) ! root fraction of each PFT in each soil layer:- real(r8), allocatable :: rootr_ft(:,:) ! fraction of water taken from each PFT and soil layer:- - real(r8) :: btran_ft(numpft_ed) ! btran calculated seperately for each PFT:- + real(r8) :: btran_ft(maxpft) ! btran calculated seperately for each PFT:- ! DISTURBANCE real(r8) :: disturbance_rates(n_dist_types) ! disturbance rate from 1) mortality and 2) fire: fraction/day @@ -399,8 +374,8 @@ module EDTypesMod ! Pools of litter (non respiring) real(r8) :: cwd_ag(ncwd) ! above ground coarse wood debris litter that does not respire. KgC/m2 real(r8) :: cwd_bg(ncwd) ! below ground coarse wood debris litter that does not respire. KgC/m2 - real(r8) :: leaf_litter(numpft_ed) ! above ground leaf litter that does not respire. KgC/m2 - real(r8) :: root_litter(numpft_ed) ! below ground fine root litter that does not respire. KgC/m2 + real(r8) :: leaf_litter(maxpft) ! above ground leaf litter that does not respire. KgC/m2 + real(r8) :: root_litter(maxpft) ! below ground fine root litter that does not respire. KgC/m2 ! Fluxes of litter (non respiring) real(r8) :: fragmentation_scaler ! Scale rate of litter fragmentation. 0 to 1. @@ -410,18 +385,18 @@ module EDTypesMod real(r8) :: cwd_bg_out(ncwd) ! Flux out of BG CWD into BG litter KgC/m2/ - real(r8) :: leaf_litter_in(numpft_ed) ! Flux in to AG leaf litter from leaf turnover and mortality KgC/m2/y - real(r8) :: leaf_litter_out(numpft_ed) ! Flux out of AG leaf litter from fragmentation KgC/m2/y - real(r8) :: root_litter_in(numpft_ed) ! Flux in to BG root litter from leaf turnover and mortality KgC/m2/y - real(r8) :: root_litter_out(numpft_ed) ! Flux out of BG root from fragmentation KgC/m2/y + real(r8) :: leaf_litter_in(maxpft) ! Flux in to AG leaf litter from leaf turnover and mortality KgC/m2/y + real(r8) :: leaf_litter_out(maxpft) ! Flux out of AG leaf litter from fragmentation KgC/m2/y + real(r8) :: root_litter_in(maxpft) ! Flux in to BG root litter from leaf turnover and mortality KgC/m2/y + real(r8) :: root_litter_out(maxpft) ! Flux out of BG root from fragmentation KgC/m2/y ! Derivatives of litter (non respiring) real(r8) :: dcwd_AG_dt(ncwd) ! rate of change of above ground CWD in each size class: KgC/m2/year. real(r8) :: dcwd_BG_dt(ncwd) ! rate of change of below ground CWD in each size class: KgC/m2/year. - real(r8) :: dleaf_litter_dt(numpft_ed) ! rate of change of leaf litter in each size class: KgC/m2/year. - real(r8) :: droot_litter_dt(numpft_ed) ! rate of change of root litter in each size class: KgC/m2/year. + real(r8) :: dleaf_litter_dt(maxpft) ! rate of change of leaf litter in each size class: KgC/m2/year. + real(r8) :: droot_litter_dt(maxpft) ! rate of change of root litter in each size class: KgC/m2/year. - real(r8) :: repro(numpft_ed) ! allocation to reproduction per PFT : KgC/m2 + real(r8) :: repro(maxpft) ! allocation to reproduction per PFT : KgC/m2 !FUEL CHARECTERISTICS real(r8) :: sum_fuel ! total ground fuel related to ros (omits 1000hr fuels): KgC/m2 @@ -531,9 +506,9 @@ module EDTypesMod real(r8) :: water_memory(numWaterMem) ! last 10 days of soil moisture memory... !SEED BANK - real(r8) :: seed_bank(numpft_ed) ! seed pool in KgC/m2/year - real(r8) :: dseed_dt(numpft_ed) - real(r8) :: seed_rain_flux(numpft_ed) ! flux of seeds from exterior KgC/m2/year (needed for C balance purposes) + real(r8) :: seed_bank(maxpft) ! seed pool in KgC/m2/year + real(r8) :: dseed_dt(maxpft) + real(r8) :: seed_rain_flux(maxpft) ! flux of seeds from exterior KgC/m2/year (needed for C balance purposes) ! FIRE real(r8) :: wind ! daily wind in m/min for Spitfire units @@ -542,7 +517,7 @@ module EDTypesMod real(r8) :: frac_burnt ! fraction of soil burnt in this day. real(r8) :: total_burn_flux_to_atm ! total carbon burnt to the atmosphere in this day. KgC/site real(r8) :: cwd_ag_burned(ncwd) - real(r8) :: leaf_litter_burned(numpft_ed) + real(r8) :: leaf_litter_burned(maxpft) ! PLANT HYDRAULICS type(ed_site_hydr_type), pointer :: si_hydr @@ -565,113 +540,8 @@ module EDTypesMod end type ed_site_type - public :: ed_hist_scpfmaps - contains - - !-------------------------------------------------------------------------------------! - subroutine ed_hist_scpfmaps - ! This subroutine allocates and populates the variables - ! that define the mapping of variables in history files in the "scpf" format - ! back to - ! its respective size-class "sc" and pft "pf" - - integer :: i - integer :: isc - integer :: ipft - integer :: icwd - integer :: ifuel - integer :: ican - integer :: ileaf - integer :: iage - - allocate( fates_hdim_levsclass(1:nlevsclass_ed )) - allocate( fates_hdim_pfmap_levscpf(1:nlevsclass_ed*maxpft)) - allocate( fates_hdim_scmap_levscpf(1:nlevsclass_ed*maxpft)) - allocate( fates_hdim_levpft(1:maxpft )) - allocate( fates_hdim_levfuel(1:NFSC )) - allocate( fates_hdim_levcwdsc(1:NCWD )) - allocate( fates_hdim_levage(1:nlevage_ed )) - - allocate( fates_hdim_levcan(nclmax)) - allocate( fates_hdim_canmap_levcnlf(nlevleaf*nclmax)) - allocate( fates_hdim_lfmap_levcnlf(nlevleaf*nclmax)) - allocate( fates_hdim_canmap_levcnlfpf(nlevleaf*nclmax*numpft_ed)) - allocate( fates_hdim_lfmap_levcnlfpf(nlevleaf*nclmax*numpft_ed)) - allocate( fates_hdim_pftmap_levcnlfpf(nlevleaf*nclmax*numpft_ed)) - allocate( fates_hdim_scmap_levscag(nlevsclass_ed * nlevage_ed )) - allocate( fates_hdim_agmap_levscag(nlevsclass_ed * nlevage_ed )) - - ! Fill the IO array of plant size classes - ! For some reason the history files did not like - ! a hard allocation of sclass_ed - fates_hdim_levsclass(:) = sclass_ed(:) - - fates_hdim_levage(:) = ageclass_ed(:) - - ! make pft array - do ipft=1,maxpft - fates_hdim_levpft(ipft) = ipft - end do - - ! make fuel array - do ifuel=1,NFSC - fates_hdim_levfuel(ifuel) = ifuel - end do - - ! make cwd array - do icwd=1,NCWD - fates_hdim_levcwdsc(icwd) = icwd - end do - - ! make canopy array - do ican = 1,nclmax - fates_hdim_levcan(ican) = ican - end do - - ! Fill the IO arrays that match pft and size class to their combined array - i=0 - do ipft=1,maxpft - do isc=1,nlevsclass_ed - i=i+1 - fates_hdim_pfmap_levscpf(i) = ipft - fates_hdim_scmap_levscpf(i) = isc - end do - end do - - i=0 - do ican=1,nclmax - do ileaf=1,nlevleaf - i=i+1 - fates_hdim_canmap_levcnlf(i) = ican - fates_hdim_lfmap_levcnlf(i) = ileaf - end do - end do - - i=0 - do iage=1,nlevage_ed - do isc=1,nlevsclass_ed - i=i+1 - fates_hdim_scmap_levscag(i) = isc - fates_hdim_agmap_levscag(i) = iage - end do - end do - - i=0 - do ipft=1,numpft_ed - do ican=1,nclmax - do ileaf=1,nlevleaf - i=i+1 - fates_hdim_canmap_levcnlfpf(i) = ican - fates_hdim_lfmap_levcnlfpf(i) = ileaf - fates_hdim_pftmap_levcnlfpf(i) = ipft - end do - end do - end do - - end subroutine ed_hist_scpfmaps - ! ===================================================================================== function get_age_class_index(age) result( patch_age_class ) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index a971ccd70a..0f93936e2b 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -12,6 +12,7 @@ module FatesHistoryInterfaceMod use FatesInterfaceMod , only : hlm_hio_ignore_val use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_use_ed_st3 + use FatesInterfaceMod , only : numpft use EDParamsMod , only : ED_val_comp_excln ! FIXME(bja, 2016-10) need to remove CLM dependancy @@ -1102,7 +1103,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) use EDtypesMod , only : ncwd use EDtypesMod , only : ican_upper use EDtypesMod , only : ican_ustory - use EDTypesMod , only : maxpft use EDTypesMod , only : get_sizeage_class_index use EDTypesMod , only : nlevleaf @@ -1720,7 +1720,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! pass the cohort termination mortality as a flux to the history, and then reset the termination mortality buffer ! note there are various ways of reporting the total mortality, so pass to these as well - do i_pft = 1, maxpft + do i_pft = 1, numpft do i_scls = 1,nlevsclass_ed i_scpf = (i_pft-1)*nlevsclass_ed + i_scls hio_m6_si_scpf(io_si,i_scpf) = (sites(s)%terminated_nindivs(i_scls,i_pft,1) + & @@ -1738,13 +1738,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) sites(s)%terminated_nindivs(:,:,:) = 0._r8 ! pass the recruitment rate as a flux to the history, and then reset the recruitment buffer - do i_pft = 1, maxpft + do i_pft = 1, numpft hio_recruitment_si_pft(io_si,i_pft) = sites(s)%recruitment_rate(i_pft) * days_per_year end do sites(s)%recruitment_rate(:) = 0._r8 ! summarize all of the mortality fluxes by PFT - do i_pft = 1, maxpft + do i_pft = 1, numpft do i_scls = 1,nlevsclass_ed i_scpf = (i_pft-1)*nlevsclass_ed + i_scls hio_mortality_si_pft(io_si,i_pft) = hio_mortality_si_pft(io_si,i_pft) + & @@ -1815,7 +1815,7 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) AREA_INV, & nlevage_ed, & nlevsclass_ed - use EDTypesMod, only : numpft_ed, nclmax, nlevleaf + use EDTypesMod , only : nclmax, nlevleaf ! ! Arguments class(fates_history_interface_type) :: this @@ -2057,7 +2057,7 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) enddo ! cohort loop ! summarize radiation profiles through the canopy - do ipft=1,numpft_ed + do ipft=1,numpft do ican=1,nclmax do ileaf=1,nlevleaf ! calculate where we are on multiplexed dimensions @@ -2171,6 +2171,7 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) use FatesHydraulicsMemMod, only : nlevsoi_hyd use EDTypesMod , only : nlevsclass_ed use EDTypesMod , only : maxpft + ! Arguments class(fates_history_interface_type) :: this @@ -2380,7 +2381,7 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) end do !patch loop if(hlm_use_ed_st3.eq.ifalse) then - do scpf=1,nlevsclass_ed*maxpft + do scpf=1,nlevsclass_ed*numpft if( abs(hio_nplant_si_scpf(io_si, scpf)-ncohort_scpf(scpf)) > 1.0E-8_r8 ) then write(fates_log(),*) 'nplant check on hio_nplant_si_scpf fails during hydraulics history updates' call endrun(msg=errMsg(sourcefile, __LINE__)) diff --git a/main/FatesIODimensionsMod.F90 b/main/FatesIODimensionsMod.F90 index 1dd5cce0b9..3d2e4bee0a 100644 --- a/main/FatesIODimensionsMod.F90 +++ b/main/FatesIODimensionsMod.F90 @@ -110,6 +110,10 @@ module FatesIODimensionsMod procedure, public :: SetThreadBounds end type fates_io_dimension_type + + + + contains ! ===================================================================================== diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 876358f560..ba4597fd8d 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -17,13 +17,14 @@ module FatesInterfaceMod use EDTypesMod , only : inir use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf - use EDTypesMod , only : numpft_ed + use EDTypesMod , only : maxpft use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue use FatesGlobals , only : fates_global_verbose use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use EDPftvarcon , only : FatesReportPFTParams + use EDPftvarcon , only : EDPftvarcon_inst use EDParamsMod , only : FatesReportParams @@ -102,6 +103,15 @@ module FatesInterfaceMod ! between the pedotransfer functions of the HLM ! and how it moves and stores water in its ! rhizosphere shells + + integer, protected :: hlm_max_patch_per_site ! The HLM needs to exchange some patch + ! level quantities with FATES + ! FATES does not dictate those allocations + ! since it happens pretty early in + ! the model initialization sequence. + ! So we want to at least query it, + ! compare it to our maxpatchpersite, + ! and gracefully halt if we are over-allocating integer, protected :: hlm_use_vertsoilc ! This flag signals whether or not the ! host model is using vertically discretized @@ -170,7 +180,30 @@ module FatesInterfaceMod ! data as some fields are arrays where each array is ! associated with one cohort - + ! ------------------------------------------------------------------------------------- + ! These vectors are used for history output mapping + ! CLM/ALM have limited support for multi-dimensional history output arrays. + ! FATES structure and composition is multi-dimensional, so we end up "multi-plexing" + ! multiple dimensions into one dimension. These new dimensions need definitions, + ! mapping to component dimensions, and definitions for those component dimensions as + ! well. + ! ------------------------------------------------------------------------------------- + + real(r8), allocatable :: fates_hdim_levsclass(:) ! plant size class lower bound dimension + integer , allocatable :: fates_hdim_pfmap_levscpf(:) ! map of pfts into size-class x pft dimension + integer , allocatable :: fates_hdim_scmap_levscpf(:) ! map of size-class into size-class x pft dimension + real(r8), allocatable :: fates_hdim_levage(:) ! patch age lower bound dimension + integer , allocatable :: fates_hdim_levpft(:) ! plant pft dimension + integer , allocatable :: fates_hdim_levfuel(:) ! fire fuel class dimension + integer , allocatable :: fates_hdim_levcwdsc(:) ! cwd class dimension + integer , allocatable :: fates_hdim_levcan(:) ! canopy-layer dimension + integer , allocatable :: fates_hdim_canmap_levcnlf(:) ! canopy-layer map into the canopy-layer x leaf-layer dim + integer , allocatable :: fates_hdim_lfmap_levcnlf(:) ! leaf-layer map into the can-layer x leaf-layer dimension + integer , allocatable :: fates_hdim_canmap_levcnlfpf(:) ! can-layer map into the can-layer x pft x leaf-layer dim + integer , allocatable :: fates_hdim_lfmap_levcnlfpf(:) ! leaf-layer map into the can-layer x pft x leaf-layer dim + integer , allocatable :: fates_hdim_pftmap_levcnlfpf(:) ! pft map into the canopy-layer x pft x leaf-layer dim + integer , allocatable :: fates_hdim_scmap_levscag(:) ! map of size-class into size-class x patch age dimension + integer , allocatable :: fates_hdim_agmap_levscag(:) ! map of patch-age into size-class x patch age dimension ! ------------------------------------------------------------------------------------ ! DYNAMIC BOUNDARY CONDITIONS @@ -195,6 +228,16 @@ module FatesInterfaceMod real(r8), protected :: hlm_freq_day ! fraction of year for daily time-step ! (1/days_per_year_, this is a frequency + + ! ------------------------------------------------------------------------------------- + ! + ! Constant parameters that are dictated by the fates parameter file + ! + ! ------------------------------------------------------------------------------------- + + integer, protected :: numpft ! The total number of PFTs defined in the simulation + + ! ------------------------------------------------------------------------------------- ! Structured Boundary Conditions (SITE/PATCH SCALE) ! For floating point arrays, it is sometimes the convention to define the arrays as @@ -209,10 +252,6 @@ module FatesInterfaceMod ! _rb means radiation band ! ------------------------------------------------------------------------------------ - - - - type, public :: bc_in_type ! The actual number of FATES' ED patches @@ -796,17 +835,65 @@ end subroutine zero_bcs ! =================================================================================== subroutine set_fates_global_elements(use_fates) + + ! -------------------------------------------------------------------------------- + ! + ! This subroutine is called directly from the HLM, and is the first FATES routine + ! that is called. + ! + ! This subroutine MUST BE CALLED AFTER the FATES parameter file has been read in, + ! and the EDPftvarcon_inst structure has been made. + ! This subroutine must ALSO BE CALLED BEFORE the history file dimensions + ! are set. + ! + ! This routine requires no information from the HLM. This routine is responsible + ! for generating the globals that are required by the HLM that are entirely + ! FATES derived. + ! + ! -------------------------------------------------------------------------------- + + implicit none logical,intent(in) :: use_fates ! Is fates turned on? if (use_fates) then + ! Identify the number of PFTs by evaluating a pft array + ! Using wood density as that is not expected to be deprecated any time soon + + if(lbound(EDPftvarcon_inst%wood_density(:),dim=1) .eq. 0 ) then + numpft = size(EDPftvarcon_inst%wood_density,dim=1)-1 + elseif(lbound(EDPftvarcon_inst%wood_density(:),dim=1) .eq. 1 ) then + numpft = size(EDPftvarcon_inst%wood_density,dim=1) + else + write(fates_log(), *) 'While assessing the number of FATES PFTs,' + write(fates_log(), *) 'it was found that the lower bound was neither 0 or 1?' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if(numpft>maxpft) then + write(fates_log(), *) 'The number of PFTs dictated by the FATES parameter file' + write(fates_log(), *) 'is larger than the maximum allowed. Increase the FATES parameter constant' + write(fates_log(), *) 'FatesInterfaceMod.F90:maxpft accordingly' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + + ! These values are used to define the restart file allocations and general structure + ! of memory for the cohort arrays + fates_maxElementsPerPatch = max(maxCohortsPerPatch, & - numpft_ed * nclmax * nlevleaf) - + numpft * nclmax * nlevleaf) + fates_maxElementsPerSite = maxPatchesPerSite * fates_maxElementsPerPatch + + ! Set Various Mapping Arrays used in history output as well + ! These will not be used if use_ed or use_fates is false + call fates_history_maps() + + else ! If we are not using FATES, the cohort dimension is still ! going to be initialized, lets set it to the smallest value @@ -822,6 +909,123 @@ subroutine set_fates_global_elements(use_fates) end subroutine set_fates_global_elements + !============================================================================================== + + subroutine fates_history_maps + + use EDTypesMod, only : nlevsclass_ed + use EDTypesMod, only : NFSC + use EDTypesMod, only : NCWD + use EDTypesMod, only : nlevage_ed + use EDTypesMod, only : nlevsclass_ed + use EDTypesMod, only : nclmax + use EDTypesMod, only : nlevleaf + use EDTypesMod, only : sclass_ed + use EDTypesMod, only : ageclass_ed + + ! ------------------------------------------------------------------------------------------ + ! This subroutine allocates and populates the variables + ! that define the mapping of variables in history files in multiplexed dimensions liked + ! the "scpf" format + ! back to + ! their respective single component dimensions, like size-class "sc" and pft "pf" + ! ------------------------------------------------------------------------------------------ + + integer :: i + integer :: isc + integer :: ipft + integer :: icwd + integer :: ifuel + integer :: ican + integer :: ileaf + integer :: iage + + allocate( fates_hdim_levsclass(1:nlevsclass_ed )) + allocate( fates_hdim_pfmap_levscpf(1:nlevsclass_ed*numpft)) + allocate( fates_hdim_scmap_levscpf(1:nlevsclass_ed*numpft)) + allocate( fates_hdim_levpft(1:numpft )) + allocate( fates_hdim_levfuel(1:NFSC )) + allocate( fates_hdim_levcwdsc(1:NCWD )) + allocate( fates_hdim_levage(1:nlevage_ed )) + + allocate( fates_hdim_levcan(nclmax)) + allocate( fates_hdim_canmap_levcnlf(nlevleaf*nclmax)) + allocate( fates_hdim_lfmap_levcnlf(nlevleaf*nclmax)) + allocate( fates_hdim_canmap_levcnlfpf(nlevleaf*nclmax*numpft)) + allocate( fates_hdim_lfmap_levcnlfpf(nlevleaf*nclmax*numpft)) + allocate( fates_hdim_pftmap_levcnlfpf(nlevleaf*nclmax*numpft)) + allocate( fates_hdim_scmap_levscag(nlevsclass_ed * nlevage_ed )) + allocate( fates_hdim_agmap_levscag(nlevsclass_ed * nlevage_ed )) + + ! Fill the IO array of plant size classes + ! For some reason the history files did not like + ! a hard allocation of sclass_ed + fates_hdim_levsclass(:) = sclass_ed(:) + + fates_hdim_levage(:) = ageclass_ed(:) + + ! make pft array + do ipft=1,numpft + fates_hdim_levpft(ipft) = ipft + end do + + ! make fuel array + do ifuel=1,NFSC + fates_hdim_levfuel(ifuel) = ifuel + end do + + ! make cwd array + do icwd=1,NCWD + fates_hdim_levcwdsc(icwd) = icwd + end do + + ! make canopy array + do ican = 1,nclmax + fates_hdim_levcan(ican) = ican + end do + + ! Fill the IO arrays that match pft and size class to their combined array + i=0 + do ipft=1,numpft + do isc=1,nlevsclass_ed + i=i+1 + fates_hdim_pfmap_levscpf(i) = ipft + fates_hdim_scmap_levscpf(i) = isc + end do + end do + + i=0 + do ican=1,nclmax + do ileaf=1,nlevleaf + i=i+1 + fates_hdim_canmap_levcnlf(i) = ican + fates_hdim_lfmap_levcnlf(i) = ileaf + end do + end do + + i=0 + do iage=1,nlevage_ed + do isc=1,nlevsclass_ed + i=i+1 + fates_hdim_scmap_levscag(i) = isc + fates_hdim_agmap_levscag(i) = iage + end do + end do + + i=0 + do ipft=1,numpft + do ican=1,nclmax + do ileaf=1,nlevleaf + i=i+1 + fates_hdim_canmap_levcnlfpf(i) = ican + fates_hdim_lfmap_levcnlfpf(i) = ileaf + fates_hdim_pftmap_levcnlfpf(i) = ipft + end do + end do + end do + + end subroutine fates_history_maps + ! =================================================================================== subroutine SetFatesTime(current_year_in, current_month_in, & @@ -912,6 +1116,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_hio_ignore_val = unset_double hlm_masterproc = unset_int hlm_ipedof = unset_int + hlm_max_patch_per_site = unset_int hlm_use_vertsoilc = unset_int hlm_use_spitfire = unset_int hlm_use_planthydro = unset_int @@ -1064,6 +1269,21 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if(hlm_max_patch_per_site .eq. unset_int ) then + if (fates_global_verbose()) then + write(fates_log(), *) 'the number of patch-space per site unset: hlm_max_patch_per_site, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + elseif(hlm_max_patch_per_site < maxPatchesPerSite ) then + if (fates_global_verbose()) then + write(fates_log(), *) 'FATES is trying to allocate space for more patches per site, than the HLM has space for.' + write(fates_log(), *) 'hlm_max_patch_per_site (HLM side): ', hlm_max_patch_per_site + write(fates_log(), *) 'maxPatchesPerSite (FATES side): ', maxPatchesPerSite + write(fates_log(), *) + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if(hlm_use_vertsoilc .eq. unset_int) then if (fates_global_verbose()) then write(fates_log(), *) 'switch for the HLMs soil carbon discretization unset: hlm_use_vertsoilc, exiting' @@ -1148,6 +1368,12 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering hlm_ipedof = ',ival,' to FATES' end if + case('max_patch_per_site') + hlm_max_patch_per_site = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_max_patch_per_site = ',ival,' to FATES' + end if + case('use_vertsoilc') hlm_use_vertsoilc = ival if (fates_global_verbose()) then diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 634a0e9ad9..52a75e480c 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -34,7 +34,6 @@ module FatesInventoryInitMod use EDTypesMod , only : ed_cohort_type use EDTypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst - use EDEcophysConType , only : EDecophyscon implicit none private @@ -76,7 +75,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) use shr_file_mod, only : shr_file_getUnit use shr_file_mod, only : shr_file_freeUnit use EDTypesMod, only : nclmax - use EDTypesMod, only : numpft_ed use EDTypesMod, only : maxpft use EDTypesMod, only : ncwd use EDParamsMod, only : ED_val_maxspread @@ -252,12 +250,12 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) spread_init(:) = ED_val_maxspread cwd_ag_init(:) = 0.0_r8 cwd_bg_init(:) = 0.0_r8 - leaf_litter_init(1:numpft_ed) = 0.0_r8 - root_litter_init(1:numpft_ed) = 0.0_r8 + leaf_litter_init(:) = 0.0_r8 + root_litter_init(:) = 0.0_r8 call create_patch(sites(s), newpatch, age_init, area_init, spread_init, & cwd_ag_init, cwd_bg_init, & - leaf_litter_init(1:numpft_ed), root_litter_init(1:numpft_ed) ) + leaf_litter_init, root_litter_init ) if( inv_format_list(invsite) == 1 ) then call set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) @@ -644,7 +642,6 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name use EDTypesMod, only: get_age_class_index use EDtypesMod, only: AREA - use EDTypesMod, only: numpft_ed use EDTypesMod, only: ncwd use SFParamsMod , only : SF_val_CWD_frac @@ -719,10 +716,10 @@ subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name newpatch%cwd_bg(icwd) = 0.0_r8 end do - do ipft = 1, numpft_ed - newpatch%leaf_litter(ipft) = 0.0_r8 - newpatch%root_litter(ipft) = 0.0_r8 - end do + + newpatch%leaf_litter(:) = 0.0_r8 + newpatch%root_litter(:) = 0.0_r8 + return end subroutine set_inventory_edpatch_type1 @@ -754,11 +751,11 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! avgRG (cm/yr?) Average Radial Growth (NOT USED) ! -------------------------------------------------------------------------------------------- - use EDTypesMod , only : numpft_ed use EDGrowthFunctionsMod, only : hite use EDGrowthFunctionsMod, only : bleaf use EDGrowthFunctionsMod, only : bdead use EDCohortDynamicsMod , only : create_cohort + use FatesInterfaceMod , only : numpft ! Arguments type(ed_site_type),intent(inout), target :: csite ! current site @@ -824,10 +821,10 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & ! pft, nplant and dbh are the critical ones in this format specification ! ------------------------------------------------------------------------------------------- - if (c_pft > numpft_ed ) then + if (c_pft > numpft ) then write(fates_log(), *) 'inventory pft: ',c_pft write(fates_log(), *) 'An inventory cohort file specified a pft index' - write(fates_log(), *) 'greater than the maximum specified pfts ed_numpft' + write(fates_log(), *) 'greater than the maximum specified pfts numpft' call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -877,13 +874,13 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & temp_cohort%b = temp_cohort%balive + temp_cohort%bdead if( EDPftvarcon_inst%evergreen(c_pft) == 1) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(c_pft) + temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(c_pft) temp_cohort%laimemory = 0._r8 cstatus = 2 endif if( EDPftvarcon_inst%season_decid(c_pft) == 1 ) then !for dorment places - temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(c_pft) !stored carbon in new seedlings. + temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(c_pft) !stored carbon in new seedlings. if(csite%status == 2)then temp_cohort%laimemory = 0.0_r8 else @@ -895,7 +892,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & endif if ( EDPftvarcon_inst%stress_decid(c_pft) == 1 ) then - temp_cohort%bstore = Bleaf(temp_cohort) * EDecophyscon%cushion(c_pft) + temp_cohort%bstore = Bleaf(temp_cohort) * EDPftvarcon_inst%cushion(c_pft) temp_cohort%laimemory = Bleaf(temp_cohort) temp_cohort%balive = temp_cohort%balive - temp_cohort%laimemory cstatus = csite%dstatus diff --git a/main/FatesParameterDerivedMod.F90 b/main/FatesParameterDerivedMod.F90 index bd177af2b5..ee1ebbd8b7 100644 --- a/main/FatesParameterDerivedMod.F90 +++ b/main/FatesParameterDerivedMod.F90 @@ -36,28 +36,28 @@ module FatesParameterDerivedMod contains - subroutine InitAllocate(this,maxpft) + subroutine InitAllocate(this,numpft) class(param_derived_type), intent(inout) :: this - integer, intent(in) :: maxpft + integer, intent(in) :: numpft - allocate(this%vcmax25top(maxpft)) - allocate(this%jmax25top(maxpft)) - allocate(this%tpu25top(maxpft)) - allocate(this%kp25top(maxpft)) - allocate(this%lmr25top(maxpft)) + allocate(this%vcmax25top(numpft)) + allocate(this%jmax25top(numpft)) + allocate(this%tpu25top(numpft)) + allocate(this%kp25top(numpft)) + allocate(this%lmr25top(numpft)) return end subroutine InitAllocate ! ===================================================================================== - subroutine Init(this,maxpft) + subroutine Init(this,numpft) use EDPftvarcon, only: EDPftvarcon_inst class(param_derived_type), intent(inout) :: this - integer, intent(in) :: maxpft + integer, intent(in) :: numpft ! local variables integer :: ft ! pft index @@ -70,9 +70,9 @@ subroutine Init(this,maxpft) ! projected area basis [m^2/gC] leafcn => EDPftvarcon_inst%leafcn ) ! leaf C:N (gC/gN) - call this%InitAllocate(maxpft) + call this%InitAllocate(numpft) - do ft = 1,maxpft + do ft = 1,numpft ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) lnc = 1._r8 / (slatop(ft) * leafcn(ft)) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 12b1d6bc9b..347bfc9052 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -932,7 +932,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) use EDTypesMod, only : nclmax use EDTypesMod, only : nlevleaf use FatesInterfaceMod, only : fates_maxElementsPerPatch - use EDTypesMod, only : numpft_ed + use FatesInterfaceMod, only : numpft use EDTypesMod, only : ed_site_type use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type @@ -1083,7 +1083,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_pa_sunz = io_idx_co_1st ! write seed_bank info(site-level, but PFT-resolved) - do i = 1,numpft_ed + do i = 1,numpft rio_seed_bank_sift(io_idx_co_1st+i-1) = sites(s)%seed_bank(i) end do @@ -1184,9 +1184,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! ! deal with patch level fields of arrays here ! - ! these are arrays of length numpft_ed, each patch contains one + ! these are arrays of length numpft, each patch contains one ! vector so we increment - do i = 1,numpft_ed + do i = 1,numpft rio_leaf_litter_paft(io_idx_pa_pft) = cpatch%leaf_litter(i) rio_root_litter_paft(io_idx_pa_pft) = cpatch%root_litter(i) rio_leaf_litter_in_paft(io_idx_pa_pft) = cpatch%leaf_litter_in(i) @@ -1207,10 +1207,10 @@ subroutine set_restart_vectors(this,nc,nsites,sites) if ( DEBUG ) write(fates_log(),*) 'CLTV io_idx_pa_sunz 1 ',io_idx_pa_sunz - if ( DEBUG ) write(fates_log(),*) 'CLTV 1186 ',nlevleaf,numpft_ed,nclmax + if ( DEBUG ) write(fates_log(),*) 'CLTV 1186 ',nlevleaf,numpft,nclmax - do k = 1,nlevleaf ! nlevleaf currently 40 - do j = 1,numpft_ed ! numpft_ed currently 2 + do k = 1,nlevleaf ! nlevleaf currently 40 + do j = 1,numpft ! dependent on parameter file do i = 1,nclmax ! nclmax currently 2 rio_fsun_paclftls(io_idx_pa_sunz) = cpatch%f_sun(i,j,k) rio_fabd_sun_z_paclftls(io_idx_pa_sunz) = cpatch%fabd_sun_z(i,j,k) @@ -1229,9 +1229,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! by the maximum number of cohorts per patch io_idx_co_1st = io_idx_co_1st + fates_maxElementsPerPatch - ! reset counters so that they are all advanced evenly. Currently - ! the offset is 10, the max of numpft_ed, ncwd, nclmax, - ! io_idx_si_wmem and the number of allowed cohorts per patch + ! reset counters so that they are all advanced evenly. io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st io_idx_pa_cl = io_idx_co_1st @@ -1308,7 +1306,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) use EDTypesMod, only : nlevleaf use EDTypesMod, only : nclmax use FatesInterfaceMod, only : fates_maxElementsPerPatch - use EDTypesMod, only : numpft_ed + use EDTypesMod, only : maxpft use EDTypesMod, only : area use EDPatchDynamicsMod, only : zero_patch use EDGrowthFunctionsMod, only : Dbh @@ -1332,8 +1330,8 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) real(r8) :: cwd_ag_local(ncwd) real(r8) :: cwd_bg_local(ncwd) real(r8) :: spread_local(nclmax) - real(r8) :: leaf_litter_local(numpft_ed) - real(r8) :: root_litter_local(numpft_ed) + real(r8) :: leaf_litter_local(maxpft) + real(r8) :: root_litter_local(maxpft) real(r8) :: patch_age integer :: cohortstatus integer :: s ! site index @@ -1503,10 +1501,10 @@ subroutine get_restart_vectors(this, nc, nsites, sites) use EDTypesMod, only : ed_site_type use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type - use EDTypesMod, only : numpft_ed use EDTypesMod, only : ncwd use EDTypesMod, only : nlevleaf use EDTypesMod, only : nclmax + use FatesInterfaceMod, only : numpft use FatesInterfaceMod, only : fates_maxElementsPerPatch use EDTypesMod, only : numWaterMem @@ -1644,7 +1642,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_wmem = io_idx_co_1st ! read seed_bank info(site-level, but PFT-resolved) - do i = 1,numpft_ed + do i = 1,numpft sites(s)%seed_bank(i) = rio_seed_bank_sift(io_idx_co_1st+i-1) enddo @@ -1745,10 +1743,10 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! ! deal with patch level fields of arrays here ! - ! these are arrays of length numpft_ed, each patch contains one + ! these are arrays of length numpft, each patch contains one ! vector so we increment - do i = 1,numpft_ed + do i = 1,numpft cpatch%leaf_litter(i) = rio_leaf_litter_paft(io_idx_pa_pft) cpatch%root_litter(i) = rio_root_litter_paft(io_idx_pa_pft) cpatch%leaf_litter_in(i) = rio_leaf_litter_in_paft(io_idx_pa_pft) @@ -1770,7 +1768,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) if ( DEBUG ) write(fates_log(),*) 'CVTL io_idx_pa_sunz 1 ',io_idx_pa_sunz do k = 1,nlevleaf ! nlevleaf currently 40 - do j = 1,numpft_ed ! numpft_ed currently 2 + do j = 1,numpft do i = 1,nclmax ! nclmax currently 2 cpatch%f_sun(i,j,k) = rio_fsun_paclftls(io_idx_pa_sunz) cpatch%fabd_sun_z(i,j,k) = rio_fabd_sun_z_paclftls(io_idx_pa_sunz)