search

ESCOMP / CTSM

Community Terrestrial Systems Model (includes the Community Land Model of CESM)
http://www.cesm.ucar.edu/models/cesm2.0/land/
Other
301 stars 306 forks source link

Perform gridcell-level carbon and nitrogen balance checks bracketing the entire run loop #314

Closed billsacks closed 4 years ago

billsacks commented 6 years ago

This is similar to #201 but for carbon and nitrogen rather than water and energy.

To summarize, we want an ordering like this:

The gridcell-level balance checks need to include some extra terms that exist at the gridcell level rather than the column-level:

We should try to introduce the gridcell-level checks without having a lot of duplication between these and the col-level balance checks. To do this, probably save shared fluxes (i.e., inputs & outputs) in the balance check object... then can compute col-level inputs & outputs for the col-level check, and average them to gridcell level for the grc-level check.

billsacks commented 6 years ago

cc @bishtgautam

bishtgautam commented 6 years ago

In CStateUpdateDynPatch(), isn't dzsoi_decomp(j) missing when updating of decomp_cpools_vr_col? The updated code should be:

cs_soil%decomp_cpools_vr_col(c,j,i_met_lit) = &
cs_soil%decomp_cpools_vr_col(c,j,i_met_lit) + &
cf_veg%dwt_frootc_to_litr_met_c_col(c,j) * dt / dzsoi_decomp(j)
billsacks commented 6 years ago

@bishtgautam Sorry, it's been a while since I've looked at this. Can you please describe why you think this should be changed?

bishtgautam commented 6 years ago

When trying to understand how updates to decomp_cpools_vr_col(:,:,:) modify the soil carbon pool totals, I called %Summary() before and after %dwt_* fluxes were accounted for. I found that if I don't account for dzsoi_decomp(j), the change reported by %Summary() does not match the sum of dwt fluxes. The %Summary() uses decomp_cpools_col(c,j) += decomp_cpools_vr_col(c,j,l) * dzsoi_decomp(j)

Additionally, the tridiagonal solve in SoilBiogeochemLittVertTranspMod.F90#L412 uses the source term contribution as source * dzsoi_decomp /dtime (i.e. accounts fordzsoi_decomp).

(Note: I used ELM and not CTSM, but the two models are essentially performing the same computation)

billsacks commented 6 years ago

There is very similar update code in CNCStateUpdate2Mod.F90 and CNCStateUpdate3Mod.F90, and none of that update code uses dzsoi_decomp. I'm guessing that I based the code in CStateUpdateDynPatch on that existing code. Can you describe why you think the code in CStateUpdateDynPatch needs the dzsoi_decomp factor even though the code in CNCStateUpdate2Mod and CNCStateUpdate3Mod does not have this factor?

I think of dzsoi_decomp as being needed when you're computing vertical integrals, as in the case of the Summary routine - though I realize that doesn't explain why it's needed in the tridiagonal solve.

bishtgautam commented 6 years ago

You are correct that CNCStateUpdate2Mod.F90 and CNCStateUpdate3Mod.F90 use a similar equation to update decomp_cpools_vr_col without accounting for dzsoi_decomp. I'm probably interpreting things incorrectly and need to understand the code better.

billsacks commented 6 years ago

@bishtgautam it's certainly possible that there is an issue here: I haven't looked back closely enough to really be able to defend the current implementation. So please feel free to add more if you think there's truly an issue.

bishtgautam commented 6 years ago

Here is an update on my first attempt at adding a grid-level nutrient balance check before and after calling dynSubgrid_driver() in ELM.

subroutine clm_driver()
...
do nc = 1,nclumps
   call Cflux_vars%ZeroDWT()
   call Nflux_vars%ZeroDWT()
   call Pflux_vars%ZeroDWT()
   call Cstate_vars%Summary()
   call Nstate_vars%Summary()
   call Pstate_vars%Summary()
   call BeginGridCBalanceBeforeDynSubgridDriver()
   call BeginGridNBalanceBeforeDynSubgridDriver()
   call BeginGridPBalanceBeforeDynSubgridDriver()
enddo

call dynSubgrid_driver()

do nc = 1,nclumps
   call Cstate_vars%Summary()
   call Nstate_vars%Summary()
   call Pstate_vars%Summary()
   call EndGridCBalanceAfterDynSubgridDriver()
   call EndGridNBalanceAfterDynSubgridDriver()
   call EndGridPBalanceAfterDynSubgridDriver()
enddo
end subroutine clm_driver
subroutine BeginGrid<CNP>BalanceBeforeDynSubgridDriver()
   call c2g(totcolc_col, begbal_grc)
end subroutine BeginGrid<CNP>BalanceBeforeDynSubgridDriver
subroutine EndGrid<CNP>BalanceAfterDynSubgridDriver()
   call c2g(totcolc_col, endbal_grc)

   do g = begg,endg
      grc_pinputs = &
         dwt_seedp_to_leaf_grc(g)     + &
        dwt_seedp_to_deadstem_grc(g)

      grc_poutputs = &
         dwt_conv_pflux_grc(g)        + &
         dwt_prod10p_gain_grc(g)      + &
         dwt_prod100p_gain_grc(g)

       errpb_grc(g) = (grc_pinputs - grc_poutputs)*dt - (endbal_grc(g) - begbal_grc(g))
   enddo
end subroutine EndGrid<CNP>BalanceAfterDynSubgridDriver

Note: At the moment, I haven't accounted for this%cropseedc_deficit_patch in ELM's pseudo-code because ELM presently only supports use_crop = .false.

bishtgautam commented 6 years ago

Forgot to mention in the previous comment that ELM does not use dribblers and all changes are applied instantaneously.

bishtgautam commented 6 years ago

DynamicColumnAdjustments() is also not included yet in ELM.

billsacks commented 5 years ago

Tagging this high priority to avoid further problems like #675 .

billsacks commented 5 years ago

It looks like @bishtgautam has put in place these grid cell level checks in e3sm master.

@bishtgautam Can you say why you ended up doing the gridcell-level CN balance checks immediately after the dyn subgrid stuff rather than near the end of the time loop (as is done for the gridcell-level water balance checks)? Do you remember if there was a problem with the latter for CN balance checks?

bishtgautam commented 5 years ago

@billsacks, The pseudo code I mentioned in https://github.com/ESCOMP/ctsm/issues/314#issuecomment-384101074 was only to perform grid-level checks before and after dyn subgrid. The aim was to check if dyn subgrid is conserving CN at grid-level.

Here is a possible extension of my current code to introduce new BeginGridCBalance()/EndGridCBalance() that check the conservation of Carbon at grid-level within the a CTSM time step

  1. BeginGridCBalance() is equivalent to calling BeginGridCBalanceBeforeDynSubgridDriver()
  2. EndGridCBalance(), which is called at the very end of the time loop: 2.1. Aggregates inputs/outputs for C at grid-level due to dyn subgrid driver (i.e. the code I added in EndGridCBalanceAfterDynSubgridDriver()) 2.2. Aggregates inputs/outputs for C at grid-level after dyn subgrid driver (i.e. the code in CBalanceCheck()#L185.
billsacks commented 5 years ago

Thanks @bishtgautam !

Also, a suggestion from @dlawrenncar - If the performance cost of the column-level checks is significant, we could do something like: Initialize both the gridcell and column-level balance checks. Then, at the end of the time step, just do the gridcell-level check; only do the column-level check if this gridcell-level check fails. The rationale is: the gridcell-level check should probably be sufficient to catch balance errors; the main point of the column-level check would then be to give additional information about the column with the problem.

Whether or not we do that, I realized that: We should either do the column-level checks first or, if we do the gridcell-level checks first (e.g., as per @dlawrenncar 's suggestion), then we should wait to abort until after we've done the column-level checks, so that we print information from the column-level checks.