update_ocn_f problem

[dabail10]

Not sure when this code came in, but there is a bug here. When update_ocn_f is true, we do need to compute the freshwater and salt due to frazil formation even for ktherm == 2. This was supposed to be "somewhere" else, but I don't see where it is done.

 if (update_ocn_f) then
         if (ktherm <= 1) then
            dfresh = -rhoi*vi0new/dt 
            dfsalt = ice_ref_salinity*p001*dfresh
            fresh  = fresh + dfresh
            fsalt  = fsalt + dfsalt
         ! elseif (ktherm == 2) the fluxes are added elsewhere
         endif
      else ! update_ocn_f = false
         if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
            vi0tmp = fnew*dt / (rhoi*Lfresh)
            dfresh = -rhoi*(vi0new - vi0tmp)/dt
            dfsalt = ice_ref_salinity*p001*dfresh
            fresh  = fresh + dfresh
            fsalt  = fsalt + dfsalt
            frazil_diag = frazil - vi0tmp
         ! elseif ktherm==1 do nothing
         endif
      endif

[eclare108213]

I think it's in ice_comp_mct.F, at least in E3SM. The logic associated with update_ocn_f is extremely convoluted and confusing. I'd like to fix it as part of the Icepack merge effort - will need to discuss how to implement it broadly, since it likely will affect many/most coupled model configurations.

See lines 1973ff and 2555-2556 in
https://github.com/E3SM-Project/E3SM/blob/302618e817d9c4299a213b62b49707fcf2c5b01e/components/mpas-seaice/driver/ice_comp_mct.F.

These changes were made at the same time as removing the temperature limiting mentioned in
#385.

[dabail10]

What was the solution in E3SM here? I had to remove the if block in the update_ocn_f = .true. code. This was done in #418 . So, it is answer changing as the "bug" here was the the fresh water and salt flux were not computed for frazil formation here. This is done somewhere else in E3SM.

[eclare108213]

Here is the solution currently proposed, which fixes a conservation issue in E3SM using the newer Icepack code and is backward compatible for all other codes, as far as we know. It adds a new namelist flag, cpl_frazil, with default combinations listed below for some models. This is being implemented in Icepack now -- if anyone has concerns, please speak up soon.

   dfresh = c0
   dfsalt = c0
   if (cpl_frazil == 'external') then
      ! do nothing here, calculations are in the coupler or elsewhere
   else
      if (update_ocn_f) then
         dfresh = -rhoi*vi0new/dt
      elseif (cpl_frazil == 'fresh_ice_correction' .and. ktherm == 2) then
         ! correct frazil fluxes for mushy
         vi0tmp = fnew*dt / (rhoi*Lfresh) ! ocn/cpl assumes frazil volume is pure, fresh ice
         dfresh = -rhoi*(vi0new - vi0tmp)/dt
         frazil_diag = frazil - vi0tmp
!     else 
!        do nothing - other correction options could be implemented in the future             
      endif

      if (saltflux_option == 'prognostic') then
         dfsalt = Si0new*p001*dfresh
      else
         dfsalt = ice_ref_salinity*p001*dfresh
      endif
      fresh  = fresh + dfresh
      fsalt  = fsalt + dfsalt
      
   endif

Icepack default:

cpl_frazil = 'fresh_ice_correction'
update_ocn_f = .false.
ktherm = 2

E3SM:

cpl_frazil = 'external'
update_ocn_f = .true. or .false. or neither (icepack will assign its default value)
ktherm = 2

CESM, RASM:

cpl_frazil = 'fresh_ice_correction'
update_ocn_f = .false.
ktherm = 2

[dabail10]

The default for CESM2+ is:

update_ocn_f = .true.
saltflux_option = 'prognostic'

[eclare108213]

And so I think you'd also need cpl_frazil = 'internal' (or something other than 'external' or 'fresh_ice_correction') in the revised version.

Edit: Actually, you could just not set cpl_frazil and it would work, although it might be confusing to anyone looking at namelist settings in log files.

[dabail10]

Agreed.

[eclare108213]

See #458
Can we close this issue now?

[apcraig]

I think we need an update to CICE? Don't we want to modify the implementation so it sets these parameters at initialization and drops the argument passed at runtime (even though it continues to be allowed). At some point, we'll want to remove that optional argument, maybe when some other Icepack interfaces change. We'll want to remember to do that.

[apcraig]

I'm working on the CICE update. I'm a little confused about the CICE diagnostics and history code. @dabail10, can you grep for update_ocn_f in

cicecore/cicedyn/analysis/ice_history.F90
cicecore/cicedyn/general/ice_flux.F90

and let me know what you think. I think the following changes are needed but not sure. I guess I'm wondering whether there is a better way to implement this so we don't need any "if" tests here. Why can't CICE/Icepack set variables in the physics that can be used directly in history/diagnostics without "if" checks in history/diagnostics. For some cases, those variables would probably be 0. The problem now is that anytime there are changes in the frazil ice in Icepack, it looks like those changes have to be "duplicated" in the CICE history and diagnostics code, that's not particularly friendly. If we need to add some additional arguments to Icepack interfaces to generalize, I favor that over what we're doing now.

--- a/cicecore/cicedyn/analysis/ice_diagnostics.F90
-         if (ktherm == 2 .and. .not.update_ocn_f) &
+         if (.not. update_ocn_f .and. ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
             work1(:,:,:) = (frazil(:,:,:)-frazil_diag(:,:,:))*rhoi/dt
+         endif

--- a/cicecore/cicedyn/analysis/ice_history.F90
 !                Add in frazil flux
                  if (.not. update_ocn_f) then
-                 if ( ktherm == 2) then
-                    dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
-                 else
-                    dfresh = -rhoi*frazil(i,j,iblk)/dt
-                 endif
+                    if ( ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
+                       dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
+                    else
+                       dfresh = -rhoi*frazil(i,j,iblk)/dt
+                    endif
                  endif

 !                Add in frazil flux
                  if (.not. update_ocn_f) then
-                 if ( ktherm == 2) then
-                    dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
-                 else
-                    dfresh = -rhoi*frazil(i,j,iblk)/dt
-                 endif
+                    if ( ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
+                       dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
+                    else
+                       dfresh = -rhoi*frazil(i,j,iblk)/dt
+                    endif

[dabail10]

I agree with this logic. The frazil_diag should only be used when ktherm ==2 and cpl_frazil == 'fresh_ice_correction'. We can rename if this makes sense.

[apcraig]

Thanks @dabail10, is there a way we can avoid all the if checks in the diagnostics and history modules? Like if frazil and frazil_diag were set to 0 unless they need to have a value. Could we just have a generic

dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt

without any if checks anywhere and if frazil and/or frazil_diag were 0 except under certain conditions, we'd get the right results all the time? Is that already happening in the code via Icepack? I noticed that while frazil is set to 0 at initialization, it's NOT reset to 0 everytime icepack is called. That seems like it could lead to some confusion.

[apcraig]

My hope here is to NOT have code like this in CICE that has to be "kept up" with changes in Icepack. Icepack should provide the "diagnostics" with all the logic needed to make sure they are "correct" for all conditions, then CICE just uses it. Thoughts?

[dabail10]

This is a nice idea. Just resetting frazil_diag to zero each time. In terms of the diagnostics code, this is tricky. For CICE, these are total budgets for the NH and SH. Whereas icepack is point by point. Not sure how to take advantage of this.

[apcraig]

@dabail10, are the diagnostics and history self consistent? In diagnostics, we have

         work1(:,:,:) = frazil(:,:,:)*rhoi/dt
         if (.not. update_ocn_f .and. ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
            work1(:,:,:) = (frazil(:,:,:)-frazil_diag(:,:,:))*rhoi/dt
         endif

in history, we have

                 if (.not. update_ocn_f) then
                    if ( ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
                       dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
                    else
                       dfresh = -rhoi*frazil(i,j,iblk)/dt
                    endif
                 endif

In diagnostics, "work1" is always defined based on frazil with one special case. In history, dfresh is zero except in the two special cases. Just asking because I might expect the logic to be the same but it's not.

[apcraig]

Looking at the Icepack implementation now. I'm not sure the diagnostics/history are consistent with Icepack either which looks like this

      if (cpl_frazil == 'external') then
         ! do nothing here, calculations are in the coupler or elsewhere                                                                                        
      else
         if (update_ocn_f) then
            dfresh = -rhoi*vi0new/dt
         elseif (cpl_frazil == 'fresh_ice_correction' .and. ktherm == 2) then
            ! correct frazil fluxes for mushy                                                                                                                   
            vi0tmp = fnew*dt / (rhoi*Lfresh) ! ocn/cpl assumes frazil volume is pure, fresh ice                                                                 
            dfresh = -rhoi*(vi0new - vi0tmp)/dt
            frazil_diag = frazil - vi0tmp
!        else                                                                                                                                                   
!           do nothing - other correction options could be implemented in the future                                                                            
         endif

So that sort of looks like diagnostics (except the update_ocn_f check is different).

[apcraig]

Also, there is an adjustment to frazil in Icepack due to fsd AFTER the definition of frazil and dfresh above. It seems like we are computing the history/diagnostics in CICE based on frazil without taking into account fsd changes correctly.

[dabail10]

This is not exactly right. If you look at the code in add_new_ice, vi0tmp is only computed when ktherm == 2 and the cpl_frazil is 'fresh_ice_correction'. So, this is mean to "correct" the freshwater and salt that is sent to the coupler. So, frazil_diag = frazil - vi0tmp is the "excess" frazil created by the mushy layer physics (ktherm == 2). So, the history variable output and diagnostic budgets should only include whatever is actually computed in CICE. Maybe things are inconsistent here now that I look at it. Marika and I went over it though and the salt and freshwater budgets do balance.

[dabail10]

Just seeing the FSD comment. I had never looked at the frazil adjustment in the FSD, so this could definitely be a problem.

[apcraig]

OK, there is also the challenge that CICE doesn't really test all the update_ocn_f modes, and we have not run with FSD all that much. I continue to hope that Icepack can provide a diagnostics (which will sometimes be zero) that can be added to the CICE diagnostics and history appropriately. Hopefully, that's just one diagnostic. We could redefine the frazil_diag term to be the "whole thing", so diagnostics would look like

work1(:,:,:) = frazil_diag(:,:,:)*rhoi/dt

and history would be

dfresh = frazil_diag(i,j,iblk)*rhoi/dt

That would be great if possible. Then all the logic and complexity would live in Icepack and we wouldn't have to worry about this stuff on the CICE side. Icepack would just have to be aware that frazil_diag would have to be kept up to date.

Again, not sure if that's possible, just proposing it for now.

[dabail10]

Good thoughts. I am running with the FSD in CESM now. However, I have update_ocn_f = .true. as we are doing all of the freshwater and salt fluxes in CICE when coupled to MOM6. There are still runs coupled to POP with the CESM3 development code, so I think we still need something.

[apcraig]

Good points. I'm not even sure what the diagnostics and history should do IF some of this stuff is computed outside CICE. Are those coupled fields taken into account in CICE/Icepack and then do they break conservation? I guess when things are computed "external", then there is no frazil ice generated in Icepack? So maybe that's OK and self consistent??

[dabail10]

Not exactly. The subroutine add_new_ice always computes the frazil volume based on the heat flux from the ocean. This is all about whether the salt and fresh water fluxes are computed in Icepack or not. So, yes the budgets are impacted by this. If the fresh water and salt are exchanged with the ocean then they are added to the budgets. I kind of like using "frazil_diag" here as zero or non-zero. Then we could always add frazil_diag in the budgets.

[apcraig]

Thanks @dabail10, how should we proceed? Should I try to update the code in Icepack and CICE to handle the frazil_diag change? Or should I just update the code in CICE to keep the history/diagnostics consistent with what we have now (which may have some problems) and then we can tackle the frazil_diag separately. I'm updating some the ocn_update_f implementation in CICE at the moment, but it could be a two-step process. But happy to try to do it all now with some suggestions/guidance. Also, we need to decide how to test and validate, does some of that have to be done in a coupled system?

[dabail10]

I will prototype something this week and make sure it meets the needs of the coupled models. I am out of the office M-W, but back on Thursday.

[apcraig]

Thanks @dabail10. In the mean time, I have created a PR for the initial changes. See CICE-Consortium/CICE#881

[apcraig]

I am going to close CICE-Consortium/CICE#881, we should continue discussion/design, and then create a new PR (maybe in both Icepack and CICE) to address this issue.

[apcraig]

I created a PR, CICE-Consortium/CICE#889, in CICE to update the basic implementation. I would still like to see the Icepack/CICE frazil diagnostics updated as proposed above if possible. @dabail10, let me know if I can help with the implementation or testing.

[dabail10]

Has this been fixed?

[eclare108213]

I think the fundamental issue has been fixed, but there are a number of other things in the conversation above that still need to be looked at. Maybe pull those out as separate issue(s) and then we can close this one.