Changeset 3302 for palm/trunk/SOURCE/ocean_mod.f90
 Timestamp:
 Oct 3, 2018 2:39:40 AM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/ocean_mod.f90
r3294 r3302 25 25 !  26 26 ! $Id$ 27 ! Craik Leibovich force (Stokes drift) + wave breaking effect added 28 ! 29 ! 3294 20181001 02:37:10Z raasch 27 30 ! initial revision 28 ! 29 ! 30 ! 31 ! 31 32 ! 32 33 ! Authors: … … 36 37 ! Description: 37 38 !  38 !> This module contains relevant code for PALM's ocean option, e.g. the39 !> This module contains relevant code for PALM's ocean mode, e.g. the 39 40 !> 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 42 42 !! 43 43 MODULE ocean_mod … … 66 66 67 67 INTEGER(iwp) :: ibc_sa_t !< integer flag for bc_sa_t 68 INTEGER(iwp) :: iran_ocean = 1234567 !< random number used for wave breaking effects 68 69 69 70 INTEGER(iwp) :: sa_vertical_gradient_level_ind(10) = 9999 !< grid index values of sa_vertical_gradient_level(s) 70 71 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 74 76 REAL(wp) :: prho_reference !< reference state of potential density at ocean surface 75 77 REAL(wp) :: sa_surface = 35.0_wp !< surface salinity, namelist parameter 76 78 REAL(wp) :: sa_vertical_gradient(10) = 0.0_wp !< namelist parameter 77 79 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 80 84 81 85 REAL(wp), DIMENSION(12), PARAMETER :: nom = & … … 100 104 SAVE 101 105 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, & 103 107 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 106 109 107 110 … … 139 142 END INTERFACE ocean_data_output_3d 140 143 144 INTERFACE ocean_header 145 MODULE PROCEDURE ocean_header 146 END INTERFACE ocean_header 147 141 148 INTERFACE ocean_init 142 149 MODULE PROCEDURE ocean_init … … 180 187 END INTERFACE ocean_3d_data_averaging 181 188 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 183 198 184 199 PRIVATE … … 188 203 ocean_check_data_output_pr, ocean_check_parameters, & 189 204 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 194 210 195 211 … … 484 500 485 501 486 NAMELIST /ocean_parameters/ bc_sa_t, bottom_salinityflux, 487 langmuir_number, sa_surface, sa_vertical_gradient,&488 s a_vertical_gradient_level, stokes_force, stokes_waveheight,&489 stokes_wavelength, top_salinityflux, wall_salinityflux502 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 490 506 491 507 ! … … 580 596 'top_salinityflux /= 0.0' 581 597 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 582 615 ENDIF 583 616 … … 1045 1078 ! Description: 1046 1079 !  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 1097 1 FORMAT (//' Ocean settings:'/ & 1098 ' '/) 1099 2 FORMAT (' > CraikLeibovich vortex force and Stokes drift switched', & 1100 ' on'/ & 1101 ' waveheight: ',F4.1,' m wavelength: ',F6.1,' m') 1102 3 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 !  1047 1112 !> Allocate arrays and assign pointers. 1048 1113 !! … … 1093 1158 1094 1159 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 1096 1162 1097 1163 USE basic_constants_and_equations_mod, & 1098 1164 ONLY: g 1099 1165 1166 USE basic_constants_and_equations_mod, & 1167 ONLY: pi 1168 1100 1169 USE control_parameters, & 1101 1170 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 1103 1173 1104 1174 USE indices, & … … 1107 1177 USE kinds 1108 1178 1109 USE pegrid 1179 USE pegrid, & 1180 ONLY: myid 1110 1181 1111 1182 USE statistics, & … … 1119 1190 INTEGER(iwp) :: n !< loop index 1120 1191 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 1124 1199 1125 1200 REAL(wp), DIMENSION(nzb:nzt+1) :: rho_ocean_init !< local array for initial density … … 1270 1345 ENDIF 1271 1346 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 1272 1411 1273 1412 END SUBROUTINE ocean_init … … 1599 1738 READ ( 13 ) bottom_salinityflux 1600 1739 1601 CASE ( 'langmuir_number' )1602 READ ( 13 ) langmuir_number1603 1604 1740 CASE ( 'sa_init' ) 1605 1741 READ ( 13 ) sa_init … … 1614 1750 READ ( 13 ) sa_vertical_gradient_level 1615 1751 1616 CASE ( 'stokes_force' )1617 READ ( 13 ) stokes_force1618 1619 1752 CASE ( 'stokes_waveheight' ) 1620 1753 READ ( 13 ) stokes_waveheight … … 1628 1761 CASE ( 'wall_salinityflux' ) 1629 1762 READ ( 13 ) wall_salinityflux 1763 1764 CASE ( 'wave_breaking' ) 1765 READ ( 13 ) wave_breaking 1630 1766 1631 1767 CASE DEFAULT … … 1731 1867 WRITE ( 14 ) bottom_salinityflux 1732 1868 1733 CALL wrd_write_string( 'langmuir_number' )1734 WRITE ( 14 ) langmuir_number1735 1736 1869 CALL wrd_write_string( 'sa_init' ) 1737 1870 WRITE ( 14 ) sa_init … … 1746 1879 WRITE ( 14 ) sa_vertical_gradient_level 1747 1880 1748 CALL wrd_write_string( 'stokes_force' )1749 WRITE ( 14 ) stokes_force1750 1751 1881 CALL wrd_write_string( 'stokes_waveheight' ) 1752 1882 WRITE ( 14 ) stokes_waveheight … … 1760 1890 CALL wrd_write_string( 'wall_salinityflux' ) 1761 1891 WRITE ( 14 ) wall_salinityflux 1892 1893 CALL wrd_write_string( 'wave_breaking' ) 1894 WRITE ( 14 ) wave_breaking 1762 1895 1763 1896 END SUBROUTINE ocean_wrd_global … … 1792 1925 1793 1926 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 ! ucomponent 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,i1) & 1969 + v(k,j,i)  v(k,j,i1) ) * ddx & 1970  0.5 * ( u(k,j+1,i)  u(k,j1,i) ) * ddy & 1971 ) & 1972 + f * v_stokes_zu(k) 1973 ENDDO 1974 ENDDO 1975 ENDDO 1976 1977 ! 1978 ! vcomponent 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,i1) ) * ddx & 1985  0.5 * ( u(k,j,i)  u(k,j1,i) & 1986 + u(k,j,i+1)  u(k,j1,i+1) ) * ddy & 1987 ) & 1988  f * u_stokes_zu(k) 1989 ENDDO 1990 ENDDO 1991 ENDDO 1992 1993 ! 1994 ! wcomponent 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,i1) & 2004 ) * ddx ) & 2005  v_stokes_zw(k) * ( & 2006 0.5 * ( w(k,j+1,i)  w(k,j1,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 ! ucomponent 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,i1) & 2067 + v(k,j,i)  v(k,j,i1) ) * ddx & 2068  0.5 * ( u(k,j+1,i)  u(k,j1,i) ) * ddy & 2069 ) & 2070 + f * v_stokes_zu(k) 2071 ENDDO 2072 ! 2073 ! vcomponent 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,i1) ) * ddx & 2078  0.5 * ( u(k,j,i)  u(k,j1,i) & 2079 + u(k,j,i+1)  u(k,j1,i+1) ) * ddy & 2080 ) & 2081  f * u_stokes_zu(k) 2082 ENDDO 2083 2084 ! 2085 ! wcomponent 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,i1) & 2093 ) * ddx ) & 2094  v_stokes_zw(k) * ( & 2095 0.5 * ( w(k,j+1,i)  w(k,j1,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 ! ucomponent 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 ! vcomponent 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/vcomponent 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 1794 2237 END MODULE ocean_mod
Note: See TracChangeset
for help on using the changeset viewer.