Merge branch 'main' into main-to-gfdl_2024-05-31

Resolved conflicts in src/core/MOM_checksum_packages.F90 with
collision between switch to unscale= and the addition of new
checksums with previous scale= argument.

src/core/MOM.F90

@@ -3299,17 +3299,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
 
   CS%mixedlayer_restrat = mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, &
                                                   CS%mixedlayer_restrat_CSp, restart_CSp)
-  if (CS%mixedlayer_restrat) then
-    if (GV%Boussinesq .and. associated(CS%visc%h_ML)) then
-      ! This is here to allow for a transition of restart files between model versions.
-      call get_param(param_file, "MOM", "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
-                     default=.false., do_not_log=.true.)
-      if (MLE_use_PBL_MLD .and. .not.query_initialized(CS%visc%h_ML, "h_ML", restart_CSp) .and. &
-          associated(CS%visc%MLD)) then
-        do j=js,je ; do i=is,ie ; CS%visc%h_ML(i,j) = GV%Z_to_H * CS%visc%MLD(i,j) ; enddo ; enddo
-      endif
+
+  if (GV%Boussinesq .and. associated(CS%visc%h_ML)) then
+    ! This is here to allow for a transition of restart files between model versions.
+    call get_param(param_file, "MOM", "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
+                   default=.false., do_not_log=.true.)
+    if (MLE_use_PBL_MLD .and. .not.query_initialized(CS%visc%h_ML, "h_ML", restart_CSp) .and. &
+        associated(CS%visc%MLD)) then
+      do j=js,je ; do i=is,ie ; CS%visc%h_ML(i,j) = GV%Z_to_H * CS%visc%MLD(i,j) ; enddo ; enddo
     endif
+  endif
 
+  if (CS%mixedlayer_restrat) then
     if (.not.(bulkmixedlayer .or. CS%use_ALE_algorithm)) &
       call MOM_error(FATAL, "MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
     ! When DIABATIC_FIRST=False and using CS%visc%ML in mixedlayer_restrat we need to update after a restart
@@ -4056,7 +4057,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
     endif
   endif
 
-  if (CS%debug) call MOM_surface_chksum("Post extract_sfc", sfc_state, G, US, haloshift=0)
+  if (CS%debug) call MOM_surface_chksum("Post extract_sfc", sfc_state, G, US, haloshift=0, symmetric=.true.)
 
   ! Rotate sfc_state back onto the input grid, sfc_state_in
   if (CS%rotate_index) then

src/core/MOM_checksum_packages.F90

@@ -167,7 +167,7 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, US, haloshift, symmetric)
   logical :: sym
 
   sym = .false. ; if (present(symmetric)) sym = symmetric
-  hs = 1 ; if (present(haloshift)) hs = haloshift
+  hs = 0 ; if (present(haloshift)) hs = haloshift
 
   if (allocated(sfc_state%SST)) call hchksum(sfc_state%SST, mesg//" SST", G%HI, haloshift=hs, &
                                              unscale=US%C_to_degC)
@@ -182,6 +182,14 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, US, haloshift, symmetric)
                   unscale=US%L_T_to_m_s)
   if (allocated(sfc_state%frazil)) call hchksum(sfc_state%frazil, mesg//" frazil", G%HI, &
                                                 haloshift=hs, unscale=US%Q_to_J_kg*US%RZ_to_kg_m2)
+  if (allocated(sfc_state%melt_potential)) call hchksum(sfc_state%melt_potential, mesg//" melt_potential", &
+                      G%HI, haloshift=hs, unscale=US%Q_to_J_kg*US%RZ_to_kg_m2)
+  if (allocated(sfc_state%ocean_mass)) call hchksum(sfc_state%ocean_mass, mesg//" ocean_mass", &
+                      G%HI, haloshift=hs, unscale=US%RZ_to_kg_m2)
+  if (allocated(sfc_state%ocean_heat)) call hchksum(sfc_state%ocean_heat, mesg//" ocean_heat", &
+                      G%HI, haloshift=hs, unscale=US%C_to_degC*US%RZ_to_kg_m2)
+  if (allocated(sfc_state%ocean_salt)) call hchksum(sfc_state%ocean_salt, mesg//" ocean_salt", &
+                      G%HI, haloshift=hs, unscale=US%S_to_ppt*US%RZ_to_kg_m2)
 
 end subroutine MOM_surface_chksum
 

src/equation_of_state/MOM_EOS.F90

@@ -1526,6 +1526,7 @@ subroutine EOS_init(param_file, EOS, US)
                  "If true, use a bug in the calculation of the second derivatives of density "//&
                  "with temperature and with temperature and pressure that causes some terms "//&
                  "to be only 2/3 of what they should be.", default=.false.)
+    call EOS_manual_init(EOS, form_of_EOS=EOS_WRIGHT, use_Wright_2nd_deriv_bug=EOS%use_Wright_2nd_deriv_bug)
   endif
 
   EOS_quad_default = .not.((EOS%form_of_EOS == EOS_LINEAR) .or. &
@@ -1645,6 +1646,8 @@ subroutine EOS_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Co
     select type (t => EOS%type)
       type is (linear_EOS)
         call t%set_params_linear(Rho_T0_S0, dRho_dT, dRho_dS)
+      type is (buggy_Wright_EOS)
+        call t%set_params_buggy_Wright(use_Wright_2nd_deriv_bug)
     end select
   endif
   if (present(form_of_TFreeze))  EOS%form_of_TFreeze = form_of_TFreeze

src/equation_of_state/MOM_EOS_Wright.F90

@@ -12,6 +12,7 @@ module MOM_EOS_Wright
 public buggy_Wright_EOS
 public int_density_dz_wright, int_spec_vol_dp_wright
 public avg_spec_vol_buggy_Wright
+public set_params_buggy_Wright
 
 !>@{ Parameters in the Wright equation of state using the reduced range formula, which is a fit to the UNESCO
 !    equation of state for the restricted range: -2 < theta < 30 [degC], 28 < S < 38 [PSU], 0  < p < 5e7 [Pa].
@@ -39,6 +40,8 @@ module MOM_EOS_Wright
 !> The EOS_base implementation of the Wright 1997 equation of state with some bugs
 type, extends (EOS_base) :: buggy_Wright_EOS
 
+  real :: three = 3.0 !< A constant that can be adjusted to recreate some bugs [nondim]
+
 contains
   !> Implementation of the in-situ density as an elemental function [kg m-3]
   procedure :: density_elem => density_elem_buggy_Wright
@@ -59,6 +62,9 @@ module MOM_EOS_Wright
   !> Implementation of the range query function
   procedure :: EOS_fit_range => EOS_fit_range_buggy_Wright
 
+  !> Instance specific function to set internal parameters
+  procedure :: set_params_buggy_Wright => set_params_buggy_Wright
+
   !> Local implementation of generic calculate_density_array for efficiency
   procedure :: calculate_density_array => calculate_density_array_buggy_Wright
   !> Local implementation of generic calculate_spec_vol_array for efficiency
@@ -228,13 +234,16 @@ elemental subroutine calculate_density_second_derivs_elem_buggy_Wright(this, T,
   real :: z11    ! A local work variable [Pa m2 s-2 PSU-1] = [kg m s-4 PSU-1]
   real :: z2_2   ! A local work variable [m4 s-4]
   real :: z2_3   ! A local work variable [m6 s-6]
+  real :: six    ! A constant that can be adjusted from 6. to 4. to recreate a bug [nondim]
 
   ! Based on the above expression with common terms factored, there probably exists a more numerically stable
   ! and/or efficient expression
 
+  six = 2.0*this%three  ! When recreating a bug from the original version of this routine, six = 4.
+
   z0 = T*(b1 + b5*S + T*(b2 + b3*T))
   z1 = (b0 + pressure + b4*S + z0)
-  z3 = (b1 + b5*S + T*(2.*b2 + 2.*b3*T)) ! BUG: This should be z3 = b1 + b5*S + T*(2.*b2 + 3.*b3*T)
+  z3 = (b1 + b5*S + T*(2.*b2 + this%three*b3*T))  ! When recreating a bug here this%three = 2.
   z4 = (c0 + c4*S + T*(c1 + c5*S + T*(c2 + c3*T)))
   z5 = (b1 + b5*S + T*(b2 + b3*T) + T*(b2 + 2.*b3*T))
   z6 = c1 + c5*S + T*(c2 + c3*T) + T*(c2 + 2.*c3*T)
@@ -249,8 +258,7 @@ elemental subroutine calculate_density_second_derivs_elem_buggy_Wright(this, T,
 
   drho_ds_ds = (z10*(c4 + c5*T) - a2*z10*z1 - z10*z7)/z2_2 - (2.*(c4 + c5*T + z9*z10 + a2*z1)*z11)/z2_3
   drho_ds_dt = (z10*z6 - z1*(c5 + a2*z5) + b5*z4 - z5*z7)/z2_2 - (2.*(z6 + z9*z5 + a1*z1)*z11)/z2_3
-  ! BUG: In the following line: (2.*b2 + 4.*b3*T) should be (2.*b2 + 6.*b3*T)
-  drho_dt_dt = (z3*z6 - z1*(2.*c2 + 6.*c3*T + a1*z5) + (2.*b2 + 4.*b3*T)*z4 - z5*z8)/z2_2 - &
+  drho_dt_dt = (z3*z6 - z1*(2.*c2 + 6.*c3*T + a1*z5) + (2.*b2 + six*b3*T)*z4 - z5*z8)/z2_2 - &
                   (2.*(z6 + z9*z5 + a1*z1)*(z3*z4 - z1*z8))/z2_3
   drho_ds_dp = (-c4 - c5*T - 2.*a2*z1)/z2_2 - (2.*z9*z11)/z2_3
   drho_dt_dp = (-c1 - c5*S - T*(2.*c2 + 3.*c3*T) - 2.*a1*z1)/z2_2 - (2.*z9*(z3*z4 - z1*z8))/z2_3
@@ -933,6 +941,23 @@ subroutine calculate_spec_vol_array_buggy_Wright(this, T, S, pressure, specvol,
 end subroutine calculate_spec_vol_array_buggy_Wright
 
 
+!> Set coefficients that can correct bugs un the buggy Wright equation of state.
+subroutine set_params_buggy_Wright(this, use_Wright_2nd_deriv_bug)
+  class(buggy_Wright_EOS), intent(inout) :: this !< This EOS
+  logical, optional, intent(in)    :: use_Wright_2nd_deriv_bug !< If true, use a buggy
+                           !! buggy version of the calculations of the second
+                           !! derivative of density with temperature and with temperature and
+                           !! pressure.  This bug is corrected in the default version.
+
+  this%three = 3.0
+  if (present(use_Wright_2nd_deriv_bug)) then
+    if (use_Wright_2nd_deriv_bug) then ; this%three = 2.0
+    else  ; this%three = 3.0 ; endif
+  endif
+
+end subroutine set_params_buggy_Wright
+
+
 !> \namespace mom_eos_wright
 !!
 !! \section section_EOS_Wright Wright equation of state

src/framework/MOM_diag_mediator.F90

@@ -14,7 +14,6 @@ module MOM_diag_mediator
 use MOM_diag_manager_infra, only : send_data_infra, MOM_diag_field_add_attribute, EAST, NORTH
 use MOM_diag_manager_infra, only : register_diag_field_infra, register_static_field_infra
 use MOM_diag_manager_infra, only : get_MOM_diag_field_id, DIAG_FIELD_NOT_FOUND
-use MOM_diag_manager_infra, only : diag_send_complete_infra
 use MOM_diag_remap,       only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask
 use MOM_diag_remap,       only : diag_remap_init, diag_remap_end, diag_remap_do_remap
 use MOM_diag_remap,       only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field
@@ -2079,10 +2078,8 @@ end subroutine enable_averages
 subroutine disable_averaging(diag_cs)
   type(diag_ctrl), intent(inout) :: diag_CS !< Structure used to regulate diagnostic output
 
-  call diag_send_complete_infra()
   diag_cs%time_int = 0.0
   diag_cs%ave_enabled = .false.
-
 end subroutine disable_averaging
 
 !> Call this subroutine to determine whether the averaging is

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -676,6 +676,14 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim
     call cpu_clock_begin(id_clock_kpp)
     ! total vertical viscosity in the interior is represented via visc%Kv_shear
 
+    ! NOTE: The following do not require initialization, but their checksums do
+    !   require initialization, and past versions were initialized to zero.
+    KPP_NLTheat(:,:,:) = 0.
+    KPP_NLTscalar(:,:,:) = 0.
+    KPP_buoy_flux(:,:,:) = 0.
+    KPP_temp_flux(:,:) = 0.
+    KPP_salt_flux(:,:) = 0.
+
     ! KPP needs the surface buoyancy flux but does not update state variables.
     ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
     ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux
@@ -1285,6 +1293,14 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
       visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
     enddo ; enddo ; enddo
 
+    ! NOTE: The following do not require initialization, but their checksums do
+    !   require initialization, and past versions were initialized to zero.
+    KPP_NLTheat(:,:,:) = 0.
+    KPP_NLTscalar(:,:,:) = 0.
+    KPP_buoy_flux(:,:,:) = 0.
+    KPP_temp_flux(:,:) = 0.
+    KPP_salt_flux(:,:) = 0.
+
     ! KPP needs the surface buoyancy flux but does not update state variables.
     ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
     ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux
@@ -1890,6 +1906,15 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e
 
   if (CS%useKPP) then
     call cpu_clock_begin(id_clock_kpp)
+
+    ! NOTE: The following do not require initialization, but their checksums do
+    !   require initialization, and past versions were initialized to zero.
+    KPP_NLTheat(:,:,:) = 0.
+    KPP_NLTscalar(:,:,:) = 0.
+    KPP_buoy_flux(:,:,:) = 0.
+    KPP_temp_flux(:,:) = 0.
+    KPP_salt_flux(:,:) = 0.
+
     ! KPP needs the surface buoyancy flux but does not update state variables.
     ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
     ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -1180,7 +1180,8 @@ subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV)
 
   ! Each cell extends from x=-1/2 to 1/2, and has a topography
   ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12.
-  crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3
+  !crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3
+  crv_3 = (Dp + Dm - (2.0*D_vel)) ; crv = 3.0*crv_3
   slope = Dp - Dm
 
   ! Calculate the volume above which the entire cell is open and the volume at which the
@@ -1207,10 +1208,13 @@ subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV)
       !   vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L)).
       if (a2x48_apb3*vol_below(K) < 1e-8) then ! Could be 1e-7?
         ! There is a very good approximation here for massless layers.
-        L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0)
+        !L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0)
+        L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + (ax2_3apb*L0))
       else
+        !L(K) = apb_4a * (1.0 - &
+        !         2.0 * cos(C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3))
         L(K) = apb_4a * (1.0 - &
-                 2.0 * cos(C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3))
+                 2.0 * cos(C1_3*acos((a2x48_apb3*vol_below(K)) - 1.0) - C2pi_3))
       endif
       ! To check the answers.
       ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K)

src/tracer/MOM_tracer_advect.F90

@@ -34,6 +34,8 @@ module MOM_tracer_advect
   logical :: debug                 !< If true, write verbose checksums for debugging purposes.
   logical :: usePPM                !< If true, use PPM instead of PLM
   logical :: useHuynh              !< If true, use the Huynh scheme for PPM interface values
+  logical :: useHuynhStencilBug = .false. !< If true, use the incorrect stencil width.
+                                   !! This is provided for compatibility with legacy simuations.
   type(group_pass_type) :: pass_uhr_vhr_t_hprev !< A structure used for group passes
 end type tracer_advect_CS
 
@@ -94,6 +96,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
                                        ! can be simply discarded [H L2 ~> m3 or kg].
 
   real :: landvolfill                   ! An arbitrary? nonzero cell volume [H L2 ~> m3 or kg].
+  logical :: use_PPM_stencil            ! If true, use the correct PPM stencil width.
   real :: Idt                           ! 1/dt [T-1 ~> s-1].
   logical :: domore_u(SZJ_(G),SZK_(GV))  ! domore_u and domore_v indicate whether there is more
   logical :: domore_v(SZJB_(G),SZK_(GV)) ! advection to be done in the corresponding row or column.
@@ -112,7 +115,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
   landvolfill = 1.0e-20         ! This is arbitrary, but must be positive.
-  stencil = 2                   ! The scheme's stencil; 2 for PLM and PPM:H3
+  stencil = 2                   ! The scheme's stencil; 2 for PLM
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_advect: "// &
        "tracer_advect_init must be called before advect_tracer.")
@@ -123,8 +126,8 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
   x_first = (MOD(G%first_direction,2) == 0)
 
   ! increase stencil size for Colella & Woodward PPM
-! if (CS%usePPM .and. .not. CS%useHuynh) stencil = 3
-  if (CS%usePPM) stencil = 3
+  use_PPM_stencil = CS%usePPM .and. .not. CS%useHuynhStencilBug
+  if (use_PPM_stencil) stencil = 3
 
   ntr = Reg%ntr
   Idt = 1.0 / dt
@@ -1130,6 +1133,16 @@ subroutine tracer_advect_init(Time, G, US, param_file, diag, CS)
            "Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg))
   end select
 
+  if (CS%useHuynh) then
+    call get_param(param_file, mdl, "USE_HUYNH_STENCIL_BUG", &
+        CS%useHuynhStencilBug, &
+        desc="If true, use a stencil width of 2 in PPM:H3 tracer advection. " &
+        // "This is incorrect and will produce regressions in certain " &
+        // "configurations, but may be required to reproduce results in " &
+        // "legacy simulations.", &
+        default=.false.)
+  endif
+
   id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=CLOCK_MODULE)
   id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=CLOCK_ROUTINE)
   id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=CLOCK_ROUTINE)