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 054d42a084..0ac99a3940 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -299,6 +299,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., & press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, & cfc=CS%use_CFC, hevap=CS%enthalpy_cpl) + call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) @@ -704,6 +705,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) if (CS%rigid_sea_ice) then call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed) @@ -864,6 +866,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)) enddo ; enddo call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) else ! C-grid wind stresses. diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 9fb1a6b356..df28dc0338 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -65,10 +65,13 @@ 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 +use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix 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 +use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS +use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS +use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member implicit none ; private @@ -134,6 +137,8 @@ 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] @@ -172,10 +177,12 @@ module MOM_dynamics_split_RK2 !! Euler (1) [nondim]. 0 is often used. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: debug_OBC !< If true, do debugging calls for open boundary conditions. + logical :: fpmix !< If true, applies profiles of momentum flux magnitude and direction. 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 @@ -346,6 +353,11 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_OBC ! The starting meridional velocities, which are ! 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 + ! [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(:,:) :: & @@ -369,9 +381,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1]. h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth from Cvmix real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. real :: Idt_bc ! Inverse of the baroclinic timestep [T-1 ~> s-1] - logical :: dyn_p_surf logical :: BT_cont_BT_thick ! If true, use the BT_cont_type to estimate the ! relative weightings of the layers in calculating @@ -659,10 +671,40 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (CS%debug) then call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, scale=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 vertvisc_coef(up, vp, h, forces, visc, dt_pred, G, GV, US, CS%vertvisc_CSp, & CS%OBC, VarMix) call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, & GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) + + if (CS%fpmix) then + hbl(:,:) = 0.0 + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) + if (ASSOCIATED(CS%energetic_PBL_CSp)) & + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) + call vertFPmix(up, vp, uold, vold, hbl, h, forces, & + dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) + call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, & + 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)") if (G%nonblocking_updates) then call cpu_clock_end(id_clock_vertvisc) @@ -880,9 +922,35 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! u <- u + dt d/dz visc d/dz u ! 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(I,j,k) + enddo + enddo + do J = Jsq, Jeq + do i = is, ie + vold(i,J,k) = v(i,J,k) + enddo + enddo + enddo + endif + call vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix) call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves) + + if (CS%fpmix) then + call vertFPmix(u, v, uold, vold, hbl, h, forces, dt, & + G, GV, US, CS%vertvisc_CSp, CS%OBC) + call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & + CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) + endif + if (G%nonblocking_updates) then call cpu_clock_end(id_clock_vertvisc) call start_group_pass(CS%pass_uv, G%Domain, clock=id_clock_pass) @@ -963,6 +1031,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s CS%CAu_pred_stored = .false. endif + if (CS%fpmix) then + if (CS%id_uold > 0) call post_data(CS%id_uold , uold, CS%diag) + if (CS%id_vold > 0) call post_data(CS%id_vold , vold, CS%diag) + endif ! The time-averaged free surface height has already been set by the last call to btstep. @@ -1298,6 +1370,9 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param "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, apply profiles of momentum flux magnitude and "//& + " direction", default=.false.) call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", CS%remap_aux, & "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//& "variables that are needed to reproduce across restarts, similarly to "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index a36fec3bb5..9cc4261e5f 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -67,6 +67,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 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 @@ -226,7 +227,9 @@ 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] + 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) real, pointer, dimension(:,:) :: p_surf_full => NULL() @@ -362,8 +365,8 @@ module MOM_forcing_type integer :: id_taux = -1 integer :: id_tauy = -1 integer :: id_ustar = -1 + integer :: id_omega_w2x = -1 integer :: id_tau_mag = -1 - integer :: id_psurf = -1 integer :: id_TKE_tidal = -1 integer :: id_buoy = -1 @@ -1328,6 +1331,9 @@ 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') + if (present(use_berg_fluxes)) then if (use_berg_fluxes) then handles%id_ustar_berg = register_diag_field('ocean_model', 'ustar_berg', diag%axesT1, Time, & @@ -2165,6 +2171,11 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) 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) @@ -2302,6 +2313,11 @@ subroutine copy_back_forcing_fields(fluxes, forces, G) 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) @@ -2950,6 +2966,9 @@ 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_ustar_berg > 0) .and. associated(fluxes%ustar_berg)) & call post_data(handles%id_ustar_berg, fluxes%ustar_berg, diag) @@ -3275,6 +3294,7 @@ end subroutine myAlloc subroutine deallocate_forcing_type(fluxes) type(forcing), intent(inout) :: fluxes !< Forcing fields structure + if (associated(fluxes%omega_w2x)) deallocate(fluxes%omega_w2x) if (associated(fluxes%ustar)) deallocate(fluxes%ustar) if (associated(fluxes%ustar_gustless)) deallocate(fluxes%ustar_gustless) if (associated(fluxes%tau_mag)) deallocate(fluxes%tau_mag) @@ -3334,6 +3354,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%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_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 47d4dffef6..8200201ee3 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -36,7 +36,8 @@ module MOM_set_visc #include public set_viscous_BBL, set_viscous_ML, set_visc_init, set_visc_end -public set_visc_register_restarts, remap_vertvisc_aux_vars +public set_visc_register_restarts, set_u_at_v, set_v_at_u +public remap_vertvisc_aux_vars ! 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 diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 80fff62f21..f513f50158 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -3,6 +3,7 @@ module MOM_vert_friction ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domains, only : pass_var, To_All, Omit_corners +use MOM_domains, only : pass_vector, Scalar_Pair use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : post_product_u, post_product_sum_u use MOM_diag_mediator, only : post_product_v, post_product_sum_v @@ -27,6 +28,7 @@ module MOM_vert_friction use MOM_variables, only : ocean_internal_state use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS +use MOM_set_visc, only : set_v_at_u, set_u_at_v use MOM_lateral_mixing_coeffs, only : VarMix_CS implicit none ; private @@ -35,6 +37,7 @@ module MOM_vert_friction public vertvisc, vertvisc_remnant, vertvisc_coef public vertvisc_limit_vel, vertvisc_init, vertvisc_end public updateCFLtruncationValue +public vertFPmix ! 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 @@ -166,6 +169,9 @@ 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_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 @@ -184,6 +190,381 @@ 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) + 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] + 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 + + ! 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, omega_tmp !< constants and dummy variables + real :: du, dv, depth, sigma, Wind_x, Wind_y !< intermediate variables + real :: taux, tauy, tauxDG, tauyDG, tauxDGup, tauyDGup, ustar2, tauh !< intermediate variables + real :: tauNLup, tauNLdn, tauNL_CG, tauNL_DG, tauNL_X, tauNL_Y, tau_MAG !< intermediate variables + real :: omega_w2s, omega_tau2s, omega_s2x, omega_tau2x, omega_tau2w, omega_s2w !< intermediate angles + integer :: kblmin, kbld, kp1, k, nz !< 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. + + 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_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 + 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 + endif + enddo + enddo + 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 + endif + enddo + enddo + + if (CS%debug) then + 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) + 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 + + 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 + + ! 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 + + 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 ) + 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 + 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 + enddo + do k= kbld+1, nz + ui(I,j,k) = uold(I,j,k) + uold(I,j,k) = 0.0 + enddo + endif + 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) + 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 + 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 + + 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. +real function G_sig(sigma) + real , intent(in) :: sigma !< non-dimensional normalized boundary layer depth [m] + + ! local variables + real :: p1, c2, c3 !< parameters used to fit and match empirycal shape-functions. + + ! 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, @@ -2487,6 +2868,21 @@ 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=GV%H_to_MKS) + 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') + 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) if (CS%id_du_dt_visc > 0) call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)