Ensure non-negative tracer values and conservation by virtual flux correction

[matsbn]

Tracers are enforced to be non-negative and and any tracer quantity added to ensure non-negative values are compensation for by a global virtual flux correction.

The simulation n1850.ne30_tn14.hybrid_wavparam_vfcor_fatessp.20250204 was started as a hybrid simulation from year 29 of NorESMhub/noresm3_dev_simulations#68 to test the impact of this PR. Diagnostics can be found here: https://ns2345k.web.sigma2.no/datalake/diagnostics/noresm/matsbn/n1850.ne30_tn14.hybrid_wavparam_vfcor_fatessp.20250204/

Some erratic biogeochemistry global flux time-series (net downward CO2 flux) and global averages (DIC and alkalinity) seems to be mitigated in the new simulation. It should be noted that recent changes to iHAMOCC included in the code base used for 1850.ne30_tn14.hybrid_wavparam_vfcor_fatessp.20250204 could also contribute to this improvement. Then there are difference to the evolution of global averages of oxygen, phosphate and nitrate that would be useful with comments on.

Global average salinity has a stronger increase in n1850.ne30_tn14.hybrid_wavparam_vfcor_fatessp.20250204 compared to NorESMhub/noresm3_dev_simulations#68. Apparently and unintentionally, CICE started from the initial condition set in user_nl_cice, and not the year 29 restart of NorESMhub/noresm3_dev_simulations#68. The resulting sea-ice volume build-up might explain the salinity evolution. A new simulation is underway to make sure CICE starts from the year 29 restart of NorESMhub/noresm3_dev_simulations#68.

[TomasTorsvik]

@matsbn - Thanks, the changes in oxygen, phosphate and nitrate are likely due to PR #464. I tested the impact of this change, in
NorESMhub/noresm3_dev_simulations#77
NorESMhub/noresm3_dev_simulations#74
See in particular
https://ns2345k.web.sigma2.no/datalake/diagnostics/noresm/tomast/n1850.ne30pg3_tn14.2_5_alpha08d.BLOMv1.7.0.hybrid_wavparam_fatessp.20250124/HAMOCC_DIAG/ts1to28-n1850.ne30_tn14.hybrid_fatessp.20241204/set1/set1_ann_avg_2models.png

[JorgSchwinger]

A general question about the time-stepping: Now the salt and tracer correction is calculated but applied one time-step later, correct? This is why it has to be saved in the restart files. When in this sequence is HAMOCC called? Is this after the correction has been applied? HAMOCC cannot be called with negative tracer values. In fact we have the routine trc_limitc, which is called before HAMOCC and makes sure no negative values occur. Assuming that the virtual fluxes are the only potential source of negative values, this routine should then be removed?

[matsbn]

A general question about the time-stepping: Now the salt and tracer correction is calculated but applied one time-step later, correct? This is why it has to be saved in the restart files. When in this sequence is HAMOCC called? Is this after the correction has been applied? HAMOCC cannot be called with negative tracer values. In fact we have the routine trc_limitc, which is called before HAMOCC and makes sure no negative values occur. Assuming that the virtual fluxes are the only potential source of negative values, this routine should then be removed?

Yes, the tracer added to ensure non-negative values is compensated as a correction of the virtual tracer flux applied the next time-level. Ensuring non-negative tracer values (of the 'k+nn' time-level) happens in routine 'mxlayr' when isopycnic coordinate is used and in routine 'ale_vdifft' otherwise. These routines are called before routine 'updtrc' is called, so there should be no negative tracer values entering HAMOCC (assuming 'k+nn' time-level is used). So I think the limiting in 'trc_limitc' can now be skipped.

[JorgSchwinger]

Regarding the trc_limitc routine, I would suggest merging this now and then removing the routine later (we can first issue an error if negative tracer values are encountered to make sure everything works as expected). I can do this in a subsequent PR.

[matsbn]

Regarding the trc_limitc routine, I would suggest merging this now and then removing the routine later (we can first issue an error if negative tracer values are encountered to make sure everything works as expected). I can do this in a subsequent PR.

Sounds like a good plan. A simulation NorESMhub/noresm3_dev_simulations#86 with correct CICE restart file for a hybrid run of NorESMhub/noresm3_dev_simulations#68 has soon completed 7+7 years. After the first 7 years, the salinity evolution is very close to the reference case: https://ns2345k.web.sigma2.no/datalake/diagnostics/noresm/matsbn/n1850.ne30_tn14.hybrid_wavparam_vfcor_fatessp.20250205/BLOM_DIAG/yrs29to35-n1850.ne30_tn14.hybrid_wavparam_fatessp.20241226/set1/set1_ann_salnga_2models.png.

Let me verify that the last 7 years also look reasonable, then I'll merge.

[matsbn]

The simulation n1850.ne30_tn14.hybrid_wavparam_vfcor_fatessp.20250205 (NorESMhub/noresm3_dev_simulations#86) completed 13 years (crash in final year that seems unrelated to the code changes).

Diagnostics can be found here: https://ns2345k.web.sigma2.no/datalake/diagnostics/noresm/matsbn/n1850.ne30_tn14.hybrid_wavparam_vfcor_fatessp.20250205/

Time series of global average salinity and iHAMOCC fields are shown below. Salinity evolution seems fine and also iHAMOCC averaged variables with the additional iHAMOCC changes taken into account.