Return bottom flux #377
Return bottom flux #377
Conversation
Currently this array is just set to zero, but eventually terms that are added to the kmt level of interior_tendencies(:,:) will be returned to the GCM as a flux instead
Instead, these terms are converted to fluxes and stored in bottom_fluxes; sed_denitrif affects no3, alk, and alk_alt_co2 while other_remin affects o2. For testing in the stand-alone test suite, I then convert these bottom_fluxes back to tendencies and add them to the interior_tendencies array. This change led to an update in the Jint_Ntot conservation computation as well (since sed_denitrif is not included in the tendency, it does not need to be added back into the total budget). I have verified that this introduces round-off level changes in Jint_Ntot and some of the tendency diagnostics, as well as a larger-than-roundoff change in o2_consumption (since this term no longer includes sed_denitrif and other_remin). I suspect the CI will fail because of this change to o2_consumption, but I'll wait until I have properly computed the rest of the bottom_flux terms before updating the baselines.
|
The CI failure in 75a1998 is expected: Details from the failing test: The change in |
This only modifies the stand-alone driver; as with the interior tendencies and the surface fluxes, it is up to the GCM to add these fields to output.
|
After adding bottom fluxes to driver output in 963df13, I figured the easiest way to test this would be to write a python script... I don't think it's necessary to show the full script, but the important lines are o2_diff = (ds.O2_CONSUMPTION - ds_baseline.O2_CONSUMPTION).isel(num_cols=colid, num_levels=kmt).values
btf_term = (ds.BTF_O2.isel(num_cols=colid)/(100*ds_grid.delta_z.isel(zt=kmt))).valuesAnd the output does show that the changes in |
This commit results in large changes in PON_REMIN_NH4 and PON_REMIN_DONr, but I have confirmed that those changes are offset (within machine precision) by the bottom flux values. This also results in an update to the Ntot conservation check.
|
For 838c8de, I expect As before, these are the only test failures and I have added those variables to my script to verify the changes are offset by the new |
This commit results in large changes in CaCO3_REMIN and CaCO3_ALT_CO2_REMIN that are offset (within machine precision) by the bottom flux values. This also results in an update to the Ctot conservation check.
Rather than change the kmt layer of the term we are integrating, I adjust the kmt level of the terms being summed for the comparison check. I think it looks a little cleaner, and it avoids dividing by delta_z(kmt) when adding a term to work (which is then multiplied by delta_z).
|
For 5fa8cec, we add From my tool's output (switching to column 2, because that's the only one with non-zero changes in I also verified that the absolute difference in So these failures are all still expected |
This commit results in a large change in SiO2_REMIN that is offset (within machine precision) by the bottom flux value. This also results in an update to the Sitot conservation check.
|
Continuing the pattern, 227b0cc adds But I am convinced this change is being offset by the |
This commit removes the bottom flux terms from POC%remin, which affects the bottom flux for the above three variables. This commit has been harder to verify because BTF_O2 and BTF_DIC were already non-zero -- the changes to O2_CONSUMPTION and POC_REMIN_DIC were not completely offset by bottom fluxes. I was able to compare POC_REMIN_DIC to output from the previous commit, but I'm not sure how to convince myself that O2_CONSUMPTION is still correct.
This commit results in large changes to POP_REMIN_PO4 and POP_REMIN_DOPr that are offset (within machine precision) by the bottom flux values. This also results in an update to the Ptot conservation check.
Despite only using the variable when POC%to_floor is positive, it was computed in one (POC%to_floor > 0) block and used in another, which raised a warning in some versions of gfortran. Computing it at the top of the (k == kmt) block should avoid that warning.
I had missed the relationship between lig_loss and POC_remin in previous commits
|
I think the bulk of this PR is done, but it looks like there is some trouble-shooting to do. Here's a full list of the errors from Travis: I don't know what to make of I'll look through the changes and see if anything catches my eye. Once I get to the bottom of this, I'll comment on whether the rest of the differences in the variables are expected. |
This was causing a problem in my netcdf comparison tool, because it was treating the max value of some fields as 1e37 and reporting a tiny relative error when in fact the fields were changing quite a bit.
Still need to figure out why J_O2 is off by so much
With this commit, J_O2 is no longer significantly different from the baseline
|
As of 4a5ab9d the variables that fail are all expected: I've updated my script to verify these changes are all reflected in So I think this is ready for review. @klindsay28, I'll send you a calendar invite -- besides going over the Fortran changes, I can also walk you through the python script that generated the output above; while the comparison is clear in most variables, the terms marked "offset due to changes in BTF in the new model" require multiplying by some scale factors and / or taking the difference of terms. |
|
It doesn't seem like ciso tracers are being handled. |
Unfortunately we don't have input data to run the call_compute_subroutines test with CISO enabled, so this will need to be tested in POP... but I think I updated the correct bottom flux terms and accounted for them correctly when doing the conservation checks in the diagnostics.
Good catch! I forgot about them because the datasets for our stand-alone test doesn't include CISO data. I think e6e3930 handles those tracers; I'll test it out in POP tomorrow. The POP test will abort if I botched the conservation checks, but I don't have a good way to test for correctness... maybe I should run POP for a day with nstep output? I could compare results from before and after pulling down that last commit. |
|
In subroutine |
Good find! I think this will be tricky to implement, though, because the oxygen bottom flux includes Maybe the best solution is to introduce a |
|
Introducing a I'm struck by how complicated it has become to implement an idea that initially seemed so simple. I guess I should know by now that things are rarely as simple as they initially seem. |
|
Creating And the same variables when we compute I think the issue is that previously we scaled denitrif(k) = work * ((DOC_remin(k) + DOCr_remin(k) + POC_remin(k) * (c1 - POCremin_refract) &
- other_remin(k)) / denitrif_C_N - sed_denitrif(k))
! scale down denitrif if computed rate would consume all NO3 in 10 days
if (NO3_loc(k) < ((c10*spd)*(denitrif(k)+sed_denitrif(k)))) then
work = NO3_loc(k) / ((c10*spd)*(denitrif(k)+sed_denitrif(k)))
denitrif(k) = denitrif(k) * work
sed_denitrif(k) = sed_denitrif(k) * work
end ifIs it reasonable to chalk the bigger-than-roundoff differences up to changes in |
For ciso, set bottom_fluxes in marbl_ciso_interior_tendency_compute() before setting the tendencies. For non-ciso, introduce compute_bottom_fluxes() which gets called before compute_local_tendencies() (ciso doesn't have quite the modularity of the non-iso module, but I think this split fits the existing code in both places)
|
Okay, these are the differences between 3c145c7 (introduction of These are not the only |
The right-hand side of conservation integrals is now to_floor instead of sed_loss (this branch had changed it to sed_loss + (to_floor - sed_loss) and this commit just does the obvious algebraic cancellation). Some additional clean-up renamed variables involved in the vertical integral computation (giving work and work2 more descriptive names) and also cleaning up comments in interior_tendency_mod.
The call_compute_subroutine test did not deallocate the global memory allocated for the bottom fluxes
|
Regarding the hypothesis that the remaining differences are due to the downscaling of |
|
I suggest adding a column where the downscaling of |
|
I'm not immediately seeing a robust way to modify the code in This block of code was added to avoid total depletion of N at the sea floor, which can happen without the downscaling. I'm not sure how to include this safety rail with the loss term due to |
This is significantly easier than the tests I was thinking of, I've started running it now
I definitely think we need a column like that in the test. I'd also like to include the carbon isotopes in testing -- maybe it makes the most sense to keep the current test as is and also introduce a second initial condition file that contains a few different columns as well as CISO initial conditions / forcing fields? As helpful as that would be for this PR, I'd rather finish this work and then add the new test columns... I think there was some testing to make sure I got the data out of POP correctly that would side-track me from the testing I'm doing right now.
One possibility would be to scale the bottom fluxes of NO3 and NH4 by the same factor as the tendencies (e.g. if the scaling term is returned by |
Agreed |
There is a POC_remin term in denitrif that needs to be reflected in bottom_fluxes(no3_ind) that I was missing. It is weighted by min(max(((parm_o2_min + parm_o2_min_delta) - O2_loc) / parm_o2_min_delta, c0), c1) So I assume it wasn't caught in the stand-alone test due to O2_loc(kmt) being laeger than (parm_o2_min + parm_o2_min_delta) in every column?
|
The good news is that testing with the threshold check commented out helped me find a bug (well, the good part is that 55a06ab fixes said bug). With the threshold check uncommented, though, there are still some bigger-than-roundoff errors: The |
This term is then passed to compute_bottom_fluxes and used to scale the POC terms in bottom_fluxes(o2_ind) and bottom_fluxes(no3_ind) to account for the change to denitrif if NO3 concentrations in the bottom level are below a certain threshold.
|
In 7678ac2, I scaled the but it let me do one additional test with the following diff in 410c410
< dissolved_organic_matter%DOCr_remin(:), &
---
> dissolved_organic_matter%DOCr_remin(:), c1 / domain%delta_z(kmt), POC, &
3356c3356
< POC_remin, other_remin, sed_denitrif, denitrif, denitrif_coeff_kmt)
---
> dzr_kmt, POC, POC_remin, other_remin, sed_denitrif, denitrif, denitrif_coeff_kmt)
3367a3368,3369
> real(r8), intent(in) :: dzr_kmt
> type(column_sinking_particle_type), intent(in) :: POC
3378c3380
< real(r8) :: work
---
> real(r8) :: work, work2
3391a3394,3398
> work2 = (c10*spd)*(denitrif(k)+sed_denitrif(k))
> if ((k == kmt) .and. (POC%to_floor > c0)) then
> work2 = work2 + (c10*spd) * work * (POC%to_floor - POC%sed_loss(k)) * (c1 - POCremin_refract) * dzr_kmt / &
> denitrif_C_N
> end if
3393,3394c3400,3401
< if (NO3_loc(k) < ((c10*spd)*(denitrif(k)+sed_denitrif(k)))) then
< work = NO3_loc(k) / ((c10*spd)*(denitrif(k)+sed_denitrif(k)))
---
> if (NO3_loc(k) < work2) then
> work = NO3_loc(k) / work2For context, the first block is adding arguments ( Basically, I added the So I am convinced that everything removed from |
For the call_compute_subroutines test, I had been adding bottom_flux / dz to the tendency to ensure the tendencies were within round-off of the previous values. That was just for testing purposes, and now that this PR is done with testing it makes sense to generate fresh baselines where the bottom fluxes and tendencies are output separately.
Previous commit included the wrong baseline file (it still had BTF added to J_)
There were two lines that I added while testing out some mods, and then commented out because I wasn't sure if I needed them... I don't need them, so rather than leave as comments I am removing them.
Small cleanup to give the variable a more accurate / descriptive name
|
Superseded by #381, which is a much better approach to the problem |
Create a new array
bottom_fluxesin themarbl_interface_classthat represents per-tracer flux into the bottom of the column.Opening this as a draft PR so I can see CI test results (and better share code snippets), but it's not ready for review / merging yet.
Fixes #360