Index: palm/trunk/SOURCE/ocean_mod.f90
===================================================================
--- palm/trunk/SOURCE/ocean_mod.f90	(revision 3294)
+++ palm/trunk/SOURCE/ocean_mod.f90	(revision 3302)
@@ -25,8 +25,9 @@
 ! -----------------
 ! $Id$
+! Craik Leibovich force (Stokes drift) + wave breaking effect added
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! initial revision
-! 
-! 
-! 
+!
 !
 ! Authors:
@@ -36,8 +37,7 @@
 ! Description:
 ! ------------
-!> This module contains relevant code for PALM's ocean option, e.g. the
+!> This module contains relevant code for PALM's ocean mode, e.g. the
 !> prognostic equation for salinity, the equation of state for seawater,
-!> and the Craik Leibovich force (Stokes force)
-!>
+!> the Craik Leibovich force (Stokes force), and wave breaking effects
 !------------------------------------------------------------------------------!
  MODULE ocean_mod
@@ -66,16 +66,20 @@
 
     INTEGER(iwp) ::  ibc_sa_t   !< integer flag for bc_sa_t
+    INTEGER(iwp) ::  iran_ocean = -1234567  !< random number used for wave breaking effects
 
     INTEGER(iwp) ::  sa_vertical_gradient_level_ind(10) = -9999  !< grid index values of sa_vertical_gradient_level(s)
 
-    LOGICAL ::  stokes_force = .FALSE.              !< switch to switch on the Stokes force
-
-    REAL(wp) ::  langmuir_number = 0.34_wp    !< turbulent Langmuir number, default from Li et al., 2005
+    LOGICAL ::  stokes_force = .FALSE.        !< switch to switch on the Stokes force
+    LOGICAL ::  wave_breaking = .FALSE.       !< switch to switch on wave breaking effects
+
+    REAL(wp) ::  alpha_wave_breaking = 3.0_wp !< coefficient for wave breaking generated turbulence from Noh et al. (2004), JPO
     REAL(wp) ::  prho_reference               !< reference state of potential density at ocean surface
     REAL(wp) ::  sa_surface = 35.0_wp         !< surface salinity, namelist parameter
     REAL(wp) ::  sa_vertical_gradient(10) = 0.0_wp               !< namelist parameter
     REAL(wp) ::  sa_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
-    REAL(wp) ::  stokes_waveheight =  1.0_wp  !< typical wave height assumed for Stokes velocity
-    REAL(wp) ::  stokes_wavelength = 40.0_wp  !< typical wavelength assumed for Stokes velocity
+    REAL(wp) ::  stokes_waveheight = 0.0_wp  !< wave height assumed for Stokes drift velocity
+    REAL(wp) ::  stokes_wavelength = 0.0_wp  !< wavelength assumed for Stokes drift velocity
+    REAL(wp) ::  timescale_wave_breaking     !< time scale of random forcing
+    REAL(wp) ::  u_star_wave_breaking        !< to store the absolute value of friction velocity at the ocean surface
 
     REAL(wp), DIMENSION(12), PARAMETER ::  nom =                               &
@@ -100,8 +104,7 @@
     SAVE
 
-    PUBLIC ::  bc_sa_t, ibc_sa_t, langmuir_number, prho_reference, sa_surface, &
+    PUBLIC ::  bc_sa_t, ibc_sa_t, prho_reference, sa_surface,                  &
                sa_vertical_gradient, sa_vertical_gradient_level,               &
-               sa_vertical_gradient_level_ind, stokes_force,                   &
-               stokes_waveheight, stokes_wavelength
+               sa_vertical_gradient_level_ind, stokes_force, wave_breaking
 
 
@@ -139,4 +142,8 @@
     END INTERFACE ocean_data_output_3d
 
+    INTERFACE ocean_header
+       MODULE PROCEDURE ocean_header
+    END INTERFACE ocean_header
+
     INTERFACE ocean_init
        MODULE PROCEDURE ocean_init
@@ -180,5 +187,13 @@
     END INTERFACE ocean_3d_data_averaging
 
-    SAVE
+    INTERFACE stokes_drift_terms
+       MODULE PROCEDURE stokes_drift_terms
+       MODULE PROCEDURE stokes_drift_terms_ij
+    END INTERFACE stokes_drift_terms
+
+    INTERFACE wave_breaking_term
+       MODULE PROCEDURE wave_breaking_term
+       MODULE PROCEDURE wave_breaking_term_ij
+    END INTERFACE wave_breaking_term
 
     PRIVATE
@@ -188,8 +203,9 @@
            ocean_check_data_output_pr, ocean_check_parameters,                 &
            ocean_data_output_2d, ocean_data_output_3d,                         &
-           ocean_define_netcdf_grid, ocean_init, ocean_init_arrays,            &
-           ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel,      &
-           ocean_rrd_global, ocean_rrd_local, ocean_wrd_global,                &
-           ocean_wrd_local, ocean_3d_data_averaging
+           ocean_define_netcdf_grid, ocean_header, ocean_init,                 &
+           ocean_init_arrays, ocean_parin, ocean_prognostic_equations,         &
+           ocean_swap_timelevel, ocean_rrd_global, ocean_rrd_local,            &
+           ocean_wrd_global, ocean_wrd_local, ocean_3d_data_averaging,         &
+           stokes_drift_terms, wave_breaking_term
 
 
@@ -484,8 +500,8 @@
 
 
-    NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux,                 &
-             langmuir_number, sa_surface, sa_vertical_gradient,                &
-             sa_vertical_gradient_level, stokes_force, stokes_waveheight,      &
-             stokes_wavelength, top_salinityflux, wall_salinityflux
+    NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux, sa_surface,     &
+             sa_vertical_gradient, sa_vertical_gradient_level,                 &
+             stokes_waveheight, stokes_wavelength, top_salinityflux,           &
+             wall_salinityflux, wave_breaking
 
 !
@@ -580,4 +596,21 @@
                         'top_salinityflux /= 0.0'
        CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
+!-- Check if Stokes force is to be used
+    IF ( stokes_waveheight /= 0.0_wp  .AND.  stokes_wavelength /= 0.0_wp )  THEN
+       stokes_force = .TRUE.
+    ELSE
+       IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength > 0.0_wp ) &
+            .OR.                                                               &
+            ( stokes_waveheight > 0.0_wp .AND. stokes_wavelength <= 0.0_wp ) &
+            .OR.                                                               &
+            ( stokes_waveheight < 0.0_wp .AND. stokes_wavelength < 0.0_wp  ) ) &
+       THEN
+          message_string = 'wrong settings for stokes_wavelength and/or ' //   &
+                           'stokes_waveheight'
+          CALL message( 'ocean_check_parameters', 'PA0460', 1, 2, 0, 6, 0 )
+       ENDIF
     ENDIF
 
@@ -1045,4 +1078,36 @@
 ! Description:
 ! ------------
+!> Header output for ocean parameters
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_header( io )
+
+
+    IMPLICIT NONE
+
+    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
+
+!
+!-- Write ocean header
+    WRITE( io, 1 )
+    IF ( stokes_force  )  WRITE( io, 2 )  stokes_waveheight, stokes_wavelength
+    IF ( wave_breaking )  THEN
+       WRITE( io, 3 )  alpha_wave_breaking, timescale_wave_breaking
+    ENDIF
+
+1   FORMAT (//' Ocean settings:'/                                              &
+              ' ------------------------------------------'/)
+2   FORMAT ('    --> Craik-Leibovich vortex force and Stokes drift switched',  &
+                     ' on'/                                                    &
+            '        waveheight: ',F4.1,' m   wavelength: ',F6.1,' m')
+3   FORMAT ('    --> wave breaking generated turbulence switched on'/          &
+            '        alpha:    ',F4.1/                                         &
+            '        timescale:',F5.1,' s')
+
+ END SUBROUTINE ocean_header
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Allocate arrays and assign pointers.
 !------------------------------------------------------------------------------!
@@ -1093,12 +1158,17 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzu, hyp, pt_init, ref_state, zu, zw
+        ONLY:  dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw,    &
+               v_stokes_zu, v_stokes_zw, zu, zw
 
     USE basic_constants_and_equations_mod,                                     &
         ONLY:  g
 
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  pi
+
     USE control_parameters,                                                    &
         ONLY:  initializing_actions, molecular_viscosity, rho_surface,         &
-               rho_reference, surface_pressure, use_single_reference_value
+               rho_reference, surface_pressure, top_momentumflux_u,            &
+               top_momentumflux_v, use_single_reference_value
 
     USE indices,                                                               &
@@ -1107,5 +1177,6 @@
     USE kinds
 
-    USE pegrid
+    USE pegrid,                                                                &
+        ONLY:  myid
 
     USE statistics,                                                            &
@@ -1119,7 +1190,11 @@
     INTEGER(iwp) ::  n  !< loop index
 
-    REAL(wp)     ::  dum   !< dummy argument
-    REAL(wp)     ::  pt_l  !< local scalar for pt used in equation of state function
-    REAL(wp)     ::  sa_l  !< local scalar for sa used in equation of state function
+    REAL(wp) ::  alpha !< angle of surface stress
+    REAL(wp) ::  dum   !< dummy argument
+    REAL(wp) ::  pt_l  !< local scalar for pt used in equation of state function
+    REAL(wp) ::  sa_l  !< local scalar for sa used in equation of state function
+    REAL(wp) ::  velocity_amplitude  !< local scalar for amplitude of Stokes drift velocity
+    REAL(wp) ::  x     !< temporary variable to store surface stress along x
+    REAL(wp) ::  y     !< temporary variable to store surface stress along y
 
     REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init  !< local array for initial density
@@ -1270,4 +1345,68 @@
     ENDIF
 
+!
+!-- Calculate the Stokes drift velocity profile
+    IF ( stokes_force )  THEN
+
+!
+!--    First, calculate angle of surface stress
+       x = -top_momentumflux_u
+       y = -top_momentumflux_v
+       IF ( x == 0.0_wp )  THEN
+          IF ( y > 0.0_wp )  THEN
+             alpha = pi / 2.0_wp
+          ELSEIF ( y < 0.0_wp )  THEN
+             alpha = 3.0_wp * pi / 2.0_wp
+          ENDIF
+       ELSE
+          IF ( x < 0.0_wp )  THEN
+             alpha = ATAN( y / x ) + pi
+          ELSE
+             IF ( y < 0.0_wp )  THEN
+                alpha = ATAN( y / x ) + 2.0_wp * pi
+             ELSE
+                alpha = ATAN( y / x )
+             ENDIF
+          ENDIF
+       ENDIF
+
+       velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 *&
+                            SQRT( g * stokes_wavelength / ( 2.0_wp * pi ) )
+
+       DO  k = nzb, nzt
+          u_stokes_zu(k) = velocity_amplitude * COS( alpha ) *                 &
+                           EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
+          u_stokes_zw(k) = velocity_amplitude * COS( alpha ) *                 &
+                           EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
+          v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) *                 &
+                           EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
+          v_stokes_zw(k) = velocity_amplitude * SIN( alpha ) *                 &
+                           EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
+       ENDDO
+       u_stokes_zu(nzt+1) = u_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
+       u_stokes_zw(nzt+1) = u_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
+       v_stokes_zu(nzt+1) = v_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
+       v_stokes_zw(nzt+1) = v_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
+
+    ENDIF
+
+!
+!-- Wave breaking effects
+    IF ( wave_breaking )  THEN
+!
+!--    Calculate friction velocity at ocean surface
+       u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 +              &
+                                          top_momentumflux_v**2 ) )
+!
+!--    Set the time scale of random forcing. The vertical grid spacing at the
+!--    ocean surface is assumed as the length scale of turbulence.
+!--    Formula follows Noh et al. (2004), JPO
+       timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking /     &
+                                 u_star_wave_breaking
+!
+!--    Set random number seeds differently on the processor cores in order to
+!--    create different random number sequences
+       iran_ocean = iran_ocean + myid
+    ENDIF
 
  END SUBROUTINE ocean_init
@@ -1599,7 +1738,4 @@
           READ ( 13 )  bottom_salinityflux
 
-       CASE ( 'langmuir_number' )
-          READ ( 13 )  langmuir_number
-
        CASE ( 'sa_init' )
           READ ( 13 )  sa_init
@@ -1614,7 +1750,4 @@
           READ ( 13 )  sa_vertical_gradient_level
 
-       CASE ( 'stokes_force' )
-          READ ( 13 )  stokes_force
-
        CASE ( 'stokes_waveheight' )
           READ ( 13 )  stokes_waveheight
@@ -1628,4 +1761,7 @@
        CASE ( 'wall_salinityflux' )
           READ ( 13 )  wall_salinityflux
+
+       CASE ( 'wave_breaking' )
+          READ ( 13 )  wave_breaking
 
        CASE DEFAULT
@@ -1731,7 +1867,4 @@
     WRITE ( 14 )  bottom_salinityflux
 
-    CALL wrd_write_string( 'langmuir_number' )
-    WRITE ( 14 )  langmuir_number
-
     CALL wrd_write_string( 'sa_init' )
     WRITE ( 14 )  sa_init
@@ -1746,7 +1879,4 @@
     WRITE ( 14 )  sa_vertical_gradient_level
 
-    CALL wrd_write_string( 'stokes_force' )
-    WRITE ( 14 )  stokes_force
-
     CALL wrd_write_string( 'stokes_waveheight' )
     WRITE ( 14 )  stokes_waveheight
@@ -1760,4 +1890,7 @@
     CALL wrd_write_string( 'wall_salinityflux' )
     WRITE ( 14 )  wall_salinityflux
+
+    CALL wrd_write_string( 'wave_breaking' )
+    WRITE ( 14 )  wave_breaking
 
  END SUBROUTINE ocean_wrd_global
@@ -1792,3 +1925,313 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates the Craik Leibovich vortex force and the additional
+!> effect of the Stokes drift on the Coriolis force
+!> Call for all gridpoints.
+!------------------------------------------------------------------------------!
+ SUBROUTINE stokes_drift_terms( component )
+
+    USE arrays_3d,                                                             &
+        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
+               v_stokes_zw, w, tend
+
+    USE control_parameters,                                                    &
+        ONLY:  f, fs, message_string
+
+    USE grid_variables,                                                        &
+        ONLY:  ddx, ddy
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nys, nysv, nyn, nzb, nzt
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  component  !< component of momentum equation
+    INTEGER(iwp) ::  i          !< loop index along x
+    INTEGER(iwp) ::  j          !< loop index along y
+    INTEGER(iwp) ::  k          !< loop index along z
+
+
+!
+!-- Compute Stokes terms for the respective velocity components
+    SELECT CASE ( component )
+
+!
+!--    u-component
+       CASE ( 1 )
+          DO  i = nxl, nxr
+             DO  j = nysv, nyn
+                DO  k = nzb+1, nzt
+                   tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (              &
+                                   0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)           &
+                                         + v(k,j,i)   - v(k,j,i-1)   ) * ddx   &
+                                 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy   &
+                                                                )              &
+                                 + f * v_stokes_zu(k)
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--    v-component
+       CASE ( 2 )
+          DO  i = nxl, nxr
+             DO  j = nysv, nyn
+                DO  k = nzb+1, nzt
+                   tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (              &
+                                   0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx   &
+                                 - 0.5 * ( u(k,j,i) - u(k,j-1,i)               &
+                                         + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
+                                                                )              &
+                                 - f * u_stokes_zu(k)
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--    w-component
+       CASE ( 3 )
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
+                                             0.5 * ( u(k+1,j,i) - u(k,j,i)     &
+                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
+                                                   ) * ddzu(k+1)               &
+                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
+                                                   ) * ddx      )              &
+                                             - v_stokes_zw(k) * (              &
+                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
+                                                   ) * ddy                     &
+                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
+                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
+                                                   ) * ddzu(k)  )              &
+                                           + fs * u_stokes_zw(k)
+                ENDDO
+             ENDDO
+          ENDDO
+
+       CASE DEFAULT
+          WRITE( message_string, * ) 'wrong component of Stokes force: ',      &
+                                     component
+          CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
+
+    END SELECT
+
+ END SUBROUTINE stokes_drift_terms
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates the Craik Leibovich vortex force and the additional
+!> effect of the Stokes drift on the Coriolis force
+!> Call for gridpoints i,j.
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE stokes_drift_terms_ij( i, j, component )
+
+    USE arrays_3d,                                                             &
+        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
+               v_stokes_zw, w, tend
+
+    USE control_parameters,                                                    &
+        ONLY:  f, fs, message_string
+
+    USE grid_variables,                                                        &
+        ONLY:  ddx, ddy
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nys, nysv, nyn, nzb, nzt
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  component  !< component of momentum equation
+    INTEGER(iwp) ::  i          !< loop index along x
+    INTEGER(iwp) ::  j          !< loop index along y
+    INTEGER(iwp) ::  k          !< loop incex along z
+
+
+!
+!-- Compute Stokes terms for the respective velocity components
+    SELECT CASE ( component )
+
+!
+!--    u-component
+       CASE ( 1 )
+          DO  k = nzb+1, nzt
+             tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (                    &
+                                     0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)         &
+                                           + v(k,j,i)   - v(k,j,i-1)   ) * ddx &
+                                   - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy &
+                                                          )                    &
+                                       + f * v_stokes_zu(k)
+          ENDDO
+!
+!--    v-component
+       CASE ( 2 )
+          DO  k = nzb+1, nzt
+             tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (                    &
+                                     0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx &
+                                   - 0.5 * ( u(k,j,i) - u(k,j-1,i)             &
+                                           + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
+                                                          )                    &
+                                       - f * u_stokes_zu(k)
+          ENDDO
+
+!
+!--    w-component
+       CASE ( 3 )
+          DO  k = nzb+1, nzt
+             tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
+                                     0.5 * ( u(k+1,j,i) - u(k,j,i)     &
+                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
+                                                   ) * ddzu(k+1)               &
+                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
+                                                   ) * ddx )                   &
+                                       - v_stokes_zw(k) * (                    &
+                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
+                                                   ) * ddy                     &
+                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
+                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
+                                                   ) * ddzu(k)  )              &
+                                       + fs * u_stokes_zw(k)
+          ENDDO
+
+       CASE DEFAULT
+          WRITE( message_string, * ) ' wrong component: ', component
+          CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
+
+    END SELECT
+
+ END SUBROUTINE stokes_drift_terms_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates turbulence generated by wave breaking near the ocean
+!> surface, following a parameterization given in Noh et al. (2004), JPO
+!> Call for all gridpoints.
+!> TODO: so far, this routine only works if the model time step has about the
+!>       same value as the time scale of wave breaking!
+!------------------------------------------------------------------------------!
+ SUBROUTINE wave_breaking_term( component )
+
+    USE arrays_3d,                                                             &
+        ONLY:  u_p, v_p
+
+    USE control_parameters,                                                    &
+        ONLY:  dt_3d, message_string
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzt
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  component  !< component of momentum equation
+    INTEGER(iwp) ::  i          !< loop index along x
+    INTEGER(iwp) ::  j          !< loop index along y
+
+    REAL(wp) ::  random_gauss  !< function that creates a random number with a
+                               !< Gaussian distribution
+
+
+!
+!-- Compute wave breaking terms for the respective velocity components.
+!-- Velocities are directly manipulated, since this is not a real force
+    SELECT CASE ( component )
+
+!
+!--    u-component
+       CASE ( 1 )
+          DO  i = nxlu, nxr
+             DO  j = nys, nyn
+                u_p(nzt,j,i) = u_p(nzt,j,i) +                                  &
+                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
+                               * alpha_wave_breaking * u_star_wave_breaking    &
+                               / timescale_wave_breaking * dt_3d
+             ENDDO
+          ENDDO
+!
+!--    v-component
+       CASE ( 2 )
+          DO  i = nxl, nxr
+             DO  j = nysv, nyn
+                v_p(nzt,j,i) = v_p(nzt,j,i) +                                  &
+                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
+                               * alpha_wave_breaking * u_star_wave_breaking    &
+                               / timescale_wave_breaking * dt_3d
+             ENDDO
+          ENDDO
+
+       CASE DEFAULT
+          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
+                                     component
+          CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
+
+    END SELECT
+
+ END SUBROUTINE wave_breaking_term
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates turbulence generated by wave breaking near the ocean
+!> surface, following a parameterization given in Noh et al. (2004), JPO
+!> Call for gridpoint i,j.
+!> TODO: so far, this routine only works if the model time step has about the
+!>       same value as the time scale of wave breaking!
+!------------------------------------------------------------------------------!
+ SUBROUTINE wave_breaking_term_ij( i, j, component )
+
+    USE arrays_3d,                                                             &
+        ONLY:  u_p, v_p
+
+    USE control_parameters,                                                    &
+        ONLY:  dt_3d, message_string
+
+    USE indices,                                                               &
+        ONLY:  nzt
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  component  !< component of momentum equation
+    INTEGER(iwp) ::  i          !< loop index along x
+    INTEGER(iwp) ::  j          !< loop index along y
+
+    REAL(wp) ::  random_gauss  !< function that creates a random number with a
+                               !< Gaussian distribution
+
+!
+!-- Compute wave breaking terms for the respective velocity components
+    SELECT CASE ( component )
+
+!
+!--    u-/v-component
+       CASE ( 1 )
+          u_p(nzt,j,i) = u_p(nzt,j,i) +                                        &
+                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
+                         * alpha_wave_breaking * u_star_wave_breaking          &
+                         / timescale_wave_breaking * dt_3d
+
+       CASE ( 2 )
+          v_p(nzt,j,i) = v_p(nzt,j,i) +                                        &
+                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
+                         * alpha_wave_breaking * u_star_wave_breaking          &
+                         / timescale_wave_breaking * dt_3d
+
+       CASE DEFAULT
+          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
+                                     component
+          CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
+
+    END SELECT
+
+ END SUBROUTINE wave_breaking_term_ij
+
+
  END MODULE ocean_mod
