Fixes for sitemptop, sitempbot, and sitempsnic.

[dabail10]

For detailed information about submitting Pull Requests (PRs) to the CICE-Consortium,
please refer to: https://github.com/CICE-Consortium/About-Us/wiki/Resource-Index#information-for-developers

PR checklist

• Short (1 sentence) summary of your PR:
  This is the first in a series of fixes for SIMIP variables. For now, it is just the sitemp variables.
• Developer(s):
  dabail10 (D. Bailey)
• 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.
  https://github.com/CICE-Consortium/Test-Results/wiki/cice_by_hash_forks#ae5b8fd81f597724f1d804902cbb7df56f83f725
• How much do the PR code changes differ from the unmodified code?
  • bit for bit except gnu and cray due to changes in Icepack
  • 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 document the changes in detail, including why the changes are made. This will become part of the PR commit log.

This PR will also require an icepack commit later on. For now, this is the changes for CICE only. Start of addressing issue #1038 and #1033

Changes in CICE:

• move conversion to K into the define_hist_field call
• fix the accumulation of sitemptop, sitempsnic, and sitempbot to reflect the "intrinsic" versus "extrinsic" definitions.
• redo the avg_ice_present division by aice, so that conb is added at the end
• also undo the ravgct multiplication and division for the case of avg_ice_present
• add cell_methods as "area: time: mean where sea ice (mask=siconc)"
• Update icepack hash to bug fix version of "Tsnice".

[NickSzapiro-NOAA]

This isn't bit-for-bit in UFS for sitempsnic and sitempbot. Not that it should be...I just see that box checked.

More substantially, I don't know that we should separate the CICE and Icepack changes for these variables. But maybe I just don't see your plan

[dabail10]

This isn't bit-for-bit in UFS for sitempsnic and sitempbot. Not that it should be...I just see that box checked.

More substantially, I don't know that we should separate the CICE and Icepack changes for these variables. But maybe I just don't see your plan

I guess I think of bfb meaning only the internal variables and not the history variables. In terms of CICE and Icepack, this is the workflow where we have to do each separately and then they come in together.

[anton-seaice]

This looks ok - but CICE-Consortium/Icepack#542 has changed Tsnice (in Celsius) to already be aice weighted, so in combination with this change , it would get weighted by aice twice (sortof)

[dabail10]

This looks ok - but CICE-Consortium/Icepack#542 has changed Tsnice (in Celsius) to already be aice weighted, so in combination with this change , it would get weighted by aice twice (sortof)

I would be happy to have a discussion about this over Google Meet. Icepack is producing the gridcell Tsnice. So, it is:

Tsnice = sum(aicen(n)*Tsnice).

Technically this can be n = 0,ncat where n=0 is the open water category and Tsnice is zero here. So, this is like dividing by sum(aicen(n)), n = 0,ncat = 1. Then in ice_history.F90 we do sum(aice(t)*Tsnice(t)) / sum(aice(t)) over each timestep in the averaging interval. The sum(aice(t)) division is done later on. So this code:

                    if (tmask(i,j,iblk)) then
                          a2D(i,j,n,iblk) = &
                          a2D(i,j,n,iblk)*avgct(ns)*ravgip(i,j)
                    endif

[eclare108213]

The category and time averaging are confusing enough for everyone (including me) that I think it's worth writing it out mathematically and making that available in the docs. I'll work on that.

My question for @dabail10 is whether the CMIP data request asks for Tsnice over the grid cell or over just the ice. I think it's the latter, and the confusing part is why you're only computing a grid-cell value in Icepack.

I think of the time averaging this way:

Average volume V_avg = sum(V) / M, where the sum is over all time steps in the averaging period and M is the number of time steps with ice present.
Average content, e.g. salt_avg = sum(V*S) / M where S is salinity i.e. salt tracer
Average tracer, e.g. salinity_avg = sum(V*S) / sum(V) = salt_avg / V_avg

The quantities being averaged in time here have already been integrated or averaged over categories, and it's always better to use 'content' quantities for conservation and/or accuracy. I think that's what @dabail10 is doing but I want to write it out with units to check. The temperatures are surface quantities so the integrals would be over the surfaces rather than the volumes.

Let me know if you see errors in my thinking here!

[dabail10]

The category and time averaging are confusing enough for everyone (including me) that I think it's worth writing it out mathematically and making that available in the docs. I'll work on that.

My question for @dabail10 is whether the CMIP data request asks for Tsnice over the grid cell or over just the ice. I think it's the latter, and the confusing part is why you're only computing a grid-cell value in Icepack.

I think of the time averaging this way:

Average volume V_avg = sum(V) / M, where the sum is over all time steps in the averaging period and M is the number of time steps with ice present. Average content, e.g. salt_avg = sum(V*S) / M where S is salinity i.e. salt tracer Average tracer, e.g. salinity_avg = sum(V*S) / sum(V) = salt_avg / V_avg

The quantities being averaged in time here have already been integrated or averaged over categories, and it's always better to use 'content' quantities for conservation and/or accuracy. I think that's what @dabail10 is doing but I want to write it out with units to check. The temperatures are surface quantities so the integrals would be over the surfaces rather than the volumes.

Let me know if you see errors in my thinking here!

The SIMIP request only asks for the temperature within / over / under the ice. However, you point is valid in terms of the definition of intensive versus extensive. I think the point is we want to avoid averaging in zeroes that are not physical. So, the aggregate Tsnice is:

Tsnice = sum(aicen(n)*Tisnicen(n)), for n=1,5

From the Notz et al. 2016 paper:

"For all variables that are proportional to area fraction, i.e. extensive variables such as volume, mass, or area fraction, a zero should always be averaged in for all time steps where no sea ice is present. This is because the extensive variables naturally approach zero as area fraction approaches zero."

So in this sense, this definition of Tsnice is indeed extensive. It becomes intensive if we do:

Tsnice_cell = sum(aicen(n)*Tsnicen(n)) / aice(n), n=1,5

this is a weighted mean of Tsnicen, again only over the ice. However, it does not go to zero as aice goes to zero. So, I agree that I should just do a straight time average of Tsnice as it is defined. This is similar for trcr(:,:,nt_Tsfc) and Tbot. For Tbot though since it is all the same for each category, then Tbot_cell = Tbot.

[dabail10]

Here are some images with the updated history averaging and fix to Tsnice in Icepack.

[dabail10]

@eclare108213 @anton-seaice @NickSzapiro-NOAA Hopefully I have addressed your concerns on this as well.

[eclare108213]

One question about average_ice_present...

[eclare108213]

Can you explain why avg_ice_present should be false? Since these temperature variables are supposed to be only over the ice area, it seems to me that avg_ice_present should be true. If it is false, then the accumulated field (e.g. a2D) is not multiplied by ravgip, so the field written out is the average over all space and time for that accumulation interval, right?

[dabail10]

I keep going around and around on this one. The variable sitemptop is trcr(:,:,nt_Tsfc), so the category aggregate surface temperature. As aicen goes to zero, then trcr(:,:,nt_Tsfc) goes to zero and hence is "extensive". However, this is only true in Celsius and 0C is actually a valid temperature. I would like others interpretation here.

[eclare108213]

The "extensive" and "intensive" nomenclature is nonintuitive in my opinion. Since nearly every quantity in the ice model is aggregated over categories, ice area is always going to be a factor and those will all go to zero as area goes to zero. The exceptions will be tracers (or similar) for which you want the "actual" value over the ice. I would think that skin temperature would be one of those, and so the aggregated value ought to be divided by ice area. Can you work through an example case that is easy to understand, by hand, following the code? E.g. let's say the ice area changes from 0 to 0.5 to 1 over 3 time steps, and the temperature changes from nothing to something... What should the time averaged value be over the average ice, and does the code get it right? If not, is the problem with avg_ice_present or somewhere else?

[anton-seaice]

However, this is only true in Celsius and 0C is actually a valid temperature.

I think this is kind of coincidental. (We report temperature in K after all)

It clear in Notz 2016 / the CMIP7 data request that temperatures should be masked / averaged only where there is sea ice

[anton-seaice]

Does removing the avg_ice_present and removing if (aice(i,j,iblk) > puny) later on mean that this field is not masked over land anymore ?

[anton-seaice]

Nevermind, tmask is applied seperately

[anton-seaice]

Sorry @dabail10 - I am still confused, well probably need to have a meeting !

In my very simply example, with one ice thickness category:

If ncat = 1 and aicen = aice = 0.5
And there is no snow
And surface temperature is 2°C

Then icepack produces a grid cell mean temperature (Tsnice) of 1°C (as sum(aicen*Tsf))

If we assume that there is no change during a day, then the implementation in this PR calculates the sitemptop daily average as 274.15K (for either setting of avg_ice_present)

But the correct answer is 275.15K ? (As its a diagnostic of surface temperature over sea ice)

[anton-seaice]

I realised in our discussion this morning we forgot to talk about the right step the C->K conversion should happen (:

[eclare108213]

the right step the C->K conversion should happen

I think that should be last, to just do the calculation once rather than over and over during the time accumulation.

[dabail10]

the right step the C->K conversion should happen

I think that should be last, to just do the calculation once rather than over and over during the time accumulation.

I always forget that we can add a conversion factor in the variable definition. This is added at the end.

[anton-seaice]

Is it correct to use trcr(:,:,nt_Tsfc,:) here?

I believe trcr(:,:,nt_Tsfc,:) is a grid cell average, where the open water area is Tfrz.

But what we want is a ice area aggregate only of Tsfc?

e.g.

sum(aicen(n)*trcrn(nt_Tsfc,n)) for n=1,ncat

without including the open water temperature

[dabail10]

This is a good question. My assumption was that all trcr tracers were just sum(aicen*trcrn) from n=1,5. Is this not the case?

[anton-seaice]

I found this quite hard to figure out, but I think Tf is averaged in.

There is this note for the Tsfc diagnostic:

CICE/cicecore/cicedyn/analysis/ice_history.F90

Lines 757 to 760 in bd30d63

	call define_hist_field(n_Tsfc,"Tsfc","C",tstr2D, tcstr, &
	"snow/ice surface temperature", &
	"averaged with Tf if no ice is present", c1, c0, &
	ns1, f_Tsfc)

And this special handling for Tsfc:

            if (aicen > puny) then
               trcrn(it) = atrcrn(it) / aicen
            else
               trcrn(it) = c0
               if (it == nt_Tsfc) then
                  trcrn(it) = Tf  ! surface temperature
               endif

from https://github.com/CICE-Consortium/Icepack/blob/7ee2266089f80effb5fee8335c4a1056d7da11e7/columnphysics/icepack_tracers.F90#L1288

[dabail10]

Right, so this means trcrn(:,:,nt_Tsfc) is only set to Tf when aicen(n) < puny. So, this would not be averaged in. Here is the code in icepack_aggregate. I suppose this might still average values of Tf when aicen(n) < puny. @eclare108213 is this correct?

      do n = 1, ncat

            aice = aice + aicen(n)
            vice = vice + vicen(n)
            vsno = vsno + vsnon(n)

         do it = 1, ntrcr
            atrcrn = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
                                + trcr_base(it,2) * vicen(n) &
                                + trcr_base(it,3) * vsnon(n))
            if (n_trcr_strata(it) > 0) then  ! additional tracer layers
               do itl = 1, n_trcr_strata(it)
                  ntr = nt_strata(it,itl)
                  atrcrn = atrcrn * trcrn(ntr,n)
               enddo
            endif
            atrcr(it) = atrcr(it) + atrcrn
         enddo                  ! ntrcr
      enddo                     ! ncat

[anton-seaice]

Oh ok - not as problematic as I was thinking. Probably ok then !

[anton-seaice]

the right step the C->K conversion should happen

I think that should be last, to just do the calculation once rather than over and over during the time accumulation.

I always forget that we can add a conversion factor in the variable definition. This is added at the end.

I just realised we might need to change the order of operations to do this, the masking for "avg_ice_present" being true occurs after "conb" (the constant to add to a diagnostic) is added

see

CICE/cicecore/cicedyn/analysis/ice_history.F90

Line 3720 in bd30d63

if (avail_hist_fields(n)%avg_ice_present) then

[dabail10]

I remember this code originally. I think I need to change the logic. The following two lines are accumulating normally:

                    a2D(i,j,n,iblk) = avail_hist_fields(n)%cona*a2D(i,j,n,iblk) &
                                   * ravgct + avail_hist_fields(n)%conb

where the multiplication of ravgct is a division by the number of steps in the averaging period. However, the next set of lines in the avg_ice_present if block are:

                          a2D(i,j,n,iblk) = &
                          a2D(i,j,n,iblk)*avgct(ns)*ravgip(i,j)

Are undoing the multiplication by ravgct by mutiplying by avgct(ns) and then dividing by the accumulated ice area (ravgip). I think this is a bug. I should probably only do the first part when avg_ice_present is .false. and then change the second line to:

                      a2D(i,j,n,iblk) = avail_hist_fields(n)%cona*a2D(i,j,n,iblk) &
                                   * ravgip(i,j) + avail_hist_fields(n)%conb

this was how the pond variables were coded. Again, I think I did this. Not sure why I did it this way.

[dabail10]

Ok. I am happy with the averaging when ice present for sitemptop, sitempsnic, and sitempbot. I also added the updated cell_methods for variables when avg_ice_present = .true. I will work on the aice, aice_mid, aice_init stuff in the next PR.

[apcraig]

I think everything is committed now.

The last push seems to be from 3 days ago if I'm looking push history above correctly, you might check whether you need to push to the repo.

[dabail10]

I think everything is committed now.

The last push seems to be from 3 days ago if I'm looking push history above correctly, you might check whether you need to push to the repo.

Yes, this was it.

58de585

[apcraig]

I ran

--suite quick_suite,base_suite,io_suite,unittest_suite -m derecho -e intel,gnu

Looks like gnu is changing answers. For example,

PASS derecho_gnu_smoke_gx3_8x2_diag1_run5day build
PASS derecho_gnu_smoke_gx3_8x2_diag1_run5day run 10.33 1.43 1.66
PASS derecho_gnu_smoke_gx3_8x2_diag1_run5day test 
PASS derecho_gnu_smoke_gx3_8x2_diag1_run5day generate cice.dh02
FAIL derecho_gnu_smoke_gx3_8x2_diag1_run5day complog cice.bd30d63f31.251101 different-data
FAIL derecho_gnu_smoke_gx3_8x2_diag1_run5day compare cice.bd30d63f31.251101 10.22 1.43 1.70 different-data

Can that be explained? I am going to continue testing, see what else I come up with.

[apcraig]

The cray compiler also shows answers changes.

The nvhpc compiler fails to run cases,

derecho_nvhpc_smoke_gx3_8x2_diag1_run5day: ./cice: error while loading shared libraries: libhugetlbfs.so: cannot open shared object file: No such file or directory

I wonder if that's related to the "upgrade" this week. Weekly testing will give us more information.

[apcraig]

The weekend testing of Icepack shows that both the gnu and cray compiler change answers with the related Icepack change, CICE-Consortium/Icepack#542. See https://github.com/CICE-Consortium/Test-Results/wiki/icepack_by_hash#4954a6f9033f78e5c32bf33780384cbf2d0843e6. Both debug and non-debug tests change answers which suggests that even at -O0, the results change. On the other hand, the intel compiler is bit-for-bit. Anyway, that suggests the related Icepack changes are likely leading to the non-bit-for-bit results in this CICE PR. That Icepack PR indicated bit-for-bit, but testing seems to have been a little inadequate. This PR is showing the same behavior, Intel bit-for-bit, but cray and gnu not bit-for-bit with debug and non-debug tests.

So if we accept that the Icepack results change answers for gnu and cray compilers with CICE-Consortium/Icepack#542, then these results are probably OK.

[apcraig]

The nvhpc compiler issues are unrelated to this PR. I tested last weeks Icepack main which worked last week and now it's not running. I will reach out to CISL and update the port in a separate PR if needed. We can ignore the nvhpc issues for now.

[apcraig]

@dabail10, I'm ready to approve this PR. Could you update the PR checklist at the bottom and summarize the changes so I can make sure the merge has a comprehensive commit message.

[dabail10]

@dabail10, I'm ready to approve this PR. Could you update the PR checklist at the bottom and summarize the changes so I can make sure the merge has a comprehensive commit message.

I guess I am not sure what you mean here. Do I need to do some more testing on Icepack to figure out the change? Do I need to run the QC test?

[apcraig]

I think testing is OK at this point. What would be helpful is a fairly comprehensive summary of the current changes written into the bottom of the PR checklist (you can edit that box). Then I'll cut and paste that into the commit when I squash and merge.

[eclare108213]

It's a little concerning to me that the answers change even without any optimization. The Icepack changes are minimal -- I can see why Tsnice itself would not be BFB, but it's only diagnostic in both Icepack and CICE, so shouldn't be changing the BFBness in the test suites. Is this a problem with the code, or is it a problem with the gnu and cray compilers? Since intel is BFB, I'd suspect the latter, but at -O0?

[dabail10]

It's a little concerning to me that the answers change even without any optimization. The Icepack changes are minimal -- I can see why Tsnice itself would not be BFB, but it's only diagnostic in both Icepack and CICE, so shouldn't be changing the BFBness in the test suites. Is this a problem with the code, or is it a problem with the gnu and cray compilers? Since intel is BFB, I'd suspect the latter, but at -O0?

I have rechecked the code and the only changes are related to sitemptop, sitempsnic, and sitempbot which are all diagnostic. There was a derecho update last week that looks to have changed answers.

[dabail10]

I think testing is OK at this point. What would be helpful is a fairly comprehensive summary of the current changes written into the bottom of the PR checklist (you can edit that box). Then I'll cut and paste that into the commit when I squash and merge.

Ok. I have updated the summary.

[apcraig]

The answer changing is a little surprising. I hope it's not from the upgrade, we are specifying the modules. I will run last weeks tag against last weeks test results to check. That would clarify several things though. The answer changing results for debug tests is surprising too. Usually if it's compiler optimization, it doesn't affect -O0 compilations. I'll run some more tests.

[apcraig]

I have confirmed the answer changes we're seeing on Derecho in the last week (since the upgrade) are the result of the upgrade. I am following up with CISL, and I'm retesting this PR against this past weekend's results to confirm it's bit-for-bit.

[apcraig]

OK, I retested against the after update results with

-m derecho -e intel,gnu,cray --suite quick_suite,base_suite,io_suite,unittest_suite --bcmp cice.8e3ef7c4cb.251108

And everything looks fine. Results are bit-for-bit. I continue to discuss with CISL the answer changes on Derecho with gnu and cray with the upgrade last week.

I think this can be merged. I will merge this week unless I hear other comments.

[eclare108213]

@anton-seaice Have your questions and concerns above been adequately addressed? Some of them were postponed to the next PR. I mainly want to make sure that the answers for all three variables now reflect what you'd expect as a CMIP/SIMIP user.

[anton-seaice]

@anton-seaice Have your questions and concerns above been adequately addressed? Some of them were postponed to the next PR. I mainly want to make sure that the answers for all three variables now reflect what you'd expect as a CMIP/SIMIP user.

I think what this change does is ok by the definition:

sitempbot, sitemptop and sitempsnic are defined in Kelvin, and are being calculated as a mean weighted by sea ice area fraction. Because the output is by aice>puny, users can be reasonably confident its a area-weighted mean. Cell_methods are defined correctly for these variables (although they seem to be still in a state of flux in a general sense in the data request, so it's probably hard to get them correct straight out of the model).

I guess users may want to know that sitempbot and sitempsnic are not tracers - so are approximate?

In a general sense, applying consistent notes and cell methods to the output would help but its not really specific to this change. (e.g. the metadata for all variables could state what time averaging & spatial averaging was applied). Similarly - internally labelling variables consistently based on their area fraction would help developers a lot (but also unrelated to this change) (e.g. Tsnice is a grid cell aggregate and Tbot is an true temperature). The other unresolved part is obviously applying aice weightings consistent - Tsnice/Tbot are from the middle of the timestep, but then the area weighting in time uses aice at the end of the timestep.

Anyway - I don't think i've said anything new - happy to merge!

[eclare108213]

This work has highlighted how much more there is to do, to improve the history output as well as more general conservation issues.
@dabail10 thanks for the huge effort to get these issues fixed and for your patience with the rest of us.
@anton-seaice and @nick thanks for your thoughtful reviews.
@apcraig thanks for your extensive testing.
There will likely be more refactoring, but for now I think this PR is ready to merge.