*Updates in FPMix and Stokes Most

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -298,7 +298,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., &
                                press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, &
                                cfc=CS%use_CFC, hevap=CS%enthalpy_cpl, tau_mag=.true.)
-    !call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed)
+    call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed)
 
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
@@ -704,7 +704,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
     call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed)
     call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed)
-    !call safe_alloc_ptr(forces%omega_w2x,isd,ied,jsd,jed)
+    call safe_alloc_ptr(forces%omega_w2x,isd,ied,jsd,jed)
 
     if (CS%rigid_sea_ice) then
       call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed)
@@ -865,7 +865,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
       forces%tau_mag(i,j) = gustiness + G%mask2dT(i,j) * sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)
       forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * &
                                sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
-      !forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j))
+      forces%omega_w2x(i,j) = atan(tauy_at_h(i,j), taux_at_h(i,j))
     enddo ; enddo
     call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1)
   else ! C-grid wind stresses.

src/core/MOM.F90

@@ -77,6 +77,7 @@ module MOM
 use MOM_dynamics_split_RK2,    only : step_MOM_dyn_split_RK2, register_restarts_dyn_split_RK2
 use MOM_dynamics_split_RK2,    only : initialize_dyn_split_RK2, end_dyn_split_RK2
 use MOM_dynamics_split_RK2,    only : MOM_dyn_split_RK2_CS, remap_dyn_split_rk2_aux_vars
+use MOM_dynamics_split_RK2,    only : init_dyn_split_RK2_diabatic
 use MOM_dynamics_unsplit_RK2,  only : step_MOM_dyn_unsplit_RK2, register_restarts_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : MOM_dyn_unsplit_RK2_CS
@@ -2069,6 +2070,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   logical :: symmetric         ! If true, use symmetric memory allocation.
   logical :: save_IC           ! If true, save the initial conditions.
   logical :: do_unit_tests     ! If true, call unit tests.
+  logical :: fpmix             ! Needed to decide if BLD should be passed to RK2.
   logical :: test_grid_copy = .false.
 
   logical :: bulkmixedlayer    ! If true, a refined bulk mixed layer scheme is used
@@ -2167,6 +2169,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
                  default=.false.)
   endif
 
+  ! FPMIX is needed to decide if boundary layer depth should be passed to RK2
+  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 .and. .not. CS%split)  then
+    call MOM_error(FATAL, "initialize_MOM: "//&
+       "FPMIX=True only works when SPLIT=True.")
+  endif
+
   call get_param(param_file, "MOM", "BOUSSINESQ", Boussinesq, &
                  "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.)
   call get_param(param_file, "MOM", "SEMI_BOUSSINESQ", semi_Boussinesq, &
@@ -3241,6 +3253,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
                               CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%int_tide_CSp)
   endif
 
+  ! GMM
+  if (CS%split .and. fpmix) then
+    call init_dyn_split_RK2_diabatic(CS%diabatic_CSp, CS%dyn_split_RK2_CSp)
+  endif
+
   if (associated(CS%sponge_CSp)) &
     call init_sponge_diags(Time, G, GV, US, diag, CS%sponge_CSp)
 
@@ -3533,6 +3550,7 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp)
   ! hML is needed when using the ice shelf module
   call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., &
                  do_not_log=.true.)
+
   if (use_ice_shelf .and. associated(CS%Hml)) then
     call register_restart_field(CS%Hml, "hML", .false., restart_CSp, &
                                 "Mixed layer thickness", "m", conversion=US%Z_to_m)

src/core/MOM_dynamics_split_RK2.F90

@@ -180,11 +180,10 @@ module MOM_dynamics_split_RK2
                      !! Euler (1) [nondim].  0 is often used.
   logical :: debug   !< If true, write verbose checksums for debugging purposes.
   logical :: debug_OBC !< If true, do debugging calls for open boundary conditions.
-  logical :: fpmix = .false.                 !< If true, applies profiles of momentum flux magnitude and direction.
+  logical :: fpmix     !< If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.
   logical :: module_is_initialized = .false. !< Record whether this module has been initialized.
 
   !>@{ Diagnostic IDs
-  integer :: id_uold   = -1, id_vold   = -1
   integer :: id_uh     = -1, id_vh     = -1
   integer :: id_umo    = -1, id_vmo    = -1
   integer :: id_umo_2d = -1, id_vmo_2d = -1
@@ -271,6 +270,7 @@ module MOM_dynamics_split_RK2
 public register_restarts_dyn_split_RK2
 public initialize_dyn_split_RK2
 public remap_dyn_split_RK2_aux_vars
+public init_dyn_split_RK2_diabatic
 public end_dyn_split_RK2
 
 !>@{ CPU time clock IDs
@@ -397,7 +397,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   logical :: Use_Stokes_PGF ! If true, add Stokes PGF to hydrostatic PGF
   !---For group halo pass
   logical :: showCallTree, sym
-
+  logical :: lFPpost        ! post diagnostics from vertFPmix when fpmix=true
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
   integer :: cont_stencil, obc_stencil
 
@@ -712,17 +712,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
   call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt_pred, G, GV, US, CS%vertvisc_CSp, &
                      CS%OBC, VarMix)
-  call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, &
-                GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
 
   if (CS%fpmix) then
     hbl(:,:) = 0.0
     if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H)
     if (ASSOCIATED(CS%energetic_PBL_CSp)) &
       call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H)
-    call vertFPmix(up, vp, uold, vold, hbl, h, forces, &
-                   dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC)
-    call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, &
+
+    ! lFPpost must be false in the predictor step to avoid averaging into the diagnostics
+    lFPpost = .false.
+    call vertFPmix(up, vp, uold, vold, hbl, h, forces, dt_pred, lFPpost,  &
+                    G, GV, US, CS%vertvisc_CSp, CS%OBC, waves=waves)
+    call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, &
+                    GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, fpmix=CS%fpmix, waves=waves)
+  else
+    call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, &
                 GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
   endif
 
@@ -964,12 +968,15 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
 
   call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
   call vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
-  call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
-                CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves)
 
   if (CS%fpmix) then
-    call vertFPmix(u, v, uold, vold, hbl, h, forces, dt, &
-                   G, GV, US, CS%vertvisc_CSp, CS%OBC)
+    lFPpost = .true.
+    call vertFPmix(u, v, uold, vold, hbl, h, forces, dt, lFPpost, &
+                   G, GV, US, CS%vertvisc_CSp, CS%OBC, Waves=Waves)
+    call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
+         CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, fpmix=CS%fpmix, waves=waves)
+
+  else
     call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
                   CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
   endif
@@ -1054,11 +1061,6 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     CS%CAu_pred_stored = .false.
   endif
 
-  if (CS%fpmix) then
-    if (CS%id_uold > 0) call post_data(CS%id_uold, uold, CS%diag)
-    if (CS%id_vold > 0) call post_data(CS%id_vold, vold, CS%diag)
-  endif
-
   ! The time-averaged free surface height has already been set by the last call to btstep.
 
   ! Deallocate this memory to avoid a memory leak. ### We should revisit how this array is declared. -RWH
@@ -1292,6 +1294,16 @@ subroutine remap_dyn_split_RK2_aux_vars(G, GV, CS, h_old_u, h_old_v, h_new_u, h_
 
 end subroutine remap_dyn_split_RK2_aux_vars
 
+!> Initializes aspects of the dyn_split_RK2 that depend on diabatic processes.
+subroutine init_dyn_split_RK2_diabatic(diabatic_CSp, CS)
+  type(diabatic_CS),                intent(in) :: diabatic_CSp !< diabatic structure
+  type(MOM_dyn_split_RK2_CS),       pointer    :: CS           !< module control structure
+
+  call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp)
+  call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp)
+
+end subroutine init_dyn_split_RK2_diabatic
+
 !> This subroutine initializes all of the variables that are used by this
 !! dynamic core, including diagnostics and the cpu clocks.
 subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, &
@@ -1403,8 +1415,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
                  "timestep for use in the predictor step of the next split RK2 timestep.", &
                  default=.true.)
   call get_param(param_file, mdl, "FPMIX", CS%fpmix, &
-                 "If true, apply profiles of momentum flux magnitude and "//&
-                 " direction", default=.false.)
+                 "If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", &
+                 default=.false.)
   call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", CS%remap_aux, &
                  "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&
                  "variables that are needed to reproduce across restarts, similarly to "//&
@@ -1483,7 +1495,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
                           CS%SAL_CSp, CS%tides_CSp)
   call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp)
   call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, &
-                     ntrunc, CS%vertvisc_CSp)
+                     ntrunc, CS%vertvisc_CSp, CS%fpmix)
   CS%set_visc_CSp => set_visc
   call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, &
                                 activate=is_new_run(restart_CS) )

src/core/MOM_forcing_type.F90

@@ -74,7 +74,7 @@ module MOM_forcing_type
 
   ! surface stress components and turbulent velocity scale
   real, pointer, dimension(:,:) :: &
-    !omega_w2x     => NULL(), & !< the counter-clockwise angle of the wind stress with respect
+    omega_w2x     => NULL(), & !< the counter-clockwise angle of the wind stress with respect
     ustar         => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1].
     tau_mag       => NULL(), & !< Magnitude of the wind stress averaged over tracer cells,
                                !! including any contributions from sub-gridscale variability
@@ -239,8 +239,8 @@ module MOM_forcing_type
     tau_mag => NULL(), & !< Magnitude of the wind stress averaged over tracer cells, including any
                        !! contributions from sub-gridscale variability or gustiness [R L Z T-2 ~> Pa]
     ustar => NULL(), & !< surface friction velocity scale [Z T-1 ~> m s-1].
-    net_mass_src => NULL() !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1]
-    !omega_w2x    => NULL()    !< the counter-clockwise angle of the wind stress with respect
+    net_mass_src => NULL(), & !< The net mass source to the ocean [R Z T-1 ~> kg m-2 s-1]
+    omega_w2x    => NULL()    !< the counter-clockwise angle of the wind stress with respect
                               !! to the horizontal abscissa (x-coordinate) at tracer points [rad].
 
   ! applied surface pressure from other component models (e.g., atmos, sea ice, land ice)
@@ -378,7 +378,7 @@ module MOM_forcing_type
   integer :: id_taux  = -1
   integer :: id_tauy  = -1
   integer :: id_ustar = -1
-  !integer :: id_omega_w2x = -1
+  integer :: id_omega_w2x = -1
   integer :: id_tau_mag = -1
   integer :: id_psurf     = -1
   integer :: id_TKE_tidal = -1
@@ -1506,8 +1506,8 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
       'Surface friction velocity = [(gustiness + tau_magnitude)/rho0]^(1/2)', &
       'm s-1', conversion=US%Z_to_m*US%s_to_T)
 
-  !handles%id_omega_w2x = register_diag_field('ocean_model', 'omega_w2x', diag%axesT1, Time, &
-  !    'Counter-clockwise angle of the wind stress from the horizontal axis.', 'rad')
+  handles%id_omega_w2x = register_diag_field('ocean_model', 'omega_w2x', diag%axesT1, Time, &
+      'Counter-clockwise angle of the wind stress from the horizontal axis.', 'rad')
 
   if (present(use_berg_fluxes)) then
     if (use_berg_fluxes) then
@@ -2370,11 +2370,11 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres)
       fluxes%ustar(i,j) = forces%ustar(i,j)
     enddo ; enddo
   endif
-  !if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then
-  !  do j=js,je ; do i=is,ie
-  !    fluxes%omega_w2x(i,j) = forces%omega_w2x(i,j)
-  !  enddo ; enddo
-  !endif
+  if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then
+    do j=js,je ; do i=is,ie
+      fluxes%omega_w2x(i,j) = forces%omega_w2x(i,j)
+    enddo ; enddo
+  endif
   if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then
     do j=js,je ; do i=is,ie
       fluxes%tau_mag(i,j) = forces%tau_mag(i,j)
@@ -2516,11 +2516,11 @@ subroutine copy_back_forcing_fields(fluxes, forces, G)
       forces%ustar(i,j) = fluxes%ustar(i,j)
     enddo ; enddo
   endif
-  !if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then
-  !  do j=js,je ; do i=is,ie
-  !    forces%omega_w2x(i,j) = fluxes%omega_w2x(i,j)
-  !  enddo ; enddo
-  !endif
+  if (associated(forces%omega_w2x) .and. associated(fluxes%omega_w2x)) then
+    do j=js,je ; do i=is,ie
+      forces%omega_w2x(i,j) = fluxes%omega_w2x(i,j)
+    enddo ; enddo
+  endif
   if (associated(forces%tau_mag) .and. associated(fluxes%tau_mag)) then
     do j=js,je ; do i=is,ie
       forces%tau_mag(i,j) = fluxes%tau_mag(i,j)
@@ -3172,8 +3172,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h
     if ((handles%id_ustar > 0) .and. associated(fluxes%ustar)) &
       call post_data(handles%id_ustar, fluxes%ustar, diag)
 
-    !if ((handles%id_omega_w2x > 0) .and. associated(fluxes%omega_w2x)) &
-    !  call post_data(handles%id_omega_w2x, fluxes%omega_w2x, diag)
+    if ((handles%id_omega_w2x > 0) .and. associated(fluxes%omega_w2x)) &
+      call post_data(handles%id_omega_w2x, fluxes%omega_w2x, diag)
 
     if ((handles%id_ustar_berg > 0) .and. associated(fluxes%ustar_berg)) &
       call post_data(handles%id_ustar_berg, fluxes%ustar_berg, diag)
@@ -3512,7 +3512,7 @@ end subroutine myAlloc
 subroutine deallocate_forcing_type(fluxes)
   type(forcing), intent(inout) :: fluxes !< Forcing fields structure
 
-  !if (associated(fluxes%omega_w2x))            deallocate(fluxes%omega_w2x)
+  if (associated(fluxes%omega_w2x))            deallocate(fluxes%omega_w2x)
   if (associated(fluxes%ustar))                deallocate(fluxes%ustar)
   if (associated(fluxes%ustar_gustless))       deallocate(fluxes%ustar_gustless)
   if (associated(fluxes%tau_mag))              deallocate(fluxes%tau_mag)
@@ -3572,7 +3572,7 @@ end subroutine deallocate_forcing_type
 subroutine deallocate_mech_forcing(forces)
   type(mech_forcing), intent(inout) :: forces  !< Forcing fields structure
 
-  !if (associated(forces%omega_w2x))      deallocate(forces%omega_w2x)
+  if (associated(forces%omega_w2x))      deallocate(forces%omega_w2x)
   if (associated(forces%taux))           deallocate(forces%taux)
   if (associated(forces%tauy))           deallocate(forces%tauy)
   if (associated(forces%ustar))          deallocate(forces%ustar)