Update grid averaging for tmass, aice, uvelT, vvelT

[apcraig]

PR checklist

• Short (1 sentence) summary of your PR:
  Update grid averaging for tmass, aice, uvelT, vvelT
• Developer(s):
  apcraig
• Suggest PR reviewers from list in the column to the right.
• Please copy the PR test results link or provide a summary of testing completed below.
  Everything runs, changes most results, https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#658799554aafd24f9c2fd575b68b47b281b03027. QC tests pass relative to current main (cice.tr502 is current main, #6587995

PASS cheyenne_intel_qcchk_gx1_144x1_medium_qc_qcchk compare cice.tr502 -1 -1 -1
PASS cheyenne_pgi_qcchk_gx1_144x1_medium_qc_qcchk compare cice.tr502 -1 -1 -1
PASS cheyenne_gnu_qcchk_gx1_144x1_medium_qc_qcchk compare cice.tr502 -1 -1 -1

• How much do the PR code changes differ from the unmodified code?

  • bit for bit
  • different at roundoff level
  • more substantial

• Does this PR create or have dependencies on Icepack or any other models?

  • Yes
  • No

• Does this PR update the Icepack submodule? If so, the Icepack submodule must point to a hash on Icepack's main branch.

  • Yes
  • No

• Does this PR add any new test cases?

  • Yes
  • No

• Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/. A test build of the technical docs will be performed as part of the PR testing.)

  • Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • No

• Please provide any additional information or relevant details below:

• Update tmass and aice T2U mapping, switch from "F" to "S", F was backwards compatible but not correct (changes answers), closes change to_ugrid calculation of mass, aice to use T2US instead of T2UF (check and fix everywhere, eventually) #736

• Update ocean forcing T2U averaging in ocn_data_ncar_init, change "F" to "A".

• Update uvelT, vvelT averaging in step_therm1, change from 4 point average to U2TA (changes answers for highfreq=.true.), closes change ice_step_therm1 to compute uvel, vvel at cell centers using new averaging routine #735

• Remove history grids not needed (i.e. ustr3Dz), see C-grid Code Cleanup  #660

[eclare108213]

The change here for tmass is right: to get the correct mass at the U point, we need to compute
sum(rhoi*vice*aice*T cell area)/sum(aice*T cell area)
over the T cells [ignoring snow for now], since vice has units of volume per unit area.

But I do not think this is correct for aice_init. Since ice concentration is assumed to be uniformly distributed over each grid cell, the concentration at U should be
sum(aice_init*T cell area)/sum(T cell area)
but I think the calculation here is
sum(aice_init*aice*T cell area)/sum(aice*T cell area).

[apcraig]

Just to clarify, The "S" mapping is masked area weighted while the "A" mapping is just area weighted. When you use "aice" above, is that the ice fraction? There are no ice fraction weighted averages built into CICE. If we want to do that, we can use user defined averaging weights. Again,

"S":  sum(var*Tmask*Tarea)/sum(Tmask*Tarea)
"A":  sum(var*Tarea)/sum(Tarea)
"F":  sum(var*Tarea_var)/Tarea_dst

Where Tarea and Tmask are the static grid cell areas and masks computed at initialization. If we do want to have fraction weighted mapping, we can do that, but it's not the standard "S", "A", or "F". To do that we'd do

mask_s(:,:,;) = 1
wght_s(:,:,:) = Tarea(:,:,:) * aice(:,:,:)
call grid_average_X2Y('A', tmass, 'T', wght_s, mask_s, umass, 'U')

and that would do

umass = sum(tmass*Tarea*aice)/sum(Tarea*aice)

At this point, we NEVER weigh the grid averaging by the ice fraction anywhere. Is that something we should be doing?

[eclare108213]

Sorry, I misunderstood "area weighted" in the doc to mean ice-area-weighted, but it clearly says grid cell area. My mistake. Let me think about this some more - what I wrote above is not quite correct.

[eclare108213]

Correction to my comment above:

To get the correct mass at the U point, we need to compute
sum(rhoi*vice*Tarea)/sum(Tarea)
over the T cells [ignoring snow for now], since vice has units of volume per unit area and already incorporates aice (i.e. vice = hi * aice). So the tmass line is correct.

Since ice concentration is assumed to be uniformly distributed over each grid cell, the initial concentration at U should be
sum(aice_init*Tarea)/sum(Tarea)
and that is what's being done in the S mapping.

So I agree with both of these modifications. Sorry for the confusion!

[apcraig]

The "S" mapping includes the mask, the "A" mapping does not. I think we need to change "S" to "A" for umass and aiU if we want an unmasked weighted average as written above. This only matters for gridcells near land. So should the average be an average of only ocean cells ("S" mapping) or an average of all ocean/land cells ("A" mapping)?

One other thing that I just want reiterate. Maybe there is a place for fraction weighted mapping in CICE. We do that in the coupler when we map ice coupling fields to other grids. We fraction area weight the mapping so a gridcell with fraction=0.9 has much greater weight than a gridcell with fraction=0.2 when interpolating from ice to atmosphere for example. You can also think of it as a fraction weighted merge. How we do that depends whether it's a flux (which has to be strictly conserved and has a particular mapping requirement) or a state (which could be done many ways). Maybe we need to do the same in CICE in some places for some fields?? I guess it partly depends whether we can define an ice fraction on the grid we're mapping from. We'd have to think about it a bit.

[eclare108213]

I think S is correct in this case, with the mask included, although it would be good for @JFLemieux73 @dabail10 to confirm that it's consistent with other assumptions in the C-grid approach. On a B-grid, the velocity being computed would be zero if any of the T cells is land (i.e. the point is moot).

[apcraig]

Great, thanks for the clarification.

[eclare108213]

I also wondered whether ice area weighting could/should be included as part of the averaging infrastructure. For coupling, the actual averaging happens in the coupler with other components' data, so we can't control the entire average anyhow. Once the model starts being coupled with different grids, we should check whether we are calculating the coupling fields appropriately, but I think you've already thought through all of that. In the ice model, we "aggregate" across ice thickness categories, but that's only within a particular cell, not for shifting things around horizontally (there are vertical grid remappings). The tracer loading/unloading also involves a lot of multiplying/dividing by "weights" (including ice area and other tracers). Neither of these cases involve grid shifts. So: I don't think we should worry about it unless we come across a particular case in which we need it.