kinetic energy diagnostic mixes T and U variables for B grid

[eclare108213]

No description provided.

[apcraig]

I see this in ice_diagnostics.F90. I assume what we need to do is compute work1 on either the T or U grid consistently? It seems like it would be better to compute the KE on the U grid, moving the T variables to the U grid? For the C grid, would we compute the "u" part of KE on the E face and the "v" part of KE on the N face then add them? Or is it better to just average the U,V to T and compute KE on T?

      ! total ice-snow kinetic energy                                                                                                                                                                         
      !$OMP PARALLEL DO PRIVATE(iblk,i,j)                                                                                                                                                                     
      do iblk = 1, nblocks
         do j = 1, ny_block
         do i = 1, nx_block
            work1(i,j,iblk) = p5 &
                           * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) &
                           * (uvel(i,j,iblk)**2 + vvel(i,j,iblk)**2)
         enddo
         enddo
      enddo
      ! Eventually do energy diagnostic on T points.                                                                                                                                                          
!     if (grid_ice == 'CD') then                                                                                                                                                                              
!     !$OMP PARALLEL DO PRIVATE(iblk,i,j)                                                                                                                                                                     
!     do iblk = 1, nblocks                                                                                                                                                                                    
!        do j = 1, ny_block                                                                                                                                                                                   
!        do i = 1, nx_block                                                                                                                                                                                   
!           call grid_average_X2Y('E2TS',uvelE,uvelT)                                                                                                                                                         
!           call grid_average_X2Y('N2TS',vvelN,vvelT)                                                                                                                                                         
!           work1(i,j,iblk) = p5 &                                                                                                                                                                            
!                          * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) &                                                                                                                                    
!                          * (uvelT(i,j,iblk)*uvelT(i,j,iblk) &                                                                                                                                               
!                          +  vvelT(i,j,iblk)*vvelT(i,j,iblk))                                                                                                                                                
!        enddo                                                                                                                                                                                                
!        enddo                                                                                                                                                                                                
!     enddo                                                                                                                                                                                                   
!     endif                                                                                                                                                                                                   
!     !$OMP END PARALLEL DO                                                                                                                                                                                   
      ketotn = c0
      ketots = c0
      ketotn = global_sum(work1, distrb_info, field_loc_center, tarean)
      ketots = global_sum(work1, distrb_info, field_loc_center, tareas)

[apcraig]

See also #663. I'm going to close that now as it seems redundant.

[apcraig]

Would this also address tasks in #655? Or are there still more things to do there. And if so, should I include them now?

[eclare108213]

MPAS-Ocean also mixes cell-centered and vertex values in its diagnostic KE calculation. See https://github.com/E3SM-Project/E3SM/blob/master/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F. That code is pretty complicated, but I think the take-home point is that it uses velocities from all of the edges (with appropriate area weighting) around a cell to calculate a cell-centered value. Then it adds additional information from vertices, which are also calculated using the nearest edge velocities. I'm not sure we need the second step, but we could follow the guidance in the first step for getting a cell-centered value, and use that as the primary diagnostic. Would the diagnostic be needed at other cell locations?

[eclare108213]

Would this also address tasks in #655? Or are there still more things to do there. And if so, should I include them now?

I do not think this fully addresses #655, but it might be close. We need to review the documentation and double-check that the diagnostics are all properly computed (if that's done, please check the box in #655).

[apcraig]

When we average velocity to T gridpoints for diagostics, should it be via an "S" or an "A" average? "S" just averages unmasked values. "A" averages all values including "0s" on land points.

A quick test with a B-grid case, averaging velocity to the T grid with "S", leads to differences in KE of up to ~20% and rms ice speed of up to ~10%. Arctic seems to increase and Antarctic decreases which is also odd. That is surprisingly large all things considered. I don't think "A" will change the results all that much, it will just reduce velocity at the land boundary. Is 10-20% changes what we should expect by making this change?? Sample output from a few timesteps (< is original, > is new computation). Otherwise, all diagnostics, results, and restarts are bit-for-bit.

< tot kinetic energy (J) =    1.05704983559942594E+14   2.17837451740381656E+14
< rms ice speed    (m/s) =        0.12141714861097812       0.13593151480976387
---
> tot kinetic energy (J) =    1.18221584250556125E+14   1.86410938012182125E+14
> rms ice speed    (m/s) =        0.12840463052301013       0.12574466800508025
1954,1955c1954,1955
< tot kinetic energy (J) =    1.03116916338079937E+14   2.10713254698175312E+14
< rms ice speed    (m/s) =        0.11989095333854706       0.13373704987236470
---
> tot kinetic energy (J) =    1.16166934063425047E+14   1.81509234400792750E+14
> rms ice speed    (m/s) =        0.12725144384255277       0.12412383165292681
2036,2037c2036,2037
< tot kinetic energy (J) =    1.02387446679585031E+14   2.07076425696231437E+14
< rms ice speed    (m/s) =        0.11943572466777044       0.13262459221801276
---
> tot kinetic energy (J) =    1.15960377431561109E+14   1.79024040035439594E+14
> rms ice speed    (m/s) =        0.12710589783025453       0.12331457439291411
2118,2119c2118,2119
< tot kinetic energy (J) =    1.03040335377542250E+14   2.07422253853324656E+14
< rms ice speed    (m/s) =        0.11978526883994178       0.13277740763082344
---
> tot kinetic energy (J) =    1.17193638553626375E+14   1.79488965646758125E+14
> rms ice speed    (m/s) =        0.12774732039487935       0.12351377290518413

[JFLemieux73]

Hi Tony,

I think we should use A averaging. When you say 'including "0s" on land points' you mean land-ocean boundary points that have have u,v=0 because of the boundary conditions?

[apcraig]

@JFLemieux73, correct. There are 0s on neighbor points because nothing is computed there. The "A" averaging includes those. The "S" averaging just averages points that are "active". I will switch to "A" averaging. This will reduce the averaged velocity near land boundaries, but otherwise will be the same. It will be interesting to see how this changes the values.

[apcraig]

Output with "A" mapping results (< is original same as above, > is with "A" averaging of u/v). This makes a relatively big difference. All values are now reduced, now up to 25% or even a bit more. I still find this surprising.

< tot kinetic energy (J) =    1.05704983559942594E+14   2.17837451740381656E+14
< rms ice speed    (m/s) =        0.12141714861097812       0.13593151480976387
---
> tot kinetic energy (J) =    9.16786980527338437E+13   1.62428433854588812E+14
> rms ice speed    (m/s) =        0.11307497923412818       0.11737751507186428
1954,1955c1954,1955
< tot kinetic energy (J) =    1.03116916338079937E+14   2.10713254698175312E+14
< rms ice speed    (m/s) =        0.11989095333854706       0.13373704987236470
---
> tot kinetic energy (J) =    8.99214739711308750E+13   1.57901911038422937E+14
> rms ice speed    (m/s) =        0.11195749263441422       0.11577092525216240
2036,2037c2036,2037
< tot kinetic energy (J) =    1.02387446679585031E+14   2.07076425696231437E+14
< rms ice speed    (m/s) =        0.11943572466777044       0.13262459221801276
---
> tot kinetic energy (J) =    8.95471057261088750E+13   1.55700009765719500E+14
> rms ice speed    (m/s) =        0.11169575531991531       0.11500137988797084
2118,2119c2118,2119
< tot kinetic energy (J) =    1.03040335377542250E+14   2.07422253853324656E+14
< rms ice speed    (m/s) =        0.11978526883994178       0.13277740763082344
---
> tot kinetic energy (J) =    9.01660479834729062E+13   1.56250972499639781E+14
> rms ice speed    (m/s) =        0.11205243472702531       0.11524122740140137

[JFLemieux73]

I think the differences come from the fact that you interpolate first and then calculate u2 + v2. I would try to calculate u2 + v2 first and then do the interpolation.

[JFLemieux73]

u2 = u**2...

[apcraig]

How would averaging uu+vv work on the C-grid since u and v are at different locations? I think we want to compute uu and vv separately and map that then combine them? That should work fine for all grids. So the difference would be mapping u/v or mapping u2/v2.

A quick test, averaging uu+vv for the B grid produces much closer results. This is not intuitive to me, interesting.

< tot kinetic energy (J) =    1.05704983559942594E+14   2.17837451740381656E+14
< rms ice speed    (m/s) =        0.12141714861097812       0.13593151480976387
---
> tot kinetic energy (J) =    1.05318296876123922E+14   1.88866104358449875E+14
> rms ice speed    (m/s) =        0.12119486291635775       0.12657003316437856
1954,1955c1954,1955
< tot kinetic energy (J) =    1.03116916338079937E+14   2.10713254698175312E+14
< rms ice speed    (m/s) =        0.11989095333854706       0.13373704987236470
---
> tot kinetic energy (J) =    1.02887899940024375E+14   1.83150175609628500E+14
> rms ice speed    (m/s) =        0.11975774407812759       0.12468364238710151
2036,2037c2036,2037
< tot kinetic energy (J) =    1.02387446679585031E+14   2.07076425696231437E+14
< rms ice speed    (m/s) =        0.11943572466777044       0.13262459221801276
---
> tot kinetic energy (J) =    1.02271550092365047E+14   1.80396438860153156E+14
> rms ice speed    (m/s) =        0.11936810840871509       0.12378633693330938
2118,2119c2118,2119
< tot kinetic energy (J) =    1.03040335377542250E+14   2.07422253853324656E+14
< rms ice speed    (m/s) =        0.11978526883994178       0.13277740763082344
---
> tot kinetic energy (J) =    1.03024821182457625E+14   1.81072573160082594E+14
> rms ice speed    (m/s) =        0.11977625080841743       0.12405744922111049

[apcraig]

There is also a diagnostic for the maximum speed. Should that also be interpolated to the T grid? Currently, we are just evaluating them at the U, N, or E grid points, with the C grid incorrectly combining uvelE^2 + vvelN^2. I guess I'll interpolate the squared velocity to the T point here as well then taking the sqrt.

[apcraig]

For max ice speed, I'm averaging uu and vv to the T grid then adding them and taking the sqrt. For many timesteps, the max speed is reduced by a few percent. For some timesteps, the max speed is reduced by 20%. This is different than the ke or rms ice speed as it represents the value on one grid point, not aggregated over the grid. I think it makes sense that the max speed is reduced when interpolating and that the amount of reduction can vary over time quite a bit. Is this what we want or do we just want to leave it like it was. Some output (< old, > new),

6554c6554
< max ice speed    (m/s) =        0.70513280849380267       0.30503580050552248
---
> max ice speed    (m/s) =        0.66869569762562509       0.29630303924568568

7622c7622
< max ice speed    (m/s) =        0.72585606435592342       0.32659921016949545
---
> max ice speed    (m/s) =        0.56588394034081257       0.28828361334048946

10658c10658
< max ice speed    (m/s) =        0.38427407167820776       0.31749894854027921
---
> max ice speed    (m/s) =        0.36867488792117187       0.27899196551749433

[eclare108213]

I don't see the point of averaging in order to get a maximum ice speed diagnostic. Just use the raw data.

[eclare108213]

Well, I guess the point is that u and v are in different locations on the C-grid, and speed needs both of them. What you're doing seems reasonable, and this is just a diagnostic... A true max would probably interpolate u values to v locations, calculate speed and take the max, then interpolate v values to u locations, calculate speed and take the max, then take the maximum of the two maxes. Not sure it's worth the trouble.

[apcraig]

I agree, I think max speed might as well be the raw data. For the C-grid, what you think about using

sqrt(uvelE*uvelE + vvelN*vvelN)?

If you think we shouldn't do that, I'll do the u to v and v to u calculation as you propose. I do think we should either average all (B, C, CD) max speeds to the T grid or none (as much as possible) for some consistency.

[eclare108213]

For the KE:
A mathematical truth is that
sum(ui^2 + vi^2) <= (sum(ui))^2 + (sum(vi))^2
because of the cross terms when squaring sums on the right hand side. That's what @JFLemieux73 is getting at.
But I think the right hand side is more correct in this case, which is averaging u and v individually to T and then calculating KE there. So I think your initial, larger values are both understandable and more correct. Does that make sense? Feel free to disagree, since I'm not 100% convinced, myself!

[apcraig]

I'm not going to argue for one implementation over another. Just let me know how you want me to do this and I will.

[eclare108213]

I don't really like
sqrt(uvelE*uvelE + vvelN*vvelN)
for the C grid locations. If we're going to fix it, I think we should do it right.

[apcraig]

On KE, I guess I do feel interpolating U and V to the T grid then computing KE is more straight-forward to justify.

[eclare108213]

Yes. Conceptually, if we wanted to sample KE anywhere in the domain (regardless of whether it's a cell center or a node or somewhere in between, or out in the real world), we would sample the u, v and m values from continuous fields and then calculate the KE. Like you say, interpolating squares or other mid-calculation quantities is harder to justify.

[JFLemieux73]

I agree with this. We should interpolate u and v to the T point and then calculate KE. It is simpler like that. I was just trying to explain the differences between the new approach and the old one. Should we also simply get max speed from u and v interpolated at the T point?