diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index ca0dfe0145..0d06e0e33a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2354,7 +2354,9 @@ 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 :: nlVstress ! Needed to decide if BLD should be passed to RK2. + logical :: fpmix ! The value of the deprecated FPMIX runtime parameter, now used + ! to determine the default for its replacement, NL_VSTRESS. logical :: test_grid_copy = .false. logical :: bulkmixedlayer ! If true, a refined bulk mixed layer scheme is used @@ -2464,13 +2466,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 + ! NL_VSTRESS is needed to decide if boundary layer depth should be passed to RK2. + ! The deprecated FPMIX determines its default value, in case NL_VSTRESS is not set. 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 get_param(param_file, '', "NL_VSTRESS", nlVstress, & + "If true, add non-local momentum flux increments.", & + default=fpmix, do_not_log=.true.) + if (nlVstress .and. .not. CS%split) then call MOM_error(FATAL, "initialize_MOM: "//& - "FPMIX=True only works when SPLIT=True.") + "NL_VSTRESS=True only works when SPLIT=True.") endif call openParameterBlock(param_file, 'KPP', do_not_log=.true.) call get_param(param_file, '', 'STOKES_MOST', CS%StokesMOST, & @@ -3677,7 +3682,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%vertex_shear = kappa_shear_at_vertex(param_file) ! GMM, the following is needed to get BLDs into the dynamics module - if (CS%split .and. fpmix) then + if (CS%split .and. nlVstress) then call init_dyn_split_RK2_diabatic(CS%diabatic_CSp, CS%dyn_split_RK2_CSp) endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 8588a0c41e..7fdd76e720 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -29,7 +29,7 @@ module MOM_dynamics_split_RK2 use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint -use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : get_param, read_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_io, only : vardesc, var_desc, EAST_FACE, NORTH_FACE use MOM_restart, only : register_restart_field, register_restart_pair @@ -76,7 +76,7 @@ module MOM_dynamics_split_RK2 use MOM_unit_scaling, only : unit_scale_type use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant use MOM_vert_friction, only : vertvisc_init, vertvisc_end, vertvisc_CS -use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix +use MOM_vert_friction, only : updateCFLtruncationValue, vertNLstress use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units use MOM_wave_interface, only : wave_parameters_CS, Stokes_PGF @@ -191,7 +191,7 @@ module MOM_dynamics_split_RK2 logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debug_OBC !< If true, do additional calls resetting values to help debug the correctness !! of the open boundary condition code. - logical :: fpmix !< If true, add non-local momentum flux increments and diffuse down the Eulerian gradient. + logical :: nlVstress !< If true, add non-local momentum flux increments. logical :: module_is_initialized = .false. !< Record whether this module has been initialized. logical :: visc_rem_dt_bug = .true. !< If true, recover a bug that uses dt_pred rather than dt for vertvisc_rem !! at the end of predictor. @@ -375,11 +375,6 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_OBC ! The starting meridional velocities, which are ! saved for use in the radiation open boundary condition code [L T-1 ~> m s-1] - ! GMM, TODO: make these allocatable? - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for fpmix - ! [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold ! v-velocity before vert_visc is applied, for fpmix - ! [L T-1 ~> m s-1] real :: pres_to_eta ! A factor that converts pressures to the units of eta ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1] real, pointer, dimension(:,:) :: & @@ -414,8 +409,9 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f logical :: Use_Stokes_PGF ! If true, add Stokes PGF to hydrostatic PGF !---For group halo pass logical :: showCallTree, sym - logical :: lFPpost ! Used to only post diagnostics in vertFPmix when fpmix=true and + logical :: lNLpost ! Used to only post diagnostics in vertNLstress when nlVstress=true and ! in the corrector step (not the predict) + logical :: waves_ok ! True if the Waves control structure is present and associated integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: cont_stencil, obc_stencil, vel_stencil integer :: cor_stencil @@ -431,6 +427,15 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("step_MOM_dyn_split_RK2(), MOM_dynamics_split_RK2.F90") + ! The nonlocal stress increments require wave (Stokes drift) information, so fail early + ! if an associated Waves control structure was not provided. + if (CS%nlVstress) then + waves_ok = .false. + if (present(Waves)) then ; if (associated(Waves)) waves_ok = .true. ; endif + if (.not. waves_ok) call MOM_error(FATAL, "MOM_dynamics_split_RK2, step_MOM_dyn_split_RK2: "//& + "NL_VSTRESS=True requires an associated Waves control structure (e.g., USE_WAVES=True).") + endif + !$OMP parallel do default(shared) do k=1,nz do j=G%jsd,G%jed ; do i=G%isdB,G%iedB ; up(i,j,k) = 0.0 ; enddo ; enddo @@ -731,43 +736,23 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s) endif - if (CS%fpmix) then - uold(:,:,:) = 0.0 - vold(:,:,:) = 0.0 - do k = 1, nz - do j = js , je - do I = Isq, Ieq - uold(I,j,k) = up(I,j,k) - enddo - enddo - do J = Jsq, Jeq - do i = is, ie - vold(i,J,k) = vp(i,J,k) - enddo - enddo - enddo - endif - 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) - if (CS%fpmix) then + if (CS%nlVstress) 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) - ! 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, CS%Cemp_NL, & - 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) + ! lNLpost must be false in the predictor step to avoid averaging into the diagnostics + lNLpost = .false. + call vertNLstress(up, vp, hbl, h, forces, dt_pred, lNLpost, CS%Cemp_NL, & + G, GV, US, CS%vertvisc_CSp, CS%OBC, waves=waves) endif + 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 (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)") if (G%nonblocking_updates) then @@ -993,37 +978,16 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f ! u_av <- u_av + dt d/dz visc d/dz u_av call cpu_clock_begin(id_clock_vertvisc) - if (CS%fpmix) then - uold(:,:,:) = 0.0 - vold(:,:,:) = 0.0 - do k = 1, nz - do j = js , je - do I = Isq, Ieq - uold(I,j,k) = u_inst(I,j,k) - enddo - enddo - do J = Jsq, Jeq - do i = is, ie - vold(i,J,k) = v_inst(i,J,k) - enddo - enddo - enddo - endif - call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix) - if (CS%fpmix) then - lFPpost = .true. - call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, lFPpost, CS%Cemp_NL, & - G, GV, US, CS%vertvisc_CSp, CS%OBC, Waves=Waves) - call vertvisc(u_inst, v_inst, 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_inst, v_inst, 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%nlVstress) then + lNLpost = .true. + call vertNLstress(u_inst, v_inst, hbl, h, forces, dt, lNLpost, CS%Cemp_NL, & + G, GV, US, CS%vertvisc_CSp, CS%OBC, Waves=Waves) endif + call vertvisc(u_inst, v_inst, 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 (G%nonblocking_updates) then call cpu_clock_end(id_clock_vertvisc) @@ -1424,6 +1388,9 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to ! recreate the bugs, or if false bugs are only used if actively selected. logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG. + logical :: fpmix ! The value of the FPMIX runtime parameter, which is + ! adopted as the default for its replacement, NL_VSTRESS. + logical :: fpmix_set ! True if FPMIX was explicitly set in the input parameter files. integer :: cor_stencil integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz @@ -1489,10 +1456,17 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p "If true, calculate the Coriolis accelerations at the end of each "//& "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, add non-local momentum flux increments and diffuse down the Eulerian gradient.", & - default=.false.) - if (CS%fpmix) then + ! FPMIX has been renamed NL_VSTRESS. To allow a gradual transition, an explicitly set + ! FPMIX still works, acting as the default for NL_VSTRESS, but it triggers a warning. + fpmix = .false. ; fpmix_set = .false. + call read_param(param_file, "FPMIX", fpmix, set=fpmix_set) + if (fpmix_set) call MOM_error(WARNING, "MOM_dynamics_split_RK2, initialize_dyn_split_RK2: "//& + "The FPMIX runtime parameter has been renamed NL_VSTRESS. FPMIX still works, but will be "//& + "obsoleted in future. Please use NL_VSTRESS instead.") + call get_param(param_file, mdl, "NL_VSTRESS", CS%nlVstress, & + "If true, add non-local momentum flux increments.", & + default=fpmix) + if (CS%nlVstress) then call get_param(param_file, "MOM", "CEMP_NL", CS%Cemp_NL, & "Empirical coefficient of non-local momentum mixing.", & units="nondim", default=3.6) @@ -1597,7 +1571,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p 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, CS%fpmix) + ntrunc, CS%vertvisc_CSp, CS%nlVstress) CS%set_visc_CSp => set_visc call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, activate=is_new_run(restart_CS) ) diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index dcdfd9b834..60aad4711b 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -74,7 +74,7 @@ module MOM_dynamics_split_RK2b use MOM_unit_scaling, only : unit_scale_type use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant use MOM_vert_friction, only : vertvisc_init, vertvisc_end, vertvisc_CS -use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix +use MOM_vert_friction, only : updateCFLtruncationValue, vertNLstress use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units use MOM_wave_interface, only : wave_parameters_CS, Stokes_PGF @@ -181,7 +181,7 @@ module MOM_dynamics_split_RK2b logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debug_OBC !< If true, do additional calls resetting values to help verify the correctness !! of the open boundary condition code. - logical :: fpmix = .false. !< If true, applies profiles of momentum flux magnitude and direction. + logical :: nlVstress = .false. !< If true, add non-local momentum flux increments. logical :: module_is_initialized = .false. !< Record whether this module has been initialized. logical :: visc_rem_dt_bug = .true. !< If true, recover a bug that uses dt_pred rather than dt for vertvisc_rem !! at the end of predictor. @@ -371,9 +371,9 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1] ! GMM, TODO: make these allocatable? - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for fpmix + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for nlVstress ! [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold ! v-velocity before vert_visc is applied, for fpmix + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold ! v-velocity before vert_visc is applied, for nlVstress ! [L T-1 ~> m s-1] real :: pres_to_eta ! A factor that converts pressures to the units of eta ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1] @@ -747,7 +747,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s) endif - ! if (CS%fpmix) then + ! if (CS%nlVstress) then ! uold(:,:,:) = 0.0 ! vold(:,:,:) = 0.0 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -764,10 +764,10 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc 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 + ! if (CS%nlVstress) then ! hbl(:,:) = 0.0 ! if (associated(visc%h_ML)) hbl(:,:) = visc%h_ML(:,:) - ! call vertFPmix(up, vp, uold, vold, hbl, h, forces, & + ! call vertNLstress(up, vp, 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, & ! GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) @@ -980,7 +980,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc ! u_av <- u_av + dt d/dz visc d/dz u_av call cpu_clock_begin(id_clock_vertvisc) - ! if (CS%fpmix) then + ! if (CS%nlVstress) then ! uold(:,:,:) = 0.0 ! vold(:,:,:) = 0.0 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -996,8 +996,8 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc call vertvisc(u_inst, v_inst, 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_inst, v_inst, uold, vold, hbl, h, forces, dt, & + ! if (CS%nlVstress) then + ! call vertNLstress(u_inst, v_inst, hbl, h, forces, dt, & ! G, GV, US, CS%vertvisc_CSp, CS%OBC) ! call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & ! CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) @@ -1063,7 +1063,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, h, CS%OBC) endif - ! if (CS%fpmix) then + ! if (CS%nlVstress) 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 @@ -1388,10 +1388,10 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, "predictor step. This should make little difference in the "//& "deep ocean but appears to help for vanished layers. If false, "//& "uses the same mass source as from the predictor step.", default=.true.) - ! call get_param(param_file, mdl, "FPMIX", CS%fpmix, & - ! "If true, apply profiles of momentum flux magnitude and direction.", & + ! call get_param(param_file, mdl, "NL_VSTRESS", CS%nlVstress, & + ! "If true, add non-local momentum flux increments.", & ! default=.false.) - CS%fpmix = .false. + CS%nlVstress = .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 "//& diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index c72d25a260..192fada40c 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -131,6 +131,8 @@ module MOM_CVMix_KPP real :: KPP_ER_Cb !< Entrainment Rule TKE buoyancy production weight [nondim] real :: KPP_ER_Cs !< Entrainment Rule TKE Stokes production weight [nondim] real :: KPP_ER_Cu !< Entrainment Rule TKE shear production weight [nondim] + logical :: KPP_USE_ER !< If true and STOKES_MOST is enabled, use the Entrainment + !! Rule to diagnose entraining boundary layer depths logical :: STOKES_MIXING !< Flag if model is mixing down Stokes gradient !! This is relevant for which current to use in RiB logical :: OBL_depth_bounds_bug !< If true, limit the KPP boundary layer depth relative to @@ -222,7 +224,7 @@ subroutine register_KPP_restarts(G, param_file, restart_CSp, CS) type(KPP_CS), pointer :: CS !< module control structure character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module - logical :: use_kpp, fpmix + logical :: use_kpp, nlVstress, fpmix if (associated(CS)) call MOM_error(FATAL, 'MOM_CVMix_KPP, register_KPP_restarts: '// & 'Control structure has already been initialized') @@ -233,11 +235,14 @@ subroutine register_KPP_restarts(G, param_file, restart_CSp, 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 + ! NL_VSTRESS is needed to decide if boundary layer depth should be added to restart file. + ! The FPMIX now determines its default value, in case NL_VSTRESS is not set. 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) + call get_param(param_file, '', "NL_VSTRESS", nlVstress, & + "If true, add non-local momentum flux increments.", & + default=fpmix, do_not_log=.true.) + if (nlVstress) call register_restart_field(CS%OBLdepth, 'KPP_OBLdepth', .false., restart_CSp) end subroutine register_KPP_restarts @@ -542,6 +547,10 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 'Parameter for Stokes MOST convection entrainment (unresolved shear)', & units="nondim", default=1.6) + call get_param(paramFile, mdl, "KPP_USE_ER", CS%KPP_USE_ER, & + 'If true and STOKES_MOST is enabled, use the Entrainment Rule to '//& + 'diagnose entraining boundary layer depths.', & + default=.true., do_not_log=.not.CS%StokesMOST) call get_param(paramFile, mdl, "KPP_ER_Cb", CS%KPP_ER_Cb, & 'Entrainment Rule TKE buoyancy production weight', & units="nondim", default=0.96) @@ -1466,7 +1475,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer else ERdepth = 0.0 - if ( CS%StokesMOST .and. (surfBuoy_NS < 0.0) ) then + if ( (CS%StokesMOST .and. CS%KPP_USE_ER) .and. (surfBuoy_NS < 0.0) ) then ! Search for Entrainment rule depth (ER_depth) call CVMix_kpp_compute_ER_depth( & z_inter, & ! (in) Interface heights <= 0 [m] @@ -1486,7 +1495,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl endif endif - ! Original Richardson Number method (always the case with CS%StokesMOST=False) + ! Original Richardson Number method (always the case when the Entrainment Rule is not in use) if (CS%ERdepth(i,j) == 0.) then Vt_layer = 1.0 ! CS%surf_layer_ext call CVMix_kpp_compute_turbulent_scales( & ! 1d_OBL diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 2bfada3bf1..f7ea29fe99 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -2714,7 +2714,8 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C ! Local variables logical :: use_kappa_shear, KS_at_vertex logical :: adiabatic, useKPP, useEPBL, use_ideal_age - logical :: do_brine_plume, use_hor_bnd_diff, use_neutral_diffusion, use_fpmix, use_StokesMOST + logical :: do_brine_plume, use_hor_bnd_diff, use_neutral_diffusion, use_nlVstress, use_StokesMOST + logical :: use_fpmix ! This now determines the default for use_nlVstress. logical :: use_CVMix_shear, MLE_use_PBL_MLD, MLE_use_Bodner, use_CVMix_conv integer :: isd, ied, jsd, jed, nz real :: hfreeze !< If hfreeze > 0 [Z ~> m], melt potential will be computed. @@ -2797,6 +2798,8 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C default=.false., do_not_log=.true.) call get_param(param_file, mdl, "FPMIX", use_fpmix, & default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "NL_VSTRESS", use_nlVstress, & + default=use_fpmix, do_not_log=.true.) call openParameterBlock(param_file, 'KPP', do_not_log=.true.) call get_param(param_file, mdl, 'STOKES_MOST', use_StokesMOST, & default=.false., do_not_log=.true.) @@ -2814,7 +2817,7 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C if (use_StokesMOST .and. MLE_use_Bodner) then call safe_alloc_ptr(visc%Lam2, isd, ied, jsd, jed) endif - if ((hfreeze >= 0.0) .or. MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. & + if ((hfreeze >= 0.0) .or. MLE_use_PBL_MLD .or. do_brine_plume .or. use_nlVstress .or. & use_neutral_diffusion .or. use_hor_bnd_diff .or. use_ideal_age) then call safe_alloc_ptr(visc%h_ML, isd, ied, jsd, jed) endif @@ -2827,7 +2830,7 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C call register_restart_field(visc%Lam2, "Lam2", .false., restart_CS, & "(Langmuir Number)^-2", units="" ) endif - if (MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. & + if (MLE_use_PBL_MLD .or. do_brine_plume .or. use_nlVstress .or. & use_neutral_diffusion .or. use_hor_bnd_diff .or. MLE_use_Bodner) then call register_restart_field(visc%h_ML, "h_ML", .false., restart_CS, & "Instantaneous active mixing layer thickness", & diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 9eba7e3471..c94a449012 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -43,7 +43,7 @@ module MOM_vert_friction public vertvisc, vertvisc_remnant, vertvisc_coef public vertvisc_limit_vel, vertvisc_init, vertvisc_end public updateCFLtruncationValue -public vertFPmix +public vertNLstress ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -160,10 +160,6 @@ module MOM_vert_friction !! u-accelerations are written if velocity truncations occur. character(len=200) :: v_trunc_file !< The complete path to a file in which a column of !! v-accelerations are written if velocity truncations occur. - logical :: StokesMixing !< If true, do Stokes drift mixing via the Lagrangian current - !! (Eulerian plus Stokes drift). False by default and set - !! via STOKES_MIXING_COMBINED. - type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. real, allocatable, dimension(:,:) :: kappa_gl90_2d !< 2D kappa_gl90 at h-points [L2 H Z-1 T-1 ~> m2 s-1 or Pa s] @@ -198,29 +194,24 @@ module MOM_vert_friction contains !> Add nonlocal stress increments to ui^n and vi^n. -subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, lpost, Cemp_NL, G, GV, US, CS, OBC, Waves) +subroutine vertNLstress(ui, vi, hbl_h, h, forces, dt, lpost, Cemp_NL, G, GV, US, CS, OBC, Waves) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: ui !< Zonal velocity after vertvisc [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(inout) :: vi !< Meridional velocity after vertvisc [L T-1 ~> m s-1] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: uold !< Old Zonal velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(inout) :: vold !< Old Meridional velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: hbl_h !< boundary layer depth [H ~> m] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces real, intent(in) :: dt !< Time increment [T ~> s] real, intent(in) :: Cemp_NL !< empirical coefficient of non-local momentum mixing [nondim] - logical, intent(in) :: lpost !< Compute and make available FPMix diagnostics + logical, intent(in) :: lpost !< Compute and make available vertNLstress diagnostics type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure - type(wave_parameters_CS), & - optional, pointer :: Waves !< Container for wave/Stokes information + type(wave_parameters_CS), pointer :: Waves !< Container for wave/Stokes information ! local variables real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth (u-pts) [H ~> m] @@ -256,6 +247,9 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, lpost, Cemp_NL, G is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke + if (.not.associated(Waves)) call MOM_error(FATAL, "MOM_vert_friction, vertNLstress: "//& + "an associated Waves control structure is required (e.g., USE_WAVES=True).") + pi = 4. * atan2(1.,1.) Irho0 = 1.0 / GV%Rho0 @@ -418,7 +412,7 @@ subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, lpost, Cemp_NL, G endif -end subroutine vertFPmix +end subroutine vertNLstress !> Compute coupling coefficient associated with vertical viscosity parameterization as in Greatbatch and Lamb @@ -542,7 +536,7 @@ end subroutine find_coupling_coef_gl90 !! There is an additional stress term on the right-hand side !! if DIRECT_STRESS is true, applied to the surface layer. subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & - taux_bot, tauy_bot, fpmix, Waves) + taux_bot, tauy_bot, Waves) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -566,7 +560,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real, dimension(SZI_(G),SZJB_(G)), & optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to !! rock [R L Z T-2 ~> Pa] - logical, optional, intent(in) :: fpmix !< fpmix along Eulerian shear type(wave_parameters_CS), & optional, pointer :: Waves !< Container for wave/Stokes information @@ -611,8 +604,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real, allocatable, dimension(:,:,:) :: KE_v ! The area integral of a KE term in a layer at v-points ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3] - logical :: DoStokesMixing - logical :: lfpmix + logical :: mixEulerianShear ! If true, remove the Stokes drift from the velocities before + ! the implicit vertical viscosity solver and restore it + ! afterward, so that the mixing acts on the Eulerian shear. integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec @@ -642,30 +636,19 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & accel_underflow = CS%vel_underflow * Idt - !Check if Stokes mixing allowed if requested (present and associated) - DoStokesMixing=.false. - if (CS%StokesMixing) then - if (present(Waves)) DoStokesMixing = associated(Waves) - if (.not. DoStokesMixing) & - call MOM_error(FATAL,"Stokes Mixing called without allocated"//& - "Waves Control Structure") - endif - lfpmix = .false. - if ( present(fpmix) ) lfpmix = fpmix + ! Determine whether the implicit vertical mixing should act on the Eulerian shear, + ! which is obtained by removing the Stokes drift from the model (Lagrangian mean) + ! velocities before the solver and restoring it afterward. + mixEulerianShear = .false. + if (present(Waves)) then ; if (associated(Waves)) then + if (Waves%UseWaves) mixEulerianShear = .not. Waves%LagrangianMixing + endif ; endif ! Update the zonal velocity component using a modification of a standard ! tridiagonal solver. - ! WGL: Brandon Reichl says the following is obsolete. u(I,j,k) already - ! includes Stokes. - ! When mixing down Eulerian current + Stokes drift add before calling solver - if (DoStokesMixing) then - do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then - u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) - endif ; enddo ; enddo ; enddo - endif - - if (lfpmix) then + ! When mixing down the Eulerian shear, remove the Stokes drift before calling the solver + if (mixEulerianShear) then do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) endif ; enddo ; enddo ; enddo @@ -863,14 +846,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & endif endif - ! When mixing down Eulerian current + Stokes drift subtract after calling solver - if (DoStokesMixing) then - do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then - u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) - endif ; enddo ; enddo ; enddo - endif - - if (lfpmix) then + ! When mixing down the Eulerian shear, restore the Stokes drift after calling the solver + if (mixEulerianShear) then do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) endif ; enddo ; enddo ; enddo @@ -878,14 +855,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & ! == Now work on the meridional velocity component. - ! When mixing down Eulerian current + Stokes drift add before calling solver - if (DoStokesMixing) then - do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then - v(i,j,k) = v(i,j,k) + Waves%Us_y(i,j,k) - endif ; enddo ; enddo ; enddo - endif - - if (lfpmix) then + ! When mixing down the Eulerian shear, remove the Stokes drift before calling the solver + if (mixEulerianShear) then do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then v(i,j,k) = v(i,j,k) - Waves%Us_y(i,j,k) endif ; enddo ; enddo ; enddo @@ -1055,14 +1026,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & endif endif - ! When mixing down Eulerian current + Stokes drift subtract after calling solver - if (DoStokesMixing) then - do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then - v(i,J,k) = v(i,J,k) - Waves%Us_y(i,J,k) - endif ; enddo ; enddo ; enddo - endif - - if (lfpmix) then + ! When mixing down the Eulerian shear, restore the Stokes drift after calling the solver + if (mixEulerianShear) then do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then v(i,J,k) = v(i,J,k) + Waves%Us_y(i,J,k) endif ; enddo ; enddo ; enddo @@ -2712,7 +2677,7 @@ end subroutine vertvisc_limit_vel !> Initialize the vertical friction module subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & - ntrunc, CS, fpmix) + ntrunc, CS, nlVstress) type(ocean_internal_state), & target, intent(in) :: MIS !< The "MOM Internal State", a set of pointers !! to the fields and accelerations that make @@ -2727,7 +2692,9 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & type(directories), intent(in) :: dirs !< Relevant directory paths integer, target, intent(inout) :: ntrunc !< Number of velocity truncations type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure - logical, optional, intent(in) :: fpmix !< Nonlocal momentum mixing + logical, optional, intent(in) :: nlVstress !< If true, nonlocal momentum flux + !! increments are in use, so register + !! the related diagnostics ! Local variables @@ -2735,7 +2702,9 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & real :: Kv_back_z ! A background kinematic viscosity [Z2 T-1 ~> m2 s-1] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz - logical :: lfpmix + logical :: l_nlVstress + logical :: StokesMixing ! The value of STOKES_MIXING_COMBINED, an option whose partial + ! implementation has been removed; enabling it is/was a fatal error. character(len=200) :: kappa_gl90_file, inputdir, kdgl90_varname ! This include declares and sets the variable "version". # include "version_variable.h" @@ -2760,8 +2729,8 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & CS%diag => diag ; CS%ntrunc => ntrunc ; ntrunc = 0 - lfpmix = .false. - if (present(fpmix)) lfpmix = fpmix + l_nlVstress = .false. + if (present(nlVstress)) l_nlVstress = nlVstress ! Default, read and log parameters call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.) @@ -2982,24 +2951,19 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & "The start value of the truncation CFL number used when "//& "ramping up CFL_TRUNC.", & units="nondim", default=0.) - call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", CS%StokesMixing, & + call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", StokesMixing, & "Flag to use Stokes drift Mixing via the Lagrangian "//& " current (Eulerian plus Stokes drift). "//& " Still needs work and testing, so not recommended for use.",& default=.false.) - !BGR 04/04/2018{ - ! StokesMixing is required for MOM6 for some Langmuir mixing parameterization. - ! The code used here has not been developed for vanishing layers or in - ! conjunction with any bottom friction. Therefore, the following line is - ! added so this functionality cannot be used without user intervention in - ! the code. This will prevent general use of this functionality until proper - ! care is given to the previously mentioned issues. Comment out the following - ! MOM_error to use, but do so at your own risk and with these points in mind. - !} - if (CS%StokesMixing) then - call MOM_error(FATAL, "Stokes mixing requires user intervention in the code.\n"//& - " Model now exiting. See MOM_vert_friction.F90 for \n"//& - " details (search 'BGR 04/04/2018' to locate comment).") + ! The partial implementation of STOKES_MIXING_COMBINED was incomplete (it had not been + ! developed for vanishing layers or in conjunction with any bottom friction) and had been + ! deliberately unreachable since 2018, with a FATAL error here that could only be bypassed + ! by editing the source code. The unreachable code has been removed from vertvisc. + if (StokesMixing) then + call MOM_error(FATAL, "STOKES_MIXING_COMBINED is not implemented. Its partial \n"//& + " implementation was unusable and has been removed from \n"//& + " MOM_vert_friction.F90.") endif call get_param(param_file, mdl, "VEL_UNDERFLOW", CS%vel_underflow, & "A negligibly small velocity magnitude below which velocity "//& @@ -3059,7 +3023,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', & thickness_units, conversion=US%Z_to_m) - if (lfpmix) then + if (l_nlVstress) then CS%id_uE_h = register_diag_field('ocean_model', 'uE_h' , CS%diag%axesTL, & Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) CS%id_vE_h = register_diag_field('ocean_model', 'vE_h' , CS%diag%axesTL, & diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index d7ce6cb95e..5c38e71b1b 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -120,6 +120,13 @@ module MOM_wave_interface integer, public :: NumBands = 0 !< Number of wavenumber/frequency partitions !! Must match the number of bands provided !! via either coupling or file. + logical, public :: LagrangianMixing = .true. !< True if the implicit vertical mixing of + !! momentum acts on the model (Lagrangian mean) current, + !! as in Reichl et al., 2016 KPP-LT approach. If false + !! and waves are in use, the Stokes drift is removed + !! before the implicit vertical viscosity solver and + !! restored afterward, so that the mixing acts on the + !! Eulerian shear. ! The remainder of this control structure is private integer :: WaveMethod = -99 !< Options for including wave information @@ -129,11 +136,6 @@ module MOM_wave_interface !! 2 - DHH85 !! 3 - LF17 !! -99 - No waves computed, but empirical Langmuir number used. - logical :: LagrangianMixing !< This feature is in development and not ready - !! True if Stokes drift is present and mixing - !! should be applied to Lagrangian current - !! (mean current + Stokes drift). - !! See Reichl et al., 2016 KPP-LT approach logical :: StokesMixing !< This feature is in development and not ready. !! True if vertical mixing of momentum !! should be applied directly to Stokes current @@ -299,6 +301,8 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags logical :: use_waves logical :: StatisticalWaves + logical :: fpmix ! The value of the FPMIX runtime parameter, which determines + ! the default for LAGRANGIAN_MIXING. (FPMIX will be obsolete in the future.) ! Dummy Check if (.not. associated(CS)) then @@ -361,13 +365,19 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) ! Wave modified physics ! Presently these are all in research mode + ! The FPMIX parameter historically enabled both the nonlocal momentum flux + ! increments (now NL_VSTRESS) and the mixing of the Eulerian shear in vertvisc (now + ! LAGRANGIAN_MIXING=False). To reproduce the behavior of existing configurations that + ! set FPMIX, it is honored here as setting the default of LAGRANGIAN_MIXING; an + ! explicitly set LAGRANGIAN_MIXING takes precedence. + fpmix = .false. + call get_param(param_file, mdl, "FPMIX", fpmix, default=.false., do_not_log=.true.) call get_param(param_file, mdl, "LAGRANGIAN_MIXING", CS%LagrangianMixing, & - "Flag to use Lagrangian Mixing of momentum", default=.false., & - do_not_log=.not.use_waves) - if (CS%LagrangianMixing) then - ! Force Code Intervention - call MOM_error(FATAL,"Should you be enabling Lagrangian Mixing? Code not ready.") - endif + "If true, the implicit vertical mixing of momentum acts on the model "//& + "(Lagrangian mean) current. If false and waves are in use, the Stokes "//& + "drift is removed before the implicit vertical viscosity solver and "//& + "restored afterward, so that the mixing acts on the Eulerian shear.", & + default=.not.fpmix, do_not_log=.not.use_waves) call get_param(param_file, mdl, "STOKES_MIXING", CS%StokesMixing, & "Flag to use Stokes Mixing of momentum", default=.false., & do_not_log=.not.use_waves)