Check strocnxT implementation

[apcraig]

I'm looking at the implementation of strocnxT, strocnyT in CICE. This is defined as the ice/ocn stress on the T grid and is something that is passed out of the model for coupling. It is computed by taking the strocnxU, strocnyU, normalizing the ice fraction on the U grid (aiU) and then averaging to the T grid. Normalization by ice fraction is needed for the coupling field, many coupling fields are normalized the same way, see subroutine scale_fluxes. However, strocnxT and strocnyT are computed in the dynamics, not in scale fluxes. I'm planning to clean up the implementation, see #660, and move the computation into scale fluxes.

I do have one question. I think there may be a bug. strocnxT and strocnyT are also passed into icepack_step_therm1. They are used in the frzmlt_bottom_lateral subroutine in computation of ustar,

         ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2                                                                                                                                           
         ustar = sqrt (sqrt(strocnxT**2+strocnyT**2)/rhow)

My question is whether strocnxT and strocnyT here should be the stresses normalized by the fraction or not. In the current implementation, they are. But I'm pretty sure that's wrong. Thoughts? @eclare108213 @dabail10

FYI, strocnxT and strocnyT are computed as follows,

      ! strocn computed on U, N, E as needed. Map strocn U divided by aiU to T                                                                                                                                
      ! TODO: This should be done elsewhere as part of generalization?                                                                                                                                        
      ! TODO: Rename strocn[x,y]T since it's different than strocn[x,y][U,N,E]                                                                                                                                
      ! conservation requires aiU be divided before averaging                                                                                                                                                 
      work1 = c0
      work2 = c0
      !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) SCHEDULE(runtime)                                                                                                                                                
      do iblk = 1, nblocks
      do ij = 1, icellu(iblk)
         i = indxui(ij,iblk)
         j = indxuj(ij,iblk)
         work1(i,j,iblk) = strocnxU(i,j,iblk)/aiU(i,j,iblk)
         work2(i,j,iblk) = strocnyU(i,j,iblk)/aiU(i,j,iblk)
      enddo
      enddo
      !$OMP END PARALLEL DO                                                                                                                                                                                   
      call ice_HaloUpdate (work1,              halo_info, &
                           field_loc_NEcorner, field_type_vector)
      call ice_HaloUpdate (work2,              halo_info, &
                           field_loc_NEcorner, field_type_vector)
      call grid_average_X2Y('F', work1, 'U', strocnxT, 'T')    ! shift                                                                                                                                        
      call grid_average_X2Y('F', work2, 'U', strocnyT, 'T')

[apcraig]

@eclare108213 @dabail10 I'd like to work on this task, could someone have a quick look at my question above if you have a chance. Thanks.

[eclare108213]

This is a little complicated. The calculation is in the dynamics rather than scale_fluxes for a couple of reasons (probably). Here's what I remember, with the caveat that a lot of things have changed and so some of this might not still be true.

First, just to clarify the language: The units of ice volume are m^3 of ice per m^2 of grid cell area = grid-cell-avg ice thickness in m. Dividing by aice, which is unitless, gives us m^3 of ice per m^2 of ice area = "actual" ice thickness averaged over the total ice area, in m.

Second, the calculation of strocnx,y is (or was) split -- part of it done in the momentum equation and the rest of it after the subcycling, for efficiency. That required extreme care to make sure the terms were consistent, and I have a vague recollection that this changed in the code to not be split. So now the scaling might not need to be done in the dynamics anymore. However -

There was a huge community argument 20 years ago about whether the stress terms in the momentum equation should be multiplied by the ice area or not. In the original (Hibler) momentum equations, they were not. The Connolley et al (2004) paper finally cleared that up (answer: yes). At the time I carefully traced all of the area-weighted scalings through the code to make sure they were working as expected. But that was for the dynamics.

You're asking about the thermodynamics. The question in that case is whether the calculations should involve the entire grid cell or just the portion under the ice. Generally speaking, the assumption is that the vertical thermo calculations are performed for a point in the grid cell, or equivalently as if ice [in a given thickness category] covers the whole grid cell (i.e. divided by the ice [category] area, which I think is what you mean by 'normalized'), and then the results are scaled as needed using ice area. But that means that the input fluxes for the calculations also have to be appropriately scaled for the calculations.

In the dynamics, we assume that the ice volume spread evenly over the whole grid cell, which means multiplying stresses by aice to get the grid-cell average. In the thermo, we assume the ice area is concentrated into one point, so for example we use the actual thickness, hi = vice/aice, rather than the grid-cell-averaged thickness vice. And so dividing the thermo forcing fluxes by aice (if that's what you mean by 'normalizing') might be okay, conceptually. How the scaling translates from strocn into ustar and then into the actual heat flux might need to be reviewed, though.

Does that make sense? And if so, do you still think there's a bug? If the answer to that is yes, then I'll look at the code more closely.

Sorry this is a long-winded answer that doesn't completely answer your question, but hopefully the background is helpful. Thanks for bringing this up -- it's good to review questionable things like this, just in case.

[dabail10]

I guess had not realized that ustar was computed from strocn before. I believe you are correct that we should not be dividing by ice area here. I had always though this was a coupler field and that was it. I see that this was in CICE5 as well, so I think it has been there awhile.

[apcraig]

Thanks for the info. I don't have a good sense of what's correct in the thermodynamics. It sounds a bit complicated? I agree with @dabail10 in the sense that it looks to me like strocnx,yT was coded as a diagnostic and coupling field, and the normalization (divide by aiU, required for coupling) is happening in the dynamics because aiU is readily available there. But then somehow strocnx,yT got passed into Icepack with the normalization applied. That's sort of what perked my curiosity/concern.

If @dabail10 is correct, then we should not divide strocnx,yT by aiU when we pass it into Icepack. If that's true, it'll change answers and we can move the "normalization" by aiU into scale fluxes quite cleanly. I would add aiU to the ice_state module data to support that. Are we comfortable that's correct?

[eclare108213]

Think about it this way, rather than in terms of whether it's a coupling field or not: the stress terms in the dynamics are calculated and used with aice as a multiplicative factor in the equations (see dyn_finish in ice_dyn_shared.F90). Dividing by aiU then takes the stress term back to its "normal" form without the aice factor.