Merge tag 'dev/ncar_241101' into running_means

config_src/external/MARBL/marbl_interface.F90

@@ -36,7 +36,8 @@ module marbl_interface
        contains
         procedure, public  :: put_setting                !< dummy put_setting routine
         procedure, public  :: get_setting                !< dummy get_setting routine
-        procedure, public  :: init                       !< dummy routine
+        procedure, public  :: init                       !< dummy init routine
+        procedure, public  :: compute_totChl             !< dummy routine to compute total Chlorophyll
         procedure, public  :: surface_flux_compute       !< dummy surface flux routine
         procedure, public  :: interior_tendency_compute  !< dummy interior tendency routine
         procedure, public  :: add_output_for_GCM         !< dummy add_output_for_GCM routine
@@ -89,6 +90,15 @@ subroutine init(self,                  &
         call MOM_error(FATAL, error_msg)
     end subroutine init
 
+    !> Dummy version of MARBL's compute_totChl() function
+    subroutine compute_totChl(self)
+
+      class(marbl_interface_class), intent(inout) :: self
+
+      call MOM_error(FATAL, error_msg)
+
+    end subroutine compute_totChl
+
     !> Dummy version of MARBL's surface_flux_compute() function
     subroutine surface_flux_compute(self)
 

config_src/external/MARBL/marbl_interface_public_types.F90

@@ -77,6 +77,7 @@ module marbl_interface_public_types
         ! that needs to be passed to the GCM / flux coupler.
         ! Data must be accessed via the marbl_output_for_GCM_type
         ! data structure.
+        character(len=0) :: short_name  !< dummy name
         real, allocatable, dimension(:)   :: forcing_field_0d  !< dummy forcing_field_0d
         real, allocatable, dimension(:,:) :: forcing_field_1d  !< forcing_field_1d
     end type marbl_single_output_type

src/ALE/MOM_ALE.F90

@@ -380,12 +380,12 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
       'Rate of change in half rho0 times depth integral of squared zonal'//&
       ' velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&
       ' this measures the change before the KE-conserving correction is applied.', &
-      'W m-2', conversion=US%RZ3_T3_to_W_m2 * US%L_to_Z**2)
+      'W m-2', conversion=GV%H_to_kg_m2 * US%L_T_to_m_s**2 * US%s_to_T)
   CS%id_remap_delta_integ_v2 = register_diag_field('ocean_model', 'ale_v2', diag%axesCv1, Time, &
       'Rate of change in half rho0 times depth integral of squared meridional'//&
       ' velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&
       ' this measures the change before the KE-conserving correction is applied.', &
-      'W m-2', conversion=US%RZ3_T3_to_W_m2 * US%L_to_Z**2)
+      'W m-2', conversion=GV%H_to_kg_m2 * US%L_T_to_m_s**2 * US%s_to_T)
 
 end subroutine ALE_register_diags
 
@@ -1172,7 +1172,11 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
         ke_c_tgt = ke_c_tgt + h2(k) * (u_tgt(k) - u_bt)**2
       enddo
       ! Next rescale baroclinic component on target grid to conserve ke
-      rescale_coef = min(1.25, sqrt(ke_c_src / (ke_c_tgt + 1.E-19)))
+      if (ke_c_src < 1.5625 * ke_c_tgt) then
+        rescale_coef = sqrt(ke_c_src / ke_c_tgt)
+      else
+        rescale_coef = 1.25
+      endif
       do k=1,nz
         u_tgt(k) = u_bt + rescale_coef * (u_tgt(k) - u_bt)
       enddo
@@ -1182,7 +1186,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
       do k=1,nz
         u2h_tot = u2h_tot + h2(k) * (u_tgt(k)**2)
       enddo
-      du2h_tot(I,j) = GV%H_to_RZ * u2h_tot * I_dt
+      du2h_tot(I,j) = u2h_tot * I_dt
     endif
 
     if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) &
@@ -1240,7 +1244,11 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
         ke_c_tgt = ke_c_tgt + h2(k) * (v_tgt(k) - v_bt)**2
       enddo
       ! Next rescale baroclinic component on target grid to conserve ke
-      rescale_coef = min(1.25, sqrt(ke_c_src / (ke_c_tgt + 1.E-19)))
+      if (ke_c_src < 1.5625 * ke_c_tgt) then
+        rescale_coef = sqrt(ke_c_src / ke_c_tgt)
+      else
+        rescale_coef = 1.25
+      endif
       do k=1,nz
         v_tgt(k) = v_bt + rescale_coef * (v_tgt(k) - v_bt)
       enddo
@@ -1250,7 +1258,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
       do k=1,nz
         v2h_tot = v2h_tot + h2(k) * (v_tgt(k)**2)
       enddo
-      dv2h_tot(I,j) = GV%H_to_RZ * v2h_tot * I_dt
+      dv2h_tot(I,j) = v2h_tot * I_dt
     endif
 
     if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then

src/core/MOM.F90

@@ -2895,7 +2895,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   endif
 
   if (.not. CS%adiabatic) then
-    call register_diabatic_restarts(G, US, param_file, CS%int_tide_CSp, restart_CSp)
+    call register_diabatic_restarts(G, US, param_file, CS%int_tide_CSp, restart_CSp, CS%diabatic_CSp)
   endif
 
   call callTree_waypoint("restart registration complete (initialize_MOM)")

src/parameterizations/vertical/MOM_CVMix_KPP.F90

@@ -13,6 +13,7 @@ module MOM_CVMix_KPP
 use MOM_file_parser,    only : openParameterBlock, closeParameterBlock
 use MOM_grid,           only : ocean_grid_type, isPointInCell
 use MOM_interface_heights, only : thickness_to_dz
+use MOM_restart,        only : MOM_restart_CS, register_restart_field
 use MOM_unit_scaling,   only : unit_scale_type
 use MOM_variables,      only : thermo_var_ptrs
 use MOM_verticalGrid,   only : verticalGrid_type
@@ -36,6 +37,7 @@ module MOM_CVMix_KPP
 
 #include "MOM_memory.h"
 
+public :: register_KPP_restarts
 public :: KPP_init
 public :: KPP_compute_BLD
 public :: KPP_calculate
@@ -112,6 +114,7 @@ module MOM_CVMix_KPP
   logical :: LT_K_Enhancement          !< Flags if enhancing mixing coefficients due to LT
   integer :: LT_K_Shape                !< Integer for constant or shape function enhancement
   integer :: LT_K_Method               !< Integer for mixing coefficients LT method
+  real    :: KPP_CVt2                  !< Parameter for Stokes MOST convection entrainment
   real    :: KPP_K_ENH_FAC             !< Factor to multiply by K if Method is CONSTANT [nondim]
   logical :: LT_Vt2_Enhancement        !< Flags if enhancing Vt2 due to LT
   integer :: LT_VT2_METHOD             !< Integer for Vt2 LT method
@@ -152,7 +155,7 @@ module MOM_CVMix_KPP
   !>@}
 
   ! Diagnostics arrays
-  real, allocatable, dimension(:,:)   :: OBLdepth  !< Depth (positive) of ocean boundary layer (OBL) [Z ~> m]
+  real, pointer,     dimension(:,:)   :: OBLdepth  !< Depth (positive) of ocean boundary layer (OBL) [Z ~> m]
   real, allocatable, dimension(:,:)   :: OBLdepth_original  !< Depth (positive) of OBL [Z ~> m] without smoothing
   real, allocatable, dimension(:,:)   :: StokesParXI !< Stokes similarity parameter
   real, allocatable, dimension(:,:)   :: Lam2        !< La^(-2) = Ustk0/u*
@@ -188,6 +191,33 @@ module MOM_CVMix_KPP
 
 contains
 
+!> Routine to register restarts, pass-through to children modules
+subroutine register_KPP_restarts(G, param_file, restart_CSp, CS)
+  type(ocean_grid_type), intent(in)    :: G           !< The ocean's grid structure
+  type(param_file_type), intent(in)    :: param_file  !< A structure to parse for run-time parameters
+  type(MOM_restart_CS),  pointer       :: restart_CSp  !< MOM restart control structure
+  type(KPP_CS),         pointer        :: CS           !< module control structure
+
+  character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
+  logical :: use_kpp, fpmix
+
+  if (associated(CS)) call MOM_error(FATAL, 'MOM_CVMix_KPP, register_KPP_restarts: '// &
+           'Control structure has already been initialized')
+  call get_param(param_file, mdl, "USE_KPP", use_kpp, default=.false., do_not_log=.true.)
+  ! Forego remainder of initialization if not using this scheme
+  if (.not. use_kpp) return
+  allocate(CS)
+
+  allocate(CS%OBLdepth(SZI_(G),SZJ_(G)), source=0.0)
+
+  ! FPMIX is needed to decide if boundary layer depth should be added to restart file
+  call get_param(param_file, '', "FPMIX", fpmix, &
+                 "If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", &
+                 default=.false., do_not_log=.true.)
+    if (fpmix) call register_restart_field(CS%OBLdepth, 'KPP_OBLdepth', .false., restart_CSp)
+
+end subroutine register_KPP_restarts
+
 !> Initialize the CVMix KPP module and set up diagnostics
 !! Returns True if KPP is to be used, False otherwise.
 logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
@@ -213,9 +243,6 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
   logical :: CS_IS_ONE=.false.         !< Logical for setting Cs based on Non-local
   logical :: lnoDGat1=.false.          !< True => G'(1) = 0 (shape function)
                                        !! False => compute G'(1) as in LMD94
-  if (associated(CS)) call MOM_error(FATAL, 'MOM_CVMix_KPP, KPP_init: '// &
-           'Control structure has already been initialized')
-
   ! Read parameters
   call get_param(paramFile, mdl, "USE_KPP", KPP_init, default=.false., do_not_log=.true.)
   call log_version(paramFile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // &
@@ -226,7 +253,6 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
                  default=.false.)
   ! Forego remainder of initialization if not using this scheme
   if (.not. KPP_init) return
-  allocate(CS)
 
   call get_param(paramFile, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
@@ -489,6 +515,10 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
                    units="m", default=1.0, scale=US%m_to_Z)
   endif
 
+  call get_param(paramFile, mdl, "KPP_CVt2", CS%KPP_CVt2, &
+                 'Parameter for Stokes MOST convection entrainment', &
+                 units="nondim", default=1.6)
+
   call get_param(paramFile, mdl, "ANSWER_DATE", CS%answer_date, &
                  "The vintage of the order of arithmetic in the CVMix KPP calculations.  Values "//&
                  "below 20240501 recover the answers from early in 2024, while higher values "//&
@@ -504,6 +534,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
                        minVtsqr=US%L_T_to_m_s**2*CS%minVtsqr, &
                        vonKarman=CS%vonKarman,             &
                        surf_layer_ext=CS%surf_layer_ext,   &
+                       CVt2=CS%KPP_CVt2,                   &
                        interp_type=CS%interpType,          &
                        interp_type2=CS%interpType2,        &
                        lEkman=CS%computeEkman,             &
@@ -599,7 +630,6 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
       'Surface-layer Langmuir number computed in [CVMix] KPP','nondim')
 
   allocate( CS%N( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
-  allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ), source=0. )
   allocate( CS%StokesParXI( SZI_(G), SZJ_(G) ), source=0. )
   allocate( CS%Lam2    ( SZI_(G), SZJ_(G) ), source=0. )
   allocate( CS%kOBL( SZI_(G), SZJ_(G) ), source=0. )
@@ -1028,6 +1058,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
   real :: hcorr         ! A cumulative correction arising from inflation of vanished layers [Z ~> m]
 
   ! For Langmuir Calculations
+  real :: Vt_layer     ! non-dimensional extent contribution to unresolved shear
   real :: LangEnhVt2   ! Langmuir enhancement for unresolved shear [nondim]
   real, dimension(GV%ke) :: U_H, V_H ! Velocities at tracer points [L T-1 ~> m s-1]
   real :: MLD_guess    ! A guess at the mixed layer depth for calculating the Langmuir number [Z ~> m]
@@ -1085,7 +1116,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
   !$OMP                           uS_SLD, vS_SLD, uS_SLC, vS_SLC, uSbar_SLD, vSbar_SLD, &
   !$OMP                           StokesXI, StokesXI_1d, StokesVt_1d, kbl) &
   !$OMP                           shared(G, GV, CS, US, uStar, h, dz, buoy_scale, buoyFlux, &
-  !$OMP                           Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult)
+  !$OMP                           Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult, Vt_layer)
   do j = G%jsc, G%jec
     do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then
 
@@ -1099,7 +1130,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
           vE_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)-Waves%US_y(i,J,k)-Waves%US_y(i,J-1,k))
         enddo
       endif
-
       ! things independent of position within the column
       Coriolis = 0.25*US%s_to_T*( (G%CoriolisBu(i,j)   + G%CoriolisBu(i-1,j-1)) + &
                                   (G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) )
@@ -1137,9 +1167,11 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
           endif
         enddo
 
-         if (CS%StokesMOST) then
+        if (CS%StokesMOST) then
+          ! if k=1, want buoyFlux(i,j,1) - buoyFlux(i,j,2), otherwise
+          ! subtract average of buoyFlux(i,j,k) and buoyFlux(i,j,k+1)
           surfBuoyFlux   = buoy_scale * &
-                          (buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,k)+buoyFlux(i,j,k+1)) )
+                          (buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,max(2,k))+buoyFlux(i,j,k+1)) )
           surfBuoyFlux2(k) = surfBuoyFlux
           call Compute_StokesDrift(i,j, iFaceHeight(k),iFaceHeight(k+1), &
               uS_Hi(k+1), vS_Hi(k+1), uS_H(k), vS_H(k), uSbar_H(k), vSbar_H(k), Waves)
@@ -1149,10 +1181,10 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
                surfFricVel,waves%omega_w2x(i,j), uE_H, vE_H, uS_Hi, vS_Hi, uSbar_H, vSbar_H, uS_SLD,&
                vS_SLD, uSbar_SLD, vSbar_SLD, StokesXI, CVMix_kpp_params_user=CS%KPP_params )
           StokesXI_1d(k) = StokesXI
-          StokesVt_1d(k) = StokesXI
+          StokesVt_1d(k) = 0.0  ! StokesXI
 
           ! average temperature, salinity, u and v over surface layer starting at ksfc
-          delH      = SLdepth_0d + iFaceHeight(ksfc-1)
+          delH      = SLdepth_0d + iFaceHeight(ksfc)
           surfHtemp = Temp(i,j,ksfc) * delH
           surfHsalt = Salt(i,j,ksfc) * delH
           surfHu    = (uE_H(ksfc) + uSbar_SLD)  * delH
@@ -1171,12 +1203,12 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
           surfSalt = surfHsalt * I_hTot
           surfU    = surfHu    * I_hTot
           surfV    = surfHv    * I_hTot
+
           Uk       = uE_H(k) + uS_H(k) - surfU
           Vk       = vE_H(k) + vS_H(k) - surfV
 
         else   !not StokesMOST
           StokesXI_1d(k) = 0.0
-
           ! average temperature, salinity, u and v over surface layer
           ! use C-grid average to get u and v on T-points.
           surfHtemp = 0.0
@@ -1205,13 +1237,20 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
             endif
 
           enddo
+          !I_hTot = 1./hTot
+          !surfTemp = surfHtemp * I_hTot
+          !surfSalt = surfHsalt * I_hTot
+          !surfU    = surfHu    * I_hTot
+          !surfV    = surfHv    * I_hTot
+          !surfUs   = surfHus   * I_hTot
+          !surfVs   = surfHvs   * I_hTot
+
           surfTemp = surfHtemp / hTot
           surfSalt = surfHsalt / hTot
           surfU    = surfHu    / hTot
           surfV    = surfHv    / hTot
           surfUs   = surfHus   / hTot
           surfVs   = surfHvs   / hTot
-
           ! vertical shear between present layer and surface layer averaged surfU and surfV.
           ! C-grid average to get Uk and Vk on T-points.
           Uk         = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU
@@ -1296,12 +1335,13 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
         z_inter(K) = US%Z_to_m*iFaceHeight(K)
       enddo
 
-      ! turbulent velocity scales w_s and w_m computed at the cell centers.
-      ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales
-      ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass
-      ! sigma=CS%surf_layer_ext for this calculation.
+      ! CVMix_kpp_compute_turbulent_scales_1d_OBL computes w_s velocity scale at cell centers for
+      ! CVmix_kpp_compute_bulk_Richardson call to CVmix_kpp_compute_unresolved_shear
+      ! at sigma=Vt_layer (CS%surf_layer_ext or 1.0) for this calculation.
+      ! StokesVt_1d controls Stokes enhancement  (= 0 for none)
+      Vt_layer = 1.0      ! CS%surf_layer_ext
       call CVMix_kpp_compute_turbulent_scales( &                             ! 1d_OBL
-              CS%surf_layer_ext, & ! (in)  Normalized surface layer depth; sigma = CS%surf_layer_ext
+              Vt_layer,          & ! (in)  Boundary layer extent contributing to unresolved shear
               OBL_depth,         & ! (in)  OBL depth [m]
               surfBuoyFlux2,     & ! (in)  Buoyancy flux at surface [m2 s-3]
               surfFricVel,       & ! (in)  Turbulent friction velocity at surface [m s-1]
@@ -1341,7 +1381,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
                   N_iface=N_col,                     & ! Buoyancy frequency [s-1]
                   EFactor=LangEnhVT2,                & ! Langmuir enhancement factor [nondim]
                   LaSL=CS%La_SL(i,j),                & ! surface layer averaged Langmuir number [nondim]
-                  bfsfc=surfBuoyFlux2,                & ! surface buoyancy flux [m2 s-3]
+                  bfsfc=surfBuoyFlux2,               & ! surface buoyancy flux [m2 s-3]
                   uStar=surfFricVel,                 & ! surface friction velocity [m s-1]
                   CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters
 
@@ -1368,7 +1408,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
       CS%OBLdepth(i,j) = US%m_to_Z * KPP_OBL_depth
 
     if (CS%StokesMOST) then
-      kbl = nint(CS%kOBL(i,j))
+      kbl = int(CS%kOBL(i,j))
       SLdepth_0d = CS%surf_layer_ext*CS%OBLdepth(i,j)
       surfBuoyFlux = surfBuoyFlux2(kbl)
         ! find ksfc for cell where "surface layer" sits