diff --git a/src/biogeochem/CNBalanceCheckMod.F90 b/src/biogeochem/CNBalanceCheckMod.F90 index 7b88422843..acaac60347 100644 --- a/src/biogeochem/CNBalanceCheckMod.F90 +++ b/src/biogeochem/CNBalanceCheckMod.F90 @@ -204,6 +204,9 @@ subroutine BeginCNColumnBalance(this, bounds, num_soilc, filter_soilc, & c = filter_soilc(fc) col_begcb(c) = totcolc(c) + if(use_fates_bgc)then + col_begcb(c)= soilbiogeochem_carbonstate_inst%totsomc_col(c) + endif col_begnb(c) = totcoln(c) end do @@ -304,17 +307,18 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & if( col%is_fates(c) ) then - + col_endcb(c) = soilbiogeochem_carbonstate_inst%totsomc_col(c) s = clm_fates%f2hmap(ic)%hsites(c) + ! this balance check determines the column level balance of the soil pool + ! when FATES is on. col_cinputs = fates_litter_flux(c) - - ! calculate total column-level outputs - ! fates has already exported burn losses and fluxes to the atm - ! So they are irrelevant here + + ! calculate total column-level outputs from soil pool ! (gC/m2/s) total heterotrophic respiration + ! what about leaching? + col_coutputs = soilbiogeochem_carbonflux_inst%hr_col(c) - else ! calculate total column-level inputs @@ -349,13 +353,14 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & err_found = .true. err_index = c end if - if (abs(col_errcb(c)) > this%cwarning) then - write(iulog,*) 'cbalance warning at c =', c, col_errcb(c), col_endcb(c) + if (abs(col_errcb(c)) > this%cwarning .and. .not. use_fates_bgc) then + write(iulog,*) 'cbalance warning at c =', c, col_errcb(c), col_endcb(c) + write(iulog,*) 'flux, dstock',(col_cinputs - col_coutputs)*dt, (col_endcb(c) - col_begcb(c)) end if - end do ! end of columns loop +! end do ! end of columns loop - if (err_found) then + if (err_found .and. .not. use_fates_bgc) then c = err_index write(iulog,*)'column cbalance error = ', col_errcb(c), c write(iulog,*)'is fates column? = ', col%is_fates(c) @@ -381,8 +386,13 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & write(iulog,*)'hr = ',soilbiogeochem_carbonflux_inst%hr_col(c)*dt end if write(iulog,*)'-1*som_c_leached = ',som_c_leached(c)*dt - call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, msg=errMsg(sourcefile, __LINE__)) + if(.not. use_fates_bgc)then + call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, msg=errMsg(sourcefile, __LINE__)) + endif end if + end do ! end of columns loop + + ! Repeat error check at the gridcell level call c2g( bounds = bounds, & @@ -390,17 +400,25 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & garr = totgrcc(bounds%begg:bounds%endg), & c2l_scale_type = 'unity', & l2g_scale_type = 'unity') + call c2g( bounds = bounds, & carr = som_c_leached(bounds%begc:bounds%endc), & garr = som_c_leached_grc(bounds%begg:bounds%endg), & c2l_scale_type = 'unity', & l2g_scale_type = 'unity') + if(use_fates_bgc)then + call c2g( bounds = bounds, & + carr = soilbiogeochem_carbonflux_inst%fates_nbp_col(bounds%begc:bounds%endc), & + garr = soilbiogeochem_carbonflux_inst%fates_nbp_grc(bounds%begg:bounds%endg), & + c2l_scale_type = 'unity', & + l2g_scale_type = 'unity') + end if + err_found = .false. do g = bounds%begg, bounds%endg ! calculate gridcell-level carbon storage for mass conservation check - ! Notes: - ! totgrcc = totcolc = totc_p2c_col(c) + soilbiogeochem_cwdc_col(c) + soilbiogeochem_totlitc_col(c) + soilbiogeochem_totmicc_col(c) + soilbiogeochem_totsomc_col(c) + soilbiogeochem_ctrunc_col(c) + ! Notes: totgrcc = totcolc = totc_p2c_col(c) + soilbiogeochem_cwdc_col(c) + soilbiogeochem_totlitc_col(c) + soilbiogeochem_totmicc_col(c) + soilbiogeochem_totsomc_col(c) + soilbiogeochem_ctrunc_col(c) ! totc_p2c_col = totc_patch = totvegc_patch(p) + xsmrpool_patch(p) + ctrunc_patch(p) + cropseedc_deficit_patch(p) ! Not including seedc_grc in grc_begcb and grc_endcb because ! seedc_grc forms out of thin air, for now, and equals @@ -439,12 +457,18 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & else - ! Totally punt on this for now. We just don't track these gridscale variables yet (RGK) - grc_cinputs = 0._r8 - grc_endcb(g) = grc_begcb(g) - grc_coutputs = 0._r8 - grc_errcb(g) = 0._r8 + ! FATES gridcell level balance check. + ! Inputs are calculated from the fates_nbp + ! Outputs the same as for cn mode. + ! Total land carbon is informed by FATES carbon. + grc_endcb(g) = totgrcc(g) + tot_woodprod_grc(g) + cropprod1_grc(g) + grc_cinputs = soilbiogeochem_carbonflux_inst%fates_nbp_grc(g) + + grc_coutputs = - som_c_leached_grc(g) + + grc_errcb(g) = (grc_cinputs - grc_coutputs) * dt - & + (grc_endcb(g) - grc_begcb(g)) end if ! check for significant errors @@ -452,13 +476,12 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & err_found = .true. err_index = g end if - if (abs(grc_errcb(g)) > this%cwarning) then - write(iulog,*) 'cbalance warning at g =', g, grc_errcb(g), grc_endcb(g) + if (abs(grc_errcb(g)) > this%cwarning .and. .not. use_fates_bgc ) then + write(iulog,*) 'cbal warning:', g, grc_errcb(g), grc_endcb(g) end if - end do ! end of gridcell loop - if (err_found) then - g = err_index + if (err_found .and. .not. use_fates_bgc) then + !g = err_index write(iulog,*)'gridcell cbalance error =', grc_errcb(g), g write(iulog,*)'latdeg, londeg =', grc%latdeg(g), grc%londeg(g) write(iulog,*)'begcb =', grc_begcb(g) @@ -472,7 +495,10 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & write(iulog,*)'-1*som_c_leached_grc = ', som_c_leached_grc(g) * dt call endrun(subgrid_index=g, subgrid_level=subgrid_level_gridcell, msg=errMsg(sourcefile, __LINE__)) end if - + + ! moving the update until after the error report. + grc_endcb(g) = grc_begcb(g) + end do ! end associate end subroutine CBalanceCheck @@ -649,7 +675,7 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, & err_index = c end if - if (abs(col_errnb(c)) > this%nwarning) then + if (abs(col_errnb(c)) > this%nwarning ) then write(iulog,*) 'nbalance warning at c =', c, col_errnb(c), col_endnb(c) write(iulog,*)'inputs,ffix,nfix,ndep = ',ffix_to_sminn(c)*dt,nfix_to_sminn(c)*dt,ndep_to_sminn(c)*dt write(iulog,*)'outputs,lch,roff,dnit = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90 index 0274fcc87f..13990c1e2b 100644 --- a/src/biogeochem/CNDriverMod.F90 +++ b/src/biogeochem/CNDriverMod.F90 @@ -860,6 +860,7 @@ subroutine CNDriverNoLeaching(bounds, end if if_bgc_vegp2: if(num_bgc_vegp>0)then + call c_products_inst%UpdateProducts(bounds, & num_bgc_vegp, filter_bgc_vegp, & dwt_wood_product_gain_patch = cnveg_carbonflux_inst%dwt_wood_productc_gain_patch(begp:endp), & @@ -903,11 +904,13 @@ subroutine CNDriverNoLeaching(bounds, if (use_c14) call c14_products_inst%ComputeProductSummaryVars(bounds) call n_products_inst%ComputeProductSummaryVars(bounds) - call c_products_inst%ComputeSummaryVars(bounds) if (use_c13) call c13_products_inst%ComputeSummaryVars(bounds) if (use_c14) call c14_products_inst%ComputeSummaryVars(bounds) call n_products_inst%ComputeSummaryVars(bounds) + + ! unsure where else to put this - it can't go in the summary as that is all column variables. + soilbiogeochem_carbonflux_inst%fates_product_loss_grc(bounds%begg:bounds%endg)=c_products_inst%product_loss_grc(bounds%begg:bounds%endg) call t_stopf('CNWoodProducts') @@ -1261,7 +1264,7 @@ subroutine CNDriverSummarizeFluxes(bounds, & ! ---------------------------------------------- ! soilbiogeochem carbon/nitrogen flux summary ! ---------------------------------------------- - + call soilbiogeochem_carbonflux_inst%Summary(bounds, num_bgc_soilc, filter_bgc_soilc, & num_bgc_vegp, filter_bgc_vegp, & soilbiogeochem_carbonflux_inst%decomp_cascade_ctransfer_col(begc:endc,1:ndecomp_cascade_transitions), & diff --git a/src/biogeochem/CNVegetationFacade.F90 b/src/biogeochem/CNVegetationFacade.F90 index dd20ce50a2..d01d9b3906 100644 --- a/src/biogeochem/CNVegetationFacade.F90 +++ b/src/biogeochem/CNVegetationFacade.F90 @@ -250,6 +250,7 @@ subroutine Init(this, bounds, NLFilename, nskip_steps, params_ncid) call this%CNReadNML( NLFilename ) ! MUST be called first as passes down control information to others end if + if(use_cn.or.use_fates_bgc)then call this%cnveg_carbonstate_inst%Init(bounds, carbon_type='c12', ratio=1._r8, & NLFilename=NLFilename, dribble_crophrv_xsmrpool_2atm=this%dribble_crophrv_xsmrpool_2atm, & @@ -912,7 +913,7 @@ subroutine InitGridcellBalance(this, bounds, num_allc, filter_allc, & garr = soilbiogeochem_carbonstate_inst%totc_grc(bounds%begg:bounds%endg), & c2l_scale_type = 'unity', & l2g_scale_type = 'unity') - + ! total gridcell nitrogen (TOTGRIDCELLN) call c2g( bounds = bounds, & carr = soilbiogeochem_nitrogenstate_inst%totn_col(bounds%begc:bounds%endc), & @@ -1061,7 +1062,8 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, & soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & - soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst) + soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst,c_products_inst,& + clm_fates) ! ! !DESCRIPTION: ! Do the main science for CN vegetation that needs to be done after hydrology-drainage @@ -1100,7 +1102,9 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, & type(soilbiogeochem_carbonstate_type) , intent(inout) :: c14_soilbiogeochem_carbonstate_inst type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst - ! + type(cn_products_type) , intent(inout) :: c_products_inst + type(hlm_fates_interface_type) , intent(inout) :: clm_fates + ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'EcosystemDynamicsPostDrainage' @@ -1140,8 +1144,11 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, & c14_soilbiogeochem_carbonstate_inst,soilbiogeochem_nitrogenstate_inst) call t_stopf('SoilBiogeochemPrecisionControl') + if(use_fates_bgc)then + call clm_fates%wrap_AtmosphericCarbonFluxes(bounds,soilbiogeochem_carbonflux_inst,soilbiogeochem_carbonstate_inst) + endif + ! Call to all CN summary routines - call CNDriverSummarizeStates(bounds, & num_allc, filter_allc, & num_bgc_soilc, filter_bgc_soilc, & @@ -1172,6 +1179,8 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, & soilbiogeochem_nitrogenstate_inst, & soilbiogeochem_nitrogenflux_inst) + + ! On the radiation time step, use C state variables to calculate ! vegetation structure (LAI, SAI, height) if(num_bgc_vegp>0)then @@ -1181,6 +1190,7 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, & crop_inst, this%cnveg_carbonstate_inst, canopystate_inst) end if end if + end subroutine EcosystemDynamicsPostDrainage @@ -1217,6 +1227,7 @@ subroutine BalanceCheck(this, bounds, num_bgc_soilc, filter_bgc_soilc, & character(len=*), parameter :: subname = 'BalanceCheck' !----------------------------------------------------------------------- + DA_nstep = get_nstep_since_startup_or_lastDA_restart_or_pause() if (DA_nstep <= skip_steps )then if (masterproc) then diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index 7ed9f82586..ccebb88bbb 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -84,6 +84,7 @@ module clm_driver use clm_instMod use SoilMoistureStreamMod , only : PrescribedSoilMoistureInterp, PrescribedSoilMoistureAdvance use SoilBiogeochemDecompCascadeConType , only : no_soil_decomp, decomp_method + use CNProductsMod , only : cn_products_type ! ! !PUBLIC TYPES: implicit none @@ -128,6 +129,9 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro logical, intent(in) :: nlend ! true => end of run on this step character(len=*),intent(in) :: rdate ! restart file time stamp for name + !where are all the other instances defined? I can't figure this out! + type(cn_products_type) :: c_products_inst + ! Whether we're running with a prognostic ROF component. This shouldn't change from ! timestep to timestep, but we pass it into the driver loop because it isn't available ! in initialization. @@ -1041,6 +1045,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call t_stopf('ecosysdyn') end if + ! Prescribed biogeography - prescribed canopy structure, some prognostic carbon fluxes if (((.not. use_cn) .and. (.not. use_fates) .and. (doalb))) then @@ -1104,7 +1109,6 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call t_stopf('hydro2_drainage') - if (use_cn .or. use_fates_bgc) then call t_startf('EcosysDynPostDrainage') call bgc_vegetation_inst%EcosystemDynamicsPostDrainage(bounds_clump, & @@ -1120,7 +1124,8 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & - soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst) + soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst, c_products_inst,& + clm_fates) call t_stopf('EcosysDynPostDrainage') end if @@ -1148,7 +1153,9 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, & water_inst%wateratm2lndbulk_inst, canopystate_inst, soilbiogeochem_carbonflux_inst, & frictionvel_inst, soil_water_retention_curve) - + if(use_fates)then +! call clm_fates%wrap_AtmosphericCarbonFluxes(nc,bounds_proc,soilbiogeochem_carbonflux_inst,c_products_inst) + endif ! TODO(wjs, 2016-04-01) I think this setFilters call should be replaced by a ! call to reweight_wrapup, if it's needed at all. call setFilters( bounds_clump, glc_behavior ) @@ -1304,9 +1311,13 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro ! based on similar pgi compiler bugs that we have run into before). Also note that I ! don't have explicit bounds on the left-hand-side of this assignment: excluding these ! explicit bounds seemed to be needed to get around other compiler bugs. - allocate(net_carbon_exchange_grc(bounds_proc%begg:bounds_proc%endg)) - net_carbon_exchange_grc = bgc_vegetation_inst%get_net_carbon_exchange_grc(bounds_proc) + allocate(net_carbon_exchange_grc(bounds_proc%begg:bounds_proc%endg)) + if(use_fates_bgc)then + net_carbon_exchange_grc = soilbiogeochem_carbonflux_inst%fates_nbp_grc(bounds_proc%begg:bounds_proc%endg) + else + net_carbon_exchange_grc = bgc_vegetation_inst%get_net_carbon_exchange_grc(bounds_proc) + endif call lnd2atm(bounds_proc, & atm2lnd_inst, surfalb_inst, temperature_inst, frictionvel_inst, & water_inst, & diff --git a/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 b/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 index 23f24e44d5..4a5c7061b6 100644 --- a/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 @@ -6,7 +6,7 @@ module SoilBiogeochemCarbonFluxType use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, ndecomp_cascade_outtransitions use clm_varpar , only : nlevdecomp_full, nlevgrnd, nlevdecomp, nlevsoi, ndecomp_pools_vr, i_cwdl2 use clm_varcon , only : spval, ispval, dzsoi_decomp - use clm_varctl , only : use_fates,use_cn + use clm_varctl , only : use_fates,use_fates_bgc,use_cn use pftconMod , only : pftcon use landunit_varcon , only : istsoil, istcrop, istdlak use ch4varcon , only : allowlakeprod @@ -52,7 +52,17 @@ module SoilBiogeochemCarbonFluxType ! nitrif_denitrif real(r8), pointer :: phr_vr_col (:,:) ! (gC/m3/s) potential hr (not N-limited) real(r8), pointer :: fphr_col (:,:) ! fraction of potential heterotrophic respiration - + !-----FATES compposite fluxes ----------! + real(r8), pointer :: fates_npp_col (:) ! (gC/m2/s) net ecosystem exchange when FATES is on + real(r8), pointer :: fates_nee_col (:) ! (gC/m2/s) net ecosystem exchange when FATES is on + real(r8), pointer :: fates_nep_col (:) ! (gC/m2/s) net ecosystem productivity when FATES is on + + real(r8), pointer :: fates_nbp_col (:) ! (gC/m2/s) net biome productivity when FATES is on column level + real(r8), pointer :: fates_nbp_grc (:) ! (gC/m2/s) net biome productivity when FATES is on gridcell level + real(r8), pointer :: fates_fire_grazing_col (:) ! (gC/m2/s) fire plus grazing fluxes to the atmosphere (happen the -next- day as fire & grazing are calcualted at midnight. + real(r8), pointer :: fates_product_loss_grc (:) ! (gC/m2/s) total loss from product pools for calcualtion of fates_nbp + + ! ----- Hetertrophic Respiration fluxes --------! real(r8), pointer :: hr_col (:) ! (gC/m2/s) total heterotrophic respiration real(r8), pointer :: michr_col (:) ! (gC/m2/s) microbial heterotrophic respiration: donor-pool based definition, so expect it to be zero with MIMICS; microbial decomposition is responsible for heterotrophic respiration of donor pools (litter and soil), but in the accounting we assign it to the donor pool for consistency with CENTURY real(r8), pointer :: cwdhr_col (:) ! (gC/m2/s) coarse woody debris heterotrophic respiration: donor-pool based definition @@ -121,12 +131,14 @@ subroutine InitAllocate(this, bounds) ! !LOCAL VARIABLES: integer :: begp,endp ! Begin and end patch integer :: begc,endc ! Begin and end column + integer :: begg,endg integer :: Ntrans,Ntrans_diag ! N trans size for matrix solution !------------------------------------------------------------------------ begp = bounds%begp; endp = bounds%endp begc = bounds%begc; endc = bounds%endc - + begg = bounds%begg; endg = bounds%endg + allocate(this%t_scalar_col (begc:endc,1:nlevdecomp_full)); this%t_scalar_col (:,:) =spval allocate(this%w_scalar_col (begc:endc,1:nlevdecomp_full)); this%w_scalar_col (:,:) =spval allocate(this%o_scalar_col (begc:endc,1:nlevdecomp_full)); this%o_scalar_col (:,:) =spval @@ -171,6 +183,15 @@ subroutine InitAllocate(this, bounds) allocate(this%decomp_cpools_transport_tendency_col(begc:endc,1:nlevdecomp_full,1:ndecomp_pools)) this%decomp_cpools_transport_tendency_col(:,:,:)= nan + if(use_fates_bgc)then + allocate(this%fates_npp_col (begc:endc)) ; this%fates_npp_col (:) = nan + allocate(this%fates_nee_col (begc:endc)) ; this%fates_nee_col (:) = nan + allocate(this%fates_nep_col (begc:endc)) ; this%fates_nep_col (:) = nan + allocate(this%fates_nbp_col (begc:endc)) ; this%fates_nbp_col (:) = nan + allocate(this%fates_nbp_grc (begg:endg)) ; this%fates_nbp_grc (:) = nan + allocate(this%fates_fire_grazing_col (begc:endc)) ; this%fates_fire_grazing_col (:) = nan + allocate(this%fates_product_loss_grc (begg:endg)) ; this%fates_product_loss_grc (:) = nan + endif allocate(this%hr_col (begc:endc)) ; this%hr_col (:) = nan allocate(this%michr_col (begc:endc)) ; this%michr_col (:) = nan @@ -235,6 +256,7 @@ subroutine InitHistory(this, bounds, carbon_type) character(10) :: active integer :: begp,endp integer :: begc,endc + integer :: begg,endg character(24) :: fieldname character(100) :: longname real(r8), pointer :: data1dptr(:) ! temp. pointer for slicing larger arrays @@ -243,7 +265,8 @@ subroutine InitHistory(this, bounds, carbon_type) begp = bounds%begp; endp = bounds%endp begc = bounds%begc; endc = bounds%endc - + begg = bounds%begg; endg = bounds%endg + if (nlevdecomp > 1) then vr_suffix = "_vr" else @@ -255,7 +278,20 @@ subroutine InitHistory(this, bounds, carbon_type) !------------------------------- ! add history fields for all CLAMP CN variables - + if (use_fates_bgc)then + call hist_addfld1d (fname='FATES_NEE', units='gC/m^2/s', & + avgflag='A', long_name='FATES net ecosystem productivity', & + ptr_col=this%fates_nee_col) + + call hist_addfld1d (fname='FATES_NBP', units='gC/m^2/s', & + avgflag='A', long_name='FATES net biome productivity', & + ptr_gcell=this%fates_nbp_grc) + ! this is writing the COLUMN level NBP, but we want to write out the gridcell fluxes. + + call hist_addfld1d (fname='FATES_FIRE_GRAZING', units='gC/m^2/s', & + avgflag='A', long_name='FATES fire and grazing fluxes (NBP-NEP) ', & + ptr_col=this%fates_fire_grazing_col) + endif if (carbon_type == 'c12') then this%hr_col(begc:endc) = spval @@ -306,7 +342,18 @@ subroutine InitHistory(this, bounds, carbon_type) avgflag='A', long_name=longname, & ptr_col=data2dptr, default='inactive') end do - + + if(use_fates_bgc)then + this%fates_npp_col(begc:endc) = spval + this%fates_nee_col(begc:endc) = spval + this%fates_nep_col(begc:endc) = spval + this%fates_nbp_col(begc:endc) = spval + this%fates_fire_grazing_col(begc:endc) = spval + this%fates_nbp_grc(begg:endg) = spval + this%fates_product_loss_grc(begg:endg) = spval + endif + + this%decomp_cascade_hr_col(begc:endc,:) = spval this%decomp_cascade_hr_vr_col(begc:endc,:,:) = spval this%decomp_cascade_ctransfer_col(begc:endc,:) = spval @@ -787,6 +834,13 @@ subroutine SetValues ( this, num_column, filter_column, value_column) do fi = 1,num_column i = filter_column(fi) + if(use_fates_bgc)then + this%fates_npp_col(i) = value_column + this%fates_nee_col(i) = value_column + this%fates_nep_col(i) = value_column + this%fates_nbp_col(i) = value_column + this%fates_fire_grazing_col(i) = value_column + endif this%hr_col(i) = value_column this%somc_fire_col(i) = value_column this%som_c_leached_col(i) = value_column @@ -795,8 +849,9 @@ subroutine SetValues ( this, num_column, filter_column, value_column) this%cwdhr_col(i) = value_column this%michr_col(i) = value_column this%soilc_change_col(i) = value_column - end do + end do + end subroutine SetValues !----------------------------------------------------------------------- @@ -810,7 +865,7 @@ subroutine Summary(this, bounds, & ! On the radiation time step, carbon summary calculations ! ! !USES: - use subgridAveMod, only: p2c + use subgridAveMod, only: p2c , c2g use CNSharedParamsMod, only: CNParamsShareInst ! !ARGUMENTS: @@ -824,11 +879,12 @@ subroutine Summary(this, bounds, & real(r8), intent(in), optional :: soilbiogeochem_cwdn_col(bounds%begc:) real(r8), intent(in), optional :: soilbiogeochem_decomp_cascade_ctransfer_col(bounds%begc:,1:) + real(r8), intent(in), optional :: leafc_to_litter_patch(:) real(r8), intent(in), optional :: frootc_to_litter_patch(:) ! ! !LOCAL VARIABLES: - integer :: c,j,k,l,p + integer :: c,j,k,l,p,g integer :: fc, fp real(r8) :: ligninNratio_cwd ! lignin to N ratio of CWD real(r8) :: ligninNratio_leaf_patch(bounds%begp:bounds%endp) ! lignin to N ratio of leaves, patch level @@ -1007,14 +1063,40 @@ subroutine Summary(this, bounds, & max(1.0e-3_r8, leafc_to_litter_col(c) + & frootc_to_litter_col(c) + & soilbiogeochem_decomp_cascade_ctransfer_col(c,i_cwdl2)) - !else - ! For FATES: - ! this array is currently updated here: - ! clmfates_interfaceMod.F90:wrap_update_hlmfates_dyn() + else + ! fates does not work with MIMICS yet! end if end do end if if_mimics + + do fc = 1,num_bgc_soilc + c = filter_bgc_soilc(fc) + + ! For FATES, we update the coposite flux pools here with the infomationwe + ! recevied from FATES on column level NPP and grazing/fire fluxes. + if(col%is_fates(c)) then + ! These are all gc/,2/s and instantaneous. + this%fates_nep_col(c) = this%fates_npp_col(c) - this%hr_col(c) + this%fates_nbp_col(c) = this%fates_nep_col(c) - this%fates_fire_grazing_col(c) + ! do not add product loss here as it is a gridcell variable + + end if + end do + + ! create the gridcell level NBP variable + call c2g( bounds = bounds, & + carr = this%fates_nbp_col(bounds%begc:bounds%endc), & + garr = this%fates_nbp_grc(bounds%begg:bounds%endg), & + c2l_scale_type = 'unity', & + l2g_scale_type = 'unity') + + ! remove the gridcell level product_loss flux are the gridcell level. + ! this is a FATES variable that is passed into this module after the product pool summary. + do g = bounds%begg,bounds%endg + this%fates_nbp_grc(g) = this%fates_nbp_grc(g) - this%fates_product_loss_grc(g) + enddo + end subroutine Summary end module SoilBiogeochemCarbonFluxType diff --git a/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 b/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 index 768811925d..da65e1f012 100644 --- a/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 @@ -1,5 +1,5 @@ module SoilBiogeochemCarbonStateType - + use shr_kind_mod , only : r8 => shr_kind_r8 use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use shr_log_mod , only : errMsg => shr_log_errMsg @@ -53,8 +53,8 @@ module SoilBiogeochemCarbonStateType real(r8), pointer :: totc_col (:) ! (gC/m2) total column carbon, incl veg and cpool real(r8), pointer :: totecosysc_col (:) ! (gC/m2) total ecosystem carbon, incl veg but excl cpool real(r8), pointer :: totc_grc (:) ! (gC/m2) total gridcell carbon - - + real(r8), pointer :: fates_total_carbon_col (:) ! (gC/m2) total carbon contain in FATES state variables (vegetation, seeds, litter) for gridcell balance check. + ! Matrix-cn real(r8), pointer :: matrix_cap_decomp_cpools_col (:,:) ! (gC/m2) C capacity in decomposing (litter, cwd, soil) N pools in dimension (col,npools) real(r8), pointer :: matrix_cap_decomp_cpools_vr_col (:,:,:) ! (gC/m3) vertically-resolved C capacity in decomposing (litter, cwd, soil) pools in dimension(col,nlev,npools) @@ -177,7 +177,10 @@ subroutine InitAllocate(this, bounds) allocate(this%totc_col (begc:endc)) ; this%totc_col (:) = nan allocate(this%totecosysc_col (begc:endc)) ; this%totecosysc_col (:) = nan allocate(this%totc_grc (begg:endg)) ; this%totc_grc (:) = nan - + if(use_fates_bgc)then + allocate(this%fates_total_carbon_col (begc:endc)) ; this%fates_total_carbon_col (:) = nan + endif + this%restart_file_spinup_state = huge(1) end subroutine InitAllocate @@ -288,6 +291,12 @@ subroutine InitHistory(this, bounds, carbon_type) avgflag='A', long_name='total soil organic matter carbon', & ptr_col=this%totsomc_col) + if(use_fates_bgc)then + call hist_addfld1d (fname='FATES_TOTAL_CARBON', units='gC/m^2', & + avgflag='A', long_name='FATES total carbon stock (biomass + litter + seeds)', & + ptr_col=this%fates_total_carbon_col) + endif + if ( nlevdecomp_full > 1 ) then this%totlitc_1m_col(begc:endc) = spval call hist_addfld1d (fname='TOTLITC_1m', units='gC/m^2', & @@ -686,6 +695,10 @@ subroutine InitCold(this, bounds, ratio, c12_soilbiogeochem_carbonstate_inst) this%totc_col(c) = 0._r8 this%totecosysc_col(c) = 0._r8 end if + + if(use_fates_bgc)then + this%fates_total_carbon_col(c) = 0._r8 + endif end do @@ -1330,6 +1343,8 @@ subroutine SetValues ( this, num_column, filter_column, value_column) i = filter_column(fi) if ( .not. col%is_fates(i) ) then this%cwdc_col(i) = value_column + else + this%fates_total_carbon_col(i) = value_column end if this%ctrunc_col(i) = value_column this%totmicc_col(i) = value_column @@ -1609,6 +1624,7 @@ subroutine Summary(this, bounds, num_allc, filter_allc, num_bgc_soilc, filter_bg else num_local = num_allc end if + do fc = 1,num_local if(use_fates_bgc) then c = filter_bgc_soilc(fc) @@ -1616,9 +1632,13 @@ subroutine Summary(this, bounds, num_allc, filter_allc, num_bgc_soilc, filter_bg c = filter_allc(fc) end if if(col%is_fates(c)) then - totvegc_col = 0._r8 - ecovegc_col = 0._r8 + totvegc_col = this%fates_total_carbon_col(c) + ecovegc_col = this%fates_total_carbon_col(c) else + if(use_fates_bgc)then + totvegc_col = 0._r8 + ecovegc_col = 0._r8 + endif do l = 1, ndecomp_pools if ( decomp_cascade_con%is_cwd(l) ) then this%cwdc_col(c) = this%cwdc_col(c) + this%decomp_cpools_col(c,l) @@ -1627,7 +1647,6 @@ subroutine Summary(this, bounds, num_allc, filter_allc, num_bgc_soilc, filter_bg totvegc_col = cnveg_carbonstate_inst%totc_p2c_col(c) ecovegc_col = cnveg_carbonstate_inst%totvegc_col(c) end if - ! total ecosystem carbon, including veg but excluding cpool (TOTECOSYSC) this%totecosysc_col(c) = & this%cwdc_col(c) + & diff --git a/src/utils/clmfates_interfaceMod.F90 b/src/utils/clmfates_interfaceMod.F90 index efeee5e1f0..d74c7b77b8 100644 --- a/src/utils/clmfates_interfaceMod.F90 +++ b/src/utils/clmfates_interfaceMod.F90 @@ -1,3 +1,4 @@ + module CLMFatesInterfaceMod ! ------------------------------------------------------------------------------------- @@ -95,6 +96,7 @@ module CLMFatesInterfaceMod use atm2lndType , only : atm2lnd_type use SurfaceAlbedoType , only : surfalb_type use SolarAbsorbedType , only : solarabs_type + use CNVegCarbonFluxType , only : cnveg_carbonflux_type use SoilBiogeochemCarbonFluxType, only : soilbiogeochem_carbonflux_type use SoilBiogeochemCarbonStateType, only : soilbiogeochem_carbonstate_type use SoilBiogeochemNitrogenFluxType, only : soilbiogeochem_nitrogenflux_type @@ -255,6 +257,7 @@ module CLMFatesInterfaceMod procedure, public :: wrap_btran procedure, public :: wrap_drydep procedure, public :: wrap_photosynthesis + procedure, public :: wrap_atmosphericCarbonFluxes procedure, public :: wrap_accumulatefluxes procedure, public :: prep_canopyfluxes procedure, public :: wrap_canopy_radiation @@ -2917,7 +2920,64 @@ subroutine wrap_WoodProducts(this, bounds_clump, num_soilc, filter_soilc, & return end subroutine wrap_WoodProducts - ! ====================================================================================== +! ====================================================================================== + + subroutine wrap_atmosphericCarbonFluxes(this,bounds_clump,soilbiogeochem_carbonflux_inst,soilbiogeochem_carbonstate_inst) + + ! summarize the high-level fluxes that integrate information from both + ! FATES and outside-of-FATES decomposition and product decay code. + use subgridAveMod , only: c2g + use FatesConstantsMod , only : g_per_kg + + ! !ARGUMENTS: + type(bounds_type), intent(in) :: bounds_clump + class(hlm_fates_interface_type), intent(inout) :: this + integer :: nc + type(soilbiogeochem_carbonflux_type), intent(in) :: soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst + integer :: g,s,c + + ! NPP, grazing and fire and total carbon are the outputs from FATES that we pass to the HLM here + ! The composite fluxes are calculated using the HR and Product pools flues from the HLM in the + ! summary routines (which happen later). + ! Total carbon is needed for gridcell balance check purposes. + ! All fluxes here are in + + associate(& + fates_npp => soilbiogeochem_carbonflux_inst%fates_npp_col , & + fire_grazing => soilbiogeochem_carbonflux_inst%fates_fire_grazing_col , & + product_closs => soilbiogeochem_carbonflux_inst%fates_product_loss_grc , & + fates_total_carbon => soilbiogeochem_carbonstate_inst%fates_total_carbon_col) + + nc = bounds_clump%clump_index + + do s = 1, this%fates(nc)%nsites + c = this%f2hmap(nc)%fcolumn(s) + g = col%gridcell(c) + + ! Instantaneous site level NPP is calculate in FATES in AccumulateFluxes_ED + ! NEP, NBP are calculated in the summary after hr_col is determineds. + fates_npp(c) = this%fates(nc)%bc_out(s)%npp_site + + ! Convert yesteday's dribbled fluxes from fire and grazing from kgC/m2/s to g/m2/s + fire_grazing(c) = this%fates(nc)%bc_out(s)%grazing_closs_to_atm_si*g_per_kg & + + this%fates(nc)%bc_out(s)%fire_closs_to_atm_si*g_per_kg + + ! Total carbon in all FATES-side pools (biomass, litter & seeds). + fates_total_carbon(c) = this%fates(nc)%bc_out(s)%fates_total_carbon_site ! gC/m2 + + ! Add the instantaneous amount of carbon in the accumulated NPP pool, which at this + ! model timestep has not been allocated to a FATES biomass pool + ! (but will be at the end of the day) + + end do + + end associate + return + end subroutine wrap_atmosphericCarbonFluxes + + + ! ====================================================================================== subroutine wrap_canopy_radiation(this, bounds_clump, nc, fcansno, surfalb_inst) @@ -3528,6 +3588,7 @@ subroutine init_history_io(this,bounds_proc) select case(trim(ioname)) case(site_r8) + call hist_addfld1d(fname=trim(vname),units=trim(vunits), & avgflag=trim(vavgflag),long_name=trim(vlong), & ptr_col=fates_hist%hvars(ivar)%r81d, & @@ -3632,10 +3693,10 @@ subroutine ComputeRootSoilFlux(this, bounds_clump, num_filterc, filterc, & integer :: s integer :: c integer :: l - integer :: nc + integer :: n integer :: num_filter_fates integer :: nlevsoil - + integer :: nc if( .not. use_fates_planthydro ) return