From 15deea43ecb97a3920b3a06162363f1bac2b8ba5 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Sun, 8 Sep 2024 23:21:05 -0600 Subject: [PATCH] *Updates in FPMix and Stokes Most (#283) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This PR revises the formulation of the legacy K-profile parameterization (KPP) ocean boundary layer scheme. It incorporates: 1. a non-local momentum flux—the Flux-profile parameterization (`FPMIX`); when the local shear is not aligned with the wind, this scheme adds a non-local momentum flux in the direction of the wind; and 2. mixing with and without waves following the Monin-Obukhov Similarity Theory expanded to include Stokes drift (`STOKES_MOST`). This option provides the transition from waveless to ocean surface waves in any stage of growth or decay. **Summary:** * Uncomment omega w2x entries; * Simplify the nonlocal increments in `vertFPMix`; * In the call to `CVmix_kpp_compute_unresolved_shear`, passes the 2D surface buoyancy flux (`surfBuoyFlux2`) instead of the 1D version (`surfBuoyFlux`), which is preferable. **This is answer changing**; * Remove `uold` and `vold` diagnostics. These were used in an alternative time-stepping scheme that is now obsolete; * Pass boundary layer depths to the RK2 and add consistency check to make sure `FPMix` is always used with `SPLIT`; * Add the capability to mix down the Eulerian gradient instead of the Lagrangian; * Make a minimum set of `FPMix` diagnostics available. This PR relies on https://github.com/CVMix/CVMix-src/pull/94/. New diagnostics: ``` "StokesXI" ! modules: {ocean_model,ocean_model_d2} ! long_name: Stokes Similarity Parameter ! units: nondim ! cell_methods: xh:mean yh:mean area:mean "Lam2" ! modules: {ocean_model,ocean_model_d2} ! long_name: Ustk0_ustar ! units: nondim ! cell_methods: xh:mean yh:mean area:mean "uE_h" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: x-zonal Eulerian ! units: m s-1 ! cell_methods: xh:mean yh:mean zl:mean area:mean ! variants: {uE_h,uE_h_xyave} "vE_h" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: y-merid Eulerian ! units: m s-1 ! cell_methods: xh:mean yh:mean zl:mean area:mean ! variants: {vE_h,vE_h_xyave} "uInc_h" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: x-zonal Eulerian ! units: m s-1 ! cell_methods: xh:mean yh:mean zl:mean area:mean ! variants: {uInc_h,uInc_h_xyave} "vInc_h" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: x-zonal Eulerian ! units: m s-1 ! cell_methods: xh:mean yh:mean zl:mean area:mean ! variants: {vInc_h,vInc_h_xyave} "uStk" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: x-FP du increment ! units: m s-1 ! cell_methods: xh:mean yh:mean zl:mean area:mean ! variants: {uStk,uStk_xyave} "vStk" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: y-FP dv increment ! units: m s-1 ! cell_methods: xh:mean yh:mean zl:mean area:mean ! variants: {vStk,vStk_xyave} "Omega_tau2s" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: Stress direction from shear ! units: radians ! cell_methods: xh:mean yh:mean zi:point area:mean ! variants: {Omega_tau2s,Omega_tau2s_xyave} "Omega_tau2w" ! modules: {ocean_model,ocean_model_z,ocean_model_rho2,ocean_model_d2,ocean_model_z_d2,ocean_model_rho2_d2} ! long_name: Stress direction from wind ! units: radians ! cell_methods: xh:mean yh:mean zi:point area:mean ! variants: {Omega_tau2w,Omega_tau2w_xyave} "uStk0" ! modules: {ocean_model,ocean_model_d2} ! long_name: Zonal Surface Stokes ! units: m s-1 ! cell_methods: xh:mean yh:mean area:mean "vStk0" ! modules: {ocean_model,ocean_model_d2} ! long_name: Merid Surface Stokes ! units: m s-1 ! cell_methods: xh:mean yh:mean area:mean ``` --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 6 +- pkg/CVMix-src | 2 +- src/core/MOM.F90 | 19 +- src/core/MOM_dynamics_split_RK2.F90 | 88 ++- src/core/MOM_forcing_type.F90 | 40 +- .../vertical/MOM_CVMix_KPP.F90 | 288 +++++++-- .../vertical/MOM_vert_friction.F90 | 583 +++++++----------- src/user/MOM_wave_interface.F90 | 2 +- 8 files changed, 552 insertions(+), 476 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 9f409a1af9..897491711f 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -325,7 +325,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, press=.true., fix_accum_bug=.not.CS%ustar_gustless_bug, & cfc=CS%use_CFC, marbl=CS%use_marbl_tracers, hevap=CS%enthalpy_cpl, & tau_mag=.true., ice_ncat=IOB%ice_ncat) - !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) call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed) @@ -762,7 +762,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) @@ -923,7 +923,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. diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 52aac958e0..3ec78bac83 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 52aac958e05cdb2471dc73f9ef7fb4e816c550f2 +Subproject commit 3ec78bac8306ef2e61a33e0c4beafa0875a2c787 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2cbbf6bcc7..1c94cf18fb 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -78,6 +78,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_split_RK2b, only : step_MOM_dyn_split_RK2b, register_restarts_dyn_split_RK2b use MOM_dynamics_split_RK2b, only : initialize_dyn_split_RK2b, end_dyn_split_RK2b use MOM_dynamics_split_RK2b, only : MOM_dyn_split_RK2b_CS, remap_dyn_split_RK2b_aux_vars @@ -2102,6 +2103,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 @@ -2206,6 +2208,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, & @@ -3334,6 +3346,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, the following is needed to get BLDs into the dynamics module + 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) @@ -3609,7 +3626,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) then + 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) endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 35e9c2722e..217ec42c20 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -12,6 +12,7 @@ module MOM_dynamics_split_RK2 use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member use MOM_diag_mediator, only : diag_mediator_init, enable_averages use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : post_product_u, post_product_sum_u @@ -45,7 +46,9 @@ module MOM_dynamics_split_RK2 use MOM_continuity, only : continuity_init, continuity_stencil use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_CS use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS use MOM_debugging, only : check_redundant +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS, hor_visc_vel_stencil @@ -137,14 +140,16 @@ module MOM_dynamics_split_RK2 real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure !! anomaly in each layer due to free surface height !! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]. + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to ge + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure - real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean - !! to the seafloor [R L Z T-2 ~> Pa] - real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean - !! to the seafloor [R L Z T-2 ~> Pa] - type(BT_cont_type), pointer :: BT_cont => NULL() !< A structure with elements that describe the - !! effective summed open face areas as a function - !! of barotropic flow. + real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean + !! to the seafloor [R L Z T-2 ~> Pa] + real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean + !! to the seafloor [R L Z T-2 ~> Pa] + type(BT_cont_type), pointer :: BT_cont => NULL() !< A structure with elements that describe the + !! effective summed open face areas as a function + !! of barotropic flow. ! This is to allow the previous, velocity-based coupling with between the ! baroclinic and barotropic modes. @@ -174,13 +179,13 @@ module MOM_dynamics_split_RK2 !! the extent to which the treatment of gravity waves !! is forward-backward (0) or simulated backward !! Euler (1) [nondim]. 0 is often used. - logical :: debug !< If true, write verbose checksums for debugging purposes. + real :: Cemp_NL !< Empirical coefficient of non-local momentum mixing [nondim] + 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 @@ -267,6 +272,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 @@ -394,7 +400,8 @@ 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 + ! in the corrector step (not the predict) integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: cont_stencil, obc_stencil, vel_stencil @@ -710,16 +717,22 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f 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(visc%h_ML)) hbl(:,:) = visc%h_ML(:,:) - 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, & - GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) + 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) endif if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)") @@ -960,12 +973,15 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f 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) - 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, & - G, GV, US, CS%vertvisc_CSp, CS%OBC) + 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) endif @@ -1052,11 +1068,6 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f 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 @@ -1290,6 +1301,17 @@ 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. +!! Needed when BLDs are used in the dynamics. +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, tv, uh, vh, eta, Time, G, GV, US, param_file, & @@ -1402,8 +1424,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p "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.) + if (CS%fpmix) then + call get_param(param_file, "MOM", "CEMP_NL", CS%Cemp_NL, & + "Empirical coefficient of non-local momentum mixing.", & + units="nondim", default=3.6) + endif 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 "//& @@ -1480,7 +1507,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & 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) + call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, & + 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) ) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index bef29ffe86..d25df710ce 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -81,7 +81,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 @@ -263,8 +263,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) @@ -408,7 +408,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 @@ -1577,8 +1577,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 @@ -2509,11 +2509,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) @@ -2661,11 +2661,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) @@ -3367,8 +3367,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) @@ -3748,7 +3748,7 @@ end subroutine myAlloc_3d 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) @@ -3821,7 +3821,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) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index fec7219e7d..816d5d7498 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -30,6 +30,7 @@ module MOM_CVMix_KPP use CVMix_kpp, only : CVMix_kpp_compute_unresolved_shear use CVMix_kpp, only : CVMix_kpp_params_type use CVMix_kpp, only : CVMix_kpp_compute_kOBL_depth +use CVMix_kpp, only : CVMix_kpp_compute_StokesXi implicit none ; private @@ -82,6 +83,7 @@ module MOM_CVMix_KPP logical :: enhance_diffusion !< If True, add enhanced diffusivity at base of boundary layer. character(len=32) :: interpType !< Type of interpolation to compute bulk Richardson number character(len=32) :: interpType2 !< Type of interpolation to compute diff and visc at OBL_depth + logical :: StokesMOST !< If True, use Stokes similarity package logical :: computeEkman !< If True, compute Ekman depth limit for OBLdepth logical :: computeMoninObukhov !< If True, compute Monin-Obukhov limit for OBLdepth logical :: passiveMode !< If True, makes KPP passive meaning it does NOT alter the diffusivity @@ -145,11 +147,15 @@ module MOM_CVMix_KPP integer :: id_EnhW = -1 integer :: id_La_SL = -1 integer :: id_OBLdepth_original = -1 + integer :: id_StokesXI = -1 + integer :: id_Lam2 = -1 !>@} ! Diagnostics arrays real, allocatable, 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* real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent [nondim] real, allocatable, dimension(:,:) :: OBLdepthprev !< previous Depth (positive) of OBL [Z ~> m] real, allocatable, dimension(:,:) :: La_SL !< Langmuir number used in KPP [nondim] @@ -272,6 +278,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 'Type of interpolation to compute diff and visc at OBL_depth.\n'// & 'Allowed types are: linear, quadratic, cubic or LMD94.', & default='LMD94') + call get_param(paramFile, mdl, 'STOKES_MOST', CS%StokesMOST, & + 'If True, use Stokes Similarity package.', & + default=.False.) call get_param(paramFile, mdl, 'COMPUTE_EKMAN', CS%computeEkman, & 'If True, limit OBL depth to be no deeper than Ekman depth.', & default=.False.) @@ -498,6 +507,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) interp_type=CS%interpType, & interp_type2=CS%interpType2, & lEkman=CS%computeEkman, & + lStokesMOST=CS%StokesMOST, & lMonOb=CS%computeMoninObukhov, & MatchTechnique=CS%MatchTechnique, & lenhanced_diff=CS%enhance_diffusion,& @@ -524,6 +534,12 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', & cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme') endif + if( CS%StokesMOST ) then + CS%id_StokesXI = register_diag_field('ocean_model', 'StokesXI', diag%axesT1, Time, & + 'Stokes Similarity Parameter', 'nondim') + CS%id_Lam2 = register_diag_field('ocean_model', 'Lam2', diag%axesT1, Time, & + 'Ustk0_ustar', 'nondim') + endif CS%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, Time, & 'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', & 'kg/m3', conversion=US%R_to_kg_m3) @@ -584,6 +600,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 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. ) allocate( CS%La_SL( SZI_(G), SZJ_(G) ), source=0. ) allocate( CS%Vt2( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. ) @@ -804,6 +822,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, & GV%ke, & ! (in) Number of levels to compute coeffs for GV%ke, & ! (in) Number of levels in array shape Langmuir_EFactor=LangEnhK,& ! Langmuir enhancement multiplier + StokesXi = CS%StokesParXI(i,j), & ! Stokes forcing parameter CVMix_kpp_params_user=CS%KPP_params ) ! safety check, Kviscosity and Kdiffusivity must be >= 0 @@ -962,7 +981,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m] - ! Variables for passing to CVMix routines, often in MKS units real, dimension( GV%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars in MKS units [m s-1] real, dimension( GV%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3] @@ -997,6 +1015,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real :: Uk, Vk ! Layer velocities relative to their averages in the surface layer [L T-1 ~> m s-1] real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth [Z ~> m] real :: hTot ! Running sum of thickness used in the surface layer average [Z ~> m] + real :: I_hTot ! The inverse of hTot [Z-1 ~> m-1] real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] real :: delH ! Thickness of a layer [Z ~> m] real :: surfTemp ! Average of temperature over the surface layer [C ~> degC] @@ -1018,6 +1037,17 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl integer :: i, j, k, km1, kk, ksfc, ktmp ! Loop indices + real, dimension(GV%ke) :: uE_H, vE_H ! Eulerian velocities h-points, centers [L T-1 ~> m s-1] + real, dimension(GV%ke) :: uS_H, vS_H ! Stokes drift components h-points, centers [L T-1 ~> m s-1] + real, dimension(GV%ke) :: uSbar_H, vSbar_H ! Cell Average Stokes drift h-points [L T-1 ~> m s-1] + real, dimension(GV%ke+1) :: uS_Hi, vS_Hi ! Stokes Drift components at interfaces [L T-1 ~> m s-1] + real :: uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD ! Stokes at/to to Surface Layer Extent + ! [L T-1 ~> m s-1] + real :: StokesXI ! Stokes similarity parameter [nondim] + real, dimension( GV%ke ) :: StokesXI_1d , StokesVt_1d ! Parameters of TKE production ratio [nondim] + real :: Llimit ! Stable boundary Layer Limit = vonk Lstar [Z ~> m] + integer :: kbl ! index of cell containing boundary layer depth + if (CS%Stokes_Mixing .and. .not.associated(Waves)) call MOM_error(FATAL, & "KPP_compute_BLD: The Waves control structure must be associated if STOKES_MIXING is True.") @@ -1046,27 +1076,36 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl !$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, & !$OMP surfBuoyFlux, U_H, V_H, Coriolis, pRef, SLdepth_0d, vt2_1d, & !$OMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, & - !$OMP surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv, & + !$OMP surfHvS, hTot, I_hTot, delH, surftemp, surfsalt, surfu, surfv, & !$OMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, N_col, & !$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_guess, LA, rho_1D, & !$OMP deltarho, deltaBuoy, N2_1d, ws_1d, LangEnhVT2,KPP_OBL_depth, z_cell, & - !$OMP z_inter, OBL_depth, BulkRi_1d, zBottomMinusOffset) & + !$OMP z_inter, OBL_depth, BulkRi_1d, zBottomMinusOffset, uE_H, vE_H, & + !$OMP uS_H, vS_H, uSbar_H, vSbar_H , uS_Hi, vS_Hi, & + !$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) do j = G%jsc, G%jec do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then do k=1,GV%ke - U_H(k) = 0.5 * (u(i,j,k)+u(i-1,j,k)) - V_H(k) = 0.5 * (v(i,j,k)+v(i,j-1,k)) + U_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)) + V_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)) enddo + if (CS%StokesMOST) then + do k=1,GV%ke + uE_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)-Waves%US_x(I,j,k)-Waves%US_x(I-1,j,k)) + 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)) ) surfFricVel = US%Z_to_m*US%s_to_T * uStar(i,j) - ! Bullk Richardson number computed for each cell in a column, + ! Bulk Richardson number computed for each cell in a column, ! assuming OBLdepth = grid cell depth. After Rib(k) is ! known for the column, then CVMix interpolates to find ! the actual OBLdepth. This approach avoids need to iterate @@ -1075,8 +1114,11 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl iFaceHeight(1) = 0.0 ! BBL is all relative to the surface pRef = 0. ; if (associated(tv%p_surf)) pRef = tv%p_surf(i,j) hcorr = 0. - do k=1,GV%ke + if (CS%StokesMOST) call Compute_StokesDrift( i, j, h(i,j,1) , iFaceHeight(1), & + uS_Hi(1), vS_Hi(1), uS_H(1), vS_H(1), uSbar_H(1), vSbar_H(1), Waves) + + do k=1,GV%ke ! cell center and cell bottom in meters (negative values in the ocean) dh = dz(i,j,k) ! Nominal thickness to use for increment dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) @@ -1095,53 +1137,99 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl endif enddo - ! average temperature, salinity, u and v over surface layer - ! use C-grid average to get u and v on T-points. - surfHtemp = 0.0 - surfHsalt = 0.0 - surfHu = 0.0 - surfHv = 0.0 - surfHuS = 0.0 - surfHvS = 0.0 - hTot = 0.0 - do ktmp = 1,ksfc - - ! SLdepth_0d can be between cell interfaces - delH = min( max(0.0, SLdepth_0d - hTot), dz(i,j,ktmp) ) - - ! surface layer thickness - hTot = hTot + delH - - ! surface averaged fields - surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH - surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH - surfHu = surfHu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delH - surfHv = surfHv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delH + if (CS%StokesMOST) then + surfBuoyFlux = buoy_scale * & + (buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,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) + call Compute_StokesDrift(i,j, iFaceHeight(ksfc) , -SLdepth_0d, & + uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD, Waves) + call cvmix_kpp_compute_StokesXi( iFaceHeight,CellHeight,ksfc ,SLdepth_0d,surfBuoyFlux, & + 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 + + ! average temperature, salinity, u and v over surface layer starting at ksfc + delH = SLdepth_0d + iFaceHeight(ksfc-1) + surfHtemp = Temp(i,j,ksfc) * delH + surfHsalt = Salt(i,j,ksfc) * delH + surfHu = (uE_H(ksfc) + uSbar_SLD) * delH + surfHv = (vE_H(ksfc) + vSbar_SLD) * delH + hTot = delH + do ktmp = 1,ksfc-1 ! if ksfc >=2 + delH = h(i,j,ktmp)*GV%H_to_Z + hTot = hTot + delH + surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH + surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH + surfHu = surfHu + (uE_H(ktmp) + uSbar_H(ktmp)) * delH + surfHv = surfHv + (vE_H(ktmp) + vSbar_H(ktmp)) * delH + enddo + I_hTot = 1./hTot + surfTemp = surfHtemp * I_hTot + 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 + surfHsalt = 0.0 + surfHu = 0.0 + surfHv = 0.0 + surfHuS = 0.0 + surfHvS = 0.0 + hTot = 0.0 + do ktmp = 1,ksfc + + ! SLdepth_0d can be between cell interfaces + delH = min( max(0.0, SLdepth_0d - hTot), dz(i,j,ktmp) ) + + ! surface layer thickness + hTot = hTot + delH + + ! surface averaged fields + surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH + surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH + surfHu = surfHu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delH + surfHv = surfHv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delH + if (CS%Stokes_Mixing) then + surfHus = surfHus + 0.5*(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) * delH + surfHvs = surfHvs + 0.5*(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) * delH + endif + + enddo + 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 + Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV + if (CS%Stokes_Mixing) then - surfHus = surfHus + 0.5*(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) * delH - surfHvs = surfHvs + 0.5*(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) * delH + ! If momentum is mixed down the Stokes drift gradient, then + ! the Stokes drift must be included in the bulk Richardson number + ! calculation. + Uk = Uk + (0.5*(Waves%Us_x(i,j,k)+Waves%US_x(i-1,j,k)) - surfUs ) + Vk = Vk + (0.5*(Waves%Us_y(i,j,k)+Waves%Us_y(i,j-1,k)) - surfVs ) endif - enddo - 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 - Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV - - if (CS%Stokes_Mixing) then - ! If momentum is mixed down the Stokes drift gradient, then - ! the Stokes drift must be included in the bulk Richardson number - ! calculation. - Uk = Uk + (0.5*(Waves%Us_x(i,j,k)+Waves%US_x(i-1,j,k)) - surfUs ) - Vk = Vk + (0.5*(Waves%Us_y(i,j,k)+Waves%Us_y(i,j-1,k)) - surfVs ) - endif + ! this difference accounts for penetrating SW + surfBuoyFlux = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,k+1)) + surfBuoyFlux2(k) = surfBuoyFlux + + endif ! StokesMOST deltaU2(k) = US%L_T_to_m_s**2 * (Uk**2 + Vk**2) @@ -1165,9 +1253,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! iterate pRef for next pass through k-loop. pRef = pRef + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k) - ! this difference accounts for penetrating SW - surfBuoyFlux2(k) = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,k+1)) - enddo ! k-loop finishes if ( (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) .and. .not. present(lamult)) then @@ -1215,11 +1300,12 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! 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. - call CVMix_kpp_compute_turbulent_scales( & + call CVMix_kpp_compute_turbulent_scales( & ! 1d_OBL CS%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext 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] + xi=StokesVt_1d, & ! (in) Stokes similarity parameter-->1/CHI(xi) enhance of Vt w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1] CVMix_kpp_params_user=CS%KPP_params ) @@ -1255,10 +1341,17 @@ 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=surfBuoyFlux, & ! 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 +! ! A hack to avoid KPP reaching the bottom. It was needed during development +! ! because KPP was unable to handle vanishingly small layers near the bottom. +! if (CS%deepOBLoffset>0.) then +! zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(CS%deepOBLoffset, -0.1*iFaceHeight(GV%ke+1)) +! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset ) +! endif + zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(GV%ke+1)) call CVMix_kpp_compute_OBL_depth( & BulkRi_1d, & ! (in) Bulk Richardson number @@ -1267,11 +1360,39 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent zt_cntr=z_cell, & ! (in) Height of cell centers [m] surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] - surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] + surf_buoy=surfBuoyFlux2, & ! (in) Buoyancy flux at surface [m2 s-3] Coriolis=Coriolis, & ! (in) Coriolis parameter [s-1] + Xi = StokesXI_1d, & ! (in) Stokes similarity parameter Lmob limit (1-Xi) + zBottom = zBottomMinusOffset, & ! (in) Numerical limit on OBLdepth CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters CS%OBLdepth(i,j) = US%m_to_Z * KPP_OBL_depth + if (CS%StokesMOST) then + kbl = nint(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 + ksfc = kbl + do ktmp = 1, kbl + if (-1.0*iFaceHeight(ktmp+1) >= SLdepth_0d) then + ksfc = ktmp + exit + endif + enddo + + call Compute_StokesDrift(i,j, iFaceHeight(ksfc) , -SLdepth_0d, & + uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD, Waves) + call cvmix_kpp_compute_StokesXi( iFaceHeight,CellHeight,ksfc ,SLdepth_0d, & + surfBuoyFlux, 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 ) + CS%StokesParXI(i,j) = StokesXI + CS%Lam2(i,j) = sqrt(US_Hi(1)**2+VS_Hi(1)**2) / MAX(surfFricVel,0.0002) + + else !.not Stokes_MOST + CS%StokesParXI(i,j) = 10.0 + CS%Lam2(i,j) = sqrt(US_Hi(1)**2+VS_Hi(1)**2) / MAX(surfFricVel,0.0002) + ! A hack to avoid KPP reaching the bottom. It was needed during development ! because KPP was unable to handle vanishingly small layers near the bottom. if (CS%deepOBLoffset>0.) then @@ -1285,6 +1406,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(GV%ke+1) ) ! no deeper than bottom CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + endif !Stokes_MOST + ! compute unresolved squared velocity for diagnostics if (CS%id_Vt2 > 0) then Vt2_1d(:) = CVmix_kpp_compute_unresolved_shear( & @@ -1293,7 +1416,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl N_iface=N_col, & ! Buoyancy frequency at interface [s-1] EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim] LaSL=CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim] - bfsfc=surfBuoyFlux, & ! 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 CS%Vt2(i,j,:) = US%m_to_Z*US%T_to_s * Vt2_1d(:) @@ -1307,6 +1430,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl US%Z_to_m*CS%OBLdepth(i,j), & ! (in) OBL depth [m] surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] + xi=StokesXI, & ! (in) Stokes similarity parameter-->1/CHI(xi) enhance w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1] CVMix_kpp_params_user=CS%KPP_params) ! KPP parameters CS%Ws(i,j,:) = US%m_to_Z*US%T_to_s*Ws_1d(:) @@ -1342,6 +1466,11 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (CS%id_La_SL > 0) call post_data(CS%id_La_SL, CS%La_SL, CS%diag) if (CS%id_Vt2 > 0) call post_data(CS%id_Vt2, CS%Vt2, CS%diag) + if (CS%StokesMOST) then + if (CS%id_StokesXI > 0) call post_data(CS%id_StokesXI, CS%StokesParXI, CS%diag) + if (CS%id_Lam2 > 0) call post_data(CS%id_Lam2 , CS%Lam2 , CS%diag) + endif + ! BLD smoothing: if (CS%n_smooth > 0) call KPP_smooth_BLD(CS, G, GV, US, dz) @@ -1594,6 +1723,49 @@ subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, end subroutine KPP_NonLocalTransport_saln +!> Compute Stokes Drift components at zbot < ztop <= 0 and at k=0.5*(ztop+zbot) and +!! average components from ztop to zbot <= 0 +subroutine Compute_StokesDrift(i ,j, ztop, zbot, uS_i, vS_i, uS_k, vS_k, uSbar, vSbar, waves) + + type(wave_parameters_CS), pointer :: waves !< Wave CS for Langmuir turbulence + real, intent(in) :: ztop !< cell top + real, intent(in) :: zbot !< cell bottom + real, intent(inout) :: uS_i !< Stokes u velocity at zbot interface + real, intent(inout) :: vS_i !< Stokes v velocity at zbot interface + real, intent(inout) :: uS_k !< Stokes u velocity at zk center + real, intent(inout) :: vS_k !< Stokes v at zk =0.5(ztop+zbot) + real, intent(inout) :: uSbar !< mean Stokes u (ztop to zbot) + real, intent(inout) :: vSbar !< mean Stokes v (ztop to zbot) + integer, intent(in) :: i !< Meridional index of H-point + integer, intent(in) :: j !< Zonal index of H-point + + ! local variables + integer :: b !< wavenumber band index + real :: fexp !< an exponential function + real :: WaveNum !< Wavenumber + + uS_i = 0.0 + vS_i = 0.0 + uS_k = 0.0 + vS_k = 0.0 + uSbar = 0.0 + vSbar = 0.0 + do b = 1, waves%NumBands + WaveNum = waves%WaveNum_Cen(b) + fexp = exp(2. * WaveNum * zbot) + uS_i = uS_i + waves%Ustk_Hb(i,j,b) * fexp + vS_i = vS_i + waves%Vstk_Hb(i,j,b) * fexp + fexp = exp( WaveNum * (ztop + zbot) ) + uS_k = uS_k+ waves%Ustk_Hb(i,j,b) * fexp + vS_k = vS_k+ waves%Vstk_Hb(i,j,b) * fexp + fexp = exp(2. * WaveNum * ztop) - exp(2. * WaveNum * zbot) + uSbar = uSbar + 0.5 * waves%Ustk_Hb(i,j,b) * fexp / WaveNum + vSbar = vSbar + 0.5 * waves%Vstk_Hb(i,j,b) * fexp / WaveNum + enddo + uSbar = uSbar / (ztop-zbot) + vSbar = vSbar / (ztop-zbot) + +end subroutine Compute_StokesDrift !> Clear pointers, deallocate memory subroutine KPP_end(CS) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index c26ee4ac75..3f968b2101 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -31,6 +31,8 @@ module MOM_vert_friction use MOM_set_visc, only : set_v_at_u, set_u_at_v use MOM_lateral_mixing_coeffs, only : VarMix_CS +use CVMix_kpp, only : cvmix_kpp_composite_Gshape + implicit none ; private #include @@ -170,9 +172,11 @@ module MOM_vert_friction integer :: id_au_vv = -1, id_av_vv = -1, id_au_gl90_vv = -1, id_av_gl90_vv = -1 integer :: id_du_dt_str = -1, id_dv_dt_str = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 - integer :: id_FPw2x = -1 !W id_FPhbl_u = -1, id_FPhbl_v = -1 - integer :: id_tauFP_u = -1, id_tauFP_v = -1 !W, id_FPtau2x_u = -1, id_FPtau2x_v = -1 - integer :: id_FPtau2s_u = -1, id_FPtau2s_v = -1, id_FPtau2w_u = -1, id_FPtau2w_v = -1 + integer :: id_Omega_w2x = -1, id_FPtau2s = -1 , id_FPtau2w = -1 + integer :: id_uE_h = -1, id_vE_h = -1 + integer :: id_uStk = -1, id_vStk = -1 + integer :: id_uStk0 = -1, id_vStk0 = -1 + integer :: id_uInc_h= -1, id_vInc_h= -1 integer :: id_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 integer :: id_Kv_gl90_u = -1, id_Kv_gl90_v = -1 @@ -191,392 +195,211 @@ module MOM_vert_friction contains -!> Add nonlocal stress increments to u^n (uold) and v^n (vold) using ui and vi. -subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, G, GV, US, CS, OBC) +!> 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) 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] + 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] - 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 + 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 + 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 ! local variables - real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth at u-pts [H ~> m] - real, dimension(SZI_(G),SZJB_(G)) :: hbl_v !< boundary layer depth at v-pts [H ~> m] - integer, dimension(SZIB_(G),SZJ_(G)) :: kbl_u !< index of the BLD at u-pts [nondim] - integer, dimension(SZI_(G),SZJB_(G)) :: kbl_v !< index of the BLD at v-pts [nondim] - real, dimension(SZIB_(G),SZJ_(G)) :: ustar2_u !< ustar squared at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G)) :: ustar2_v !< ustar squared at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G)) :: taux_u !< zonal wind stress at u-pts [R L Z T-2 ~> Pa] - real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< meridional wind stress at v-pts [R L Z T-2 ~> Pa] - !real, dimension(SZIB_(G),SZJ_(G)) :: omega_w2x_u !< angle between wind and x-axis at u-pts [rad] - !real, dimension(SZI_(G),SZJB_(G)) :: omega_w2x_v !< angle between wind and y-axis at v-pts [rad] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tau_u !< kinematic zonal mtm flux at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tau_v !< kinematic mer. mtm flux at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauxDG_u !< downgradient zonal mtm flux at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: tauyDG_u !< downgradient meri mtm flux at u-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauxDG_v !< downgradient zonal mtm flux at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: tauyDG_v !< downgradient meri mtm flux at v-pts [L2 T-2 ~> m2 s-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2s_u !< angle between mtm flux and vert shear at u-pts [rad] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2s_v !< angle between mtm flux and vert shear at v-pts [rad] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: omega_tau2w_u !< angle between mtm flux and wind at u-pts [rad] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: omega_tau2w_v !< angle between mtm flux and wind at v-pts [rad] - - real :: pi, Cemp_CG, tmp, cos_tmp, sin_tmp !< constants and dummy variables [nondim] - real :: omega_tmp !< A dummy angle [radians] - real :: du, dv !< Velocity increments [L T-1 ~> m s-1] - real :: depth !< Cumulative layer thicknesses [H ~> m or kg m=2] - real :: sigma !< Fractional depth in the mixed layer [nondim] - real :: Wind_x, Wind_y !< intermediate wind stress componenents [L2 T-2 ~> m2 s-2] - real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables [L2 T-2 ~> m2 s-2] - real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables [L2 T-2 ~> m2 s-2] - real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles [radians] - integer :: kblmin, kbld, kp1, k, nz !< vertical indices - integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq ! horizontal indices + real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth (u-pts) [H ~> m] + real, dimension(SZI_(G),SZJB_(G)) :: hbl_v !< boundary layer depth (v-pts) [H ~> m] + real, dimension(SZIB_(G),SZJ_(G)) :: taux_u !< kinematic zonal wind stress (u-pts) [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< kinematic merid wind stress (v-pts) [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G)) :: uS0 !< surface zonal Stokes drift [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G)) :: vS0 !< surface zonal Stokes drift [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uE_u !< zonal Eulerian u-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uE_h !< zonal Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vE_v !< merid Eulerian v-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vE_h !< merid Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uInc_u !< zonal Eulerian u-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uInc_h !< zonal Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vInc_v !< merid Eulerian v-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vInc_h !< merid Eulerian h-pts [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uStk !< zonal Stokes Drift (h-pts) [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vStk !< merid Stokes Drift (h-pts) [L T-1 ~> m s-1] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2s !< angle stress to shear (h-pts) [rad] + real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2w !< angle stress to wind (h-pts) [rad] + real :: pi, tmp_u, tmp_v, omega_tmp, Irho0, fexp !< constants and dummy variables + real :: sigma,Gat1,Gsig,dGdsig !< Shape parameters + real :: du, dv, depth, Wind_x, Wind_y !< intermediate variables + real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, ustar2min, tauh !< intermediate variables + real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables + real :: omega_w2s, omega_s2x, omega_tau2x, omega_s2w , omega_e2x, omega_l2x !< intermediate angles + integer :: b, kbld, kp1, k, nz !< band and vertical indices + integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq !< horizontal indices 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 pi = 4. * atan2(1.,1.) - Cemp_CG = 3.6 - kblmin = 1 - taux_u(:,:) = 0. - tauy_v(:,:) = 0. + Irho0 = 1.0 / GV%Rho0 - do j = js,je - do I = Isq,Ieq - taux_u(I,j) = forces%taux(I,j) / GV%H_to_RZ !W rho0=1035. - enddo - enddo - - do J = Jsq,Jeq - do i = is,ie - tauy_v(i,J) = forces%tauy(i,J) / GV%H_to_RZ - enddo - enddo + call pass_var(hbl_h , G%Domain, halo=1) - call pass_var( hbl_h ,G%Domain, halo=1 ) - call pass_vector(taux_u , tauy_v, G%Domain, To_All ) - ustar2_u(:,:) = 0. - ustar2_v(:,:) = 0. - hbl_u(:,:) = 0. - hbl_v(:,:) = 0. - kbl_u(:,:) = 0 - kbl_v(:,:) = 0 - !omega_w2x_u(:,:) = 0.0 - !omega_w2x_v(:,:) = 0.0 - tauxDG_u(:,:,:) = 0.0 - tauyDG_v(:,:,:) = 0.0 + ! u-points do j = js,je do I = Isq,Ieq - if( (G%mask2dCu(I,j) > 0.5) ) then - tmp = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) ) - hbl_u(I,j) = (G%mask2dT(i,j)* hbl_h(i,j) + G%mask2dT(i+1,j) * hbl_h(i+1,j)) /tmp - tmp = MAX(1.0, (G%mask2dCv(i,j) + G%mask2dCv(i,j-1) + G%mask2dCv(i+1,j) + G%mask2dCv(i+1,j-1) ) ) - tauy = ( G%mask2dCv(i ,j )*tauy_v(i ,j ) + G%mask2dCv(i ,j-1)*tauy_v(i ,j-1) & - + G%mask2dCv(i+1,j )*tauy_v(i+1,j ) + G%mask2dCv(i+1,j-1)*tauy_v(i+1,j-1) ) / tmp - ustar2_u(I,j) = sqrt( taux_u(I,j)*taux_u(I,j) + tauy*tauy ) - !omega_w2x_u(I,j) = atan2( tauy , taux_u(I,j) ) - tauxDG_u(I,j,1) = taux_u(I,j) - depth = 0.0 - do k = 1, nz - depth = depth + CS%h_u(I,j,k) - if( (depth >= hbl_u(I,j)) .and. (kbl_u(I,j) == 0 ) .and. (k > (kblmin-1)) ) then - kbl_u(I,j) = k - hbl_u(I,j) = depth - endif - enddo + taux_u(I,j) = forces%taux(I,j) * Irho0 + if ( (G%mask2dCu(I,j) > 0.5) ) then + ! h to u-pts + tmp_u = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) ) + hbl_u(I,j) = (G%mask2dT(i,j)* hbl_h(i,j) + G%mask2dT(i+1,j) * hbl_h(i+1,j)) / tmp_u endif + depth = 0. + Gat1 = 0. + do k=1, nz + ! cell center + depth = depth + 0.5*CS%h_u(I,j,k) + uE_u(I,j,k) = ui(I,j,k) - waves%Us_x(I,j,k) + if ( depth < hbl_u(I,j) ) then + sigma = depth / hbl_u(i,j) + ! cell bottom + depth = depth + 0.5*CS%h_u(I,j,k) + call cvmix_kpp_composite_Gshape(sigma,Gat1,Gsig,dGdsig) + ! nonlocal boundary-layer increment + uInc_u(I,j,k) = dt * Cemp_NL * taux_u(I,j) * dGdsig / hbl_u(I,j) + ui(I,j,k) = ui(I,j,k) + uInc_u(I,j,k) + else + uInc_u(I,j,k) = 0.0 + endif + enddo enddo enddo + + ! v-points do J = Jsq,Jeq do i = is,ie - if( (G%mask2dCv(i,J) > 0.5) ) then - tmp = max( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1))) - hbl_v(i,J) = (G%mask2dT(i,j) * hbl_h(i,J) + G%mask2dT(i,j+1) * hbl_h(i,j+1)) /tmp - tmp = max(1.0, (G%mask2dCu(i,j) + G%mask2dCu(i,j+1) + G%mask2dCu(i-1,j) + G%mask2dCu(i-1,j+1))) - taux = ( G%mask2dCu(i ,j) * taux_u(i ,j) + G%mask2dCu(i ,j+1) * taux_u(i ,j+1) & - + G%mask2dCu(i-1,j) * taux_u(i-1,j) + G%mask2dCu(i-1,j+1) * taux_u(i-1,j+1)) / tmp - ustar2_v(i,J) = sqrt(tauy_v(i,J)*tauy_v(i,J) + taux*taux) - !omega_w2x_v(i,J) = atan2( tauy_v(i,J), taux ) - tauyDG_v(i,J,1) = tauy_v(i,J) - depth = 0.0 - do k = 1, nz - depth = depth + CS%h_v(i,J,k) - if( (depth >= hbl_v(i,J)) .and. (kbl_v(i,J) == 0) .and. (k > (kblmin-1))) then - kbl_v(i,J) = k - hbl_v(i,J) = depth - endif - enddo + tauy_v(i,J) = forces%tauy(i,J) * Irho0 + if ( (G%mask2dCv(i,J) > 0.5) ) then + ! h to v-pts + tmp_v = max( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1))) + hbl_v(i,J) = (G%mask2dT(i,j) * hbl_h(i,J) + G%mask2dT(i,j+1) * hbl_h(i,j+1)) / tmp_v endif - enddo - enddo - - if (CS%debug) then - !### These checksum calls are missing necessary dimensional scaling factors. - call uvchksum("surface tau[xy]_[uv] ", taux_u, tauy_v, G%HI, haloshift=1, scalar_pair=.true.) - call uvchksum("ustar2", ustar2_u, ustar2_v, G%HI, haloshift=0, scalar_pair=.true.) - call uvchksum(" hbl", hbl_u , hbl_v , G%HI, haloshift=0, scalar_pair=.true.) - endif - - ! Compute downgradient stresses - do k = 1, nz - kp1 = min( k+1 , nz) - do j = js ,je - do I = Isq , Ieq - tauxDG_u(I,j,k+1) = CS%a_u(I,j,kp1) * (ui(I,j,k) - ui(I,j,kp1)) - enddo - enddo - do J = Jsq , Jeq - do i = is , ie - tauyDG_v(i,J,k+1) = CS%a_v(i,J,kp1) * (vi(i,J,k) - vi(i,J,kp1)) - enddo - enddo - enddo - - call pass_vector(tauxDG_u, tauyDG_v , G%Domain, To_All) - call pass_vector(ui,vi, G%Domain, To_All) - tauxDG_v(:,:,:) = 0. - tauyDG_u(:,:,:) = 0. - - ! Thickness weighted interpolations - do k = 1, nz - ! v to u points - do j = js , je - do I = Isq, Ieq - tauyDG_u(I,j,k) = set_v_at_u(tauyDG_v, h, G, GV, I, j, k, G%mask2dCv, OBC) - enddo - enddo - ! u to v points - do J = Jsq, Jeq - do i = is, ie - tauxDG_v(i,J,k) = set_u_at_v(tauxDG_u, h, G, GV, i, J, k, G%mask2dCu, OBC) + depth = 0. + Gat1 = 0. + do k=1, nz + ! cell center + depth = depth + 0.5* CS%h_v(i,J,k) + vE_v(i,J,k) = vi(i,J,k) - waves%Us_y(i,J,k) + if ( depth < hbl_v(i,J) ) then + sigma = depth / hbl_v(i,J) + ! cell bottom + depth = depth + 0.5* CS%h_v(i,J,k) + call cvmix_kpp_composite_Gshape(sigma,Gat1,Gsig,dGdsig) + ! nonlocal boundary-layer increment + vInc_v(i,J,k) = dt * Cemp_NL * tauy_v(i,J) * dGdsig / hbl_v(i,J) + vi(i,J,k) = vi(i,J,k) + vInc_v(i,J,k) + else + vInc_v(i,J,k) = 0.0 + endif enddo enddo enddo - if (CS%debug) then - call uvchksum(" tauyDG_u tauxDG_v",tauyDG_u,tauxDG_v, G%HI, haloshift=0, scalar_pair=.true.) - endif - ! compute angles, tau2x_[u,v], tau2w_[u,v], tau2s_[u,v], s2w_[u,v] and stress mag tau_[u,v] - omega_tau2w_u(:,:,:) = 0.0 - omega_tau2w_v(:,:,:) = 0.0 - omega_tau2s_u(:,:,:) = 0.0 - omega_tau2s_v(:,:,:) = 0.0 - tau_u(:,:,:) = 0.0 - tau_v(:,:,:) = 0.0 - - ! stress magnitude tau_[uv] & direction Omega_tau2(w,s,x)_[uv] - do j = js,je - do I = Isq,Ieq - if( (G%mask2dCu(I,j) > 0.5) ) then - ! SURFACE - tauyDG_u(I,j,1) = ustar2_u(I,j) !* cos(omega_w2x_u(I,j)) - tau_u(I,j,1) = ustar2_u(I,j) - Omega_tau2w_u(I,j,1) = 0.0 - Omega_tau2s_u(I,j,1) = 0.0 + ! Compute and store diagnostics, only during the corrector step. + if (lpost) then + call pass_vector(uE_u , vE_v , G%Domain, To_All) + call pass_vector(uInc_u, vInc_v , G%Domain, To_All) + uStk = 0.0 + vStk = 0.0 + uS0 = 0.0 + vS0 = 0.0 + + do j = js,je + do i = is,ie + if (G%mask2dT(i,j) > 0.5) then + ! u to h-pts + tmp_u = max( 1.0 ,(G%mask2dCu(i,j) + G%mask2dCu(i-1,j))) + ! v to h-pts + tmp_v = max( 1.0 ,(G%mask2dCv(i,j) + G%mask2dCv(i,j-1))) + do k = 1,nz + uE_h(i,j,k) = (G%mask2dCu(i,j) * uE_u(i,j,k) + G%mask2dCu(i-1,j) * uE_u(i-1,j,k)) / tmp_u + uInc_h(i,j,k) = (G%mask2dCu(i,j) * uInc_u(i,j,k) + G%mask2dCu(i-1,j) * uInc_u(i-1,j,k)) / tmp_u + vE_h(i,j,k) = (G%mask2dCv(i,j) * vE_v(i,j,k) + G%mask2dCv(i,j-1) * vE_v(i,j-1,k)) / tmp_v + vInc_h(i,j,k) = (G%mask2dCv(i,j) * vInc_v(i,j,k) + G%mask2dCv(i,j-1) * vInc_v(i,j-1,k)) / tmp_v + enddo + ! Wind, Stress and Shear align at surface + Omega_tau2w(i,j,1) = 0.0 + Omega_tau2s(i,j,1) = 0.0 + do k = 1,nz + kp1 = min( nz , k+1) + du = uE_h(i,j,k) - uE_h(i,j,kp1) + dv = vE_h(i,j,k) - vE_h(i,j,kp1) + omega_s2x = atan2( dv , du ) + + du = du + uInc_h(i,j,k) - uInc_h(i,j,kp1) + dv = dv + vInc_h(i,j,k) - vInc_h(i,j,kp1) + omega_tau2x = atan2( dv , du ) + + omega_tmp = omega_tau2x - forces%omega_w2x(i,j) + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2w(i,j,kp1) = omega_tmp + + omega_tmp = omega_tau2x - omega_s2x + if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi + if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi + Omega_tau2s(i,j,kp1) = omega_tmp - do k=1,nz - kp1 = MIN(k+1 , nz) - tau_u(I,j,k+1) = sqrt( tauxDG_u(I,j,k+1)*tauxDG_u(I,j,k+1) + tauyDG_u(I,j,k+1)*tauyDG_u(I,j,k+1)) - Omega_tau2x = atan2( tauyDG_u(I,j,k+1) , tauxDG_u(I,j,k+1) ) - omega_tmp = Omega_tau2x !- omega_w2x_u(I,j) - if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi - if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi - Omega_tau2w_u(I,j,k+1) = omega_tmp - Omega_tau2s_u(I,j,k+1) = 0.0 - enddo - endif - enddo - enddo - do J = Jsq, Jeq - do i = is, ie - if( (G%mask2dCv(i,J) > 0.5) ) then - ! SURFACE - tauxDG_v(i,J,1) = ustar2_v(i,J) !* sin(omega_w2x_v(i,J)) - tau_v(i,J,1) = ustar2_v(i,J) - Omega_tau2w_v(i,J,1) = 0.0 - Omega_tau2s_v(i,J,1) = 0.0 - - do k=1,nz-1 - kp1 = MIN(k+1 , nz) - tau_v(i,J,k+1) = sqrt ( tauxDG_v(i,J,k+1)*tauxDG_v(i,J,k+1) + tauyDG_v(i,J,k+1)*tauyDG_v(i,J,k+1) ) - omega_tau2x = atan2( tauyDG_v(i,J,k+1) , tauxDG_v(i,J,k+1) ) - omega_tmp = omega_tau2x !- omega_w2x_v(i,J) - if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi - if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi - Omega_tau2w_v(i,J,k+1) = omega_tmp - Omega_tau2s_v(i,J,k+1) = 0.0 - enddo - endif - enddo - enddo + enddo + endif - ! Parameterized stress orientation from the wind at interfaces (tau2x) - ! and centers (tau2x) OVERWRITE to kbl-interface above hbl - do j = js,je - do I = Isq,Ieq - if( (G%mask2dCu(I,j) > 0.5) ) then - kbld = min( (kbl_u(I,j)) , (nz-2) ) - if ( tau_u(I,j,kbld+2) > tau_u(I,j,kbld+1) ) kbld = kbld + 1 - - !### This expression is dimensionally inconsistent. - tauh = tau_u(I,j,kbld+1) + GV%H_subroundoff - ! surface boundary conditions - depth = 0. - tauNLup = 0.0 - do k=1, kbld - depth = depth + CS%h_u(I,j,k) - sigma = min( 1.0 , depth / hbl_u(i,j) ) - - ! linear stress mag - tau_MAG = (ustar2_u(I,j) * (1.-sigma) ) + (tauh * sigma ) - !### The following expressions are dimensionally inconsistent. - cos_tmp = tauxDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) - sin_tmp = tauyDG_u(I,j,k+1) / (tau_u(I,j,k+1) + GV%H_subroundoff) - - ! rotate to wind coordinates - Wind_x = ustar2_u(I,j) !* cos(omega_w2x_u(I,j)) - Wind_y = ustar2_u(I,j) !* sin(omega_w2x_u(I,j)) - tauNL_DG = (Wind_x * cos_tmp + Wind_y * sin_tmp) - tauNL_CG = (Wind_y * cos_tmp - Wind_x * sin_tmp) - omega_w2s = atan2(tauNL_CG, tauNL_DG) - omega_s2w = 0.0-omega_w2s - tauNL_CG = Cemp_CG * G_sig(sigma) * tauNL_CG - tau_MAG = max(tau_MAG, tauNL_CG) - tauNL_DG = sqrt(tau_MAG*tau_MAG - tauNL_CG*tauNL_CG) - tau_u(I,j,k+1) - - ! back to x,y coordinates - tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp) - tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp) - tauNLdn = tauNL_X - - ! nonlocal increment and update to uold - !### The following expression is dimensionally inconsistent and missing parentheses. - du = (tauNLup - tauNLdn) * (dt/CS%h_u(I,j,k) + GV%H_subroundoff) - ui(I,j,k) = uold(I,j,k) + du - uold(I,j,k) = du - tauNLup = tauNLdn - - ! diagnostics - Omega_tau2s_u(I,j,k+1) = atan2(tauNL_CG , (tau_u(I,j,k+1)+tauNL_DG)) - tau_u(I,j,k+1) = sqrt((tauxDG_u(I,j,k+1) + tauNL_X)**2 + (tauyDG_u(I,j,k+1) + tauNL_Y)**2) - omega_tau2x = atan2((tauyDG_u(I,j,k+1) + tauNL_Y), (tauxDG_u(I,j,k+1) + tauNL_X)) - omega_tau2w = omega_tau2x !- omega_w2x_u(I,j) - if (omega_tau2w >= pi ) omega_tau2w = omega_tau2w - 2.*pi - if (omega_tau2w <= (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi - Omega_tau2w_u(I,j,k+1) = omega_tau2w + ! Stokes drift + do b=1,waves%NumBands + uS0(i,j) = uS0(i,j) + waves%UStk_Hb(i,j,b) ! or forces%UStkb(i,j,b) + vS0(i,j) = vS0(i,j) + waves%VStk_Hb(i,j,b) ! or forces%VStkb(i,j,b) enddo - do k= kbld+1, nz - ui(I,j,k) = uold(I,j,k) - uold(I,j,k) = 0.0 + depth = 0.0 + do k = 1,nz + do b = 1, waves%NumBands + ! cell center + fexp = exp(-2. * waves%WaveNum_Cen(b) * (depth+0.5*h(i,j,k)) ) + uStk(i,j,k) = uStk(i,j,k) + waves%UStk_Hb(i,j,b) * fexp + vStk(i,j,k) = vStk(i,j,k) + waves%VStk_Hb(i,j,b) * fexp + enddo + ! cell bottom + depth = depth + h(i,j,k) enddo - endif + enddo enddo - enddo - ! v-point dv increment - do J = Jsq,Jeq - do i = is,ie - if( (G%mask2dCv(i,J) > 0.5) ) then - kbld = min((kbl_v(i,J)), (nz-2)) - if (tau_v(i,J,kbld+2) > tau_v(i,J,kbld+1)) kbld = kbld + 1 - tauh = tau_v(i,J,kbld+1) - - !surface boundary conditions - depth = 0. - tauNLup = 0.0 - do k=1, kbld - depth = depth + CS%h_v(i,J,k) - sigma = min(1.0, depth/ hbl_v(I,J)) - - ! linear stress - tau_MAG = (ustar2_v(i,J) * (1.-sigma)) + (tauh * sigma) - !### The following expressions are dimensionally inconsistent. - cos_tmp = tauxDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) - sin_tmp = tauyDG_v(i,J,k+1) / (tau_v(i,J,k+1) + GV%H_subroundoff) - - ! rotate into wind coordinate - Wind_x = ustar2_v(i,J) !* cos(omega_w2x_v(i,J)) - Wind_y = ustar2_v(i,J) !* sin(omega_w2x_v(i,J)) - tauNL_DG = (Wind_x * cos_tmp + Wind_y * sin_tmp) - tauNL_CG = (Wind_y * cos_tmp - Wind_x * sin_tmp) - omega_w2s = atan2(tauNL_CG , tauNL_DG) - omega_s2w = 0.0 - omega_w2s - tauNL_CG = Cemp_CG * G_sig(sigma) * tauNL_CG - tau_MAG = max( tau_MAG , tauNL_CG ) - tauNL_DG = 0.0 - tau_v(i,J,k+1) + sqrt(tau_MAG*tau_MAG - tauNL_CG*tauNL_CG) - - ! back to x,y coordinate - tauNL_X = (tauNL_DG * cos_tmp - tauNL_CG * sin_tmp) - tauNL_Y = (tauNL_DG * sin_tmp + tauNL_CG * cos_tmp) - tauNLdn = tauNL_Y - !### The following expression is dimensionally inconsistent, [L T-1] vs. [L2 H-1 T-1] on the right, - ! and it is inconsistent with the counterpart expression for du. - dv = (tauNLup - tauNLdn) * (dt/(CS%h_v(i,J,k)) ) - vi(i,J,k) = vold(i,J,k) + dv - vold(i,J,k) = dv - tauNLup = tauNLdn - - ! diagnostics - Omega_tau2s_v(i,J,k+1) = atan2(tauNL_CG, tau_v(i,J,k+1) + tauNL_DG) - tau_v(i,J,k+1) = sqrt((tauxDG_v(i,J,k+1) + tauNL_X)**2 + (tauyDG_v(i,J,k+1) + tauNL_Y)**2) - !omega_tau2x = atan2((tauyDG_v(i,J,k+1) + tauNL_Y) , (tauxDG_v(i,J,k+1) + tauNL_X)) - !omega_tau2w = omega_tau2x - omega_w2x_v(i,J) - if (omega_tau2w > pi) omega_tau2w = omega_tau2w - 2.*pi - if (omega_tau2w .le. (0.-pi) ) omega_tau2w = omega_tau2w + 2.*pi - Omega_tau2w_v(i,J,k+1) = omega_tau2w - enddo + ! post FPmix diagnostics + if (CS%id_uE_h > 0) call post_data(CS%id_uE_h , uE_h , CS%diag) + if (CS%id_vE_h > 0) call post_data(CS%id_vE_h , vE_h , CS%diag) + if (CS%id_uInc_h > 0) call post_data(CS%id_uInc_h , uInc_h , CS%diag) + if (CS%id_vInc_h > 0) call post_data(CS%id_vInc_h , vInc_h , CS%diag) + if (CS%id_FPtau2s > 0) call post_data(CS%id_FPtau2s, Omega_tau2s, CS%diag) + if (CS%id_FPtau2w > 0) call post_data(CS%id_FPtau2w, Omega_tau2w, CS%diag) + if (CS%id_uStk0 > 0) call post_data(CS%id_uStk0 , uS0 , CS%diag) + if (CS%id_vStk0 > 0) call post_data(CS%id_vStk0 , vS0 , CS%diag) + if (CS%id_uStk > 0) call post_data(CS%id_uStk , uStk , CS%diag) + if (CS%id_vStk > 0) call post_data(CS%id_vStk , vStk , CS%diag) + if (CS%id_Omega_w2x > 0) call post_data(CS%id_Omega_w2x, forces%omega_w2x, CS%diag) - do k= kbld+1, nz - vi(i,J,k) = vold(i,J,k) - vold(i,J,k) = 0.0 - enddo - endif - enddo - enddo - - if (CS%debug) then - call uvchksum("FP-tau_[uv] ", tau_u, tau_v, G%HI, haloshift=0, scalar_pair=.true.) endif - if (CS%id_tauFP_u > 0) call post_data(CS%id_tauFP_u, tau_u, CS%diag) - if (CS%id_tauFP_v > 0) call post_data(CS%id_tauFP_v, tau_v, CS%diag) - if (CS%id_FPtau2s_u > 0) call post_data(CS%id_FPtau2s_u, omega_tau2s_u, CS%diag) - if (CS%id_FPtau2s_v > 0) call post_data(CS%id_FPtau2s_v, omega_tau2s_v, CS%diag) - if (CS%id_FPtau2w_u > 0) call post_data(CS%id_FPtau2w_u, omega_tau2w_u, CS%diag) - if (CS%id_FPtau2w_v > 0) call post_data(CS%id_FPtau2w_v, omega_tau2w_v, CS%diag) - !if (CS%id_FPw2x > 0) call post_data(CS%id_FPw2x, forces%omega_w2x , CS%diag) - end subroutine vertFPmix -!> Returns the empirical shape-function given sigma [nondim] -real function G_sig(sigma) - real , intent(in) :: sigma !< Normalized boundary layer depth [nondim] - - ! local variables - real :: p1, c2, c3 !< Parameters used to fit and match empirical shape-functions [nondim] - - ! parabola - p1 = 0.287 - ! cubic function - c2 = 1.74392 - c3 = 2.58538 - G_sig = min( p1 * (1.-sigma)*(1.-sigma) , sigma * (1. + sigma * (c2*sigma - c3) ) ) -end function G_sig - !> Compute coupling coefficient associated with vertical viscosity parameterization as in Greatbatch and Lamb !! (1990), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme !! redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization, @@ -701,7 +524,7 @@ end subroutine find_coupling_coef_gl90 !! 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, Waves) + taux_bot, tauy_bot, fpmix, 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 @@ -725,6 +548,7 @@ 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 @@ -765,6 +589,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & logical :: do_i(SZIB_(G)) logical :: DoStokesMixing + logical :: lfpmix 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 @@ -802,6 +627,8 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & call MOM_error(FATAL,"Stokes Mixing called without allocated"//& "Waves Control Structure") endif + lfpmix = .false. + if ( present(fpmix) ) lfpmix = fpmix do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo @@ -814,11 +641,17 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & do j=G%jsc,G%jec do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0.0) ; enddo + ! 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 I=Isq,Ieq if (do_i(I)) u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) + enddo ; enddo ; endif + if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = u(I,j,k) enddo ; enddo ; endif @@ -976,6 +809,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (do_i(I)) u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) + enddo ; enddo ; endif + enddo ! end u-component j loop ! Now work on the meridional velocity component. @@ -991,6 +828,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (do_i(i)) v(i,j,k) = v(i,j,k) + Waves%Us_y(i,j,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do i=is,ie + if (do_i(i)) v(i,j,k) = v(i,j,k) - Waves%Us_y(i,j,k) + enddo ; enddo ; endif + if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = v(i,J,k) enddo ; enddo ; endif @@ -1119,6 +960,10 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (do_i(i)) v(i,J,k) = v(i,J,k) - Waves%Us_y(i,J,k) enddo ; enddo ; endif + if ( lfpmix ) then ; do k=1,nz ; do i=is,ie + if (do_i(i)) v(i,J,k) = v(i,J,k) + Waves%Us_y(i,J,k) + enddo ; enddo ; endif + enddo ! end of v-component J loop ! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3]. @@ -2618,7 +2463,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) + ntrunc, CS, fpmix) type(ocean_internal_state), & target, intent(in) :: MIS !< The "MOM Internal State", a set of pointers !! to the fields and accelerations that make @@ -2633,6 +2478,7 @@ 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 ! Local variables @@ -2640,6 +2486,7 @@ 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 character(len=200) :: kappa_gl90_file, inputdir, kdgl90_varname ! This include declares and sets the variable "version". # include "version_variable.h" @@ -2664,6 +2511,9 @@ 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 + ! Default, read and log parameters call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & @@ -2966,20 +2816,29 @@ 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) - CS%id_FPw2x = register_diag_field('ocean_model', 'FPw2x', diag%axesT1, Time, & - 'Wind direction from x-axis','radians') - CS%id_tauFP_u = register_diag_field('ocean_model', 'tauFP_u', diag%axesCui, Time, & - 'Stress Mag Profile (u-points)', 'm2 s-2') - CS%id_tauFP_v = register_diag_field('ocean_model', 'tauFP_v', diag%axesCvi, Time, & - 'Stress Mag Profile (v-points)', 'm2 s-2') - CS%id_FPtau2s_u = register_diag_field('ocean_model', 'FPtau2s_u', diag%axesCui, Time, & - 'stress from shear direction (u-points)', 'radians ') - CS%id_FPtau2s_v = register_diag_field('ocean_model', 'FPtau2s_v', diag%axesCvi, Time, & - 'stress from shear direction (v-points)', 'radians') - CS%id_FPtau2w_u = register_diag_field('ocean_model', 'FPtau2w_u', diag%axesCui, Time, & - 'stress from wind direction (u-points)', 'radians') - CS%id_FPtau2w_v = register_diag_field('ocean_model', 'FPtau2w_v', diag%axesCvi, Time, & - 'stress from wind direction (v-points)', 'radians') + if (lfpmix) 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, & + Time, 'y-merid Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_uInc_h = register_diag_field('ocean_model','uInc_h',CS%diag%axesTL, & + Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_vInc_h = register_diag_field('ocean_model','vInc_h',CS%diag%axesTL, & + Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_uStk = register_diag_field('ocean_model', 'uStk' , CS%diag%axesTL, & + Time, 'x-FP du increment' , 'm s-1', conversion=US%L_T_to_m_s) + CS%id_vStk = register_diag_field('ocean_model', 'vStk' , CS%diag%axesTL, & + Time, 'y-FP dv increment' , 'm s-1', conversion=US%L_T_to_m_s) + + CS%id_FPtau2s = register_diag_field('ocean_model','Omega_tau2s',CS%diag%axesTi, & + Time, 'Stress direction from shear','radians') + CS%id_FPtau2w = register_diag_field('ocean_model','Omega_tau2w',CS%diag%axesTi, & + Time, 'Stress direction from wind','radians') + CS%id_uStk0 = register_diag_field('ocean_model', 'uStk0' , diag%axesT1, & + Time, 'Zonal Surface Stokes', 'm s-1', conversion=US%L_T_to_m_s) + CS%id_vStk0 = register_diag_field('ocean_model', 'vStk0' , diag%axesT1, & + Time, 'Merid Surface Stokes', 'm s-1', conversion=US%L_T_to_m_s) + endif CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, Time, & 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 656ff5b569..3744469891 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -715,7 +715,7 @@ subroutine Update_Surface_Waves(G, GV, US, Time_present, dt, CS, forces) enddo do j=G%jsc,G%jec do i=G%isc,G%iec - !CS%Omega_w2x(i,j) = forces%omega_w2x(i,j) + CS%Omega_w2x(i,j) = forces%omega_w2x(i,j) do b=1,CS%NumBands CS%UStk_Hb(i,j,b) = forces%UStkb(i,j,b) CS%VStk_Hb(i,j,b) = forces%VStkb(i,j,b)