Changeset 3302 for palm/trunk
- Timestamp:
- Oct 3, 2018 2:39:40 AM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 9 edited
-
advec_s_pw.f90 (modified) (7 diffs)
-
advec_ws.f90 (modified) (20 diffs)
-
check_parameters.f90 (modified) (2 diffs)
-
header.f90 (modified) (3 diffs)
-
init_3d_model.f90 (modified) (2 diffs)
-
modules.f90 (modified) (2 diffs)
-
ocean_mod.f90 (modified) (21 diffs)
-
prognostic_equations.f90 (modified) (12 diffs)
-
subsidence_mod.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_s_pw.f90
r2718 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Stokes drift velocity added 28 ! 29 ! 2718 2018-01-02 08:49:38Z maronga 27 30 ! Corrected "Former revisions" section 28 31 ! … … 104 107 105 108 USE arrays_3d, & 106 ONLY: dd2zu, tend, u, v, w109 ONLY: dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w 107 110 108 111 USE control_parameters, & … … 124 127 INTEGER(iwp) :: j !< 125 128 INTEGER(iwp) :: k !< 129 130 REAL(wp) :: gu !< local additional advective velocity 131 REAL(wp) :: gv !< local additional advective velocity 126 132 127 133 #if defined( __nopointer ) … … 135 141 DO j = nys, nyn 136 142 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 ) ) 148 167 ENDDO 149 168 ENDDO … … 161 180 162 181 USE arrays_3d, & 163 ONLY: dd2zu, tend, u, v, w182 ONLY: dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w 164 183 165 184 USE control_parameters, & … … 180 199 INTEGER(iwp) :: j !< 181 200 INTEGER(iwp) :: k !< 201 202 REAL(wp) :: gu !< local additional advective velocity 203 REAL(wp) :: gv !< local additional advective velocity 182 204 183 205 #if defined( __nopointer ) … … 189 211 190 212 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 ) ) 202 237 ENDDO 203 238 -
palm/trunk/SOURCE/advec_ws.f90
r3298 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Stokes drift velocity added in scalar advection 28 ! 29 ! 3298 2018-10-02 12:21:11Z kanani 27 30 ! Add formatted note from ketelsen about statistics for chemistry 28 31 ! … … 1147 1150 1148 1151 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 1150 1154 1151 1155 USE control_parameters, & … … 1227 1231 ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1) 1228 1232 1229 v_comp = v(k,j,i) - v_gtrans 1233 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) 1230 1234 swap_flux_y_local(k,tn) = v_comp * ( & 1231 1235 ( 37.0_wp * ibit5 * adv_sca_5 & … … 1263 1267 DO k = nzb_max+1, nzt 1264 1268 1265 v_comp = v(k,j,i) - v_gtrans 1269 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) 1266 1270 swap_flux_y_local(k,tn) = v_comp * ( & 1267 1271 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & … … 1288 1292 ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1) 1289 1293 1290 u_comp = u(k,j,i) - u_gtrans 1294 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) 1291 1295 swap_flux_x_local(k,j,tn) = u_comp * ( & 1292 1296 ( 37.0_wp * ibit2 * adv_sca_5 & … … 1323 1327 DO k = nzb_max+1, nzt 1324 1328 1325 u_comp = u(k,j,i) - u_gtrans 1329 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) 1326 1330 swap_flux_x_local(k,j,tn) = u_comp * ( & 1327 1331 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & … … 1355 1359 ibit0 = IBITS(advc_flags_1(k,j,i),0,1) 1356 1360 1357 u_comp = u(k,j,i+1) - u_gtrans 1361 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) 1358 1362 flux_r(k) = u_comp * ( & 1359 1363 ( 37.0_wp * ibit2 * adv_sca_5 & … … 1390 1394 ibit3 = IBITS(advc_flags_1(k,j,i),3,1) 1391 1395 1392 v_comp = v(k,j+1,i) - v_gtrans 1396 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) 1393 1397 flux_n(k) = v_comp * ( & 1394 1398 ( 37.0_wp * ibit5 * adv_sca_5 & … … 1514 1518 DO k = nzb_max+1, nzt 1515 1519 1516 u_comp = u(k,j,i+1) - u_gtrans 1520 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) 1517 1521 flux_r(k) = u_comp * ( & 1518 1522 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & … … 1524 1528 + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 1525 1529 1526 v_comp = v(k,j+1,i) - v_gtrans 1530 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) 1527 1531 flux_n(k) = v_comp * ( & 1528 1532 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & … … 3192 3196 SUBROUTINE advec_s_ws( sk, sk_char ) 3193 3197 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, & 3198 3203 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 3199 3204 3200 USE grid_variables, &3205 USE grid_variables, & 3201 3206 ONLY: ddx, ddy 3202 3207 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, & 3205 3210 nzt, advc_flags_1 3206 3211 3207 3212 USE kinds 3208 3213 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, & 3212 3217 sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep 3213 3218 … … 3272 3277 ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1) 3273 3278 3274 u_comp = u(k,j,i) - u_gtrans 3279 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) 3275 3280 swap_flux_x_local(k,j) = u_comp * ( & 3276 3281 ( 37.0_wp * ibit2 * adv_sca_5 & … … 3307 3312 DO k = nzb_max+1, nzt 3308 3313 3309 u_comp = u(k,j,i) - u_gtrans 3314 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) 3310 3315 swap_flux_x_local(k,j) = u_comp * ( & 3311 3316 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & … … 3333 3338 ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1) 3334 3339 3335 v_comp = v(k,j,i) - v_gtrans 3340 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) 3336 3341 swap_flux_y_local(k) = v_comp * ( & 3337 3342 ( 37.0_wp * ibit5 * adv_sca_5 & … … 3369 3374 DO k = nzb_max+1, nzt 3370 3375 3371 v_comp = v(k,j,i) - v_gtrans 3376 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) 3372 3377 swap_flux_y_local(k) = v_comp * ( & 3373 3378 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & … … 3375 3380 + ( sk(k,j+2,i) + sk(k,j-3,i) ) & 3376 3381 ) * adv_sca_5 3377 swap_diss_y_local(k) = -ABS( v_comp ) * (&3382 swap_diss_y_local(k) = -ABS( v_comp ) * ( & 3378 3383 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & 3379 3384 - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & … … 3396 3401 ibit0 = IBITS(advc_flags_1(k,j,i),0,1) 3397 3402 3398 u_comp = u(k,j,i+1) - u_gtrans 3403 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) 3399 3404 flux_r(k) = u_comp * ( & 3400 3405 ( 37.0_wp * ibit2 * adv_sca_5 & … … 3431 3436 ibit3 = IBITS(advc_flags_1(k,j,i),3,1) 3432 3437 3433 v_comp = v(k,j+1,i) - v_gtrans 3438 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) 3434 3439 flux_n(k) = v_comp * ( & 3435 3440 ( 37.0_wp * ibit5 * adv_sca_5 & … … 3548 3553 DO k = nzb_max+1, nzt 3549 3554 3550 u_comp = u(k,j,i+1) - u_gtrans 3555 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) 3551 3556 flux_r(k) = u_comp * ( & 3552 3557 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & … … 3558 3563 + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 3559 3564 3560 v_comp = v(k,j+1,i) - v_gtrans 3565 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) 3561 3566 flux_n(k) = v_comp * ( & 3562 3567 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & -
palm/trunk/SOURCE/check_parameters.f90
r3301 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! small code rearrangements 28 ! 29 ! 3301 2018-10-02 17:10:50Z gronemeier 27 30 ! corrected former revision section 28 31 ! … … 1468 1471 1469 1472 !-- 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 1471 1483 CALL tcm_check_parameters 1472 CALL stg_check_parameters1473 IF ( plant_canopy ) CALL pcm_check_parameters1474 IF ( calculate_spectra ) CALL spectra_check_parameters1475 IF ( land_surface ) CALL lsm_check_parameters1476 IF ( bulk_cloud_model ) CALL bcm_check_parameters1477 1484 IF ( urban_surface ) CALL usm_check_parameters 1478 1485 IF ( wind_turbine ) CALL wtm_check_parameters 1479 IF ( radiation ) CALL radiation_check_parameters1480 IF ( gust_module_enabled ) CALL gust_check_parameters1481 IF ( large_scale_forcing .OR. nudging ) CALL lsf_nudging_check_parameters1482 1486 1483 1487 ! -
palm/trunk/SOURCE/header.f90
r3298 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! call of ocean_header 28 ! 29 ! 3298 2018-10-02 12:21:11Z kanani 27 30 ! Minor formatting (kanani) 28 31 ! Add chemistry header (basit) … … 437 440 438 441 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 441 445 442 446 USE particle_attributes, & … … 1120 1124 ENDIF 1121 1125 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 ) 1131 1134 1132 1135 IF ( air_chemistry ) CALL chem_header ( io ) -
palm/trunk/SOURCE/init_3d_model.f90
r3298 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! allocate and set stokes drift velocity profiles 28 ! 29 ! 3298 2018-10-02 12:21:11Z kanani 27 30 ! Minor formatting (kanani) 28 31 ! Added CALL to the chem_emissions module (Russo) … … 788 791 #endif 789 792 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 790 803 791 804 ! -
palm/trunk/SOURCE/modules.f90
r3298 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! +u/v_stokes_zu, u/v_stokes_zw 28 ! 29 ! 3298 2018-10-02 12:21:11Z kanani 27 30 ! Minor formatting/clean-up (kanani) 28 31 ! Added some variables for time treatment (Russo) … … 739 742 REAL(wp), DIMENSION(:), ALLOCATABLE :: ug !< geostrophic wind component in x-direction 740 743 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 741 746 REAL(wp), DIMENSION(:), ALLOCATABLE :: vg !< geostrophic wind component in y-direction 742 747 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 743 750 REAL(wp), DIMENSION(:), ALLOCATABLE :: w_subs !< subsidence/ascent velocity 744 751 REAL(wp), DIMENSION(:), ALLOCATABLE :: zu !< vertical grid coordinate of u-grid (in m) -
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 2018-10-01 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 (' --> Craik-Leibovich 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 !-- 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 1794 2237 END MODULE ocean_mod -
palm/trunk/SOURCE/prognostic_equations.f90
r3298 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Stokes drift + wave breaking term added 28 ! 29 ! 3298 2018-10-02 12:21:11Z kanani 27 30 ! Code added for decycling chemistry (basit) 28 31 ! … … 392 395 393 396 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 395 399 396 400 USE plant_canopy_model_mod, & … … 595 599 596 600 ! 601 !-- Effect of Stokes drift (in ocean mode only) 602 IF ( stokes_force ) CALL stokes_drift_terms( i, j, 1 ) 603 604 ! 597 605 !-- Forces by wind turbines 598 606 IF ( wind_turbine ) CALL wtm_tendencies( i, j, 1 ) … … 612 620 ) 613 621 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 614 630 615 631 ! … … 666 682 667 683 ! 684 !-- Effect of Stokes drift (in ocean mode only) 685 IF ( stokes_force ) CALL stokes_drift_terms( i, j, 2 ) 686 687 ! 668 688 !-- Forces by wind turbines 669 689 IF ( wind_turbine ) CALL wtm_tendencies( i, j, 2 ) … … 682 702 ) 683 703 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 684 712 685 713 ! … … 731 759 !-- Drag by plant canopy 732 760 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 ) 733 765 734 766 ! … … 1383 1415 1384 1416 ! 1417 !-- Effect of Stokes drift (in ocean mode only) 1418 IF ( stokes_force ) CALL stokes_drift_terms( 1 ) 1419 1420 ! 1385 1421 !-- Forces by wind turbines 1386 1422 IF ( wind_turbine ) CALL wtm_tendencies( 1 ) … … 1405 1441 1406 1442 ! 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 ! 1407 1451 !-- Calculate tendencies for the next Runge-Kutta step 1408 1452 IF ( timestep_scheme(1:5) == 'runge' ) THEN … … 1466 1510 !-- Nudging 1467 1511 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 ) 1468 1516 1469 1517 ! … … 1490 1538 1491 1539 ! 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 ! 1492 1548 !-- Calculate tendencies for the next Runge-Kutta step 1493 1549 IF ( timestep_scheme(1:5) == 'runge' ) THEN … … 1547 1603 !-- Drag by plant canopy 1548 1604 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 ) 1549 1609 1550 1610 ! -
palm/trunk/SOURCE/subsidence_mod.f90
r3294 r3302 25 25 ! ----------------- 26 26 ! $Id$ 27 ! small message change 28 ! 29 ! 3294 2018-10-01 02:37:10Z raasch 27 30 ! ocean renamed ocean_mode 28 31 ! … … 151 154 152 155 IF ( ocean_mode ) THEN 153 message_string = ' Applying large scale vertical motion is not ' // &156 message_string = 'applying large scale vertical motion is not ' // & 154 157 'allowed for ocean mode' 155 158 CALL message( 'init_w_subsidence', 'PA0324', 2, 2, 0, 6, 0 )
Note: See TracChangeset
for help on using the changeset viewer.
