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