Updating name of Stokes PGF routine to be more descriptive.

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

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