From 77084ecc92bb7a26fb4d92a57dd889418b8246d0 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Fri, 21 May 2021 12:29:08 -0400 Subject: [PATCH 1/4] Adding Stokes forces implementation and diags. - Add Craik-Leibovich terms, vorticity and dynamic pressure components. - Add exploratory Stokes ddt version in MOM_wave_interface and MOM.F90 - Add tendency diagnostics for Stokes terms w/ control and option for diagnostic only Stokes force computation in MOM_CoriolisAdv. --- src/core/MOM.F90 | 15 ++- src/core/MOM_CoriolisAdv.F90 | 97 +++++++++++++- src/core/MOM_dynamics_split_RK2.F90 | 58 +++++++- src/user/MOM_wave_interface.F90 | 199 +++++++++++++++++++++++++++- 4 files changed, 354 insertions(+), 15 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4659b685e5..dddfea228d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -646,12 +646,21 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (CS%UseWaves) then ! Update wave information, which is presently kept static over each call to step_mom call enable_averages(time_interval, Time_start + real_to_time(US%T_to_s*time_interval), CS%diag) - call Update_Stokes_Drift(G, GV, US, Waves, h, forces%ustar) + call Update_Stokes_Drift(G, GV, US, Waves, h, forces%ustar, dt) + if (Waves%Stokes_DDT) then + u(:,:,:) = u(:,:,:) + Waves%ddt_us_x(:,:,:)*dt + v(:,:,:) = v(:,:,:) + Waves%ddt_us_y(:,:,:)*dt + endif call disable_averaging(CS%diag) endif else ! not do_dyn. - if (CS%UseWaves) & ! Diagnostics are not enabled in this call. - call Update_Stokes_Drift(G, GV, US, Waves, h, fluxes%ustar) + if (CS%UseWaves) then ! Diagnostics are not enabled in this call. + call Update_Stokes_Drift(G, GV, US, Waves, h, fluxes%ustar, dt) + if (Waves%Stokes_DDT) then + u(:,:,:) = u(:,:,:) + Waves%ddt_us_x(:,:,:)*dt + v(:,:,:) = v(:,:,:) + Waves%ddt_us_y(:,:,:)*dt + endif + endif endif if (CS%debug) then diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 231b6ed058..a8ab9ed5bc 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -16,6 +16,7 @@ module MOM_CoriolisAdv use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : accel_diag_ptrs use MOM_verticalGrid, only : verticalGrid_type +use MOM_wave_interface, only : wave_parameters_CS implicit none ; private @@ -82,6 +83,7 @@ module MOM_CoriolisAdv integer :: id_h_gKEu = -1, id_h_gKEv = -1 integer :: id_h_rvxu = -1, id_h_rvxv = -1 integer :: id_intz_rvxu_2d = -1, id_intz_rvxv_2d = -1 + integer :: id_CAuS = -1, id_CAvS = -1 !>@} end type CoriolisAdv_CS @@ -117,7 +119,7 @@ module MOM_CoriolisAdv contains !> Calculates the Coriolis and momentum advection contributions to the acceleration. -subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) +subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, Waves) type(ocean_grid_type), intent(in) :: G !< Ocen grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1] @@ -135,10 +137,12 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) type(accel_diag_ptrs), intent(inout) :: AD !< Storage for acceleration diagnostics type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(CoriolisAdv_CS), pointer :: CS !< Control structure for MOM_CoriolisAdv - + type(Wave_parameters_CS), optional, pointer :: Waves !< An optional pointer to Stokes drift CS + ! Local variables real, dimension(SZIB_(G),SZJB_(G)) :: & q, & ! Layer potential vorticity [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1]. + qS, & ! Layer Stokes vorticity [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1]. Ih_q, & ! The inverse of thickness interpolated to q points [H-1 ~> m-1 or m2 kg-1]. Area_q ! The sum of the ocean areas at the 4 adjacent thickness points [L2 ~> m2]. @@ -172,8 +176,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) ! discretization [H-1 s-1 ~> m-1 s-1 or m2 kg-1 s-1]. real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! Contributions to the circulation around q-points [L2 T-1 ~> m2 s-1] + dvSdx, duSdy, & ! idem. for Stokes drift [L2 T-1 ~> m2 s-1] rel_vort, & ! Relative vorticity at q-points [T-1 ~> s-1]. abs_vort, & ! Absolute vorticity at q-points [T-1 ~> s-1]. + stk_vort, & ! Stokes vorticity at q-points [T-1 ~> s-1]. q2, & ! Relative vorticity over thickness [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1]. max_fvq, & ! The maximum of the adjacent values of (-u) times absolute vorticity [L T-2 ~> m s-2]. min_fvq, & ! The minimum of the adjacent values of (-u) times absolute vorticity [L T-2 ~> m s-2]. @@ -182,6 +188,8 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: & PV, & ! A diagnostic array of the potential vorticities [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1]. RV ! A diagnostic array of the relative vorticities [T-1 ~> s-1]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & CAuS ! + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & CAvS ! real :: fv1, fv2, fv3, fv4 ! (f+rv)*v [L T-2 ~> m s-2]. real :: fu1, fu2, fu3, fu4 ! -(f+rv)*u [L T-2 ~> m s-2]. real :: max_fv, max_fu ! The maximum or minimum of the neighboring Coriolis @@ -220,7 +228,8 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) real :: UHeff, VHeff ! More temporary variables [H L2 T-1 ~> m3 s-1 or kg s-1]. real :: QUHeff,QVHeff ! More temporary variables [H L2 T-1 s-1 ~> m3 s-2 or kg s-2]. integer :: i, j, k, n, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz - + logical :: Stokes_VF, Passive_Stokes_VF + ! Diagnostics for fractional thickness-weighted terms real, allocatable, dimension(:,:) :: & hf_gKEu_2d, hf_gKEv_2d, & ! Depth sum of hf_gKEu, hf_gKEv [L T-2 ~> m s-2]. @@ -283,6 +292,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) Area_q(i,j) = (Area_h(i,j) + Area_h(i+1,j+1)) + & (Area_h(i+1,j) + Area_h(i,j+1)) enddo ; enddo + + Stokes_VF = present(Waves) + if (Stokes_VF) Stokes_VF = associated(Waves) + if (Stokes_VF) Stokes_VF = Waves%Stokes_VF !$OMP parallel do default(private) shared(u,v,h,uh,vh,CAu,CAv,G,GV,CS,AD,Area_h,Area_q,& !$OMP RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,vol_neglect,h_tiny,OBC,eps_vel) @@ -293,6 +306,34 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) ! vorticity is second order accurate everywhere with free slip b.c.s, ! but only first order accurate at boundaries with no slip b.c.s. ! First calculate the contributions to the circulation around the q-point. + if (Stokes_VF) then + if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + dvSdx(I,J) = ((-Waves%us_y(i+1,J,k))*G%dyCv(i+1,J) - & + (-Waves%us_y(i,J,k))*G%dyCv(i,J)) + duSdy(I,J) = ((-Waves%us_x(I,j+1,k))*G%dxCu(I,j+1) - & + (-Waves%us_x(I,j,k))*G%dxCu(I,j)) + enddo; enddo + endif + if (Passive_Stokes_VF) then + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + dvdx(I,J) = (v(i+1,J,k)*G%dyCv(i+1,J) - v(i,J,k)*G%dyCv(i,J)) + dudy(I,J) = (u(I,j+1,k)*G%dxCu(I,j+1) - u(I,j,k)*G%dxCu(I,j)) + enddo; enddo + else + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + dvdx(I,J) = ((v(i+1,J,k)-Waves%us_y(i+1,J,k))*G%dyCv(i+1,J) - & + (v(i,J,k)-Waves%us_y(i,J,k))*G%dyCv(i,J)) + dudy(I,J) = ((u(I,j+1,k)-Waves%us_x(I,j+1,k))*G%dxCu(I,j+1) - & + (u(I,j,k)-Waves%us_x(I,j,k))*G%dxCu(I,j)) + enddo; enddo + endif + else + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + dvdx(I,J) = (v(i+1,J,k)*G%dyCv(i+1,J) - v(i,J,k)*G%dyCv(i,J)) + dudy(I,J) = (u(I,j+1,k)*G%dxCu(I,j+1) - u(I,j,k)*G%dxCu(I,j)) + enddo; enddo + endif do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 dvdx(I,J) = (v(i+1,J,k)*G%dyCv(i+1,J) - v(i,J,k)*G%dyCv(i,J)) dudy(I,J) = (u(I,j+1,k)*G%dxCu(I,j+1) - u(I,j,k)*G%dxCu(I,j)) @@ -440,11 +481,21 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) if (CS%no_slip) then do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 rel_vort(I,J) = (2.0 - G%mask2dBu(I,J)) * (dvdx(I,J) - dudy(I,J)) * G%IareaBu(I,J) - enddo ; enddo + enddo; enddo + if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + stk_vort(I,J) = (2.0 - G%mask2dBu(I,J)) * (dvSdx(I,J) - duSdy(I,J)) * G%IareaBu(I,J) + enddo; enddo + endif else do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 rel_vort(I,J) = G%mask2dBu(I,J) * (dvdx(I,J) - dudy(I,J)) * G%IareaBu(I,J) - enddo ; enddo + enddo; enddo + if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + stk_vort(I,J) = (2.0 - G%mask2dBu(I,J)) * (dvSdx(I,J) - duSdy(I,J)) * G%IareaBu(I,J) + enddo; enddo + endif endif do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 @@ -455,7 +506,13 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) hArea_q = (hArea_u(I,j) + hArea_u(I,j+1)) + (hArea_v(i,J) + hArea_v(i+1,J)) Ih_q(I,J) = Area_q(I,J) / (hArea_q + vol_neglect) q(I,J) = abs_vort(I,J) * Ih_q(I,J) - enddo ; enddo + enddo; enddo + + if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 + qS(I,J) = stk_vort(I,J) * Ih_q(I,J) + enddo; enddo + endif if (CS%id_rv > 0) then do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 @@ -679,6 +736,15 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) (ep_u(i,j)*uh(I-1,j,k) - ep_u(i+1,j)*uh(I+1,j,k)) * G%IdxCu(I,j) enddo ; enddo ; endif + if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then + ! Computing the diagnostic Stokes contribution to CAu + do j=js,je ; do I=Isq,Ieq + CAuS(I,j,k) = 0.25 * & + (qS(I,J) * (vh(i+1,J,k) + vh(i,J,k)) + & + qS(I,J-1) * (vh(i,J-1,k) + vh(i+1,J-1,k))) * G%IdxCu(I,j) + enddo ; enddo + endif + if (CS%bound_Coriolis) then do j=js,je ; do I=Isq,Ieq fv1 = abs_vort(I,J) * v(i+1,J,k) @@ -792,6 +858,15 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) (ep_v(i,j)*vh(i,J-1,k) - ep_v(i,j+1)*vh(i,J+1,k)) * G%IdyCv(i,J) enddo ; enddo ; endif + if (CS%id_CAuS>0 .or. CS%id_CAvS>0) then + ! Computing the diagnostic Stokes contribution to CAv + do J=Jsq,Jeq ; do i=is,ie + CAvS(I,j,k) = 0.25 * & + (qS(I,J) * (uh(I,j+1,k) + uh(I,j,k)) + & + qS(I,J-1) * (uh(I-1,j,k) + uh(I-1,j+1,k))) * G%IdyCv(i,J) + enddo; enddo + endif + if (CS%bound_Coriolis) then do J=Jsq,Jeq ; do i=is,ie fu1 = -abs_vort(I,J) * u(I,j+1,k) @@ -868,6 +943,8 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) if (CS%id_gKEv>0) call post_data(CS%id_gKEv, AD%gradKEv, CS%diag) if (CS%id_rvxu > 0) call post_data(CS%id_rvxu, AD%rv_x_u, CS%diag) if (CS%id_rvxv > 0) call post_data(CS%id_rvxv, AD%rv_x_v, CS%diag) + if (CS%id_CAuS > 0) call post_data(CS%id_CAuS, CAuS, CS%diag) + if (CS%id_CAvS > 0) call post_data(CS%id_CAvS, CAvS, CS%diag) ! Diagnostics for terms multiplied by fractional thicknesses @@ -1274,6 +1351,14 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) 'Zonal Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_rvxv > 0) call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) + CS%id_CAuS = register_diag_field('ocean_model', 'CAuS', diag%axesCuL, Time, & + 'Zonal Acceleration from Stokes Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2) + ! add to AD + + CS%id_CAvS = register_diag_field('ocean_model', 'CAvS', diag%axesCvL, Time, & + 'Meridional Acceleration from Stokes Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2) + ! add to AD + !CS%id_hf_gKEu = register_diag_field('ocean_model', 'hf_gKEu', diag%axesCuL, Time, & ! 'Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', & ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 07a302b2b0..34539dbea7 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -60,7 +60,7 @@ module MOM_dynamics_split_RK2 use MOM_vert_friction, only : updateCFLtruncationValue 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 +use MOM_wave_interface, only: wave_parameters_CS, Stokes_PGF implicit none ; private @@ -71,11 +71,13 @@ module MOM_dynamics_split_RK2 real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: & CAu, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2] PFu, & !< PFu = -dM/dx [L T-2 ~> m s-2] + PFu_Stokes, & !< PFu_Stokes = -d/dx int_r (u_L*duS/dr) [L T-2 ~> m s-2] diffu !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2] real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: & CAv, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2] PFv, & !< PFv = -dM/dy [L T-2 ~> m s-2] + PFv_Stokes, & !< PFv_Stokes = -d/dy int_r (v_L*dvS/dr) [L T-2 ~> m s-2] diffv !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: visc_rem_u @@ -360,6 +362,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s logical :: BT_cont_BT_thick ! If true, use the BT_cont_type to estimate the ! relative weightings of the layers in calculating ! the barotropic accelerations. + logical :: Use_Stokes_PGF ! If true, add Stokes PGF to hydrostatic PGF !---For group halo pass logical :: showCallTree, sym @@ -454,6 +457,30 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s eta_PF_start(i,j) = CS%eta_PF(i,j) - pres_to_eta * (p_surf_begin(i,j) - p_surf_end(i,j)) enddo ; enddo endif + ! Stokes shear force contribution to pressure gradient + Use_Stokes_PGF = present(Waves) + if (Use_Stokes_PGF) then + Use_Stokes_PGF = associated(Waves) + if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF + if (Use_Stokes_PGF) then + call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + ! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv + ! will therefore report the sum total PGF and we avoid other + ! modifications in the code. The PFu_Stokes can be output within the waves routines. + if (.not.Waves%Passive_Stokes_PGF) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + CS%PFu(I,j,k) = CS%PFu(I,j,k) + CS%PFu_Stokes(I,j,k) + enddo ; enddo + enddo + do k=1,nz + do J=Jsq,Jeq ; do i=is,ie + CS%PFv(i,J,k) = CS%PFv(i,J,k) + CS%PFv_Stokes(i,J,k) + enddo ; enddo + enddo + endif + endif + endif call cpu_clock_end(id_clock_pres) call disable_averaging(CS%diag) if (showCallTree) call callTree_wayPoint("done with PressureForce (step_MOM_dyn_split_RK2)") @@ -470,7 +497,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av call cpu_clock_begin(id_clock_Cor) call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, Gv, US, CS%CoriolisAdv_CSp) + G, Gv, US, CS%CoriolisAdv_CSp, Waves=Waves) call cpu_clock_end(id_clock_Cor) if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)") @@ -692,6 +719,27 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s call cpu_clock_begin(id_clock_pres) call PressureForce(hp, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, & CS%ALE_CSp, p_surf, CS%pbce, CS%eta_PF) + ! Stokes shear force contribution to pressure gradient + Use_Stokes_PGF = present(Waves) + if (Use_Stokes_PGF) then + Use_Stokes_PGF = associated(Waves) + if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF + if (Use_Stokes_PGF) then + call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + if (.not.Waves%Passive_Stokes_PGF) then + do k=1,nz + do j=js,je ; do I=Isq,Ieq + CS%PFu(I,j,k) = CS%PFu(I,j,k) + CS%PFu_Stokes(I,j,k) + enddo ; enddo + enddo + do k=1,nz + do J=Jsq,Jeq ; do i=is,ie + CS%PFv(i,J,k) = CS%PFv(i,J,k) + CS%PFv_Stokes(i,J,k) + enddo ; enddo + enddo + endif + endif + endif call cpu_clock_end(id_clock_pres) if (showCallTree) call callTree_wayPoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)") endif @@ -726,7 +774,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av call cpu_clock_begin(id_clock_Cor) call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, GV, US, CS%CoriolisAdv_CSp) + G, GV, US, CS%CoriolisAdv_CSp, Waves=Waves) call cpu_clock_end(id_clock_Cor) if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)") @@ -1317,7 +1365,9 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param ALLOC_(CS%u_accel_bt(IsdB:IedB,jsd:jed,nz)) ; CS%u_accel_bt(:,:,:) = 0.0 ALLOC_(CS%v_accel_bt(isd:ied,JsdB:JedB,nz)) ; CS%v_accel_bt(:,:,:) = 0.0 - + ALLOC_(CS%PFu_Stokes(IsdB:IedB,jsd:jed,nz)) ; CS%PFu_Stokes(:,:,:) = 0.0 + ALLOC_(CS%PFv_Stokes(isd:ied,JsdB:JedB,nz)) ; CS%PFv_Stokes(:,:,:) = 0.0 + MIS%diffu => CS%diffu MIS%diffv => CS%diffv MIS%PFu => CS%PFu diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index b9563f9369..fb188884e8 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -31,6 +31,7 @@ module MOM_wave_interface ! called in step_mom. public get_Langmuir_Number ! Public interface to compute Langmuir number called from ! ePBL or KPP routines. +public Stokes_PGF ! Public interface to compute Stokes-shear modifications to pressure gradient force public StokesMixing ! NOT READY - Public interface to add down-Stokes gradient ! momentum mixing (e.g. the approach of Harcourt 2013/2015) public CoriolisStokes ! NOT READY - Public interface to add Coriolis-Stokes acceleration @@ -50,6 +51,13 @@ module MOM_wave_interface ! Main surface wave options and publicly visible variables logical, public :: UseWaves = .false. !< Flag to enable surface gravity wave feature + logical, public :: Stokes_VF !< Developmental: + !! True if Stokes vortex force is used + logical, public :: Stokes_PGF !< Developmental: + !! True if Stokes shear pressure Gradient force is used + logical, public :: Passive_Stokes_PGF !< Keeps Stokes_PGF on, but doesn't affect dynamics + logical, public :: Stokes_DDT !< Developmental: + !! True if Stokes d/dt is used real, allocatable, dimension(:,:,:), public :: & Us_x !< 3d zonal Stokes drift profile [L T-1 ~> m s-1] !! Horizontal -> U points @@ -58,9 +66,18 @@ module MOM_wave_interface Us_y !< 3d meridional Stokes drift profile [L T-1 ~> m s-1] !! Horizontal -> V points !! Vertical -> Mid-points + real, allocatable, dimension(:,:,:), public :: & + ddt_Us_x !< 3d zonal Stokes drift profile [m s-1] + !! Horizontal -> U points + !! Vertical -> Mid-points + real, allocatable, dimension(:,:,:), public :: & + ddt_Us_y !< 3d meridional Stokes drift profile [m s-1] + !! Horizontal -> V points + !! Vertical -> Mid-points real, allocatable, dimension(:,:,:), public :: & KvS !< Viscosity for Stokes Drift shear [Z2 T-1 ~> m2 s-1] + integer, public :: id_PFu_Stokes = -1 , id_PFv_Stokes = -1 ! The remainder of this control structure is private logical :: LagrangianMixing !< This feature is in development and not ready !! True if Stokes drift is present and mixing @@ -131,6 +148,7 @@ module MOM_wave_interface !>@{ Diagnostic handles integer :: id_surfacestokes_x = -1 , id_surfacestokes_y = -1 integer :: id_3dstokes_x = -1 , id_3dstokes_y = -1 + integer :: id_ddt_3dstokes_x = -1 , id_ddt_3dstokes_y = -1 integer :: id_La_turb = -1 !>@} @@ -263,6 +281,18 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) ! Force Code Intervention call MOM_error(FATAL,"Should you be enabling Coriolis-Stokes? Code not ready.") endif + call get_param(param_file, mdl, "STOKES_VF", CS%Stokes_VF, & + "Flag to use Stokes vortex force", units="", & + Default=.false.) + call get_param(param_file, mdl, "STOKES_PGF", CS%Stokes_PGF, & + "Flag to use Stokes pressure gradient force", units="", & + Default=.false.) + call get_param(param_file, mdl, "PASSIVE_STOKES_PGF", CS%Passive_Stokes_PGF, & + "Flag to make Stokes pressure gradient force diagnostic only.", units="", & + Default=.false.) + call get_param(param_file, mdl, "STOKES_DDT", CS%Stokes_DDT, & + "Flag to use Stokes d/dt", units="", & + Default=.false.) CS%g_Earth = US%L_to_Z**2*GV%g_Earth ! Get Wave Method and write to integer WaveMethod @@ -395,6 +425,12 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) CS%Us_x(:,:,:) = 0.0 allocate(CS%Us_y(G%isd:G%Ied,G%jsdB:G%jedB,GV%ke)) CS%Us_y(:,:,:) = 0.0 + if (CS%Stokes_DDT) then + allocate(CS%ddt_Us_x(G%isdB:G%IedB,G%jsd:G%jed,G%ke)) + CS%ddt_Us_x(:,:,:) = 0.0 + allocate(CS%ddt_Us_y(G%isd:G%Ied,G%jsdB:G%jedB,G%ke)) + CS%ddt_Us_y(:,:,:) = 0.0 + endif ! b. Surface Values allocate(CS%US0_x(G%isdB:G%iedB,G%jsd:G%jed)) CS%US0_x(:,:) = 0.0 @@ -420,6 +456,14 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) CS%diag%axesCvL,Time,'3d Stokes drift (y)', 'm s-1', conversion=US%L_T_to_m_s) CS%id_3dstokes_x = register_diag_field('ocean_model','3d_stokes_x', & CS%diag%axesCuL,Time,'3d Stokes drift (x)', 'm s-1', conversion=US%L_T_to_m_s) + CS%id_ddt_3dstokes_y = register_diag_field('ocean_model','dvdt_Stokes', & + CS%diag%axesCvL,Time,'d/dt Stokes drift (meridional)','m s-2') + CS%id_ddt_3dstokes_x = register_diag_field('ocean_model','dudt_Stokes', & + CS%diag%axesCuL,Time,'d/dt Stokes drift (zonal)','m s-2') + CS%id_PFv_Stokes = register_diag_field('ocean_model','PFv_Stokes', & + CS%diag%axesCvL,Time,'PF from Stokes drift (meridional)','m s-2')!Needs conversion + CS%id_PFu_Stokes = register_diag_field('ocean_model','PFu_Stokes', & + CS%diag%axesCuL,Time,'PF from Stokes drift (zonal)','m s-2')!Needs conversion CS%id_La_turb = register_diag_field('ocean_model','La_turbulent',& CS%diag%axesT1,Time,'Surface (turbulent) Langmuir number','nondim') @@ -524,7 +568,7 @@ end subroutine Update_Surface_Waves !> Constructs the Stokes Drift profile on the model grid based on !! desired coupling options -subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) +subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt) type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure type(ocean_grid_type), intent(inout) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -533,6 +577,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) intent(in) :: h !< Thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: ustar !< Wind friction velocity [Z T-1 ~> m s-1]. + real, intent(in) :: dt !< Time-step for computing Stokes-tendency ! Local Variables real :: Top, MidPoint, Bottom ! Positions within the layer [Z ~> m] real :: one_cm ! One centimeter in the units of wavelengths [Z ~> m] @@ -544,10 +589,15 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) real :: UStokes ! A Stokes drift velocity [L T-1 ~> m s-1] real :: La ! The local Langmuir number [nondim] integer :: ii, jj, kk, b, iim1, jjm1 - + real :: idt ! 1 divided by the time step + one_cm = 0.01*US%m_to_Z min_level_thick_avg = 1.e-3*US%m_to_Z + idt = 1.0/dt + CS%ddt_us_x(:,:,:) = CS%US_x(:,:,:) + CS%ddt_us_y(:,:,:) = CS%US_y(:,:,:) + ! 1. If Test Profile Option is chosen ! Computing mid-point value from surface value and decay wavelength if (WaveMethod==TESTPROF) then @@ -772,6 +822,9 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo enddo + CS%ddt_us_x(:,:,:) = (CS%US_x(:,:,:) - CS%ddt_us_x(:,:,:)) * idt + CS%ddt_us_y(:,:,:) = (CS%US_y(:,:,:) - CS%ddt_us_y(:,:,:)) * idt + ! Output any desired quantities if (CS%id_surfacestokes_y>0) & call post_data(CS%id_surfacestokes_y, CS%us0_y, CS%diag) @@ -781,6 +834,10 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) call post_data(CS%id_3dstokes_y, CS%us_y, CS%diag) if (CS%id_3dstokes_x>0) & call post_data(CS%id_3dstokes_x, CS%us_x, CS%diag) + if (CS%id_ddt_3dstokes_x>0) & + call post_data(CS%id_ddt_3dstokes_x, CS%ddt_us_x, CS%diag) + if (CS%id_ddt_3dstokes_y>0) & + call post_data(CS%id_ddt_3dstokes_y, CS%ddt_us_y, CS%diag) if (CS%id_La_turb>0) & call post_data(CS%id_La_turb, CS%La_turb, CS%diag) @@ -1323,6 +1380,144 @@ subroutine StokesMixing(G, GV, dt, h, u, v, Waves ) end subroutine StokesMixing +!> Computes tendency due to Stokes pressure gradient force +subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) + type(ocean_grid_type), & + intent(in) :: G !< Ocean grid + type(verticalGrid_type), & + intent(in) :: GV !< Ocean vertical grid + real, dimension(SZI_(G),SZJ_(G),SZK_(G)),& + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & + intent(in) :: u !< Velocity i-component [m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(in) :: v !< Velocity j-component [m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & + intent(out) :: PFu_Stokes !< PGF Stokes-shear i-component [L T-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(out) :: PFv_Stokes !< PGF Stokes-shear j-component [m s-1] + type(Wave_parameters_CS), & + pointer :: CS !< Surface wave related control structure. + + ! Local variables + real :: P_Stokes_l, P_Stokes_r ! Contribution of Stokes shear to pressure (left/right index) [L2 T-2 ~> m2 s-2] + real :: u_l, u_r, v_l, v_r ! Velocity components + real :: dUs_dz_l, dUs_dz_r ! Vertical derivative of zonal Stokes drift (left/right index) [T-1 ~> s-1] + real :: dVs_dz_l, dVs_dz_r ! Vertical derivative of meridional Stokes drift (left/right index) [T-1 ~> s-1] + real :: z_top_l, z_top_r ! The height of the top of the cell (left/right index) [Z ~> m]. + real :: z_mid_l, z_mid_r ! The height of the middle of the cell (left/right index) [Z ~> m]. + real :: h_l, h_r ! The thickness of the cell (left/right index) [Z ~> m]. + real :: wavenum,TwoKexpL, TwoKexpR !TMP DELETE THIS + integer :: i,j,k + + ! Comput the Stokes contribution to the pressure gradient force + + PFu_Stokes(:,:,:) = 0.0 + PFv_Stokes(:,:,:) = 0.0 + + wavenum = (2.*3.14)/50. + do j = G%jsc, G%jec ; do I = G%iscB, G%iecB + if (G%mask2dCu(I,j)>0.5) then + z_top_l = 0.0 + z_top_r = 0.0 + P_Stokes_l = 0.0 + P_Stokes_r = 0.0 + do k = 1, G%ke + h_l = h(i,j,k) + h_r = h(i+1,j,k) + z_mid_l = z_top_l - 0.5*h_l + z_mid_r = z_top_r - 0.5*h_r + TwoKexpL = (2.*wavenum)*exp(2*wavenum*z_mid_l) + TwoKexpR = (2.*wavenum)*exp(2*wavenum*z_mid_r) + !UL -> I-1 & I, j + !UR -> I & I+1, j + !VL -> i, J-1 & J + !VR -> i+1, J-1 & J + dUs_dz_l = TwoKexpL*0.5 * & + (CS%Us0_x(I-1,j)*G%mask2dCu(I-1,j) + & + CS%Us0_x(I,j)*G%mask2dCu(I,j)) + dUs_dz_r = TwoKexpR*0.5 * & + (CS%Us0_x(I,j)*G%mask2dCu(I,j) + & + CS%Us0_x(I+1,j)*G%mask2dCu(I+1,j)) + dVs_dz_l = TwoKexpL*0.5 * & + (CS%Us0_y(i,J-1)*G%mask2dCv(i,J-1) + & + CS%Us0_y(i,J)*G%mask2dCv(i,J)) + dVs_dz_r = TwoKexpR*0.5 * & + (CS%Us0_y(i+1,J-1)*G%mask2dCv(i+1,J-1) + & + CS%Us0_y(i+1,J)*G%mask2dCv(i+1,J)) + u_l = 0.5*(u(I-1,j,k)*G%mask2dCu(I-1,j) + & + u(I,j,k)*G%mask2dCu(I,j)) + u_r = 0.5*(u(I,j,k)*G%mask2dCu(I,j) + & + u(I+1,j,k)*G%mask2dCu(I+1,j)) + v_l = 0.5*(v(i,J-1,k)*G%mask2dCv(i,J-1) + & + v(i,J,k)*G%mask2dCv(i,J)) + v_r = 0.5*(v(i+1,J-1,k)*G%mask2dCv(i+1,J-1) + & + v(i+1,J,k)*G%mask2dCv(i+1,J)) + if (G%mask2dT(i,j)>0.5) & + P_Stokes_l = P_Stokes_l + h_l*(dUs_dz_l*u_l+dVs_dz_l*v_l) + if (G%mask2dT(i+1,j)>0.5) & + P_Stokes_r = P_Stokes_r + h_r*(dUs_dz_r*u_r+dVs_dz_r*v_r) + PFu_Stokes(I,j,k) = (P_Stokes_r - P_Stokes_l)*G%IdxCu(I,j) + z_top_l = z_top_l - h_l + z_top_r = z_top_r - h_r + enddo + endif + enddo ; enddo + do J = G%jscB, G%jecB ; do i = G%isc, G%iec + if (G%mask2dCv(i,J)>0.5) then + z_top_l = 0.0 + z_top_r = 0.0 + P_Stokes_l = 0.0 + P_Stokes_r = 0.0 + do k = 1, G%ke + h_l = h(i,j,k) + h_r = h(i,j+1,k) + z_mid_l = z_top_l - 0.5*h_l + z_mid_r = z_top_r - 0.5*h_r + TwoKexpL = (2.*wavenum)*exp(2*wavenum*z_mid_l) + TwoKexpR = (2.*wavenum)*exp(2*wavenum*z_mid_r) + !UL -> I-1 & I, j + !UR -> I-1 & I, j+1 + !VL -> i, J & J-1 + !VR -> i, J & J+1 + dUs_dz_l = TwoKexpL*0.5 * & + (CS%Us0_x(I-1,j)*G%mask2dCu(I-1,j) + & + CS%Us0_x(I,j)*G%mask2dCu(I,j)) + dUs_dz_r = TwoKexpR*0.5 * & + (CS%Us0_x(I-1,j+1)*G%mask2dCu(I-1,j+1) + & + CS%Us0_x(I,j+1)*G%mask2dCu(I,j+1)) + dVs_dz_l = TwoKexpL*0.5 * & + (CS%Us0_y(i,J-1)*G%mask2dCv(i,J-1) + & + CS%Us0_y(i,J)*G%mask2dCv(i,J)) + dVs_dz_r = TwoKexpR*0.5 * & + (CS%Us0_y(i,J)*G%mask2dCv(i,J) + & + CS%Us0_y(i,J+1)*G%mask2dCv(i,J+1)) + u_l = 0.5*(u(I-1,j,k)*G%mask2dCu(I-1,j) + & + u(I,j,k)*G%mask2dCu(I,j)) + u_r = 0.5*(u(I-1,j+1,k)*G%mask2dCu(I-1,j+1) + & + u(I,j+1,k)*G%mask2dCu(I,j+1)) + v_l = 0.5*(v(i,J-1,k)*G%mask2dCv(i,J-1) + & + v(i,J,k)*G%mask2dCv(i,J)) + v_r = 0.5*(v(i,J,k)*G%mask2dCv(i,J) + & + v(i,J+1,k)*G%mask2dCv(i,J+1)) + if (G%mask2dT(i,j)>0.5) & + P_Stokes_l = P_Stokes_l + h_l*(dUs_dz_l*u_l+dVs_dz_l*v_l) + if (G%mask2dT(i,j+1)>0.5) & + P_Stokes_r = P_Stokes_r + h_r*(dUs_dz_r*u_r+dVs_dz_r*v_r) + PFv_Stokes(i,J,k) = (P_Stokes_r - P_Stokes_l)*G%IdyCv(i,J) + z_top_l = z_top_l - h_l + z_top_r = z_top_r - h_r + enddo + endif + enddo ; enddo + + if (CS%id_PFv_Stokes>0) & + call post_data(CS%id_PFv_Stokes, PFv_Stokes, CS%diag) + if (CS%id_PFu_Stokes>0) & + call post_data(CS%id_PFu_Stokes, PFu_Stokes, CS%diag) + +end subroutine Stokes_PGF + !> Solver to add Coriolis-Stokes to model !! Still in development and not meant for general use. !! Can be activated (with code intervention) for LES comparison From 9fbdf11ec3c5bf0b698aed7b6eca674b388d3950 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Fri, 4 Jun 2021 13:47:21 -0400 Subject: [PATCH 2/4] Finish incomplete merge in MOM_wave_interface. --- src/user/MOM_wave_interface.F90 | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 611387c41e..71f29eb852 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -460,19 +460,15 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) CS%diag%axesCvL,Time,'3d Stokes drift (y)', 'm s-1', conversion=US%L_T_to_m_s) CS%id_3dstokes_x = register_diag_field('ocean_model','3d_stokes_x', & CS%diag%axesCuL,Time,'3d Stokes drift (x)', 'm s-1', conversion=US%L_T_to_m_s) -<<<<<<< HEAD CS%id_ddt_3dstokes_y = register_diag_field('ocean_model','dvdt_Stokes', & - CS%diag%axesCvL,Time,'d/dt Stokes drift (meridional)','m s-2') + CS%diag%axesCvL,Time,'d/dt Stokes drift (meridional)','m s-2')!Needs conversion CS%id_ddt_3dstokes_x = register_diag_field('ocean_model','dudt_Stokes', & - CS%diag%axesCuL,Time,'d/dt Stokes drift (zonal)','m s-2') + CS%diag%axesCuL,Time,'d/dt Stokes drift (zonal)','m s-2')!Needs conversion CS%id_PFv_Stokes = register_diag_field('ocean_model','PFv_Stokes', & CS%diag%axesCvL,Time,'PF from Stokes drift (meridional)','m s-2')!Needs conversion CS%id_PFu_Stokes = register_diag_field('ocean_model','PFu_Stokes', & CS%diag%axesCuL,Time,'PF from Stokes drift (zonal)','m s-2')!Needs conversion - CS%id_La_turb = register_diag_field('ocean_model','La_turbulent',& -======= CS%id_La_turb = register_diag_field('ocean_model','La_turbulent', & ->>>>>>> dev/gfdl CS%diag%axesT1,Time,'Surface (turbulent) Langmuir number','nondim') end subroutine MOM_wave_interface_init From a6de11c92a09abfc748e4970aa3aa09e7f0c98f5 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Mon, 7 Jun 2021 13:07:33 -0400 Subject: [PATCH 3/4] Remote application of Stokes tendency on thermodynamic step - This only makes sense to add on the dynamics timestep. --- src/core/MOM.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index dddfea228d..f7167a1a52 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -656,10 +656,10 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS else ! not do_dyn. if (CS%UseWaves) then ! Diagnostics are not enabled in this call. call Update_Stokes_Drift(G, GV, US, Waves, h, fluxes%ustar, dt) - if (Waves%Stokes_DDT) then - u(:,:,:) = u(:,:,:) + Waves%ddt_us_x(:,:,:)*dt - v(:,:,:) = v(:,:,:) + Waves%ddt_us_y(:,:,:)*dt - endif + !if (Waves%Stokes_DDT) then + ! u(:,:,:) = u(:,:,:) + Waves%ddt_us_x(:,:,:)*dt + ! v(:,:,:) = v(:,:,:) + Waves%ddt_us_y(:,:,:)*dt + !endif endif endif From c5f1d442e162e28e0f9c7918a2f90d6557aa6637 Mon Sep 17 00:00:00 2001 From: "brandon.reichl" Date: Mon, 7 Jun 2021 13:08:43 -0400 Subject: [PATCH 4/4] Updating name of Stokes PGF routine to be more descriptive. --- src/core/MOM_dynamics_split_RK2.F90 | 6 +- src/user/MOM_wave_interface.F90 | 255 ++++++++++++++-------------- 2 files changed, 135 insertions(+), 126 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 34539dbea7..bffd07f66c 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -60,7 +60,7 @@ module MOM_dynamics_split_RK2 use MOM_vert_friction, only : updateCFLtruncationValue 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_wave_interface, only: wave_parameters_CS, Stokes_PGF_Add_FD implicit none ; private @@ -463,7 +463,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s Use_Stokes_PGF = associated(Waves) if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF if (Use_Stokes_PGF) then - call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + call Stokes_PGF_Add_FD(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) ! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv ! will therefore report the sum total PGF and we avoid other ! modifications in the code. The PFu_Stokes can be output within the waves routines. @@ -725,7 +725,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s Use_Stokes_PGF = associated(Waves) if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF if (Use_Stokes_PGF) then - call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + call Stokes_PGF_Add_FD(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) if (.not.Waves%Passive_Stokes_PGF) then do k=1,nz do j=js,je ; do I=Isq,Ieq diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 71f29eb852..ca93ae0a93 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -30,7 +30,8 @@ module MOM_wave_interface ! called in step_mom. public get_Langmuir_Number ! Public interface to compute Langmuir number called from ! ePBL or KPP routines. -public Stokes_PGF ! Public interface to compute Stokes-shear modifications to pressure gradient force +public Stokes_PGF_Add_FD ! Public interface to compute Stokes-shear modifications to pressure gradient force + ! using an additive, finite difference method public StokesMixing ! NOT READY - Public interface to add down-Stokes gradient ! momentum mixing (e.g. the approach of Harcourt 2013/2015) public CoriolisStokes ! NOT READY - Public interface to add Coriolis-Stokes acceleration @@ -54,7 +55,7 @@ module MOM_wave_interface !! True if Stokes vortex force is used logical, public :: Stokes_PGF !< Developmental: !! True if Stokes shear pressure Gradient force is used - logical, public :: Passive_Stokes_PGF !< Keeps Stokes_PGF on, but doesn't affect dynamics + logical, public :: Passive_Stokes_PGF !< Keeps Stokes_PGF on, but doesn't affect dynamics logical, public :: Stokes_DDT !< Developmental: !! True if Stokes d/dt is used real, allocatable, dimension(:,:,:), public :: & @@ -461,9 +462,9 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) CS%id_3dstokes_x = register_diag_field('ocean_model','3d_stokes_x', & CS%diag%axesCuL,Time,'3d Stokes drift (x)', 'm s-1', conversion=US%L_T_to_m_s) CS%id_ddt_3dstokes_y = register_diag_field('ocean_model','dvdt_Stokes', & - CS%diag%axesCvL,Time,'d/dt Stokes drift (meridional)','m s-2')!Needs conversion + CS%diag%axesCvL,Time,'d/dt Stokes drift (meridional)','m s-2') CS%id_ddt_3dstokes_x = register_diag_field('ocean_model','dudt_Stokes', & - CS%diag%axesCuL,Time,'d/dt Stokes drift (zonal)','m s-2')!Needs conversion + CS%diag%axesCuL,Time,'d/dt Stokes drift (zonal)','m s-2') CS%id_PFv_Stokes = register_diag_field('ocean_model','PFv_Stokes', & CS%diag%axesCvL,Time,'PF from Stokes drift (meridional)','m s-2')!Needs conversion CS%id_PFu_Stokes = register_diag_field('ocean_model','PFu_Stokes', & @@ -565,14 +566,15 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt) real :: La ! The local Langmuir number [nondim] integer :: ii, jj, kk, b, iim1, jjm1 real :: idt ! 1 divided by the time step - + one_cm = 0.01*US%m_to_Z min_level_thick_avg = 1.e-3*US%m_to_Z idt = 1.0/dt + ! Getting Stokes drift profile from previous step CS%ddt_us_x(:,:,:) = CS%US_x(:,:,:) CS%ddt_us_y(:,:,:) = CS%US_y(:,:,:) - + ! 1. If Test Profile Option is chosen ! Computing mid-point value from surface value and decay wavelength if (CS%WaveMethod==TESTPROF) then @@ -810,9 +812,11 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt) enddo enddo + ! Finding tendency of Stokes drift over the time step to apply + ! as an acceleration to the models current. CS%ddt_us_x(:,:,:) = (CS%US_x(:,:,:) - CS%ddt_us_x(:,:,:)) * idt CS%ddt_us_y(:,:,:) = (CS%US_y(:,:,:) - CS%ddt_us_y(:,:,:)) * idt - + ! Output any desired quantities if (CS%id_surfacestokes_y>0) & call post_data(CS%id_surfacestokes_y, CS%us0_y, CS%diag) @@ -1422,8 +1426,56 @@ subroutine StokesMixing(G, GV, dt, h, u, v, Waves ) end subroutine StokesMixing -!> Computes tendency due to Stokes pressure gradient force -subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) +!> Solver to add Coriolis-Stokes to model +!! Still in development and not meant for general use. +!! Can be activated (with code intervention) for LES comparison +!! CHECK THAT RIGHT TIMESTEP IS PASSED IF YOU USE THIS** +!! +!! Not accessed in the standard code. +subroutine CoriolisStokes(G, GV, dt, h, u, v, WAVES, US) + type(ocean_grid_type), & + intent(in) :: G !< Ocean grid + type(verticalGrid_type), & + intent(in) :: GV !< Ocean vertical grid + real, intent(in) :: dt !< Time step of MOM6 [T ~> s] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: u !< Velocity i-component [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: v !< Velocity j-component [L T-1 ~> m s-1] + type(Wave_parameters_CS), & + pointer :: Waves !< Surface wave related control structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + ! Local variables + real :: DVel ! A rescaled velocity change [L T-2 ~> m s-2] + integer :: i,j,k + + do k = 1, GV%ke +do j = G%jsc, G%jec + do I = G%iscB, G%iecB + DVel = 0.25*(WAVES%us_y(i,j+1,k)+WAVES%us_y(i-1,j+1,k))*G%CoriolisBu(i,j+1) + & + 0.25*(WAVES%us_y(i,j,k)+WAVES%us_y(i-1,j,k))*G%CoriolisBu(i,j) + u(I,j,k) = u(I,j,k) + DVEL*dt + enddo + enddo + enddo + + do k = 1, GV%ke + do J = G%jscB, G%jecB + do i = G%isc, G%iec + DVel = 0.25*(WAVES%us_x(i+1,j,k)+WAVES%us_x(i+1,j-1,k))*G%CoriolisBu(i+1,j) + & + 0.25*(WAVES%us_x(i,j,k)+WAVES%us_x(i,j-1,k))*G%CoriolisBu(i,j) + v(i,J,k) = v(i,j,k) - DVEL*dt + enddo + enddo + enddo +end subroutine CoriolisStokes + + +!> Computes tendency due to Stokes pressure gradient force using an +!! additive finite difference method +subroutine Stokes_PGF_Add_FD(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) type(ocean_grid_type), & intent(in) :: G !< Ocean grid type(verticalGrid_type), & @@ -1449,17 +1501,16 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) real :: z_top_l, z_top_r ! The height of the top of the cell (left/right index) [Z ~> m]. real :: z_mid_l, z_mid_r ! The height of the middle of the cell (left/right index) [Z ~> m]. real :: h_l, h_r ! The thickness of the cell (left/right index) [Z ~> m]. - real :: wavenum,TwoKexpL, TwoKexpR !TMP DELETE THIS - integer :: i,j,k - - ! Comput the Stokes contribution to the pressure gradient force + real :: TwoKexpL, TwoKexpR + integer :: i,j,k,l + ! Compute the Stokes contribution to the pressure gradient force PFu_Stokes(:,:,:) = 0.0 PFv_Stokes(:,:,:) = 0.0 - wavenum = (2.*3.14)/50. do j = G%jsc, G%jec ; do I = G%iscB, G%iecB if (G%mask2dCu(I,j)>0.5) then + z_top_l = 0.0 z_top_r = 0.0 P_Stokes_l = 0.0 @@ -1469,37 +1520,39 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) h_r = h(i+1,j,k) z_mid_l = z_top_l - 0.5*h_l z_mid_r = z_top_r - 0.5*h_r - TwoKexpL = (2.*wavenum)*exp(2*wavenum*z_mid_l) - TwoKexpR = (2.*wavenum)*exp(2*wavenum*z_mid_r) - !UL -> I-1 & I, j - !UR -> I & I+1, j - !VL -> i, J-1 & J - !VR -> i+1, J-1 & J - dUs_dz_l = TwoKexpL*0.5 * & - (CS%Us0_x(I-1,j)*G%mask2dCu(I-1,j) + & - CS%Us0_x(I,j)*G%mask2dCu(I,j)) - dUs_dz_r = TwoKexpR*0.5 * & - (CS%Us0_x(I,j)*G%mask2dCu(I,j) + & - CS%Us0_x(I+1,j)*G%mask2dCu(I+1,j)) - dVs_dz_l = TwoKexpL*0.5 * & - (CS%Us0_y(i,J-1)*G%mask2dCv(i,J-1) + & - CS%Us0_y(i,J)*G%mask2dCv(i,J)) - dVs_dz_r = TwoKexpR*0.5 * & - (CS%Us0_y(i+1,J-1)*G%mask2dCv(i+1,J-1) + & - CS%Us0_y(i+1,J)*G%mask2dCv(i+1,J)) - u_l = 0.5*(u(I-1,j,k)*G%mask2dCu(I-1,j) + & - u(I,j,k)*G%mask2dCu(I,j)) - u_r = 0.5*(u(I,j,k)*G%mask2dCu(I,j) + & - u(I+1,j,k)*G%mask2dCu(I+1,j)) - v_l = 0.5*(v(i,J-1,k)*G%mask2dCv(i,J-1) + & - v(i,J,k)*G%mask2dCv(i,J)) - v_r = 0.5*(v(i+1,J-1,k)*G%mask2dCv(i+1,J-1) + & - v(i+1,J,k)*G%mask2dCv(i+1,J)) - if (G%mask2dT(i,j)>0.5) & - P_Stokes_l = P_Stokes_l + h_l*(dUs_dz_l*u_l+dVs_dz_l*v_l) - if (G%mask2dT(i+1,j)>0.5) & - P_Stokes_r = P_Stokes_r + h_r*(dUs_dz_r*u_r+dVs_dz_r*v_r) - PFu_Stokes(I,j,k) = (P_Stokes_r - P_Stokes_l)*G%IdxCu(I,j) + do l = 1, CS%numbands + TwoKexpL = (2.*CS%WaveNum_Cen(l))*exp(2*CS%WaveNum_Cen(l)*z_mid_l) + TwoKexpR = (2.*CS%WaveNum_Cen(l))*exp(2*CS%WaveNum_Cen(l)*z_mid_r) + !UL -> I-1 & I, j + !UR -> I & I+1, j + !VL -> i, J-1 & J + !VR -> i+1, J-1 & J + dUs_dz_l = TwoKexpL*0.5 * & + (CS%Stkx0(I-1,j,l)*G%mask2dCu(I-1,j) + & + CS%Stkx0(I,j,l)*G%mask2dCu(I,j)) + dUs_dz_r = TwoKexpR*0.5 * & + (CS%Stkx0(I,j,l)*G%mask2dCu(I,j) + & + CS%Stkx0(I+1,j,l)*G%mask2dCu(I+1,j)) + dVs_dz_l = TwoKexpL*0.5 * & + (CS%Stky0(i,J-1,l)*G%mask2dCv(i,J-1) + & + CS%Stky0(i,J,l)*G%mask2dCv(i,J)) + dVs_dz_r = TwoKexpR*0.5 * & + (CS%Stky0(i+1,J-1,l)*G%mask2dCv(i+1,J-1) + & + CS%Stky0(i+1,J,l)*G%mask2dCv(i+1,J)) + u_l = 0.5*(u(I-1,j,k)*G%mask2dCu(I-1,j) + & + u(I,j,k)*G%mask2dCu(I,j)) + u_r = 0.5*(u(I,j,k)*G%mask2dCu(I,j) + & + u(I+1,j,k)*G%mask2dCu(I+1,j)) + v_l = 0.5*(v(i,J-1,k)*G%mask2dCv(i,J-1) + & + v(i,J,k)*G%mask2dCv(i,J)) + v_r = 0.5*(v(i+1,J-1,k)*G%mask2dCv(i+1,J-1) + & + v(i+1,J,k)*G%mask2dCv(i+1,J)) + if (G%mask2dT(i,j)>0.5) & + P_Stokes_l = P_Stokes_l + h_l*(dUs_dz_l*u_l+dVs_dz_l*v_l) + if (G%mask2dT(i+1,j)>0.5) & + P_Stokes_r = P_Stokes_r + h_r*(dUs_dz_r*u_r+dVs_dz_r*v_r) + PFu_Stokes(I,j,k) = PFu_Stokes(I,j,k) + (P_Stokes_r - P_Stokes_l)*G%IdxCu(I,j) + enddo z_top_l = z_top_l - h_l z_top_r = z_top_r - h_r enddo @@ -1516,37 +1569,39 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) h_r = h(i,j+1,k) z_mid_l = z_top_l - 0.5*h_l z_mid_r = z_top_r - 0.5*h_r - TwoKexpL = (2.*wavenum)*exp(2*wavenum*z_mid_l) - TwoKexpR = (2.*wavenum)*exp(2*wavenum*z_mid_r) - !UL -> I-1 & I, j - !UR -> I-1 & I, j+1 - !VL -> i, J & J-1 - !VR -> i, J & J+1 - dUs_dz_l = TwoKexpL*0.5 * & - (CS%Us0_x(I-1,j)*G%mask2dCu(I-1,j) + & - CS%Us0_x(I,j)*G%mask2dCu(I,j)) - dUs_dz_r = TwoKexpR*0.5 * & - (CS%Us0_x(I-1,j+1)*G%mask2dCu(I-1,j+1) + & - CS%Us0_x(I,j+1)*G%mask2dCu(I,j+1)) - dVs_dz_l = TwoKexpL*0.5 * & - (CS%Us0_y(i,J-1)*G%mask2dCv(i,J-1) + & - CS%Us0_y(i,J)*G%mask2dCv(i,J)) - dVs_dz_r = TwoKexpR*0.5 * & - (CS%Us0_y(i,J)*G%mask2dCv(i,J) + & - CS%Us0_y(i,J+1)*G%mask2dCv(i,J+1)) - u_l = 0.5*(u(I-1,j,k)*G%mask2dCu(I-1,j) + & - u(I,j,k)*G%mask2dCu(I,j)) - u_r = 0.5*(u(I-1,j+1,k)*G%mask2dCu(I-1,j+1) + & - u(I,j+1,k)*G%mask2dCu(I,j+1)) - v_l = 0.5*(v(i,J-1,k)*G%mask2dCv(i,J-1) + & - v(i,J,k)*G%mask2dCv(i,J)) - v_r = 0.5*(v(i,J,k)*G%mask2dCv(i,J) + & - v(i,J+1,k)*G%mask2dCv(i,J+1)) - if (G%mask2dT(i,j)>0.5) & - P_Stokes_l = P_Stokes_l + h_l*(dUs_dz_l*u_l+dVs_dz_l*v_l) - if (G%mask2dT(i,j+1)>0.5) & - P_Stokes_r = P_Stokes_r + h_r*(dUs_dz_r*u_r+dVs_dz_r*v_r) - PFv_Stokes(i,J,k) = (P_Stokes_r - P_Stokes_l)*G%IdyCv(i,J) + do l = 1, CS%numbands + TwoKexpL = (2.*CS%WaveNum_Cen(l))*exp(2*CS%WaveNum_Cen(l)*z_mid_l) + TwoKexpR = (2.*CS%WaveNum_Cen(l))*exp(2*CS%WaveNum_Cen(l)*z_mid_r) + !UL -> I-1 & I, j + !UR -> I-1 & I, j+1 + !VL -> i, J & J-1 + !VR -> i, J & J+1 + dUs_dz_l = TwoKexpL*0.5 * & + (CS%Stkx0(I-1,j,l)*G%mask2dCu(I-1,j) + & + CS%Stkx0(I,j,l)*G%mask2dCu(I,j)) + dUs_dz_r = TwoKexpR*0.5 * & + (CS%Stkx0(I-1,j+1,l)*G%mask2dCu(I-1,j+1) + & + CS%Stkx0(I,j+1,l)*G%mask2dCu(I,j+1)) + dVs_dz_l = TwoKexpL*0.5 * & + (CS%Stky0(i,J-1,l)*G%mask2dCv(i,J-1) + & + CS%Stky0(i,J,l)*G%mask2dCv(i,J)) + dVs_dz_r = TwoKexpR*0.5 * & + (CS%Stky0(i,J,l)*G%mask2dCv(i,J) + & + CS%Stky0(i,J+1,l)*G%mask2dCv(i,J+1)) + u_l = 0.5*(u(I-1,j,k)*G%mask2dCu(I-1,j) + & + u(I,j,k)*G%mask2dCu(I,j)) + u_r = 0.5*(u(I-1,j+1,k)*G%mask2dCu(I-1,j+1) + & + u(I,j+1,k)*G%mask2dCu(I,j+1)) + v_l = 0.5*(v(i,J-1,k)*G%mask2dCv(i,J-1) + & + v(i,J,k)*G%mask2dCv(i,J)) + v_r = 0.5*(v(i,J,k)*G%mask2dCv(i,J) + & + v(i,J+1,k)*G%mask2dCv(i,J+1)) + if (G%mask2dT(i,j)>0.5) & + P_Stokes_l = P_Stokes_l + h_l*(dUs_dz_l*u_l+dVs_dz_l*v_l) + if (G%mask2dT(i,j+1)>0.5) & + P_Stokes_r = P_Stokes_r + h_r*(dUs_dz_r*u_r+dVs_dz_r*v_r) + PFv_Stokes(i,J,k) = PFv_Stokes(i,J,k) + (P_Stokes_r - P_Stokes_l)*G%IdyCv(i,J) + enddo z_top_l = z_top_l - h_l z_top_r = z_top_r - h_r enddo @@ -1558,53 +1613,7 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) if (CS%id_PFu_Stokes>0) & call post_data(CS%id_PFu_Stokes, PFu_Stokes, CS%diag) -end subroutine Stokes_PGF - -!> Solver to add Coriolis-Stokes to model -!! Still in development and not meant for general use. -!! Can be activated (with code intervention) for LES comparison -!! CHECK THAT RIGHT TIMESTEP IS PASSED IF YOU USE THIS** -!! -!! Not accessed in the standard code. -subroutine CoriolisStokes(G, GV, dt, h, u, v, WAVES, US) - type(ocean_grid_type), & - intent(in) :: G !< Ocean grid - type(verticalGrid_type), & - intent(in) :: GV !< Ocean vertical grid - real, intent(in) :: dt !< Time step of MOM6 [T ~> s] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: u !< Velocity i-component [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(inout) :: v !< Velocity j-component [L T-1 ~> m s-1] - type(Wave_parameters_CS), & - pointer :: Waves !< Surface wave related control structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - ! Local variables - real :: DVel ! A rescaled velocity change [L T-2 ~> m s-2] - integer :: i,j,k - - do k = 1, GV%ke - do j = G%jsc, G%jec - do I = G%iscB, G%iecB - DVel = 0.25*(WAVES%us_y(i,j+1,k)+WAVES%us_y(i-1,j+1,k))*G%CoriolisBu(i,j+1) + & - 0.25*(WAVES%us_y(i,j,k)+WAVES%us_y(i-1,j,k))*G%CoriolisBu(i,j) - u(I,j,k) = u(I,j,k) + DVEL*dt - enddo - enddo - enddo - - do k = 1, GV%ke - do J = G%jscB, G%jecB - do i = G%isc, G%iec - DVel = 0.25*(WAVES%us_x(i+1,j,k)+WAVES%us_x(i+1,j-1,k))*G%CoriolisBu(i+1,j) + & - 0.25*(WAVES%us_x(i,j,k)+WAVES%us_x(i,j-1,k))*G%CoriolisBu(i,j) - v(i,J,k) = v(i,j,k) - DVEL*dt - enddo - enddo - enddo -end subroutine CoriolisStokes +end subroutine Stokes_PGF_Add_FD !> Computes wind speed from ustar_air based on COARE 3.5 Cd relationship !! Probably doesn't belong in this module, but it is used here to estimate