diag bugs with AusCOM driver (affecting ACCESS-OM2)

[aekiss]

Thanks to Stewart Allen for flagging this issue (see Slack discussion https://arccss.slack.com/archives/C6PP0GU9Y/p1627269245007400).

These diagnostics are identical in access-om2, and shouldn't be:

• fresh and fresh_ai
• fsalt and fsalt_ai
• fhocn and fhocn_ai
• fswthru and fswthru_ai

This issue probably also affects these diagnostics and their _ai counterparts: alvdr, alidr, alvdf, alidf, fNO, fNH, fN, fSil, but I haven't checked.

I did a test run in /home/156/aek156/payu/1deg_jra55_iaf_cice_diag_test, which gives these test results for equality of some diagnostics and their _ai counterparts:

import xarray as xr

ds = xr.open_dataset('/scratch/v45/aek156/access-om2/archive/1deg_jra55_iaf_cice_diag_test/output132/ice/OUTPUT/iceh.1968-02.nc')
allvars = list(ds.variables.keys())
vs = [v for v in allvars if v[:-2]+'_ai_m' in allvars]

for v in vs:
    print(v, ds[v].equals(ds[v[:-2]+'_ai_m']))

prints

snow_m False
rain_m False
fswabs_m False
flat_m False
fsens_m False
flwup_m False
evap_m False
fresh_m True
fsalt_m True
fhocn_m True
fswthru_m True

[aekiss]

The CICE5.1 manual p75 says:

Several history variables are available in two forms, a value representing an average over the sea ice fraction of the grid cell, and another that is multiplied by ai, representing an average over the grid cell area. Our naming convention attaches the suffix “_ai” to the grid-cell-mean variable names.

(This is a confusing convention. My interpretation is: if a variable name doesn't end in _ai, we need to check if an _ai version is a namelist option to decide how to interpret it. If there's an _ai version, then the non-ai variable is an average over the ice area. If there's no _ai version, the non-ai variable is a grid-cell mean.)

[aekiss]

In the AusCOM driver the _ai forms of these variables are first initialised to the non-ai values:
https://github.com/COSIMA/cice5/blob/2572851/drivers/auscom/CICE_RunMod.F90#L560-L577

alvdf
alidf
alvdr
alidr
fresh
fsalt
fhocn
fswthru
flux_bio

In some drivers the subroutine scale_fluxes is then called to divide the following non-ai values by aice (the ice area fraction), so the non-ai variables are ice area averages and the _ai variables are grid-cell means, as exected: https://github.com/COSIMA/cice5/blob/2572851/source/ice_flux.F90#L833-L1001
Note that this is a superset of the list of variables that were initialised above.

strairxT
strairyT
fsens
flat
fswabs
flwout
evap
Tref
Qref
Uref
fresh
fsalt
fhocn
fswthru
alvdr
alidr
alvdf
alidf
flux_bio
faero_ocn
fsurf
fcondtop

But (as @russfiedler pointed out) scale_fluxes isn't called in the AusCOM driver used by ACCESS-OM2: https://github.com/COSIMA/cice5/blob/2572851/drivers/auscom/CICE_RunMod.F90#L591-L626, so for the variables in the first list the _ai and non-ai versions remain identical.
(This has apparently been the case since 2014: https://github.com/COSIMA/cice5/blob/9a52804/drivers/auscom/CICE_RunMod.F90#L658-L693.)

Consequently (as confirmed above), these non-ai diagnostics and their _ai counterparts are both grid cell averages (i.e. the non-ai values are incorrect):

• fresh, fsalt, fhocn, fswthru: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1447-L1463

This is probably also the case for these, which output the _ai and non-ai variables without recalculation (though I haven't confirmed this using actual outputs):

• alvdr, alidr, alvdf, alidf: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1381-L1396
• fNO, fNH, fN, fSil (using flux_bio): https://github.com/COSIMA/cice5/blob/2572851/source/ice_history_bgc.F90#L614-L640

[aekiss]

In contrast, for the following diagnostics the _ai version is recalculated by multiplying the non-ai version by aice, so the values differ (see test above).
However, since the non-ai values have not been divided by aice in scale_fluxes, it's possible (I haven't checked) that in these cases the non-ai diagnostics are actually grid cell averages.
If so, both _ai and non-ai variables are incorrect (the non-ai ones are actually _ai, and the _ai ones are just plain wrong).

• fsens, flat, flwup (using flout), evap: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1407-L1422
• fswabs: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1366-L1373

[aekiss]

Only _ai versions are output for these. They might actually be grid averages multiplied by aice, i.e. wrong (see above):

• fsurf, fcondtop: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1516-L1520

[aekiss]

These are output only as non-_ai versions. I guess we'd expect them to be grid cell averages?

• Tref, Qref: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1426-L1429
• faero_ocn: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history_bgc.F90#L548-L553

[aekiss]

The remaining variables that didn't go through scale_fluxes are

• strairxT, strairyT, Uref

I'm not sure whether they will affect any diagnostics.

[hakaseh]

In the AusCOM driver the _ai forms of these variables are first initialised to the non-ai values: https://github.com/COSIMA/cice5/blob/2572851/drivers/auscom/CICE_RunMod.F90#L560-L577

alvdf
alidf
alvdr
alidr
fresh
fsalt
fhocn
fswthru
flux_bio

In some drivers the subroutine scale_fluxes is then called to divide the following non-ai values by aice (the ice area fraction), so the non-ai variables are ice area averages and the _ai variables are grid-cell means, as exected: https://github.com/COSIMA/cice5/blob/2572851/source/ice_flux.F90#L833-L1001 Note that this is a superset of the list of variables that were initialised above.

strairxT
strairyT
fsens
flat
fswabs
flwout
evap
Tref
Qref
Uref
fresh
fsalt
fhocn
fswthru
alvdr
alidr
alvdf
alidf
flux_bio
faero_ocn
fsurf
fcondtop

But (as @russfiedler pointed out) scale_fluxes isn't called in the AusCOM driver used by ACCESS-OM2: https://github.com/COSIMA/cice5/blob/2572851/drivers/auscom/CICE_RunMod.F90#L591-L626, so for the variables in the first list the _ai and non-ai versions remain identical. (This has apparently been the case since 2014: https://github.com/COSIMA/cice5/blob/9a52804/drivers/auscom/CICE_RunMod.F90#L658-L693.)

Consequently (as confirmed above), these non-ai diagnostics and their _ai counterparts are both grid cell averages (i.e. the non-ai values are incorrect):

• fresh, fsalt, fhocn, fswthru: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1447-L1463

This is probably also the case for these, which output the _ai and non-ai variables without recalculation (though I haven't confirmed this using actual outputs):

• alvdr, alidr, alvdf, alidf: https://github.com/COSIMA/cice5/blob/2572851/source/ice_history.F90#L1381-L1396
• fNO, fNH, fN, fSil (using flux_bio): https://github.com/COSIMA/cice5/blob/2572851/source/ice_history_bgc.F90#L614-L640

I checked and confirm that BGC fluxes are identical between _ai and non-_ai (for example, fNO vs fNO_ai).

My question is whether they are really grid cell average or not. I ask this because flux_bio is multiplied by ice area (aicen in the below code) when updating the ocean surface BGC concentration:

cice5/source/ice_algae.F90

Lines 467 to 470 in 2572851

	! Uncomment to update ocean concentration
	! Currently not coupled with ocean biogeochemistry
	!ocean_bio(i,j,nn) = ocean_bio(i,j,nn) &
	! + flux_bio(i,j,nn)/hmix(i,j)*aicen(i,j)

This multiplication by ice area suggests that flux_bio is not grid cell average, and therefore needs to be multiplied by ice area to get grid cell average to update ocean surface concentration.

[hakaseh]

Just to follow up on my previous comment.
I contacted Nicole Jeffery (the main developer of sea-ice BGC in CICE), and she indicated that flux_bio_ai is the quantity before divided by aice in scale_fluxes. So in our auscom case, flux_bio and flux_bio_ai are identical and they are grid cell averages.