hfds is broken

[StephenGriffies]

The CMIP6 diagnostic hfds is broken. This diagnostic measures the net heat through the ocean surface. The problem is related to incorrectly diagnosed heat associated with water transfer across the ocean surface.

[StephenGriffies]

The following documents investigations from 19May from Alistair and Steve.

extractFluxes1d is called inside of calculateBuoyancyFlux1d. This
call modifies the diagnostics heat_content_massout and
heat_content_massin.

calculateBuoyancyFlux is called from inside of MOM_diabatic.F90 on
line 570, with KPP option true. It is also called inside of
MOM_diabatic.F90 if CS%use_energetic_PBL=.true. The EBPL call is
presumably not needed. Indeed, the all to calculateBuoyancyFlux passes
some KPP flags.

We should perhaps remove the calculateBuoyancyFlux
call when using EPBL. This hypothesis needs to be pursued.

[StephenGriffies]

Update on this issue. There may be a code bug with EPBL. At least, I am puzzled why it calls extractfluxes1d twice.

Details

From inside of MOM_diabatic, the EPBL logic first calls applyBoundaryFluxesInOut (which is inside MOM_diabatic_aux). Next, it calls calculateBuoyancyFlux2d (which is inside MOM_forcing_type). Each of these routines in turn calls extractFluxes1d (which is inside of MOM_forcing_type). The two calls to extractFluxes1d is a problem for diagnostics and potential for the prognostic solution.

Notes:

1/ Each call to extractFluxes1d reinitializes the flux diagnostics, with the second initialization leaving some of the diagnostics spuriously set to zero. That is the root problem with the hfds diagnostic that started this thread.

2/ Two calls to extractFluxes1d is inefficient at best, and potentially a code bug at worse.

3/ KPP calls calculateBuoyancyFlux2d. Yet it does not call applyBoundaryFluxesInOut. So the KPP method only calls extractFluxes1d, meaning the flux diagnostics are fine.

4/ Inside of MOM_diabatic_driver, the EPBL call to calculateBuoyancyFlux2d makes use of a KPP type to pass three boundary fluxes through the subroutine:

CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux

Use of the KPP type for EPBL is confusing, since EPBL has nothing to do with KPP.

As I am concerned that there is more than a mere diagnostic bug going on here, I need advise from the EPBL developers to help resolve the above questions.

[Hallberg-NOAA]

Steve and I looked through this, and it definitely appears that the heat flux associated with river or calving heat fluxes are being double incremented. One solution might be to optionally have applyBoundaryFluxesInOut return the surface buoyancy flux, which would eliminate the need for the subsequent call to calculateBuoyancyFlux2d. A second option would be to add an optional flag to calculateBuoyancyFlux2d to avoid the side-effect of incrementing the accumulated heat input.

Regardless of which option we choose, there is a second mystery in that a 3-d buoyancy flux array is being passed into a 2-d buoyancy flux argument to energetic_PBL without triggering a compile-time error and without an explanatory comment. I know that such a thing is done routinely in C, but in F90 I think that such array reconfiguration is a bad idea and should be avoided in MOM6 code.

[adcroft]

I've been looking at this too and thought I could correct the diagnostics by re-ordering the call to the buoyancy calculation as follows:

--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -746,13 +746,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, G, GV, CS)
         h_prebound(i,j,k) = h(i,j,k)
     enddo ; enddo ; enddo
     if (CS%use_energetic_PBL) then
+      call calculateBuoyancyFlux2d(G, GV, fluxes, CS%optics, h, tv%T, tv%S, tv, &
+           CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux)
+
       call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, &
                           h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, &
                           CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS)
 
-      call calculateBuoyancyFlux2d(G, GV, fluxes, CS%optics, h, tv%T, tv%S, tv, &
-           CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux)
-

but this changed answers - which I think is because h, tv%T and tv%S are correctly different before and after applyBoundaryFluxesInOut(). CS%KPP_buoy_flux depends on those state variables which is why ePBL does something different.

@breichl I think the buoyancy call should be before the apply-fluxes so that the buoyancy flux is more consistent with the actual forcing used in ePBL?

@Hallberg-NOAA's suggestion of having applyBoundaryFluxesInOut() return the buoyancy flux would be the cleaner way to do things except the buoyancy flux calculation adds up fluxes differently for reasons I forget but related to assumptions in KPP. So taking this approach will also change answers but I agree we should do it this way.