diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index dd3b26fd28..caf6861577 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -504,10 +504,10 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) allocate(copyc) call InitPRTCohort(copyc) - call copy_cohort(currentCohort, copyc) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(currentSite,copyc) endif + call copy_cohort(currentCohort, copyc) newarea = currentCohort%c_area - cc_loss copyc%n = currentCohort%n*newarea/currentCohort%c_area @@ -851,11 +851,11 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) allocate(copyc) call InitPRTCohort(copyc) - call copy_cohort(currentCohort, copyc) !makes an identical copy... if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(CurrentSite,copyc) endif - + call copy_cohort(currentCohort, copyc) !makes an identical copy... + newarea = currentCohort%c_area - cc_gain !new area of existing cohort call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & @@ -1422,11 +1422,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! is obscured by snow. layer_top_hite = currentCohort%hite - & - ( dble(iv-1.0)/currentCohort%NV * currentCohort%hite * & + ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite * & EDPftvarcon_inst%crown(currentCohort%pft) ) layer_bottom_hite = currentCohort%hite - & - ( dble(iv)/currentCohort%NV * currentCohort%hite * & + ( real(iv,r8)/currentCohort%NV * currentCohort%hite * & EDPftvarcon_inst%crown(currentCohort%pft) ) fraction_exposed = 1.0_r8 @@ -1449,7 +1449,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) if(iv==currentCohort%NV) then remainder = (currentCohort%treelai + currentCohort%treesai) - & - (dinc_ed*dble(currentCohort%nv-1.0_r8)) + (dinc_ed*real(currentCohort%nv-1,r8)) if(remainder > dinc_ed )then write(fates_log(), *)'ED: issue with remainder', & currentCohort%treelai,currentCohort%treesai,dinc_ed, & diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 0c535751bb..d09b90393e 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -372,6 +372,7 @@ subroutine nan_cohort(cc_p) currentCohort%NV = fates_unset_int ! Number of leaf layers: - currentCohort%status_coh = fates_unset_int ! growth status of plant (2 = leaves on , 1 = leaves off) currentCohort%size_class = fates_unset_int ! size class index + currentCohort%size_class_lasttimestep = fates_unset_int ! size class index currentCohort%size_by_pft_class = fates_unset_int ! size by pft classification index currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) @@ -476,6 +477,8 @@ subroutine zero_cohort(cc_p) currentcohort%ts_net_uptake(:) = 0._r8 currentcohort%seed_prod = 0._r8 currentcohort%fraction_crown_burned = 0._r8 + currentCohort%size_class = 1 + currentCohort%size_class_lasttimestep = 0 currentcohort%npp_acc_hold = 0._r8 currentcohort%gpp_acc_hold = 0._r8 currentcohort%dmort = 0._r8 @@ -489,7 +492,6 @@ subroutine zero_cohort(cc_p) currentcohort%prom_weight = 0._r8 currentcohort%crownfire_mort = 0._r8 currentcohort%cambial_mort = 0._r8 - end subroutine zero_cohort @@ -728,6 +730,9 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) real(r8) :: diff real(r8) :: dynamic_fusion_tolerance + integer :: largersc, smallersc, sc_i ! indices for tracking the growth flux caused by fusion + real(r8) :: larger_n, smaller_n + logical, parameter :: fuse_debug = .false. ! This debug is over-verbose ! and gets its own flag @@ -847,6 +852,42 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%canopy_layer_yesterday = (currentCohort%n*currentCohort%canopy_layer_yesterday + & nextc%n*nextc%canopy_layer_yesterday)/newn + ! keep track of the size class bins so that we can monitor growth fluxes + ! compare the values. if they are the same, then nothing needs to be done. if not, track the diagnostic flux + if (currentCohort%size_class_lasttimestep .ne. nextc%size_class_lasttimestep ) then + ! + ! keep track of which was which, irresespective of which cohort they were in + if (currentCohort%size_class_lasttimestep .gt. nextc%size_class_lasttimestep) then + largersc = currentCohort%size_class_lasttimestep + smallersc = nextc%size_class_lasttimestep + larger_n = currentCohort%n + smaller_n = nextc%n + else + largersc = nextc%size_class_lasttimestep + smallersc = currentCohort%size_class_lasttimestep + larger_n = nextc%n + smaller_n = currentCohort%n + endif + ! + ! it is possible that fusion has caused cohorts separated by at least two size bin deltas to join. + ! so slightly complicated to keep track of because the resulting cohort could be in one of the old bins or in between + ! structure as a loop to handle the general case + ! + ! first the positive growth case + do sc_i = smallersc + 1, currentCohort%size_class + currentSite%growthflux_fusion(sc_i, currentCohort%pft) = & + currentSite%growthflux_fusion(sc_i, currentCohort%pft) + smaller_n + end do + ! + ! next the negative growth case + do sc_i = currentCohort%size_class + 1, largersc + currentSite%growthflux_fusion(sc_i, currentCohort%pft) = & + currentSite%growthflux_fusion(sc_i, currentCohort%pft) - larger_n + end do + ! now that we've tracked the change flux. reset the memory of the prior timestep + currentCohort%size_class_lasttimestep = currentCohort%size_class + endif + ! Flux and biophysics variables have not been calculated for recruits we just default to ! their initization values, which should be the same for eahc @@ -1215,6 +1256,7 @@ subroutine copy_cohort( currentCohort,copyc ) n%excl_weight = o%excl_weight n%prom_weight = o%prom_weight n%size_class = o%size_class + n%size_class_lasttimestep = o%size_class_lasttimestep n%size_by_pft_class = o%size_by_pft_class ! This transfers the PRT objects over. @@ -1292,6 +1334,7 @@ subroutine copy_cohort( currentCohort,copyc ) ! indices for binning n%size_class = o%size_class + n%size_class_lasttimestep = o%size_class_lasttimestep n%size_by_pft_class = o%size_by_pft_class !Pointers diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index ad81833e0c..2c4cf1b7e7 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -350,12 +350,17 @@ subroutine spawn_patches( currentSite, bc_in) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate + ! Only create new patches that have non-negligible amount of land + if((currentPatch%area*currentPatch%disturbance_rate) > nearzero ) then + site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate + end if + currentPatch => currentPatch%older enddo ! end loop over patches. sum area disturbed for all patches. - if (site_areadis > 0.0_r8) then + if (site_areadis > nearzero) then + cwd_ag_local = 0.0_r8 cwd_bg_local = 0.0_r8 leaf_litter_local = 0.0_r8 @@ -377,6 +382,8 @@ subroutine spawn_patches( currentSite, bc_in) ! This is the amount of patch area that is disturbed, and donated by the donor patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate + if (patch_site_areadis > nearzero) then + call average_patch_properties(currentPatch, new_patch, patch_site_areadis) if (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. & @@ -543,7 +550,7 @@ subroutine spawn_patches( currentSite, bc_in) nc%n * currentCohort%crownfire_mort / hlm_freq_day currentSite%fmort_carbonflux(levcan) = currentSite%fmort_carbonflux(levcan) + & (nc%n * currentCohort%fire_mort) * & - currentCohort%b_total() * g_per_kg * days_per_sec * ha_per_m2 + total_c * g_per_kg * days_per_sec * ha_per_m2 ! loss of individual from fire in new patch. nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) @@ -690,10 +697,6 @@ subroutine spawn_patches( currentSite, bc_in) enddo ! currentCohort call sort_cohorts(currentPatch) - !zero disturbance accumulators - currentPatch%disturbance_rate = 0._r8 - currentPatch%disturbance_rates = 0._r8 - !update area of donor patch currentPatch%area = currentPatch%area - patch_site_areadis @@ -706,6 +709,12 @@ subroutine spawn_patches( currentSite, bc_in) call terminate_cohorts(currentSite, currentPatch, 2) call sort_cohorts(currentPatch) + end if ! if (patch_site_areadis > nearzero) then + + !zero disturbance rate trackers + currentPatch%disturbance_rate = 0._r8 + currentPatch%disturbance_rates = 0._r8 + currentPatch => currentPatch%younger enddo ! currentPatch patch loop. @@ -731,9 +740,11 @@ subroutine spawn_patches( currentSite, bc_in) endif !end new_patch area + call check_patch_area(currentSite) call set_patchno(currentSite) + return end subroutine spawn_patches ! ============================================================================ @@ -1153,14 +1164,20 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat !not right to recalcualte dmort here. canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * hlm_freq_day * fates_mortality_disturbance_fraction) - - canopy_mortality_woody_litter(p)= canopy_mortality_woody_litter(p) + & canopy_dead*(struct_c + sapw_c) canopy_mortality_leaf_litter(p) = canopy_mortality_leaf_litter(p) + & canopy_dead*leaf_c + + ! Some plants upon death will transfer storage carbon to seed production + ! Storage carbon that is not transferred to seeds goes to root litter flux + canopy_mortality_root_litter(p) = canopy_mortality_root_litter(p) + & - canopy_dead*(fnrt_c + store_c) + canopy_dead*(fnrt_c + store_c*(1.0_r8-EDPftvarcon_inst%allom_frbstor_repro(p)) ) + + currentSite%seed_bank(p) = currentSite%seed_bank(p) + & + canopy_dead * store_c * EDPftvarcon_inst%allom_frbstor_repro(p)/AREA + if( hlm_use_planthydro == itrue ) then call AccumulateMortalityWaterStorage(currentSite,currentCohort, canopy_dead) @@ -1324,6 +1341,8 @@ subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_l new_patch%frac_burnt = 0._r8 new_patch%total_tree_area = 0.0_r8 new_patch%NCL_p = 1 + + end subroutine create_patch @@ -1448,6 +1467,12 @@ subroutine zero_patch(cp_p) currentPatch%c_stomata = 0.0_r8 ! This is calculated immediately before use currentPatch%c_lblayer = 0.0_r8 + currentPatch%solar_zenith_flag = .false. + currentPatch%solar_zenith_angle = nan + + currentPatch%gnd_alb_dir(:) = nan + currentPatch%gnd_alb_dif(:) = nan + end subroutine zero_patch ! ============================================================================ diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index bd92830686..ad639950aa 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -172,7 +172,7 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in ) ! update fragmenting pool fluxes call cwd_input( currentSite, currentPatch) call cwd_out( currentSite, currentPatch, bc_in) - + do p = 1,numpft currentSite%dseed_dt(p) = currentSite%dseed_dt(p) + & (currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - & @@ -589,7 +589,7 @@ subroutine phenology( currentSite, bc_in ) if ((t >= currentSite%dleafondate - 30.and.t <= currentSite%dleafondate + 30).or.(t > 360 - 15.and. & currentSite%dleafondate < 15))then ! are we in the window? ! TODO: CHANGE THIS MATH, MOVE THE DENOMENATOR OUTSIDE OF THE SUM (rgk 01-2017) - if (sum(currentSite%water_memory(1:numWaterMem)/dble(numWaterMem)) & + if (sum(currentSite%water_memory(1:numWaterMem)/real(numWaterMem,r8)) & >= ED_val_phen_drought_threshold.and.currentSite%dstatus == 1.and.t >= 10)then ! leave some minimum time between leaf off and leaf on to prevent 'flickering'. if (timesincedleafoff > ED_val_phen_doff_time)then @@ -792,6 +792,7 @@ subroutine seeds_in( currentSite, cp_pnt ) type(ed_cohort_type), pointer :: currentCohort integer :: p logical :: pft_present(maxpft) + real(r8) :: store_c_to_repro ! carbon sent from storage to reproduction upon death [kg/plant] real(r8) :: npfts_present !---------------------------------------------------------------------- @@ -822,10 +823,18 @@ subroutine seeds_in( currentSite, cp_pnt ) currentPatch => cp_pnt currentCohort => currentPatch%tallest do while (associated(currentCohort)) + + ! a certain fraction of bstore goes to clonal reproduction when plants die + store_c_to_repro = currentCohort%prt%GetState(store_organ,all_carbon_elements) * & + EDPftvarcon_inst%allom_frbstor_repro(currentCohort%pft) + 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) + + currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + & + (currentCohort%seed_prod * currentCohort%n - & + currentCohort%dndt*store_c_to_repro) & + /(currentPatch%area * npfts_present) endif end do currentCohort => currentCohort%shorter @@ -836,8 +845,15 @@ subroutine seeds_in( currentSite, cp_pnt ) currentCohort => currentPatch%tallest do while (associated(currentCohort)) p = currentCohort%pft + + ! a certain fraction of bstore goes to clonal reproduction when plants die + store_c_to_repro = currentCohort%prt%GetState(store_organ,all_carbon_elements) * & + EDPftvarcon_inst%allom_frbstor_repro(p) + currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + & - currentCohort%seed_prod * currentCohort%n/currentPatch%area + (currentCohort%seed_prod * currentCohort%n - & + currentCohort%dndt*store_c_to_repro)/currentPatch%area + currentCohort => currentCohort%shorter enddo !cohort loop @@ -913,13 +929,18 @@ subroutine seed_germination( currentSite, currentPatch ) do p = 1,numpft currentPatch%seed_germination(p) = min(currentSite%seed_bank(p) * & - EDPftvarcon_inst%germination_timescale(p),max_germination) + EDPftvarcon_inst%germination_timescale(p),max_germination) + !set the germination only under the growing season...c.xu + if (EDPftvarcon_inst%season_decid(p) == 1.and.currentSite%status == 1)then + currentPatch%seed_germination(p) = 0.0_r8 + endif + if (EDPftvarcon_inst%stress_decid(p) == 1.and.currentSite%dstatus == 1)then + currentPatch%seed_germination(p) = 0.0_r8 + endif enddo end subroutine seed_germination - ! ============================================================================ - subroutine recruitment( currentSite, currentPatch, bc_in ) ! ! !DESCRIPTION: @@ -1134,9 +1155,8 @@ subroutine CWD_Input( currentSite, currentPatch) ! the litter flux has already been counted since it captured ! the losses of live trees and those flagged for death - currentPatch%root_litter_in(pft) = currentPatch%root_litter_in(pft) + & - (fnrt_c + store_c ) * dead_n + (fnrt_c + store_c*(1._r8-EDPftvarcon_inst%allom_frbstor_repro(pft)) ) * dead_n ! Update diagnostics that track resource management currentSite%resources_management%delta_litter_stock = & diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 81881a190a..7bb2f0ec9c 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2156,7 +2156,7 @@ subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h ) ! of the diameter increment counter = 0 step_frac = step_frac0 - do while( (bdead-bt_dead) > calloc_abs_error ) + do while( (bdead-bt_dead) > calloc_abs_error .and. dbt_dead_dd>0.0_r8) ! vulnerable to div0 dd = step_frac*(bdead-bt_dead)/dbt_dead_dd diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 7f046dfe68..2b10f18899 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -25,6 +25,7 @@ module EDBtranMod public :: btran_ed public :: get_active_suction_layers + public :: check_layer_water contains @@ -197,7 +198,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) cpatch%rootr_ft(ft,j) * pftgs(ft)/sum_pftgs else bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & - cpatch%rootr_ft(ft,j) * 1._r8/dble(numpft) + cpatch%rootr_ft(ft,j) * 1._r8/real(numpft,r8) end if enddo enddo diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index 38c087f942..d9e4c6b3f2 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -14,6 +14,8 @@ module EDSurfaceRadiationMod use EDTypesMod , only : maxPatchesPerSite use EDTypesMod , only : maxpft use FatesConstantsMod , only : r8 => fates_r8 + use FatesConstantsMod , only : itrue + use FatesConstantsMod , only : pi_const use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : bc_out_type use FatesInterfaceMod , only : hlm_numSWb @@ -29,7 +31,8 @@ module EDSurfaceRadiationMod use EDTypesMod , only : ipar use EDCanopyStructureMod, only: calc_areaindex use FatesGlobals , only : fates_log - + use FatesGlobals, only : endrun => fates_endrun + ! CIME globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -37,6 +40,7 @@ module EDSurfaceRadiationMod private public :: ED_Norman_Radiation ! Surface albedo and two-stream fluxes + public :: PatchNormanRadiation public :: ED_SunShadeFracs logical :: debug = .false. ! for debugging this module @@ -66,84 +70,18 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) ! !LOCAL VARIABLES: - ! ============================================================================ - ! ED/NORMAN RADIATION DECS - ! ============================================================================ - type (ed_patch_type) , pointer :: currentPatch - integer :: radtype, L, ft, j, ifp - integer :: iter ! Iteration index - integer :: irep ! Flag to exit iteration loop - 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,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) - real(r8) :: weighted_fsun(nclmax) - real(r8) :: weighted_dif_ratio(nclmax,maxSWb) - real(r8) :: weighted_dif_down(nclmax) - real(r8) :: weighted_dif_up(nclmax) - 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,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(nclmax) - - integer :: fp,iv,s ! array indices - integer :: ib ! waveband number - real(r8) :: cosz ! 0.001 <= coszen <= 1.000 - real(r8) :: chil(maxPatchesPerSite) ! -0.4 <= xl <= 0.6 - real(r8) :: gdir(maxPatchesPerSite) ! leaf projection in solar direction (0 to 1) + integer :: s ! site loop counter + integer :: ifp ! patch loop counter + integer :: ib ! radiation broad band counter + type(ed_patch_type), pointer :: currentPatch ! patch pointer !----------------------------------------------------------------------- - - associate(& - rhol => EDPftvarcon_inst%rhol , & ! Input: [real(r8) (:) ] leaf reflectance: 1=vis, 2=nir - rhos => EDPftvarcon_inst%rhos , & ! Input: [real(r8) (:) ] stem reflectance: 1=vis, 2=nir - taul => EDPftvarcon_inst%taul , & ! Input: [real(r8) (:) ] leaf transmittance: 1=vis, 2=nir - taus => EDPftvarcon_inst%taus , & ! Input: [real(r8) (:) ] stem transmittance: 1=vis, 2=nir - xl => EDPftvarcon_inst%xl , & ! Input: [real(r8) (:) ] ecophys const - leaf/stem orientation index - clumping_index => EDPftvarcon_inst%clumping_index) ! Input: [real(r8) (:) ] ecophys const - leaf/stem clumping index - -! albd => surfalb_inst%albd_patch , & ! Output: [real(r8) (:,:) ] surface albedo (direct) (USED IN LND2ATM,BALANCE_CHECK) -! albi => surfalb_inst%albi_patch , & ! Output: [real(r8) (:,:) ] surface albedo (diffuse) (LND2ATM,BALANCE_CHECK) -! fabd => surfalb_inst%fabd_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by canopy per unit direct flux (BALANCE_CHECK) -! fabi => surfalb_inst%fabi_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by canopy per unit diffuse flux (BALANCE_CHECK) -! ftdd => surfalb_inst%ftdd_patch , & ! Output: [real(r8) (:,:) ] down direct flux below canopy per unit direct flx (BALANCE_CHECK) -! ftid => surfalb_inst%ftid_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit direct flx (BALANCE_CHECK) -! ftii => surfalb_inst%ftii_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit diffuse flx (BALANCE_CHECK) - - ! ------------------------------------------------------------------------------- - ! TODO (mv, 2014-10-29) the filter here is different than below - ! this is needed to have the VOC's be bfb - this needs to be - ! re-examined int he future - ! RGK,2016-08-06: FATES is still incompatible with VOC emission module - ! ------------------------------------------------------------------------------- + ! ------------------------------------------------------------------------------- + ! TODO (mv, 2014-10-29) the filter here is different than below + ! this is needed to have the VOC's be bfb - this needs to be + ! re-examined int he future + ! RGK,2016-08-06: FATES is still incompatible with VOC emission module + ! ------------------------------------------------------------------------------- do s = 1, nsites @@ -167,24 +105,17 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) currentPatch%nrmlzd_parprof_dir_z(:,:,:) = 0._r8 currentPatch%nrmlzd_parprof_dif_z(:,:,:) = 0._r8 - if(bc_in(s)%filter_vegzen_pa(ifp))then + currentPatch%solar_zenith_flag = bc_in(s)%filter_vegzen_pa(ifp) + currentPatch%solar_zenith_angle = bc_in(s)%coszen_pa(ifp) + currentPatch%gnd_alb_dif(1:hlm_numSWb) = bc_in(s)%albgr_dif_rb(1:hlm_numSWb) + currentPatch%gnd_alb_dir(1:hlm_numSWb) = bc_in(s)%albgr_dir_rb(1:hlm_numSWb) + + if(currentPatch%solar_zenith_flag )then - weighted_dir_tr(:) = 0._r8 - weighted_dif_down(:) = 0._r8 - weighted_dif_up(:) = 0._r8 bc_out(s)%albd_parb(ifp,:) = 0._r8 ! output HLM bc_out(s)%albi_parb(ifp,:) = 0._r8 ! output HLM bc_out(s)%fabi_parb(ifp,:) = 0._r8 ! output HLM bc_out(s)%fabd_parb(ifp,:) = 0._r8 ! output HLM - tr_dir_z(:,:,:) = 0._r8 - tr_dif_z(:,:,:) = 0._r8 - ftweight(:,:,:) = 0._r8 - lai_change(:,:,:) = 0._r8 - Dif_up(:,:,:) = 0._r8 - Dif_dn(:,:,:) = 0._r8 - refl_dif(:,:,:,:) = 0.0_r8 - tran_dif(:,:,:,:) = 0.0_r8 - dif_ratio(:,:,:,:) = 0.0_r8 bc_out(s)%ftdd_parb(ifp,:) = 1._r8 ! output HLM bc_out(s)%ftid_parb(ifp,:) = 1._r8 ! output HLM bc_out(s)%ftii_parb(ifp,:) = 1._r8 ! output HLM @@ -201,750 +132,904 @@ subroutine ED_Norman_Radiation (nsites, sites, bc_in, bc_out ) bc_out(s)%ftid_parb(ifp,ib)= 1.0_r8 bc_out(s)%ftii_parb(ifp,ib)= 1.0_r8 enddo + else + + call PatchNormanRadiation (currentPatch, & + bc_out(s)%albd_parb(ifp,:), & + bc_out(s)%albi_parb(ifp,:), & + bc_out(s)%fabd_parb(ifp,:), & + bc_out(s)%fabi_parb(ifp,:), & + bc_out(s)%ftdd_parb(ifp,:), & + bc_out(s)%ftid_parb(ifp,:), & + bc_out(s)%ftii_parb(ifp,:)) - ! Is this pft/canopy layer combination present in this patch? - do L = 1,nclmax - do ft = 1,numpft - currentPatch%canopy_mask(L,ft) = 0 - do iv = 1, currentPatch%nrad(L,ft) - if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then - currentPatch%canopy_mask(L,ft) = 1 - !I think 'present' is only used here... - endif - end do !iv - end do !ft - end do !L - - !do this once for one unit of diffuse, and once for one unit of direct radiation - do radtype = 1, n_rad_stream_types - do ib = 1,hlm_numSWb - if (radtype == idirect) then - ! Set the hypothetical driving radiation. We do this once for a single unit of direct and - ! once for a single unit of diffuse radiation. - forc_dir(ifp,ib) = 1.00_r8 - forc_dif(ifp,ib) = 0.00_r8 - else !dif - forc_dir(ifp,ib) = 0.00_r8 - forc_dif(ifp,ib) = 1.00_r8 - end if - end do !ib - - !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 - 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) - end do !iv - end do !ft1 - end do !L - if (sum(ftweight(1,:,1))<0.999_r8)then - write(fates_log(),*) 'canopy not full',ftweight(1,:,1) - endif - if (sum(ftweight(1,:,1))>1.0001_r8)then - write(fates_log(),*) 'canopy too full',ftweight(1,:,1) - endif - - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! 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 - 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 - chil(ifp) = 0.01_r8 - end if - phi1b(ifp,ft) = 0.5_r8 - 0.633_r8*chil(ifp) - 0.330_r8*chil(ifp)*chil(ifp) - phi2b(ifp,ft) = 0.877_r8 * (1._r8 - 2._r8*phi1b(ifp,ft)) !0 = horiz leaves, 1 - vert leaves. - gdir(ifp) = phi1b(ifp,ft) + phi2b(ifp,ft) * sin(sb) - !how much direct light penetrates a singleunit of lai? - k_dir(ft) = clumping_index(ft) * gdir(ifp) / sin(sb) - end do !FT - - do L = 1,currentPatch%NCL_p !start at the top canopy layer (1 is the top layer.) - weighted_dir_tr(L) = 0.0_r8 - 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 - if (currentPatch%canopy_mask(L,ft) == 1)then !only do calculation if there are the appropriate leaves. - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Diffuse transmittance, tr_dif, do each layer with thickness elai_z. - ! Estimated do nine sky angles in increments of 10 degrees - ! PFT specific... - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - tr_dif_z(L,ft,:) = 0._r8 - do iv = 1,currentPatch%nrad(L,ft) - do j = 1,9 - angle = (5._r8 + (j - 1) * 10._r8) * 3.142 / 180._r8 - gdir(ifp) = phi1b(ifp,ft) + phi2b(ifp,ft) * sin(angle) - tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) + exp(-clumping_index(ft) * & - gdir(ifp) / sin(angle) * & - (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv))) * & - sin(angle)*cos(angle) - end do - - tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) * 2._r8 * (10.00*pi/180._r8) - - end do - - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Direct beam transmittance, tr_dir_z, uses cumulative LAI above layer J to give - ! unscattered direct beam onto layer J. do each PFT section. - ! This is just an decay curve based on k_dir. (leaf & sun angle) - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - if (L==1)then - tr_dir_z(L,ft,1) = 1._r8 - else - tr_dir_z(L,ft,1) = weighted_dir_tr(L-1) - endif - laisum = 0.00_r8 - !total direct beam getting to the bottom of the top canopy. - do iv = 1,currentPatch%nrad(L,ft) - laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) - lai_change(L,ft,iv) = 0.0_r8 - if (( ftweight(L,ft,iv+1) > 0.0_r8 ) .and. ( ftweight(L,ft,iv+1) < ftweight(L,ft,iv) ))then - !where there is a partly empty leaf layer, some fluxes go straight through. - lai_change(L,ft,iv) = ftweight(L,ft,iv)-ftweight(L,ft,iv+1) - endif - if (ftweight(L,ft,iv+1) - ftweight(L,ft,iv) > 1.e-10_r8)then - write(fates_log(),*) 'lower layer has more coverage. This is wrong' , & - ftweight(L,ft,iv),ftweight(L,ft,iv+1),ftweight(L,ft,iv+1)-ftweight(L,ft,iv) - endif - - !n.b. in theory lai_change could be calculated daily in the ED code. - !This is light coming striaght through the canopy. - if (L==1)then - tr_dir_z(L,ft,iv+1) = exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - else - tr_dir_z(L,ft,iv+1) = weighted_dir_tr(L-1)*exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - endif - - if (iv == 1)then - !this is the top layer. - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & - ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) - else - !the lai_change(iv) affects the light incident on layer iv+2 not iv+1 - ! light coming from the layer above (iv-1) goes through iv and onto iv+1. - if (lai_change(L,ft,iv-1) > 0.0_r8)then - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv)* & - lai_change(L,ft,iv-1) / ftweight(L,ft,1) - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv-1)* & - (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) - else - !account fot the light that comes striaght down from unfilled layers above. - tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & - ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) - endif - endif - end do - - !add up all the weighted contributions from the different PFT columns. - weighted_dir_tr(L) = weighted_dir_tr(L) + tr_dir_z(L,ft,currentPatch%nrad(L,ft)+1)*ftweight(L,ft,1) + + endif ! is there vegetation? + + end if ! if the vegetation and zenith filter is active + currentPatch => currentPatch%younger + end do ! Loop linked-list patches + enddo ! Loop Sites + + return + end subroutine ED_Norman_Radiation + - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Sunlit and shaded fraction of leaf layer - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - - !laisum = 0._r8 - do iv = 1,currentPatch%nrad(L,ft) - ! Cumulative leaf area. Original code uses cumulative lai do layer. - ! Now use cumulative lai at center of layer. - ! Same as tr_dir_z calcualtions, but in the middle of the layer? FIX(RF,032414)-WHY? - if (iv == 1) then - laisum = 0.5_r8 * (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv)) - else - laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) - end if - - - if (L == 1)then !top canopy layer - currentPatch%f_sun(L,ft,iv) = exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - else - currentPatch%f_sun(L,ft,iv) = weighted_fsun(L-1)* exp(-k_dir(ft) * laisum)* & - (ftweight(L,ft,iv)/ftweight(L,ft,1)) - endif - - if ( iv > 1 ) then ! becasue we are looking at this layer (not the next) - ! we only ever add fluxes if iv>1 - if (lai_change(L,ft,iv-1) > 0.0_r8)then - currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & - currentPatch%f_sun(L,ft,iv) * & - lai_change(L,ft,iv-1)/ftweight(L,ft,1) - currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & - currentPatch%f_sun(L,ft,iv-1) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) - else - currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & - currentPatch%f_sun(L,ft,iv-1) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - endif - endif - - end do !iv - weighted_fsun(L) = weighted_fsun(L) + currentPatch%f_sun(L,ft,currentPatch%nrad(L,ft))* & - ftweight(L,ft,1) - - ! instance where the first layer ftweight is used a proxy for the whole column. FTWA - ! this is possibly a source of slight error. If we use the ftweight at the top of the PFT column, - ! then we willl underestimate fsun, but if we use ftweight at the bottom of the column, we will - ! underestimate it. Really, we should be tracking the release of direct light from the column as it tapers - ! towards the ground. Is that necessary to get energy closure? It would be quite hard... - endif !present. - end do!pft loop - end do !L - - do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up. - do ft = 1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - !==============================================================================! - ! Iterative solution do scattering - !==============================================================================! + ! ====================================================================================== + + subroutine PatchNormanRadiation (currentPatch, & + albd_parb_out, & ! (ifp,ib) + albi_parb_out, & ! (ifp,ib) + fabd_parb_out, & ! (ifp,ib) + fabi_parb_out, & ! (ifp,ib) + ftdd_parb_out, & ! (ifp,ib) + ftid_parb_out, & ! (ifp,ib) + ftii_parb_out) ! (ifp,ib) + + ! ----------------------------------------------------------------------------------- + ! + ! This routine performs the Norman Radiation scattering for each patch. + ! + ! ----------------------------------------------------------------------------------- + + ! + ! !USES: + use EDPftvarcon , only : EDPftvarcon_inst + use EDtypesMod , only : ed_patch_type + + ! ----------------------------------------------------------------------------------- + ! !ARGUMENTS: + ! ----------------------------------------------------------------------------------- + + type(ed_patch_type), intent(inout), target :: currentPatch + real(r8), intent(inout) :: albd_parb_out(hlm_numSWb) + real(r8), intent(inout) :: albi_parb_out(hlm_numSWb) + real(r8), intent(inout) :: fabd_parb_out(hlm_numSWb) + real(r8), intent(inout) :: fabi_parb_out(hlm_numSWb) + real(r8), intent(inout) :: ftdd_parb_out(hlm_numSWb) + real(r8), intent(inout) :: ftid_parb_out(hlm_numSWb) + real(r8), intent(inout) :: ftii_parb_out(hlm_numSWb) + + ! Locals + ! ----------------------------------------------------------------------------------- + + integer :: radtype, L, ft, j + integer :: iter ! Iteration index + integer :: irep ! Flag to exit iteration loop + 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,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) :: weighted_dir_tr(nclmax) + real(r8) :: weighted_fsun(nclmax) + real(r8) :: weighted_dif_ratio(nclmax,maxSWb) + real(r8) :: weighted_dif_down(nclmax) + real(r8) :: weighted_dif_up(nclmax) + 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(maxpft) ! Radiation transmitted to the soil surface. + real(r8) :: phi2b(maxpft) + real(r8) :: laisum ! cumulative lai+sai for canopy layer (at middle of layer) + real(r8) :: angle + + real(r8),parameter :: tolerance = 0.000000001_r8 + + + 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(nclmax) + + integer :: fp,iv,s ! array indices + integer :: ib ! waveband number + real(r8) :: cosz ! 0.001 <= coszen <= 1.000 + real(r8) :: chil + real(r8) :: gdir + + + real(r8), parameter :: forc_dir(n_rad_stream_types) = (/ 1.0_r8, 0.0_r8 /) ! These are binary switches used + real(r8), parameter :: forc_dif(n_rad_stream_types) = (/ 0.0_r8, 1.0_r8 /) ! to turn off and on radiation streams + + + + associate(& + rhol => EDPftvarcon_inst%rhol , & ! Input: [real(r8) (:) ] leaf reflectance: 1=vis, 2=nir + rhos => EDPftvarcon_inst%rhos , & ! Input: [real(r8) (:) ] stem reflectance: 1=vis, 2=nir + taul => EDPftvarcon_inst%taul , & ! Input: [real(r8) (:) ] leaf transmittance: 1=vis, 2=nir + taus => EDPftvarcon_inst%taus , & ! Input: [real(r8) (:) ] stem transmittance: 1=vis, 2=nir + xl => EDPftvarcon_inst%xl , & ! Input: [real(r8) (:) ] ecophys const - leaf/stem orientation index + clumping_index => EDPftvarcon_inst%clumping_index) + + + ! Initialize local arrays + + weighted_dir_tr(:) = 0._r8 + weighted_dif_down(:) = 0._r8 + weighted_dif_up(:) = 0._r8 + + tr_dir_z(:,:,:) = 0._r8 + tr_dif_z(:,:,:) = 0._r8 + lai_change(:,:,:) = 0._r8 + Dif_up(:,:,:) = 0._r8 + Dif_dn(:,:,:) = 0._r8 + refl_dif(:,:,:,:) = 0.0_r8 + tran_dif(:,:,:,:) = 0.0_r8 + dif_ratio(:,:,:,:) = 0.0_r8 + + + ! Initialize the ouput arrays + ! --------------------------------------------------------------------------------- + albd_parb_out(1:hlm_numSWb) = 0.0_r8 + albi_parb_out(1:hlm_numSWb) = 0.0_r8 + fabd_parb_out(1:hlm_numSWb) = 0.0_r8 + fabi_parb_out(1:hlm_numSWb) = 0.0_r8 + ftdd_parb_out(1:hlm_numSWb) = 1.0_r8 + ftid_parb_out(1:hlm_numSWb) = 1.0_r8 + ftii_parb_out(1:hlm_numSWb) = 1.0_r8 + + ! Is this pft/canopy layer combination present in this patch? + + do L = 1,nclmax + do ft = 1,numpft + currentPatch%canopy_mask(L,ft) = 0 + do iv = 1, currentPatch%nrad(L,ft) + if (currentPatch%canopy_area_profile(L,ft,iv) > 0._r8)then + currentPatch%canopy_mask(L,ft) = 1 + !I think 'present' is only used here... + endif + end do !iv + end do !ft + end do !L + + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Direct beam extinction coefficient, k_dir. PFT specific. + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + cosz = max(0.001_r8, currentPatch%solar_zenith_angle ) !copied from previous radiation code... + do ft = 1,numpft + sb = (90._r8 - (acos(cosz)*180._r8/pi_const)) * (pi_const / 180._r8) + chil = xl(ft) !min(max(xl(ft), -0.4_r8), 0.6_r8 ) + if ( abs(chil) <= 0.01_r8) then + chil = 0.01_r8 + end if + phi1b(ft) = 0.5_r8 - 0.633_r8*chil - 0.330_r8*chil*chil + phi2b(ft) = 0.877_r8 * (1._r8 - 2._r8*phi1b(ft)) !0 = horiz leaves, 1 - vert leaves. + gdir = phi1b(ft) + phi2b(ft) * sin(sb) + !how much direct light penetrates a singleunit of lai? + k_dir(ft) = clumping_index(ft) * gdir / sin(sb) + end do !FT + + + + + !do this once for one unit of diffuse, and once for one unit of direct radiation + do radtype = 1, n_rad_stream_types + + ! Extract information that needs to be provided by ED into local array. + ! RGK: NOT SURE WHY WE NEED FTWEIGHT ... + ! ------------------------------------------------------------------------------ + + ftweight(:,:,:) = 0._r8 + do L = 1,currentPatch%NCL_p + 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) + end do !iv + end do !ft1 + end do !L + if (sum(ftweight(1,:,1))<0.999_r8)then + write(fates_log(),*) 'canopy not full',ftweight(1,:,1) + endif + if (sum(ftweight(1,:,1))>1.0001_r8)then + write(fates_log(),*) 'canopy too full',ftweight(1,:,1) + endif + + do L = 1,currentPatch%NCL_p !start at the top canopy layer (1 is the top layer.) + + weighted_dir_tr(L) = 0.0_r8 + 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 + + if (currentPatch%canopy_mask(L,ft) == 1)then !only do calculation if there are the appropriate leaves. + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Diffuse transmittance, tr_dif, do each layer with thickness elai_z. + ! Estimated do nine sky angles in increments of 10 degrees + ! PFT specific... + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + tr_dif_z(L,ft,:) = 0._r8 + do iv = 1,currentPatch%nrad(L,ft) + do j = 1,9 + angle = (5._r8 + real(j - 1,r8) * 10._r8) * pi_const / 180._r8 + gdir = phi1b(ft) + phi2b(ft) * sin(angle) + tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) + exp(-clumping_index(ft) * & + gdir / sin(angle) * & + (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv))) * & + sin(angle)*cos(angle) + end do + + tr_dif_z(L,ft,iv) = tr_dif_z(L,ft,iv) * 2._r8 * (10._r8 * pi_const / 180._r8) + + end do + + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Direct beam transmittance, tr_dir_z, uses cumulative LAI above layer J to give + ! unscattered direct beam onto layer J. do each PFT section. + ! This is just an decay curve based on k_dir. (leaf & sun angle) + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + if (L==1)then + tr_dir_z(L,ft,1) = 1._r8 + else + tr_dir_z(L,ft,1) = weighted_dir_tr(L-1) + endif + laisum = 0.00_r8 + !total direct beam getting to the bottom of the top canopy. + do iv = 1,currentPatch%nrad(L,ft) + laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) + lai_change(L,ft,iv) = 0.0_r8 + if (( ftweight(L,ft,iv+1) > 0.0_r8 ) .and. ( ftweight(L,ft,iv+1) < ftweight(L,ft,iv) ))then + !where there is a partly empty leaf layer, some fluxes go straight through. + lai_change(L,ft,iv) = ftweight(L,ft,iv)-ftweight(L,ft,iv+1) + endif + if (ftweight(L,ft,iv+1) - ftweight(L,ft,iv) > 1.e-10_r8)then + write(fates_log(),*) 'lower layer has more coverage. This is wrong' , & + ftweight(L,ft,iv),ftweight(L,ft,iv+1),ftweight(L,ft,iv+1)-ftweight(L,ft,iv) + endif + + !n.b. in theory lai_change could be calculated daily in the ED code. + !This is light coming striaght through the canopy. + if (L==1)then + tr_dir_z(L,ft,iv+1) = exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + else + tr_dir_z(L,ft,iv+1) = weighted_dir_tr(L-1)*exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + endif + + if (iv == 1)then + !this is the top layer. + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & + ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) + else + !the lai_change(iv) affects the light incident on layer iv+2 not iv+1 + ! light coming from the layer above (iv-1) goes through iv and onto iv+1. + if (lai_change(L,ft,iv-1) > 0.0_r8)then + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv)* & + lai_change(L,ft,iv-1) / ftweight(L,ft,1) + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv-1)* & + (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) + else + !account fot the light that comes striaght down from unfilled layers above. + tr_dir_z(L,ft,iv+1) = tr_dir_z(L,ft,iv+1) + tr_dir_z(L,ft,iv) * & + ((ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1)) + endif + endif + + end do + + !add up all the weighted contributions from the different PFT columns. + weighted_dir_tr(L) = weighted_dir_tr(L) + tr_dir_z(L,ft,currentPatch%nrad(L,ft)+1)*ftweight(L,ft,1) + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Sunlit and shaded fraction of leaf layer + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - do ib = 1,hlm_numSWb !vis, nir - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Leaf scattering coefficient and terms do diffuse radiation reflected - ! and transmitted by a layer - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - f_not_abs(ft,ib) = rhol(ft,ib) + taul(ft,ib) !leaf level fraction NOT absorbed. - !tr_dif_z is a term that uses the LAI in each layer, whereas rhol and taul do not, - !because they are properties of leaf surfaces and not of the leaf matrix. - do iv = 1,currentPatch%nrad(L,ft) - !How much diffuse light is intercepted and then reflected? - refl_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * rhol(ft,ib) - !How much diffuse light in this layer is transmitted? - tran_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * & - taul(ft,ib) + tr_dif_z(L,ft,iv) - end do + !laisum = 0._r8 + do iv = 1,currentPatch%nrad(L,ft) + ! Cumulative leaf area. Original code uses cumulative lai do layer. + ! Now use cumulative lai at center of layer. + ! Same as tr_dir_z calcualtions, but in the middle of the layer? FIX(RF,032414)-WHY? + if (iv == 1) then + laisum = 0.5_r8 * (currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv)) + else + laisum = laisum + currentPatch%elai_profile(L,ft,iv)+currentPatch%esai_profile(L,ft,iv) + end if + + + if (L == 1)then !top canopy layer + currentPatch%f_sun(L,ft,iv) = exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + else + currentPatch%f_sun(L,ft,iv) = weighted_fsun(L-1)* exp(-k_dir(ft) * laisum)* & + (ftweight(L,ft,iv)/ftweight(L,ft,1)) + endif + + if ( iv > 1 ) then ! becasue we are looking at this layer (not the next) + ! we only ever add fluxes if iv>1 + if (lai_change(L,ft,iv-1) > 0.0_r8)then + currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & + currentPatch%f_sun(L,ft,iv) * & + lai_change(L,ft,iv-1)/ftweight(L,ft,1) + currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & + currentPatch%f_sun(L,ft,iv-1) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ftweight(L,ft,1) + else + currentPatch%f_sun(L,ft,iv) = currentPatch%f_sun(L,ft,iv) + & + currentPatch%f_sun(L,ft,iv-1) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + endif + endif - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Ratio of upward to forward diffuse fluxes, dif_ratio - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! Soil diffuse reflectance (ratio of down to up radiation). - iv = currentPatch%nrad(L,ft) + 1 - if (L == currentPatch%NCL_p)then !nearest the soil - dif_ratio(L,ft,iv,ib) = bc_in(s)%albgr_dif_rb(ib) - else - dif_ratio(L,ft,iv,ib) = weighted_dif_ratio(L+1,ib) - end if - ! Canopy layers, working upwardfrom soil with dif_ratio(iv+1) known - ! FIX(RF,032414) ray tracing eqution - need to find derivation of this... - ! for each unit going down, there are x units going up. - do iv = currentPatch%nrad(L,ft),1, -1 - dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv+1,ib) * & - tran_dif(L,ft,iv,ib)*tran_dif(L,ft,iv,ib) / & - (1._r8 - dif_ratio(L,ft,iv+1,ib) * refl_dif(L,ft,iv,ib)) & - + refl_dif(L,ft,iv,ib) - dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) * & - ftweight(L,ft,iv)/ftweight(L,ft,1) - dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) + dif_ratio(L,ft,iv+1,ib) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - end do - weighted_dif_ratio(L,ib) = weighted_dif_ratio(L,ib) + & - dif_ratio(L,ft,1,ib) * ftweight(L,ft,1) - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - end do!hlm_numSWb - endif ! currentPatch%canopy_mask - end do!ft - end do!L + end do !iv + + weighted_fsun(L) = weighted_fsun(L) + currentPatch%f_sun(L,ft,currentPatch%nrad(L,ft))* & + ftweight(L,ft,1) + + ! instance where the first layer ftweight is used a proxy for the whole column. FTWA + ! this is possibly a source of slight error. If we use the ftweight at the top of the PFT column, + ! then we willl underestimate fsun, but if we use ftweight at the bottom of the column, we will + ! underestimate it. Really, we should be tracking the release of direct light from the column as it tapers + ! towards the ground. Is that necessary to get energy closure? It would be quite hard... + endif !present. + end do!pft loop + end do !L + + + do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up. + do ft = 1,numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + + !==============================================================================! + ! Iterative solution do scattering + !==============================================================================! + + do ib = 1,hlm_numSWb !vis, nir + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Leaf scattering coefficient and terms do diffuse radiation reflected + ! and transmitted by a layer + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + f_not_abs(ft,ib) = rhol(ft,ib) + taul(ft,ib) !leaf level fraction NOT absorbed. + !tr_dif_z is a term that uses the LAI in each layer, whereas rhol and taul do not, + !because they are properties of leaf surfaces and not of the leaf matrix. + do iv = 1,currentPatch%nrad(L,ft) + !How much diffuse light is intercepted and then reflected? + refl_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * rhol(ft,ib) + !How much diffuse light in this layer is transmitted? + tran_dif(L,ft,iv,ib) = (1._r8 - tr_dif_z(L,ft,iv)) * & + taul(ft,ib) + tr_dif_z(L,ft,iv) + end do + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Ratio of upward to forward diffuse fluxes, dif_ratio + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! Soil diffuse reflectance (ratio of down to up radiation). + iv = currentPatch%nrad(L,ft) + 1 + if (L == currentPatch%NCL_p)then !nearest the soil + dif_ratio(L,ft,iv,ib) = currentPatch%gnd_alb_dif(ib) !bc_in(s)%albgr_dif_rb(ib) + else + dif_ratio(L,ft,iv,ib) = weighted_dif_ratio(L+1,ib) + end if + ! Canopy layers, working upwardfrom soil with dif_ratio(iv+1) known + ! FIX(RF,032414) ray tracing eqution - need to find derivation of this... + ! for each unit going down, there are x units going up. + do iv = currentPatch%nrad(L,ft),1, -1 + dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv+1,ib) * & + tran_dif(L,ft,iv,ib)*tran_dif(L,ft,iv,ib) / & + (1._r8 - dif_ratio(L,ft,iv+1,ib) * refl_dif(L,ft,iv,ib)) & + + refl_dif(L,ft,iv,ib) + dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) * & + ftweight(L,ft,iv)/ftweight(L,ft,1) + dif_ratio(L,ft,iv,ib) = dif_ratio(L,ft,iv,ib) + dif_ratio(L,ft,iv+1,ib) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + end do + weighted_dif_ratio(L,ib) = weighted_dif_ratio(L,ib) + & + dif_ratio(L,ft,1,ib) * ftweight(L,ft,1) + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + end do!hlm_numSWb + endif ! currentPatch%canopy_mask + end do!ft + end do!L - do ib = 1,hlm_numSWb - Dif_dn(:,:,:) = 0.00_r8 - 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 - if (currentPatch%canopy_mask(L,ft) == 1)then - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! First estimates do downward and upward diffuse flux - ! - ! Dif_dn = forward diffuse flux onto layer J - ! Dif_up = Upward diffuse flux above layer J - ! - ! Solved here without direct beam radiation and using dif_ratio = Dif_up / Dif_dn - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! downward diffuse flux onto the top surface of the canopy - - if (L == 1)then - Dif_dn(L,ft,1) = forc_dif(ifp,ib) - else - Dif_dn(L,ft,1) = weighted_dif_down(L-1) - end if - ! forward diffuse flux within the canopy and at soil, working forward through canopy - do iv = 1,currentPatch%nrad(L,ft) - denom = refl_dif(L,ft,iv,ib) * dif_ratio(L,ft,iv,ib) - denom = 1._r8 - denom - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) / & - denom *ftweight(L,ft,iv)/ftweight(L,ft,1) - if (iv > 1)then - if (lai_change(L,ft,iv-1) > 0.0_r8)then - !here we are thinking about whether the layer above had an laichange, - !but calculating the flux onto the layer below. - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv)* & - lai_change(L,ft,iv-1)/ftweight(L,ft,1) - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv-1)* & - (ftweight(L,ft,1)-ftweight(L,ft,iv-1)/ftweight(L,ft,1)) - else - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - endif - else - Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - endif - end do - - weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & - ftweight(L,ft,1) - - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - endif !present - end do !ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is the the (incomplete) understorey? - !Add on the radiation going through the canopy gaps. - weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1)*(1.0-sum(ftweight(L,:,1))) - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - endif - end do !L - - do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom. - weighted_dif_up(L) = 0._r8 - do ft = 1, numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - !Bounce diffuse radiation off soil surface. - iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then !is this the bottom layer ? - Dif_up(L,ft,iv) =bc_in(s)%albgr_dif_rb(ib) * Dif_dn(L,ft,iv) - else - Dif_up(L,ft,iv) = weighted_dif_up(L+1) - end if - ! Upward diffuse flux within the canopy and above the canopy, working upward through canopy - - do iv = currentPatch%nrad(L,ft), 1, -1 - if (lai_change(L,ft,iv) > 0.0_r8)then - Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * & - ftweight(L,ft,iv) / ftweight(L,ft,1) - Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & - tran_dif(L,ft,iv,ib) * lai_change(L,ft,iv)/ftweight(L,ft,1) - Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & - (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - !nb is this the right constuction? - ! the radiation that hits the empty space is not reflected. - else - Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * ftweight(L,ft,iv) - Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * (1.0_r8-ftweight(L,ft,iv)) - endif - end do + + do ib = 1,hlm_numSWb + Dif_dn(:,:,:) = 0.00_r8 + 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 + if (currentPatch%canopy_mask(L,ft) == 1)then + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! First estimates do downward and upward diffuse flux + ! + ! Dif_dn = forward diffuse flux onto layer J + ! Dif_up = Upward diffuse flux above layer J + ! + ! Solved here without direct beam radiation and using dif_ratio = Dif_up / Dif_dn + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! downward diffuse flux onto the top surface of the canopy + + if (L == 1)then + Dif_dn(L,ft,1) = forc_dif(radtype) + else + Dif_dn(L,ft,1) = weighted_dif_down(L-1) + end if + ! forward diffuse flux within the canopy and at soil, working forward through canopy + do iv = 1,currentPatch%nrad(L,ft) + denom = refl_dif(L,ft,iv,ib) * dif_ratio(L,ft,iv,ib) + denom = 1._r8 - denom + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) / & + denom *ftweight(L,ft,iv)/ftweight(L,ft,1) + if (iv > 1)then + if (lai_change(L,ft,iv-1) > 0.0_r8)then + !here we are thinking about whether the layer above had an laichange, + !but calculating the flux onto the layer below. + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv)* & + lai_change(L,ft,iv-1)/ftweight(L,ft,1) + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1)+ Dif_dn(L,ft,iv-1)* & + (ftweight(L,ft,1)-ftweight(L,ft,iv-1)/ftweight(L,ft,1)) + else + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + endif + else + Dif_dn(L,ft,iv+1) = Dif_dn(L,ft,iv+1) + Dif_dn(L,ft,iv) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + endif + end do + + weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & + ftweight(L,ft,1) - weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) - !instance where the first layer ftweight is used a proxy for the whole column. FTWA - endif !present - end do !ft - 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: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:numpft,1)))*bc_in(s)%albgr_dir_rb(ib) - endif - end do !L - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - ! 3. Iterative calculation of forward and upward diffuse fluxes, iNCL_puding - ! scattered direct beam - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - - ! Flag to exit iteration loop: 0 = exit and 1 = iterate - irep = 1 - ! Iteration loop - iter = 0 - do while(irep ==1 .and. iter<50) - - iter = iter + 1 - irep = 0 - do L = 1,currentPatch%NCL_p !working from the top down - weighted_dif_down(L) = 0._r8 - do ft =1,numpft - if (currentPatch%canopy_mask(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. - ! Note: down = forward flux onto next layer - if (L == 1)then !is this the top layer? - Dif_dn(L,ft,1) = forc_dif(ifp,ib) - else - Dif_dn(L,ft,1) = weighted_dif_down(L-1) - end if - down_rad = 0._r8 - - do iv = 1, currentPatch%nrad(L,ft) - - down_rad = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) + & - Dif_up(L,ft,iv+1) * refl_dif(L,ft,iv,ib) + & - forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) * (1.00_r8 - & - exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & - currentPatch%esai_profile(L,ft,iv)))) * taul(ft,ib) - down_rad = down_rad *(ftweight(L,ft,iv)/ftweight(L,ft,1)) - - if (iv > 1)then - if (lai_change(L,ft,iv-1) > 0.0_r8)then - down_rad = down_rad + Dif_dn(L,ft,iv) * lai_change(L,ft,iv-1)/ftweight(L,ft,1) - down_rad = down_rad + Dif_dn(L,ft,iv-1) * (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ & - ftweight(L,ft,1) - else - down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & - ftweight(L,ft,1) - endif - else - down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & - ftweight(L,ft,1) - endif - - !this is just Dif down, plus refl up, plus dir intercepted and turned into dif... , - if (abs(down_rad - Dif_dn(L,ft,iv+1)) > tolerance)then - irep = 1 - end if - Dif_dn(L,ft,iv+1) = down_rad - - end do !iv - - weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & - ftweight(L,ft,1) - - 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: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 - if (currentPatch%canopy_mask(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 - if (L==currentPatch%NCL_p)then !In the bottom canopy layer, reflect off the soil - Dif_up(L,ft,iv) = Dif_dn(L,ft,iv) *bc_in(s)%albgr_dif_rb(ib) + & - forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) *bc_in(s)%albgr_dir_rb(ib) - else !In the other canopy layers, reflect off the underlying vegetation. - Dif_up(L,ft,iv) = weighted_dif_up(L+1) - end if - - ! Upward diffuse flux within and above the canopy, working upward through canopy - ! with Dif_dn from previous interation. Note: up = upward flux above current layer - do iv = currentPatch%nrad(L,ft),1,-1 - !this is radiation up, by layer transmittance, by - - !reflection of the lower layer, - up_rad = Dif_dn(L,ft,iv) * refl_dif(L,ft,iv,ib) - up_rad = up_rad + forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) * (1.00_r8 - exp(-k_dir(ft) * & - (currentPatch%elai_profile(L,ft,iv) + currentPatch%esai_profile(L,ft,iv)))) * & - rhol(ft,ib) - up_rad = up_rad + Dif_up(L,ft,iv+1) * tran_dif(L,ft,iv,ib) - up_rad = up_rad * ftweight(L,ft,iv)/ftweight(L,ft,1) - up_rad = up_rad + Dif_up(L,ft,iv+1) *(ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) - ! THE LOWER LAYER FLUX IS HOMOGENIZED, SO WE DON"T CONSIDER THE LAI_CHANGE HERE... + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + endif !present + end do !ft + if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is the the (incomplete) understorey? + !Add on the radiation going through the canopy gaps. + weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1)*(1.0-sum(ftweight(L,:,1))) + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + endif + end do !L + + do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom. + weighted_dif_up(L) = 0._r8 + do ft = 1, numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + !Bounce diffuse radiation off soil surface. + iv = currentPatch%nrad(L,ft) + 1 + if (L==currentPatch%NCL_p)then !is this the bottom layer ? + Dif_up(L,ft,iv) = currentPatch%gnd_alb_dif(ib) * Dif_dn(L,ft,iv) + else + Dif_up(L,ft,iv) = weighted_dif_up(L+1) + end if + ! Upward diffuse flux within the canopy and above the canopy, working upward through canopy + + do iv = currentPatch%nrad(L,ft), 1, -1 + if (lai_change(L,ft,iv) > 0.0_r8)then + Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * & + ftweight(L,ft,iv) / ftweight(L,ft,1) + Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & + tran_dif(L,ft,iv,ib) * lai_change(L,ft,iv)/ftweight(L,ft,1) + Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * & + (ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + !nb is this the right constuction? + ! the radiation that hits the empty space is not reflected. + else + Dif_up(L,ft,iv) = dif_ratio(L,ft,iv,ib) * Dif_dn(L,ft,iv) * ftweight(L,ft,iv) + Dif_up(L,ft,iv) = Dif_up(L,ft,iv) + Dif_up(L,ft,iv+1) * (1.0_r8-ftweight(L,ft,iv)) + endif + end do + + weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) + !instance where the first layer ftweight is used a proxy for the whole column. FTWA + endif !present + end do !ft + 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_r8-sum(ftweight(L,1:numpft,1))) * & + weighted_dif_down(L-1) * currentPatch%gnd_alb_dif(ib) + !direct to diffuse + weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(radtype) * & + weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) * currentPatch%gnd_alb_dir(ib) + endif + end do !L + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + ! 3. Iterative calculation of forward and upward diffuse fluxes, iNCL_puding + ! scattered direct beam + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + + ! Flag to exit iteration loop: 0 = exit and 1 = iterate + irep = 1 + ! Iteration loop + iter = 0 + do while(irep ==1 .and. iter<50) + + iter = iter + 1 + irep = 0 + do L = 1,currentPatch%NCL_p !working from the top down + weighted_dif_down(L) = 0._r8 + do ft =1,numpft + if (currentPatch%canopy_mask(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. + ! Note: down = forward flux onto next layer + if (L == 1)then !is this the top layer? + Dif_dn(L,ft,1) = forc_dif(radtype) + else + Dif_dn(L,ft,1) = weighted_dif_down(L-1) + end if + down_rad = 0._r8 + + do iv = 1, currentPatch%nrad(L,ft) + + down_rad = Dif_dn(L,ft,iv) * tran_dif(L,ft,iv,ib) + & + Dif_up(L,ft,iv+1) * refl_dif(L,ft,iv,ib) + & + forc_dir(radtype) * tr_dir_z(L,ft,iv) * (1.00_r8 - & + exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & + currentPatch%esai_profile(L,ft,iv)))) * taul(ft,ib) + down_rad = down_rad *(ftweight(L,ft,iv)/ftweight(L,ft,1)) - if (abs(up_rad - Dif_up(L,ft,iv)) > tolerance) then !are we close to the tolerance level? - irep = 1 - end if - Dif_up(L,ft,iv) = up_rad + if (iv > 1)then + if (lai_change(L,ft,iv-1) > 0.0_r8)then + down_rad = down_rad + Dif_dn(L,ft,iv) * lai_change(L,ft,iv-1)/ftweight(L,ft,1) + down_rad = down_rad + Dif_dn(L,ft,iv-1) * (ftweight(L,ft,1)-ftweight(L,ft,iv-1))/ & + ftweight(L,ft,1) + else + down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & + ftweight(L,ft,1) + endif + else + down_rad = down_rad + Dif_dn(L,ft,iv) * (ftweight(L,ft,1)-ftweight(L,ft,iv))/ & + ftweight(L,ft,1) + endif - end do !iv - weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) - end if !present - end do!ft - - 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: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:numpft,1)))*bc_in(s)%albgr_dir_rb(ib) - end if - end do!L - end do ! do while over iter - - abs_rad(ib) = 0._r8 - tr_soili = 0._r8 - tr_soild = 0._r8 - - 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 - if (currentPatch%canopy_mask(L,ft) == 1)then - !==============================================================================! - ! Compute absorbed flux densities - !==============================================================================! - - ! Absorbed direct beam and diffuse do leaf layers - do iv = 1, currentPatch%nrad(L,ft) - Abs_dir_z(ft,iv) = ftweight(L,ft,iv)* forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) * & - (1.00_r8 - exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & - currentPatch%esai_profile(L,ft,iv)))) * (1.00_r8 - f_not_abs(ft,ib)) - Abs_dif_z(ft,iv) = ftweight(L,ft,iv)* ((Dif_dn(L,ft,iv) + & - Dif_up(L,ft,iv+1)) * (1.00_r8 - tr_dif_z(L,ft,iv)) * & - (1.00_r8 - f_not_abs(ft,ib))) - end do - - ! Absorbed direct beam and diffuse do soil - if (L == currentPatch%NCL_p)then - iv = currentPatch%nrad(L,ft) + 1 - Abs_dif_z(ft,iv) = ftweight(L,ft,1)*Dif_dn(L,ft,iv) * (1.0_r8 -bc_in(s)%albgr_dif_rb(ib)) - Abs_dir_z(ft,iv) = ftweight(L,ft,1)*forc_dir(ifp,ib) * & - tr_dir_z(L,ft,iv) * (1.0_r8 -bc_in(s)%albgr_dir_rb(ib)) - tr_soild = tr_soild + ftweight(L,ft,1)*forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) - tr_soili = tr_soili + ftweight(L,ft,1)*Dif_dn(L,ft,iv) - end if - ! Absorbed radiation, shaded and sunlit portions of leaf layers - !here we get one unit of diffuse radiation... how much of - !it is absorbed? - if (ib == ivis) then ! only set the absorbed PAR for the visible light band. - do iv = 1, currentPatch%nrad(L,ft) - if (radtype==idirect) then - if ( debug ) then - write(fates_log(),*) 'EDsurfAlb 730 ',Abs_dif_z(ft,iv),currentPatch%f_sun(L,ft,iv) - write(fates_log(),*) 'EDsurfAlb 731 ', currentPatch%fabd_sha_z(L,ft,iv), & - currentPatch%fabd_sun_z(L,ft,iv) - endif - currentPatch%fabd_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - (1._r8 - currentPatch%f_sun(L,ft,iv)) - currentPatch%fabd_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - currentPatch%f_sun(L,ft,iv) + & - Abs_dir_z(ft,iv) - else - currentPatch%fabi_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - (1._r8 - currentPatch%f_sun(L,ft,iv)) - currentPatch%fabi_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & - currentPatch%f_sun(L,ft,iv) - endif - if ( debug ) then - write(fates_log(),*) 'EDsurfAlb 740 ', currentPatch%fabd_sha_z(L,ft,iv), & - currentPatch%fabd_sun_z(L,ft,iv) - endif - end do - endif ! ib - - !==============================================================================! - ! Sum fluxes - !==============================================================================! - ! Solar radiation absorbed by ground - iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then - abs_rad(ib) = abs_rad(ib) + (Abs_dir_z(ft,iv) + Abs_dif_z(ft,iv)) - end if - ! Solar radiation absorbed by vegetation and sunlit/shaded leaves - do iv = 1,currentPatch%nrad(L,ft) - if (radtype == idirect)then - currentPatch%fabd(ib) = currentPatch%fabd(ib) + & - Abs_dir_z(ft,iv)+Abs_dif_z(ft,iv) - ! bc_out(s)%fabd_parb(ifp,ib) = currentPatch%fabd(ib) - else - currentPatch%fabi(ib) = currentPatch%fabi(ib) + Abs_dif_z(ft,iv) - ! bc_out(s)%fabi_parb(ifp,ib) = currentPatch%fabi(ib) - endif - end do - ! Albefor - if (L==1)then !top canopy layer. - if (radtype == idirect)then - bc_out(s)%albd_parb(ifp,ib) = bc_out(s)%albd_parb(ifp,ib) + & - Dif_up(L,ft,1) * ftweight(L,ft,1) - else - bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + & - Dif_up(L,ft,1) * ftweight(L,ft,1) - end if - end if - - ! pass normalized PAR profiles for use in diagnostic averaging for history fields - if (ib == ivis) then ! only diagnose PAR profiles for the visible band - do iv = 1, currentPatch%nrad(L,ft) - currentPatch%nrmlzd_parprof_pft_dir_z(radtype,L,ft,iv) = & - forc_dir(ifp,ib) * tr_dir_z(L,ft,iv) - currentPatch%nrmlzd_parprof_pft_dif_z(radtype,L,ft,iv) = & - Dif_dn(L,ft,iv) + Dif_up(L,ft,iv) - ! - currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) = & - currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) + & - (forc_dir(ifp,ib) * tr_dir_z(L,ft,iv)) * & - (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) - currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) = & - currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) + & - (Dif_dn(L,ft,iv) + Dif_up(L,ft,iv)) * & - (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) - end do - end if ! ib = visible - end if ! present - end do !ft - if (radtype == idirect)then - bc_out(s)%fabd_parb(ifp,ib) = currentPatch%fabd(ib) - else - bc_out(s)%fabi_parb(ifp,ib) = currentPatch%fabi(ib) - endif - - - !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: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: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 == idirect)then - currentPatch%tr_soil_dir(ib) = tr_soild - currentPatch%tr_soil_dir_dif(ib) = tr_soili - currentPatch%sabs_dir(ib) = abs_rad(ib) - bc_out(s)%ftdd_parb(ifp,ib) = tr_soild - bc_out(s)%ftid_parb(ifp,ib) = tr_soili - else - currentPatch%tr_soil_dif(ib) = tr_soili - currentPatch%sabs_dif(ib) = abs_rad(ib) - bc_out(s)%ftii_parb(ifp,ib) = tr_soili - end if - - end do!l - - - !==============================================================================! - ! Conservation check - !==============================================================================! - ! Total radiation balance: absorbed = incoming - outgoing - - if (radtype == idirect)then - error = abs(currentPatch%sabs_dir(ib) - (currentPatch%tr_soil_dir(ib) * & - (1.0_r8-bc_in(s)%albgr_dir_rb(ib)) + & - currentPatch%tr_soil_dir_dif(ib) * (1.0_r8-bc_in(s)%albgr_dif_rb(ib)))) - 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: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)) - - do ft =1,3 - iv = currentPatch%nrad(1,ft) + 1 - write(fates_log(),*) 'abs soil fluxes', Abs_dir_z(ft,iv),Abs_dif_z(ft,iv) - end do - - end if - else - if ( abs(currentPatch%sabs_dif(ib)-(currentPatch%tr_soil_dif(ib) * & - (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:numpft,1)) - endif - endif - - if (radtype == idirect)then - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabd_parb(ifp,ib) + bc_out(s)%albd_parb(ifp,ib) + currentPatch%sabs_dir(ib)) - else - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabi_parb(ifp,ib) + bc_out(s)%albi_parb(ifp,ib) + currentPatch%sabs_dif(ib)) - endif - lai_reduction(:) = 0.0_r8 - do L = 1, currentPatch%NCL_p - do ft =1,numpft - if (currentPatch%canopy_mask(L,ft) == 1)then - do iv = 1, currentPatch%nrad(L,ft) - if (lai_change(L,ft,iv) > 0.0_r8)then - lai_reduction(L) = max(lai_reduction(L),lai_change(L,ft,iv)) - endif - enddo - endif - enddo - enddo + !this is just Dif down, plus refl up, plus dir intercepted and turned into dif... , + if (abs(down_rad - Dif_dn(L,ft,iv+1)) > tolerance)then + irep = 1 + end if + Dif_dn(L,ft,iv+1) = down_rad + + end do !iv + + weighted_dif_down(L) = weighted_dif_down(L) + Dif_dn(L,ft,currentPatch%nrad(L,ft)+1) * & + ftweight(L,ft,1) + + 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_r8-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 + if (currentPatch%canopy_mask(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 + if (L==currentPatch%NCL_p)then !In the bottom canopy layer, reflect off the soil + Dif_up(L,ft,iv) = Dif_dn(L,ft,iv) * currentPatch%gnd_alb_dif(ib) + & + forc_dir(radtype) * tr_dir_z(L,ft,iv) * currentPatch%gnd_alb_dir(ib) + else !In the other canopy layers, reflect off the underlying vegetation. + Dif_up(L,ft,iv) = weighted_dif_up(L+1) + end if + + ! Upward diffuse flux within and above the canopy, working upward through canopy + ! with Dif_dn from previous interation. Note: up = upward flux above current layer + do iv = currentPatch%nrad(L,ft),1,-1 + !this is radiation up, by layer transmittance, by + + !reflection of the lower layer, + up_rad = Dif_dn(L,ft,iv) * refl_dif(L,ft,iv,ib) + up_rad = up_rad + forc_dir(radtype) * tr_dir_z(L,ft,iv) * (1.00_r8 - exp(-k_dir(ft) * & + (currentPatch%elai_profile(L,ft,iv) + currentPatch%esai_profile(L,ft,iv)))) * & + rhol(ft,ib) + up_rad = up_rad + Dif_up(L,ft,iv+1) * tran_dif(L,ft,iv,ib) + up_rad = up_rad * ftweight(L,ft,iv)/ftweight(L,ft,1) + up_rad = up_rad + Dif_up(L,ft,iv+1) *(ftweight(L,ft,1)-ftweight(L,ft,iv))/ftweight(L,ft,1) + ! THE LOWER LAYER FLUX IS HOMOGENIZED, SO WE DON"T CONSIDER THE LAI_CHANGE HERE... + + if (abs(up_rad - Dif_up(L,ft,iv)) > tolerance) then !are we close to the tolerance level? + irep = 1 + end if + Dif_up(L,ft,iv) = up_rad + + end do !iv + weighted_dif_up(L) = weighted_dif_up(L) + Dif_up(L,ft,1) * ftweight(L,ft,1) + end if !present + end do!ft + + 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:numpft,1))) * & + weighted_dif_down(L-1) * currentPatch%gnd_alb_dif(ib) + weighted_dif_up(L) = weighted_dif_up(L) + forc_dir(radtype) * & + weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1)))*currentPatch%gnd_alb_dir(ib) + end if + end do!L + end do ! do while over iter - if (radtype == idirect)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. - if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then - bc_out(s)%albd_parb(ifp,ib) = bc_out(s)%albd_parb(ifp,ib) + error - !this terms adds the error back on to the albedo. While this is partly inexcusable, it is - ! in the medium term a solution that - ! prevents the model from crashing with small and occasional energy balances issues. - ! These are extremely difficult to debug, many have been solved already, leading - ! to the complexity of this code, but where the system generates occasional errors, we - ! will deal with them for now. - end if - if (abs(error) > 0.15_r8)then - 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,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) - - bc_out(s)%albd_parb(ifp,ib) = bc_out(s)%albd_parb(ifp,ib) + error - end if - else - - if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then - bc_out(s)%albi_parb(ifp,ib) = bc_out(s)%albi_parb(ifp,ib) + error - end if - - if (abs(error) > 0.15_r8)then - 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,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,:) - write(fates_log(),*) 'ftw',sum(ftweight(1,1:numpft,1)),ftweight(1,1:numpft,1) - write(fates_log(),*) 'present',currentPatch%canopy_mask(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 - - if (radtype == idirect)then - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabd_parb(ifp,ib) + bc_out(s)%albd_parb(ifp,ib) + currentPatch%sabs_dir(ib)) - else - error = (forc_dir(ifp,ib) + forc_dif(ifp,ib)) - & - (bc_out(s)%fabi_parb(ifp,ib) + bc_out(s)%albi_parb(ifp,ib) + currentPatch%sabs_dif(ib)) - endif - - if (abs(error) > 0.00000001_r8)then - write(fates_log(),*) 'there is still error after correction',error ,ifp,ib - end if - - end if + abs_rad(ib) = 0._r8 + tr_soili = 0._r8 + tr_soild = 0._r8 + + 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 + if (currentPatch%canopy_mask(L,ft) == 1)then + !==============================================================================! + ! Compute absorbed flux densities + !==============================================================================! + + ! Absorbed direct beam and diffuse do leaf layers + do iv = 1, currentPatch%nrad(L,ft) + Abs_dir_z(ft,iv) = ftweight(L,ft,iv)* forc_dir(radtype) * tr_dir_z(L,ft,iv) * & + (1.00_r8 - exp(-k_dir(ft) * (currentPatch%elai_profile(L,ft,iv)+ & + currentPatch%esai_profile(L,ft,iv)))) * (1.00_r8 - f_not_abs(ft,ib)) + Abs_dif_z(ft,iv) = ftweight(L,ft,iv)* ((Dif_dn(L,ft,iv) + & + Dif_up(L,ft,iv+1)) * (1.00_r8 - tr_dif_z(L,ft,iv)) * & + (1.00_r8 - f_not_abs(ft,ib))) + end do + + ! Absorbed direct beam and diffuse do soil + if (L == currentPatch%NCL_p)then + iv = currentPatch%nrad(L,ft) + 1 + Abs_dif_z(ft,iv) = ftweight(L,ft,1)*Dif_dn(L,ft,iv) * (1.0_r8 - currentPatch%gnd_alb_dif(ib) ) + Abs_dir_z(ft,iv) = ftweight(L,ft,1)*forc_dir(radtype) * & + tr_dir_z(L,ft,iv) * (1.0_r8 - currentPatch%gnd_alb_dir(ib) ) + tr_soild = tr_soild + ftweight(L,ft,1)*forc_dir(radtype) * tr_dir_z(L,ft,iv) + tr_soili = tr_soili + ftweight(L,ft,1)*Dif_dn(L,ft,iv) + end if + + ! Absorbed radiation, shaded and sunlit portions of leaf layers + !here we get one unit of diffuse radiation... how much of + !it is absorbed? + if (ib == ivis) then ! only set the absorbed PAR for the visible light band. + do iv = 1, currentPatch%nrad(L,ft) + if (radtype==idirect) then + if ( debug ) then + write(fates_log(),*) 'EDsurfAlb 730 ',Abs_dif_z(ft,iv),currentPatch%f_sun(L,ft,iv) + write(fates_log(),*) 'EDsurfAlb 731 ', currentPatch%fabd_sha_z(L,ft,iv), & + currentPatch%fabd_sun_z(L,ft,iv) + endif + currentPatch%fabd_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + (1._r8 - currentPatch%f_sun(L,ft,iv)) + currentPatch%fabd_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + currentPatch%f_sun(L,ft,iv) + & + Abs_dir_z(ft,iv) + else + currentPatch%fabi_sha_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + (1._r8 - currentPatch%f_sun(L,ft,iv)) + currentPatch%fabi_sun_z(L,ft,iv) = Abs_dif_z(ft,iv) * & + currentPatch%f_sun(L,ft,iv) + endif + if ( debug ) then + write(fates_log(),*) 'EDsurfAlb 740 ', currentPatch%fabd_sha_z(L,ft,iv), & + currentPatch%fabd_sun_z(L,ft,iv) + endif + end do + endif ! ib + + + !==============================================================================! + ! Sum fluxes + !==============================================================================! + ! Solar radiation absorbed by ground + iv = currentPatch%nrad(L,ft) + 1 + if (L==currentPatch%NCL_p)then + abs_rad(ib) = abs_rad(ib) + (Abs_dir_z(ft,iv) + Abs_dif_z(ft,iv)) + end if + ! Solar radiation absorbed by vegetation and sunlit/shaded leaves + do iv = 1,currentPatch%nrad(L,ft) + if (radtype == idirect)then + currentPatch%fabd(ib) = currentPatch%fabd(ib) + & + Abs_dir_z(ft,iv)+Abs_dif_z(ft,iv) + ! bc_out(s)%fabd_parb_out(ib) = currentPatch%fabd(ib) + else + currentPatch%fabi(ib) = currentPatch%fabi(ib) + Abs_dif_z(ft,iv) + ! bc_out(s)%fabi_parb_out(ib) = currentPatch%fabi(ib) + endif + end do + + ! Albefor + if (L==1)then !top canopy layer. + if (radtype == idirect)then + albd_parb_out(ib) = albd_parb_out(ib) + & + Dif_up(L,ft,1) * ftweight(L,ft,1) + else + albi_parb_out(ib) = albi_parb_out(ib) + & + Dif_up(L,ft,1) * ftweight(L,ft,1) + end if + end if + + ! pass normalized PAR profiles for use in diagnostic averaging for history fields + if (ib == ivis) then ! only diagnose PAR profiles for the visible band + do iv = 1, currentPatch%nrad(L,ft) + currentPatch%nrmlzd_parprof_pft_dir_z(radtype,L,ft,iv) = & + forc_dir(radtype) * tr_dir_z(L,ft,iv) + currentPatch%nrmlzd_parprof_pft_dif_z(radtype,L,ft,iv) = & + Dif_dn(L,ft,iv) + Dif_up(L,ft,iv) + ! + currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) = & + currentPatch%nrmlzd_parprof_dir_z(radtype,L,iv) + & + (forc_dir(radtype) * tr_dir_z(L,ft,iv)) * & + (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) + currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) = & + currentPatch%nrmlzd_parprof_dif_z(radtype,L,iv) + & + (Dif_dn(L,ft,iv) + Dif_up(L,ft,iv)) * & + (ftweight(L,ft,iv) / sum(ftweight(L,1:numpft,iv))) + end do + end if ! ib = visible + end if ! present + end do !ft + if (radtype == idirect)then + fabd_parb_out(ib) = currentPatch%fabd(ib) + else + fabi_parb_out(ib) = currentPatch%fabi(ib) + endif + + + !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:numpft,1)))*(1.0_r8-currentPatch%gnd_alb_dif(ib) ) + abs_rad(ib) = abs_rad(ib) + forc_dir(radtype) * weighted_dir_tr(L-1) * & + (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-currentPatch%gnd_alb_dir(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(radtype) * weighted_dir_tr(L-1) * (1.0_r8-sum(ftweight(L,1:numpft,1))) + endif + + if (radtype == idirect)then + currentPatch%tr_soil_dir(ib) = tr_soild + currentPatch%tr_soil_dir_dif(ib) = tr_soili + currentPatch%sabs_dir(ib) = abs_rad(ib) + ftdd_parb_out(ib) = tr_soild + ftid_parb_out(ib) = tr_soili + else + currentPatch%tr_soil_dif(ib) = tr_soili + currentPatch%sabs_dif(ib) = abs_rad(ib) + ftii_parb_out(ib) = tr_soili + end if + + end do!l + - end do !hlm_numSWb - - enddo ! rad-type - endif ! is there vegetation? + !==============================================================================! + ! Conservation check + !==============================================================================! + ! Total radiation balance: absorbed = incoming - outgoing + + if (radtype == idirect)then + error = abs(currentPatch%sabs_dir(ib) - (currentPatch%tr_soil_dir(ib) * & + (1.0_r8-currentPatch%gnd_alb_dir(ib) ) + & + currentPatch%tr_soil_dir_dif(ib) * (1.0_r8-currentPatch%gnd_alb_dif(ib) ))) + if ( abs(error) > 0.0001)then + write(fates_log(),*)'dir ground absorption error',error,currentPatch%sabs_dir(ib), & + currentPatch%tr_soil_dir(ib)* & + (1.0_r8-currentPatch%gnd_alb_dir(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-currentPatch%gnd_alb_dir(ib) ) + + do ft =1,3 + iv = currentPatch%nrad(1,ft) + 1 + write(fates_log(),*) 'abs soil fluxes', Abs_dir_z(ft,iv),Abs_dif_z(ft,iv) + end do + + end if + else + if ( abs(currentPatch%sabs_dif(ib)-(currentPatch%tr_soil_dif(ib) * & + (1.0_r8-currentPatch%gnd_alb_dif(ib) ))) > 0.0001_r8)then + write(fates_log(),*)'dif ground absorption error',currentPatch%sabs_dif(ib) , & + (currentPatch%tr_soil_dif(ib)* & + (1.0_r8-currentPatch%gnd_alb_dif(ib) )),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) + endif + endif + + if (radtype == idirect)then + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabd_parb_out(ib) + albd_parb_out(ib) + currentPatch%sabs_dir(ib)) + else + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabi_parb_out(ib) + albi_parb_out(ib) + currentPatch%sabs_dif(ib)) + endif + lai_reduction(:) = 0.0_r8 + do L = 1, currentPatch%NCL_p + do ft =1,numpft + if (currentPatch%canopy_mask(L,ft) == 1)then + do iv = 1, currentPatch%nrad(L,ft) + if (lai_change(L,ft,iv) > 0.0_r8)then + lai_reduction(L) = max(lai_reduction(L),lai_change(L,ft,iv)) + endif + enddo + endif + enddo + enddo + + if (radtype == idirect)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. + if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then + albd_parb_out(ib) = albd_parb_out(ib) + error + !this terms adds the error back on to the albedo. While this is partly inexcusable, it is + ! in the medium term a solution that + ! prevents the model from crashing with small and occasional energy balances issues. + ! These are extremely difficult to debug, many have been solved already, leading + ! to the complexity of this code, but where the system generates occasional errors, we + ! will deal with them for now. + end if + if (abs(error) > 0.15_r8)then + write(fates_log(),*) 'Large Dir Radn consvn error',error ,ib + write(fates_log(),*) 'diags', albd_parb_out(ib), ftdd_parb_out(ib), & + ftid_parb_out(ib), fabd_parb_out(ib) + 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(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib) + + albd_parb_out(ib) = albd_parb_out(ib) + error + end if + else + + if (abs(error) > 1.e-9_r8 .and. abs(error) < 0.15_r8)then + albi_parb_out(ib) = albi_parb_out(ib) + error + end if + + if (abs(error) > 0.15_r8)then + write(fates_log(),*) '>5% Dif Radn consvn error',error ,ib + write(fates_log(),*) 'diags', albi_parb_out(ib), ftii_parb_out(ib), & + fabi_parb_out(ib) + 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(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib) + 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%canopy_mask(1,1:numpft) + write(fates_log(),*) 'CAP',currentPatch%canopy_area_profile(1,1:numpft,1) + + albi_parb_out(ib) = albi_parb_out(ib) + error + end if + + if (radtype == idirect)then + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabd_parb_out(ib) + albd_parb_out(ib) + currentPatch%sabs_dir(ib)) + else + error = (forc_dir(radtype) + forc_dif(radtype)) - & + (fabi_parb_out(ib) + albi_parb_out(ib) + currentPatch%sabs_dif(ib)) + endif + + if (abs(error) > 0.00000001_r8)then + write(fates_log(),*) 'there is still error after correction',error ,ib + end if + + end if + + end do !hlm_numSWb + + enddo ! rad-type + + + end associate + return + end subroutine PatchNormanRadiation - end if ! if the vegetation and zenith filter is active - currentPatch => currentPatch%younger - end do ! Loop linked-list patches - enddo ! Loop Sites - - end associate - return - end subroutine ED_Norman_Radiation - ! ====================================================================================== subroutine ED_SunShadeFracs(nsites, sites,bc_in,bc_out) diff --git a/biogeophys/FatesBstressMod.F90 b/biogeophys/FatesBstressMod.F90 new file mode 100644 index 0000000000..88f7d44710 --- /dev/null +++ b/biogeophys/FatesBstressMod.F90 @@ -0,0 +1,99 @@ +module FatesBstressMod + + !------------------------------------------------------------------------------------- + ! Description: calculate the stress impact on transpiration from salinity and sulphide in soils + ! note that water stress is calculated in EDBtranMod or HYDRO + ! ------------------------------------------------------------------------------------ + ! + use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm + use FatesConstantsMod , only : itrue,ifalse + use EDTypesMod , only : ed_site_type, & + ed_patch_type, & + ed_cohort_type, & + maxpft + use shr_kind_mod , only : r8 => shr_kind_r8 + use FatesInterfaceMod , only : bc_in_type, & + bc_out_type, & + numpft + use FatesInterfaceMod , only : hlm_use_planthydro + use FatesGlobals , only : fates_log + use EDBtranMod , only : check_layer_water + + implicit none + private + + public :: btran_sal_stress_fates + +contains + ! ===================================================================================== + + subroutine btran_sal_stress_fates( nsites, sites, bc_in) + + + ! --------------------------------------------------------------------------------- + ! Calculate the transpiration wetness function (BTRAN) and the root uptake + ! distribution (ROOTR). + ! Boundary conditions in: bc_in(s)%salinity_sl(j) salinity concontration[ppm] + ! Boundary conditions in: bc_in(s)%sulphide_sl(j) sulphide concontration[ppm] + ! Output cpatch%bstress_sal_ft(ft) + ! Output cpatch%bstress_sul_ft(ft) + ! --------------------------------------------------------------------------------- + + ! Arguments + + integer,intent(in) :: nsites + type(ed_site_type),intent(inout),target :: sites(nsites) + type(bc_in_type),intent(in) :: bc_in(nsites) + + ! + ! !LOCAL VARIABLES: + type(ed_patch_type),pointer :: cpatch ! Current Patch Pointer + type(ed_cohort_type),pointer :: ccohort ! Current cohort pointer + integer :: s ! site + integer :: j ! soil layer + integer :: ft ! plant functional type index + real(r8) :: salinity_node ! salinity in the soil water [ppt] + real(r8) :: rresis ! salinity limitation to transpiration independent + !------------------------------------------------------------------------------ + + do s = 1,nsites + + cpatch => sites(s)%oldest_patch + do while (associated(cpatch)) + + ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK) + + do ft = 1,numpft + cpatch%bstress_sal_ft(ft) = 0.0_r8 + do j = 1,bc_in(s)%nlevsoil + + ! Calculations are only relevant where liquid water exists + ! see clm_fates%wrap_btran for calculation with CLM/ELM + + if ( check_layer_water(bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j)) ) then + + salinity_node = bc_in(s)%salinity_sl(j) + + rresis = min( 1.244_r8/(1+exp((0.186_r8-salinity_node)/(-0.132_r8))), 1._r8) + + cpatch%bstress_sal_ft(ft) = cpatch%bstress_sal_ft(ft)+ & + cpatch%rootfr_ft(ft,j)*rresis + + end if + + end do !j + + end do !PFT + + cpatch => cpatch%younger + + end do + + end do + + + end subroutine btran_sal_stress_fates + + ! ==================================================================================== + +end module FatesBstressMod diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index ab5487d4ac..d27f6ebb63 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -402,15 +402,16 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) ! a_sapwood = a_leaf_tot / EDPftvarcon_inst%allom_latosa_int(FT)*1.e-4_r8 ! m2 sapwood = m2 leaf * cm2 sapwood/m2 leaf *1.0e-4m2 - - + ! or ... + ! a_sapwood = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * cCohort%hite ) * 1.e-4_r8 + v_sapwood = a_sapwood * z_stem ccohort_hydr%v_ag(n_hypool_leaf+1:n_hypool_ag) = v_sapwood / n_hypool_stem ! TRANSPORTING ROOT DEPTH & VOLUME !in special case where n_hypool_troot = 1, the node depth of the single troot pool !is the depth at which 50% total root distribution is attained - dcumul_rf = 1._r8/dble(n_hypool_troot) + dcumul_rf = 1._r8/real(n_hypool_troot,r8) do k=1,n_hypool_troot cumul_rf = dcumul_rf*k @@ -1026,11 +1027,11 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s csite_hydr%kmax_upper_shell(j,k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**-1._r8 + 1._r8/kmax_soil_total)**(-1._r8) csite_hydr%kmax_bound_shell(j,k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**-1._r8 + 1._r8/kmax_soil_total)**(-1._r8) csite_hydr%kmax_lower_shell(j,k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**-1._r8 + 1._r8/kmax_soil_total)**(-1._r8) end if if(j == 1) then if(csite_hydr%r_node_shell(j,k) <= csite_hydr%rs1(j)) then @@ -1041,11 +1042,11 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_1D / & log(csite_hydr%r_node_shell_1D(k)/csite_hydr%rs1(j))*hksat_s csite_hydr%kmax_upper_shell_1D(k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**-1._r8 + 1._r8/kmax_soil_total)**(-1._r8) csite_hydr%kmax_bound_shell_1D(k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**-1._r8 + 1._r8/kmax_soil_total)**(-1._r8) csite_hydr%kmax_lower_shell_1D(k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**-1._r8 + 1._r8/kmax_soil_total)**(-1._r8) end if end if else @@ -2554,7 +2555,7 @@ subroutine Hydraulics_1DSolve(cc_p, ft, z_node, v_node, ths_node, thr_node, kmax iterh1 = 0 do while( iterh1 == 0 .or. ((abs(we_local) > thresh .or. supsub_flag /= 0) .and. iterh1 < maxiter) ) dt_fac = max(imult*iterh1,1) - dt_fac2 = DBLE(dt_fac) + dt_fac2 = real(dt_fac,r8) dt_new = dtime/dt_fac2 !! restore initial states for a fresh attempt using new sub-timesteps @@ -4100,7 +4101,7 @@ subroutine swcCampbell_satfrac_from_psi(psi, psisat, B, satfrac) ! !LOCAL VARIABLES: !------------------------------------------------------------------------------ - satfrac = (psi/psisat)**(-1/B) + satfrac = (psi/psisat)**(-1.0_r8/B) end subroutine swcCampbell_satfrac_from_psi @@ -4438,7 +4439,7 @@ subroutine shellGeom(l_aroot, rs1, area, dz, r_out_shell, r_node_shell, v_shell) r_out_shell(nshell) = (pi_const*l_aroot/(area*dz))**(-0.5_r8) ! eqn(8) S98 if(nshell > 1) then do k = 1,nshell-1 - r_out_shell(k) = rs1*(r_out_shell(nshell)/rs1)**((k+0._r8)/nshell) ! eqn(7) S98 + r_out_shell(k) = rs1*(r_out_shell(nshell)/rs1)**((real(k,r8))/real(nshell,r8)) ! eqn(7) S98 enddo end if diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 56ca97ea23..4b181ed7b3 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -20,27 +20,27 @@ module FATESPlantRespPhotosynthMod ! !USES: - use FatesGlobals, only : endrun => fates_endrun - use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun + use FatesGlobals, only : fates_log use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : itrue use FatesInterfaceMod, only : hlm_use_planthydro use FatesInterfaceMod, only : hlm_parteh_mode use FatesInterfaceMod, only : numpft - use EDTypesMod, only : maxpft - use EDTypesMod, only : nlevleaf - use EDTypesMod, only : nclmax - - use PRTGenericMod, only : prt_carbon_allom_hyp - use PRTGenericMod, only : prt_cnp_flex_allom_hyp - use PRTGenericMod, only : all_carbon_elements - use PRTGenericMod, only : nitrogen_element - use PRTGenericMod, only : leaf_organ - use PRTGenericMod, only : fnrt_organ - use PRTGenericMod, only : sapw_organ - use PRTGenericMod, only : store_organ - use PRTGenericMod, only : repro_organ - use PRTGenericMod, only : struct_organ + use EDTypesMod, only : maxpft + use EDTypesMod, only : nlevleaf + use EDTypesMod, only : nclmax + use EDTypesMod, only : do_fates_salinity + use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_cnp_flex_allom_hyp + use PRTGenericMod, only : all_carbon_elements + use PRTGenericMod, only : nitrogen_element + use PRTGenericMod, only : leaf_organ + use PRTGenericMod, only : fnrt_organ + use PRTGenericMod, only : sapw_organ + use PRTGenericMod, only : store_organ + use PRTGenericMod, only : repro_organ + use PRTGenericMod, only : struct_organ ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -423,7 +423,11 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end if - + + if(do_fates_salinity)then + btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft) + endif + ! Scale for leaf nitrogen profile nscaler = exp(-kn(ft) * cumulative_lai) @@ -934,7 +938,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in if ( parsun_lsl <= 0._r8 ) then ! night time - anet_av_out = 0._r8 + anet_av_out = -lmr psn_out = 0._r8 rstoma_out = min(rsmax0, 1._r8/bbb * cf) diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 8ad11e3fbc..44f8d8d2e5 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -245,7 +245,7 @@ subroutine charecteristics_of_fuel ( currentSite ) endif ! FIX(RF,032414): needs refactoring. ! average water content !is this the correct metric? - timeav_swc = sum(currentSite%water_memory(1:numWaterMem)) / dble(numWaterMem) + timeav_swc = sum(currentSite%water_memory(1:numWaterMem)) / real(numWaterMem,r8) ! Equation B2 in Thonicke et al. 2010 ! live grass moisture content depends on upper soil layer fuel_moisture(lg_sf) = max(0.0_r8, 10.0_r8/9._r8 * timeav_swc - 1.0_r8/9.0_r8) @@ -408,8 +408,8 @@ subroutine rate_of_spread ( currentSite ) use SFParamsMod, only : SF_val_miner_total, & SF_val_part_dens, & SF_val_miner_damp, & - SF_val_fuel_energy, & - SF_val_wind_max + SF_val_fuel_energy + use FatesInterfaceMod, only : hlm_current_day, hlm_current_month type(ed_site_type), intent(in), target :: currentSite @@ -426,7 +426,6 @@ subroutine rate_of_spread ( currentSite ) real(r8) beta_ratio ! ratio of beta/beta_op real(r8) a_beta ! dummy variable for product of a* beta_ratio for react_v_opt equation real(r8) a,b,c,e ! function of fuel sav - real(r8) wind_elev_fire !wind speed (m/min) at elevevation relevant for fire logical,parameter :: debug_windspeed = .false. !for debugging @@ -493,21 +492,9 @@ subroutine rate_of_spread ( currentSite ) ! Equation A5 in Thonicke et al. 2010 ! phi_wind (unitless) - ! convert wind_elev_fire from m/min to ft/min for Rothermel ROS eqn - ! wind max per Lasslop et al 2014 to linearly reduce ROS for high wind speeds - !OLD! phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e)) - if (currentPatch%effect_wspeed .le. SF_val_wind_max) then - wind_elev_fire = currentPatch%effect_wspeed - phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e)) - if (debug_windspeed) write(fates_log(),*) 'SF wind LESS max ', currentPatch%effect_wspeed - if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day - else - !max condition 225 ft/min (FIREMIP Rabin table A10 JSBACH-Spitfire) convert to 68.577 m/min - wind_elev_fire = max(0.0_r8,(68.577_r8-0.5_r8*currentPatch%effect_wspeed)) - phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e)) - if (debug_windspeed) write(fates_log(),*) 'SF wind GREATER max ', currentPatch%effect_wspeed - if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day - endif + ! convert current_wspeed (wind at elev relevant to fire) from m/min to ft/min for Rothermel ROS eqn + phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e)) + ! ---propagating flux---- ! Equation A2 in Thonicke et al.2010 and Eq. 42 Rothermal 1972 @@ -769,7 +756,8 @@ subroutine area_burnt ( currentSite ) ! THIS SHOULD HAVE THE COLUMN AND LU AREA WEIGHT ALSO, NO? gridarea = km2_to_m2 ! 1M m2 in a km2 - !NF = number of lighting strikes per day per km2 + + ! NF = number of lighting strikes per day per km2 currentPatch%NF = ED_val_nignitions * currentPatch%area/area /365 ! If there are 15 lightening strickes per year, per km2. (approx from NASA product) @@ -784,9 +772,10 @@ subroutine area_burnt ( currentSite ) size_of_fire = ((3.1416_r8/(4.0_r8*lb))*((df+db)**2.0_r8)) !AB = daily area burnt = size fires in m2 * num ignitions * prob ignition starts fire + ! m2 per km2 per day currentPatch%AB = size_of_fire * currentPatch%NF * currentSite%FDI - patch_area_in_m2 = gridarea*currentPatch%area/area + patch_area_in_m2 = gridarea *currentPatch%area/area currentPatch%frac_burnt = currentPatch%AB / patch_area_in_m2 if(write_SF == itrue)then diff --git a/fire/SFParamsMod.F90 b/fire/SFParamsMod.F90 index a8556577fa..f40fea7a39 100644 --- a/fire/SFParamsMod.F90 +++ b/fire/SFParamsMod.F90 @@ -23,7 +23,6 @@ module SFParamsMod real(r8),protected :: SF_val_miner_damp real(r8),protected :: SF_val_max_durat real(r8),protected :: SF_val_durat_slope - real(r8),protected :: SF_val_wind_max ! Maximum wind speed expected by fire model (m/min) real(r8),protected :: SF_val_alpha_FMC(NFSC) real(r8),protected :: SF_val_CWD_frac(NCWD) real(r8),protected :: SF_val_max_decomp(NFSC) @@ -56,7 +55,6 @@ module SFParamsMod character(len=param_string_length),parameter :: SF_name_low_moisture_Slope = "fates_low_moisture_Slope" character(len=param_string_length),parameter :: SF_name_mid_moisture_Coeff = "fates_mid_moisture_Coeff" character(len=param_string_length),parameter :: SF_name_mid_moisture_Slope = "fates_mid_moisture_Slope" - character(len=param_string_length),parameter :: SF_name_wind_max = "fates_fire_wind_max" public :: SpitFireRegisterParams public :: SpitFireReceiveParams @@ -90,7 +88,6 @@ subroutine SpitFireParamsInit() SF_val_miner_damp = nan SF_val_max_durat = nan SF_val_durat_slope = nan - SF_val_wind_max = nan SF_val_CWD_frac(:) = nan @@ -150,9 +147,6 @@ subroutine SpitFireRegisterScalars(fates_params) character(len=param_string_length), parameter :: dim_names_scalar(1) = (/dimension_name_scalar/) - call fates_params%RegisterParameter(name=SF_name_wind_max, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=SF_name_fdi_a, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -191,9 +185,6 @@ subroutine SpitFireReceiveScalars(fates_params) class(fates_parameters_type), intent(inout) :: fates_params - call fates_params%RetreiveParameter(name=SF_name_wind_max, & - data=SF_val_wind_max) - call fates_params%RetreiveParameter(name=SF_name_fdi_a, & data=SF_val_fdi_a) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 7c31006a9c..8f0090d332 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -19,6 +19,8 @@ module EDInitMod use EDTypesMod , only : nuMWaterMem use EDTypesMod , only : maxpft use EDTypesMod , only : AREA + use EDTypesMod , only : init_spread_near_bare_ground + use EDTypesMod , only : init_spread_inventory use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_use_inventory_init @@ -68,13 +70,15 @@ subroutine init_site_vars( site_in ) ! !LOCAL VARIABLES: !---------------------------------------------------------------------- ! - allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,2)) + allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,1:nclmax)) allocate(site_in%demotion_rate(1:nlevsclass)) allocate(site_in%promotion_rate(1:nlevsclass)) allocate(site_in%imort_rate(1:nlevsclass,1:numpft)) allocate(site_in%fmort_rate(1:nlevsclass,1:numpft,1:nclmax)) allocate(site_in%fmort_rate_cambial(1:nlevsclass,1:numpft)) allocate(site_in%fmort_rate_crown(1:nlevsclass,1:numpft)) + allocate(site_in%growthflux_fusion(1:nlevsclass,1:numpft)) + ! end subroutine init_site_vars @@ -133,6 +137,9 @@ subroutine zero_site( site_in ) site_in%fmort_rate_cambial(:,:) = 0._r8 site_in%fmort_rate_crown(:,:) = 0._r8 + ! fusoin-induced growth flux of individuals + site_in%growthflux_fusion(:,:) = 0._r8 + ! demotion/promotion info site_in%demotion_rate(:) = 0._r8 site_in%demotion_carbonflux = 0._r8 @@ -227,7 +234,6 @@ subroutine set_site_properties( nsites, sites) sites(s)%frac_burnt = 0.0_r8 sites(s)%old_stock = 0.0_r8 - sites(s)%spread = 1.0_r8 end do return @@ -285,6 +291,13 @@ subroutine init_patches( nsites, sites, bc_in) if ( hlm_use_inventory_init.eq.itrue ) then + ! Initialize the site-level crown area spread factor (0-1) + ! It is likely that closed canopy forest inventories + ! have smaller spread factors than bare ground (they are crowded) + do s = 1, nsites + sites(s)%spread = init_spread_inventory + enddo + call initialize_sites_by_inventory(nsites,sites,bc_in) do s = 1, nsites @@ -299,6 +312,11 @@ subroutine init_patches( nsites, sites, bc_in) !FIX(SPM,032414) clean this up...inits out of this loop do s = 1, nsites + ! Initialize the site-level crown area spread factor (0-1) + ! It is likely that closed canopy forest inventories + ! have smaller spread factors than bare ground (they are crowded) + sites(s)%spread = init_spread_near_bare_ground + allocate(newp) newp%patchno = 1 diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index f1c705c3e1..4e4a1b11b5 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -330,9 +330,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year endif else - currentCohort%npp_acc_hold = currentCohort%npp_acc * dble(hlm_days_per_year) - currentCohort%gpp_acc_hold = currentCohort%gpp_acc * dble(hlm_days_per_year) - currentCohort%resp_acc_hold = currentCohort%resp_acc * dble(hlm_days_per_year) + currentCohort%npp_acc_hold = currentCohort%npp_acc * real(hlm_days_per_year,r8) + currentCohort%gpp_acc_hold = currentCohort%gpp_acc * real(hlm_days_per_year,r8) + currentCohort%resp_acc_hold = currentCohort%resp_acc * real(hlm_days_per_year,r8) endif currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n @@ -538,6 +538,12 @@ subroutine ed_total_balance_check (currentSite, call_index ) real(r8) :: error ! How much carbon did we gain or lose (should be zero!) real(r8) :: error_frac ! Error as a fraction of total biomass real(r8) :: net_flux ! Difference between recorded fluxes in and out. KgC/site + real(r8) :: leaf_c + real(r8) :: fnrt_c + real(r8) :: sapw_c + real(r8) :: store_c + real(r8) :: struct_c + real(r8) :: repro_c ! nb. There is no time associated with these variables ! because this routine can be called between any two @@ -579,7 +585,7 @@ subroutine ed_total_balance_check (currentSite, call_index ) ! burned_litter * new_patch%area !kG/site/day ! ----------------------------------------------------------------------------------- - if ( error_frac > 10e-6 ) then + if ( error_frac > 10e-6_r8 ) then write(fates_log(),*) 'carbon balance error detected' write(fates_log(),*) 'error fraction relative to biomass stock:',error_frac write(fates_log(),*) 'call index: ',call_index @@ -593,13 +599,17 @@ subroutine ed_total_balance_check (currentSite, call_index ) write(fates_log(),*) 'seeds',seed_stock write(fates_log(),*) 'previous total',currentSite%old_stock - if(debug)then + write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon + + ! If this is the first day of simulation, carbon balance reports but does not end the run + if( int(hlm_current_year*10000 + hlm_current_month*100 + hlm_current_day).ne.hlm_reference_date ) then + change_in_stock = 0.0_r8 biomass_stock = 0.0_r8 litter_stock = 0.0_r8 seed_stock = sum(currentSite%seed_bank)*AREA - currentPatch => currentSite%oldest_patch + currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) write(fates_log(),*) '---------------------------------------' write(fates_log(),*) currentPatch%area , sum(currentPatch%cwd_ag), sum(currentPatch%cwd_bg) @@ -607,19 +617,22 @@ subroutine ed_total_balance_check (currentSite, call_index ) write(fates_log(),*)'---' currentCohort => currentPatch%tallest do while(associated(currentCohort)) - write(fates_log(),*) 'structure: ',currentCohort%prt%GetState(struct_organ,all_carbon_elements) - write(fates_log(),*) 'storage: ',currentCohort%prt%GetState(store_organ,all_carbon_elements) + write(fates_log(),*) 'pft: ',currentCohort%pft + write(fates_log(),*) 'dbh: ',currentCohort%dbh + leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) + repro_c = currentCohort%prt%GetState(repro_organ,all_carbon_elements) + sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) + + write(fates_log(),*) 'lc: ',leaf_c,' dc: ',struct_c,' stc: ',store_c + write(fates_log(),*) 'fc: ',fnrt_c,' rc: ',repro_c,' sac: ',sapw_c write(fates_log(),*) 'N plant: ',currentCohort%n - currentCohort => currentCohort%shorter; - enddo !end cohort loop + currentCohort => currentCohort%shorter + enddo !end cohort loop currentPatch => currentPatch%younger - enddo !end patch loop - end if - - write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon - - ! If this is the first day of simulation, carbon balance reports but does not end the run - if( int(hlm_current_year*10000 + hlm_current_month*100 + hlm_current_day).ne.hlm_reference_date ) then + enddo !end patch loop write(fates_log(),*) 'aborting on date:',hlm_current_year,hlm_current_month,hlm_current_day call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -657,9 +670,9 @@ subroutine bypass_dynamics(currentSite) currentCohort%isnew=.false. - currentCohort%npp_acc_hold = currentCohort%npp_acc * dble(hlm_days_per_year) - currentCohort%gpp_acc_hold = currentCohort%gpp_acc * dble(hlm_days_per_year) - currentCohort%resp_acc_hold = currentCohort%resp_acc * dble(hlm_days_per_year) + currentCohort%npp_acc_hold = currentCohort%npp_acc * real(hlm_days_per_year,r8) + currentCohort%gpp_acc_hold = currentCohort%gpp_acc * real(hlm_days_per_year,r8) + currentCohort%resp_acc_hold = currentCohort%resp_acc * real(hlm_days_per_year,r8) currentCohort%npp_acc = 0.0_r8 currentCohort%gpp_acc = 0.0_r8 diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index ca2d33096a..f94cdccad2 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -43,6 +43,7 @@ module EDParamsMod real(r8),protected :: ED_val_patch_fusion_tol real(r8),protected :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry + ! two special parameters whose size is defined in the parameter file real(r8),protected,allocatable :: ED_val_history_sizeclass_bin_edges(:) real(r8),protected,allocatable :: ED_val_history_ageclass_bin_edges(:) @@ -87,7 +88,11 @@ module EDParamsMod real(r8),protected :: hydr_psicap ! sapwood water potential at which capillary reserves exhausted (MPa) character(len=param_string_length),parameter :: hydr_name_psicap = "fates_hydr_psicap" - + !Soil BGC parameters, mostly used for testing FATES when not coupled to the dynamics bgc hlm + ! ---------------------------------------------------------------------------------------------- + real(r8),protected :: bgc_soil_salinity ! site-level soil salinity for FATES when not coupled to dynamic soil BGC of salinity + character(len=param_string_length),parameter :: bgc_name_soil_salinity= "fates_soil_salinity" + ! Logging Control Parameters (ONLY RELEVANT WHEN USE_FATES_LOGGING = TRUE) ! ---------------------------------------------------------------------------------------------- @@ -155,6 +160,8 @@ subroutine FatesParamsInit() hydr_kmax_rsurf = nan hydr_psi0 = nan hydr_psicap = nan + + bgc_soil_salinity = nan logging_dbhmin = nan logging_collateral_frac = nan @@ -261,6 +268,9 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=hydr_name_psicap, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) + call fates_params%RegisterParameter(name=bgc_name_soil_salinity, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names) + call fates_params%RegisterParameter(name=logging_name_dbhmin, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -378,6 +388,9 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=hydr_name_psicap, & data=hydr_psicap) + + call fates_params%RetreiveParameter(name=bgc_name_soil_salinity, & + data=bgc_soil_salinity) call fates_params%RetreiveParameter(name=logging_name_dbhmin, & data=logging_dbhmin) @@ -450,6 +463,7 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'hydr_kmax_rsurf = ',hydr_kmax_rsurf write(fates_log(),fmt0) 'hydr_psi0 = ',hydr_psi0 write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap + write(fates_log(),fmt0) 'bgc_soil_salinity = ', bgc_soil_salinity write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin write(fates_log(),fmt0) 'logging_collateral_frac = ',logging_collateral_frac write(fates_log(),fmt0) 'logging_coll_under_frac = ',logging_coll_under_frac diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 161f0b7558..b32532f86e 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -32,55 +32,67 @@ module EDPftvarcon !ED specific variables. type, public :: EDPftvarcon_type - real(r8), allocatable :: pft_used (:) ! Switch to turn on and off PFTs - - real(r8), allocatable :: freezetol (:) ! minimum temperature tolerance - real(r8), allocatable :: wood_density (:) ! wood density g cm^-3 ... - real(r8), allocatable :: hgt_min (:) ! sapling height m + + real(r8), allocatable :: pft_used(:) ! Switch to turn on and off PFTs + real(r8), allocatable :: freezetol(:) ! minimum temperature tolerance + real(r8), allocatable :: wood_density(:) ! wood density g cm^-3 ... + real(r8), allocatable :: hgt_min(:) ! sapling height m real(r8), allocatable :: dbh_repro_threshold(:) ! diameter at which mature plants shift allocation - real(r8), allocatable :: dleaf (:) ! leaf characteristic dimension length (m) - real(r8), allocatable :: z0mr (:) ! ratio of roughness length of vegetation to height (-) - real(r8), allocatable :: displar (:) ! ratio of displacement height to canopy top height (-) - 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 + real(r8), allocatable :: dleaf(:) ! leaf characteristic dimension length (m) + real(r8), allocatable :: z0mr(:) ! ratio of roughness length of vegetation to height (-) + real(r8), allocatable :: displar(:) ! ratio of displacement height to canopy top height (-) + 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 (:) ! 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 :: seed_alloc_mature (:) ! fraction of carbon balance allocated to clonal reproduction. - real(r8), allocatable :: seed_alloc (:) ! fraction of carbon balance allocated to seeds. - real(r8), allocatable :: c2b (:) ! Carbon to biomass multiplier [kg/kgC] - real(r8), allocatable :: woody(:) - real(r8), allocatable :: stress_decid(:) - real(r8), allocatable :: season_decid(:) - real(r8), allocatable :: evergreen(:) - real(r8), allocatable :: slamax(:) - real(r8), allocatable :: slatop(:) + 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 :: roota_par(:) - real(r8), allocatable :: rootb_par(:) - real(r8), allocatable :: lf_flab(:) - real(r8), allocatable :: lf_fcel(:) - real(r8), allocatable :: lf_flig(:) - real(r8), allocatable :: fr_flab(:) - real(r8), allocatable :: fr_fcel(:) - real(r8), allocatable :: fr_flig(:) - real(r8), allocatable :: xl(:) - real(r8), allocatable :: clumping_index(:) ! factor describing how much self-occlusion - ! of leaf scattering elements decreases light interception - real(r8), allocatable :: c3psn(:) ! index defining the photosynthetic pathway C4 = 0, C3 = 1 - real(r8), allocatable :: vcmax25top(:) - real(r8), allocatable :: smpso(:) - real(r8), allocatable :: smpsc(:) + real(r8), allocatable :: seed_alloc_mature(:) ! fraction of carbon balance allocated to + ! clonal reproduction. + real(r8), allocatable :: seed_alloc(:) ! fraction of carbon balance allocated to seeds. + real(r8), allocatable :: c2b(:) ! Carbon to biomass multiplier [kg/kgC] + real(r8), allocatable :: woody(:) ! Does the plant have wood? (1=yes, 0=no) + + ! The following three PFT classes + ! are mutually exclusive + real(r8), allocatable :: stress_decid(:) ! Is the plant stress deciduous? (1=yes, 0=no) + real(r8), allocatable :: season_decid(:) ! Is the plant seasonally deciduous (1=yes, 0=no) + real(r8), allocatable :: evergreen(:) ! Is the plant an evergreen (1=yes, 0=no) + + real(r8), allocatable :: slamax(:) ! Maximum specific leaf area of plant (at bottom) [m2/gC] + real(r8), allocatable :: slatop(:) ! Specific leaf area at canopy top [m2/gC] - real(r8), allocatable :: maintresp_reduction_curvature(:) ! curvature of MR reduction as f(carbon storage), - ! 1=linear, 0=very curved - real(r8), allocatable :: maintresp_reduction_intercept(:) ! intercept of MR reduction as f(carbon storage), - ! 0=no throttling, 1=max throttling + real(r8), allocatable :: roota_par(:) ! Normalized Root profile scaling parameter A + real(r8), allocatable :: rootb_par(:) ! Normalized root profile scaling parameter B + real(r8), allocatable :: lf_flab(:) ! Leaf litter labile fraction [-] + real(r8), allocatable :: lf_fcel(:) ! Leaf litter cellulose fraction [-] + real(r8), allocatable :: lf_flig(:) ! Leaf litter lignan fraction [-] + real(r8), allocatable :: fr_flab(:) ! Fine-root litter labile fraction [-] + real(r8), allocatable :: fr_fcel(:) ! Fine-root litter cellulose fraction [-] + real(r8), allocatable :: fr_flig(:) ! Fine-root litter lignatn fraction [-] + real(r8), allocatable :: xl(:) ! Leaf-stem orientation index + real(r8), allocatable :: clumping_index(:) ! factor describing how much self-occlusion + ! of leaf scattering elements + ! decreases light interception + real(r8), allocatable :: c3psn(:) ! index defining the photosynthetic + ! pathway C4 = 0, C3 = 1 + real(r8), allocatable :: vcmax25top(:) ! maximum carboxylation rate of Rub. at 25C, + ! canopy top [umol CO2/m^2/s] + real(r8), allocatable :: smpso(:) ! Soil water potential at full stomatal opening + ! (non-HYDRO mode only) [mm] + real(r8), allocatable :: smpsc(:) ! Soil water potential at full stomatal closure + ! (non-HYDRO mode only) [mm] + + + real(r8), allocatable :: maintresp_reduction_curvature(:) ! curvature of MR reduction as f(carbon storage), + ! 1=linear, 0=very curved + real(r8), allocatable :: maintresp_reduction_intercept(:) ! intercept of MR reduction as f(carbon storage), + ! 0=no throttling, 1=max throttling real(r8), allocatable :: bmort(:) real(r8), allocatable :: mort_scalar_coldstress(:) real(r8), allocatable :: mort_scalar_cstarvation(:) @@ -147,6 +159,8 @@ module EDPftvarcon real(r8), allocatable :: allom_agb2(:) ! Parameter 2 for agb allometry real(r8), allocatable :: allom_agb3(:) ! Parameter 3 for agb allometry real(r8), allocatable :: allom_agb4(:) ! Parameter 3 for agb allometry + + real(r8), allocatable :: allom_frbstor_repro(:) ! fraction of bstrore for reproduction after mortality ! Prescribed Physiology Mode Parameters real(r8), allocatable :: prescribed_npp_canopy(:) ! this is only for the special @@ -203,7 +217,7 @@ module EDPftvarcon ! PFT Dimension real(r8), allocatable :: hydr_p_taper(:) ! xylem taper exponent - real(r8), allocatable :: hydr_rs2(:) ! absorbing root radius (mm) + real(r8), allocatable :: hydr_rs2(:) ! absorbing root radius (m) real(r8), allocatable :: hydr_srl(:) ! specific root length (m g-1) real(r8), allocatable :: hydr_rfrac_stem(:) ! fraction of total tree resistance from troot to canopy real(r8), allocatable :: hydr_avuln_gs(:) ! shape parameter for stomatal control of water vapor exiting leaf @@ -608,6 +622,10 @@ subroutine Register_PFT(this, fates_params) name = 'fates_allom_agb4' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_allom_frbstor_repro' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) name = 'fates_hydr_p_taper' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & @@ -1039,6 +1057,10 @@ subroutine Receive_PFT(this, fates_params) name = 'fates_allom_agb4' call fates_params%RetreiveParameterAllocate(name=name, & data=this%allom_agb4) + + name = 'fates_allom_frbstor_repro' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%allom_frbstor_repro) name = 'fates_hydr_p_taper' call fates_params%RetreiveParameterAllocate(name=name, & @@ -1715,7 +1737,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'allom_agb2 = ',EDPftvarcon_inst%allom_agb2 write(fates_log(),fmt0) 'allom_agb3 = ',EDPftvarcon_inst%allom_agb3 write(fates_log(),fmt0) 'allom_agb4 = ',EDPftvarcon_inst%allom_agb4 - + write(fates_log(),fmt0) 'allom_frbstor_repro = ',EDPftvarcon_inst%allom_frbstor_repro write(fates_log(),fmt0) 'hydr_p_taper = ',EDPftvarcon_inst%hydr_p_taper write(fates_log(),fmt0) 'hydr_rs2 = ',EDPftvarcon_inst%hydr_rs2 write(fates_log(),fmt0) 'hydr_srl = ',EDPftvarcon_inst%hydr_srl @@ -1939,6 +1961,22 @@ subroutine FatesCheckParams(is_master, parteh_mode) end if + ! Check if fraction of storage to reproduction is between 0-1 + ! ---------------------------------------------------------------------------------- + + if ( ( EDPftvarcon_inst%allom_frbstor_repro(ipft) < 0.0_r8 ) .or. & + ( EDPftvarcon_inst%allom_frbstor_repro(ipft) > 1.0_r8 ) ) then + + write(fates_log(),*) 'fraction of storage to reproduction' + write(fates_log(),*) ' after plants die, must be between' + write(fates_log(),*) ' 0 and 1' + write(fates_log(),*) ' PFT#: ',ipft + write(fates_log(),*) ' allom_frbstor_repro: ',EDPftvarcon_inst%allom_frbstor_repro(ipft) + write(fates_log(),*) ' Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + end if + ! Check if photosynthetic pathway is neither C3/C4 ! ---------------------------------------------------------------------------------- diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 2b94391a4c..d9c5f5c9d9 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -1,6 +1,8 @@ module EDTypesMod use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : ifalse + use FatesConstantsMod, only : itrue use FatesGlobals, only : fates_log use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use FatesHydraulicsMemMod, only : ed_cohort_hydr_type @@ -18,7 +20,9 @@ module EDTypesMod 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 + integer, parameter :: ican_ustory = 2 ! Nominal index for diagnostics that refer + ! to understory layers (all layers that + ! are not the top canopy layer) integer, parameter :: nlevleaf = 40 ! number of leaf layers in canopy layer integer, parameter :: maxpft = 15 ! maximum number of PFTs allowed @@ -79,7 +83,29 @@ module EDTypesMod ! TO-DO: THESE SHOULD BE PARAMETERS IN THE FILE OR NAMELIST - ADDING THESE ! WAS OUTSIDE THE SCOPE OF THE VERY LARGE CHANGESET WHERE THESE WERE FIRST ! INTRODUCED (RGK 03-2017) - logical, parameter :: do_ed_phenology = .true. + + logical, parameter :: do_ed_phenology = .true. + + + ! Flag to turn on/off salinity effects on the effective "btran" + ! btran stress function. + + logical, parameter :: do_fates_salinity = .false. + + + + ! This is the community level amount of spread expected in nearly-bare-ground + ! and inventory starting modes. + ! These are used to initialize only. These values will scale between + ! the PFT defined maximum and minimum crown area scaing parameters. + ! + ! A value of 1 indicates that + ! plants should have crown areas at maximum spread for their size and PFT. + ! A value of 0 means that they have the least amount of spread for their + ! size and PFT. + + real(r8), parameter :: init_spread_near_bare_ground = 1.0_r8 + real(r8), parameter :: init_spread_inventory = 0.0_r8 ! MODEL PARAMETERS @@ -121,18 +147,18 @@ module EDTypesMod integer , parameter :: N_HITE_BINS = 60 ! no. of hite bins used to distribute LAI ! COHORT TERMINATION - real(r8), parameter :: min_npm2 = 1.0E-8_r8 ! minimum cohort number density per m2 before termination - real(r8), parameter :: min_patch_area = 0.001_r8 ! smallest allowable patch area before termination - real(r8), parameter :: min_nppatch = 1.0E-11_r8 ! minimum number of cohorts per patch (min_npm2*min_patch_area) - real(r8), parameter :: min_n_safemath = 1.0E-15_r8 ! in some cases, we want to immediately remove super small - ! number densities of cohorts to prevent FPEs + real(r8), parameter :: min_npm2 = 1.0E-7_r8 ! minimum cohort number density per m2 before termination + real(r8), parameter :: min_patch_area = 0.01_r8 ! smallest allowable patch area before termination + real(r8), parameter :: min_nppatch = min_npm2*min_patch_area ! minimum number of cohorts per patch (min_npm2*min_patch_area) + real(r8), parameter :: min_n_safemath = 1.0E-12_r8 ! in some cases, we want to immediately remove super small + ! number densities of cohorts to prevent FPEs character*4 yearchar ! special mode to cause PFTs to create seed mass of all currently-existing PFTs logical, parameter :: homogenize_seed_pfts = .false. - !************************************ + !************************************ !** COHORT type structure ** !************************************ type ed_cohort_type @@ -182,6 +208,7 @@ module EDTypesMod ! type classification. We also maintain this in the main cohort memory ! because we don't want to continually re-calculate the cohort's position when ! performing size diagnostics at high-frequency calls + integer :: size_class_lasttimestep ! size class of the cohort at the end of the previous timestep (used for calculating growth flux) ! CARBON FLUXES @@ -328,6 +355,13 @@ module EDTypesMod integer :: ncan(nclmax,maxpft) ! number of total leaf layers for each canopy layer and pft !RADIATION FLUXES + + logical :: solar_zenith_flag ! integer flag specifying daylight (based on zenith angle) + real(r8) :: solar_zenith_angle ! solar zenith angle (radians) + + real(r8) :: gnd_alb_dif(maxSWb) ! ground albedo for diffuse rad, both bands (fraction) + real(r8) :: gnd_alb_dir(maxSWb) ! ground albedo for direct rad, both bands (fraction) + 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,maxpft,nlevleaf) ! shade fraction of direct light absorbed by each canopy @@ -393,7 +427,9 @@ module EDTypesMod ! 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(maxpft) ! btran calculated seperately for each PFT:- + real(r8) :: btran_ft(maxpft) ! btran calculated seperately for each PFT:- + real(r8) :: bstress_sal_ft(maxpft) ! bstress from salinity calculated seperately for each PFT:- + ! DISTURBANCE real(r8) :: disturbance_rates(n_dist_types) ! disturbance rate from 1) mortality @@ -572,20 +608,25 @@ module EDTypesMod ! TERMINATION, RECRUITMENT, DEMOTION, and DISTURBANCE - real(r8), allocatable :: terminated_nindivs(:,:,:) ! number of individuals that were in cohorts which were terminated this timestep, on size x pft x canopy array. - real(r8) :: termination_carbonflux(2) ! carbon flux from live to dead pools associated with termination mortality, per canopy level - real(r8) :: recruitment_rate(1:maxpft) ! number of individuals that were recruited into new cohorts - real(r8), allocatable :: demotion_rate(:) ! rate of individuals demoted from canopy to understory per FATES timestep - real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day] - real(r8), allocatable :: promotion_rate(:) ! rate of individuals promoted from understory to canopy per FATES timestep - real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day] + real(r8), allocatable :: terminated_nindivs(:,:,:) ! number of individuals that were in cohorts which were terminated this timestep, on size x pft x canopy array. + real(r8) :: termination_carbonflux(nclmax) ! carbon flux from live to dead pools associated with termination mortality, per canopy level + real(r8) :: recruitment_rate(1:maxpft) ! number of individuals that were recruited into new cohorts + real(r8), allocatable :: demotion_rate(:) ! rate of individuals demoted from canopy to understory per FATES timestep + real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day] + real(r8), allocatable :: promotion_rate(:) ! rate of individuals promoted from understory to canopy per FATES timestep + real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day] real(r8), allocatable :: imort_rate(:,:) ! rate of individuals killed due to impact mortality per year. on size x pft array real(r8) :: imort_carbonflux ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day] - real(r8), allocatable :: fmort_rate(:,:,:) ! rate of individuals killed due to fire mortality per year. on size x pft array - real(r8) :: fmort_carbonflux(2) ! biomass of individuals killed due to fire mortality per year. [gC/m2/sec] + + real(r8), allocatable :: fmort_rate(:,:,:) ! rate of individuals killed due to fire mortality per year. on size x pft x can-layer array + ! (1:nlevsclass,1:numpft,1:nclmax) + real(r8) :: fmort_carbonflux(nclmax) ! biomass of individuals killed due to fire mortality per year. [gC/m2/sec] real(r8), allocatable :: fmort_rate_cambial(:,:) ! rate of individuals killed due to fire mortality from cambial damage per year. on size x pft array real(r8), allocatable :: fmort_rate_crown(:,:) ! rate of individuals killed due to fire mortality from crown damage per year. on size x pft array + real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin due to fusion in a given day. on size x pft array + + ! some diagnostic-only (i.e. not resolved by ODE solver) flux of carbon to CWD and litter pools from termination and canopy mortality real(r8) :: CWD_AG_diagnostic_input_carbonflux(1:ncwd) ! diagnostic flux to AG CWD [kg C / m2 / yr] real(r8) :: CWD_BG_diagnostic_input_carbonflux(1:ncwd) ! diagnostic flux to BG CWD [kg C / m2 / yr] @@ -707,6 +748,10 @@ subroutine dump_patch(cpatch) write(fates_log(),*) 'pa%total_canopy_area = ',cpatch%total_canopy_area write(fates_log(),*) 'pa%total_tree_area = ',cpatch%total_tree_area write(fates_log(),*) 'pa%zstar = ',cpatch%zstar + write(fates_log(),*) 'pa%solar_zenith_flag = ',cpatch%solar_zenith_flag + write(fates_log(),*) 'pa%solar_zenith_angle = ',cpatch%solar_zenith_angle + write(fates_log(),*) 'pa%gnd_alb_dif = ',cpatch%gnd_alb_dif(:) + write(fates_log(),*) 'pa%gnd_alb_dir = ',cpatch%gnd_alb_dir(:) write(fates_log(),*) 'pa%c_stomata = ',cpatch%c_stomata write(fates_log(),*) 'pa%c_lblayer = ',cpatch%c_lblayer write(fates_log(),*) 'pa%disturbance_rate = ',cpatch%disturbance_rate diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 9eba8603c0..1160eb8bfd 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -8,7 +8,8 @@ module FatesHistoryInterfaceMod use FatesConstantsMod , only : calloc_abs_error use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun - + use EDTypesMod , only : nclmax + use EDTypesMod , only : ican_upper use FatesIODimensionsMod , only : fates_io_dimension_type use FatesIOVariableKindMod , only : fates_io_variable_kind_type use FatesHistoryVariableType , only : fates_history_variable_type @@ -181,6 +182,8 @@ module FatesHistoryInterfaceMod integer, private :: ih_ar_understory_si_scpf integer, private :: ih_ddbh_si_scpf + integer, private :: ih_growthflux_si_scpf + integer, private :: ih_growthflux_fusion_si_scpf integer, private :: ih_ba_si_scpf integer, private :: ih_m1_si_scpf integer, private :: ih_m2_si_scpf @@ -193,6 +196,7 @@ module FatesHistoryInterfaceMod integer, private :: ih_crownfiremort_si_scpf integer, private :: ih_cambialfiremort_si_scpf + integer, private :: ih_ar_si_scpf integer, private :: ih_ar_grow_si_scpf integer, private :: ih_ar_maint_si_scpf @@ -201,6 +205,7 @@ module FatesHistoryInterfaceMod integer, private :: ih_ar_crootm_si_scpf integer, private :: ih_ar_frootm_si_scpf + ! indices to (site x scls) variables integer, private :: ih_ba_si_scls integer, private :: ih_nplant_si_scls @@ -1270,6 +1275,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) use FatesSizeAgeTypeIndicesMod, only : get_agepft_class_index use FatesSizeAgeTypeIndicesMod, only : get_age_class_index use FatesSizeAgeTypeIndicesMod, only : get_height_index + use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index use EDTypesMod , only : nlevleaf use EDParamsMod, only : ED_val_history_height_bin_edges @@ -1408,6 +1414,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_ar_canopy_si_scpf => this%hvars(ih_ar_canopy_si_scpf)%r82d, & hio_ar_understory_si_scpf => this%hvars(ih_ar_understory_si_scpf)%r82d, & hio_ddbh_si_scpf => this%hvars(ih_ddbh_si_scpf)%r82d, & + hio_growthflux_si_scpf => this%hvars(ih_growthflux_si_scpf)%r82d, & + hio_growthflux_fusion_si_scpf => this%hvars(ih_growthflux_fusion_si_scpf)%r82d, & hio_ba_si_scpf => this%hvars(ih_ba_si_scpf)%r82d, & hio_nplant_si_scpf => this%hvars(ih_nplant_si_scpf)%r82d, & @@ -1570,6 +1578,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) do while(associated(ccohort)) ft = ccohort%pft + + call sizetype_class_index(ccohort%dbh, ccohort%pft, ccohort%size_class, ccohort%size_by_pft_class) ! Increment the number of cohorts per site hio_ncohorts_si(io_si) = hio_ncohorts_si(io_si) + 1._r8 @@ -1780,6 +1790,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! growth increment hio_ddbh_si_scpf(io_si,scpf) = hio_ddbh_si_scpf(io_si,scpf) + & ccohort%ddbhdt*ccohort%n + end if hio_agb_si_scls(io_si,scls) = hio_agb_si_scls(io_si,scls) + & @@ -1838,8 +1849,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_nplant_canopy_si_scpf(io_si,scpf) = hio_nplant_canopy_si_scpf(io_si,scpf) + ccohort%n hio_nplant_canopy_si_scls(io_si,scls) = hio_nplant_canopy_si_scls(io_si,scls) + ccohort%n - hio_lai_canopy_si_scls(io_si,scls) = hio_lai_canopy_si_scls(io_si,scls) + ccohort%n * ccohort%treelai - hio_sai_canopy_si_scls(io_si,scls) = hio_sai_canopy_si_scls(io_si,scls) + ccohort%n * ccohort%treesai + hio_lai_canopy_si_scls(io_si,scls) = hio_lai_canopy_si_scls(io_si,scls) + & + ccohort%treelai*ccohort%c_area * AREA_INV + hio_sai_canopy_si_scls(io_si,scls) = hio_sai_canopy_si_scls(io_si,scls) + & + ccohort%treesai*ccohort%c_area * AREA_INV hio_trimming_canopy_si_scls(io_si,scls) = hio_trimming_canopy_si_scls(io_si,scls) + & ccohort%n * ccohort%canopy_trim hio_crown_area_canopy_si_scls(io_si,scls) = hio_crown_area_canopy_si_scls(io_si,scls) + & @@ -1922,8 +1935,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_nplant_understory_si_scpf(io_si,scpf) = hio_nplant_understory_si_scpf(io_si,scpf) + ccohort%n hio_nplant_understory_si_scls(io_si,scls) = hio_nplant_understory_si_scls(io_si,scls) + ccohort%n - hio_lai_understory_si_scls(io_si,scls) = hio_lai_understory_si_scls(io_si,scls) + ccohort%n * ccohort%treelai - hio_sai_understory_si_scls(io_si,scls) = hio_sai_understory_si_scls(io_si,scls) + ccohort%n * ccohort%treesai + hio_lai_understory_si_scls(io_si,scls) = hio_lai_understory_si_scls(io_si,scls) + & + ccohort%treelai*ccohort%c_area * AREA_INV + hio_sai_understory_si_scls(io_si,scls) = hio_sai_understory_si_scls(io_si,scls) + & + ccohort%treelai*ccohort%c_area * AREA_INV hio_trimming_understory_si_scls(io_si,scls) = hio_trimming_understory_si_scls(io_si,scls) + & ccohort%n * ccohort%canopy_trim hio_crown_area_understory_si_scls(io_si,scls) = hio_crown_area_understory_si_scls(io_si,scls) + & @@ -1987,8 +2002,26 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! ! ccohort%canopy_layer_yesterday = real(ccohort%canopy_layer, r8) - + ! + ! growth flux of individuals into a given bin + ! track the actual growth here, the virtual growth from fusion lower down + if ( (scls - ccohort%size_class_lasttimestep ) .gt. 0) then + do i_scls = ccohort%size_class_lasttimestep + 1, scls + i_scpf = (ccohort%pft-1)*nlevsclass+i_scls + hio_growthflux_si_scpf(io_si,i_scpf) = hio_growthflux_si_scpf(io_si,i_scpf) + & + ccohort%n * days_per_year + end do + end if + ccohort%size_class_lasttimestep = scls + ! end associate + else ! i.e. cohort%isnew + ! + ! if cohort is new, track its growth flux into the first size bin + i_scpf = (ccohort%pft-1)*nlevsclass+1 + hio_growthflux_si_scpf(io_si,i_scpf) = hio_growthflux_si_scpf(io_si,i_scpf) + ccohort%n * days_per_year + ccohort%size_class_lasttimestep = 1 + ! end if ! resolve some canopy area profiles, both total and of occupied leaves @@ -2090,21 +2123,25 @@ subroutine update_history_dyn(this,nc,nsites,sites) i_scpf = (i_pft-1)*nlevsclass + i_scls ! ! termination mortality. sum of canopy and understory indices - hio_m6_si_scpf(io_si,i_scpf) = (sites(s)%terminated_nindivs(i_scls,i_pft,1) + & - sites(s)%terminated_nindivs(i_scls,i_pft,2)) * days_per_year + hio_m6_si_scpf(io_si,i_scpf) = sum(sites(s)%terminated_nindivs(i_scls,i_pft,1:nclmax)) * days_per_year + hio_m6_si_scls(io_si,i_scls) = hio_m6_si_scls(io_si,i_scls) + & - (sites(s)%terminated_nindivs(i_scls,i_pft,1) + & - sites(s)%terminated_nindivs(i_scls,i_pft,2)) * days_per_year + sum(sites(s)%terminated_nindivs(i_scls,i_pft,1:nclmax)) * days_per_year + ! ! add termination mortality to canopy and understory mortality hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + & - sites(s)%terminated_nindivs(i_scls,i_pft,1) * days_per_year + sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper) * days_per_year + hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & - sites(s)%terminated_nindivs(i_scls,i_pft,2) * days_per_year + sum(sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper+1:nclmax)) * days_per_year + hio_mortality_canopy_si_scpf(io_si,i_scpf) = hio_mortality_canopy_si_scpf(io_si,i_scpf) + & - sites(s)%terminated_nindivs(i_scls,i_pft,1) * days_per_year + sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper) * days_per_year + hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + & - sites(s)%terminated_nindivs(i_scls,i_pft,2) * days_per_year + sum(sites(s)%terminated_nindivs(i_scls,i_pft,ican_upper+1:nclmax)) * days_per_year + ! ! imort on its own hio_m4_si_scpf(io_si,i_scpf) = sites(s)%imort_rate(i_scls, i_pft) @@ -2113,7 +2150,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! add imort to other mortality terms. consider imort as understory mortality even if it happens in ! cohorts that may have been promoted as part of the patch creation, and use the pre-calculated site-level ! values to avoid biasing the results by the dramatically-reduced number densities in cohorts that are subject to imort - hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf)+ & + hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + & sites(s)%imort_rate(i_scls, i_pft) hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & sites(s)%imort_rate(i_scls, i_pft) @@ -2122,35 +2159,48 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & sites(s)%imort_rate(i_scls, i_pft) ! + ! fire mortality from the site-level diagnostic rates - hio_m5_si_scpf(io_si,i_scpf) = sum(sites(s)%fmort_rate(i_scls, i_pft,:)) - hio_m5_si_scls(io_si,i_scls) = hio_m5_si_scls(io_si,i_scls) + sum(sites(s)%fmort_rate(i_scls, i_pft,:)) + hio_m5_si_scpf(io_si,i_scpf) = sum(sites(s)%fmort_rate(i_scls, i_pft,1:nclmax)) + hio_m5_si_scls(io_si,i_scls) = hio_m5_si_scls(io_si,i_scls) + sum(sites(s)%fmort_rate(i_scls, i_pft,1:nclmax)) ! hio_crownfiremort_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_crown(i_scls, i_pft) hio_cambialfiremort_si_scpf(io_si,i_scpf) = sites(s)%fmort_rate_cambial(i_scls, i_pft) ! ! fire components of overall canopy and understory mortality hio_mortality_canopy_si_scpf(io_si,i_scpf) = hio_mortality_canopy_si_scpf(io_si,i_scpf) + & - sites(s)%fmort_rate(i_scls, i_pft,1) + sites(s)%fmort_rate(i_scls, i_pft,ican_upper) hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + & - sites(s)%fmort_rate(i_scls, i_pft,1) + sites(s)%fmort_rate(i_scls, i_pft,ican_upper) + + ! the fire mortality rates for each layer are total dead, since the usable + ! output will then normalize by the counts, we are allowed to sum over layers hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + & - sites(s)%fmort_rate(i_scls, i_pft,2) + sum(sites(s)%fmort_rate(i_scls, i_pft,ican_upper+1:nclmax)) + hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & - sites(s)%fmort_rate(i_scls, i_pft,2) + sum(sites(s)%fmort_rate(i_scls, i_pft,ican_upper+1:nclmax)) + ! ! carbon flux associated with mortality of trees dying by fire hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - sites(s)%fmort_carbonflux(1) + sites(s)%fmort_carbonflux(ican_upper) + hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sites(s)%fmort_carbonflux(2) + sum(sites(s)%fmort_carbonflux(ican_upper+1:nclmax)) + ! ! for scag variables, also treat as happening in the newly-disurbed patch - iscag = i_scls ! since imort is by definition something that only happens in newly disturbed patches, treat as such + hio_mortality_canopy_si_scag(io_si,iscag) = hio_mortality_canopy_si_scag(io_si,iscag) + & - sites(s)%fmort_rate(i_scls, i_pft,1) + sites(s)%fmort_rate(i_scls, i_pft,ican_upper) hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & - sites(s)%fmort_rate(i_scls, i_pft,2) + sum(sites(s)%fmort_rate(i_scls, i_pft,ican_upper+1:nclmax)) + + ! while in this loop, pass the fusion-induced growth rate flux to history + hio_growthflux_fusion_si_scpf(io_si,i_scpf) = hio_growthflux_fusion_si_scpf(io_si,i_scpf) + & + sites(s)%growthflux_fusion(i_scls, i_pft) * days_per_year + end do end do ! @@ -2165,6 +2215,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) sites(s)%fmort_carbonflux(:) = 0._r8 sites(s)%fmort_rate_cambial(:,:) = 0._r8 sites(s)%fmort_rate_crown(:,:) = 0._r8 + sites(s)%growthflux_fusion(:,:) = 0._r8 ! pass the recruitment rate as a flux to the history, and then reset the recruitment buffer do i_pft = 1, numpft @@ -2204,8 +2255,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & sites(s)%termination_carbonflux(ican_upper) * g_per_kg * days_per_sec * ha_per_m2 + hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sites(s)%termination_carbonflux(ican_ustory) * g_per_kg * days_per_sec * ha_per_m2 + sum(sites(s)%termination_carbonflux(ican_upper+1:nclmax)) * g_per_kg * days_per_sec * ha_per_m2 + ! and zero the site-level termination carbon flux variable sites(s)%termination_carbonflux(:) = 0._r8 ! @@ -2532,8 +2585,8 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) ! summarize radiation profiles through the canopy do ipft=1,numpft - do ican=1,nclmax - do ileaf=1,nlevleaf + do ican=1,nclmax ! cpatch%ncl_p ? + do ileaf=1,nlevleaf ! cpatch%ncan(ican,ipft) ? ! calculate where we are on multiplexed dimensions cnlfpft_indx = ileaf + (ican-1) * nlevleaf + (ipft-1) * nlevleaf * nclmax cnlf_indx = ileaf + (ican-1) * nlevleaf @@ -3791,6 +3844,16 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ddbh_si_scpf ) + call this%set_history_var(vname='GROWTHFLUX_SCPF', units = 'n/yr/ha', & + long='flux of individuals into a given size class bin via growth and recruitment',use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_growthflux_si_scpf ) + + call this%set_history_var(vname='GROWTHFLUX_FUSION_SCPF', units = 'n/yr/ha', & + long='flux of individuals into a given size class bin via fusion',use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_growthflux_fusion_si_scpf ) + call this%set_history_var(vname='DDBH_CANOPY_SCPF', units = 'cm/yr/ha', & long='diameter growth increment by pft/size',use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -4020,13 +4083,13 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_canopy_si_scls ) - call this%set_history_var(vname='LAI_CANOPY_SCLS', units = 'indiv/ha', & - long='number of canopy plants by size class', use_default='active', & + call this%set_history_var(vname='LAI_CANOPY_SCLS', units = 'm2/m2', & + long='Leaf are index (LAI) by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_lai_canopy_si_scls ) - call this%set_history_var(vname='SAI_CANOPY_SCLS', units = 'indiv/ha', & - long='number of canopy plants by size class', use_default='inactive', & + call this%set_history_var(vname='SAI_CANOPY_SCLS', units = 'm2/m2', & + long='stem area index(SAI) by size class', use_default='inactive', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_sai_canopy_si_scls ) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index e0af48ddd3..fb5002d2ac 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -18,6 +18,9 @@ module FatesInterfaceMod use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf use EDTypesMod , only : maxpft + use EDTypesMod , only : do_fates_salinity + use EDTypesMod , only : ncwd + use EDTypesMod , only : numWaterMem use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesGlobals , only : fates_global_verbose @@ -27,11 +30,13 @@ module FatesInterfaceMod use EDPftvarcon , only : FatesCheckParams use EDPftvarcon , only : EDPftvarcon_inst use EDParamsMod , only : FatesReportParams + use EDParamsMod , only : bgc_soil_salinity use PRTGenericMod , only : prt_carbon_allom_hyp use PRTGenericMod , only : prt_cnp_flex_allom_hyp use PRTAllometricCarbonMod, only : InitPRTGlobalAllometricCarbon ! use PRTAllometricCNPMod, only : InitPRTGlobalAllometricCNP + ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -403,6 +408,9 @@ module FatesInterfaceMod ! Soil suction potential of layers in each site, negative, [mm] real(r8), allocatable :: smp_sl(:) + + !soil salinity of layers in each site [ppt] + real(r8), allocatable :: salinity_sl(:) ! Effective porosity = porosity - vol_ic, of layers in each site [-] real(r8), allocatable :: eff_porosity_sl(:) @@ -581,6 +589,7 @@ module FatesInterfaceMod contains procedure, public :: zero_bcs + procedure, public :: set_bcs end type fates_interface_type @@ -667,6 +676,11 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in) allocate(bc_in%watsat_sl(nlevsoil_in)) allocate(bc_in%tempk_sl(nlevsoil_in)) allocate(bc_in%h2o_liqvol_sl(nlevsoil_in)) + + !BGC + if(do_fates_salinity) then + allocate(bc_in%salinity_sl(nlevsoil_in)) + endif ! Photosynthesis allocate(bc_in%filter_photo_pa(maxPatchesPerSite)) @@ -807,6 +821,10 @@ subroutine zero_bcs(this,s) this%bc_in(s)%tot_litc = 0.0_r8 this%bc_in(s)%snow_depth_si = 0.0_r8 this%bc_in(s)%frac_sno_eff_si = 0.0_r8 + + if(do_fates_salinity)then + this%bc_in(s)%salinity_sl(:) = 0.0_r8 + endif if (hlm_use_planthydro.eq.itrue) then @@ -864,6 +882,34 @@ subroutine zero_bcs(this,s) return end subroutine zero_bcs + + subroutine set_bcs(this,s) + + ! -------------------------------------------------------------------------------- + ! + ! This subroutine is called directly from the HLM to set boundary condition not yet + ! functional from hlm. This allows flexibility for model testing. + ! + ! This subroutine MUST BE CALLED AFTER the FATES PFT 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. + ! + ! -------------------------------------------------------------------------------- + implicit none + class(fates_interface_type), intent(inout) :: this + integer, intent(in) :: s + + ! Input boundaries + ! Warning: these "z" type variables + ! are written only once at the beginning + ! so THIS ROUTINE SHOULD NOT BE CALLED AFTER + ! INITIALIZATION + if(do_fates_salinity)then + this%bc_in(s)%salinity_sl(:) = bgc_soil_salinity + endif + + end subroutine set_bcs ! =================================================================================== @@ -924,8 +970,14 @@ subroutine set_fates_global_elements(use_fates) ! These values are used to define the restart file allocations and general structure ! of memory for the cohort arrays - fates_maxElementsPerPatch = max(maxCohortsPerPatch, & - numpft * nclmax * nlevleaf) + fates_maxElementsPerPatch = max(maxCohortsPerPatch, numpft, ncwd ) + + if (maxPatchesPerSite * fates_maxElementsPerPatch < numWaterMem) then + write(fates_log(), *) 'By using such a tiny number of maximum patches and maximum cohorts' + write(fates_log(), *) ' this could create problems for indexing in restart files' + write(fates_log(), *) ' The multiple of the two has to be greater than numWaterMem' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if fates_maxElementsPerSite = maxPatchesPerSite * fates_maxElementsPerPatch diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 480bea7c1e..5c682a9f02 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -35,6 +35,7 @@ module FatesInventoryInitMod use EDTypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst + implicit none private @@ -931,10 +932,9 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & endif ! Since spread is a canopy level calculation, we need to provide an initial guess here. - site_spread = 0.5_r8 call create_cohort(csite, cpatch, c_pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & - b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & - temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, 1, site_spread, bc_in) + b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & + temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, 1, csite%spread, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 8be110fc68..24eaa07311 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -5,12 +5,15 @@ module FatesRestartInterfaceMod use FatesConstantsMod , only : fates_avg_flag_length use FatesConstantsMod , only : fates_short_string_length use FatesConstantsMod , only : fates_long_string_length + use FatesConstantsMod , only : itrue + use FatesConstantsMod , only : ifalse use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use FatesIODimensionsMod, only : fates_io_dimension_type use FatesIOVariableKindMod, only : fates_io_variable_kind_type use FatesRestartVariableMod, only : fates_restart_variable_type use FatesInterfaceMod, only : bc_in_type + use FatesInterfaceMod, only : bc_out_type use FatesSizeAgeTypeIndicesMod, only : get_sizeage_class_index use PRTGenericMod, only : prt_global @@ -79,6 +82,7 @@ module FatesRestartInterfaceMod integer, private :: ir_canopy_layer_co integer, private :: ir_canopy_layer_yesterday_co integer, private :: ir_canopy_trim_co + integer, private :: ir_size_class_lasttimestep_co integer, private :: ir_dbh_co integer, private :: ir_height_co integer, private :: ir_laimemory_co @@ -100,6 +104,9 @@ module FatesRestartInterfaceMod integer, private :: ir_lmort_collateral_co integer, private :: ir_lmort_infra_co + ! Radiation + integer, private :: ir_solar_zenith_flag_pa + integer, private :: ir_solar_zenith_angle_pa integer, private :: ir_ddbhdt_co integer, private :: ir_resp_tstep_co @@ -108,6 +115,10 @@ module FatesRestartInterfaceMod integer, private :: ir_isnew_co integer, private :: ir_cwd_ag_pacw integer, private :: ir_cwd_bg_pacw + + integer, private :: ir_gnd_alb_dif_pasb + integer, private :: ir_gnd_alb_dir_pasb + integer, private :: ir_leaf_litter_paft integer, private :: ir_root_litter_paft integer, private :: ir_leaf_litter_in_paft @@ -117,11 +128,7 @@ module FatesRestartInterfaceMod integer, private :: ir_livegrass_pa integer, private :: ir_age_pa integer, private :: ir_area_pa - integer, private :: ir_fsun_paclftls - integer, private :: ir_fabd_sun_paclftls - integer, private :: ir_fabi_sun_paclftls - integer, private :: ir_fabd_sha_paclftls - integer, private :: ir_fabi_sha_paclftls + integer, private :: ir_watermem_siwm integer, private :: ir_prt_base ! Base index for all PRT variables @@ -190,6 +197,7 @@ module FatesRestartInterfaceMod procedure, public :: set_restart_vectors procedure, public :: create_patchcohort_structure procedure, public :: get_restart_vectors + procedure, public :: update_3dpatch_radiation ! private work functions procedure, private :: init_dim_kinds_maps procedure, private :: set_dim_indices @@ -607,6 +615,16 @@ subroutine define_restart_vars(this, initialize_variables) long_name='the number of cohorts per patch', units='unitless', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ncohort_pa ) + call this%set_restart_var(vname='fates_solar_zenith_flag_pa', vtype=cohort_int, & + long_name='switch specifying if zenith is positive', units='unitless', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_solar_zenith_flag_pa ) + + call this%set_restart_var(vname='fates_solar_zenith_angle_pa', vtype=cohort_r8, & + long_name='the angle of the solar zenith for each patch', units='radians', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_solar_zenith_angle_pa ) + + + ! 1D cohort Variables ! ----------------------------------------------------------------------------------- @@ -622,6 +640,10 @@ subroutine define_restart_vars(this, initialize_variables) long_name='ed cohort - canopy_trim', units='fraction', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_canopy_trim_co ) + call this%set_restart_var(vname='fates_size_class_lasttimestep', vtype=cohort_int, & + long_name='ed cohort - size-class last timestep', units='index', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_size_class_lasttimestep_co ) + call this%set_restart_var(vname='fates_dbh', vtype=cohort_r8, & long_name='ed cohort - diameter at breast height', units='cm', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dbh_co ) @@ -742,6 +764,16 @@ subroutine define_restart_vars(this, initialize_variables) units='kgC/m2', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cwd_bg_pacw ) + call this%set_restart_var(vname='fates_gnd_alb_dif', vtype=cohort_r8, & + long_name='ground albedo of diffuse radiation vis and ir', & + units='fraction', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_gnd_alb_dif_pasb ) + + call this%set_restart_var(vname='fates_gnd_alb_dir', vtype=cohort_r8, & + long_name='ground albedo of direct radiation vis and ir', & + units='fraction', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_gnd_alb_dir_pasb ) + call this%set_restart_var(vname='fates_leaf_litter', vtype=cohort_r8, & long_name='leaf litter, by patch x pft (non-respiring)', & units='kgC/m2', flushval = flushzero, & @@ -785,31 +817,6 @@ subroutine define_restart_vars(this, initialize_variables) long_name='are of the ED patch', units='m2', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_area_pa ) - ! These dimensions are pa "patch" cl "canopy layer" ft "functional type" ls "layer sublevel" - call this%set_restart_var(vname='fates_f_sun', vtype=cohort_r8, & - long_name='fraction of sunlit leaves, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fsun_paclftls ) - - call this%set_restart_var(vname='fates_fabd_sun_z', vtype=cohort_r8, & - long_name='sun fraction of direct light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabd_sun_paclftls ) - - call this%set_restart_var(vname='fates_fabi_sun_z', vtype=cohort_r8, & - long_name='sun fraction of indirect light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabi_sun_paclftls ) - - call this%set_restart_var(vname='fates_fabd_sha_z', vtype=cohort_r8, & - long_name='shade fraction of direct light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabd_sha_paclftls ) - - call this%set_restart_var(vname='fates_fabi_sha_z', vtype=cohort_r8, & - long_name='shade fraction of indirect light absorbed, by patch x can-layer x pft x sublayer', & - units='fraction', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fabi_sha_paclftls ) ! ! site x time level vars @@ -1020,14 +1027,13 @@ end subroutine set_restart_var subroutine set_restart_vectors(this,nc,nsites,sites) - use EDTypesMod, only : nclmax - use EDTypesMod, only : nlevleaf use FatesInterfaceMod, only : fates_maxElementsPerPatch use FatesInterfaceMod, only : numpft use EDTypesMod, only : ed_site_type use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd + use EDTypesMod, only : maxSWb use EDTypesMod, only : numWaterMem ! Arguments @@ -1052,7 +1058,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation + integer :: io_idx_pa_ib ! each SW band (vis/ir) per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site ! Some counters (for checking mostly) @@ -1096,9 +1102,12 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_seedrainflux_si => this%rvars(ir_seedrainflux_si)%r81d, & rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & + rio_solar_zenith_flag_pa => this%rvars(ir_solar_zenith_flag_pa)%int1d, & + rio_solar_zenith_angle_pa => this%rvars(ir_solar_zenith_angle_pa)%r81d, & rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%r81d, & rio_canopy_layer_yesterday_co => this%rvars(ir_canopy_layer_yesterday_co)%r81d, & rio_canopy_trim_co => this%rvars(ir_canopy_trim_co)%r81d, & + rio_size_class_lasttimestep => this%rvars(ir_size_class_lasttimestep_co)%int1d, & rio_dbh_co => this%rvars(ir_dbh_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & @@ -1109,7 +1118,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_gpp_acc_hold_co => this%rvars(ir_gpp_acc_hold_co)%r81d, & rio_resp_acc_hold_co => this%rvars(ir_resp_acc_hold_co)%r81d, & rio_npp_acc_hold_co => this%rvars(ir_npp_acc_hold_co)%r81d, & - rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & @@ -1124,6 +1132,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & rio_cwd_ag_pacw => this%rvars(ir_cwd_ag_pacw)%r81d, & rio_cwd_bg_pacw => this%rvars(ir_cwd_bg_pacw)%r81d, & + rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & + rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & rio_leaf_litter_paft => this%rvars(ir_leaf_litter_paft)%r81d, & rio_root_litter_paft => this%rvars(ir_root_litter_paft)%r81d, & rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & @@ -1133,11 +1143,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_fsun_paclftls => this%rvars(ir_fsun_paclftls)%r81d, & - rio_fabd_sun_z_paclftls => this%rvars(ir_fabd_sun_paclftls)%r81d, & - rio_fabi_sun_z_paclftls => this%rvars(ir_fabi_sun_paclftls)%r81d, & - rio_fabd_sha_z_paclftls => this%rvars(ir_fabd_sha_paclftls)%r81d, & - rio_fabi_sha_z_paclftls => this%rvars(ir_fabi_sha_paclftls)%r81d, & rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d ) totalCohorts = 0 @@ -1161,8 +1166,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st + io_idx_pa_ib = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st ! write seed_bank info(site-level, but PFT-resolved) do i = 1,numpft @@ -1232,6 +1237,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_canopy_layer_co(io_idx_co) = ccohort%canopy_layer rio_canopy_layer_yesterday_co(io_idx_co) = ccohort%canopy_layer_yesterday rio_canopy_trim_co(io_idx_co) = ccohort%canopy_trim + rio_size_class_lasttimestep(io_idx_co) = ccohort%size_class_lasttimestep rio_dbh_co(io_idx_co) = ccohort%dbh rio_height_co(io_idx_co) = ccohort%hite rio_laimemory_co(io_idx_co) = ccohort%laimemory @@ -1285,6 +1291,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! set cohorts per patch for IO rio_ncohort_pa( io_idx_co_1st ) = cohortsperpatch + ! Set zenith angle info + if ( cpatch%solar_zenith_flag ) then + rio_solar_zenith_flag_pa(io_idx_co_1st) = itrue + else + rio_solar_zenith_flag_pa(io_idx_co_1st) = ifalse + endif + rio_solar_zenith_angle_pa( io_idx_co_1st) = cpatch%solar_zenith_angle + if ( debug ) then write(fates_log(),*) 'offsetNumCohorts III ' & ,io_idx_co,cohortsperpatch @@ -1308,25 +1322,11 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 end do - if ( debug ) write(fates_log(),*) 'CLTV io_idx_pa_sunz 1 ',io_idx_pa_sunz - - if ( debug ) write(fates_log(),*) 'CLTV 1186 ',nlevleaf,numpft,nclmax - - 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) - rio_fabi_sun_z_paclftls(io_idx_pa_sunz) = cpatch%fabi_sun_z(i,j,k) - rio_fabd_sha_z_paclftls(io_idx_pa_sunz) = cpatch%fabd_sha_z(i,j,k) - rio_fabi_sha_z_paclftls(io_idx_pa_sunz) = cpatch%fabi_sha_z(i,j,k) - io_idx_pa_sunz = io_idx_pa_sunz + 1 - end do - end do + do i = 1,maxSWb + rio_gnd_alb_dif_pasb(io_idx_pa_ib) = cpatch%gnd_alb_dif(i) + rio_gnd_alb_dir_pasb(io_idx_pa_ib) = cpatch%gnd_alb_dir(i) + io_idx_pa_ib = io_idx_pa_ib + 1 end do - - if ( debug ) write(fates_log(),*) 'CLTV io_idx_pa_sunz 2 ',io_idx_pa_sunz - ! Set the first cohort index to the start of the next patch, increment ! by the maximum number of cohorts per patch @@ -1335,8 +1335,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! 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_ib = io_idx_co_1st io_idx_co = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st if ( debug ) then write(fates_log(),*) 'CLTV io_idx_co_1st ', io_idx_co_1st @@ -1408,8 +1408,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd - use EDTypesMod, only : nlevleaf - use EDTypesMod, only : nclmax + use EDTypesMod, only : maxSWb use FatesInterfaceMod, only : fates_maxElementsPerPatch use EDTypesMod, only : maxpft use EDTypesMod, only : area @@ -1605,8 +1604,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd - use EDTypesMod, only : nlevleaf - use EDTypesMod, only : nclmax + use EDTypesMod, only : maxSWb use FatesInterfaceMod, only : numpft use FatesInterfaceMod, only : fates_maxElementsPerPatch use EDTypesMod, only : numWaterMem @@ -1641,7 +1639,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation + integer :: io_idx_pa_ib ! each SW radiation band per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site ! Some counters (for checking mostly) @@ -1678,9 +1676,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_seedrainflux_si => this%rvars(ir_seedrainflux_si)%r81d, & rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & + rio_solar_zenith_flag_pa => this%rvars(ir_solar_zenith_flag_pa)%int1d, & + rio_solar_zenith_angle_pa => this%rvars(ir_solar_zenith_angle_pa)%r81d, & rio_canopy_layer_co => this%rvars(ir_canopy_layer_co)%r81d, & rio_canopy_layer_yesterday_co => this%rvars(ir_canopy_layer_yesterday_co)%r81d, & rio_canopy_trim_co => this%rvars(ir_canopy_trim_co)%r81d, & + rio_size_class_lasttimestep => this%rvars(ir_size_class_lasttimestep_co)%int1d, & rio_dbh_co => this%rvars(ir_dbh_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & @@ -1691,16 +1692,13 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_gpp_acc_hold_co => this%rvars(ir_gpp_acc_hold_co)%r81d, & rio_resp_acc_hold_co => this%rvars(ir_resp_acc_hold_co)%r81d, & rio_npp_acc_hold_co => this%rvars(ir_npp_acc_hold_co)%r81d, & - rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & - rio_lmort_direct_co => this%rvars(ir_lmort_direct_co)%r81d, & rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & rio_lmort_infra_co => this%rvars(ir_lmort_infra_co)%r81d, & - rio_ddbhdt_co => this%rvars(ir_ddbhdt_co)%r81d, & rio_resp_tstep_co => this%rvars(ir_resp_tstep_co)%r81d, & rio_pft_co => this%rvars(ir_pft_co)%int1d, & @@ -1708,6 +1706,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & rio_cwd_ag_pacw => this%rvars(ir_cwd_ag_pacw)%r81d, & rio_cwd_bg_pacw => this%rvars(ir_cwd_bg_pacw)%r81d, & + rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & + rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & rio_leaf_litter_paft => this%rvars(ir_leaf_litter_paft)%r81d, & rio_root_litter_paft => this%rvars(ir_root_litter_paft)%r81d, & rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & @@ -1717,11 +1717,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_fsun_paclftls => this%rvars(ir_fsun_paclftls)%r81d, & - rio_fabd_sun_z_paclftls => this%rvars(ir_fabd_sun_paclftls)%r81d, & - rio_fabi_sun_z_paclftls => this%rvars(ir_fabi_sun_paclftls)%r81d, & - rio_fabd_sha_z_paclftls => this%rvars(ir_fabd_sha_paclftls)%r81d, & - rio_fabi_sha_z_paclftls => this%rvars(ir_fabi_sha_paclftls)%r81d, & rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d ) totalcohorts = 0 @@ -1734,7 +1729,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st + io_idx_pa_ib = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st ! read seed_bank info(site-level, but PFT-resolved) @@ -1797,6 +1792,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%canopy_layer = rio_canopy_layer_co(io_idx_co) ccohort%canopy_layer_yesterday = rio_canopy_layer_yesterday_co(io_idx_co) ccohort%canopy_trim = rio_canopy_trim_co(io_idx_co) + ccohort%size_class_lasttimestep = rio_size_class_lasttimestep(io_idx_co) ccohort%dbh = rio_dbh_co(io_idx_co) ccohort%hite = rio_height_co(io_idx_co) ccohort%laimemory = rio_laimemory_co(io_idx_co) @@ -1846,17 +1842,22 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! ! deal with patch level fields here ! - cpatch%livegrass = rio_livegrass_pa(io_idx_co_1st) - cpatch%age = rio_age_pa(io_idx_co_1st) - cpatch%area = rio_area_pa(io_idx_co_1st) - cpatch%age_class = get_age_class_index(cpatch%age) - + cpatch%livegrass = rio_livegrass_pa(io_idx_co_1st) + cpatch%age = rio_age_pa(io_idx_co_1st) + cpatch%area = rio_area_pa(io_idx_co_1st) + cpatch%age_class = get_age_class_index(cpatch%age) + + ! Set zenith angle info + cpatch%solar_zenith_flag = ( rio_solar_zenith_flag_pa(io_idx_co_1st) .eq. itrue ) + cpatch%solar_zenith_angle = rio_solar_zenith_angle_pa(io_idx_co_1st) + ! set cohorts per patch for IO if ( debug ) then write(fates_log(),*) 'CVTL III ' & ,io_idx_co,cohortsperpatch endif + ! ! deal with patch level fields of arrays here ! @@ -1877,23 +1878,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 enddo - 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 - 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) - cpatch%fabi_sun_z(i,j,k) = rio_fabi_sun_z_paclftls(io_idx_pa_sunz) - cpatch%fabd_sha_z(i,j,k) = rio_fabd_sha_z_paclftls(io_idx_pa_sunz) - cpatch%fabi_sha_z(i,j,k) = rio_fabi_sha_z_paclftls(io_idx_pa_sunz) - io_idx_pa_sunz = io_idx_pa_sunz + 1 - end do - end do + do i = 1,maxSWb + cpatch%gnd_alb_dif(i) = rio_gnd_alb_dif_pasb(io_idx_pa_ib) + cpatch%gnd_alb_dir(i) = rio_gnd_alb_dir_pasb(io_idx_pa_ib) + io_idx_pa_ib = io_idx_pa_ib + 1 end do - - if ( debug ) write(fates_log(),*) 'CVTL io_idx_pa_sunz 2 ',io_idx_pa_sunz - + ! Now increment the position of the first cohort to that of the next ! patch @@ -1902,8 +1892,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! and max the number of allowed cohorts per patch io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st + io_idx_pa_ib = io_idx_co_1st io_idx_co = io_idx_co_1st - io_idx_pa_sunz = io_idx_co_1st if ( debug ) then write(fates_log(),*) 'CVTL io_idx_co_1st ', io_idx_co_1st @@ -1960,4 +1950,111 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end associate end subroutine get_restart_vectors + + ! ==================================================================================== + + subroutine update_3dpatch_radiation(this, nc, nsites, sites, bc_out) + + ! ------------------------------------------------------------------------- + ! This subroutine populates output boundary conditions related to radiation + ! called upon restart reads. + ! ------------------------------------------------------------------------- + + use EDTypesMod, only : ed_site_type + use EDTypesMod, only : ed_patch_type + use EDSurfaceRadiationMod, only : PatchNormanRadiation + use FatesInterfaceMod, only : hlm_numSWb + + ! !ARGUMENTS: + class(fates_restart_interface_type) , intent(inout) :: this + integer , intent(in) :: nc + integer , intent(in) :: nsites + type(ed_site_type) , intent(inout), target :: sites(nsites) + type(bc_out_type) , intent(inout) :: bc_out(nsites) + + ! locals + ! ---------------------------------------------------------------------------------- + type(ed_patch_type),pointer :: currentPatch ! current patch + integer :: s ! site counter + integer :: ib ! radiation band counter + integer :: ifp ! patch counter + + do s = 1, nsites + + ifp = 0 + currentpatch => sites(s)%oldest_patch + do while (associated(currentpatch)) + ifp = ifp+1 + + currentPatch%f_sun (:,:,:) = 0._r8 + currentPatch%fabd_sun_z (:,:,:) = 0._r8 + currentPatch%fabd_sha_z (:,:,:) = 0._r8 + currentPatch%fabi_sun_z (:,:,:) = 0._r8 + currentPatch%fabi_sha_z (:,:,:) = 0._r8 + currentPatch%fabd (:) = 0._r8 + currentPatch%fabi (:) = 0._r8 + + ! zero diagnostic radiation profiles + currentPatch%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0._r8 + currentPatch%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0._r8 + currentPatch%nrmlzd_parprof_dir_z(:,:,:) = 0._r8 + currentPatch%nrmlzd_parprof_dif_z(:,:,:) = 0._r8 + + ! ----------------------------------------------------------- + ! When calling norman radiation from the short-timestep + ! we are passing in boundary conditions to set the following + ! variables: + ! currentPatch%solar_zenith_flag (is there daylight?) + ! currentPatch%solar_zenith_angle (what is the value?) + ! ----------------------------------------------------------- + + if(currentPatch%solar_zenith_flag)then + + bc_out(s)%albd_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%albi_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%fabi_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%fabd_parb(ifp,:) = 0._r8 ! output HLM + bc_out(s)%ftdd_parb(ifp,:) = 1._r8 ! output HLM + bc_out(s)%ftid_parb(ifp,:) = 1._r8 ! output HLM + bc_out(s)%ftii_parb(ifp,:) = 1._r8 ! output HLM + + if (maxval(currentPatch%nrad(1,:))==0)then + !there are no leaf layers in this patch. it is effectively bare ground. + ! no radiation is absorbed + bc_out(s)%fabd_parb(ifp,:) = 0.0_r8 + bc_out(s)%fabi_parb(ifp,:) = 0.0_r8 + do ib = 1,hlm_numSWb + + ! REQUIRES A FIX HERE albd vs albi + + bc_out(s)%albd_parb(ifp,ib) = currentPatch%gnd_alb_dir(ib) + bc_out(s)%albd_parb(ifp,ib) = currentPatch%gnd_alb_dif(ib) + bc_out(s)%ftdd_parb(ifp,ib)= 1.0_r8 + bc_out(s)%ftid_parb(ifp,ib)= 1.0_r8 + bc_out(s)%ftii_parb(ifp,ib)= 1.0_r8 + enddo + else + + call PatchNormanRadiation (currentPatch, & + bc_out(s)%albd_parb(ifp,:), & + bc_out(s)%albi_parb(ifp,:), & + bc_out(s)%fabd_parb(ifp,:), & + bc_out(s)%fabi_parb(ifp,:), & + bc_out(s)%ftdd_parb(ifp,:), & + bc_out(s)%ftid_parb(ifp,:), & + bc_out(s)%ftii_parb(ifp,:)) + + endif ! is there vegetation? + + end if ! if the vegetation and zenith filter is active + + + currentPatch => currentPatch%younger + end do ! Loop linked-list patches + enddo ! Loop Sites + + return + end subroutine update_3dpatch_radiation + + end module FatesRestartInterfaceMod diff --git a/parameter_files/fates_params_14pfts.cdl b/parameter_files/fates_params_14pfts.cdl index e0caf701d2..cf3d4de880 100644 --- a/parameter_files/fates_params_14pfts.cdl +++ b/parameter_files/fates_params_14pfts.cdl @@ -126,6 +126,9 @@ variables: float fates_alloc_storage_cushion(fates_pft) ; fates_alloc_storage_cushion:units = "fraction" ; fates_alloc_storage_cushion:long_name = "maximum size of storage C pool, relative to maximum size of leaf C pool" ; + float fates_allom_frbstor_repro(fates_pft) ; + fates_allom_frbstor_repro:units = "fraction" ; + fates_allom_frbstor_repro:long_name = "fraction of bstore goes to reproduction after plant dies" ; float fates_allom_agb1(fates_pft) ; fates_allom_agb1:units = "variable" ; fates_allom_agb1:long_name = "Parameter 1 for agb allometry" ; @@ -315,7 +318,7 @@ variables: fates_hydr_rfrac_stem:units = "fraction" ; fates_hydr_rfrac_stem:long_name = "fraction of total tree resistance from troot to canopy" ; float fates_hydr_rs2(fates_pft) ; - fates_hydr_rs2:units = "mm" ; + fates_hydr_rs2:units = "m" ; fates_hydr_rs2:long_name = "absorbing root radius" ; float fates_hydr_srl(fates_pft) ; fates_hydr_srl:units = "m g-1" ; @@ -566,9 +569,6 @@ variables: float fates_fdi_b ; fates_fdi_b:units = "NA" ; fates_fdi_b:long_name = "spitfire parameter (unknown) " ; - float fates_fire_wind_max ; - fates_fire_wind_max:units = "m/min" ; - fates_fire_wind_max:long_name = "maximum wind speed expected by the fire model" ; float fates_fuel_energy ; fates_fuel_energy:units = "kJ/kg" ; fates_fuel_energy:long_name = "pitfire parameter, heat content of fuel" ; @@ -584,6 +584,9 @@ variables: float fates_part_dens ; fates_part_dens:units = "kg/m2" ; fates_part_dens:long_name = "spitfire parameter, oven dry particle density, Table A1 Thonicke et al 2010" ; + float fates_soil_salinity(fates_scalar) ; + fates_soil_salinity:units = "ppt" ; + fates_soil_salinity:long_name = "soil salinity used for model when not coupled to dynamic soil salinity" ; // global attributes: :history = "This file was made from FatesPFTIndexSwapper.py \n", @@ -693,6 +696,8 @@ data: "reproduction ", "structure "; + fates_allom_frbstor_repro = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; + fates_alloc_storage_cushion = 2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2 ; @@ -1181,8 +1186,6 @@ data: fates_fdi_b = 243.12 ; - fates_fire_wind_max = 45.718 ; - fates_fuel_energy = 18000 ; fates_max_durat = 240 ; @@ -1192,4 +1195,6 @@ data: fates_miner_total = 0.055 ; fates_part_dens = 513 ; + + fates_soil_salinity = 0.4; } diff --git a/parameter_files/fates_params_coastal_veg.cdl b/parameter_files/fates_params_coastal_veg.cdl new file mode 100644 index 0000000000..a6b1ac59d3 --- /dev/null +++ b/parameter_files/fates_params_coastal_veg.cdl @@ -0,0 +1,958 @@ +netcdf fates_params.2trop { +dimensions: + fates_NCWD = 4 ; + fates_history_age_bins = 7 ; + fates_history_height_bins = 6 ; + fates_history_size_bins = 13 ; + fates_hydr_organs = 4 ; + fates_litterclass = 6 ; + fates_pft = 2 ; + fates_scalar = 1 ; + fates_string_length = 60 ; + fates_variants = 2 ; +variables: + float fates_history_sizeclass_bin_edges(fates_history_size_bins) ; + fates_history_sizeclass_bin_edges:units = "cm" ; + fates_history_sizeclass_bin_edges:long_name = "Lower edges for DBH size class bins used in size-resolved cohort history output" ; + float fates_history_ageclass_bin_edges(fates_history_age_bins) ; + fates_history_ageclass_bin_edges:units = "yr" ; + fates_history_ageclass_bin_edges:long_name = "Lower edges for age class bins used in age-resolved patch history output" ; + float fates_FBD(fates_litterclass) ; + fates_FBD:units = "NA" ; + fates_FBD:long_name = "spitfire parameter related to fuel bulk density, see SFMain.F90" ; + float fates_SAV(fates_litterclass) ; + fates_SAV:units = "NA" ; + fates_SAV:long_name = "spitfire parameter related to surface area to volume ratio, see SFMain.F90" ; + float fates_alpha_FMC(fates_litterclass) ; + fates_alpha_FMC:units = "NA" ; + fates_alpha_FMC:long_name = "spitfire parameter related to fuel moisture content, Equation 6 Thonicke et al 2010" ; + float fates_history_height_bin_edges(fates_history_height_bins) ; + fates_history_height_bin_edges:units = "m" ; + fates_history_height_bin_edges:long_name = "Lower edges for height bins used in height-resolved history output" ; + float fates_low_moisture_Coeff(fates_litterclass) ; + fates_low_moisture_Coeff:units = "NA" ; + fates_low_moisture_Coeff:long_name = "spitfire parameter, equation B1 Thonicke et al 2010" ; + float fates_low_moisture_Slope(fates_litterclass) ; + fates_low_moisture_Slope:units = "NA" ; + fates_low_moisture_Slope:long_name = "spitfire parameter, equation B1 Thonicke et al 2010" ; + float fates_max_decomp(fates_litterclass) ; + fates_max_decomp:units = "kgC/m2/yr ?" ; + fates_max_decomp:long_name = "maximum rate of litter & CWD transfer from non-decomposing class into decomposing class" ; + float fates_mid_moisture(fates_litterclass) ; + fates_mid_moisture:units = "NA" ; + fates_mid_moisture:long_name = "spitfire litter moisture threshold to be considered medium dry" ; + float fates_mid_moisture_Coeff(fates_litterclass) ; + fates_mid_moisture_Coeff:units = "NA" ; + fates_mid_moisture_Coeff:long_name = "spitfire parameter, equation B1 Thonicke et al 2010" ; + float fates_mid_moisture_Slope(fates_litterclass) ; + fates_mid_moisture_Slope:units = "NA" ; + fates_mid_moisture_Slope:long_name = "spitfire parameter, equation B1 Thonicke et al 2010" ; + float fates_min_moisture(fates_litterclass) ; + fates_min_moisture:units = "NA" ; + fates_min_moisture:long_name = "spitfire litter moisture threshold to be considered very dry" ; + float fates_hydr_avuln_node(fates_hydr_organs, fates_pft) ; + fates_hydr_avuln_node:units = "unitless" ; + fates_hydr_avuln_node:long_name = "xylem vulnerability curve shape parameter" ; + float fates_hydr_epsil_node(fates_hydr_organs, fates_pft) ; + fates_hydr_epsil_node:units = "MPa" ; + fates_hydr_epsil_node:long_name = "bulk elastic modulus" ; + float fates_hydr_fcap_node(fates_hydr_organs, fates_pft) ; + fates_hydr_fcap_node:units = "unitless" ; + fates_hydr_fcap_node:long_name = "fraction of (1-resid_node) that is capillary in source" ; + float fates_hydr_kmax_node(fates_hydr_organs, fates_pft) ; + fates_hydr_kmax_node:units = "kgMPa/m/s" ; + fates_hydr_kmax_node:long_name = "maximum xylem conductivity per unit conducting xylem area" ; + float fates_hydr_p50_node(fates_hydr_organs, fates_pft) ; + fates_hydr_p50_node:units = "MPa" ; + fates_hydr_p50_node:long_name = "xylem water potential at 50% loss of conductivity" ; + float fates_hydr_pinot_node(fates_hydr_organs, fates_pft) ; + fates_hydr_pinot_node:units = "MPa" ; + fates_hydr_pinot_node:long_name = "osmotic potential at full turgor" ; + float fates_hydr_pitlp_node(fates_hydr_organs, fates_pft) ; + fates_hydr_pitlp_node:units = "MPa" ; + fates_hydr_pitlp_node:long_name = "turgor loss point" ; + float fates_hydr_resid_node(fates_hydr_organs, fates_pft) ; + fates_hydr_resid_node:units = "fraction" ; + fates_hydr_resid_node:long_name = "residual fraction" ; + float fates_hydr_thetas_node(fates_hydr_organs, fates_pft) ; + fates_hydr_thetas_node:units = "cm3/cm3" ; + fates_hydr_thetas_node:long_name = "saturated water content" ; + float fates_CWD_frac(fates_NCWD) ; + fates_CWD_frac:units = "fraction" ; + fates_CWD_frac:long_name = "fraction of woody (bdead+bsw) biomass destined for CWD pool" ; + char fates_pftname(fates_pft, fates_string_length) ; + fates_pftname:units = "unitless - string" ; + fates_pftname:long_name = "Description of plant type" ; + float fates_rootprof_beta(fates_variants, fates_pft) ; + fates_rootprof_beta:units = "unitless" ; + fates_rootprof_beta:long_name = "Rooting beta parameter, for C and N vertical discretization (NOT USED BY DEFAULT)" ; + float fates_alloc_storage_cushion(fates_pft) ; + fates_alloc_storage_cushion:units = "fraction" ; + fates_alloc_storage_cushion:long_name = "maximum size of storage C pool, relative to maximum size of leaf C pool" ; + float fates_allom_agb1(fates_pft) ; + fates_allom_agb1:units = "variable" ; + fates_allom_agb1:long_name = "Parameter 1 for agb allometry" ; + float fates_allom_agb2(fates_pft) ; + fates_allom_agb2:units = "variable" ; + fates_allom_agb2:long_name = "Parameter 2 for agb allometry" ; + float fates_allom_agb3(fates_pft) ; + fates_allom_agb3:units = "variable" ; + fates_allom_agb3:long_name = "Parameter 3 for agb allometry" ; + float fates_allom_agb4(fates_pft) ; + fates_allom_agb4:units = "variable" ; + fates_allom_agb4:long_name = "Parameter 4 for agb allometry" ; + float fates_allom_agb_frac(fates_pft) ; + fates_allom_agb_frac:units = "fraction" ; + fates_allom_agb_frac:long_name = "Fraction of woody biomass that is above ground" ; + float fates_allom_amode(fates_pft) ; + fates_allom_amode:units = "index" ; + fates_allom_amode:long_name = "AGB allometry function index" ; + float fates_allom_blca_expnt_diff(fates_pft) ; + fates_allom_blca_expnt_diff:units = "unitless" ; + fates_allom_blca_expnt_diff:long_name = "difference between allometric DBH:bleaf and DBH:crown area exponents" ; + float fates_allom_cmode(fates_pft) ; + fates_allom_cmode:units = "index" ; + fates_allom_cmode:long_name = "coarse root biomass allometry function index" ; + float fates_allom_d2bl1(fates_pft) ; + fates_allom_d2bl1:units = "variable" ; + fates_allom_d2bl1:long_name = "Parameter 1 for d2bl allometry (intercept)" ; + float fates_allom_d2bl2(fates_pft) ; + fates_allom_d2bl2:units = "variable" ; + fates_allom_d2bl2:long_name = "Parameter 2 for d2bl allometry (slope)" ; + float fates_allom_d2bl3(fates_pft) ; + fates_allom_d2bl3:units = "unitless" ; + fates_allom_d2bl3:long_name = "Parameter 3 for d2bl allometry (optional)" ; + float fates_allom_d2ca_coefficient_max(fates_pft) ; + fates_allom_d2ca_coefficient_max:units = "m2 cm^(-1/beta)" ; + fates_allom_d2ca_coefficient_max:long_name = "max (savanna) dbh to area multiplier factor where: area = n*d2ca_coeff*dbh^beta" ; + float fates_allom_d2ca_coefficient_min(fates_pft) ; + fates_allom_d2ca_coefficient_min:units = "m2 cm^(-1/beta)" ; + fates_allom_d2ca_coefficient_min:long_name = "min (forest) dbh to area multiplier factor where: area = n*d2ca_coeff*dbh^beta" ; + float fates_allom_d2h1(fates_pft) ; + fates_allom_d2h1:units = "variable" ; + fates_allom_d2h1:long_name = "Parameter 1 for d2h allometry (intercept, or c)" ; + float fates_allom_d2h2(fates_pft) ; + fates_allom_d2h2:units = "variable" ; + fates_allom_d2h2:long_name = "Parameter 2 for d2h allometry (slope, or m)" ; + float fates_allom_d2h3(fates_pft) ; + fates_allom_d2h3:units = "variable" ; + fates_allom_d2h3:long_name = "Parameter 3 for d2h allometry (optional)" ; + float fates_allom_dbh_maxheight(fates_pft) ; + fates_allom_dbh_maxheight:units = "cm" ; + fates_allom_dbh_maxheight:long_name = "the diameter (if any) corresponding to maximum height, diameters may increase beyond this" ; + float fates_allom_fmode(fates_pft) ; + fates_allom_fmode:units = "index" ; + fates_allom_fmode:long_name = "fine root biomass allometry function index" ; + float fates_allom_hmode(fates_pft) ; + fates_allom_hmode:units = "index" ; + fates_allom_hmode:long_name = "height allometry function index" ; + float fates_allom_l2fr(fates_pft) ; + fates_allom_l2fr:units = "gC/gC" ; + fates_allom_l2fr:long_name = "Allocation parameter: fine root C per leaf C" ; + float fates_allom_latosa_int(fates_pft) ; + fates_allom_latosa_int:units = "ratio" ; + fates_allom_latosa_int:long_name = "Leaf area to sap area ratio, intercept [m2/cm2]" ; + float fates_allom_latosa_slp(fates_pft) ; + fates_allom_latosa_slp:units = "unitless" ; + fates_allom_latosa_slp:long_name = "Leaf area to sap area ratio, slope (optional)" ; + float fates_allom_lmode(fates_pft) ; + fates_allom_lmode:units = "index" ; + fates_allom_lmode:long_name = "leaf biomass allometry function index" ; + float fates_allom_sai_scaler(fates_pft) ; + fates_allom_sai_scaler:units = "m2/gC" ; + fates_allom_sai_scaler:long_name = "allometric ratio of SAI to target bleaf" ; + float fates_allom_smode(fates_pft) ; + fates_allom_smode:units = "index" ; + fates_allom_smode:long_name = "sapwood allometry function index" ; + float fates_allom_stmode(fates_pft) ; + fates_allom_stmode:units = "index" ; + fates_allom_stmode:long_name = "storage allometry function index" ; + float fates_allom_frbstor_repro(fates_pft) ; + fates_allom_frbstor_repro:units = "fraction" ; + fates_allom_frbstor_repro:long_name = "fraction of bstore goes to reproduction after plant dies" ; + float fates_branch_turnover(fates_pft) ; + fates_branch_turnover:units = "yr-1" ; + fates_branch_turnover:long_name = "turnover time of branches" ; + float fates_c2b(fates_pft) ; + fates_c2b:units = "ratio" ; + fates_c2b:long_name = "Carbon to biomass multiplier of bulk structural tissues" ; + float fates_displar(fates_pft) ; + fates_displar:units = "unitless" ; + fates_displar:long_name = "Ratio of displacement height to canopy top height" ; + float fates_fire_alpha_SH(fates_pft) ; + fates_fire_alpha_SH:units = "NA" ; + fates_fire_alpha_SH:long_name = "spitfire parameter, alpha scorch height, Equation 16 Thonicke et al 2010" ; + float fates_fire_bark_scaler(fates_pft) ; + fates_fire_bark_scaler:units = "fraction" ; + fates_fire_bark_scaler:long_name = "the thickness of a cohorts bark as a fraction of its dbh" ; + float fates_fire_crown_depth_frac(fates_pft) ; + fates_fire_crown_depth_frac:units = "fraction" ; + fates_fire_crown_depth_frac:long_name = "the depth of a cohorts crown as a fraction of its height" ; + float fates_fire_crown_kill(fates_pft) ; + fates_fire_crown_kill:units = "NA" ; + fates_fire_crown_kill:long_name = "fire parameter, see equation 22 in Thonicke et al 2010" ; + float fates_fr_fcel(fates_pft) ; + fates_fr_fcel:units = "fraction" ; + fates_fr_fcel:long_name = "Fine root litter cellulose fraction" ; + float fates_fr_flab(fates_pft) ; + fates_fr_flab:units = "fraction" ; + fates_fr_flab:long_name = "Fine root litter labile fraction" ; + float fates_fr_flig(fates_pft) ; + fates_fr_flig:units = "fraction" ; + fates_fr_flig:long_name = "Fine root litter lignin fraction" ; + float fates_froot_cn_ratio(fates_pft) ; + fates_froot_cn_ratio:units = "gC/gN" ; + fates_froot_cn_ratio:long_name = "Fine root C:N" ; + float fates_grperc(fates_pft) ; + fates_grperc:units = "unitless" ; + fates_grperc:long_name = "Growth respiration factor" ; + float fates_hydr_avuln_gs(fates_pft) ; + fates_hydr_avuln_gs:units = "unitless" ; + fates_hydr_avuln_gs:long_name = "shape parameter for stomatal control of water vapor exiting leaf" ; + float fates_hydr_p50_gs(fates_pft) ; + fates_hydr_p50_gs:units = "MPa" ; + fates_hydr_p50_gs:long_name = "water potential at 50% loss of stomatal conductance" ; + float fates_hydr_p_taper(fates_pft) ; + fates_hydr_p_taper:units = "unitless" ; + fates_hydr_p_taper:long_name = "xylem taper exponent" ; + float fates_hydr_rfrac_stem(fates_pft) ; + fates_hydr_rfrac_stem:units = "fraction" ; + fates_hydr_rfrac_stem:long_name = "fraction of total tree resistance from troot to canopy" ; + float fates_hydr_rs2(fates_pft) ; + fates_hydr_rs2:units = "mm" ; + fates_hydr_rs2:long_name = "absorbing root radius" ; + float fates_hydr_srl(fates_pft) ; + fates_hydr_srl:units = "m g-1" ; + fates_hydr_srl:long_name = "specific root length" ; + float fates_leaf_BB_slope(fates_pft) ; + fates_leaf_BB_slope:units = "unitless" ; + fates_leaf_BB_slope:long_name = "stomatal slope parameter, as per Ball-Berry" ; + float fates_leaf_c3psn(fates_pft) ; + fates_leaf_c3psn:units = "flag" ; + fates_leaf_c3psn:long_name = "Photosynthetic pathway (1=c3, 0=c4)" ; + float fates_leaf_clumping_index(fates_pft) ; + fates_leaf_clumping_index:units = "fraction (0-1)" ; + fates_leaf_clumping_index:long_name = "factor describing how much self-occlusion of leaf scattering elements decreases light interception" ; + float fates_leaf_cn_ratio(fates_pft) ; + fates_leaf_cn_ratio:units = "gC/gN" ; + fates_leaf_cn_ratio:long_name = "Leaf C:N" ; + float fates_leaf_diameter(fates_pft) ; + fates_leaf_diameter:units = "m" ; + fates_leaf_diameter:long_name = "Characteristic leaf dimension" ; + float fates_leaf_jmaxha(fates_pft) ; + fates_leaf_jmaxha:units = "J/mol" ; + fates_leaf_jmaxha:long_name = "activation energy for jmax" ; + float fates_leaf_jmaxhd(fates_pft) ; + fates_leaf_jmaxhd:units = "J/mol" ; + fates_leaf_jmaxhd:long_name = "deactivation energy for jmax" ; + float fates_leaf_jmaxse(fates_pft) ; + fates_leaf_jmaxse:units = "J/mol/K" ; + fates_leaf_jmaxse:long_name = "entropy term for jmax" ; + float fates_leaf_long(fates_pft) ; + fates_leaf_long:units = "yr" ; + fates_leaf_long:long_name = "Leaf longevity (ie turnover timescale)" ; + float fates_leaf_slatop(fates_pft) ; + fates_leaf_slatop:units = "m^2/gC" ; + fates_leaf_slatop:long_name = "Specific Leaf Area (SLA) at top of canopy, projected area basis" ; + float fates_leaf_stor_priority(fates_pft) ; + fates_leaf_stor_priority:units = "unitless" ; + fates_leaf_stor_priority:long_name = "factor governing priority of replacing storage with NPP" ; + float fates_leaf_tpuha(fates_pft) ; + fates_leaf_tpuha:units = "J/mol" ; + fates_leaf_tpuha:long_name = "activation energy for tpu" ; + float fates_leaf_tpuhd(fates_pft) ; + fates_leaf_tpuhd:units = "J/mol" ; + fates_leaf_tpuhd:long_name = "deactivation energy for tpu" ; + float fates_leaf_tpuse(fates_pft) ; + fates_leaf_tpuse:units = "J/mol/K" ; + fates_leaf_tpuse:long_name = "entropy term for tpu" ; + float fates_leaf_vcmax25top(fates_pft) ; + fates_leaf_vcmax25top:units = "umol CO2/m^2/s" ; + fates_leaf_vcmax25top:long_name = "maximum carboxylation rate of Rub. at 25C, canopy top" ; + float fates_leaf_vcmaxha(fates_pft) ; + fates_leaf_vcmaxha:units = "J/mol" ; + fates_leaf_vcmaxha:long_name = "activation energy for vcmax" ; + float fates_leaf_vcmaxhd(fates_pft) ; + fates_leaf_vcmaxhd:units = "J/mol" ; + fates_leaf_vcmaxhd:long_name = "deactivation energy for vcmax" ; + float fates_leaf_vcmaxse(fates_pft) ; + fates_leaf_vcmaxse:units = "J/mol/K" ; + fates_leaf_vcmaxse:long_name = "entropy term for vcmax" ; + float fates_leaf_xl(fates_pft) ; + fates_leaf_xl:units = "unitless" ; + fates_leaf_xl:long_name = "Leaf/stem orientation index" ; + float fates_lf_fcel(fates_pft) ; + fates_lf_fcel:units = "fraction" ; + fates_lf_fcel:long_name = "Leaf litter cellulose fraction" ; + float fates_lf_flab(fates_pft) ; + fates_lf_flab:units = "fraction" ; + fates_lf_flab:long_name = "Leaf litter labile fraction" ; + float fates_lf_flig(fates_pft) ; + fates_lf_flig:units = "fraction" ; + fates_lf_flig:long_name = "Leaf litter lignin fraction" ; + float fates_maintresp_reduction_curvature(fates_pft) ; + fates_maintresp_reduction_curvature:units = "unitless (0-1)" ; + fates_maintresp_reduction_curvature:long_name = "curvature of MR reduction as f(carbon storage), 1=linear, 0=very curved" ; + float fates_maintresp_reduction_intercept(fates_pft) ; + fates_maintresp_reduction_intercept:units = "unitless (0-1)" ; + fates_maintresp_reduction_intercept:long_name = "intercept of MR reduction as f(carbon storage), 0=no throttling, 1=max throttling" ; + float fates_mort_bmort(fates_pft) ; + fates_mort_bmort:units = "1/yr" ; + fates_mort_bmort:long_name = "background mortality rate" ; + float fates_mort_freezetol(fates_pft) ; + fates_mort_freezetol:units = "NA" ; + fates_mort_freezetol:long_name = "minimum temperature tolerance (NOT USED)" ; + float fates_mort_hf_sm_threshold(fates_pft) ; + fates_mort_hf_sm_threshold:units = "unitless" ; + fates_mort_hf_sm_threshold:long_name = "soil moisture (btran units) at which drought mortality begins for non-hydraulic model" ; + float fates_mort_scalar_coldstress(fates_pft) ; + fates_mort_scalar_coldstress:units = "1/yr" ; + fates_mort_scalar_coldstress:long_name = "maximum mortality rate from cold stress" ; + float fates_mort_scalar_cstarvation(fates_pft) ; + fates_mort_scalar_cstarvation:units = "1/yr" ; + fates_mort_scalar_cstarvation:long_name = "maximum mortality rate from carbon starvation" ; + float fates_mort_scalar_hydrfailure(fates_pft) ; + fates_mort_scalar_hydrfailure:units = "1/yr" ; + fates_mort_scalar_hydrfailure:long_name = "maximum mortality rate from hydraulic failure" ; + float fates_pft_used(fates_pft) ; + fates_pft_used:units = "0 = off (dont use), 1 = on (use)" ; + fates_pft_used:long_name = "Switch to turn on and off PFTs (also see fates_initd for cold-start)" ; + float fates_phen_evergreen(fates_pft) ; + fates_phen_evergreen:units = "logical flag" ; + fates_phen_evergreen:long_name = "Binary flag for evergreen leaf habit" ; + float fates_phen_season_decid(fates_pft) ; + fates_phen_season_decid:units = "logical flag" ; + fates_phen_season_decid:long_name = "Binary flag for seasonal-deciduous leaf habit" ; + float fates_phen_stress_decid(fates_pft) ; + fates_phen_stress_decid:units = "logical flag" ; + fates_phen_stress_decid:long_name = "Binary flag for stress-deciduous leaf habit" ; + float fates_prescribed_mortality_canopy(fates_pft) ; + fates_prescribed_mortality_canopy:units = "1/yr" ; + fates_prescribed_mortality_canopy:long_name = "mortality rate of canopy trees for prescribed physiology mode" ; + float fates_prescribed_mortality_understory(fates_pft) ; + fates_prescribed_mortality_understory:units = "1/yr" ; + fates_prescribed_mortality_understory:long_name = "mortality rate of understory trees for prescribed physiology mode" ; + float fates_prescribed_npp_canopy(fates_pft) ; + fates_prescribed_npp_canopy:units = "gC / m^2 / yr" ; + fates_prescribed_npp_canopy:long_name = "NPP per unit crown area of canopy trees for prescribed physiology mode" ; + float fates_prescribed_npp_understory(fates_pft) ; + fates_prescribed_npp_understory:units = "gC / m^2 / yr" ; + fates_prescribed_npp_understory:long_name = "NPP per unit crown area of understory trees for prescribed physiology mode" ; + float fates_prescribed_recruitment(fates_pft) ; + fates_prescribed_recruitment:units = "n/yr" ; + fates_prescribed_recruitment:long_name = "recruitment rate for prescribed physiology mode" ; + float fates_recruit_hgt_min(fates_pft) ; + fates_recruit_hgt_min:units = "m" ; + fates_recruit_hgt_min:long_name = "the minimum height (ie starting height) of a newly recruited plant" ; + float fates_recruit_initd(fates_pft) ; + fates_recruit_initd:units = "stems/m2" ; + fates_recruit_initd:long_name = "initial seedling density for a cold-start near-bare-ground simulation" ; + float fates_rholnir(fates_pft) ; + fates_rholnir:units = "fraction" ; + fates_rholnir:long_name = "Leaf reflectance: near-IR" ; + float fates_rholvis(fates_pft) ; + fates_rholvis:units = "fraction" ; + fates_rholvis:long_name = "Leaf reflectance: visible" ; + float fates_rhosnir(fates_pft) ; + fates_rhosnir:units = "fraction" ; + fates_rhosnir:long_name = "Stem reflectance: near-IR" ; + float fates_rhosvis(fates_pft) ; + fates_rhosvis:units = "fraction" ; + fates_rhosvis:long_name = "Stem reflectance: visible" ; + float fates_root_long(fates_pft) ; + fates_root_long:units = "yr" ; + fates_root_long:long_name = "root longevity (alternatively, turnover time)" ; + float fates_roota_par(fates_pft) ; + fates_roota_par:units = "1/m" ; + fates_roota_par:long_name = "CLM rooting distribution parameter" ; + float fates_rootb_par(fates_pft) ; + fates_rootb_par:units = "1/m" ; + fates_rootb_par:long_name = "CLM rooting distribution parameter" ; + float fates_seed_alloc(fates_pft) ; + fates_seed_alloc:units = "fraction" ; + fates_seed_alloc:long_name = "fraction of available carbon balance allocated to seeds" ; + float fates_seed_alloc_mature(fates_pft) ; + fates_seed_alloc_mature:units = "fraction" ; + fates_seed_alloc_mature:long_name = "fraction of available carbon balance allocated to seeds in mature plants (adds to fates_seed_alloc)" ; + float fates_seed_dbh_repro_threshold(fates_pft) ; + fates_seed_dbh_repro_threshold:units = "cm" ; + fates_seed_dbh_repro_threshold:long_name = "the diameter (if any) where the plant will start extra clonal allocation to the seed pool (NOT USED YET)" ; + float fates_seed_decay_turnover(fates_pft) ; + fates_seed_decay_turnover:units = "1/yr" ; + fates_seed_decay_turnover:long_name = "turnover time for seeds with respect to germination" ; + float fates_seed_germination_timescale(fates_pft) ; + fates_seed_germination_timescale:units = "1/yr" ; + fates_seed_germination_timescale:long_name = "turnover time for seeds with respect to decay" ; + float fates_seed_rain(fates_pft) ; + fates_seed_rain:units = "KgC/m2/yr" ; + fates_seed_rain:long_name = "External seed rain from outside site (non-mass conserving)" ; + float fates_smpsc(fates_pft) ; + fates_smpsc:units = "mm" ; + fates_smpsc:long_name = "Soil water potential at full stomatal closure" ; + float fates_smpso(fates_pft) ; + fates_smpso:units = "mm" ; + fates_smpso:long_name = "Soil water potential at full stomatal opening" ; + float fates_taulnir(fates_pft) ; + fates_taulnir:units = "fraction" ; + fates_taulnir:long_name = "Leaf transmittance: near-IR" ; + float fates_taulvis(fates_pft) ; + fates_taulvis:units = "fraction" ; + fates_taulvis:long_name = "Leaf transmittance: visible" ; + float fates_tausnir(fates_pft) ; + fates_tausnir:units = "fraction" ; + fates_tausnir:long_name = "Stem transmittance: near-IR" ; + float fates_tausvis(fates_pft) ; + fates_tausvis:units = "fraction" ; + fates_tausvis:long_name = "Stem transmittance: visible" ; + float fates_trim_inc(fates_pft) ; + fates_trim_inc:units = "m2/m2" ; + fates_trim_inc:long_name = "Arbitrary incremental change in trimming function." ; + float fates_trim_limit(fates_pft) ; + fates_trim_limit:units = "m2/m2" ; + fates_trim_limit:long_name = "Arbitrary limit to reductions in leaf area with stress" ; + float fates_wood_density(fates_pft) ; + fates_wood_density:units = "g/cm3" ; + fates_wood_density:long_name = "mean density of woody tissue in plant" ; + float fates_woody(fates_pft) ; + fates_woody:units = "logical flag" ; + fates_woody:long_name = "Binary woody lifeform flag" ; + float fates_z0mr(fates_pft) ; + fates_z0mr:units = "unitless" ; + fates_z0mr:long_name = "Ratio of momentum roughness length to canopy top height" ; + float fates_base_mr_20(fates_scalar) ; + fates_base_mr_20:units = "gC/gN/s" ; + fates_base_mr_20:long_name = "Base maintenance respiration rate for plant tissues, using Ryan 1991" ; + float fates_bbopt_c3(fates_scalar) ; + fates_bbopt_c3:units = "umol H2O/m**2/s" ; + fates_bbopt_c3:long_name = "Ball-Berry minimum unstressed leaf conductance for C3" ; + float fates_bbopt_c4(fates_scalar) ; + fates_bbopt_c4:units = "umol H2O/m**2/s" ; + fates_bbopt_c4:long_name = "Ball-Berry minimum unstressed leaf conductance for C4" ; + float fates_canopy_closure_thresh(fates_scalar) ; + fates_canopy_closure_thresh:units = "unitless" ; + fates_canopy_closure_thresh:long_name = "tree canopy coverage at which crown area allometry changes from savanna to forest value" ; + float fates_cohort_fusion_tol(fates_scalar) ; + fates_cohort_fusion_tol:units = "unitless" ; + fates_cohort_fusion_tol:long_name = "minimum fraction in difference in dbh between cohorts" ; + float fates_comp_excln(fates_scalar) ; + fates_comp_excln:units = "none" ; + fates_comp_excln:long_name = "weighting factor (exponent on dbh) for canopy layer exclusion and promotion" ; + float fates_cwd_fcel(fates_scalar) ; + fates_cwd_fcel:units = "unitless" ; + fates_cwd_fcel:long_name = "Cellulose fraction for CWD" ; + float fates_cwd_flig(fates_scalar) ; + fates_cwd_flig:units = "unitless" ; + fates_cwd_flig:long_name = "Lignin fraction of coarse woody debris" ; + float fates_fire_nignitions(fates_scalar) ; + fates_fire_nignitions:units = "/m2 (?)" ; + fates_fire_nignitions:long_name = "number of daily ignitions (nfires = nignitions*FDI*area_scaling)" ; + float fates_hydr_kmax_rsurf(fates_scalar) ; + fates_hydr_kmax_rsurf:units = "kg water/m2 root area/Mpa/s" ; + fates_hydr_kmax_rsurf:long_name = "maximum conducitivity for unit root surface" ; + float fates_hydr_psi0(fates_scalar) ; + fates_hydr_psi0:units = "MPa" ; + fates_hydr_psi0:long_name = "sapwood water potential at saturation" ; + float fates_hydr_psicap(fates_scalar) ; + fates_hydr_psicap:units = "MPa" ; + fates_hydr_psicap:long_name = "sapwood water potential at which capillary reserves exhausted" ; + float fates_init_litter(fates_scalar) ; + fates_init_litter:units = "NA" ; + fates_init_litter:long_name = "Initialization value for litter pool in cold-start (NOT USED)" ; + float fates_logging_coll_under_frac(fates_scalar) ; + fates_logging_coll_under_frac:units = "fraction" ; + fates_logging_coll_under_frac:long_name = "Fraction of stems killed in the understory when logging generates disturbance" ; + float fates_logging_collateral_frac(fates_scalar) ; + fates_logging_collateral_frac:units = "fraction" ; + fates_logging_collateral_frac:long_name = "Fraction of large stems in upperstory that die from logging collateral damage" ; + float fates_logging_dbhmax_infra(fates_scalar) ; + fates_logging_dbhmax_infra:units = "cm" ; + fates_logging_dbhmax_infra:long_name = "Tree diameter, above which infrastructure from logging does not impact damage or mortality." ; + float fates_logging_dbhmin(fates_scalar) ; + fates_logging_dbhmin:units = "cm" ; + fates_logging_dbhmin:long_name = "Minimum dbh at which logging is applied" ; + float fates_logging_direct_frac(fates_scalar) ; + fates_logging_direct_frac:units = "fraction" ; + fates_logging_direct_frac:long_name = "Fraction of stems logged directly per event" ; + float fates_logging_event_code(fates_scalar) ; + fates_logging_event_code:units = "unitless" ; + fates_logging_event_code:long_name = "Integer code that options how logging events are structured" ; + float fates_logging_mechanical_frac(fates_scalar) ; + fates_logging_mechanical_frac:units = "fraction" ; + fates_logging_mechanical_frac:long_name = "Fraction of stems killed due infrastructure an other mechanical means" ; + float fates_mort_disturb_frac(fates_scalar) ; + fates_mort_disturb_frac:units = "fraction" ; + fates_mort_disturb_frac:long_name = "fraction of canopy mortality that results in disturbance (i.e. transfer of area from new to old patch)" ; + float fates_mort_understorey_death(fates_scalar) ; + fates_mort_understorey_death:units = "fraction" ; + fates_mort_understorey_death:long_name = "fraction of plants in understorey cohort impacted by overstorey tree-fall" ; + float fates_patch_fusion_tol(fates_scalar) ; + fates_patch_fusion_tol:units = "unitless" ; + fates_patch_fusion_tol:long_name = "minimum fraction in difference in profiles between patches" ; + float fates_phen_a(fates_scalar) ; + fates_phen_a:units = "none" ; + fates_phen_a:long_name = "GDD accumulation function, intercept parameter: gdd_thesh = a + b exp(c*ncd)" ; + float fates_phen_b(fates_scalar) ; + fates_phen_b:units = "none" ; + fates_phen_b:long_name = "GDD accumulation function, multiplier parameter: gdd_thesh = a + b exp(c*ncd)" ; + float fates_phen_c(fates_scalar) ; + fates_phen_c:units = "none" ; + fates_phen_c:long_name = "GDD accumulation function, exponent parameter: gdd_thesh = a + b exp(c*ncd)" ; + float fates_phen_chiltemp(fates_scalar) ; + fates_phen_chiltemp:units = "degrees C" ; + fates_phen_chiltemp:long_name = "chilling day counting threshold" ; + float fates_phen_coldtemp(fates_scalar) ; + fates_phen_coldtemp:units = "degrees C" ; + fates_phen_coldtemp:long_name = "temperature exceedance to flag a cold-day for temperature leaf drop" ; + float fates_phen_doff_time(fates_scalar) ; + fates_phen_doff_time:units = "days" ; + fates_phen_doff_time:long_name = "day threshold compared against days since leaves became off-allometry" ; + float fates_phen_drought_threshold(fates_scalar) ; + fates_phen_drought_threshold:units = "m3/m3" ; + fates_phen_drought_threshold:long_name = "liquid volume in soil layer, threashold for drought phenology" ; + float fates_phen_mindayson(fates_scalar) ; + fates_phen_mindayson:units = "days" ; + fates_phen_mindayson:long_name = "day threshold compared against days since leaves became on-allometry" ; + float fates_phen_ncolddayslim(fates_scalar) ; + fates_phen_ncolddayslim:units = "days" ; + fates_phen_ncolddayslim:long_name = "day threshold exceedance for temperature leaf-drop" ; + float fates_durat_slope ; + fates_durat_slope:units = "NA" ; + fates_durat_slope:long_name = "spitfire parameter, fire max duration slope, Equation 14 Thonicke et al 2010" ; + float fates_fdi_a ; + fates_fdi_a:units = "NA" ; + fates_fdi_a:long_name = "spitfire parameter (unknown) " ; + float fates_fdi_alpha ; + fates_fdi_alpha:units = "NA" ; + fates_fdi_alpha:long_name = "spitfire parameter, EQ 7 Venevsky et al. GCB 2002,(modified EQ 8 Thonicke et al. 2010) " ; + float fates_fdi_b ; + fates_fdi_b:units = "NA" ; + fates_fdi_b:long_name = "spitfire parameter (unknown) " ; + float fates_fuel_energy ; + fates_fuel_energy:units = "kJ/kg" ; + fates_fuel_energy:long_name = "pitfire parameter, heat content of fuel" ; + float fates_max_durat ; + fates_max_durat:units = "minutes" ; + fates_max_durat:long_name = "spitfire parameter, fire maximum duration, Equation 14 Thonicke et al 2010" ; + float fates_miner_damp ; + fates_miner_damp:units = "NA" ; + fates_miner_damp:long_name = "spitfire parameter, mineral-dampening coefficient EQ A1 Thonicke et al 2010 " ; + float fates_miner_total ; + fates_miner_total:units = "fraction" ; + fates_miner_total:long_name = "spitfire parameter, total mineral content, Table A1 Thonicke et al 2010" ; + float fates_part_dens ; + fates_part_dens:units = "kg/m2" ; + fates_part_dens:long_name = "spitfire parameter, oven dry particle density, Table A1 Thonicke et al 2010" ; + float fates_soil_salinity(fates_scalar) ; + fates_soil_salinity:units = "ppt" ; + fates_soil_salinity:long_name = "soil salinity used for model when not coupled to dynamic soil salinity" ; + +// global attributes: + :history = "This file was made from FatesPFTIndexSwapper.py \n", + " Input File = fates_params_13pfts.c180315.nc \n", + " Indices = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 13] \n", + " Wed Feb 28 13:43:44 PST 2018 Values from Jennifer Holms 13pft test file fates_params.c170929_13pfts.nc were then manually converted over into positions 2-13, position 1 kept the tropical broadleaf evergreen. Wed Mar 7 15:53:33 PST 2018 (RGK) added fates_logging_dbhmax_infra, fates_maintresp_reduction_curvature, fates_maintresp_reduction_intercept, and fates_leaf_clumping_index. CHanged fates_clone_alloc to fates_seed_alloc_mature. Updated minimum grass sizes per J Shuman notes. Set grass AGB intercepts and latosa to zero. Updated clumping to have more substanteated starter values per suggestions by Shawn Serbin.\n", + " Thu Mar 15 13:48:11 PDT 2018 Added C4 plants via suggestions by @huitang_earth.\n", + " Mon Mar 19 19:05:44 EDT 2018 forced all plants to be evergreen till carbon-imbalances are fixed.\n", + " Mon Mar 19 19:05:44 EDT 2018 added entry for fates_history_height_bin_edges.\n", + " Thu Apr 12 14:18:46 PDT 2018 updated names on CN recruitment parameters. -- fates_params_14pft.c180412.nc -- \n", + " This file was made then modified wih FatesPFTIndexSwapper.py to reduce to 2 tropical pfts again. \n", + " Indices = [1, 1]" ; +data: + + fates_history_sizeclass_bin_edges = 0, 0.05, 0.1, 0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, + 0.80, 0.90, 1.0 ; + + fates_history_ageclass_bin_edges = 0, 1, 2, 5, 10, 20, 50 ; + + fates_FBD = 4, 15.4, 16.8, 19.6, 999, 4 ; + + fates_SAV = 66, 13, 3.58, 0.98, 0.2, 66 ; + + fates_alpha_FMC = 0.0050769, 0.001, 0.0002754, 7.54e-05, 1.54e-05, 999 ; + + fates_history_height_bin_edges = 0, 0.1, 0.3, 0.5, 1, 2 ; + + fates_low_moisture_Coeff = 1.15, 1.12, 1.09, 0.98, 0.8, 1.15 ; + + fates_low_moisture_Slope = 0.62, 0.62, 0.72, 0.85, 0.8, 0.62 ; + + fates_max_decomp = 1, 0.52, 0.383, 0.383, 0.19, 999 ; + + fates_mid_moisture = 0.8, 0.72, 0.51, 0.38, 1, 0.8 ; + + fates_mid_moisture_Coeff = 3.2, 2.35, 1.47, 1.06, 0.8, 3.2 ; + + fates_mid_moisture_Slope = 3.2, 2.35, 1.47, 1.06, 0.8, 3.2 ; + + fates_min_moisture = 0.24, 0.18, 0.12, 0, 0, 0.24 ; + + fates_hydr_avuln_node = + 2, 2, + 2, 2, + 2, 2, + 2, 2 ; + + fates_hydr_epsil_node = + 12, 12, + 10, 10, + 10, 10, + 8, 8 ; + + fates_hydr_fcap_node = + 0, 0, + 0.08, 0.08, + 0.08, 0.08, + 0, 0 ; + + fates_hydr_kmax_node = + -999, -999, + 3, 3, + -999, -999, + -999, -999 ; + + fates_hydr_p50_node = + -2.25, -2.25, + -2.25, -2.25, + -2.25, -2.25, + -2.25, -2.25 ; + + fates_hydr_pinot_node = + -999, -999, + -999, -999, + -999, -999, + -999, -999 ; + + fates_hydr_pitlp_node = + -1.67, -1.67, + -1.4, -1.4, + -1.4, -1.4, + -1.2, -1.2 ; + + fates_hydr_resid_node = + 0.25, 0.25, + 0.325, 0.325, + 0.325, 0.325, + 0.15, 0.15 ; + + fates_hydr_thetas_node = + 0.65, 0.65, + 0.65, 0.65, + 0.65, 0.65, + 0.75, 0.75 ; + + fates_CWD_frac = 0.045, 0.075, 0.21, 0.67 ; + + fates_pftname = + "Spartina alterniflora (short form) ", + "Spartina alterniflora (tall form) " ; + + fates_rootprof_beta = + 0.976, 0.976, + _, _ ; + + fates_alloc_storage_cushion = 5.25, 3.75 ; + + fates_allom_agb1 = 0.001928, 0.001928; + + fates_allom_agb2 = 1.9492, 1.9492 ; + + fates_allom_agb3 = 0.0, 0.0 ; + + fates_allom_agb4 = 0.0, 0.0 ; + + fates_allom_agb_frac = 1.0, 1.0 ; + + fates_allom_amode = 1, 1; + + fates_allom_blca_expnt_diff = 0, 0 ; + + fates_allom_cmode = 1, 1 ; + + fates_allom_d2bl1 = 0.000964, 0.000964; + + fates_allom_d2bl2 = 1.9492, 1.9492; + + fates_allom_d2bl3 = 0.0, 0.0 ; + + fates_allom_d2ca_coefficient_max = 0.03, 0.03 ; + + fates_allom_d2ca_coefficient_min = 0.01, 0.01 ; + + fates_allom_d2h1 = 1.0, 1.0 ; + + fates_allom_d2h2 = 1.0, 1.0 ; + + fates_allom_d2h3 = -999.9, -999.9 ; + + fates_allom_dbh_maxheight = 1, 1 ; + + fates_allom_fmode = 1, 1 ; + + fates_allom_hmode = 3, 3 ; + + fates_allom_l2fr = 0.75, 1.25; + + fates_allom_latosa_int = 0.001, 0.001 ; + + fates_allom_latosa_slp = 0, 0 ; + + fates_allom_lmode = 1, 1 ; + + fates_allom_sai_scaler = 0.0012, 0.0012 ; + + fates_allom_smode = 1, 1 ; + + fates_allom_stmode = 1, 1 ; + + fates_allom_frbstor_repro = 1.0, 1.0; + + fates_branch_turnover = 50, 50 ; + + fates_c2b = 2, 2 ; + + fates_displar = 0.67, 0.67 ; + + fates_fire_alpha_SH = 0.2, 0.2 ; + + fates_fire_bark_scaler = 0.07, 0.07 ; + + fates_fire_crown_depth_frac = 0.5, 0.5 ; + + fates_fire_crown_kill = 0.775, 0.775 ; + + fates_fr_fcel = 0.5, 0.5 ; + + fates_fr_flab = 0.25, 0.25 ; + + fates_fr_flig = 0.25, 0.25 ; + + fates_froot_cn_ratio = 42, 42 ; + + fates_grperc = 0.11, 0.11 ; + + fates_hydr_avuln_gs = 2.5, 2.5 ; + + fates_hydr_p50_gs = -1.5, -1.5 ; + + fates_hydr_p_taper = 0.333, 0.333 ; + + fates_hydr_rfrac_stem = 0.625, 0.625 ; + + fates_hydr_rs2 = 0.0001, 0.0001 ; + + fates_hydr_srl = 25, 25 ; + + fates_leaf_BB_slope = 8, 8 ; + + fates_leaf_c3psn = 0, 0; + + fates_leaf_clumping_index = 0.85, 0.85 ; + + fates_leaf_cn_ratio = 30, 30 ; + + fates_leaf_diameter = 0.04, 0.04 ; + + fates_leaf_jmaxha = 43540, 43540 ; + + fates_leaf_jmaxhd = 152040, 152040 ; + + fates_leaf_jmaxse = 495, 495 ; + + fates_leaf_long = 1.0, 1.0 ; + + fates_leaf_slatop = 0.016, 0.016 ; + + fates_leaf_stor_priority = 0.8, 0.8 ; + + fates_leaf_tpuha = 53100, 53100 ; + + fates_leaf_tpuhd = 150650, 150650 ; + + fates_leaf_tpuse = 490, 490 ; + + fates_leaf_vcmax25top = 70, 40 ; + + fates_leaf_vcmaxha = 65330, 65330 ; + + fates_leaf_vcmaxhd = 149250, 149250 ; + + fates_leaf_vcmaxse = 485, 485 ; + + fates_leaf_xl = 0.1, 0.1 ; + + fates_lf_fcel = 0.5, 0.5 ; + + fates_lf_flab = 0.25, 0.25 ; + + fates_lf_flig = 0.25, 0.25 ; + + fates_maintresp_reduction_curvature = 0.01, 0.01 ; + + fates_maintresp_reduction_intercept = 1, 1 ; + + fates_mort_bmort = 0.9, 0.9 ; + + fates_mort_freezetol = -25, -25 ; + + fates_mort_hf_sm_threshold = 1e-06, 1e-06 ; + + fates_mort_scalar_coldstress = 3, 3 ; + + fates_mort_scalar_cstarvation = 0.6, 0.6 ; + + fates_mort_scalar_hydrfailure = 0.6, 0.6 ; + + fates_pft_used = 1, 0 ; + + fates_phen_evergreen = 0, 0; + + fates_phen_season_decid = 1, 1 ; + + fates_phen_stress_decid = 0, 0 ; + + fates_prescribed_mortality_canopy = 0.0194, 0.0194 ; + + fates_prescribed_mortality_understory = 0.025, 0.025 ; + + fates_prescribed_npp_canopy = 0.4, 0.4 ; + + fates_prescribed_npp_understory = 0.03125, 0.03125 ; + + fates_prescribed_recruitment = 0.02, 0.02 ; + + fates_recruit_hgt_min = 0.05, 0.05 ; + + fates_recruit_initd = 10000000, 0 ; + + fates_rholnir = 0.45, 0.45 ; + + fates_rholvis = 0.1, 0.1 ; + + fates_rhosnir = 0.39, 0.39 ; + + fates_rhosvis = 0.16, 0.16 ; + + fates_root_long = 0.5, 1 ; + + fates_roota_par = 7, 7 ; + + fates_rootb_par = 1, 1 ; + + fates_seed_alloc = 0.1, 0.1 ; + + fates_seed_alloc_mature = 0, 0 ; + + fates_seed_dbh_repro_threshold = 150, 150 ; + + fates_seed_decay_turnover = 0.51, 0.51 ; + + fates_seed_germination_timescale = 0.5, 0.5 ; + + fates_seed_rain = 0.28, 0.28 ; + + fates_smpsc = -255000, -255000 ; + + fates_smpso = -66000, -66000 ; + + fates_taulnir = 0.25, 0.25 ; + + fates_taulvis = 0.05, 0.05 ; + + fates_tausnir = 0.001, 0.001 ; + + fates_tausvis = 0.001, 0.001 ; + + fates_trim_inc = 0.03, 0.03 ; + + fates_trim_limit = 0.3, 0.3 ; + + fates_wood_density = 0.01, 0.01 ; + + fates_woody = 1, 1 ; + + fates_z0mr = 0.055, 0.055 ; + + fates_base_mr_20 = 2.52e-06 ; + + fates_bbopt_c3 = 10000 ; + + fates_bbopt_c4 = 40000 ; + + fates_canopy_closure_thresh = 0.8 ; + + fates_cohort_fusion_tol = 0.05 ; + + fates_comp_excln = 3 ; + + fates_cwd_fcel = 0.76 ; + + fates_cwd_flig = 0.24 ; + + fates_fire_nignitions = 15 ; + + fates_hydr_kmax_rsurf = 0.001; + + fates_hydr_psi0 = 0 ; + + fates_hydr_psicap = -0.6 ; + + fates_init_litter = 0.05 ; + + fates_logging_coll_under_frac = 0.55983 ; + + fates_logging_collateral_frac = 0.05 ; + + fates_logging_dbhmax_infra = 35 ; + + fates_logging_dbhmin = 50 ; + + fates_logging_direct_frac = 0.15 ; + + fates_logging_event_code = -30 ; + + fates_logging_mechanical_frac = 0.05 ; + + fates_mort_disturb_frac = 1 ; + + fates_mort_understorey_death = 0.09 ; + + fates_patch_fusion_tol = 0.05 ; + + fates_phen_a = -68 ; + + fates_phen_b = 638 ; + + fates_phen_c = -0.001 ; + + fates_phen_chiltemp = 5 ; + + fates_phen_coldtemp = 7.5 ; + + fates_phen_doff_time = 100 ; + + fates_phen_drought_threshold = 0.15 ; + + fates_phen_mindayson = 30 ; + + fates_phen_ncolddayslim = 5 ; + + fates_durat_slope = -11.06 ; + + fates_fdi_a = 17.62 ; + + fates_fdi_alpha = 0.00037 ; + + fates_fdi_b = 243.12 ; + + fates_fuel_energy = 18000 ; + + fates_max_durat = 240 ; + + fates_miner_damp = 0.41739 ; + + fates_miner_total = 0.055 ; + + fates_part_dens = 513 ; + + fates_soil_salinity = 0.0 ; +} diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 30a335707a..bed07aed11 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -171,6 +171,9 @@ variables: float fates_allom_stmode(fates_pft) ; fates_allom_stmode:units = "index" ; fates_allom_stmode:long_name = "storage allometry function index" ; + float fates_allom_frbstor_repro(fates_pft) ; + fates_allom_frbstor_repro:units = "fraction" ; + fates_allom_frbstor_repro:long_name = "fraction of bstore goes to reproduction after plant dies" ; float fates_branch_turnover(fates_pft) ; fates_branch_turnover:units = "yr-1" ; fates_branch_turnover:long_name = "turnover time of branches" ; @@ -258,7 +261,7 @@ variables: fates_hydr_rfrac_stem:units = "fraction" ; fates_hydr_rfrac_stem:long_name = "fraction of total tree resistance from troot to canopy" ; float fates_hydr_rs2(fates_pft) ; - fates_hydr_rs2:units = "mm" ; + fates_hydr_rs2:units = "m" ; fates_hydr_rs2:long_name = "absorbing root radius" ; float fates_hydr_srl(fates_pft) ; fates_hydr_srl:units = "m g-1" ; @@ -566,9 +569,6 @@ variables: float fates_fdi_b ; fates_fdi_b:units = "NA" ; fates_fdi_b:long_name = "spitfire parameter (unknown) " ; - float fates_fire_wind_max ; - fates_fire_wind_max:units = "m/min" ; - fates_fire_wind_max:long_name = "maximum wind speed expected by the fire model" ; float fates_fuel_energy ; fates_fuel_energy:units = "kJ/kg" ; fates_fuel_energy:long_name = "pitfire parameter, heat content of fuel" ; @@ -584,6 +584,9 @@ variables: float fates_part_dens ; fates_part_dens:units = "kg/m2" ; fates_part_dens:long_name = "spitfire parameter, oven dry particle density, Table A1 Thonicke et al 2010" ; + float fates_soil_salinity(fates_scalar) ; + fates_soil_salinity:units = "ppt" ; + fates_soil_salinity:long_name = "soil salinity used for model when not coupled to dynamic soil salinity" ; // global attributes: :history = "This file was made from FatesPFTIndexSwapper.py \n", @@ -657,10 +660,10 @@ data: -2.25, -2.25 ; fates_hydr_pinot_node = - -999, -999, - -999, -999, - -999, -999, - -999, -999 ; + -1.465984, -1.465984, + -1.228070, -1.228070, + -1.228070, -1.228070, + -1.043478, -1.043478 ; fates_hydr_pitlp_node = -1.67, -1.67, @@ -751,6 +754,8 @@ data: fates_allom_smode = 1, 1 ; fates_allom_stmode = 1, 1 ; + + fates_allom_frbstor_repro = 0.0, 0.0; fates_branch_turnover = 50, 50 ; @@ -1004,7 +1009,7 @@ data: fates_fire_nignitions = 15 ; - fates_hydr_kmax_rsurf = 0.001; + fates_hydr_kmax_rsurf = 20; fates_hydr_psi0 = 0 ; @@ -1058,8 +1063,6 @@ data: fates_fdi_b = 243.12 ; - fates_fire_wind_max = 45.718 ; - fates_fuel_energy = 18000 ; fates_max_durat = 240 ; @@ -1069,4 +1072,6 @@ data: fates_miner_total = 0.055 ; fates_part_dens = 513 ; + + fates_soil_salinity = 0.4 ; } diff --git a/tools/FatesPFTIndexSwapper.py b/tools/FatesPFTIndexSwapper.py index aa6b25eb73..3ef7deb044 100755 --- a/tools/FatesPFTIndexSwapper.py +++ b/tools/FatesPFTIndexSwapper.py @@ -22,6 +22,7 @@ # ======================================================================================= pft_dim_name = 'fates_pft' +prt_dim_name = 'fates_prt_organs' class timetype: @@ -161,19 +162,23 @@ def main(argv): # Idenfity if this variable has pft dimension pft_dim_found = -1 + prt_dim_found = -1 pft_dim_len = len(fp_in.variables.get(key).dimensions) for idim, name in enumerate(fp_in.variables.get(key).dimensions): # Manipulate data if(name==pft_dim_name): pft_dim_found = idim - + if(name==prt_dim_name): + prt_dim_found = idim + + # Copy over the input data # Tedious, but I have to permute through all combinations of dimension position if( pft_dim_len == 0 ): out_var = fp_out.createVariable(key,'f',(fp_in.variables.get(key).dimensions)) out_var.assignValue(float(fp_in.variables.get(key).data)) - elif(pft_dim_found==-1): + elif( (pft_dim_found==-1) & (prt_dim_found==-1) ): out_var = fp_out.createVariable(key,'f',(fp_in.variables.get(key).dimensions)) out_var[:] = in_var[:] elif( (pft_dim_found==0) & (pft_dim_len==1) ): # 1D fates_pft @@ -183,8 +188,9 @@ def main(argv): tmp_out[id] = fp_in.variables.get(key).data[ipft-1] out_var[:] = tmp_out - - elif( (pft_dim_found==1) & (pft_dim_len==2) ): # 2D hdyro_organ - fate_pft + # 2D hydro_organ - fates_pft + # or.. prt_organ - fates_pft + elif( (pft_dim_found==1) & (pft_dim_len==2) ): out_var = fp_out.createVariable(key,'f',(fp_in.variables.get(key).dimensions)) dim2_len = fp_in.dimensions.get(fp_in.variables.get(key).dimensions[0]) tmp_out = np.zeros([dim2_len,num_pft_out]) @@ -200,6 +206,12 @@ def main(argv): for id,ipft in enumerate(donor_pft_indices): out_var[id] = fp_in.variables.get(key).data[ipft-1] + + elif( (prt_dim_found==0) & (pft_dim_len==2) ): # fates_prt_organs - string_length + out_var = fp_out.createVariable(key,'c',(fp_in.variables.get(key).dimensions)) + out_var[:] = in_var[:] + + else: print('This variable has a dimensioning that we have not considered yet.') print('Please add this condition to the logic above this statement.')