Changeset 3302


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

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

Location:
palm/trunk/SOURCE
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_s_pw.f90

    r2718 r3302  
    2525! -----------------
    2626! $Id$
     27! Stokes drift velocity added
     28!
     29! 2718 2018-01-02 08:49:38Z maronga
    2730! Corrected "Former revisions" section
    2831!
     
    104107
    105108       USE arrays_3d,                                                          &
    106            ONLY:  dd2zu, tend, u, v, w
     109           ONLY:  dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w
    107110
    108111       USE control_parameters,                                                 &
     
    124127       INTEGER(iwp) ::  j !<
    125128       INTEGER(iwp) ::  k !<
     129
     130       REAL(wp)     ::  gu  !< local additional advective velocity
     131       REAL(wp)     ::  gv  !< local additional advective velocity
    126132
    127133#if defined( __nopointer )
     
    135141          DO  j = nys, nyn
    136142             DO  k = nzb+1, nzt
    137                 tend(k,j,i) = tend(k,j,i) +                                       &
    138             ( -0.5_wp * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
    139                         - ( u(k,j,i)   - u_gtrans ) * ( sk(k,j,i-1) - sk(k,j,i) ) &
    140                         ) * ddx                                                   &
    141               -0.5_wp * ( ( v(k,j+1,i) - v_gtrans ) * ( sk(k,j+1,i) - sk(k,j,i) ) &
    142                         - ( v(k,j,i)   - v_gtrans ) * ( sk(k,j-1,i) - sk(k,j,i) ) &
    143                         ) * ddy                                                   &
    144               -         (   w(k,j,i)                * ( sk(k+1,j,i) - sk(k,j,i) ) &
    145                         -   w(k-1,j,i)              * ( sk(k-1,j,i) - sk(k,j,i) ) &
    146                         ) * dd2zu(k)                                              &
    147             ) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     143
     144!
     145!--             Galilean transformation + Stokes drift velocity (ocean only)
     146                gu = u_gtrans - u_stokes_zu(k)
     147                gv = v_gtrans - v_stokes_zu(k)
     148
     149                tend(k,j,i) = tend(k,j,i) +                                    &
     150                                     ( -0.5_wp * ( ( u(k,j,i+1) - gu       ) * &
     151                                                   ( sk(k,j,i+1) - sk(k,j,i) ) &
     152                                                 - ( u(k,j,i)   - gu       ) * &
     153                                                   ( sk(k,j,i-1) - sk(k,j,i) ) &
     154                                                 ) * ddx                       &
     155                                       -0.5_wp * ( ( v(k,j+1,i) - gv       ) * &
     156                                                   ( sk(k,j+1,i) - sk(k,j,i) ) &
     157                                                 - ( v(k,j,i)   - gv       ) * &
     158                                                   ( sk(k,j-1,i) - sk(k,j,i) ) &
     159                                                 ) * ddy                       &
     160                                       -         (   w(k,j,i)                * &
     161                                                   ( sk(k+1,j,i) - sk(k,j,i) ) &
     162                                                 -   w(k-1,j,i)              * &
     163                                                   ( sk(k-1,j,i) - sk(k,j,i) ) &
     164                                                 ) * dd2zu(k)                  &
     165                                      ) * MERGE( 1.0_wp, 0.0_wp,               &
     166                                               BTEST( wall_flags_0(k,j,i), 0 ) )
    148167             ENDDO
    149168          ENDDO
     
    161180
    162181       USE arrays_3d,                                                          &
    163            ONLY:  dd2zu, tend, u, v, w
     182           ONLY:  dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w
    164183
    165184       USE control_parameters,                                                 &
     
    180199       INTEGER(iwp) ::  j !<
    181200       INTEGER(iwp) ::  k !<
     201
     202       REAL(wp)     ::  gu  !< local additional advective velocity
     203       REAL(wp)     ::  gv  !< local additional advective velocity
    182204
    183205#if defined( __nopointer )
     
    189211
    190212       DO  k = nzb+1, nzt
    191           tend(k,j,i) = tend(k,j,i) +                                             &
    192             ( -0.5_wp * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
    193                         - ( u(k,j,i)   - u_gtrans ) * ( sk(k,j,i-1) - sk(k,j,i) ) &
    194                         ) * ddx                                                   &
    195               -0.5_wp * ( ( v(k,j+1,i) - v_gtrans ) * ( sk(k,j+1,i) - sk(k,j,i) ) &
    196                         - ( v(k,j,i)   - v_gtrans ) * ( sk(k,j-1,i) - sk(k,j,i) ) &
    197                         ) * ddy                                                   &
    198               -         (   w(k,j,i)                * ( sk(k+1,j,i) - sk(k,j,i) ) &
    199                         -   w(k-1,j,i)              * ( sk(k-1,j,i) - sk(k,j,i) ) &
    200                         ) * dd2zu(k)                                              &
    201             ) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     213
     214!
     215!--       Galilean transformation + Stokes drift velocity (ocean only)
     216          gu = u_gtrans - u_stokes_zu(k)
     217          gv = v_gtrans - v_stokes_zu(k)
     218
     219          tend(k,j,i) = tend(k,j,i) +                                          &
     220                                    ( -0.5_wp * ( ( u(k,j,i+1) - gu        ) * &
     221                                                  ( sk(k,j,i+1) - sk(k,j,i) )  &
     222                                                - ( u(k,j,i)   - gu        ) * &
     223                                                  ( sk(k,j,i-1) - sk(k,j,i) )  &
     224                                                ) * ddx                        &
     225                                      -0.5_wp * ( ( v(k,j+1,i) - gv       ) *  &
     226                                                  ( sk(k,j+1,i) - sk(k,j,i) )  &
     227                                                - ( v(k,j,i)   - gv       ) *  &
     228                                                  ( sk(k,j-1,i) - sk(k,j,i) )  &
     229                                                ) * ddy                        &
     230                                      -         (   w(k,j,i)                *  &
     231                                                  ( sk(k+1,j,i) - sk(k,j,i) )  &
     232                                                -   w(k-1,j,i)              *  &
     233                                                  ( sk(k-1,j,i) - sk(k,j,i) )  &
     234                                                ) * dd2zu(k)                   &
     235                                    ) * MERGE( 1.0_wp, 0.0_wp,                 &
     236                                               BTEST( wall_flags_0(k,j,i), 0 ) )
    202237       ENDDO
    203238
  • palm/trunk/SOURCE/advec_ws.f90

    r3298 r3302  
    2525! -----------------
    2626! $Id$
     27! Stokes drift velocity added in scalar advection
     28!
     29! 3298 2018-10-02 12:21:11Z kanani
    2730! Add formatted note from ketelsen about statistics for chemistry
    2831!
     
    11471150
    11481151       USE arrays_3d,                                                          &
    1149            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
     1152           ONLY:  ddzw, drho_air, rho_air_zw, tend, u, u_stokes_zu, v,         &
     1153                  v_stokes_zu, w
    11501154
    11511155       USE control_parameters,                                                 &
     
    12271231             ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
    12281232
    1229              v_comp                  = v(k,j,i) - v_gtrans
     1233             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    12301234             swap_flux_y_local(k,tn) = v_comp *         (                     &
    12311235                                               ( 37.0_wp * ibit5 * adv_sca_5  &
     
    12631267          DO  k = nzb_max+1, nzt
    12641268
    1265              v_comp                  = v(k,j,i) - v_gtrans
     1269             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    12661270             swap_flux_y_local(k,tn) = v_comp * (                             &
    12671271                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )   &
     
    12881292             ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
    12891293
    1290              u_comp                     = u(k,j,i) - u_gtrans
     1294             u_comp                     = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    12911295             swap_flux_x_local(k,j,tn) = u_comp * (                           &
    12921296                                               ( 37.0_wp * ibit2 * adv_sca_5  &
     
    13231327          DO  k = nzb_max+1, nzt
    13241328
    1325              u_comp                 = u(k,j,i) - u_gtrans
     1329             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    13261330             swap_flux_x_local(k,j,tn) = u_comp * (                           &
    13271331                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
     
    13551359          ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
    13561360
    1357           u_comp    = u(k,j,i+1) - u_gtrans
     1361          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
    13581362          flux_r(k) = u_comp * (                                              &
    13591363                     ( 37.0_wp * ibit2 * adv_sca_5                            &
     
    13901394          ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
    13911395
    1392           v_comp    = v(k,j+1,i) - v_gtrans
     1396          v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
    13931397          flux_n(k) = v_comp * (                                              &
    13941398                     ( 37.0_wp * ibit5 * adv_sca_5                            &
     
    15141518       DO  k = nzb_max+1, nzt
    15151519
    1516           u_comp    = u(k,j,i+1) - u_gtrans
     1520          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
    15171521          flux_r(k) = u_comp * (                                              &
    15181522                      37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
     
    15241528                    +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    15251529
    1526           v_comp    = v(k,j+1,i) - v_gtrans
     1530          v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
    15271531          flux_n(k) = v_comp * (                                              &
    15281532                      37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
     
    31923196    SUBROUTINE advec_s_ws( sk, sk_char )
    31933197
    3194        USE arrays_3d,                                                         &
    3195            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
    3196 
    3197        USE control_parameters,                                                &
     3198       USE arrays_3d,                                                          &
     3199           ONLY:  ddzw, drho_air, rho_air_zw, tend, u, u_stokes_zu, v,         &
     3200                  v_stokes_zu, w
     3201
     3202       USE control_parameters,                                                 &
    31983203           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    31993204
    3200        USE grid_variables,                                                    &
     3205       USE grid_variables,                                                     &
    32013206           ONLY:  ddx, ddy
    32023207
    3203        USE indices,                                                           &
    3204            ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max,   &
     3208       USE indices,                                                            &
     3209           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max,    &
    32053210                  nzt, advc_flags_1
    32063211           
    32073212       USE kinds
    32083213       
    3209        USE statistics,                                                        &
    3210            ONLY:  hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,      &
    3211                   sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l,          &
     3214       USE statistics,                                                         &
     3215           ONLY:  hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,       &
     3216                  sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l,           &
    32123217                  sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep 
    32133218                 
     
    32723277             ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
    32733278
    3274              u_comp                 = u(k,j,i) - u_gtrans
     3279             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    32753280             swap_flux_x_local(k,j) = u_comp * (                              &
    32763281                                             ( 37.0_wp * ibit2 * adv_sca_5    &
     
    33073312          DO  k = nzb_max+1, nzt
    33083313
    3309              u_comp                 = u(k,j,i) - u_gtrans
     3314             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    33103315             swap_flux_x_local(k,j) = u_comp * (                              &
    33113316                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
     
    33333338             ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
    33343339
    3335              v_comp               = v(k,j,i) - v_gtrans
     3340             v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    33363341             swap_flux_y_local(k) = v_comp * (                                &
    33373342                                             ( 37.0_wp * ibit5 * adv_sca_5    &
     
    33693374          DO  k = nzb_max+1, nzt
    33703375
    3371              v_comp               = v(k,j,i) - v_gtrans
     3376             v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    33723377             swap_flux_y_local(k) = v_comp * (                               &
    33733378                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )  &
     
    33753380                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )  &
    33763381                                             ) * adv_sca_5
    3377               swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
     3382             swap_diss_y_local(k) = -ABS( v_comp ) * (                       &
    33783383                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )  &
    33793384                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )  &
     
    33963401                ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
    33973402
    3398                 u_comp    = u(k,j,i+1) - u_gtrans
     3403                u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
    33993404                flux_r(k) = u_comp * (                                        &
    34003405                          ( 37.0_wp * ibit2 * adv_sca_5                       &
     
    34313436                ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
    34323437
    3433                 v_comp    = v(k,j+1,i) - v_gtrans
     3438                v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
    34343439                flux_n(k) = v_comp * (                                        &
    34353440                          ( 37.0_wp * ibit5 * adv_sca_5                       &
     
    35483553             DO  k = nzb_max+1, nzt
    35493554
    3550                 u_comp    = u(k,j,i+1) - u_gtrans
     3555                u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
    35513556                flux_r(k) = u_comp * (                                        &
    35523557                      37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
     
    35583563                    +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    35593564
    3560                 v_comp    = v(k,j+1,i) - v_gtrans
     3565                v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
    35613566                flux_n(k) = v_comp * (                                        &
    35623567                      37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
  • palm/trunk/SOURCE/check_parameters.f90

    r3301 r3302  
    2525! -----------------
    2626! $Id$
     27! small code rearrangements
     28!
     29! 3301 2018-10-02 17:10:50Z gronemeier
    2730! corrected former revision section
    2831!
     
    14681471
    14691472!-- Check the module settings
    1470     IF ( ocean_mode )  CALL ocean_check_parameters
     1473    IF ( bulk_cloud_model )     CALL bcm_check_parameters
     1474    IF ( gust_module_enabled )  CALL gust_check_parameters
     1475    IF ( large_scale_forcing  .OR.  nudging )                                  &
     1476                                CALL lsf_nudging_check_parameters
     1477    IF ( land_surface )         CALL lsm_check_parameters
     1478    IF ( ocean_mode )           CALL ocean_check_parameters
     1479    IF ( plant_canopy )         CALL pcm_check_parameters
     1480    IF ( radiation )            CALL radiation_check_parameters
     1481    IF ( calculate_spectra )    CALL spectra_check_parameters
     1482    CALL stg_check_parameters
    14711483    CALL tcm_check_parameters
    1472     CALL stg_check_parameters
    1473     IF ( plant_canopy )         CALL pcm_check_parameters
    1474     IF ( calculate_spectra )    CALL spectra_check_parameters
    1475     IF ( land_surface )         CALL lsm_check_parameters
    1476     IF ( bulk_cloud_model )     CALL bcm_check_parameters
    14771484    IF ( urban_surface )        CALL usm_check_parameters
    14781485    IF ( wind_turbine )         CALL wtm_check_parameters
    1479     IF ( radiation )            CALL radiation_check_parameters
    1480     IF ( gust_module_enabled )  CALL gust_check_parameters
    1481     IF ( large_scale_forcing  .OR.  nudging )  CALL lsf_nudging_check_parameters
    14821486
    14831487!
  • palm/trunk/SOURCE/header.f90

    r3298 r3302  
    2525! -----------------
    2626! $Id$
     27! call of ocean_header
     28!
     29! 3298 2018-10-02 12:21:11Z kanani
    2730! Minor formatting (kanani)
    2831! Add chemistry header (basit)
     
    437440
    438441    USE ocean_mod,                                                             &
    439         ONLY:  ibc_sa_t, prho_reference, sa_surface, sa_vertical_gradient,     &
    440                sa_vertical_gradient_level, sa_vertical_gradient_level_ind
     442        ONLY:  ibc_sa_t, ocean_header, prho_reference, sa_surface,             &
     443               sa_vertical_gradient, sa_vertical_gradient_level,               &
     444               sa_vertical_gradient_level_ind
    441445
    442446    USE particle_attributes,                                                   &
     
    11201124    ENDIF
    11211125
    1122     IF ( syn_turb_gen )  CALL stg_header ( io )
    1123 
    1124     IF ( plant_canopy )  CALL pcm_header ( io )
    1125 
    1126     IF ( land_surface )  CALL lsm_header ( io )
    1127 
    1128     IF ( radiation )  CALL radiation_header ( io )
    1129 
    1130     IF ( gust_module_enabled )  CALL gust_header ( io )
     1126!
     1127!-- Header information from other modules
     1128    IF ( gust_module_enabled )  CALL gust_header( io )
     1129    IF ( land_surface        )  CALL lsm_header( io )
     1130    IF ( ocean_mode          )  CALL ocean_header( io )
     1131    IF ( plant_canopy        )  CALL pcm_header( io )
     1132    IF ( radiation           )  CALL radiation_header( io )
     1133    IF ( syn_turb_gen        )  CALL stg_header( io )
    11311134
    11321135    IF ( air_chemistry )  CALL chem_header ( io )
  • palm/trunk/SOURCE/init_3d_model.f90

    r3298 r3302  
    2525! -----------------
    2626! $Id$
     27! allocate and set stokes drift velocity profiles
     28!
     29! 3298 2018-10-02 12:21:11Z kanani
    2730! Minor formatting (kanani)
    2831! Added CALL to the chem_emissions module (Russo)
     
    788791#endif
    789792    ENDIF
     793
     794!
     795!-- Allocate and set 1d-profiles for Stokes drift velocity. It may be set to
     796!-- non-zero values later in ocean_init
     797    ALLOCATE( u_stokes_zu(nzb:nzt+1), u_stokes_zw(nzb:nzt+1),                  &
     798              v_stokes_zu(nzb:nzt+1), v_stokes_zw(nzb:nzt+1) )
     799    u_stokes_zu(:) = 0.0_wp
     800    u_stokes_zw(:) = 0.0_wp
     801    v_stokes_zu(:) = 0.0_wp
     802    v_stokes_zw(:) = 0.0_wp
    790803
    791804!
  • palm/trunk/SOURCE/modules.f90

    r3298 r3302  
    2525! -----------------
    2626! $Id$
     27! +u/v_stokes_zu, u/v_stokes_zw
     28!
     29! 3298 2018-10-02 12:21:11Z kanani
    2730! Minor formatting/clean-up (kanani)
    2831! Added some variables for time treatment (Russo)
     
    739742    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ug                     !< geostrophic wind component in x-direction
    740743    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_init                 !< initial profile of horizontal velocity component u
     744    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_stokes_zu            !< u-component of Stokes drift velocity at zu levels
     745    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_stokes_zw            !< u-component of Stokes drift velocity at zw levels
    741746    REAL(wp), DIMENSION(:), ALLOCATABLE ::  vg                     !< geostrophic wind component in y-direction
    742747    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_init                 !< initial profile of horizontal velocity component v
     748    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_stokes_zu            !< v-component of Stokes drift velocity at zu levels
     749    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_stokes_zw            !< v-component of Stokes drift velocity at zw levels
    743750    REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_subs                 !< subsidence/ascent velocity
    744751    REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu                     !< vertical grid coordinate of u-grid (in m)
  • 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
  • palm/trunk/SOURCE/prognostic_equations.f90

    r3298 r3302  
    2525! -----------------
    2626! $Id$
     27! Stokes drift + wave breaking term added
     28!
     29! 3298 2018-10-02 12:21:11Z kanani
    2730! Code added for decycling chemistry (basit)
    2831!
     
    392395
    393396    USE ocean_mod,                                                             &
    394         ONLY:  eqn_state_seawater, ocean_prognostic_equations
     397        ONLY:  ocean_prognostic_equations, stokes_drift_terms, stokes_force,   &
     398               wave_breaking, wave_breaking_term
    395399
    396400    USE plant_canopy_model_mod,                                                &
     
    595599
    596600!
     601!--          Effect of Stokes drift (in ocean mode only)
     602             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
     603
     604!
    597605!--          Forces by wind turbines
    598606             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 1 )
     
    612620                                                 )
    613621             ENDDO
     622
     623!
     624!--          Add turbulence generated by wave breaking (in ocean mode only)
     625             IF ( wave_breaking  .AND.                                         &
     626               intermediate_timestep_count == intermediate_timestep_count_max )&
     627             THEN
     628                CALL wave_breaking_term( i, j, 1 )
     629             ENDIF
    614630
    615631!
     
    666682
    667683!
     684!--          Effect of Stokes drift (in ocean mode only)
     685             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
     686
     687!
    668688!--          Forces by wind turbines
    669689             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 2 )
     
    682702                                                 )
    683703             ENDDO
     704
     705!
     706!--          Add turbulence generated by wave breaking (in ocean mode only)
     707             IF ( wave_breaking  .AND.                                         &
     708               intermediate_timestep_count == intermediate_timestep_count_max )&
     709             THEN
     710                CALL wave_breaking_term( i, j, 2 )
     711             ENDIF
    684712
    685713!
     
    731759!--       Drag by plant canopy
    732760          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
     761
     762!
     763!--       Effect of Stokes drift (in ocean mode only)
     764          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
    733765
    734766!
     
    13831415
    13841416!
     1417!-- Effect of Stokes drift (in ocean mode only)
     1418    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
     1419
     1420!
    13851421!-- Forces by wind turbines
    13861422    IF ( wind_turbine )  CALL wtm_tendencies( 1 )
     
    14051441
    14061442!
     1443!-- Add turbulence generated by wave breaking (in ocean mode only)
     1444    IF ( wave_breaking  .AND.                                                  &
     1445         intermediate_timestep_count == intermediate_timestep_count_max )      &
     1446    THEN
     1447       CALL wave_breaking_term( 1 )
     1448    ENDIF
     1449
     1450!
    14071451!-- Calculate tendencies for the next Runge-Kutta step
    14081452    IF ( timestep_scheme(1:5) == 'runge' )  THEN
     
    14661510!-- Nudging
    14671511    IF ( nudging )  CALL nudge( simulated_time, 'v' )
     1512
     1513!
     1514!-- Effect of Stokes drift (in ocean mode only)
     1515    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
    14681516
    14691517!
     
    14901538
    14911539!
     1540!-- Add turbulence generated by wave breaking (in ocean mode only)
     1541    IF ( wave_breaking  .AND.                                                  &
     1542         intermediate_timestep_count == intermediate_timestep_count_max )      &
     1543    THEN
     1544       CALL wave_breaking_term( 2 )
     1545    ENDIF
     1546
     1547!
    14921548!-- Calculate tendencies for the next Runge-Kutta step
    14931549    IF ( timestep_scheme(1:5) == 'runge' )  THEN
     
    15471603!-- Drag by plant canopy
    15481604    IF ( plant_canopy )  CALL pcm_tendency( 3 )
     1605
     1606!
     1607!-- Effect of Stokes drift (in ocean mode only)
     1608    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
    15491609
    15501610!
  • palm/trunk/SOURCE/subsidence_mod.f90

    r3294 r3302  
    2525! -----------------
    2626! $Id$
     27! small message change
     28!
     29! 3294 2018-10-01 02:37:10Z raasch
    2730! ocean renamed ocean_mode
    2831!
     
    151154
    152155       IF ( ocean_mode )  THEN
    153           message_string = 'Applying large scale vertical motion is not ' //   &
     156          message_string = 'applying large scale vertical motion is not ' //   &
    154157                           'allowed for ocean mode'
    155158          CALL message( 'init_w_subsidence', 'PA0324', 2, 2, 0, 6, 0 )
Note: See TracChangeset for help on using the changeset viewer.