Merge remote-tracking branch 'breichl/Stokes_Forces' into pstokes_wav…

…e_cpl_mehtod

src/core/MOM.F90

@@ -652,12 +652,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

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)

src/core/MOM_dynamics_split_RK2.F90

@@ -66,7 +66,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_Add_FD
 
 implicit none ; private
 
@@ -77,11 +77,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
@@ -366,6 +368,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
 
@@ -460,6 +463,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_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.
+      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)")
@@ -476,7 +503,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)")
 
@@ -698,6 +725,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_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
+              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
@@ -732,7 +780,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)")
 
@@ -1323,7 +1371,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