Some CMIP variables are computed using a mix of U and T quantities

[phil-blain]

Hi @dabail10,

I realized that some of the CMIP output variables are computed by multiplying the concentration aice (thus a T-point quantity) with variables from the dynamics which are U-point quantities. I'm wondering if we should not use grid_average_X2Y to shift the dynamical quantities to the T-point prior to computing these variables ?

Here is a diff that shows which variables I'm talking about:

diff --git a/cicecore/cicedyn/analysis/ice_history.F90 b/cicecore/cicedyn/analysis/ice_history.F90
index 570ef495..6340d503 100644
--- a/cicecore/cicedyn/analysis/ice_history.F90
+++ b/cicecore/cicedyn/analysis/ice_history.F90
@@ -2744,11 +2744,11 @@ subroutine accum_hist (dt)
 
          if (f_siu(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
-              if (aice(i,j,iblk) > puny) worka(i,j) = aice(i,j,iblk)*uvel(i,j,iblk)
+              if (aice(i,j,iblk) > puny) worka(i,j) = aice(i,j,iblk)*uvel(i,j,iblk) ! should be computed at the same location ?
            enddo
            enddo
            call accum_hist_field(n_siu, iblk, worka(:,:), a2D)
          endif
 
@@ -2816,11 +2816,11 @@ subroutine accum_hist (dt)
          if (f_sistrxdtop(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice_init(i,j,iblk) > puny) &
-                 worka(i,j) = aice(i,j,iblk)*(aice(i,j,iblk)*strairxU(i,j,iblk)/aice_init(i,j,iblk))
+                 worka(i,j) = aice(i,j,iblk)*(aice(i,j,iblk)*strairxU(i,j,iblk)/aice_init(i,j,iblk)) !idem
            enddo
            enddo
            call accum_hist_field(n_sistrxdtop, iblk, worka(:,:), a2D)
          endif
 
@@ -2838,11 +2838,11 @@ subroutine accum_hist (dt)
          if (f_sistrxubot(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) &
-                 worka(i,j) = aice(i,j,iblk)*strocnxU(i,j,iblk)
+                 worka(i,j) = aice(i,j,iblk)*strocnxU(i,j,iblk) ! idem
            enddo
            enddo
            call accum_hist_field(n_sistrxubot, iblk, worka(:,:), a2D)
          endif
 
@@ -3302,11 +3302,11 @@ subroutine accum_hist (dt)
          if (f_sidragtop(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = aice(i,j,iblk)*Cdn_atm(i,j,iblk)
+                 worka(i,j) = aice(i,j,iblk)*Cdn_atm(i,j,iblk) !idem
               endif
            enddo
            enddo
            call accum_hist_field(n_sidragtop, iblk, worka(:,:), a2D)
          endif
@@ -3327,11 +3327,11 @@ subroutine accum_hist (dt)
          if (f_siforcetiltx(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = aice(i,j,iblk)*strtltxU(i,j,iblk)
+                 worka(i,j) = aice(i,j,iblk)*strtltxU(i,j,iblk) !idem
               endif
            enddo
            enddo
            call accum_hist_field(n_siforcetiltx, iblk, worka(:,:), a2D)
          endif
@@ -3339,11 +3339,11 @@ subroutine accum_hist (dt)
          if (f_siforcetilty(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = aice(i,j,iblk)*strtltyU(i,j,iblk)
+                 worka(i,j) = aice(i,j,iblk)*strtltyU(i,j,iblk) ! idem
               endif
            enddo
            enddo
            call accum_hist_field(n_siforcetilty, iblk, worka(:,:), a2D)
          endif
@@ -3351,11 +3351,11 @@ subroutine accum_hist (dt)
          if (f_siforcecoriolx(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = aice(i,j,iblk)*fmU(i,j,iblk)*vvel(i,j,iblk)
+                 worka(i,j) = aice(i,j,iblk)*fmU(i,j,iblk)*vvel(i,j,iblk) ! idem
               endif
            enddo
            enddo
            call accum_hist_field(n_siforcecoriolx, iblk, worka(:,:), a2D)
          endif
@@ -3363,11 +3363,11 @@ subroutine accum_hist (dt)
          if (f_siforcecorioly(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = -aice(i,j,iblk)*fmU(i,j,iblk)*uvel(i,j,iblk)
+                 worka(i,j) = -aice(i,j,iblk)*fmU(i,j,iblk)*uvel(i,j,iblk) ! idem
               endif
            enddo
            enddo
            call accum_hist_field(n_siforcecorioly, iblk, worka(:,:), a2D)
          endif
@@ -3375,11 +3375,11 @@ subroutine accum_hist (dt)
          if (f_siforceintstrx(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = aice(i,j,iblk)*strintxU(i,j,iblk)
+                 worka(i,j) = aice(i,j,iblk)*strintxU(i,j,iblk) ! idem
               endif
            enddo
            enddo
            call accum_hist_field(n_siforceintstrx, iblk, worka(:,:), a2D)
          endif

I'm working on some slight refactoring of the history code so I could include such a change if you think it makes sense.

[dabail10]

This was a debate. According to the SIMIP protocol, everything was supposed to stay on the native grid. I agree that I probably should have interpolated these. I ended up leaving this as is. I am open to changing this.

[dabail10]

Note that these are sum(aice(t)*strintxU(t)) / sum(aice(t)) where t is over the averaging interval. So, this works out fine even though they are not exactly co-located. I could average aice to the B-grid points in both the numerator and denominator, but this is extra work that is really not necessary in my opinion.

[anton-seaice]

Could this issue be re-opened, I think it needs a bit more discussion.

The only note I found was this in Notz et al 2016, which is specific to velocity:

The variables that SIMIP requests for sea-ice dynamics are needed to understand how and why sea ice moves horizontally. The most important parameters there are obviously the actual velocities, which are requested on the native model grid. In addition, the sea-ice speed is requested, which allows one to account for possible back-and-forth movement of the ice during the averaging period. Because of their importance, these variables are requested at priority 1 and at daily resolution.

On this

sum(aice(t)*strintxU(t)) / sum(aice(t)) where t is over the averaging interval. So, this works out fine even though they are not exactly co-located.

This doesn't exactly work out equal. In most cases it would be similar, but in cases of strong gradients in aice it could be different. (quite different ?)

Also, if i understand correctly, these diagnostics should be different between B and C grid ?

e.g. for B-grid, should be:

if (f_siforcetilty(1:1) /= 'x') then
...
                 worka(i,j) = aiceU(i,j,iblk)*strtltyU(i,j,iblk) 

but C-grid

worka(i,j) = aiceN(i,j,iblk)*strtltyN(i,j,iblk) 

There is interpolation done between grids in other places in CICE - can it be reused for aice for these diagnostics ?

p.s. Low priority for me - it's possibly too late for CMIP7 now.

[eclare108213]

Good questions, happy to reopen this topic. Are the lines highlighted in the original issue post the only variables affected, or are there others?

[anton-seaice]

There is also this guidance document on grids, which includes this suggestion:

For model output, it is recommended that each grid-cell value be representative of the entire grid cell and carried at its cell center, but it is not forbidden for a variable to represent a value at a cell bound or corner. Analysts usually prefer all variables on a common grid, but there may be some circumstances when reporting variables on different grids could be judged by a modelling group to be of more importance than analysts' preferences.

Implying maybe interpolating the vectors to cell centre should be done (rather than) interpolating sea ice concentration to corner points.

Are the lines highlighted in the original issue post the only variables affected, or are there others?

I am not sure I know which of the native variables are on cell corners / edges well enough to answer this !

It might be all CMIP variables which depend on these native vars:

CICE/cicecore/cicedyn/analysis/ice_history.F90

Lines 340 to 372 in ce4dd0b

	if (grid_ice == 'CD' .or. grid_ice == 'C') then
	f_uvelE = f_uvel
	f_vvelE = f_vvel
	f_icespdE = f_icespd
	f_icedirE = f_icedir
	f_uvelN = f_uvel
	f_vvelN = f_vvel
	f_icespdN = f_icespd
	f_icedirN = f_icedir
	f_strairxN = f_strairx
	f_strairyN = f_strairy
	f_strairxE = f_strairx
	f_strairyE = f_strairy
	f_strocnxN = f_strocnx
	f_strocnyN = f_strocny
	f_strocnxE = f_strocnx
	f_strocnyE = f_strocny
	f_strcorxN = f_strcorx
	f_strcoryN = f_strcory
	f_strcorxE = f_strcorx
	f_strcoryE = f_strcory
	f_strintxN = f_strintx
	f_strintyN = f_strinty
	f_strintxE = f_strintx
	f_strintyE = f_strinty
	f_strtltxN = f_strtltx
	f_strtltyN = f_strtlty
	f_strtltxE = f_strtltx
	f_strtltyE = f_strtlty
	f_taubxN = f_taubx
	f_taubyN = f_tauby
	f_taubxE = f_taubx
	f_taubyE = f_tauby