- Timestamp:
- Mar 20, 2020 4:14:41 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
Makefile (modified) (5 diffs)
-
advec_ws.f90 (modified) (198 diffs)
-
time_integration.f90 (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r4458 r4466 25 25 # ----------------- 26 26 # $Id$ 27 # cpu measures in advec_ws added 28 # 29 # 4458 2020-03-11 15:37:31Z raasch 27 30 # bugfix for r4457: missing dependency added 28 # 31 # 29 32 # 4457 2020-03-11 14:20:43Z raasch 30 33 # exchange horiz has been modularized and exchange horiz 2d has been merged, dependencies updated 31 34 # accordingly 32 # 35 # 33 36 # 4453 2020-03-11 08:10:13Z raasch 34 37 # dependencies for exchange horiz modified 35 # 38 # 36 39 # 4434 2020-03-03 10:02:18Z oliver.maas 37 40 # Added output control for wind turbine model 38 # 41 # 39 42 # 4414 2020-02-19 20:16:04Z suehring 40 43 # Move dependencies for init grid from advection scheme and multigrid solver 41 44 # to module_interface 42 # 45 # 43 46 # 4411 2020-02-18 14:28:02Z maronga 44 47 # Added output routines for WTM 45 # 48 # 46 49 # 4400 2020-02-10 20:32:41Z suehring 47 50 # Add dependency for data-output module 48 # 51 # 49 52 # 4392 2020-01-31 16:14:57Z pavelkrc 50 53 # add dependency on fft for transpose 51 # 54 # 52 55 # 4347 2019-12-18 13:18:33Z suehring 53 56 # add dependency to basic_constants_and_equations_mod for dynamics_mod 54 # 57 # 55 58 # 4331 2019-12-10 18:25:02Z suehring 56 59 # Move diagnostic surface output to diagnostic_output_quantities 57 # 60 # 58 61 # 4309 2019-11-26 18:49:59Z suehring 59 62 # Add dependency to parallel random generator for synthetic turbulence generator 60 # 63 # 61 64 # 4286 2019-10-30 16:01:14Z resler 62 65 # delete boundary_conds, added missing dependencies 63 # 66 # 64 67 # 4270 2019-10-23 10:46:20Z monakurppa 65 68 # - Implement offline nesting for salsa and add dependency of nesting_offl_mod … … 68 71 # for salsa 69 72 # - Add dependency on basic_constants_and_equations_mod for salsa_mod 70 # 73 # 71 74 # 4258 2019-10-07 13:29:08Z suehring 72 75 # Add dependency of land-surface model on pmc_handle_communicator and cpu_log 73 # 76 # 74 77 # 4245 2019-09-30 08:40:37Z pavelkrc 75 78 # Remove no longer needed dependencies on surface_mod 76 # 79 # 77 80 # 4227 2019-09-10 18:04:34Z gronemeier 78 81 # Add palm_date_time_mod, remove date_and_time_mod … … 80 83 # 4223 2019-09-10 09:20:47Z gronemeier 81 84 # Corrected "Former revisions" section 82 # 85 # 83 86 # 4174 2019-08-20 12:41:13Z gronemeier 84 87 # bugfix: add missing dependencies for vdi_internal_controls 85 # 88 # 86 89 # 4173 2019-08-20 12:04:06Z gronemeier 87 90 # add vdi_internal_controls 88 # 91 # 89 92 # 4168 2019-08-16 13:50:17Z suehring 90 93 # Remove some dependencies on surface_mod that are no longer required without 91 94 # get_topography_top_index functions 92 # 95 # 93 96 # 4167 2019-08-16 11:01:48Z suehring 94 97 # Remove no longer needed dependencies on surface_mod 95 # 98 # 96 99 # 4127 2019-07-30 14:47:10Z suehring 97 # Add dependency of data_output_3d on plant_canopy_model_mod 100 # Add dependency of data_output_3d on plant_canopy_model_mod 98 101 # (merge from branch resler) 99 102 # 100 103 # 4106 2019-07-19 08:54:42Z gronemeier 101 104 # Remove dependency on pmc_interface for boundary_conds 102 # 105 # 103 106 # 4070 2019-07-03 13:51:40Z gronemeier 104 107 # Add new data output modules … … 358 361 modules.o 359 362 advec_ws.o: \ 363 cpulog_mod.o \ 360 364 exchange_horiz_mod.o \ 361 365 mod_kinds.o \ … … 825 829 modules.o \ 826 830 netcdf_data_input_mod.o \ 827 salsa_mod.o 831 salsa_mod.o 828 832 netcdf_data_input_mod.o: \ 829 833 chem_modules.o \ -
palm/trunk/SOURCE/advec_ws.f90
r4457 r4466 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - vector branch further optimized (linear dependencies along z removed and 28 ! loops are splitted) 29 ! - topography closed channel flow with symmetric boundaries also implemented 30 ! in vector branch 31 ! - some formatting adjustments made and comments added 32 ! - cpu measures for vector branch added 33 ! 34 ! 4457 2020-03-11 14:20:43Z raasch 27 35 ! use statement for exchange horiz added 28 ! 36 ! 29 37 ! 4414 2020-02-19 20:16:04Z suehring 30 38 ! Move call for initialization of control flags to ws_init 31 ! 39 ! 32 40 ! 4360 2020-01-07 11:25:50Z suehring 33 41 ! Introduction of wall_flags_total_0, which currently sets bits based on static 34 42 ! topography information used in wall_flags_static_0 35 ! 43 ! 36 44 ! 4340 2019-12-16 08:17:03Z Giersch 37 ! Topography closed channel flow with symmetric boundaries implemented 38 ! 45 ! Topography closed channel flow with symmetric boundaries implemented 46 ! 39 47 ! 4330 2019-12-10 16:16:33Z knoop 40 48 ! Bugix: removed syntax error introduced by last commit 41 ! 49 ! 42 50 ! 4329 2019-12-10 15:46:36Z motisi 43 51 ! Renamed wall_flags_0 to wall_flags_static_0 44 ! 52 ! 45 53 ! 4328 2019-12-09 18:53:04Z suehring 46 54 ! Minor formatting adjustments 47 ! 55 ! 48 56 ! 4327 2019-12-06 14:48:31Z Giersch 49 57 ! Setting of advection flags for vertical fluxes of w revised, air density for 50 58 ! vertical flux calculation of w at k=1 is considered now 51 ! 59 ! 52 60 ! 4325 2019-12-06 07:14:04Z Giersch 53 ! Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of 61 ! Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of 54 62 ! advection flags for fluxes in z-direction revised, comments extended 55 ! 63 ! 56 64 ! 4324 2019-12-06 07:11:33Z Giersch 57 65 ! Indirect indexing for calculating vertical fluxes close to boundaries is only 58 66 ! used for loop indizes where it is really necessary 59 ! 67 ! 60 68 ! 4317 2019-12-03 12:43:22Z Giersch 61 ! Comments revised/added, formatting improved, fluxes for u,v, and scalars are 69 ! Comments revised/added, formatting improved, fluxes for u,v, and scalars are 62 70 ! explicitly set to zero at nzt+1, fluxes of w-component are now calculated only 63 71 ! until nzt-1 (Prognostic equation for w-velocity component ends at nzt-1) 64 ! 72 ! 65 73 ! 4204 2019-08-30 12:30:17Z knoop 66 74 ! Bugfix: Changed sk_num initialization default to avoid implicit SAVE-Attribut 67 ! 75 ! 68 76 ! 4182 2019-08-22 15:20:23Z scharf 69 77 ! Corrected "Former revisions" section 70 ! 78 ! 71 79 ! 4110 2019-07-22 17:05:21Z suehring 72 80 ! - Separate initialization of advection flags for momentum and scalars. In this 73 ! context, resort the bits and do some minor formatting. 74 ! - Make flag initialization for scalars more flexible, introduce an 75 ! arguemnt list to indicate non-cyclic boundaries (required for decycled 81 ! context, resort the bits and do some minor formatting. 82 ! - Make flag initialization for scalars more flexible, introduce an 83 ! arguemnt list to indicate non-cyclic boundaries (required for decycled 76 84 ! scalars such as chemical species or aerosols) 77 ! - Introduce extended 'degradation zones', where horizontal advection of 78 ! passive scalars is discretized by first-order scheme at all grid points 79 ! that in the vicinity of buildings (<= 3 grid points). Even though no 80 ! building is within the numerical stencil, first-order scheme is used. 85 ! - Introduce extended 'degradation zones', where horizontal advection of 86 ! passive scalars is discretized by first-order scheme at all grid points 87 ! that in the vicinity of buildings (<= 3 grid points). Even though no 88 ! building is within the numerical stencil, first-order scheme is used. 81 89 ! At fourth and fifth grid point the order of the horizontal advection scheme 82 90 ! is successively upgraded. 83 ! These extended degradation zones are used to avoid stationary numerical 91 ! These extended degradation zones are used to avoid stationary numerical 84 92 ! oscillations, which are responsible for high concentration maxima that may 85 ! appear under shear-free stable conditions. 86 ! - Change interface for scalar advection routine. 93 ! appear under shear-free stable conditions. 94 ! - Change interface for scalar advection routine. 87 95 ! - Bugfix, avoid uninitialized value sk_num in vector version of scalar 88 96 ! advection 89 ! 97 ! 90 98 ! 4109 2019-07-22 17:00:34Z suehring 91 ! Implementation of a flux limiter according to Skamarock (2006) for the 92 ! vertical scalar advection. Please note, this is only implemented for the 99 ! Implementation of a flux limiter according to Skamarock (2006) for the 100 ! vertical scalar advection. Please note, this is only implemented for the 93 101 ! cache-optimized version at the moment. Implementation for the vector- 94 ! optimized version will follow after critical issues concerning 102 ! optimized version will follow after critical issues concerning 95 103 ! vectorization are fixed. 96 ! 104 ! 97 105 ! 3873 2019-04-08 15:44:30Z knoop 98 106 ! Moved ocean_mode specific code to ocean_mod 99 ! 107 ! 100 108 ! 3872 2019-04-08 15:03:06Z knoop 101 109 ! Moved all USE statements to module level + removed salsa dependency 102 ! 110 ! 103 111 ! 3871 2019-04-08 14:38:39Z knoop 104 112 ! Moving initialization of bcm specific flux arrays into bulk_cloud_model_mod 105 ! 113 ! 106 114 ! 3864 2019-04-05 09:01:56Z monakurppa 107 115 ! Remove tailing white spaces 108 ! 116 ! 109 117 ! 3696 2019-01-24 16:37:35Z suehring 110 118 ! Bugfix in degradation height 111 ! 119 ! 112 120 ! 3661 2019-01-08 18:22:50Z suehring 113 ! - Minor bugfix in divergence correction (only has implications at 121 ! - Minor bugfix in divergence correction (only has implications at 114 122 ! downward-facing wall surfaces) 115 123 ! - Remove setting of Neumann condition for horizontal velocity variances 116 124 ! - Split loops for tendency calculation and divergence correction in order to 117 125 ! reduce bit queries 118 ! - Introduce new parameter nzb_max_l to better control order degradation at 126 ! - Introduce new parameter nzb_max_l to better control order degradation at 119 127 ! non-cyclic boundaries 120 ! 128 ! 121 129 ! 3655 2019-01-07 16:51:22Z knoop 122 130 ! OpenACC port for SPEC … … 133 141 ! ------------ 134 142 !> Advection scheme for scalars and momentum using the flux formulation of 135 !> Wicker and Skamarock 5th order. Additionally the module contains of a 136 !> routine using for initialisation and steering of the statical evaluation. 143 !> Wicker and Skamarock 5th order. Additionally the module contains of a 144 !> routine using for initialisation and steering of the statical evaluation. 137 145 !> The computation of turbulent fluxes takes place inside the advection 138 146 !> routines. … … 140 148 !> degraded. 141 149 !> A divergence correction is applied. It is necessary for topography, since 142 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes 143 !> and could lead to numerical instabilities. 150 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes 151 !> and could lead to numerical instabilities. 144 152 !> 145 153 !> @todo Implement monotonic flux limiter also for vector version. … … 154 162 u_stokes_zu, v_stokes_zu, & 155 163 diss_l_diss, diss_l_e, diss_l_pt, diss_l_q, & 156 diss_l_s, diss_l_ sa, diss_l_u, diss_l_v, diss_l_w,&164 diss_l_s, diss_l_u, diss_l_v, diss_l_w, & 157 165 flux_l_diss, flux_l_e, flux_l_pt, flux_l_q, flux_l_s, & 158 flux_l_ sa, flux_l_u, flux_l_v, flux_l_w,&166 flux_l_u, flux_l_v, flux_l_w, & 159 167 diss_s_diss, diss_s_e, diss_s_pt, diss_s_q, diss_s_s, & 160 diss_s_ sa, diss_s_u, diss_s_v, diss_s_w,&168 diss_s_u, diss_s_v, diss_s_w, & 161 169 flux_s_diss, flux_s_e, flux_s_pt, flux_s_q, flux_s_s, & 162 flux_s_ sa, flux_s_u, flux_s_v, flux_s_w170 flux_s_u, flux_s_v, flux_s_w 163 171 164 172 USE control_parameters, & 165 ONLY: air_chemistry, & 166 bc_dirichlet_l, & 173 ONLY: bc_dirichlet_l, & 167 174 bc_dirichlet_n, & 168 175 bc_dirichlet_r, & … … 176 183 passive_scalar, & 177 184 rans_tke_e, & 178 momentum_advec, &179 salsa, &180 scalar_advec, &181 185 symmetry_flag, & 182 186 intermediate_timestep_count, & … … 186 190 ws_scheme_sca, & 187 191 dt_3d 192 193 USE cpulog, & 194 ONLY: cpu_log, & 195 log_point_s 188 196 189 197 USE exchange_horiz_mod, & … … 198 206 nxlg, & 199 207 nxlu, & 200 nxr, & 201 nxrg, & 208 nxr, & 209 nxrg, & 202 210 ny, & 203 nyn, & 204 nyng, & 211 nyn, & 212 nyng, & 205 213 nys, & 206 214 nysg, & … … 242 250 INTERFACE ws_init 243 251 MODULE PROCEDURE ws_init 244 END INTERFACE ws_init 245 252 END INTERFACE ws_init 253 246 254 INTERFACE ws_init_flags_momentum 247 255 MODULE PROCEDURE ws_init_flags_momentum 248 256 END INTERFACE ws_init_flags_momentum 249 257 250 258 INTERFACE ws_init_flags_scalar 251 259 MODULE PROCEDURE ws_init_flags_scalar … … 287 295 288 296 ! 289 !-- Set the appropriatefactors for scalar and momentum advection.297 !-- Set factors for scalar and momentum advection. 290 298 adv_sca_5 = 1.0_wp / 60.0_wp 291 299 adv_sca_3 = 1.0_wp / 12.0_wp … … 294 302 adv_mom_3 = 1.0_wp / 24.0_wp 295 303 adv_mom_1 = 1.0_wp / 4.0_wp 296 ! 304 ! 297 305 !-- Arrays needed for statical evaluation of fluxes. 298 306 IF ( ws_scheme_mom ) THEN … … 330 338 331 339 ! 332 !-- Arrays needed for reasons of speed optimization for cache version. 333 !-- For the vector version the buffer arrays are not necessary, 334 !-- because the the fluxes can swapped directly inside the loops of the 335 !-- advection routines. 340 !-- Arrays needed for reasons of speed optimization 341 IF ( ws_scheme_mom ) THEN 342 343 ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), & 344 flux_s_v(nzb+1:nzt,0:threads_per_task-1), & 345 flux_s_w(nzb+1:nzt,0:threads_per_task-1), & 346 diss_s_u(nzb+1:nzt,0:threads_per_task-1), & 347 diss_s_v(nzb+1:nzt,0:threads_per_task-1), & 348 diss_s_w(nzb+1:nzt,0:threads_per_task-1) ) 349 ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 350 flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 351 flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 352 diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 353 diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 354 diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 355 356 ENDIF 357 ! 358 !-- For the vector version the buffer arrays for scalars are not necessary, 359 !-- since internal arrays are used in the vector version. 336 360 IF ( loop_optimization /= 'vector' ) THEN 337 338 IF ( ws_scheme_mom ) THEN339 340 ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &341 flux_s_v(nzb+1:nzt,0:threads_per_task-1), &342 flux_s_w(nzb+1:nzt,0:threads_per_task-1), &343 diss_s_u(nzb+1:nzt,0:threads_per_task-1), &344 diss_s_v(nzb+1:nzt,0:threads_per_task-1), &345 diss_s_w(nzb+1:nzt,0:threads_per_task-1) )346 ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &347 flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &348 flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &349 diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &350 diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &351 diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )352 353 ENDIF354 355 361 IF ( ws_scheme_sca ) THEN 356 362 … … 358 364 flux_s_e(nzb+1:nzt,0:threads_per_task-1), & 359 365 diss_s_pt(nzb+1:nzt,0:threads_per_task-1), & 360 diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) 366 diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) 361 367 ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 362 368 flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & … … 386 392 387 393 ENDIF 388 389 394 ENDIF 390 395 ! 391 !-- Initialize the flag arrays controlling degradation near walls, i.e. 396 !-- Initialize the flag arrays controlling degradation near walls, i.e. 392 397 !-- to decrease the numerical stencil appropriately. The order of the scheme 393 398 !-- is degraded near solid walls as well as near non-cyclic inflow and outflow … … 412 417 !> Initialization of flags to control the order of the advection scheme near 413 418 !> solid walls and non-cyclic inflow boundaries, where the order is sucessively 414 !> degraded. 419 !> degraded. 415 420 !------------------------------------------------------------------------------! 416 421 SUBROUTINE ws_init_flags_momentum … … 437 442 DO j = nys, nyn 438 443 DO k = nzb+1, nzt 439 ! 440 !-- At first, set flags to WS1. 441 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 444 ! 445 !-- At first, set flags to WS1. 446 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 442 447 !-- in order to handle the left/south flux. 443 !-- near vertical walls. 448 !-- near vertical walls. 444 449 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 ) 445 450 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 ) … … 452 457 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 453 458 .AND. i == nxr ) ) & 454 THEN 455 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 ) 459 THEN 460 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 ) 456 461 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),1) .AND. & 457 462 BTEST(wall_flags_total_0(k,j,i+1),1) .OR. & … … 462 467 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 463 468 .AND. i == nxlu+1) ) & 464 THEN 465 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 ) 466 ! 467 !-- Clear flag for WS1 468 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) 469 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),1) .AND. &470 BTEST(wall_flags_total_0(k,j,i+2),1) .AND. &471 BTEST(wall_flags_total_0(k,j,i-1),1) ) &472 THEN 473 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 ) 474 ! 475 !-- Clear flag for WS1 476 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) 477 ENDIF 478 ! 479 !-- u component - y-direction 480 !-- WS1 (3), WS3 (4), WS5 (5) 469 THEN 470 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 ) 471 ! 472 !-- Clear flag for WS1 473 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) 474 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),1) .AND. & 475 BTEST(wall_flags_total_0(k,j,i+2),1) .AND. & 476 BTEST(wall_flags_total_0(k,j,i-1),1) ) & 477 THEN 478 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 ) 479 ! 480 !-- Clear flag for WS1 481 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) 482 ENDIF 483 ! 484 !-- u component - y-direction 485 !-- WS1 (3), WS3 (4), WS5 (5) 481 486 IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),1) .OR. & 482 487 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & … … 484 489 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 485 490 .AND. j == nyn ) ) & 486 THEN 487 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 ) 491 THEN 492 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 ) 488 493 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),1) .AND. & 489 494 BTEST(wall_flags_total_0(k,j+1,i),1) .OR. & … … 494 499 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 495 500 .AND. j == nyn-1 ) ) & 496 THEN 497 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 ) 498 ! 499 !-- Clear flag for WS1 500 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) 501 THEN 502 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 ) 503 ! 504 !-- Clear flag for WS1 505 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) 501 506 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),1) .AND. & 502 507 BTEST(wall_flags_total_0(k,j+2,i),1) .AND. & 503 508 BTEST(wall_flags_total_0(k,j-1,i),1) ) & 504 THEN 505 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 ) 506 ! 507 !-- Clear flag for WS1 508 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) 509 ENDIF 510 ! 511 !-- u component - z-direction. Fluxes are calculated on w-grid 512 !-- level. Boundary u-values at/within walls aren't used. 513 !-- WS1 (6), WS3 (7), WS5 (8) 514 IF ( k == nzb+1 ) THEN 515 k_mm = nzb 516 ELSE 517 k_mm = k - 2 518 ENDIF 519 IF ( k > nzt-1 ) THEN 520 k_pp = nzt+1 521 ELSE 522 k_pp = k + 2 523 ENDIF 524 IF ( k > nzt-2 ) THEN 525 k_ppp = nzt+1 526 ELSE 527 k_ppp = k + 3 528 ENDIF 529 530 flag_set = .FALSE. 531 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1) .AND. 532 BTEST(wall_flags_total_0(k,j,i),1) .AND. 533 BTEST(wall_flags_total_0(k+1,j,i),1) ) .OR. 534 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. &535 BTEST(wall_flags_total_0(k+1,j,i),1) .AND. 536 BTEST(wall_flags_total_0(k,j,i),1) ) .OR. 537 ( k == nzt .AND. symmetry_flag == 0 ) ) 538 THEN 539 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 ) 540 flag_set = .TRUE. 541 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1) .OR. 542 .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND. &543 BTEST(wall_flags_total_0(k-1,j,i),1) .AND. 544 BTEST(wall_flags_total_0(k,j,i),1) .AND. 545 BTEST(wall_flags_total_0(k+1,j,i),1) .AND. 546 BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. 547 .NOT. flag_set .OR. &548 ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) 549 THEN 550 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 ) 551 flag_set = .TRUE. 552 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1) .AND. 553 BTEST(wall_flags_total_0(k-1,j,i),1) .AND. 554 BTEST(wall_flags_total_0(k,j,i),1) .AND. 555 BTEST(wall_flags_total_0(k+1,j,i),1) .AND. 556 BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. 557 BTEST(wall_flags_total_0(k_ppp,j,i),1) .AND. 558 .NOT. flag_set ) 509 THEN 510 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 ) 511 ! 512 !-- Clear flag for WS1 513 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) 514 ENDIF 515 ! 516 !-- u component - z-direction. Fluxes are calculated on w-grid 517 !-- level. Boundary u-values at/within walls aren't used. 518 !-- WS1 (6), WS3 (7), WS5 (8) 519 IF ( k == nzb+1 ) THEN 520 k_mm = nzb 521 ELSE 522 k_mm = k - 2 523 ENDIF 524 IF ( k > nzt-1 ) THEN 525 k_pp = nzt+1 526 ELSE 527 k_pp = k + 2 528 ENDIF 529 IF ( k > nzt-2 ) THEN 530 k_ppp = nzt+1 531 ELSE 532 k_ppp = k + 3 533 ENDIF 534 535 flag_set = .FALSE. 536 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1) .AND. & 537 BTEST(wall_flags_total_0(k,j,i),1) .AND. & 538 BTEST(wall_flags_total_0(k+1,j,i),1) ) .OR. & 539 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. & 540 BTEST(wall_flags_total_0(k+1,j,i),1) .AND. & 541 BTEST(wall_flags_total_0(k,j,i),1) ) .OR. & 542 ( k == nzt .AND. symmetry_flag == 0 ) ) & 543 THEN 544 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 ) 545 flag_set = .TRUE. 546 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1) .OR. & 547 .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND.& 548 BTEST(wall_flags_total_0(k-1,j,i),1) .AND.& 549 BTEST(wall_flags_total_0(k,j,i),1) .AND.& 550 BTEST(wall_flags_total_0(k+1,j,i),1) .AND.& 551 BTEST(wall_flags_total_0(k_pp,j,i),1) .AND.& 552 .NOT. flag_set .OR.& 553 ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) & 554 THEN 555 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 ) 556 flag_set = .TRUE. 557 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1) .AND.& 558 BTEST(wall_flags_total_0(k-1,j,i),1) .AND.& 559 BTEST(wall_flags_total_0(k,j,i),1) .AND.& 560 BTEST(wall_flags_total_0(k+1,j,i),1) .AND.& 561 BTEST(wall_flags_total_0(k_pp,j,i),1) .AND.& 562 BTEST(wall_flags_total_0(k_ppp,j,i),1) .AND.& 563 .NOT. flag_set ) & 559 564 THEN 560 565 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 ) … … 570 575 DO k = nzb+1, nzt 571 576 ! 572 !-- At first, set flags to WS1. 573 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 577 !-- At first, set flags to WS1. 578 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 574 579 !-- in order to handle the left/south flux. 575 580 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 ) … … 583 588 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 584 589 .AND. i == nxr ) ) & 585 THEN 586 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 ) 587 ! 588 !-- WS3 590 THEN 591 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 ) 592 ! 593 !-- WS3 589 594 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),2) .AND. & 590 595 BTEST(wall_flags_total_0(k,j,i+1),2) ) .OR. & … … 595 600 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 596 601 .AND. i == nxlu ) ) & 597 THEN 598 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 ) 599 ! 600 !-- Clear flag for WS1 601 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) 602 THEN 603 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 ) 604 ! 605 !-- Clear flag for WS1 606 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) 602 607 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),2) .AND. & 603 608 BTEST(wall_flags_total_0(k,j,i+2),2) .AND. & 604 609 BTEST(wall_flags_total_0(k,j,i-1),2) ) & 605 THEN 606 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 ) 607 ! 608 !-- Clear flag for WS1 609 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) 610 ENDIF 611 ! 612 !-- v component - y-direction 613 !-- WS1 (12), WS3 (13), WS5 (14) 610 THEN 611 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 ) 612 ! 613 !-- Clear flag for WS1 614 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) 615 ENDIF 616 ! 617 !-- v component - y-direction 618 !-- WS1 (12), WS3 (13), WS5 (14) 614 619 IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),2) .OR. & 615 620 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & … … 617 622 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 618 623 .AND. j == nyn ) ) & 619 THEN 620 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 ) 624 THEN 625 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 ) 621 626 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),2) .AND. & 622 627 BTEST(wall_flags_total_0(k,j+1,i),2) .OR. & … … 627 632 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 628 633 .AND. j == nyn-1 ) ) & 629 THEN 630 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 ) 631 ! 632 !-- Clear flag for WS1 633 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) 634 THEN 635 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 ) 636 ! 637 !-- Clear flag for WS1 638 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) 634 639 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),2) .AND. & 635 640 BTEST(wall_flags_total_0(k,j+2,i),2) .AND. & 636 BTEST(wall_flags_total_0(k,j-1,i),2) ) & 637 THEN 641 BTEST(wall_flags_total_0(k,j-1,i),2) ) & 642 THEN 638 643 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 ) 639 644 ! … … 641 646 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) 642 647 ENDIF 643 ! 644 !-- v component - z-direction. Fluxes are calculated on w-grid 645 !-- level. Boundary v-values at/within walls aren't used. 648 ! 649 !-- v component - z-direction. Fluxes are calculated on w-grid 650 !-- level. Boundary v-values at/within walls aren't used. 646 651 !-- WS1 (15), WS3 (16), WS5 (17) 647 652 IF ( k == nzb+1 ) THEN 648 653 k_mm = nzb 649 ELSE 654 ELSE 650 655 k_mm = k - 2 651 656 ENDIF 652 657 IF ( k > nzt-1 ) THEN 653 658 k_pp = nzt+1 654 ELSE 659 ELSE 655 660 k_pp = k + 2 656 661 ENDIF 657 662 IF ( k > nzt-2 ) THEN 658 663 k_ppp = nzt+1 659 ELSE 664 ELSE 660 665 k_ppp = k + 3 661 ENDIF 662 666 ENDIF 667 663 668 flag_set = .FALSE. 664 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2) .AND.&665 BTEST(wall_flags_total_0(k,j,i),2) .AND.&666 BTEST(wall_flags_total_0(k+1,j,i),2) ) .OR.&667 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2) .AND. &668 BTEST(wall_flags_total_0(k+1,j,i),2) .AND.&669 BTEST(wall_flags_total_0(k,j,i),2) ) .OR.&670 ( k == nzt .AND. symmetry_flag == 0 ) ) 671 THEN 672 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 ) 673 flag_set = .TRUE. 674 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2) .OR. 675 .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),2) ) .AND. &676 BTEST(wall_flags_total_0(k-1,j,i),2) .AND. 677 BTEST(wall_flags_total_0(k,j,i),2) .AND. 678 BTEST(wall_flags_total_0(k+1,j,i),2) .AND. 679 BTEST(wall_flags_total_0(k_pp,j,i),2) .AND. 680 .NOT. flag_set .OR. &681 ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) 682 THEN 683 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 ) 684 flag_set = .TRUE. 685 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2) .AND.&686 BTEST(wall_flags_total_0(k-1,j,i),2) .AND.&687 BTEST(wall_flags_total_0(k,j,i),2) .AND.&688 BTEST(wall_flags_total_0(k+1,j,i),2) .AND.&689 BTEST(wall_flags_total_0(k_pp,j,i),2) .AND.&690 BTEST(wall_flags_total_0(k_ppp,j,i),2) .AND.&691 .NOT. flag_set ) 669 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2) .AND.& 670 BTEST(wall_flags_total_0(k,j,i),2) .AND.& 671 BTEST(wall_flags_total_0(k+1,j,i),2) ) .OR. & 672 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2) .AND.& 673 BTEST(wall_flags_total_0(k+1,j,i),2) .AND.& 674 BTEST(wall_flags_total_0(k,j,i),2) ) .OR. & 675 ( k == nzt .AND. symmetry_flag == 0 ) ) & 676 THEN 677 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 ) 678 flag_set = .TRUE. 679 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2) .OR. & 680 .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),2) ) .AND.& 681 BTEST(wall_flags_total_0(k-1,j,i),2) .AND.& 682 BTEST(wall_flags_total_0(k,j,i),2) .AND.& 683 BTEST(wall_flags_total_0(k+1,j,i),2) .AND.& 684 BTEST(wall_flags_total_0(k_pp,j,i),2) .AND.& 685 .NOT. flag_set .OR.& 686 ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) & 687 THEN 688 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 ) 689 flag_set = .TRUE. 690 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2) .AND.& 691 BTEST(wall_flags_total_0(k-1,j,i),2) .AND.& 692 BTEST(wall_flags_total_0(k,j,i),2) .AND.& 693 BTEST(wall_flags_total_0(k+1,j,i),2) .AND.& 694 BTEST(wall_flags_total_0(k_pp,j,i),2) .AND.& 695 BTEST(wall_flags_total_0(k_ppp,j,i),2) .AND.& 696 .NOT. flag_set ) & 692 697 THEN 693 698 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 ) … … 703 708 DO k = nzb+1, nzt 704 709 ! 705 !-- At first, set flags to WS1. 706 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 710 !-- At first, set flags to WS1. 711 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 707 712 !-- in order to handle the left/south flux. 708 713 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 ) … … 716 721 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 717 722 .AND. i == nxr ) ) & 718 THEN 719 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 ) 723 THEN 724 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 ) 720 725 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),3) .AND. & 721 726 BTEST(wall_flags_total_0(k,j,i+1),3) .OR. & … … 726 731 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 727 732 .AND. i == nxlu ) ) & 728 THEN 729 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 ) 730 ! 731 !-- Clear flag for WS1 732 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) 733 THEN 734 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 ) 735 ! 736 !-- Clear flag for WS1 737 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) 733 738 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),3) .AND. & 734 739 BTEST(wall_flags_total_0(k,j,i+2),3) .AND. & 735 740 BTEST(wall_flags_total_0(k,j,i-1),3) ) & 736 THEN 737 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 ) 738 ! 739 !-- Clear flag for WS1 740 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) 741 ENDIF 742 ! 743 !-- w component - y-direction 744 !-- WS1 (21), WS3 (22), WS5 (23) 741 THEN 742 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 ) 743 ! 744 !-- Clear flag for WS1 745 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) 746 ENDIF 747 ! 748 !-- w component - y-direction 749 !-- WS1 (21), WS3 (22), WS5 (23) 745 750 IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),3) .OR. & 746 751 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & … … 748 753 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 749 754 .AND. j == nyn ) ) & 750 THEN 751 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 ) 755 THEN 756 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 ) 752 757 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),3) .AND. & 753 758 BTEST(wall_flags_total_0(k,j+1,i),3) .OR. & … … 758 763 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 759 764 .AND. j == nyn-1 ) ) & 760 THEN 761 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 ) 762 ! 763 !-- Clear flag for WS1 764 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) 765 THEN 766 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 ) 767 ! 768 !-- Clear flag for WS1 769 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) 765 770 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),3) .AND. & 766 771 BTEST(wall_flags_total_0(k,j+2,i),3) .AND. & 767 772 BTEST(wall_flags_total_0(k,j-1,i),3) ) & 768 THEN 773 THEN 769 774 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 ) 770 775 ! … … 772 777 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) 773 778 ENDIF 774 ! 779 ! 775 780 !-- w component - z-direction. Fluxes are calculated on scalar grid 776 !-- level. Boundary w-values at walls are used. Flux at k=i is 777 !-- defined at scalar position k=i+1 with i being an integer. 781 !-- level. Boundary w-values at walls are used. Flux at k=i is 782 !-- defined at scalar position k=i+1 with i being an integer. 778 783 !-- WS1 (24), WS3 (25), WS5 (26) 779 784 IF ( k == nzb+1 ) THEN … … 784 789 IF ( k > nzt-1 ) THEN 785 790 k_pp = nzt+1 786 ELSE 791 ELSE 787 792 k_pp = k + 2 788 793 ENDIF 789 794 IF ( k > nzt-2 ) THEN 790 795 k_ppp = nzt+1 791 ELSE 796 ELSE 792 797 k_ppp = k + 3 793 ENDIF 794 795 flag_set = .FALSE. 796 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3) .AND. 797 BTEST(wall_flags_total_0(k+1,j,i),3) ) .OR. 798 ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3) .AND. 799 BTEST(wall_flags_total_0(k,j,i),3) ) .OR. &800 k == nzt -1 ) 798 ENDIF 799 800 flag_set = .FALSE. 801 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3) .AND. & 802 BTEST(wall_flags_total_0(k+1,j,i),3) ) .OR. & 803 ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & 804 BTEST(wall_flags_total_0(k,j,i),3) ) .OR. & 805 k == nzt -1 ) & 801 806 THEN 802 807 ! 803 808 !-- Please note, at k == nzb_w_inner(j,i) a flag is explicitly 804 !-- set, although this is not a prognostic level. However, 809 !-- set, although this is not a prognostic level. However, 805 810 !-- contrary to the advection of u,v and s this is necessary 806 811 !-- because flux_t(nzb_w_inner(j,i)) is used for the tendency … … 808 813 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 ) 809 814 flag_set = .TRUE. 810 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3) .AND. 811 BTEST(wall_flags_total_0(k,j,i),3) .AND. 812 BTEST(wall_flags_total_0(k+1,j,i),3) .AND. 813 BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR. 814 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. 815 BTEST(wall_flags_total_0(k+1,j,i),3) .AND. 816 BTEST(wall_flags_total_0(k,j,i),3) .AND. 817 BTEST(wall_flags_total_0(k-1,j,i),3) ) .AND. 818 .NOT. flag_set .OR. &819 k == nzt - 2 ) 820 THEN 821 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 ) 822 flag_set = .TRUE. 823 ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3) .AND. 824 BTEST(wall_flags_total_0(k,j,i),3) .AND. 825 BTEST(wall_flags_total_0(k+1,j,i),3) .AND. 826 BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. 827 .NOT. flag_set ) 828 THEN 815 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3) .AND. & 816 BTEST(wall_flags_total_0(k,j,i),3) .AND. & 817 BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & 818 BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR. & 819 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. & 820 BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & 821 BTEST(wall_flags_total_0(k,j,i),3) .AND. & 822 BTEST(wall_flags_total_0(k-1,j,i),3) ) .AND. & 823 .NOT. flag_set .OR. & 824 k == nzt - 2 ) & 825 THEN 826 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 ) 827 flag_set = .TRUE. 828 ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3) .AND. & 829 BTEST(wall_flags_total_0(k,j,i),3) .AND. & 830 BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & 831 BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. & 832 .NOT. flag_set ) & 833 THEN 829 834 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 ) 830 835 ENDIF … … 837 842 CALL exchange_horiz_int( advc_flags_m, nys, nyn, nxl, nxr, nzt, nbgp ) 838 843 ! 839 !-- Set boundary flags at inflow and outflow boundary in case of 844 !-- Set boundary flags at inflow and outflow boundary in case of 840 845 !-- non-cyclic boundary conditions. 841 846 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN … … 844 849 845 850 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 846 advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr)851 advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr) 847 852 ENDIF 848 853 … … 863 868 !> Initialization of flags to control the order of the advection scheme near 864 869 !> solid walls and non-cyclic inflow boundaries, where the order is sucessively 865 !> degraded. 870 !> degraded. 866 871 !------------------------------------------------------------------------------! 867 872 SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, & … … 885 890 LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south 886 891 887 LOGICAL, OPTIONAL :: extensive_degrad !< flag indicating that extensive degradation is required, e.g. for 892 LOGICAL, OPTIONAL :: extensive_degrad !< flag indicating that extensive degradation is required, e.g. for 888 893 !< passive scalars nearby topography along the horizontal directions, 889 !< as no monotonic limiter can be applied there 894 !< as no monotonic limiter can be applied there 890 895 ! 891 896 !-- Set flags to steer the degradation of the advection scheme in advec_ws … … 900 905 !-- scalar - x-direction 901 906 !-- WS1 (0), WS3 (1), WS5 (2) 902 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0) .OR. & 903 .NOT. BTEST(wall_flags_total_0(k,j,i+2),0) .OR. & 904 .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) ) .OR. & 905 ( non_cyclic_l .AND. i == 0 ) .OR. & 906 ( non_cyclic_r .AND. i == nx ) ) THEN 907 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) 908 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & 909 BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & 910 BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & 911 BTEST(wall_flags_total_0(k,j,i-1),0) & 912 ) .OR. & 913 ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0) .AND. & 914 BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & 915 BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & 916 BTEST(wall_flags_total_0(k,j,i-1),0) & 917 ) & 918 .OR. & 919 ( non_cyclic_r .AND. i == nx-1 ) .OR. & 920 ( non_cyclic_l .AND. i == 1 ) ) THEN 921 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) 922 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & 923 BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & 924 BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & 925 BTEST(wall_flags_total_0(k,j,i-1),0) .AND. & 926 BTEST(wall_flags_total_0(k,j,i-2),0) ) & 927 THEN 928 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 ) 929 ENDIF 930 ! 931 !-- scalar - y-direction 932 !-- WS1 (3), WS3 (4), WS5 (5) 933 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0) .OR. & 934 .NOT. BTEST(wall_flags_total_0(k,j+2,i),0) .OR. & 935 .NOT. BTEST(wall_flags_total_0(k,j-1,i),0)) .OR. & 936 ( non_cyclic_s .AND. j == 0 ) .OR. & 937 ( non_cyclic_n .AND. j == ny ) ) THEN 938 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 939 ! 940 !-- WS3 941 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & 942 BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & 943 BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & 944 BTEST(wall_flags_total_0(k,j-1,i),0) & 945 ) .OR. & 946 ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0) .AND. & 947 BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & 948 BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & 949 BTEST(wall_flags_total_0(k,j-1,i),0) & 950 ) & 951 .OR. & 952 ( non_cyclic_s .AND. j == 1 ) .OR. & 953 ( non_cyclic_n .AND. j == ny-1 ) ) THEN 954 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) 955 ! 956 !-- WS5 957 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & 958 BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & 959 BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & 960 BTEST(wall_flags_total_0(k,j-1,i),0) .AND. & 961 BTEST(wall_flags_total_0(k,j-2,i),0) ) & 962 THEN 963 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 ) 907 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0) .OR. & 908 .NOT. BTEST(wall_flags_total_0(k,j,i+2),0) .OR. & 909 .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) ) .OR. & 910 ( non_cyclic_l .AND. i == 0 ) .OR. & 911 ( non_cyclic_r .AND. i == nx ) ) THEN 912 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) 913 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & 914 BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & 915 BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & 916 BTEST(wall_flags_total_0(k,j,i-1),0) & 917 ) .OR. & 918 ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0) .AND. & 919 BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & 920 BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & 921 BTEST(wall_flags_total_0(k,j,i-1),0) & 922 ) & 923 .OR. & 924 ( non_cyclic_r .AND. i == nx-1 ) .OR. & 925 ( non_cyclic_l .AND. i == 1 ) ) THEN 926 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) 927 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & 928 BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & 929 BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & 930 BTEST(wall_flags_total_0(k,j,i-1),0) .AND. & 931 BTEST(wall_flags_total_0(k,j,i-2),0) ) & 932 THEN 933 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 ) 964 934 ENDIF 965 935 ! 966 !-- Near topography, set horizontal advection scheme to 1st order 967 !-- for passive scalars, even if only one direction may be 936 !-- scalar - y-direction 937 !-- WS1 (3), WS3 (4), WS5 (5) 938 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0) .OR. & 939 .NOT. BTEST(wall_flags_total_0(k,j+2,i),0) .OR. & 940 .NOT. BTEST(wall_flags_total_0(k,j-1,i),0)) .OR. & 941 ( non_cyclic_s .AND. j == 0 ) .OR. & 942 ( non_cyclic_n .AND. j == ny ) ) THEN 943 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 944 ! 945 !-- WS3 946 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & 947 BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & 948 BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & 949 BTEST(wall_flags_total_0(k,j-1,i),0) & 950 ) .OR. & 951 ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0) .AND. & 952 BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & 953 BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & 954 BTEST(wall_flags_total_0(k,j-1,i),0) & 955 ) & 956 .OR. & 957 ( non_cyclic_s .AND. j == 1 ) .OR. & 958 ( non_cyclic_n .AND. j == ny-1 ) ) THEN 959 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) 960 ! 961 !-- WS5 962 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & 963 BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & 964 BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & 965 BTEST(wall_flags_total_0(k,j-1,i),0) .AND. & 966 BTEST(wall_flags_total_0(k,j-2,i),0) ) & 967 THEN 968 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 ) 969 ENDIF 970 ! 971 !-- Near topography, set horizontal advection scheme to 1st order 972 !-- for passive scalars, even if only one direction may be 968 973 !-- blocked by topography. These locations will be identified 969 !-- by wall_flags_total_0 bit 31. Note, since several modules define 970 !-- advection flags but may apply different scalar boundary 971 !-- conditions, bit 31 is temporarily stored on advc_flags. 972 !-- Moreover, note that this extended degradtion for passive 973 !-- scalars is not required for the vertical direction as there 974 !-- the monotonic limiter can be applied. 974 !-- by wall_flags_total_0 bit 31. Note, since several modules define 975 !-- advection flags but may apply different scalar boundary 976 !-- conditions, bit 31 is temporarily stored on advc_flags. 977 !-- Moreover, note that this extended degradtion for passive 978 !-- scalars is not required for the vertical direction as there 979 !-- the monotonic limiter can be applied. 975 980 IF ( PRESENT( extensive_degrad ) ) THEN 976 981 IF ( extensive_degrad ) THEN … … 980 985 IF( BTEST( advc_flag(k,j,i), 31 ) ) THEN 981 986 ! 982 !-- Clear flags that might indicate higher-order 987 !-- Clear flags that might indicate higher-order 983 988 !-- advection along x- and y-direction. 984 989 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) … … 990 995 !-- x- and y-direction. 991 996 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) 992 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 997 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 993 998 ENDIF 994 999 ! 995 1000 !-- Adjacent to this extended degradation zone, successively 996 !-- upgrade the order of the scheme if this grid point isn't 997 !-- flagged with bit 31 (indicating extended degradation 998 !-- zone). 1001 !-- upgrade the order of the scheme if this grid point isn't 1002 !-- flagged with bit 31 (indicating extended degradation 1003 !-- zone). 999 1004 IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) ) THEN 1000 1005 ! 1001 1006 !-- x-direction. First, clear all previous settings, than 1002 1007 !-- set flag for 3rd-order scheme. 1003 IF ( BTEST( advc_flag(k,j,i-1), 31 ) .AND. &1008 IF ( BTEST( advc_flag(k,j,i-1), 31 ) .AND. & 1004 1009 BTEST( advc_flag(k,j,i+1), 31 ) ) THEN 1005 1010 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1006 1011 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1007 1012 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1008 1013 1009 1014 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) 1010 1015 ENDIF … … 1012 1017 !-- x-direction. First, clear all previous settings, than 1013 1018 !-- set flag for 5rd-order scheme. 1014 IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 ) .AND. &1015 BTEST( advc_flag(k,j,i-2), 31 ) .AND. &1016 .NOT. BTEST( advc_flag(k,j,i+1), 31 ) .AND. &1019 IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 ) .AND. & 1020 BTEST( advc_flag(k,j,i-2), 31 ) .AND. & 1021 .NOT. BTEST( advc_flag(k,j,i+1), 31 ) .AND. & 1017 1022 BTEST( advc_flag(k,j,i+2), 31 ) ) THEN 1018 1023 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1019 1024 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1020 1025 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1021 1026 1022 1027 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 ) 1023 1028 ENDIF … … 1025 1030 !-- y-direction. First, clear all previous settings, than 1026 1031 !-- set flag for 3rd-order scheme. 1027 IF ( BTEST( advc_flag(k,j-1,i), 31 ) .AND. &1032 IF ( BTEST( advc_flag(k,j-1,i), 31 ) .AND. & 1028 1033 BTEST( advc_flag(k,j+1,i), 31 ) ) THEN 1029 1034 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1030 1035 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1031 1036 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1032 1037 1033 1038 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) 1034 1039 ENDIF … … 1036 1041 !-- y-direction. First, clear all previous settings, than 1037 1042 !-- set flag for 5rd-order scheme. 1038 IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 ) .AND. &1039 BTEST( advc_flag(k,j-2,i), 31 ) .AND. &1040 .NOT. BTEST( advc_flag(k,j+1,i), 31 ) .AND. &1043 IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 ) .AND. & 1044 BTEST( advc_flag(k,j-2,i), 31 ) .AND. & 1045 .NOT. BTEST( advc_flag(k,j+1,i), 31 ) .AND. & 1041 1046 BTEST( advc_flag(k,j+2,i), 31 ) ) THEN 1042 1047 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1043 1048 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1044 1049 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1045 1050 1046 1051 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 ) 1047 1052 ENDIF 1048 1053 ENDIF 1049 1054 1050 1055 ENDIF 1051 1056 1052 1057 ! 1053 1058 !-- Near lateral boundary flags might be overwritten. Set 1054 !-- them again. 1059 !-- them again. 1055 1060 !-- x-direction 1056 IF ( ( non_cyclic_l .AND. i == 0 ) .OR. &1061 IF ( ( non_cyclic_l .AND. i == 0 ) .OR. & 1057 1062 ( non_cyclic_r .AND. i == nx ) ) THEN 1058 1063 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1059 1064 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1060 1065 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1061 1066 1062 1067 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) 1063 1068 ENDIF 1064 1065 IF ( ( non_cyclic_l .AND. i == 1 ) .OR. &1069 1070 IF ( ( non_cyclic_l .AND. i == 1 ) .OR. & 1066 1071 ( non_cyclic_r .AND. i == nx-1 ) ) THEN 1067 1072 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1068 1073 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1069 1074 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1070 1075 1071 1076 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) 1072 1077 ENDIF 1073 1078 ! 1074 1079 !-- y-direction 1075 IF ( ( non_cyclic_n .AND. j == 0 ) .OR. &1080 IF ( ( non_cyclic_n .AND. j == 0 ) .OR. & 1076 1081 ( non_cyclic_s .AND. j == ny ) ) THEN 1077 1082 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1078 1083 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1079 1084 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1080 1085 1081 1086 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 1082 1087 ENDIF 1083 1084 IF ( ( non_cyclic_n .AND. j == 1 ) .OR. &1088 1089 IF ( ( non_cyclic_n .AND. j == 1 ) .OR. & 1085 1090 ( non_cyclic_s .AND. j == ny-1 ) ) THEN 1086 1091 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1087 1092 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1088 1093 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1089 1094 1090 1095 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) 1091 1096 ENDIF 1092 1097 1093 1098 ENDIF 1094 1095 1096 ! 1097 !-- scalar - z-direction. Fluxes are calculated on w-grid 1098 !-- level. Boundary values at/within walls aren't used. 1099 !-- WS1 (6), WS3 (7), WS5 (8) 1100 IF ( k == nzb+1 ) THEN 1101 k_mm = nzb 1102 ELSE 1103 k_mm = k - 2 1104 ENDIF 1105 IF ( k > nzt-1 ) THEN 1106 k_pp = nzt+1 1107 ELSE 1108 k_pp = k + 2 1109 ENDIF 1110 IF ( k > nzt-2 ) THEN 1111 k_ppp = nzt+1 1112 ELSE 1113 k_ppp = k + 3 1114 ENDIF 1115 1116 flag_set = .FALSE. 1117 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),0) 1118 BTEST(wall_flags_total_0(k,j,i),0) 1119 BTEST(wall_flags_total_0(k+1,j,i),0) ) 1120 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),0) .AND. &1121 BTEST(wall_flags_total_0(k+1,j,i),0) 1122 BTEST(wall_flags_total_0(k,j,i),0) ) 1123 ( k == nzt .AND. symmetry_flag == 0 ) ) 1124 THEN 1125 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 ) 1126 flag_set = .TRUE. 1127 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),0) .OR. 1128 .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),0) ) .AND. &1129 BTEST(wall_flags_total_0(k-1,j,i),0) .AND. 1130 BTEST(wall_flags_total_0(k,j,i),0) .AND. 1131 BTEST(wall_flags_total_0(k+1,j,i),0) .AND. 1132 BTEST(wall_flags_total_0(k_pp,j,i),0) .AND. 1133 .NOT. flag_set .OR. &1134 ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) 1135 THEN 1136 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 ) 1137 flag_set = .TRUE. 1138 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),0) 1139 BTEST(wall_flags_total_0(k-1,j,i),0) 1140 BTEST(wall_flags_total_0(k,j,i),0) 1141 BTEST(wall_flags_total_0(k+1,j,i),0) 1142 BTEST(wall_flags_total_0(k_pp,j,i),0) 1143 BTEST(wall_flags_total_0(k_ppp,j,i),0) 1144 .NOT. flag_set ) 1099 1100 1101 ! 1102 !-- scalar - z-direction. Fluxes are calculated on w-grid 1103 !-- level. Boundary values at/within walls aren't used. 1104 !-- WS1 (6), WS3 (7), WS5 (8) 1105 IF ( k == nzb+1 ) THEN 1106 k_mm = nzb 1107 ELSE 1108 k_mm = k - 2 1109 ENDIF 1110 IF ( k > nzt-1 ) THEN 1111 k_pp = nzt+1 1112 ELSE 1113 k_pp = k + 2 1114 ENDIF 1115 IF ( k > nzt-2 ) THEN 1116 k_ppp = nzt+1 1117 ELSE 1118 k_ppp = k + 3 1119 ENDIF 1120 1121 flag_set = .FALSE. 1122 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),0) .AND. & 1123 BTEST(wall_flags_total_0(k,j,i),0) .AND. & 1124 BTEST(wall_flags_total_0(k+1,j,i),0) ) .OR. & 1125 ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),0) .AND. & 1126 BTEST(wall_flags_total_0(k+1,j,i),0) .AND. & 1127 BTEST(wall_flags_total_0(k,j,i),0) ) .OR. & 1128 ( k == nzt .AND. symmetry_flag == 0 ) ) & 1129 THEN 1130 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 ) 1131 flag_set = .TRUE. 1132 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),0) .OR. & 1133 .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),0) ) .AND.& 1134 BTEST(wall_flags_total_0(k-1,j,i),0) .AND.& 1135 BTEST(wall_flags_total_0(k,j,i),0) .AND.& 1136 BTEST(wall_flags_total_0(k+1,j,i),0) .AND.& 1137 BTEST(wall_flags_total_0(k_pp,j,i),0) .AND.& 1138 .NOT. flag_set .OR.& 1139 ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) & 1140 THEN 1141 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 ) 1142 flag_set = .TRUE. 1143 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),0) .AND. & 1144 BTEST(wall_flags_total_0(k-1,j,i),0) .AND. & 1145 BTEST(wall_flags_total_0(k,j,i),0) .AND. & 1146 BTEST(wall_flags_total_0(k+1,j,i),0) .AND. & 1147 BTEST(wall_flags_total_0(k_pp,j,i),0) .AND. & 1148 BTEST(wall_flags_total_0(k_ppp,j,i),0) .AND. & 1149 .NOT. flag_set ) & 1145 1150 THEN 1146 1151 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 ) … … 1156 1161 CALL exchange_horiz_int( advc_flag, nys, nyn, nxl, nxr, nzt, nbgp ) 1157 1162 ! 1158 !-- Set boundary flags at inflow and outflow boundary in case of 1163 !-- Set boundary flags at inflow and outflow boundary in case of 1159 1164 !-- non-cyclic boundary conditions. 1160 1165 IF ( non_cyclic_l ) THEN … … 1173 1178 advc_flag(:,nys-1,:) = advc_flag(:,nys,:) 1174 1179 ENDIF 1175 1176 1177 1178 END SUBROUTINE ws_init_flags_scalar 1179 1180 1181 1182 1183 END SUBROUTINE ws_init_flags_scalar 1184 1180 1185 !------------------------------------------------------------------------------! 1181 1186 ! Description: … … 1187 1192 1188 1193 ! 1189 !-- The arrays needed for statistical evaluation are set to to 0 at the 1194 !-- The arrays needed for statistical evaluation are set to to 0 at the 1190 1195 !-- beginning of prognostic_equations. 1191 1196 IF ( ws_scheme_mom ) THEN … … 1225 1230 1226 1231 1227 CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array 1228 1232 CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the 1233 !<correct dimension in the analysis array 1234 1229 1235 INTEGER(iwp) :: i !< grid index along x-direction 1230 1236 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread … … 1235 1241 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 1236 1242 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 1237 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 1243 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 1238 1244 INTEGER(iwp) :: tn !< number of OpenMP thread 1239 1245 1240 1246 INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: & 1241 1247 advc_flag !< flag array to control order of scalar advection … … 1244 1250 LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north 1245 1251 LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right 1246 LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south 1252 LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south 1247 1253 LOGICAL, OPTIONAL :: flux_limitation !< flag indicating flux limitation of the vertical advection 1248 1254 LOGICAL :: limiter !< control flag indicating the application of flux limitation … … 1251 1257 REAL(wp) :: div !< velocity diverence on scalar grid 1252 1258 REAL(wp) :: div_in !< vertical flux divergence of ingoing fluxes 1253 REAL(wp) :: div_out !< vertical flux divergence of outgoing fluxes 1259 REAL(wp) :: div_out !< vertical flux divergence of outgoing fluxes 1254 1260 REAL(wp) :: f_corr_t !< correction flux at grid-cell top, i.e. the difference between high and low-order flux 1255 1261 REAL(wp) :: f_corr_d !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux … … 1272 1278 REAL(wp) :: min_val !< maximum value of the quanitity along the numerical stencil (in vertical direction) 1273 1279 REAL(wp) :: mon !< monotone solution of the advection equation using 1st-order fluxes 1274 REAL(wp) :: u_comp !< advection velocity along x-direction 1275 REAL(wp) :: v_comp !< advection velocity along y-direction 1276 ! 1277 !-- sk is an array from parameter list. It should not be a pointer, because 1278 !-- in that case the compiler can not assume a stride 1 and cannot perform 1279 !-- a strided one vector load. Adding the CONTIGUOUS keyword makes things 1280 !-- even worse, because the compiler cannot assume strided one in the 1280 REAL(wp) :: u_comp !< advection velocity along x-direction 1281 REAL(wp) :: v_comp !< advection velocity along y-direction 1282 ! 1283 !-- sk is an array from parameter list. It should not be a pointer, because 1284 !-- in that case the compiler can not assume a stride 1 and cannot perform 1285 !-- a strided one vector load. Adding the CONTIGUOUS keyword makes things 1286 !-- even worse, because the compiler cannot assume strided one in the 1281 1287 !-- caller side. 1282 1288 REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar 1283 1289 1284 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 1285 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 1286 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box 1287 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 1288 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 1289 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box 1290 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t_1st !< discretized 1st-order flux at top of the grid box 1291 1292 REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !< discretized artificial dissipation at southward-side of the grid box 1293 REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !< discretized 6th-order flux at northward-side of the grid box 1294 1295 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !< discretized artificial dissipation at leftward-side of the grid box 1296 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !< discretized 6th-order flux at leftward-side of the grid box 1297 ! 1298 !-- Used local modified copy of nzb_max (used to degrade order of 1290 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side 1291 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side 1292 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top 1293 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side 1294 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side 1295 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top 1296 REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t_1st !< discretized 1st-order flux at top 1297 1298 REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !< discretized artificial dissipation at southward-side 1299 REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !< discretized 6th-order flux at northward-side 1300 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !< discretized artificial dissipation at leftward-side 1301 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !< discretized 6th-order flux at leftward-side 1302 ! 1303 !-- Used local modified copy of nzb_max (used to degrade order of 1299 1304 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 1300 !-- instead of the entire subdomain. This should lead to better 1305 !-- instead of the entire subdomain. This should lead to better 1301 1306 !-- load balance between boundary and non-boundary PEs. 1302 1307 IF( non_cyclic_l .AND. i <= nxl + 2 .OR. & … … 1309 1314 END IF 1310 1315 ! 1311 !-- Set control flag for flux limiter 1316 !-- Set control flag for flux limiter 1312 1317 limiter = .FALSE. 1313 1318 IF ( PRESENT( flux_limitation) ) limiter = flux_limitation … … 1324 1329 1325 1330 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) 1326 swap_flux_y_local(k,tn) = v_comp * ( &1327 ( 37.0_wp * ibit5 * adv_sca_5 &1328 + 7.0_wp * ibit4 * adv_sca_3 &1329 + ibit3 * adv_sca_1 &1330 ) * &1331 ( sk(k,j,i) + sk(k,j-1,i) ) &1332 - ( 8.0_wp * ibit5 * adv_sca_5 &1333 + ibit4 * adv_sca_3 &1334 ) * &1335 ( sk(k,j+1,i) + sk(k,j-2,i) ) &1336 + ( ibit5 * adv_sca_5 &1337 ) * &1338 ( sk(k,j+2,i) + sk(k,j-3,i) ) &1331 swap_flux_y_local(k,tn) = v_comp * ( & 1332 ( 37.0_wp * ibit5 * adv_sca_5 & 1333 + 7.0_wp * ibit4 * adv_sca_3 & 1334 + ibit3 * adv_sca_1 & 1335 ) * & 1336 ( sk(k,j,i) + sk(k,j-1,i) ) & 1337 - ( 8.0_wp * ibit5 * adv_sca_5 & 1338 + ibit4 * adv_sca_3 & 1339 ) * & 1340 ( sk(k,j+1,i) + sk(k,j-2,i) ) & 1341 + ( ibit5 * adv_sca_5 & 1342 ) * & 1343 ( sk(k,j+2,i) + sk(k,j-3,i) ) & 1339 1344 ) 1340 1345 1341 swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &1342 ( 10.0_wp * ibit5 * adv_sca_5 &1343 + 3.0_wp * ibit4 * adv_sca_3 &1344 + ibit3 * adv_sca_1 &1345 ) * &1346 ( sk(k,j,i) - sk(k,j-1,i) ) &1347 - ( 5.0_wp * ibit5 * adv_sca_5 &1348 + ibit4 * adv_sca_3 &1349 ) * &1350 ( sk(k,j+1,i) - sk(k,j-2,i) ) &1351 + ( ibit5 * adv_sca_5 &1352 ) * &1353 ( sk(k,j+2,i) - sk(k,j-3,i) ) &1346 swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & 1347 ( 10.0_wp * ibit5 * adv_sca_5 & 1348 + 3.0_wp * ibit4 * adv_sca_3 & 1349 + ibit3 * adv_sca_1 & 1350 ) * & 1351 ( sk(k,j,i) - sk(k,j-1,i) ) & 1352 - ( 5.0_wp * ibit5 * adv_sca_5 & 1353 + ibit4 * adv_sca_3 & 1354 ) * & 1355 ( sk(k,j+1,i) - sk(k,j-2,i) ) & 1356 + ( ibit5 * adv_sca_5 & 1357 ) * & 1358 ( sk(k,j+2,i) - sk(k,j-3,i) ) & 1354 1359 ) 1355 1360 … … 1360 1365 1361 1366 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) 1362 swap_flux_y_local(k,tn) = v_comp * ( &1363 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &1364 - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &1365 + ( sk(k,j+2,i) + sk(k,j-3,i) ) &1367 swap_flux_y_local(k,tn) = v_comp * ( & 1368 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & 1369 - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) & 1370 + ( sk(k,j+2,i) + sk(k,j-3,i) ) & 1366 1371 ) * adv_sca_5 1367 swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &1368 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &1369 - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &1370 + sk(k,j+2,i) - sk(k,j-3,i) &1372 swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & 1373 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & 1374 - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & 1375 + sk(k,j+2,i) - sk(k,j-3,i) & 1371 1376 ) * adv_sca_5 1372 1377 … … 1377 1382 !-- Compute leftside fluxes of the respective PE bounds. 1378 1383 IF ( i == i_omp ) THEN 1379 1384 1380 1385 DO k = nzb+1, nzb_max_l 1381 1386 … … 1384 1389 ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) 1385 1390 1386 u_comp 1387 swap_flux_x_local(k,j,tn) = u_comp * ( &1388 ( 37.0_wp * ibit2 * adv_sca_5 &1389 + 7.0_wp * ibit1 * adv_sca_3 &1390 + ibit0 * adv_sca_1 &1391 ) * &1392 ( sk(k,j,i) + sk(k,j,i-1) ) &1393 - ( 8.0_wp * ibit2 * adv_sca_5 &1394 + ibit1 * adv_sca_3 &1395 ) * &1396 ( sk(k,j,i+1) + sk(k,j,i-2) ) &1397 + ( ibit2 * adv_sca_5 &1398 ) * &1399 ( sk(k,j,i+2) + sk(k,j,i-3) ) &1391 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) 1392 swap_flux_x_local(k,j,tn) = u_comp * ( & 1393 ( 37.0_wp * ibit2 * adv_sca_5 & 1394 + 7.0_wp * ibit1 * adv_sca_3 & 1395 + ibit0 * adv_sca_1 & 1396 ) * & 1397 ( sk(k,j,i) + sk(k,j,i-1) ) & 1398 - ( 8.0_wp * ibit2 * adv_sca_5 & 1399 + ibit1 * adv_sca_3 & 1400 ) * & 1401 ( sk(k,j,i+1) + sk(k,j,i-2) ) & 1402 + ( ibit2 * adv_sca_5 & 1403 ) * & 1404 ( sk(k,j,i+2) + sk(k,j,i-3) ) & 1400 1405 ) 1401 1406 1402 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &1403 ( 10.0_wp * ibit2 * adv_sca_5 &1404 + 3.0_wp * ibit1 * adv_sca_3 &1405 + ibit0 * adv_sca_1 &1406 ) * &1407 ( sk(k,j,i) - sk(k,j,i-1) ) &1408 - ( 5.0_wp * ibit2 * adv_sca_5 &1409 + ibit1 * adv_sca_3 &1410 ) * &1411 ( sk(k,j,i+1) - sk(k,j,i-2) ) &1412 + ( ibit2 * adv_sca_5 &1413 ) * &1414 ( sk(k,j,i+2) - sk(k,j,i-3) ) &1407 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & 1408 ( 10.0_wp * ibit2 * adv_sca_5 & 1409 + 3.0_wp * ibit1 * adv_sca_3 & 1410 + ibit0 * adv_sca_1 & 1411 ) * & 1412 ( sk(k,j,i) - sk(k,j,i-1) ) & 1413 - ( 5.0_wp * ibit2 * adv_sca_5 & 1414 + ibit1 * adv_sca_3 & 1415 ) * & 1416 ( sk(k,j,i+1) - sk(k,j,i-2) ) & 1417 + ( ibit2 * adv_sca_5 & 1418 ) * & 1419 ( sk(k,j,i+2) - sk(k,j,i-3) ) & 1415 1420 ) 1416 1421 … … 1419 1424 DO k = nzb_max_l+1, nzt 1420 1425 1421 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k)1422 swap_flux_x_local(k,j,tn) = u_comp * ( &1423 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &1424 - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &1425 + ( sk(k,j,i+2) + sk(k,j,i-3) ) &1426 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) 1427 swap_flux_x_local(k,j,tn) = u_comp * ( & 1428 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & 1429 - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) & 1430 + ( sk(k,j,i+2) + sk(k,j,i-3) ) & 1426 1431 ) * adv_sca_5 1427 1432 1428 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &1429 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &1430 - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &1431 + ( sk(k,j,i+2) - sk(k,j,i-3) ) &1433 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & 1434 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) & 1435 - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) & 1436 + ( sk(k,j,i+2) - sk(k,j,i-3) ) & 1432 1437 ) * adv_sca_5 1433 1438 1434 1439 ENDDO 1435 1440 1436 1441 ENDIF 1437 ! 1438 !-- Now compute the fluxes and tendency terms for the horizontal and1439 !-- vertical parts up to the top of the highesttopography.1442 ! 1443 !-- Now compute the fluxes for the horizontal termns up to the highest 1444 !-- topography. 1440 1445 DO k = nzb+1, nzb_max_l 1441 !1442 !-- Note: It is faster to conduct all multiplications explicitly, e.g.1443 !-- * adv_sca_5 ... than to determine a factor and multiplicate the1444 !-- flux at the end.1445 1446 1446 1447 ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) … … 1449 1450 1450 1451 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) 1451 flux_r(k) = u_comp * ( &1452 ( 37.0_wp * ibit2 * adv_sca_5 &1453 + 7.0_wp * ibit1 * adv_sca_3 &1454 + ibit0 * adv_sca_1 &1455 ) * &1456 ( sk(k,j,i+1) + sk(k,j,i) ) &1457 - ( 8.0_wp * ibit2 * adv_sca_5 &1458 + ibit1 * adv_sca_3 &1459 ) * &1460 ( sk(k,j,i+2) + sk(k,j,i-1) ) &1461 + ( ibit2 * adv_sca_5 &1462 ) * &1463 ( sk(k,j,i+3) + sk(k,j,i-2) ) &1452 flux_r(k) = u_comp * ( & 1453 ( 37.0_wp * ibit2 * adv_sca_5 & 1454 + 7.0_wp * ibit1 * adv_sca_3 & 1455 + ibit0 * adv_sca_1 & 1456 ) * & 1457 ( sk(k,j,i+1) + sk(k,j,i) ) & 1458 - ( 8.0_wp * ibit2 * adv_sca_5 & 1459 + ibit1 * adv_sca_3 & 1460 ) * & 1461 ( sk(k,j,i+2) + sk(k,j,i-1) ) & 1462 + ( ibit2 * adv_sca_5 & 1463 ) * & 1464 ( sk(k,j,i+3) + sk(k,j,i-2) ) & 1464 1465 ) 1465 1466 1466 diss_r(k) = -ABS( u_comp ) * ( &1467 ( 10.0_wp * ibit2 * adv_sca_5 &1468 + 3.0_wp * ibit1 * adv_sca_3 &1469 + ibit0 * adv_sca_1 &1470 ) * &1471 ( sk(k,j,i+1) - sk(k,j,i) ) &1472 - ( 5.0_wp * ibit2 * adv_sca_5 &1473 + ibit1 * adv_sca_3 &1474 ) * &1475 ( sk(k,j,i+2) - sk(k,j,i-1) ) &1476 + ( ibit2 * adv_sca_5 &1477 ) * &1478 ( sk(k,j,i+3) - sk(k,j,i-2) ) &1467 diss_r(k) = -ABS( u_comp ) * ( & 1468 ( 10.0_wp * ibit2 * adv_sca_5 & 1469 + 3.0_wp * ibit1 * adv_sca_3 & 1470 + ibit0 * adv_sca_1 & 1471 ) * & 1472 ( sk(k,j,i+1) - sk(k,j,i) ) & 1473 - ( 5.0_wp * ibit2 * adv_sca_5 & 1474 + ibit1 * adv_sca_3 & 1475 ) * & 1476 ( sk(k,j,i+2) - sk(k,j,i-1) ) & 1477 + ( ibit2 * adv_sca_5 & 1478 ) * & 1479 ( sk(k,j,i+3) - sk(k,j,i-2) ) & 1479 1480 ) 1480 1481 … … 1515 1516 ENDDO 1516 1517 ! 1517 !-- Now compute the fluxes and tendency terms for the horizontal and1518 !-- vertical parts above the top of the highest topography. No degradation1519 !-- for the horizontal parts, but for the vertical it is stell needed.1518 !-- Now compute the fluxes for the horizontal terms above the topography 1519 !-- where no degradation along the horizontal parts is necessary (except 1520 !-- for the non-cyclic lateral boundaries treated by nzb_max_l). 1520 1521 DO k = nzb_max_l+1, nzt 1521 1522 … … 1542 1543 ENDDO 1543 1544 ! 1544 !-- Now, compute vertical fluxes. Split loop into a part treating the 1545 !-- Now, compute vertical fluxes. Split loop into a part treating the 1545 1546 !-- lowest grid points with indirect indexing, a main loop without 1546 1547 !-- indirect indexing, and a loop for the uppermost grip points with 1547 1548 !-- indirect indexing. This allows better vectorization for the main loop. 1548 !-- First, compute the flux at model surface, which need has to be 1549 !-- First, compute the flux at model surface, which need has to be 1549 1550 !-- calculated explicetely for the tendency at 1550 1551 !-- the first w-level. For topography wall this is done implicitely by … … 1552 1553 flux_t(nzb) = 0.0_wp 1553 1554 diss_t(nzb) = 0.0_wp 1554 1555 1555 1556 DO k = nzb+1, nzb+1 1556 1557 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) … … 1563 1564 k_pp = k + 2 * ( 1 - ibit6 ) 1564 1565 k_mm = k - 2 * ibit8 1565 1566 1567 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &1568 ( 37.0_wp * ibit8 * adv_sca_5 &1569 + 7.0_wp * ibit7 * adv_sca_3 &1570 + ibit6 * adv_sca_1 &1571 ) * &1572 ( sk(k+1,j,i) + sk(k,j,i) ) &1573 - ( 8.0_wp * ibit8 * adv_sca_5 &1574 + ibit7 * adv_sca_3 &1575 ) * &1576 ( sk(k_pp,j,i) + sk(k-1,j,i) ) &1577 + ( ibit8 * adv_sca_5 &1578 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &1579 )1580 1581 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &1582 ( 10.0_wp * ibit8 * adv_sca_5 &1583 + 3.0_wp * ibit7 * adv_sca_3 &1584 + ibit6 * adv_sca_1 &1585 ) * &1586 ( sk(k+1,j,i) - sk(k,j,i) ) &1587 - ( 5.0_wp * ibit8 * adv_sca_5 &1588 + ibit7 * adv_sca_3 &1589 ) * &1590 ( sk(k_pp,j,i) - sk(k-1,j,i) ) &1591 + ( ibit8 * adv_sca_5 &1592 ) * &1593 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &1594 )1595 ENDDO1596 1597 DO k = nzb+2, nzt-21598 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )1599 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )1600 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )1601 1602 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &1603 ( 37.0_wp * ibit8 * adv_sca_5 &1604 + 7.0_wp * ibit7 * adv_sca_3 &1605 + ibit6 * adv_sca_1 &1606 ) * &1607 ( sk(k+1,j,i) + sk(k,j,i) ) &1608 - ( 8.0_wp * ibit8 * adv_sca_5 &1609 + ibit7 * adv_sca_3 &1610 ) * &1611 ( sk(k+2,j,i) + sk(k-1,j,i) ) &1612 + ( ibit8 * adv_sca_5 &1613 ) * ( sk(k+3,j,i)+ sk(k-2,j,i) ) &1614 )1615 1616 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &1617 ( 10.0_wp * ibit8 * adv_sca_5 &1618 + 3.0_wp * ibit7 * adv_sca_3 &1619 + ibit6 * adv_sca_1 &1620 ) * &1621 ( sk(k+1,j,i) - sk(k,j,i) ) &1622 - ( 5.0_wp * ibit8 * adv_sca_5 &1623 + ibit7 * adv_sca_3 &1624 ) * &1625 ( sk(k+2,j,i) - sk(k-1,j,i) ) &1626 + ( ibit8 * adv_sca_5 &1627 ) * &1628 ( sk(k+3,j,i) - sk(k-2,j,i) ) &1629 )1630 ENDDO1631 1632 DO k = nzt-1, nzt-symmetry_flag1633 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )1634 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )1635 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )1636 !1637 !-- k index has to be modified near bottom and top, else array1638 !-- subscripts will be exceeded.1639 k_ppp = k + 3 * ibit81640 k_pp = k + 2 * ( 1 - ibit6 )1641 k_mm = k - 2 * ibit81642 1643 1566 1644 1567 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & … … 1671 1594 ) 1672 1595 ENDDO 1673 1596 1597 DO k = nzb+2, nzt-2 1598 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 1599 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 1600 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 1601 1602 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 1603 ( 37.0_wp * ibit8 * adv_sca_5 & 1604 + 7.0_wp * ibit7 * adv_sca_3 & 1605 + ibit6 * adv_sca_1 & 1606 ) * & 1607 ( sk(k+1,j,i) + sk(k,j,i) ) & 1608 - ( 8.0_wp * ibit8 * adv_sca_5 & 1609 + ibit7 * adv_sca_3 & 1610 ) * & 1611 ( sk(k+2,j,i) + sk(k-1,j,i) ) & 1612 + ( ibit8 * adv_sca_5 & 1613 ) * ( sk(k+3,j,i)+ sk(k-2,j,i) ) & 1614 ) 1615 1616 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 1617 ( 10.0_wp * ibit8 * adv_sca_5 & 1618 + 3.0_wp * ibit7 * adv_sca_3 & 1619 + ibit6 * adv_sca_1 & 1620 ) * & 1621 ( sk(k+1,j,i) - sk(k,j,i) ) & 1622 - ( 5.0_wp * ibit8 * adv_sca_5 & 1623 + ibit7 * adv_sca_3 & 1624 ) * & 1625 ( sk(k+2,j,i) - sk(k-1,j,i) ) & 1626 + ( ibit8 * adv_sca_5 & 1627 ) * & 1628 ( sk(k+3,j,i) - sk(k-2,j,i) ) & 1629 ) 1630 ENDDO 1631 1632 DO k = nzt-1, nzt-symmetry_flag 1633 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 1634 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 1635 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 1636 ! 1637 !-- k index has to be modified near bottom and top, else array 1638 !-- subscripts will be exceeded. 1639 k_ppp = k + 3 * ibit8 1640 k_pp = k + 2 * ( 1 - ibit6 ) 1641 k_mm = k - 2 * ibit8 1642 1643 1644 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 1645 ( 37.0_wp * ibit8 * adv_sca_5 & 1646 + 7.0_wp * ibit7 * adv_sca_3 & 1647 + ibit6 * adv_sca_1 & 1648 ) * & 1649 ( sk(k+1,j,i) + sk(k,j,i) ) & 1650 - ( 8.0_wp * ibit8 * adv_sca_5 & 1651 + ibit7 * adv_sca_3 & 1652 ) * & 1653 ( sk(k_pp,j,i) + sk(k-1,j,i) ) & 1654 + ( ibit8 * adv_sca_5 & 1655 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & 1656 ) 1657 1658 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 1659 ( 10.0_wp * ibit8 * adv_sca_5 & 1660 + 3.0_wp * ibit7 * adv_sca_3 & 1661 + ibit6 * adv_sca_1 & 1662 ) * & 1663 ( sk(k+1,j,i) - sk(k,j,i) ) & 1664 - ( 5.0_wp * ibit8 * adv_sca_5 & 1665 + ibit7 * adv_sca_3 & 1666 ) * & 1667 ( sk(k_pp,j,i) - sk(k-1,j,i) ) & 1668 + ( ibit8 * adv_sca_5 & 1669 ) * & 1670 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1671 ) 1672 ENDDO 1673 1674 1674 ! 1675 1675 !-- Set resolved/turbulent flux at model top to zero (w-level). In case that 1676 1676 !-- a symmetric behavior between bottom and top shall be guaranteed (closed 1677 !-- channel flow), the flux at nzt is also set to zero. 1677 !-- channel flow), the flux at nzt is also set to zero. 1678 1678 IF ( symmetry_flag == 1 ) THEN 1679 1679 flux_t(nzt) = 0.0_wp … … 1682 1682 flux_t(nzt+1) = 0.0_wp 1683 1683 diss_t(nzt+1) = 0.0_wp 1684 1685 1684 1685 1686 1686 IF ( limiter ) THEN 1687 1687 ! … … 1695 1695 ! 1696 1696 !-- In flux limitation the total flux will be corrected. For the sake 1697 !-- of cleariness the higher-order advective and disspative fluxes 1698 !-- will be merged onto flux_t. 1697 !-- of cleariness the higher-order advective and disspative fluxes 1698 !-- will be merged onto flux_t. 1699 1699 flux_t(k) = flux_t(k) + diss_t(k) 1700 1700 diss_t(k) = 0.0_wp 1701 1701 ENDDO 1702 1702 ! 1703 !-- Flux limitation of vertical fluxes according to Skamarock (2006). 1703 !-- Flux limitation of vertical fluxes according to Skamarock (2006). 1704 1704 !-- Please note, as flux limitation implies linear dependencies of fluxes, 1705 1705 !-- flux limitation is only made for the vertical advection term. Limitation 1706 !-- of the horizontal terms cannot be parallelized. 1706 !-- of the horizontal terms cannot be parallelized. 1707 1707 !-- Due to the linear dependency, the following loop will not be vectorized. 1708 1708 !-- Further, note that the flux limiter is only applied within the urban 1709 !-- layer, i.e up to the topography top. 1709 !-- layer, i.e up to the topography top. 1710 1710 DO k = nzb+1, nzb_max_l 1711 1711 ! 1712 1712 !-- Compute one-dimensional divergence along the vertical direction, 1713 !-- which is used to correct the advection discretization. This is 1713 !-- which is used to correct the advection discretization. This is 1714 1714 !-- necessary as in one-dimensional space the advection velocity 1715 !-- should actually be constant. 1715 !-- should actually be constant. 1716 1716 div = ( w(k,j,i) * rho_air_zw(k) & 1717 - w(k-1,j,i) * rho_air_zw(k-1) & 1717 - w(k-1,j,i) * rho_air_zw(k-1) & 1718 1718 ) * drho_air(k) * ddzw(k) 1719 1719 ! 1720 !-- Compute monotone solution of the advection equation from 1720 !-- Compute monotone solution of the advection equation from 1721 1721 !-- 1st-order fluxes. Please note, the advection equation is corrected 1722 1722 !-- by the divergence term (in 1D the advective flow should be divergence 1723 !-- free). Moreover, please note, as time-increment the full timestep 1724 !-- is used, even though a Runge-Kutta scheme will be used. However, 1725 !-- the length of the actual time increment is not important at all 1726 !-- since it cancels out later when the fluxes are limited. 1723 !-- free). Moreover, please note, as time-increment the full timestep 1724 !-- is used, even though a Runge-Kutta scheme will be used. However, 1725 !-- the length of the actual time increment is not important at all 1726 !-- since it cancels out later when the fluxes are limited. 1727 1727 mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) ) & 1728 1728 * drho_air(k) * ddzw(k) & 1729 1729 + div * sk(k,j,i) & 1730 ) * dt_3d 1730 ) * dt_3d 1731 1731 ! 1732 1732 !-- Determine minimum and maximum values along the numerical stencil. 1733 1733 k_mmm = MAX( k - 3, nzb + 1 ) 1734 k_ppp = MIN( k + 3, nzt + 1 ) 1734 k_ppp = MIN( k + 3, nzt + 1 ) 1735 1735 1736 1736 min_val = MINVAL( sk(k_mmm:k_ppp,j,i) ) 1737 1737 max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) ) 1738 1738 ! 1739 !-- Compute difference between high- and low-order fluxes, which may 1739 !-- Compute difference between high- and low-order fluxes, which may 1740 1740 !-- act as correction fluxes 1741 1741 f_corr_t = flux_t(k) - flux_t_1st(k) 1742 1742 f_corr_d = flux_t(k-1) - flux_t_1st(k-1) 1743 1743 ! 1744 !-- Determine outgoing fluxes, i.e. the part of the fluxes which can 1744 !-- Determine outgoing fluxes, i.e. the part of the fluxes which can 1745 1745 !-- decrease the value within the grid box 1746 1746 f_corr_t_out = MAX( 0.0_wp, f_corr_t ) 1747 1747 f_corr_d_out = MIN( 0.0_wp, f_corr_d ) 1748 1748 ! 1749 !-- Determine ingoing fluxes, i.e. the part of the fluxes which can 1750 !-- increase the value within the grid box 1749 !-- Determine ingoing fluxes, i.e. the part of the fluxes which can 1750 !-- increase the value within the grid box 1751 1751 f_corr_t_in = MIN( 0.0_wp, f_corr_t) 1752 1752 f_corr_d_in = MAX( 0.0_wp, f_corr_d) … … 1762 1762 !-- Check if outgoing fluxes can lead to undershoots, i.e. values smaller 1763 1763 !-- than the minimum value within the numerical stencil. If so, limit 1764 !-- them. 1764 !-- them. 1765 1765 IF ( mon - min_val < - div_out .AND. ABS( div_out ) > 0.0_wp ) & 1766 1766 THEN … … 1772 1772 !-- Check if ingoing fluxes can lead to overshoots, i.e. values larger 1773 1773 !-- than the maximum value within the numerical stencil. If so, limit 1774 !-- them. 1774 !-- them. 1775 1775 IF ( mon - max_val > - div_in .AND. ABS( div_in ) > 0.0_wp ) & 1776 1776 THEN … … 1780 1780 ENDIF 1781 1781 ! 1782 !-- Finally add the limited fluxes to the original ones. If no 1783 !-- flux limitation was done, the fluxes equal the original ones. 1782 !-- Finally add the limited fluxes to the original ones. If no 1783 !-- flux limitation was done, the fluxes equal the original ones. 1784 1784 flux_t(k) = flux_t_1st(k) + f_corr_t_out + f_corr_t_in 1785 1785 flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in 1786 1786 ENDDO 1787 1787 ENDIF 1788 1788 ! 1789 !-- Now compute the tendency term including divergence correction. 1789 1790 DO k = nzb+1, nzb_max_l 1790 1791 1791 1792 flux_d = flux_t(k-1) 1792 1793 diss_d = diss_t(k-1) 1793 1794 1794 1795 ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) 1795 1796 ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) 1796 1797 ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) 1797 1798 1798 1799 ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) 1799 1800 ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) 1800 1801 ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) 1801 1802 1802 1803 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 1803 1804 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) … … 1807 1808 !-- correction is needed to overcome numerical instabilities introduced 1808 1809 !-- by a not sufficient reduction of divergences near topography. 1809 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &1810 - u(k,j,i) * ( &1811 REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) &1812 + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) &1813 + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) &1814 ) &1815 ) * ddx &1816 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &1817 - v(k,j,i) * ( &1818 REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) &1819 + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) &1820 + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) &1821 ) &1822 ) * ddy &1823 + ( w(k,j,i) * rho_air_zw(k) * &1824 ( ibit6 + ibit7 + ibit8 ) &1825 - w(k-1,j,i) * rho_air_zw(k-1) * &1826 ( &1827 REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp ) &1828 + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp ) &1829 + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp ) &1830 ) &1810 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 1811 - u(k,j,i) * ( & 1812 REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) & 1813 + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) & 1814 + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) & 1815 ) & 1816 ) * ddx & 1817 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 1818 - v(k,j,i) * ( & 1819 REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) & 1820 + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) & 1821 + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) & 1822 ) & 1823 ) * ddy & 1824 + ( w(k,j,i) * rho_air_zw(k) * & 1825 ( ibit6 + ibit7 + ibit8 ) & 1826 - w(k-1,j,i) * rho_air_zw(k-1) * & 1827 ( & 1828 REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp ) & 1829 + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp ) & 1830 + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp ) & 1831 ) & 1831 1832 ) * drho_air(k) * ddzw(k) 1832 1833 1833 tend(k,j,i) = tend(k,j,i) - ( &1834 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &1835 swap_diss_x_local(k,j,tn) ) * ddx &1836 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &1837 swap_diss_y_local(k,tn) ) * ddy &1838 + ( ( flux_t(k) + diss_t(k) ) - &1839 ( flux_d + diss_d ) &1840 ) * drho_air(k) * ddzw(k) &1834 tend(k,j,i) = tend(k,j,i) - ( & 1835 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & 1836 swap_diss_x_local(k,j,tn) ) * ddx & 1837 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & 1838 swap_diss_y_local(k,tn) ) * ddy & 1839 + ( ( flux_t(k) + diss_t(k) ) - & 1840 ( flux_d + diss_d ) & 1841 ) * drho_air(k) * ddzw(k) & 1841 1842 ) + sk(k,j,i) * div 1842 1843 … … 1848 1849 1849 1850 ENDDO 1850 1851 1851 1852 DO k = nzb_max_l+1, nzt 1852 1853 … … 1857 1858 !-- correction is needed to overcome numerical instabilities introduced 1858 1859 !-- by a not sufficient reduction of divergences near topography. 1859 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &1860 + ( v(k,j+1,i) - v(k,j,i) ) * ddy &1861 + ( w(k,j,i) * rho_air_zw(k) &1862 - w(k-1,j,i) * rho_air_zw(k-1) &1860 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 1861 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 1862 + ( w(k,j,i) * rho_air_zw(k) & 1863 - w(k-1,j,i) * rho_air_zw(k-1) & 1863 1864 ) * drho_air(k) * ddzw(k) 1864 1865 1865 tend(k,j,i) = tend(k,j,i) - ( &1866 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &1867 swap_diss_x_local(k,j,tn) ) * ddx &1868 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &1869 swap_diss_y_local(k,tn) ) * ddy &1870 + ( ( flux_t(k) + diss_t(k) ) - &1871 ( flux_d + diss_d ) &1872 ) * drho_air(k) * ddzw(k) &1866 tend(k,j,i) = tend(k,j,i) - ( & 1867 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & 1868 swap_diss_x_local(k,j,tn) ) * ddx & 1869 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & 1870 swap_diss_y_local(k,tn) ) * ddy & 1871 + ( ( flux_t(k) + diss_t(k) ) - & 1872 ( flux_d + diss_d ) & 1873 ) * drho_air(k) * ddzw(k) & 1873 1874 ) + sk(k,j,i) * div 1874 1875 … … 1895 1896 ) * weight_substep(intermediate_timestep_count) 1896 1897 ENDDO 1897 1898 1898 1899 CASE ( 'sa' ) 1899 1900 … … 1906 1907 ) * weight_substep(intermediate_timestep_count) 1907 1908 ENDDO 1908 1909 1909 1910 CASE ( 'q' ) 1910 1911 … … 1988 1989 !kk Has to be implemented for kpp chemistry 1989 1990 1990 1991 1991 END SELECT 1992 1992 … … 2011 2011 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 2012 2012 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 2013 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 2013 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 2014 2014 INTEGER(iwp) :: tn !< number of OpenMP thread 2015 2015 2016 2016 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction 2017 2017 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction … … 2029 2029 REAL(wp) :: gv !< Galilei-transformation velocity along y 2030 2030 REAL(wp) :: u_comp_l !< advection velocity along x at leftmost grid point on subdomain 2031 2031 2032 2032 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 2033 2033 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box … … 2040 2040 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z 2041 2041 ! 2042 !-- Used local modified copy of nzb_max (used to degrade order of 2042 !-- Used local modified copy of nzb_max (used to degrade order of 2043 2043 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 2044 !-- instead of the entire subdomain. This should lead to better 2044 !-- instead of the entire subdomain. This should lead to better 2045 2045 !-- load balance between boundary and non-boundary PEs. 2046 2046 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & … … 2052 2052 nzb_max_l = nzb_max 2053 2053 END IF 2054 2054 2055 2055 gu = 2.0_wp * u_gtrans 2056 2056 gv = 2.0_wp * v_gtrans … … 2058 2058 !-- Compute southside fluxes for the respective boundary of PE 2059 2059 IF ( j == nys ) THEN 2060 2060 2061 2061 DO k = nzb+1, nzb_max_l 2062 2062 … … 2104 2104 37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) & 2105 2105 - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) & 2106 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 2106 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 2107 2107 diss_s_u(k,tn) = - ABS(v_comp(k)) * ( & 2108 2108 10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) & … … 2111 2111 2112 2112 ENDDO 2113 2113 2114 2114 ENDIF 2115 2115 ! 2116 2116 !-- Compute leftside fluxes for the respective boundary of PE 2117 2117 IF ( i == i_omp .OR. i == nxlu ) THEN 2118 2118 2119 2119 DO k = nzb+1, nzb_max_l 2120 2120 … … 2137 2137 ) * & 2138 2138 ( u(k,j,i+2) + u(k,j,i-3) ) & 2139 ) 2140 2139 ) 2140 2141 2141 diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( & 2142 2142 ( 10.0_wp * ibit2 * adv_mom_5 & … … 2169 2169 2170 2170 ENDDO 2171 2171 2172 2172 ENDIF 2173 2173 ! … … 2188 2188 ( u(k,j,i+1) + u(k,j,i) ) & 2189 2189 - ( 8.0_wp * ibit2 * adv_mom_5 & 2190 + ibit1 * adv_mom_3 & 2190 + ibit1 * adv_mom_3 & 2191 2191 ) * & 2192 2192 ( u(k,j,i+2) + u(k,j,i-1) ) & … … 2194 2194 ) * & 2195 2195 ( u(k,j,i+3) + u(k,j,i-2) ) & 2196 ) 2197 2196 ) 2197 2198 2198 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2199 2199 ( 10.0_wp * ibit2 * adv_mom_5 & … … 2229 2229 ) * & 2230 2230 ( u(k,j+3,i) + u(k,j-2,i) ) & 2231 ) 2232 2231 ) 2232 2233 2233 diss_n(k) = - ABS ( v_comp(k) ) * ( & 2234 2234 ( 10.0_wp * ibit5 * adv_mom_5 & … … 2253 2253 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) & 2254 2254 - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) & 2255 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 2255 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 2256 2256 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2257 2257 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) & 2258 2258 - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) & 2259 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2260 2261 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 2259 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2260 2261 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 2262 2262 flux_n(k) = v_comp(k) * ( & 2263 2263 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & 2264 2264 - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & 2265 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2265 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2266 2266 diss_n(k) = - ABS( v_comp(k) ) * ( & 2267 2267 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & … … 2271 2271 ENDDO 2272 2272 ! 2273 !-- Now, compute vertical fluxes. Split loop into a part treating the 2273 !-- Now, compute vertical fluxes. Split loop into a part treating the 2274 2274 !-- lowest grid points with indirect indexing, a main loop without 2275 2275 !-- indirect indexing, and a loop for the uppermost grip points with 2276 2276 !-- indirect indexing. This allows better vectorization for the main loop. 2277 !-- First, compute the flux at model surface, which need has to be 2277 !-- First, compute the flux at model surface, which need has to be 2278 2278 !-- calculated explicitly for the tendency at 2279 2279 !-- the first w-level. For topography wall this is done implicitely by … … 2282 2282 diss_t(nzb) = 0.0_wp 2283 2283 w_comp(nzb) = 0.0_wp 2284 2284 2285 2285 DO k = nzb+1, nzb+1 2286 2286 ! … … 2309 2309 ) * & 2310 2310 ( u(k_ppp,j,i) + u(k_mm,j,i) ) & 2311 ) 2312 2311 ) 2312 2313 2313 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2314 2314 ( 10.0_wp * ibit8 * adv_mom_5 & … … 2326 2326 ) 2327 2327 ENDDO 2328 2328 2329 2329 DO k = nzb+2, nzt-2 2330 2330 … … 2344 2344 ) * & 2345 2345 ( u(k+2,j,i) + u(k-1,j,i) ) & 2346 + ( ibit8 * adv_mom_5 & 2346 + ( ibit8 * adv_mom_5 & 2347 2347 ) * & 2348 2348 ( u(k+3,j,i) + u(k-2,j,i) ) & … … 2364 2364 ) 2365 2365 ENDDO 2366 2366 2367 2367 DO k = nzt-1, nzt-symmetry_flag 2368 2368 ! … … 2408 2408 ) 2409 2409 ENDDO 2410 2410 2411 2411 ! 2412 2412 !-- Set resolved/turbulent flux at model top to zero (w-level). In case that 2413 2413 !-- a symmetric behavior between bottom and top shall be guaranteed (closed 2414 !-- channel flow), the flux at nzt is also set to zero. 2414 !-- channel flow), the flux at nzt is also set to zero. 2415 2415 IF ( symmetry_flag == 1 ) THEN 2416 2416 flux_t(nzt) = 0.0_wp … … 2421 2421 diss_t(nzt+1) = 0.0_wp 2422 2422 w_comp(nzt+1) = 0.0_wp 2423 2423 2424 2424 DO k = nzb+1, nzb_max_l 2425 2425 2426 2426 flux_d = flux_t(k-1) 2427 2427 diss_d = diss_t(k-1) 2428 2428 2429 2429 ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp ) 2430 2430 ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp ) 2431 2431 ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp ) 2432 2432 2433 2433 ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp ) 2434 2434 ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp ) 2435 2435 ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp ) 2436 2436 2437 2437 ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) 2438 2438 ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) … … 2464 2464 + REAL( IBITS(advc_flags_m(k-1,j,i),7,1), KIND = wp ) & 2465 2465 + REAL( IBITS(advc_flags_m(k-1,j,i),8,1), KIND = wp ) & 2466 ) & 2466 ) & 2467 2467 ) * drho_air(k) * ddzw(k) & 2468 ) * 0.5_wp 2469 2468 ) * 0.5_wp 2469 2470 2470 tend(k,j,i) = tend(k,j,i) - ( & 2471 2471 ( flux_r(k) + diss_r(k) & … … 2504 2504 ) * weight_substep(intermediate_timestep_count) 2505 2505 ENDDO 2506 2506 2507 2507 DO k = nzb_max_l+1, nzt 2508 2508 … … 2516 2516 + ( v_comp(k) + gv - ( v(k,j,i) + v(k,j,i-1) ) ) * ddy & 2517 2517 + ( w_comp(k) * rho_air_zw(k) & 2518 - w_comp(k-1) * rho_air_zw(k-1) & 2518 - w_comp(k-1) * rho_air_zw(k-1) & 2519 2519 ) * drho_air(k) * ddzw(k) & 2520 2520 ) * 0.5_wp … … 2528 2528 - ( flux_d + diss_d ) & 2529 2529 ) * drho_air(k) * ddzw(k) & 2530 2530 ) + div * u(k,j,i) 2531 2531 2532 2532 flux_l_u(k,j,tn) = flux_r(k) … … 2578 2578 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 2579 2579 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 2580 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 2580 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 2581 2581 INTEGER(iwp) :: tn !< number of OpenMP thread 2582 2582 2583 2583 REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction 2584 2584 REAL(wp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction 2585 2585 REAL(wp) :: ibit11 !< flag indicating 5th-order scheme along x-direction 2586 REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction 2586 REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction 2587 2587 REAL(wp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction 2588 2588 REAL(wp) :: ibit14 !< flag indicating 3rd-order scheme along y-direction … … 2596 2596 REAL(wp) :: gv !< Galilei-transformation velocity along y 2597 2597 REAL(wp) :: v_comp_l !< advection velocity along y on leftmost grid point on subdomain 2598 2598 2599 2599 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 2600 2600 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box … … 2607 2607 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z 2608 2608 ! 2609 !-- Used local modified copy of nzb_max (used to degrade order of 2609 !-- Used local modified copy of nzb_max (used to degrade order of 2610 2610 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 2611 !-- instead of the entire subdomain. This should lead to better 2611 !-- instead of the entire subdomain. This should lead to better 2612 2612 !-- load balance between boundary and non-boundary PEs. 2613 2613 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & … … 2619 2619 nzb_max_l = nzb_max 2620 2620 END IF 2621 2621 2622 2622 gu = 2.0_wp * u_gtrans 2623 2623 gv = 2.0_wp * v_gtrans 2624 2624 2625 ! 2625 ! 2626 2626 !-- Compute leftside fluxes for the respective boundary. 2627 2627 IF ( i == i_omp ) THEN … … 2634 2634 2635 2635 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu 2636 flux_l_v(k,j,tn) = u_comp(k) * ( &2637 ( 37.0_wp * ibit11 * adv_mom_5 &2638 + 7.0_wp * ibit10 * adv_mom_3 &2639 + ibit9 * adv_mom_1 &2640 ) * &2641 ( v(k,j,i) + v(k,j,i-1) ) &2642 - ( 8.0_wp * ibit11 * adv_mom_5 &2643 + ibit10 * adv_mom_3 &2644 ) * &2645 ( v(k,j,i+1) + v(k,j,i-2) ) &2646 + ( ibit11 * adv_mom_5 &2647 ) * &2648 ( v(k,j,i+2) + v(k,j,i-3) ) &2636 flux_l_v(k,j,tn) = u_comp(k) * ( & 2637 ( 37.0_wp * ibit11 * adv_mom_5 & 2638 + 7.0_wp * ibit10 * adv_mom_3 & 2639 + ibit9 * adv_mom_1 & 2640 ) * & 2641 ( v(k,j,i) + v(k,j,i-1) ) & 2642 - ( 8.0_wp * ibit11 * adv_mom_5 & 2643 + ibit10 * adv_mom_3 & 2644 ) * & 2645 ( v(k,j,i+1) + v(k,j,i-2) ) & 2646 + ( ibit11 * adv_mom_5 & 2647 ) * & 2648 ( v(k,j,i+2) + v(k,j,i-3) ) & 2649 2649 ) 2650 2650 2651 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( &2652 ( 10.0_wp * ibit11 * adv_mom_5 &2653 + 3.0_wp * ibit10 * adv_mom_3 &2654 + ibit9 * adv_mom_1 &2655 ) * &2656 ( v(k,j,i) - v(k,j,i-1) ) &2657 - ( 5.0_wp * ibit11 * adv_mom_5 &2658 + ibit10 * adv_mom_3 &2659 ) * &2660 ( v(k,j,i+1) - v(k,j,i-2) ) &2661 + ( ibit11 * adv_mom_5 &2662 ) * &2663 ( v(k,j,i+2) - v(k,j,i-3) ) &2651 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & 2652 ( 10.0_wp * ibit11 * adv_mom_5 & 2653 + 3.0_wp * ibit10 * adv_mom_3 & 2654 + ibit9 * adv_mom_1 & 2655 ) * & 2656 ( v(k,j,i) - v(k,j,i-1) ) & 2657 - ( 5.0_wp * ibit11 * adv_mom_5 & 2658 + ibit10 * adv_mom_3 & 2659 ) * & 2660 ( v(k,j,i+1) - v(k,j,i-2) ) & 2661 + ( ibit11 * adv_mom_5 & 2662 ) * & 2663 ( v(k,j,i+2) - v(k,j,i-3) ) & 2664 2664 ) 2665 2665 … … 2669 2669 2670 2670 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu 2671 flux_l_v(k,j,tn) = u_comp(k) * ( &2672 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) &2673 - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) &2671 flux_l_v(k,j,tn) = u_comp(k) * ( & 2672 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) & 2673 - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) & 2674 2674 + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 2675 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( &2676 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) &2677 - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) &2675 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & 2676 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) & 2677 - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) & 2678 2678 + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 2679 2679 2680 2680 ENDDO 2681 2681 2682 2682 ENDIF 2683 2683 ! 2684 2684 !-- Compute southside fluxes for the respective boundary. 2685 2685 IF ( j == nysv ) THEN 2686 2686 2687 2687 DO k = nzb+1, nzb_max_l 2688 2688 … … 2692 2692 2693 2693 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv 2694 flux_s_v(k,tn) = v_comp_l * ( &2695 ( 37.0_wp * ibit14 * adv_mom_5 &2696 + 7.0_wp * ibit13 * adv_mom_3 &2697 + ibit12 * adv_mom_1 &2698 ) * &2699 ( v(k,j,i) + v(k,j-1,i) ) &2700 - ( 8.0_wp * ibit14 * adv_mom_5 &2701 + ibit13 * adv_mom_3 &2702 ) * &2703 ( v(k,j+1,i) + v(k,j-2,i) ) &2704 + ( ibit14 * adv_mom_5 &2705 ) * &2706 ( v(k,j+2,i) + v(k,j-3,i) ) &2694 flux_s_v(k,tn) = v_comp_l * ( & 2695 ( 37.0_wp * ibit14 * adv_mom_5 & 2696 + 7.0_wp * ibit13 * adv_mom_3 & 2697 + ibit12 * adv_mom_1 & 2698 ) * & 2699 ( v(k,j,i) + v(k,j-1,i) ) & 2700 - ( 8.0_wp * ibit14 * adv_mom_5 & 2701 + ibit13 * adv_mom_3 & 2702 ) * & 2703 ( v(k,j+1,i) + v(k,j-2,i) ) & 2704 + ( ibit14 * adv_mom_5 & 2705 ) * & 2706 ( v(k,j+2,i) + v(k,j-3,i) ) & 2707 2707 ) 2708 2708 2709 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &2710 ( 10.0_wp * ibit14 * adv_mom_5 &2711 + 3.0_wp * ibit13 * adv_mom_3 &2712 + ibit12 * adv_mom_1 &2713 ) * &2714 ( v(k,j,i) - v(k,j-1,i) ) &2715 - ( 5.0_wp * ibit14 * adv_mom_5 &2716 + ibit13 * adv_mom_3 &2717 ) * &2718 ( v(k,j+1,i) - v(k,j-2,i) ) &2719 + ( ibit14 * adv_mom_5 &2720 ) * &2721 ( v(k,j+2,i) - v(k,j-3,i) ) &2709 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & 2710 ( 10.0_wp * ibit14 * adv_mom_5 & 2711 + 3.0_wp * ibit13 * adv_mom_3 & 2712 + ibit12 * adv_mom_1 & 2713 ) * & 2714 ( v(k,j,i) - v(k,j-1,i) ) & 2715 - ( 5.0_wp * ibit14 * adv_mom_5 & 2716 + ibit13 * adv_mom_3 & 2717 ) * & 2718 ( v(k,j+1,i) - v(k,j-2,i) ) & 2719 + ( ibit14 * adv_mom_5 & 2720 ) * & 2721 ( v(k,j+2,i) - v(k,j-3,i) ) & 2722 2722 ) 2723 2723 … … 2727 2727 2728 2728 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv 2729 flux_s_v(k,tn) = v_comp_l * ( &2730 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) &2731 - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) &2729 flux_s_v(k,tn) = v_comp_l * ( & 2730 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) & 2731 - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) & 2732 2732 + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 2733 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &2734 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) &2735 - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) &2733 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & 2734 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) & 2735 - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) & 2736 2736 + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 2737 2737 2738 2738 ENDDO 2739 2739 2740 2740 ENDIF 2741 2741 ! … … 2747 2747 ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp ) 2748 2748 ibit9 = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp ) 2749 2749 2750 2750 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu 2751 flux_r(k) = u_comp(k) * ( &2752 ( 37.0_wp * ibit11 * adv_mom_5 &2753 + 7.0_wp * ibit10 * adv_mom_3 &2754 + ibit9 * adv_mom_1 &2755 ) * &2756 ( v(k,j,i+1) + v(k,j,i) ) &2757 - ( 8.0_wp * ibit11 * adv_mom_5 &2758 + ibit10 * adv_mom_3 &2759 ) * &2760 ( v(k,j,i+2) + v(k,j,i-1) ) &2761 + ( ibit11 * adv_mom_5 &2762 ) * &2763 ( v(k,j,i+3) + v(k,j,i-2) ) &2751 flux_r(k) = u_comp(k) * ( & 2752 ( 37.0_wp * ibit11 * adv_mom_5 & 2753 + 7.0_wp * ibit10 * adv_mom_3 & 2754 + ibit9 * adv_mom_1 & 2755 ) * & 2756 ( v(k,j,i+1) + v(k,j,i) ) & 2757 - ( 8.0_wp * ibit11 * adv_mom_5 & 2758 + ibit10 * adv_mom_3 & 2759 ) * & 2760 ( v(k,j,i+2) + v(k,j,i-1) ) & 2761 + ( ibit11 * adv_mom_5 & 2762 ) * & 2763 ( v(k,j,i+3) + v(k,j,i-2) ) & 2764 2764 ) 2765 2765 2766 diss_r(k) = - ABS( u_comp(k) ) * ( &2767 ( 10.0_wp * ibit11 * adv_mom_5 &2768 + 3.0_wp * ibit10 * adv_mom_3 &2769 + ibit9 * adv_mom_1 &2770 ) * &2771 ( v(k,j,i+1) - v(k,j,i) ) &2772 - ( 5.0_wp * ibit11 * adv_mom_5 &2773 + ibit10 * adv_mom_3 &2774 ) * &2775 ( v(k,j,i+2) - v(k,j,i-1) ) &2776 + ( ibit11 * adv_mom_5 &2777 ) * &2778 ( v(k,j,i+3) - v(k,j,i-2) ) &2766 diss_r(k) = - ABS( u_comp(k) ) * ( & 2767 ( 10.0_wp * ibit11 * adv_mom_5 & 2768 + 3.0_wp * ibit10 * adv_mom_3 & 2769 + ibit9 * adv_mom_1 & 2770 ) * & 2771 ( v(k,j,i+1) - v(k,j,i) ) & 2772 - ( 5.0_wp * ibit11 * adv_mom_5 & 2773 + ibit10 * adv_mom_3 & 2774 ) * & 2775 ( v(k,j,i+2) - v(k,j,i-1) ) & 2776 + ( ibit11 * adv_mom_5 & 2777 ) * & 2778 ( v(k,j,i+3) - v(k,j,i-2) ) & 2779 2779 ) 2780 2780 … … 2785 2785 2786 2786 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2787 flux_n(k) = ( v_comp(k) - gv ) * ( &2788 ( 37.0_wp * ibit14 * adv_mom_5 &2789 + 7.0_wp * ibit13 * adv_mom_3 &2790 + ibit12 * adv_mom_1 &2791 ) * &2792 ( v(k,j+1,i) + v(k,j,i) ) &2793 - ( 8.0_wp * ibit14 * adv_mom_5 &2794 + ibit13 * adv_mom_3 &2795 ) * &2796 ( v(k,j+2,i) + v(k,j-1,i) ) &2797 + ( ibit14 * adv_mom_5 &2798 ) * &2799 ( v(k,j+3,i) + v(k,j-2,i) ) &2787 flux_n(k) = ( v_comp(k) - gv ) * ( & 2788 ( 37.0_wp * ibit14 * adv_mom_5 & 2789 + 7.0_wp * ibit13 * adv_mom_3 & 2790 + ibit12 * adv_mom_1 & 2791 ) * & 2792 ( v(k,j+1,i) + v(k,j,i) ) & 2793 - ( 8.0_wp * ibit14 * adv_mom_5 & 2794 + ibit13 * adv_mom_3 & 2795 ) * & 2796 ( v(k,j+2,i) + v(k,j-1,i) ) & 2797 + ( ibit14 * adv_mom_5 & 2798 ) * & 2799 ( v(k,j+3,i) + v(k,j-2,i) ) & 2800 2800 ) 2801 2801 2802 diss_n(k) = - ABS( v_comp(k) - gv ) * ( &2803 ( 10.0_wp * ibit14 * adv_mom_5 &2804 + 3.0_wp * ibit13 * adv_mom_3 &2805 + ibit12 * adv_mom_1 &2806 ) * &2807 ( v(k,j+1,i) - v(k,j,i) ) &2808 - ( 5.0_wp * ibit14 * adv_mom_5 &2809 + ibit13 * adv_mom_3 &2810 ) * &2811 ( v(k,j+2,i) - v(k,j-1,i) ) &2812 + ( ibit14 * adv_mom_5 &2813 ) * &2814 ( v(k,j+3,i) - v(k,j-2,i) ) &2802 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 2803 ( 10.0_wp * ibit14 * adv_mom_5 & 2804 + 3.0_wp * ibit13 * adv_mom_3 & 2805 + ibit12 * adv_mom_1 & 2806 ) * & 2807 ( v(k,j+1,i) - v(k,j,i) ) & 2808 - ( 5.0_wp * ibit14 * adv_mom_5 & 2809 + ibit13 * adv_mom_3 & 2810 ) * & 2811 ( v(k,j+2,i) - v(k,j-1,i) ) & 2812 + ( ibit14 * adv_mom_5 & 2813 ) * & 2814 ( v(k,j+3,i) - v(k,j-2,i) ) & 2815 2815 ) 2816 2816 ENDDO … … 2819 2819 2820 2820 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu 2821 flux_r(k) = u_comp(k) * ( &2822 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &2823 - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &2821 flux_r(k) = u_comp(k) * ( & 2822 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) & 2823 - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) & 2824 2824 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 2825 2825 2826 diss_r(k) = - ABS( u_comp(k) ) * ( &2827 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &2828 - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &2826 diss_r(k) = - ABS( u_comp(k) ) * ( & 2827 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) & 2828 - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) & 2829 2829 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 2830 2830 2831 2831 2832 2832 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2833 flux_n(k) = ( v_comp(k) - gv ) * ( &2834 37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &2835 - 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &2833 flux_n(k) = ( v_comp(k) - gv ) * ( & 2834 37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) & 2835 - 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) & 2836 2836 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 2837 2837 2838 diss_n(k) = - ABS( v_comp(k) - gv ) * ( &2839 10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &2840 - 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &2838 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 2839 10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) & 2840 - 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) & 2841 2841 + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 2842 2842 ENDDO 2843 2843 ! 2844 !-- Now, compute vertical fluxes. Split loop into a part treating the 2844 !-- Now, compute vertical fluxes. Split loop into a part treating the 2845 2845 !-- lowest grid points with indirect indexing, a main loop without 2846 2846 !-- indirect indexing, and a loop for the uppermost grip points with 2847 2847 !-- indirect indexing. This allows better vectorization for the main loop. 2848 !-- First, compute the flux at model surface, which need has to be 2848 !-- First, compute the flux at model surface, which need has to be 2849 2849 !-- calculated explicitly for the tendency at 2850 2850 !-- the first w-level. For topography wall this is done implicitely by … … 2853 2853 diss_t(nzb) = 0.0_wp 2854 2854 w_comp(nzb) = 0.0_wp 2855 2855 2856 2856 DO k = nzb+1, nzb+1 2857 2857 ! … … 2867 2867 2868 2868 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2869 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2870 ( 37.0_wp * ibit17 * adv_mom_5 &2871 + 7.0_wp * ibit16 * adv_mom_3 &2872 + ibit15 * adv_mom_1 &2873 ) * &2874 ( v(k+1,j,i) + v(k,j,i) ) &2875 - ( 8.0_wp * ibit17 * adv_mom_5 &2876 + ibit16 * adv_mom_3 &2877 ) * &2878 ( v(k_pp,j,i) + v(k-1,j,i) ) &2879 + ( ibit17 * adv_mom_5 &2880 ) * &2881 ( v(k_ppp,j,i) + v(k_mm,j,i) ) &2869 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2870 ( 37.0_wp * ibit17 * adv_mom_5 & 2871 + 7.0_wp * ibit16 * adv_mom_3 & 2872 + ibit15 * adv_mom_1 & 2873 ) * & 2874 ( v(k+1,j,i) + v(k,j,i) ) & 2875 - ( 8.0_wp * ibit17 * adv_mom_5 & 2876 + ibit16 * adv_mom_3 & 2877 ) * & 2878 ( v(k_pp,j,i) + v(k-1,j,i) ) & 2879 + ( ibit17 * adv_mom_5 & 2880 ) * & 2881 ( v(k_ppp,j,i) + v(k_mm,j,i) ) & 2882 2882 ) 2883 2883 2884 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2885 ( 10.0_wp * ibit17 * adv_mom_5 &2886 + 3.0_wp * ibit16 * adv_mom_3 &2887 + ibit15 * adv_mom_1 &2888 ) * &2889 ( v(k+1,j,i) - v(k,j,i) ) &2890 - ( 5.0_wp * ibit17 * adv_mom_5 &2891 + ibit16 * adv_mom_3 &2892 ) * &2893 ( v(k_pp,j,i) - v(k-1,j,i) ) &2894 + ( ibit17 * adv_mom_5 &2895 ) * &2896 ( v(k_ppp,j,i) - v(k_mm,j,i) ) &2884 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2885 ( 10.0_wp * ibit17 * adv_mom_5 & 2886 + 3.0_wp * ibit16 * adv_mom_3 & 2887 + ibit15 * adv_mom_1 & 2888 ) * & 2889 ( v(k+1,j,i) - v(k,j,i) ) & 2890 - ( 5.0_wp * ibit17 * adv_mom_5 & 2891 + ibit16 * adv_mom_3 & 2892 ) * & 2893 ( v(k_pp,j,i) - v(k-1,j,i) ) & 2894 + ( ibit17 * adv_mom_5 & 2895 ) * & 2896 ( v(k_ppp,j,i) - v(k_mm,j,i) ) & 2897 2897 ) 2898 2898 ENDDO 2899 2899 2900 2900 DO k = nzb+2, nzt-2 2901 2901 … … 2905 2905 2906 2906 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2907 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2908 ( 37.0_wp * ibit17 * adv_mom_5 &2909 + 7.0_wp * ibit16 * adv_mom_3 &2910 + ibit15 * adv_mom_1 &2911 ) * &2912 ( v(k+1,j,i) + v(k,j,i) ) &2913 - ( 8.0_wp * ibit17 * adv_mom_5 &2914 + ibit16 * adv_mom_3 &2915 ) * &2916 ( v(k+2,j,i) + v(k-1,j,i) ) &2917 + ( ibit17 * adv_mom_5 &2918 ) * &2919 ( v(k+3,j,i) + v(k-2,j,i) ) &2907 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2908 ( 37.0_wp * ibit17 * adv_mom_5 & 2909 + 7.0_wp * ibit16 * adv_mom_3 & 2910 + ibit15 * adv_mom_1 & 2911 ) * & 2912 ( v(k+1,j,i) + v(k,j,i) ) & 2913 - ( 8.0_wp * ibit17 * adv_mom_5 & 2914 + ibit16 * adv_mom_3 & 2915 ) * & 2916 ( v(k+2,j,i) + v(k-1,j,i) ) & 2917 + ( ibit17 * adv_mom_5 & 2918 ) * & 2919 ( v(k+3,j,i) + v(k-2,j,i) ) & 2920 2920 ) 2921 2921 2922 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2923 ( 10.0_wp * ibit17 * adv_mom_5 &2924 + 3.0_wp * ibit16 * adv_mom_3 &2925 + ibit15 * adv_mom_1 &2926 ) * &2927 ( v(k+1,j,i) - v(k,j,i) ) &2928 - ( 5.0_wp * ibit17 * adv_mom_5 &2929 + ibit16 * adv_mom_3 &2930 ) * &2922 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2923 ( 10.0_wp * ibit17 * adv_mom_5 & 2924 + 3.0_wp * ibit16 * adv_mom_3 & 2925 + ibit15 * adv_mom_1 & 2926 ) * & 2927 ( v(k+1,j,i) - v(k,j,i) ) & 2928 - ( 5.0_wp * ibit17 * adv_mom_5 & 2929 + ibit16 * adv_mom_3 & 2930 ) * & 2931 2931 ( v(k+2,j,i) - v(k-1,j,i) ) & 2932 + ( ibit17 * adv_mom_5 &2933 ) * &2934 ( v(k+3,j,i) - v(k-2,j,i) ) &2932 + ( ibit17 * adv_mom_5 & 2933 ) * & 2934 ( v(k+3,j,i) - v(k-2,j,i) ) & 2935 2935 ) 2936 2936 ENDDO 2937 2937 2938 2938 DO k = nzt-1, nzt-symmetry_flag 2939 2939 ! … … 2949 2949 2950 2950 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2951 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2952 ( 37.0_wp * ibit17 * adv_mom_5 &2953 + 7.0_wp * ibit16 * adv_mom_3 &2954 + ibit15 * adv_mom_1 &2955 ) * &2956 ( v(k+1,j,i) + v(k,j,i) ) &2957 - ( 8.0_wp * ibit17 * adv_mom_5 &2958 + ibit16 * adv_mom_3 &2959 ) * &2960 ( v(k_pp,j,i) + v(k-1,j,i) ) &2961 + ( ibit17 * adv_mom_5 &2962 ) * &2963 ( v(k_ppp,j,i) + v(k_mm,j,i) ) &2951 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2952 ( 37.0_wp * ibit17 * adv_mom_5 & 2953 + 7.0_wp * ibit16 * adv_mom_3 & 2954 + ibit15 * adv_mom_1 & 2955 ) * & 2956 ( v(k+1,j,i) + v(k,j,i) ) & 2957 - ( 8.0_wp * ibit17 * adv_mom_5 & 2958 + ibit16 * adv_mom_3 & 2959 ) * & 2960 ( v(k_pp,j,i) + v(k-1,j,i) ) & 2961 + ( ibit17 * adv_mom_5 & 2962 ) * & 2963 ( v(k_ppp,j,i) + v(k_mm,j,i) ) & 2964 2964 ) 2965 2965 2966 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2967 ( 10.0_wp * ibit17 * adv_mom_5 &2968 + 3.0_wp * ibit16 * adv_mom_3 &2969 + ibit15 * adv_mom_1 &2970 ) * &2971 ( v(k+1,j,i) - v(k,j,i) ) &2972 - ( 5.0_wp * ibit17 * adv_mom_5 &2973 + ibit16 * adv_mom_3 &2974 ) * &2975 ( v(k_pp,j,i) - v(k-1,j,i) ) &2976 + ( ibit17 * adv_mom_5 &2977 ) * &2978 ( v(k_ppp,j,i) - v(k_mm,j,i) ) &2966 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2967 ( 10.0_wp * ibit17 * adv_mom_5 & 2968 + 3.0_wp * ibit16 * adv_mom_3 & 2969 + ibit15 * adv_mom_1 & 2970 ) * & 2971 ( v(k+1,j,i) - v(k,j,i) ) & 2972 - ( 5.0_wp * ibit17 * adv_mom_5 & 2973 + ibit16 * adv_mom_3 & 2974 ) * & 2975 ( v(k_pp,j,i) - v(k-1,j,i) ) & 2976 + ( ibit17 * adv_mom_5 & 2977 ) * & 2978 ( v(k_ppp,j,i) - v(k_mm,j,i) ) & 2979 2979 ) 2980 2980 ENDDO 2981 2981 2982 2982 ! 2983 2983 !-- Set resolved/turbulent flux at model top to zero (w-level). In case that 2984 2984 !-- a symmetric behavior between bottom and top shall be guaranteed (closed 2985 !-- channel flow), the flux at nzt is also set to zero. 2985 !-- channel flow), the flux at nzt is also set to zero. 2986 2986 IF ( symmetry_flag == 1 ) THEN 2987 2987 flux_t(nzt) = 0.0_wp … … 2992 2992 diss_t(nzt+1) = 0.0_wp 2993 2993 w_comp(nzt+1) = 0.0_wp 2994 2994 2995 2995 DO k = nzb+1, nzb_max_l 2996 2996 2997 2997 flux_d = flux_t(k-1) 2998 2998 diss_d = diss_t(k-1) 2999 2999 3000 3000 ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp ) 3001 3001 ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp ) 3002 3002 ibit9 = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp ) 3003 3003 3004 3004 ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp ) 3005 3005 ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp ) 3006 3006 ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp ) 3007 3007 3008 3008 ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp ) 3009 3009 ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp ) … … 3013 3013 !-- correction is needed to overcome numerical instabilities introduced 3014 3014 !-- by a not sufficient reduction of divergences near topography. 3015 div = ( ( ( u_comp(k) + gu ) &3016 * ( ibit9 + ibit10 + ibit11 ) &3017 - ( u(k,j-1,i) + u(k,j,i) ) &3018 * ( &3019 REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp ) &3020 + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp ) &3021 + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp ) &3022 ) &3023 ) * ddx &3024 + ( v_comp(k) &3025 * ( ibit12 + ibit13 + ibit14 ) &3026 - ( v(k,j,i) + v(k,j-1,i) ) &3027 * ( &3028 REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp ) &3029 + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp ) &3030 + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp ) &3031 ) &3032 ) * ddy &3033 + ( w_comp(k) * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) &3034 - w_comp(k-1) * rho_air_zw(k-1) &3035 * ( &3036 REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp ) &3037 + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp ) &3038 + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp ) &3039 ) &3040 ) * drho_air(k) * ddzw(k) &3015 div = ( ( ( u_comp(k) + gu ) & 3016 * ( ibit9 + ibit10 + ibit11 ) & 3017 - ( u(k,j-1,i) + u(k,j,i) ) & 3018 * ( & 3019 REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp ) & 3020 + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp ) & 3021 + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp ) & 3022 ) & 3023 ) * ddx & 3024 + ( v_comp(k) & 3025 * ( ibit12 + ibit13 + ibit14 ) & 3026 - ( v(k,j,i) + v(k,j-1,i) ) & 3027 * ( & 3028 REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp ) & 3029 + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp ) & 3030 + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp ) & 3031 ) & 3032 ) * ddy & 3033 + ( w_comp(k) * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) & 3034 - w_comp(k-1) * rho_air_zw(k-1) & 3035 * ( & 3036 REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp ) & 3037 + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp ) & 3038 + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp ) & 3039 ) & 3040 ) * drho_air(k) * ddzw(k) & 3041 3041 ) * 0.5_wp 3042 3042 3043 tend(k,j,i) = tend(k,j,i) - ( &3044 ( flux_r(k) + diss_r(k) &3045 - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &3046 + ( flux_n(k) + diss_n(k) &3047 - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &3048 + ( ( flux_t(k) + diss_t(k) ) &3049 - ( flux_d + diss_d ) &3050 ) * drho_air(k) * ddzw(k) &3043 tend(k,j,i) = tend(k,j,i) - ( & 3044 ( flux_r(k) + diss_r(k) & 3045 - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx & 3046 + ( flux_n(k) + diss_n(k) & 3047 - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & 3048 + ( ( flux_t(k) + diss_t(k) ) & 3049 - ( flux_d + diss_d ) & 3050 ) * drho_air(k) * ddzw(k) & 3051 3051 ) + v(k,j,i) * div 3052 3052 … … 3078 3078 3079 3079 ENDDO 3080 3080 3081 3081 DO k = nzb_max_l+1, nzt 3082 3082 … … 3152 3152 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 3153 3153 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 3154 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 3154 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 3155 3155 INTEGER(iwp) :: tn !< number of OpenMP thread 3156 3156 3157 3157 REAL(wp) :: ibit18 !< flag indicating 1st-order scheme along x-direction 3158 3158 REAL(wp) :: ibit19 !< flag indicating 3rd-order scheme along x-direction … … 3169 3169 REAL(wp) :: gu !< Galilei-transformation velocity along x 3170 3170 REAL(wp) :: gv !< Galilei-transformation velocity along y 3171 3171 3172 3172 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 3173 3173 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box … … 3180 3180 REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z 3181 3181 ! 3182 !-- Used local modified copy of nzb_max (used to degrade order of 3182 !-- Used local modified copy of nzb_max (used to degrade order of 3183 3183 !-- discretization) at non-cyclic boundaries. Modify only at relevant points 3184 !-- instead of the entire subdomain. This should lead to better 3184 !-- instead of the entire subdomain. This should lead to better 3185 3185 !-- load balance between boundary and non-boundary PEs. 3186 3186 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & … … 3192 3192 nzb_max_l = nzb_max 3193 3193 END IF 3194 3194 3195 3195 gu = 2.0_wp * u_gtrans 3196 3196 gv = 2.0_wp * v_gtrans … … 3205 3205 3206 3206 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv 3207 flux_s_w(k,tn) = v_comp(k) * ( &3208 ( 37.0_wp * ibit23 * adv_mom_5 &3209 + 7.0_wp * ibit22 * adv_mom_3 &3210 + ibit21 * adv_mom_1 &3211 ) * &3212 ( w(k,j,i) + w(k,j-1,i) ) &3213 - ( 8.0_wp * ibit23 * adv_mom_5 &3214 + ibit22 * adv_mom_3 &3215 ) * &3216 ( w(k,j+1,i) + w(k,j-2,i) ) &3217 + ( ibit23 * adv_mom_5 &3218 ) * &3219 ( w(k,j+2,i) + w(k,j-3,i) ) &3207 flux_s_w(k,tn) = v_comp(k) * ( & 3208 ( 37.0_wp * ibit23 * adv_mom_5 & 3209 + 7.0_wp * ibit22 * adv_mom_3 & 3210 + ibit21 * adv_mom_1 & 3211 ) * & 3212 ( w(k,j,i) + w(k,j-1,i) ) & 3213 - ( 8.0_wp * ibit23 * adv_mom_5 & 3214 + ibit22 * adv_mom_3 & 3215 ) * & 3216 ( w(k,j+1,i) + w(k,j-2,i) ) & 3217 + ( ibit23 * adv_mom_5 & 3218 ) * & 3219 ( w(k,j+2,i) + w(k,j-3,i) ) & 3220 3220 ) 3221 3221 3222 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( &3223 ( 10.0_wp * ibit23 * adv_mom_5 &3224 + 3.0_wp * ibit22 * adv_mom_3 &3225 + ibit21 * adv_mom_1 &3226 ) * &3227 ( w(k,j,i) - w(k,j-1,i) ) &3228 - ( 5.0_wp * ibit23 * adv_mom_5 &3229 + ibit22 * adv_mom_3 &3230 ) * &3231 ( w(k,j+1,i) - w(k,j-2,i) ) &3232 + ( ibit23 * adv_mom_5 &3233 ) * &3234 ( w(k,j+2,i) - w(k,j-3,i) ) &3222 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( & 3223 ( 10.0_wp * ibit23 * adv_mom_5 & 3224 + 3.0_wp * ibit22 * adv_mom_3 & 3225 + ibit21 * adv_mom_1 & 3226 ) * & 3227 ( w(k,j,i) - w(k,j-1,i) ) & 3228 - ( 5.0_wp * ibit23 * adv_mom_5 & 3229 + ibit22 * adv_mom_3 & 3230 ) * & 3231 ( w(k,j+1,i) - w(k,j-2,i) ) & 3232 + ( ibit23 * adv_mom_5 & 3233 ) * & 3234 ( w(k,j+2,i) - w(k,j-3,i) ) & 3235 3235 ) 3236 3236 … … 3240 3240 3241 3241 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv 3242 flux_s_w(k,tn) = v_comp(k) * ( &3243 37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) &3244 - 8.0_wp * ( w(k,j+1,i) + w(k,j-2,i) ) &3242 flux_s_w(k,tn) = v_comp(k) * ( & 3243 37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) & 3244 - 8.0_wp * ( w(k,j+1,i) + w(k,j-2,i) ) & 3245 3245 + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 3246 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( &3247 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) &3248 - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) &3246 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( & 3247 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) & 3248 - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) & 3249 3249 + ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5 3250 3250 … … 3263 3263 3264 3264 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu 3265 flux_l_w(k,j,tn) = u_comp(k) * ( &3266 ( 37.0_wp * ibit20 * adv_mom_5 &3267 + 7.0_wp * ibit19 * adv_mom_3 &3268 + ibit18 * adv_mom_1 &3269 ) * &3270 ( w(k,j,i) + w(k,j,i-1) ) &3271 - ( 8.0_wp * ibit20 * adv_mom_5 &3272 + ibit19 * adv_mom_3 &3273 ) * &3274 ( w(k,j,i+1) + w(k,j,i-2) ) &3275 + ( ibit20 * adv_mom_5 &3276 ) * &3277 ( w(k,j,i+2) + w(k,j,i-3) ) &3265 flux_l_w(k,j,tn) = u_comp(k) * ( & 3266 ( 37.0_wp * ibit20 * adv_mom_5 & 3267 + 7.0_wp * ibit19 * adv_mom_3 & 3268 + ibit18 * adv_mom_1 & 3269 ) * & 3270 ( w(k,j,i) + w(k,j,i-1) ) & 3271 - ( 8.0_wp * ibit20 * adv_mom_5 & 3272 + ibit19 * adv_mom_3 & 3273 ) * & 3274 ( w(k,j,i+1) + w(k,j,i-2) ) & 3275 + ( ibit20 * adv_mom_5 & 3276 ) * & 3277 ( w(k,j,i+2) + w(k,j,i-3) ) & 3278 3278 ) 3279 3279 3280 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( &3281 ( 10.0_wp * ibit20 * adv_mom_5 &3282 + 3.0_wp * ibit19 * adv_mom_3 &3283 + ibit18 * adv_mom_1 &3284 ) * &3285 ( w(k,j,i) - w(k,j,i-1) ) &3286 - ( 5.0_wp * ibit20 * adv_mom_5 &3287 + ibit19 * adv_mom_3 &3288 ) * &3289 ( w(k,j,i+1) - w(k,j,i-2) ) &3290 + ( ibit20 * adv_mom_5 &3291 ) * &3292 ( w(k,j,i+2) - w(k,j,i-3) ) &3280 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( & 3281 ( 10.0_wp * ibit20 * adv_mom_5 & 3282 + 3.0_wp * ibit19 * adv_mom_3 & 3283 + ibit18 * adv_mom_1 & 3284 ) * & 3285 ( w(k,j,i) - w(k,j,i-1) ) & 3286 - ( 5.0_wp * ibit20 * adv_mom_5 & 3287 + ibit19 * adv_mom_3 & 3288 ) * & 3289 ( w(k,j,i+1) - w(k,j,i-2) ) & 3290 + ( ibit20 * adv_mom_5 & 3291 ) * & 3292 ( w(k,j,i+2) - w(k,j,i-3) ) & 3293 3293 ) 3294 3294 … … 3298 3298 3299 3299 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu 3300 flux_l_w(k,j,tn) = u_comp(k) * ( &3301 37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) &3302 - 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) &3300 flux_l_w(k,j,tn) = u_comp(k) * ( & 3301 37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) & 3302 - 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) & 3303 3303 + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 3304 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( &3305 10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) &3306 - 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) &3307 + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 3304 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( & 3305 10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) & 3306 - 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) & 3307 + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 3308 3308 3309 3309 ENDDO … … 3320 3320 3321 3321 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 3322 flux_r(k) = u_comp(k) * ( &3323 ( 37.0_wp * ibit20 * adv_mom_5 &3324 + 7.0_wp * ibit19 * adv_mom_3 &3325 + ibit18 * adv_mom_1 &3326 ) * &3327 ( w(k,j,i+1) + w(k,j,i) ) &3328 - ( 8.0_wp * ibit20 * adv_mom_5 &3329 + ibit19 * adv_mom_3 &3330 ) * &3331 ( w(k,j,i+2) + w(k,j,i-1) ) &3332 + ( ibit20 * adv_mom_5 &3333 ) * &3334 ( w(k,j,i+3) + w(k,j,i-2) ) &3322 flux_r(k) = u_comp(k) * ( & 3323 ( 37.0_wp * ibit20 * adv_mom_5 & 3324 + 7.0_wp * ibit19 * adv_mom_3 & 3325 + ibit18 * adv_mom_1 & 3326 ) * & 3327 ( w(k,j,i+1) + w(k,j,i) ) & 3328 - ( 8.0_wp * ibit20 * adv_mom_5 & 3329 + ibit19 * adv_mom_3 & 3330 ) * & 3331 ( w(k,j,i+2) + w(k,j,i-1) ) & 3332 + ( ibit20 * adv_mom_5 & 3333 ) * & 3334 ( w(k,j,i+3) + w(k,j,i-2) ) & 3335 3335 ) 3336 3336 3337 diss_r(k) = - ABS( u_comp(k) ) * ( &3338 ( 10.0_wp * ibit20 * adv_mom_5 &3339 + 3.0_wp * ibit19 * adv_mom_3 &3340 + ibit18 * adv_mom_1 &3341 ) * &3342 ( w(k,j,i+1) - w(k,j,i) ) &3343 - ( 5.0_wp * ibit20 * adv_mom_5 &3344 + ibit19 * adv_mom_3 &3345 ) * &3346 ( w(k,j,i+2) - w(k,j,i-1) ) &3347 + ( ibit20 * adv_mom_5 &3348 ) * &3349 ( w(k,j,i+3) - w(k,j,i-2) ) &3337 diss_r(k) = - ABS( u_comp(k) ) * ( & 3338 ( 10.0_wp * ibit20 * adv_mom_5 & 3339 + 3.0_wp * ibit19 * adv_mom_3 & 3340 + ibit18 * adv_mom_1 & 3341 ) * & 3342 ( w(k,j,i+1) - w(k,j,i) ) & 3343 - ( 5.0_wp * ibit20 * adv_mom_5 & 3344 + ibit19 * adv_mom_3 & 3345 ) * & 3346 ( w(k,j,i+2) - w(k,j,i-1) ) & 3347 + ( ibit20 * adv_mom_5 & 3348 ) * & 3349 ( w(k,j,i+3) - w(k,j,i-2) ) & 3350 3350 ) 3351 3351 … … 3355 3355 3356 3356 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 3357 flux_n(k) = v_comp(k) * ( &3358 ( 37.0_wp * ibit23 * adv_mom_5 &3359 + 7.0_wp * ibit22 * adv_mom_3 &3360 + ibit21 * adv_mom_1 &3361 ) * &3362 ( w(k,j+1,i) + w(k,j,i) ) &3363 - ( 8.0_wp * ibit23 * adv_mom_5 &3364 + ibit22 * adv_mom_3 &3365 ) * &3366 ( w(k,j+2,i) + w(k,j-1,i) ) &3367 + ( ibit23 * adv_mom_5 &3368 ) * &3369 ( w(k,j+3,i) + w(k,j-2,i) ) &3357 flux_n(k) = v_comp(k) * ( & 3358 ( 37.0_wp * ibit23 * adv_mom_5 & 3359 + 7.0_wp * ibit22 * adv_mom_3 & 3360 + ibit21 * adv_mom_1 & 3361 ) * & 3362 ( w(k,j+1,i) + w(k,j,i) ) & 3363 - ( 8.0_wp * ibit23 * adv_mom_5 & 3364 + ibit22 * adv_mom_3 & 3365 ) * & 3366 ( w(k,j+2,i) + w(k,j-1,i) ) & 3367 + ( ibit23 * adv_mom_5 & 3368 ) * & 3369 ( w(k,j+3,i) + w(k,j-2,i) ) & 3370 3370 ) 3371 3371 3372 diss_n(k) = - ABS( v_comp(k) ) * ( &3373 ( 10.0_wp * ibit23 * adv_mom_5 &3374 + 3.0_wp * ibit22 * adv_mom_3 &3375 + ibit21 * adv_mom_1 &3376 ) * &3377 ( w(k,j+1,i) - w(k,j,i) ) &3378 - ( 5.0_wp * ibit23 * adv_mom_5 &3379 + ibit22 * adv_mom_3 &3380 ) * &3381 ( w(k,j+2,i) - w(k,j-1,i) ) &3382 + ( ibit23 * adv_mom_5 &3383 ) * &3384 ( w(k,j+3,i) - w(k,j-2,i) ) &3372 diss_n(k) = - ABS( v_comp(k) ) * ( & 3373 ( 10.0_wp * ibit23 * adv_mom_5 & 3374 + 3.0_wp * ibit22 * adv_mom_3 & 3375 + ibit21 * adv_mom_1 & 3376 ) * & 3377 ( w(k,j+1,i) - w(k,j,i) ) & 3378 - ( 5.0_wp * ibit23 * adv_mom_5 & 3379 + ibit22 * adv_mom_3 & 3380 ) * & 3381 ( w(k,j+2,i) - w(k,j-1,i) ) & 3382 + ( ibit23 * adv_mom_5 & 3383 ) * & 3384 ( w(k,j+3,i) - w(k,j-2,i) ) & 3385 3385 ) 3386 3386 ENDDO … … 3389 3389 3390 3390 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 3391 flux_r(k) = u_comp(k) * ( &3392 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) &3393 - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) &3391 flux_r(k) = u_comp(k) * ( & 3392 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) & 3393 - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) & 3394 3394 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 3395 3395 3396 diss_r(k) = - ABS( u_comp(k) ) * ( &3397 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) &3398 - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) &3396 diss_r(k) = - ABS( u_comp(k) ) * ( & 3397 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) & 3398 - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) & 3399 3399 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 3400 3400 3401 3401 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 3402 flux_n(k) = v_comp(k) * ( &3403 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) &3404 - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) &3402 flux_n(k) = v_comp(k) * ( & 3403 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) & 3404 - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) & 3405 3405 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 3406 3406 3407 diss_n(k) = - ABS( v_comp(k) ) * ( &3408 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) &3409 - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) &3407 diss_n(k) = - ABS( v_comp(k) ) * ( & 3408 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) & 3409 - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) & 3410 3410 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 3411 3411 ENDDO 3412 3412 3413 3413 ! 3414 !-- Now, compute vertical fluxes. Split loop into a part treating the 3414 !-- Now, compute vertical fluxes. Split loop into a part treating the 3415 3415 !-- lowest grid points with indirect indexing, a main loop without 3416 3416 !-- indirect indexing, and a loop for the uppermost grip points with 3417 3417 !-- indirect indexing. This allows better vectorization for the main loop. 3418 !-- First, compute the flux at model surface, which need has to be 3418 !-- First, compute the flux at model surface, which need has to be 3419 3419 !-- calculated explicitly for the tendency at 3420 3420 !-- the first w-level. For topography wall this is done implicitely by 3421 !-- advc_flags_m. 3421 !-- advc_flags_m. First, compute flux at lowest level, located at z=dz/2. 3422 3422 k = nzb + 1 3423 3423 w_comp(k) = w(k,j,i) + w(k-1,j,i) … … 3426 3426 diss_t(0) = -ABS(w_comp(k)) * rho_air(k) & 3427 3427 * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 3428 3428 3429 3429 DO k = nzb+1, nzb+1 3430 3430 ! … … 3440 3440 3441 3441 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3442 flux_t(k) = w_comp(k) * rho_air(k+1) * ( &3443 ( 37.0_wp * ibit26 * adv_mom_5 &3444 + 7.0_wp * ibit25 * adv_mom_3 &3445 + ibit24 * adv_mom_1 &3446 ) * &3447 ( w(k+1,j,i) + w(k,j,i) ) &3448 - ( 8.0_wp * ibit26 * adv_mom_5 &3449 + ibit25 * adv_mom_3 &3450 ) * &3451 ( w(k_pp,j,i) + w(k-1,j,i) ) &3452 + ( ibit26 * adv_mom_5 &3453 ) * &3454 ( w(k_ppp,j,i) + w(k_mm,j,i) ) &3442 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 3443 ( 37.0_wp * ibit26 * adv_mom_5 & 3444 + 7.0_wp * ibit25 * adv_mom_3 & 3445 + ibit24 * adv_mom_1 & 3446 ) * & 3447 ( w(k+1,j,i) + w(k,j,i) ) & 3448 - ( 8.0_wp * ibit26 * adv_mom_5 & 3449 + ibit25 * adv_mom_3 & 3450 ) * & 3451 ( w(k_pp,j,i) + w(k-1,j,i) ) & 3452 + ( ibit26 * adv_mom_5 & 3453 ) * & 3454 ( w(k_ppp,j,i) + w(k_mm,j,i) ) & 3455 3455 ) 3456 3456 3457 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( &3458 ( 10.0_wp * ibit26 * adv_mom_5 &3459 + 3.0_wp * ibit25 * adv_mom_3 &3460 + ibit24 * adv_mom_1 &3461 ) * &3462 ( w(k+1,j,i) - w(k,j,i) ) &3463 - ( 5.0_wp * ibit26 * adv_mom_5 &3464 + ibit25 * adv_mom_3 &3465 ) * &3466 ( w(k_pp,j,i) - w(k-1,j,i) ) &3467 + ( ibit26 * adv_mom_5 &3468 ) * &3469 ( w(k_ppp,j,i) - w(k_mm,j,i) ) &3457 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 3458 ( 10.0_wp * ibit26 * adv_mom_5 & 3459 + 3.0_wp * ibit25 * adv_mom_3 & 3460 + ibit24 * adv_mom_1 & 3461 ) * & 3462 ( w(k+1,j,i) - w(k,j,i) ) & 3463 - ( 5.0_wp * ibit26 * adv_mom_5 & 3464 + ibit25 * adv_mom_3 & 3465 ) * & 3466 ( w(k_pp,j,i) - w(k-1,j,i) ) & 3467 + ( ibit26 * adv_mom_5 & 3468 ) * & 3469 ( w(k_ppp,j,i) - w(k_mm,j,i) ) & 3470 3470 ) 3471 3471 ENDDO 3472 3472 3473 3473 DO k = nzb+2, nzt-2 3474 3474 3475 3475 ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp ) 3476 3476 ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp ) … … 3478 3478 3479 3479 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3480 flux_t(k) = w_comp(k) * rho_air(k+1) * ( &3481 ( 37.0_wp * ibit26 * adv_mom_5 &3482 + 7.0_wp * ibit25 * adv_mom_3 &3483 +