Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Oct 3, 2018 2:39:40 AM (6 years ago)
Author:
raasch
Message:

Craik-Leibovich force and wave breaking effects added to ocean mode, header output for ocean mode

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/ocean_mod.f90

    r3294 r3302  
    2525! -----------------
    2626! $Id$
     27! Craik Leibovich force (Stokes drift) + wave breaking effect added
     28!
     29! 3294 2018-10-01 02:37:10Z raasch
    2730! initial revision
    28 !
    29 !
    30 !
     31!
    3132!
    3233! Authors:
    … …  
    3637! Description:
    3738! ------------
    38 !> This module contains relevant code for PALM's ocean option, e.g. the
     39!> This module contains relevant code for PALM's ocean mode, e.g. the
    3940!> prognostic equation for salinity, the equation of state for seawater,
    40 !> and the Craik Leibovich force (Stokes force)
    41 !>
     41!> the Craik Leibovich force (Stokes force), and wave breaking effects
    4242!------------------------------------------------------------------------------!
    4343 MODULE ocean_mod
    … …  
    6666
    6767    INTEGER(iwp) ::  ibc_sa_t   !< integer flag for bc_sa_t
     68    INTEGER(iwp) ::  iran_ocean = -1234567  !< random number used for wave breaking effects
    6869
    6970    INTEGER(iwp) ::  sa_vertical_gradient_level_ind(10) = -9999  !< grid index values of sa_vertical_gradient_level(s)
    7071
    71     LOGICAL ::  stokes_force = .FALSE.              !< switch to switch on the Stokes force
    72 
    73     REAL(wp) ::  langmuir_number = 0.34_wp    !< turbulent Langmuir number, default from Li et al., 2005
     72    LOGICAL ::  stokes_force = .FALSE.        !< switch to switch on the Stokes force
     73    LOGICAL ::  wave_breaking = .FALSE.       !< switch to switch on wave breaking effects
     74
     75    REAL(wp) ::  alpha_wave_breaking = 3.0_wp !< coefficient for wave breaking generated turbulence from Noh et al. (2004), JPO
    7476    REAL(wp) ::  prho_reference               !< reference state of potential density at ocean surface
    7577    REAL(wp) ::  sa_surface = 35.0_wp         !< surface salinity, namelist parameter
    7678    REAL(wp) ::  sa_vertical_gradient(10) = 0.0_wp               !< namelist parameter
    7779    REAL(wp) ::  sa_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
    78     REAL(wp) ::  stokes_waveheight =  1.0_wp  !< typical wave height assumed for Stokes velocity
    79     REAL(wp) ::  stokes_wavelength = 40.0_wp  !< typical wavelength assumed for Stokes velocity
     80    REAL(wp) ::  stokes_waveheight = 0.0_wp  !< wave height assumed for Stokes drift velocity
     81    REAL(wp) ::  stokes_wavelength = 0.0_wp  !< wavelength assumed for Stokes drift velocity
     82    REAL(wp) ::  timescale_wave_breaking     !< time scale of random forcing
     83    REAL(wp) ::  u_star_wave_breaking        !< to store the absolute value of friction velocity at the ocean surface
    8084
    8185    REAL(wp), DIMENSION(12), PARAMETER ::  nom =                               &
    … …  
    100104    SAVE
    101105
    102     PUBLIC ::  bc_sa_t, ibc_sa_t, langmuir_number, prho_reference, sa_surface, &
     106    PUBLIC ::  bc_sa_t, ibc_sa_t, prho_reference, sa_surface,                  &
    103107               sa_vertical_gradient, sa_vertical_gradient_level,               &
    104                sa_vertical_gradient_level_ind, stokes_force,                   &
    105                stokes_waveheight, stokes_wavelength
     108               sa_vertical_gradient_level_ind, stokes_force, wave_breaking
    106109
    107110
    … …  
    139142    END INTERFACE ocean_data_output_3d
    140143
     144    INTERFACE ocean_header
     145       MODULE PROCEDURE ocean_header
     146    END INTERFACE ocean_header
     147
    141148    INTERFACE ocean_init
    142149       MODULE PROCEDURE ocean_init
    … …  
    180187    END INTERFACE ocean_3d_data_averaging
    181188
    182     SAVE
     189    INTERFACE stokes_drift_terms
     190       MODULE PROCEDURE stokes_drift_terms
     191       MODULE PROCEDURE stokes_drift_terms_ij
     192    END INTERFACE stokes_drift_terms
     193
     194    INTERFACE wave_breaking_term
     195       MODULE PROCEDURE wave_breaking_term
     196       MODULE PROCEDURE wave_breaking_term_ij
     197    END INTERFACE wave_breaking_term
    183198
    184199    PRIVATE
    … …  
    188203           ocean_check_data_output_pr, ocean_check_parameters,                 &
    189204           ocean_data_output_2d, ocean_data_output_3d,                         &
    190            ocean_define_netcdf_grid, ocean_init, ocean_init_arrays,            &
    191            ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel,      &
    192            ocean_rrd_global, ocean_rrd_local, ocean_wrd_global,                &
    193            ocean_wrd_local, ocean_3d_data_averaging
     205           ocean_define_netcdf_grid, ocean_header, ocean_init,                 &
     206           ocean_init_arrays, ocean_parin, ocean_prognostic_equations,         &
     207           ocean_swap_timelevel, ocean_rrd_global, ocean_rrd_local,            &
     208           ocean_wrd_global, ocean_wrd_local, ocean_3d_data_averaging,         &
     209           stokes_drift_terms, wave_breaking_term
    194210
    195211
    … …  
    484500
    485501
    486     NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux,                 &
    487              langmuir_number, sa_surface, sa_vertical_gradient,                &
    488              sa_vertical_gradient_level, stokes_force, stokes_waveheight,      &
    489              stokes_wavelength, top_salinityflux, wall_salinityflux
     502    NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux, sa_surface,     &
     503             sa_vertical_gradient, sa_vertical_gradient_level,                 &
     504             stokes_waveheight, stokes_wavelength, top_salinityflux,           &
     505             wall_salinityflux, wave_breaking
    490506
    491507!
    … …  
    580596                        'top_salinityflux /= 0.0'
    581597       CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
     598    ENDIF
     599
     600!
     601!-- Check if Stokes force is to be used
     602    IF ( stokes_waveheight /= 0.0_wp  .AND.  stokes_wavelength /= 0.0_wp )  THEN
     603       stokes_force = .TRUE.
     604    ELSE
     605       IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength > 0.0_wp ) &
     606            .OR.                                                               &
     607            ( stokes_waveheight > 0.0_wp .AND. stokes_wavelength <= 0.0_wp ) &
     608            .OR.                                                               &
     609            ( stokes_waveheight < 0.0_wp .AND. stokes_wavelength < 0.0_wp  ) ) &
     610       THEN
     611          message_string = 'wrong settings for stokes_wavelength and/or ' //   &
     612                           'stokes_waveheight'
     613          CALL message( 'ocean_check_parameters', 'PA0460', 1, 2, 0, 6, 0 )
     614       ENDIF
    582615    ENDIF
    583616
    … …  
    10451078! Description:
    10461079! ------------
     1080!> Header output for ocean parameters
     1081!------------------------------------------------------------------------------!
     1082 SUBROUTINE ocean_header( io )
     1083
     1084
     1085    IMPLICIT NONE
     1086
     1087    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
     1088
     1089!
     1090!-- Write ocean header
     1091    WRITE( io, 1 )
     1092    IF ( stokes_force  )  WRITE( io, 2 )  stokes_waveheight, stokes_wavelength
     1093    IF ( wave_breaking )  THEN
     1094       WRITE( io, 3 )  alpha_wave_breaking, timescale_wave_breaking
     1095    ENDIF
     1096
     10971   FORMAT (//' Ocean settings:'/                                              &
     1098              ' ------------------------------------------'/)
     10992   FORMAT ('    --> Craik-Leibovich vortex force and Stokes drift switched',  &
     1100                     ' on'/                                                    &
     1101            '        waveheight: ',F4.1,' m   wavelength: ',F6.1,' m')
     11023   FORMAT ('    --> wave breaking generated turbulence switched on'/          &
     1103            '        alpha:    ',F4.1/                                         &
     1104            '        timescale:',F5.1,' s')
     1105
     1106 END SUBROUTINE ocean_header
     1107
     1108
     1109!------------------------------------------------------------------------------!
     1110! Description:
     1111! ------------
    10471112!> Allocate arrays and assign pointers.
    10481113!------------------------------------------------------------------------------!
    … …  
    10931158
    10941159    USE arrays_3d,                                                             &
    1095         ONLY:  dzu, hyp, pt_init, ref_state, zu, zw
     1160        ONLY:  dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw,    &
     1161               v_stokes_zu, v_stokes_zw, zu, zw
    10961162
    10971163    USE basic_constants_and_equations_mod,                                     &
    10981164        ONLY:  g
    10991165
     1166    USE basic_constants_and_equations_mod,                                     &
     1167        ONLY:  pi
     1168
    11001169    USE control_parameters,                                                    &
    11011170        ONLY:  initializing_actions, molecular_viscosity, rho_surface,         &
    1102                rho_reference, surface_pressure, use_single_reference_value
     1171               rho_reference, surface_pressure, top_momentumflux_u,            &
     1172               top_momentumflux_v, use_single_reference_value
    11031173
    11041174    USE indices,                                                               &
    … …  
    11071177    USE kinds
    11081178
    1109     USE pegrid
     1179    USE pegrid,                                                                &
     1180        ONLY:  myid
    11101181
    11111182    USE statistics,                                                            &
    … …  
    11191190    INTEGER(iwp) ::  n  !< loop index
    11201191
    1121     REAL(wp)     ::  dum   !< dummy argument
    1122     REAL(wp)     ::  pt_l  !< local scalar for pt used in equation of state function
    1123     REAL(wp)     ::  sa_l  !< local scalar for sa used in equation of state function
     1192    REAL(wp) ::  alpha !< angle of surface stress
     1193    REAL(wp) ::  dum   !< dummy argument
     1194    REAL(wp) ::  pt_l  !< local scalar for pt used in equation of state function
     1195    REAL(wp) ::  sa_l  !< local scalar for sa used in equation of state function
     1196    REAL(wp) ::  velocity_amplitude  !< local scalar for amplitude of Stokes drift velocity
     1197    REAL(wp) ::  x     !< temporary variable to store surface stress along x
     1198    REAL(wp) ::  y     !< temporary variable to store surface stress along y
    11241199
    11251200    REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init  !< local array for initial density
    … …  
    12701345    ENDIF
    12711346
     1347!
     1348!-- Calculate the Stokes drift velocity profile
     1349    IF ( stokes_force )  THEN
     1350
     1351!
     1352!--    First, calculate angle of surface stress
     1353       x = -top_momentumflux_u
     1354       y = -top_momentumflux_v
     1355       IF ( x == 0.0_wp )  THEN
     1356          IF ( y > 0.0_wp )  THEN
     1357             alpha = pi / 2.0_wp
     1358          ELSEIF ( y < 0.0_wp )  THEN
     1359             alpha = 3.0_wp * pi / 2.0_wp
     1360          ENDIF
     1361       ELSE
     1362          IF ( x < 0.0_wp )  THEN
     1363             alpha = ATAN( y / x ) + pi
     1364          ELSE
     1365             IF ( y < 0.0_wp )  THEN
     1366                alpha = ATAN( y / x ) + 2.0_wp * pi
     1367             ELSE
     1368                alpha = ATAN( y / x )
     1369             ENDIF
     1370          ENDIF
     1371       ENDIF
     1372
     1373       velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 *&
     1374                            SQRT( g * stokes_wavelength / ( 2.0_wp * pi ) )
     1375
     1376       DO  k = nzb, nzt
     1377          u_stokes_zu(k) = velocity_amplitude * COS( alpha ) *                 &
     1378                           EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
     1379          u_stokes_zw(k) = velocity_amplitude * COS( alpha ) *                 &
     1380                           EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
     1381          v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) *                 &
     1382                           EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
     1383          v_stokes_zw(k) = velocity_amplitude * SIN( alpha ) *                 &
     1384                           EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
     1385       ENDDO
     1386       u_stokes_zu(nzt+1) = u_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
     1387       u_stokes_zw(nzt+1) = u_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
     1388       v_stokes_zu(nzt+1) = v_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
     1389       v_stokes_zw(nzt+1) = v_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
     1390
     1391    ENDIF
     1392
     1393!
     1394!-- Wave breaking effects
     1395    IF ( wave_breaking )  THEN
     1396!
     1397!--    Calculate friction velocity at ocean surface
     1398       u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 +              &
     1399                                          top_momentumflux_v**2 ) )
     1400!
     1401!--    Set the time scale of random forcing. The vertical grid spacing at the
     1402!--    ocean surface is assumed as the length scale of turbulence.
     1403!--    Formula follows Noh et al. (2004), JPO
     1404       timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking /     &
     1405                                 u_star_wave_breaking
     1406!
     1407!--    Set random number seeds differently on the processor cores in order to
     1408!--    create different random number sequences
     1409       iran_ocean = iran_ocean + myid
     1410    ENDIF
    12721411
    12731412 END SUBROUTINE ocean_init
    … …  
    15991738          READ ( 13 )  bottom_salinityflux
    16001739
    1601        CASE ( 'langmuir_number' )
    1602           READ ( 13 )  langmuir_number
    1603 
    16041740       CASE ( 'sa_init' )
    16051741          READ ( 13 )  sa_init
    … …  
    16141750          READ ( 13 )  sa_vertical_gradient_level
    16151751
    1616        CASE ( 'stokes_force' )
    1617           READ ( 13 )  stokes_force
    1618 
    16191752       CASE ( 'stokes_waveheight' )
    16201753          READ ( 13 )  stokes_waveheight
    … …  
    16281761       CASE ( 'wall_salinityflux' )
    16291762          READ ( 13 )  wall_salinityflux
     1763
     1764       CASE ( 'wave_breaking' )
     1765          READ ( 13 )  wave_breaking
    16301766
    16311767       CASE DEFAULT
    … …  
    17311867    WRITE ( 14 )  bottom_salinityflux
    17321868
    1733     CALL wrd_write_string( 'langmuir_number' )
    1734     WRITE ( 14 )  langmuir_number
    1735 
    17361869    CALL wrd_write_string( 'sa_init' )
    17371870    WRITE ( 14 )  sa_init
    … …  
    17461879    WRITE ( 14 )  sa_vertical_gradient_level
    17471880
    1748     CALL wrd_write_string( 'stokes_force' )
    1749     WRITE ( 14 )  stokes_force
    1750 
    17511881    CALL wrd_write_string( 'stokes_waveheight' )
    17521882    WRITE ( 14 )  stokes_waveheight
    … …  
    17601890    CALL wrd_write_string( 'wall_salinityflux' )
    17611891    WRITE ( 14 )  wall_salinityflux
     1892
     1893    CALL wrd_write_string( 'wave_breaking' )
     1894    WRITE ( 14 )  wave_breaking
    17621895
    17631896 END SUBROUTINE ocean_wrd_global
    … …  
    17921925
    17931926
     1927!------------------------------------------------------------------------------!
     1928! Description:
     1929! ------------
     1930!> This routine calculates the Craik Leibovich vortex force and the additional
     1931!> effect of the Stokes drift on the Coriolis force
     1932!> Call for all gridpoints.
     1933!------------------------------------------------------------------------------!
     1934 SUBROUTINE stokes_drift_terms( component )
     1935
     1936    USE arrays_3d,                                                             &
     1937        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
     1938               v_stokes_zw, w, tend
     1939
     1940    USE control_parameters,                                                    &
     1941        ONLY:  f, fs, message_string
     1942
     1943    USE grid_variables,                                                        &
     1944        ONLY:  ddx, ddy
     1945
     1946    USE indices,                                                               &
     1947        ONLY:  nxl, nxr, nys, nysv, nyn, nzb, nzt
     1948
     1949    IMPLICIT NONE
     1950
     1951    INTEGER(iwp) ::  component  !< component of momentum equation
     1952    INTEGER(iwp) ::  i          !< loop index along x
     1953    INTEGER(iwp) ::  j          !< loop index along y
     1954    INTEGER(iwp) ::  k          !< loop index along z
     1955
     1956
     1957!
     1958!-- Compute Stokes terms for the respective velocity components
     1959    SELECT CASE ( component )
     1960
     1961!
     1962!--    u-component
     1963       CASE ( 1 )
     1964          DO  i = nxl, nxr
     1965             DO  j = nysv, nyn
     1966                DO  k = nzb+1, nzt
     1967                   tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (              &
     1968                                   0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)           &
     1969                                         + v(k,j,i)   - v(k,j,i-1)   ) * ddx   &
     1970                                 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy   &
     1971                                                                )              &
     1972                                 + f * v_stokes_zu(k)
     1973                ENDDO
     1974             ENDDO
     1975          ENDDO
     1976
     1977!
     1978!--    v-component
     1979       CASE ( 2 )
     1980          DO  i = nxl, nxr
     1981             DO  j = nysv, nyn
     1982                DO  k = nzb+1, nzt
     1983                   tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (              &
     1984                                   0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx   &
     1985                                 - 0.5 * ( u(k,j,i) - u(k,j-1,i)               &
     1986                                         + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
     1987                                                                )              &
     1988                                 - f * u_stokes_zu(k)
     1989                ENDDO
     1990             ENDDO
     1991          ENDDO
     1992
     1993!
     1994!--    w-component
     1995       CASE ( 3 )
     1996          DO  i = nxl, nxr
     1997             DO  j = nys, nyn
     1998                DO  k = nzb+1, nzt
     1999                   tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
     2000                                             0.5 * ( u(k+1,j,i) - u(k,j,i)     &
     2001                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
     2002                                                   ) * ddzu(k+1)               &
     2003                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
     2004                                                   ) * ddx      )              &
     2005                                             - v_stokes_zw(k) * (              &
     2006                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
     2007                                                   ) * ddy                     &
     2008                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
     2009                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
     2010                                                   ) * ddzu(k)  )              &
     2011                                           + fs * u_stokes_zw(k)
     2012                ENDDO
     2013             ENDDO
     2014          ENDDO
     2015
     2016       CASE DEFAULT
     2017          WRITE( message_string, * ) 'wrong component of Stokes force: ',      &
     2018                                     component
     2019          CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
     2020
     2021    END SELECT
     2022
     2023 END SUBROUTINE stokes_drift_terms
     2024
     2025
     2026!------------------------------------------------------------------------------!
     2027! Description:
     2028! ------------
     2029!> This routine calculates the Craik Leibovich vortex force and the additional
     2030!> effect of the Stokes drift on the Coriolis force
     2031!> Call for gridpoints i,j.
     2032!------------------------------------------------------------------------------!
     2033
     2034 SUBROUTINE stokes_drift_terms_ij( i, j, component )
     2035
     2036    USE arrays_3d,                                                             &
     2037        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
     2038               v_stokes_zw, w, tend
     2039
     2040    USE control_parameters,                                                    &
     2041        ONLY:  f, fs, message_string
     2042
     2043    USE grid_variables,                                                        &
     2044        ONLY:  ddx, ddy
     2045
     2046    USE indices,                                                               &
     2047        ONLY:  nxl, nxr, nys, nysv, nyn, nzb, nzt
     2048
     2049    IMPLICIT NONE
     2050
     2051    INTEGER(iwp) ::  component  !< component of momentum equation
     2052    INTEGER(iwp) ::  i          !< loop index along x
     2053    INTEGER(iwp) ::  j          !< loop index along y
     2054    INTEGER(iwp) ::  k          !< loop incex along z
     2055
     2056
     2057!
     2058!-- Compute Stokes terms for the respective velocity components
     2059    SELECT CASE ( component )
     2060
     2061!
     2062!--    u-component
     2063       CASE ( 1 )
     2064          DO  k = nzb+1, nzt
     2065             tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (                    &
     2066                                     0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)         &
     2067                                           + v(k,j,i)   - v(k,j,i-1)   ) * ddx &
     2068                                   - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy &
     2069                                                          )                    &
     2070                                       + f * v_stokes_zu(k)
     2071          ENDDO
     2072!
     2073!--    v-component
     2074       CASE ( 2 )
     2075          DO  k = nzb+1, nzt
     2076             tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (                    &
     2077                                     0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx &
     2078                                   - 0.5 * ( u(k,j,i) - u(k,j-1,i)             &
     2079                                           + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
     2080                                                          )                    &
     2081                                       - f * u_stokes_zu(k)
     2082          ENDDO
     2083
     2084!
     2085!--    w-component
     2086       CASE ( 3 )
     2087          DO  k = nzb+1, nzt
     2088             tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
     2089                                     0.5 * ( u(k+1,j,i) - u(k,j,i)     &
     2090                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
     2091                                                   ) * ddzu(k+1)               &
     2092                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
     2093                                                   ) * ddx )                   &
     2094                                       - v_stokes_zw(k) * (                    &
     2095                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
     2096                                                   ) * ddy                     &
     2097                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
     2098                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
     2099                                                   ) * ddzu(k)  )              &
     2100                                       + fs * u_stokes_zw(k)
     2101          ENDDO
     2102
     2103       CASE DEFAULT
     2104          WRITE( message_string, * ) ' wrong component: ', component
     2105          CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
     2106
     2107    END SELECT
     2108
     2109 END SUBROUTINE stokes_drift_terms_ij
     2110
     2111
     2112!------------------------------------------------------------------------------!
     2113! Description:
     2114! ------------
     2115!> This routine calculates turbulence generated by wave breaking near the ocean
     2116!> surface, following a parameterization given in Noh et al. (2004), JPO
     2117!> Call for all gridpoints.
     2118!> TODO: so far, this routine only works if the model time step has about the
     2119!>       same value as the time scale of wave breaking!
     2120!------------------------------------------------------------------------------!
     2121 SUBROUTINE wave_breaking_term( component )
     2122
     2123    USE arrays_3d,                                                             &
     2124        ONLY:  u_p, v_p
     2125
     2126    USE control_parameters,                                                    &
     2127        ONLY:  dt_3d, message_string
     2128
     2129    USE indices,                                                               &
     2130        ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzt
     2131
     2132    IMPLICIT NONE
     2133
     2134    INTEGER(iwp) ::  component  !< component of momentum equation
     2135    INTEGER(iwp) ::  i          !< loop index along x
     2136    INTEGER(iwp) ::  j          !< loop index along y
     2137
     2138    REAL(wp) ::  random_gauss  !< function that creates a random number with a
     2139                               !< Gaussian distribution
     2140
     2141
     2142!
     2143!-- Compute wave breaking terms for the respective velocity components.
     2144!-- Velocities are directly manipulated, since this is not a real force
     2145    SELECT CASE ( component )
     2146
     2147!
     2148!--    u-component
     2149       CASE ( 1 )
     2150          DO  i = nxlu, nxr
     2151             DO  j = nys, nyn
     2152                u_p(nzt,j,i) = u_p(nzt,j,i) +                                  &
     2153                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
     2154                               * alpha_wave_breaking * u_star_wave_breaking    &
     2155                               / timescale_wave_breaking * dt_3d
     2156             ENDDO
     2157          ENDDO
     2158!
     2159!--    v-component
     2160       CASE ( 2 )
     2161          DO  i = nxl, nxr
     2162             DO  j = nysv, nyn
     2163                v_p(nzt,j,i) = v_p(nzt,j,i) +                                  &
     2164                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
     2165                               * alpha_wave_breaking * u_star_wave_breaking    &
     2166                               / timescale_wave_breaking * dt_3d
     2167             ENDDO
     2168          ENDDO
     2169
     2170       CASE DEFAULT
     2171          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
     2172                                     component
     2173          CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
     2174
     2175    END SELECT
     2176
     2177 END SUBROUTINE wave_breaking_term
     2178
     2179
     2180!------------------------------------------------------------------------------!
     2181! Description:
     2182! ------------
     2183!> This routine calculates turbulence generated by wave breaking near the ocean
     2184!> surface, following a parameterization given in Noh et al. (2004), JPO
     2185!> Call for gridpoint i,j.
     2186!> TODO: so far, this routine only works if the model time step has about the
     2187!>       same value as the time scale of wave breaking!
     2188!------------------------------------------------------------------------------!
     2189 SUBROUTINE wave_breaking_term_ij( i, j, component )
     2190
     2191    USE arrays_3d,                                                             &
     2192        ONLY:  u_p, v_p
     2193
     2194    USE control_parameters,                                                    &
     2195        ONLY:  dt_3d, message_string
     2196
     2197    USE indices,                                                               &
     2198        ONLY:  nzt
     2199
     2200    IMPLICIT NONE
     2201
     2202    INTEGER(iwp) ::  component  !< component of momentum equation
     2203    INTEGER(iwp) ::  i          !< loop index along x
     2204    INTEGER(iwp) ::  j          !< loop index along y
     2205
     2206    REAL(wp) ::  random_gauss  !< function that creates a random number with a
     2207                               !< Gaussian distribution
     2208
     2209!
     2210!-- Compute wave breaking terms for the respective velocity components
     2211    SELECT CASE ( component )
     2212
     2213!
     2214!--    u-/v-component
     2215       CASE ( 1 )
     2216          u_p(nzt,j,i) = u_p(nzt,j,i) +                                        &
     2217                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
     2218                         * alpha_wave_breaking * u_star_wave_breaking          &
     2219                         / timescale_wave_breaking * dt_3d
     2220
     2221       CASE ( 2 )
     2222          v_p(nzt,j,i) = v_p(nzt,j,i) +                                        &
     2223                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
     2224                         * alpha_wave_breaking * u_star_wave_breaking          &
     2225                         / timescale_wave_breaking * dt_3d
     2226
     2227       CASE DEFAULT
     2228          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
     2229                                     component
     2230          CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
     2231
     2232    END SELECT
     2233
     2234 END SUBROUTINE wave_breaking_term_ij
     2235
     2236
    17942237 END MODULE ocean_mod
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |