Changeset 4649 for palm/trunk/SOURCE/pres.f90
- Timestamp:
- Aug 25, 2020 12:11:17 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/pres.f90 (modified) (41 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pres.f90
r4457 r4649 1 1 !> @file pres.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 5 ! PALM is free software: you can redistribute it and/or modify it under the 6 ! terms of the GNU General Public License as published by the Free Software 7 ! Foundation, either version 3 of the License, or (at your option) any later 8 ! version. 9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY 11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 12 ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 13 ! 14 ! You should have received a copy of the GNU General Public License along with 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------! 17 !--------------------------------------------------------------------------------------------------! 18 ! 19 19 ! 20 20 ! Current revisions: 21 ! ----------------- -22 ! 23 ! 21 ! ----------------- 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! File re-formatted to follow the PALM coding standard 28 ! 29 ! 30 ! 4457 2020-03-11 14:20:43Z raasch 27 31 ! use statement for exchange horiz added 28 ! 32 ! 29 33 ! 4360 2020-01-07 11:25:50Z suehring 30 34 ! Introduction of wall_flags_total_0, which currently sets bits based on static 31 35 ! topography information used in wall_flags_static_0 32 ! 36 ! 33 37 ! 4329 2019-12-10 15:46:36Z motisi 34 38 ! Renamed wall_flags_0 to wall_flags_static_0 35 ! 39 ! 36 40 ! 4182 2019-08-22 15:20:23Z scharf 37 41 ! Corrected "Former revisions" section 38 ! 42 ! 39 43 ! 4015 2019-06-05 13:25:35Z raasch 40 44 ! variable child_domain_nvn eliminated 41 ! 45 ! 42 46 ! 3849 2019-04-01 16:35:16Z knoop 43 47 ! OpenACC port for SPEC … … 47 51 ! 48 52 ! 53 !--------------------------------------------------------------------------------------------------! 49 54 ! Description: 50 55 ! ------------ 51 !> Compute the divergence of the provisional velocity field. Solve the Poisson 52 !> equation for the perturbation pressure. Compute the final velocities using53 !> this perturbation pressure. Compute theremaining divergence.54 !------------------------------------------------------------------------------ !56 !> Compute the divergence of the provisional velocity field. Solve the Poisson equation for the 57 !> perturbation pressure. Compute the final velocities using this perturbation pressure. Compute the 58 !> remaining divergence. 59 !--------------------------------------------------------------------------------------------------! 55 60 SUBROUTINE pres 56 57 58 USE arrays_3d, & 59 ONLY: d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw, & 60 tend, u, v, w 61 62 USE control_parameters, & 63 ONLY: bc_lr_cyc, bc_ns_cyc, bc_radiation_l, bc_radiation_n, & 64 bc_radiation_r, bc_radiation_s, child_domain, & 65 conserve_volume_flow, coupling_mode, & 66 dt_3d, gathered_size, ibc_p_b, ibc_p_t, & 67 intermediate_timestep_count, intermediate_timestep_count_max, & 68 mg_switch_to_pe0_level, nesting_offline, & 69 psolver, subdomain_size, & 70 topography, volume_flow, volume_flow_area, volume_flow_initial 71 72 USE cpulog, & 73 ONLY: cpu_log, log_point, log_point_s 74 75 USE exchange_horiz_mod, & 61 62 63 USE arrays_3d, & 64 ONLY: d, & 65 ddzu, & 66 ddzu_pres, & 67 ddzw, & 68 dzw, & 69 p, & 70 p_loc, & 71 rho_air, & 72 rho_air_zw, & 73 tend, & 74 u, & 75 v, & 76 w 77 78 USE control_parameters, & 79 ONLY: bc_lr_cyc, & 80 bc_ns_cyc, & 81 bc_radiation_l, & 82 bc_radiation_n, & 83 bc_radiation_r, & 84 bc_radiation_s, & 85 child_domain, & 86 conserve_volume_flow, & 87 coupling_mode, & 88 dt_3d, & 89 gathered_size, & 90 ibc_p_b, & 91 ibc_p_t, & 92 intermediate_timestep_count, & 93 intermediate_timestep_count_max, & 94 mg_switch_to_pe0_level, & 95 nesting_offline, & 96 psolver, & 97 subdomain_size, & 98 topography, & 99 volume_flow, & 100 volume_flow_area, & 101 volume_flow_initial 102 103 USE cpulog, & 104 ONLY: cpu_log, & 105 log_point, & 106 log_point_s 107 108 USE exchange_horiz_mod, & 76 109 ONLY: exchange_horiz 77 110 78 USE grid_variables, & 79 ONLY: ddx, ddy 80 81 USE indices, & 82 ONLY: nbgp, ngp_2dh_outer, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, & 83 ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzt, nzt_mg, & 111 USE grid_variables, & 112 ONLY: ddx, & 113 ddy 114 115 USE indices, & 116 ONLY: nbgp, & 117 ngp_2dh_outer, & 118 nx, & 119 nxl, & 120 nxlg, & 121 nxl_mg, & 122 nxr, & 123 nxrg, & 124 nxr_mg, & 125 ny, & 126 nys, & 127 nysg, & 128 nys_mg, & 129 nyn, & 130 nyng, & 131 nyn_mg, & 132 nzb, & 133 nzt, & 134 nzt_mg, & 84 135 wall_flags_total_0 85 136 … … 87 138 88 139 USE pegrid 89 90 USE pmc_interface, &91 ONLY: nesting_mode 92 93 USE poisfft_mod, &140 141 USE pmc_interface, & 142 ONLY: nesting_mode 143 144 USE poisfft_mod, & 94 145 ONLY: poisfft 95 146 … … 98 149 USE poismg_noopt_mod 99 150 100 USE statistics, & 101 ONLY: statistic_regions, sums_divnew_l, sums_divold_l, weight_pres, & 151 USE statistics, & 152 ONLY: statistic_regions, & 153 sums_divnew_l, & 154 sums_divold_l, & 155 weight_pres, & 102 156 weight_substep 103 157 104 USE surface_mod, &158 USE surface_mod, & 105 159 ONLY : bc_h 106 160 107 161 IMPLICIT NONE 108 162 109 INTEGER(iwp) :: i 110 INTEGER(iwp) :: j 111 INTEGER(iwp) :: k 112 INTEGER(iwp) :: m 113 114 REAL(wp) :: ddt_3d!<115 REAL(wp) :: d_weight_pres!<116 REAL(wp) :: localsum!<117 REAL(wp) :: threadsum!<118 REAL(wp) :: weight_pres_l!<119 REAL(wp) :: weight_substep_l!<163 INTEGER(iwp) :: i !< 164 INTEGER(iwp) :: j !< 165 INTEGER(iwp) :: k !< 166 INTEGER(iwp) :: m !< 167 168 REAL(wp) :: ddt_3d !< 169 REAL(wp) :: d_weight_pres !< 170 REAL(wp) :: localsum !< 171 REAL(wp) :: threadsum !< 172 REAL(wp) :: weight_pres_l !< 173 REAL(wp) :: weight_substep_l !< 120 174 121 175 REAL(wp), DIMENSION(1:3) :: volume_flow_l !< … … 144 198 ! 145 199 !-- Multigrid method expects array d to have one ghost layer. 146 !-- 200 !-- 147 201 IF ( psolver(1:9) == 'multigrid' ) THEN 148 202 149 203 DEALLOCATE( d ) 150 ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 151 152 ! 153 !-- Since p is later used to hold the weighted average of the substeps, it 154 !-- cannot be used in the iterative solver. Therefore, its initial value is155 !-- stored on p_loc, which is then iterativelyadvanced in every substep.204 ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 205 206 ! 207 !-- Since p is later used to hold the weighted average of the substeps, it cannot be used in the 208 !-- iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively 209 !-- advanced in every substep. 156 210 IF ( intermediate_timestep_count <= 1 ) THEN 157 211 DO i = nxl-1, nxr+1 … … 163 217 ENDDO 164 218 ENDIF 165 219 166 220 ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN 167 221 168 222 ! 169 !-- Since p is later used to hold the weighted average of the substeps, it 170 !-- cannot be used in the iterative solver. Therefore, its initial value is171 !-- stored on p_loc, which is then iterativelyadvanced in every substep.223 !-- Since p is later used to hold the weighted average of the substeps, it cannot be used in the 224 !-- iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively 225 !-- advanced in every substep. 172 226 p_loc = p 173 227 … … 175 229 176 230 ! 177 !-- Conserve the volume flow at the outflow in case of non-cyclic lateral 178 !-- boundary conditions 179 !-- WARNING: so far, this conservation does not work at the left/south 180 !-- boundary if the topography at the inflow differs from that at the 181 !-- outflow! For this case, volume_flow_area needs adjustment! 231 !-- Conserve the volume flow at the outflow in case of non-cyclic lateral boundary conditions 232 !-- WARNING: so far, this conservation does not work at the left/south boundary if the topography at 233 !-- the inflow differs from that at the outflow! For this case, volume_flow_area needs 234 !-- adjustment! 182 235 ! 183 236 !-- Left/right 184 IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. & 185 bc_radiation_r ) ) THEN 237 IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. bc_radiation_r ) ) THEN 186 238 187 239 volume_flow(1) = 0.0_wp … … 198 250 !-- Sum up the volume flow through the south/north boundary 199 251 DO k = nzb+1, nzt 200 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) & 201 * MERGE( 1.0_wp, 0.0_wp, & 202 BTEST( wall_flags_total_0(k,j,i), 1 ) & 203 ) 204 ENDDO 205 ENDDO 206 207 #if defined( __parallel ) 252 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) & 253 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 254 ENDDO 255 ENDDO 256 257 #if defined( __parallel ) 208 258 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 209 CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, & 210 MPI_SUM, comm1dy, ierr ) 259 CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, MPI_SUM, comm1dy, ierr ) 211 260 #else 212 volume_flow = volume_flow_l 261 volume_flow = volume_flow_l 213 262 #endif 214 volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) & 215 / volume_flow_area(1) 263 volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) / volume_flow_area(1) 216 264 217 265 DO j = nysg, nyng 218 266 DO k = nzb+1, nzt 219 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) & 220 * MERGE( 1.0_wp, 0.0_wp, & 221 BTEST( wall_flags_total_0(k,j,i), 1 ) & 222 ) 267 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) & 268 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 223 269 ENDDO 224 270 ENDDO … … 243 289 !-- Sum up the volume flow through the south/north boundary 244 290 DO k = nzb+1, nzt 245 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) & 246 * MERGE( 1.0_wp, 0.0_wp, & 247 BTEST( wall_flags_total_0(k,j,i), 2 ) & 248 ) 249 ENDDO 250 ENDDO 251 252 #if defined( __parallel ) 291 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) & 292 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 293 ENDDO 294 ENDDO 295 296 #if defined( __parallel ) 253 297 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 254 CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, & 255 MPI_SUM, comm1dx, ierr ) 298 CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, MPI_SUM, comm1dx, ierr ) 256 299 #else 257 volume_flow = volume_flow_l 300 volume_flow = volume_flow_l 258 301 #endif 259 volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) & 260 / volume_flow_area(2) 302 volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) / volume_flow_area(2) 261 303 262 304 DO i = nxlg, nxrg 263 305 DO k = nzb+1, nzt 264 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) & 265 * MERGE( 1.0_wp, 0.0_wp, & 266 BTEST( wall_flags_total_0(k,j,i), 2 ) & 267 ) 306 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) & 307 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 268 308 ENDDO 269 309 ENDDO … … 273 313 ! 274 314 !-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top 275 !-- boundary. With Neumann conditions at both vertical boundaries, the solver cannot remove 276 !-- mean vertical velocities. They should be removed, because incompressibility requires that277 !-- the vertical gradient of vertical velocity is zero. Since w=0 at the solid surface, it must be278 !-- zeroeverywhere.315 !-- boundary. With Neumann conditions at both vertical boundaries, the solver cannot remove mean 316 !-- vertical velocities. They should be removed, because incompressibility requires that the 317 !-- vertical gradient of vertical velocity is zero. Since w=0 at the solid surface, it must be zero 318 !-- everywhere. 279 319 !-- This must not be done in case of a 3d-nesting child domain, because a mean vertical velocity 280 320 !-- can physically exist in such a domain. … … 282 322 !-- caused by horizontal divergence/convergence of the large scale flow that is prescribed at the 283 323 !-- side boundaries. 284 !-- The removal cannot be done before the first initial time step because ngp_2dh_outer 285 !-- is not yetknown then.324 !-- The removal cannot be done before the first initial time step because ngp_2dh_outer is not yet 325 !-- known then. 286 326 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline & 287 327 .AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' ) & 288 .AND. intermediate_timestep_count /= 0 ) & 289 THEN 328 .AND. intermediate_timestep_count /= 0 ) THEN 290 329 w_l = 0.0_wp; w_l_l = 0.0_wp 291 330 DO i = nxl, nxr … … 293 332 DO k = nzb+1, nzt 294 333 w_l_l(k) = w_l_l(k) + w(k,j,i) & 295 296 ENDDO 297 ENDDO 298 ENDDO 299 #if defined( __parallel ) 334 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 335 ENDDO 336 ENDDO 337 ENDDO 338 #if defined( __parallel ) 300 339 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 301 340 CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, ierr ) … … 310 349 DO k = nzb+1, nzt 311 350 w(k,j,i) = w(k,j,i) - w_l(k) & 312 351 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 313 352 ENDDO 314 353 ENDDO … … 351 390 DO j = nys, nyn 352 391 DO k = nzb+1, nzt 353 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 354 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 355 ( w(k,j,i) * rho_air_zw(k) - & 356 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 357 ) * ddt_3d * d_weight_pres & 358 * MERGE( 1.0_wp, 0.0_wp, & 359 BTEST( wall_flags_total_0(k,j,i), 0 ) & 360 ) 392 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 393 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 394 ( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) & 395 * ddzw(k) ) * ddt_3d * d_weight_pres & 396 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 361 397 ENDDO 362 398 ! 363 399 !-- Compute possible PE-sum of divergences for flow_statistics 364 400 DO k = nzb+1, nzt 365 threadsum = threadsum + ABS( d(k,j,i) ) & 366 * MERGE( 1.0_wp, 0.0_wp, & 367 BTEST( wall_flags_total_0(k,j,i), 0 ) & 368 ) 401 threadsum = threadsum + ABS( d(k,j,i) ) & 402 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 369 403 ENDDO 370 404 … … 372 406 ENDDO 373 407 374 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &408 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. & 375 409 intermediate_timestep_count == 0 ) THEN 376 410 localsum = localsum + threadsum * dt_3d * weight_pres_l … … 387 421 DO j = nys, nyn 388 422 DO k = 1, nzt 389 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 390 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 391 ( w(k,j,i) * rho_air_zw(k) - & 392 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 393 ) * ddt_3d * d_weight_pres & 394 * MERGE( 1.0_wp, 0.0_wp, & 395 BTEST( wall_flags_total_0(k,j,i), 0 ) & 396 ) 423 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 424 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 425 ( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) & 426 * ddzw(k) ) * ddt_3d * d_weight_pres & 427 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 397 428 ENDDO 398 429 ENDDO … … 401 432 402 433 ! 403 !-- Compute possible PE-sum of divergences for flow_statistics. Carry out 404 !-- computation only at lastRunge-Kutta substep.405 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &434 !-- Compute possible PE-sum of divergences for flow_statistics. Carry out computation only at last 435 !-- Runge-Kutta substep. 436 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. & 406 437 intermediate_timestep_count == 0 ) THEN 407 438 !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum) … … 423 454 424 455 ! 425 !-- For completeness, set the divergence sum of all statistic regions to those 426 !-- of the total domain 427 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. & 456 !-- For completeness, set the divergence sum of all statistic regions to those of the total domain 457 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. & 428 458 intermediate_timestep_count == 0 ) THEN 429 459 sums_divold_l(0:statistic_regions) = localsum … … 441 471 442 472 ! 443 !-- Store computed perturbation pressure and set boundary condition in 444 !-- z-direction 473 !-- Store computed perturbation pressure and set boundary condition in z-direction 445 474 !$OMP PARALLEL DO PRIVATE (i,j,k) 446 475 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & … … 456 485 ! 457 486 !-- Bottom boundary: 458 !-- This condition is only required for internal output. The pressure 459 !-- gradient(dp(nzb+1)-dp(nzb))/dz is not used anywhere else.487 !-- This condition is only required for internal output. The pressure gradient 488 !-- (dp(nzb+1)-dp(nzb))/dz is not used anywhere else. 460 489 IF ( ibc_p_b == 1 ) THEN 461 490 ! 462 !-- Neumann (dp/dz = 0). Using surfa e data type, first for non-natural463 !-- surfaces, then fornatural and urban surfaces491 !-- Neumann (dp/dz = 0). Using surface data type, first for non-natural surfaces, then for 492 !-- natural and urban surfaces 464 493 !-- Upward facing 465 494 !$OMP PARALLEL DO PRIVATE( i, j, k ) … … 486 515 ELSE 487 516 ! 488 !-- Dirichlet. Using surface data type, first for non-natural 489 !-- surfaces, then for natural andurban surfaces517 !-- Dirichlet. Using surface data type, first for non-natural surfaces, then for natural and 518 !-- urban surfaces 490 519 !-- Upward facing 491 520 !$OMP PARALLEL DO PRIVATE( i, j, k ) … … 537 566 !-- Exchange boundaries for p 538 567 CALL exchange_horiz( tend, nbgp ) 539 568 540 569 ELSEIF ( psolver == 'sor' ) THEN 541 570 542 571 ! 543 !-- Solve Poisson equation for perturbation pressure using SOR-Red/Black 544 !-- scheme 572 !-- Solve Poisson equation for perturbation pressure using SOR-Red/Black scheme 545 573 CALL sor( d, ddzu_pres, ddzw, p_loc ) 546 574 tend = p_loc … … 549 577 550 578 ! 551 !-- Solve Poisson equation for perturbation pressure using Multigrid scheme, 552 !-- array tend is used to store the residuals. 553 554 !-- If the number of grid points of the gathered grid, which is collected 555 !-- on PE0, is larger than the number of grid points of an PE, than array 556 !-- tend will be enlarged. 579 !-- Solve Poisson equation for perturbation pressure using Multigrid scheme, array tend is used 580 !-- to store the residuals. 581 582 !-- If the number of grid points of the gathered grid, which is collected on PE0, is larger than 583 !-- the number of grid points of an PE, than array tend will be enlarged. 557 584 IF ( gathered_size > subdomain_size ) THEN 558 585 DEALLOCATE( tend ) 559 ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( & 560 mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,& 561 nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( & 562 mg_switch_to_pe0_level)+1) ) 586 ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( & 587 mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1, & 588 nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(mg_switch_to_pe0_level)+1) ) 563 589 ENDIF 564 590 … … 575 601 576 602 ! 577 !-- Restore perturbation pressure on tend because this array is used 578 !-- further below to correct thevelocity fields603 !-- Restore perturbation pressure on tend because this array is used further below to correct the 604 !-- velocity fields 579 605 DO i = nxl-1, nxr+1 580 606 DO j = nys-1, nyn+1 … … 598 624 DO j = nys-1, nyn+1 599 625 DO k = nzb, nzt+1 600 p(k,j,i) = tend(k,j,i) * & 601 weight_substep_l 626 p(k,j,i) = tend(k,j,i) * weight_substep_l 602 627 ENDDO 603 628 ENDDO … … 613 638 DO j = nys-1, nyn+1 614 639 DO k = nzb, nzt+1 615 p(k,j,i) = p(k,j,i) + tend(k,j,i) * & 616 weight_substep_l 640 p(k,j,i) = p(k,j,i) + tend(k,j,i) * weight_substep_l 617 641 ENDDO 618 642 ENDDO … … 621 645 622 646 ENDIF 623 647 624 648 ! 625 649 !-- SOR-method needs ghost layers for the next timestep … … 627 651 628 652 ! 629 !-- Correction of the provisional velocities with the current perturbation 630 !-- pressure just computed 653 !-- Correction of the provisional velocities with the current perturbation pressure just computed 631 654 IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN 632 655 volume_flow_l(1) = 0.0_wp … … 634 657 ENDIF 635 658 ! 636 !-- Add pressure gradients to the velocity components. Note, the loops are 637 !-- running over the entire model domain, even though, in case of non-cyclic 638 !-- boundaries u- and v-component are not prognostic at i=0 and j=0, 639 !-- respectiveley. However, in case of Dirichlet boundary conditions for the 640 !-- velocities, zero-gradient conditions for the pressure are set, so that 641 !-- no modification is imposed at the boundaries. 659 !-- Add pressure gradients to the velocity components. Note, the loops are running over the entire 660 !-- model domain, even though, in case of non-cyclic boundaries u- and v-component are not 661 !-- prognostic at i=0 and j=0, respectiveley. However, in case of Dirichlet boundary conditions for 662 !-- the velocities, zero-gradient conditions for the pressure are set, so that no modification is 663 !-- imposed at the boundaries. 642 664 !$OMP PARALLEL PRIVATE (i,j,k) 643 665 !$OMP DO 644 666 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) & 645 667 !$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0) 646 DO i = nxl, nxr 668 DO i = nxl, nxr 647 669 DO j = nys, nyn 648 670 649 671 DO k = nzb+1, nzt 650 w(k,j,i) = w(k,j,i) - dt_3d * & 651 ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) & 652 * weight_pres_l & 653 * MERGE( 1.0_wp, 0.0_wp, & 654 BTEST( wall_flags_total_0(k,j,i), 3 ) & 655 ) 656 ENDDO 657 658 DO k = nzb+1, nzt 659 u(k,j,i) = u(k,j,i) - dt_3d * & 660 ( tend(k,j,i) - tend(k,j,i-1) ) * ddx & 661 * weight_pres_l & 662 * MERGE( 1.0_wp, 0.0_wp, & 663 BTEST( wall_flags_total_0(k,j,i), 1 ) & 664 ) 665 ENDDO 666 667 DO k = nzb+1, nzt 668 v(k,j,i) = v(k,j,i) - dt_3d * & 669 ( tend(k,j,i) - tend(k,j-1,i) ) * ddy & 670 * weight_pres_l & 671 * MERGE( 1.0_wp, 0.0_wp, & 672 BTEST( wall_flags_total_0(k,j,i), 2 ) & 673 ) 674 ENDDO 672 w(k,j,i) = w(k,j,i) - dt_3d * ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) & 673 * weight_pres_l & 674 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 675 ENDDO 676 677 DO k = nzb+1, nzt 678 u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx * weight_pres_l & 679 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 680 ENDDO 681 682 DO k = nzb+1, nzt 683 v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy * weight_pres_l & 684 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 685 ENDDO 675 686 676 687 ENDDO … … 679 690 680 691 ! 681 !-- The vertical velocity is not set to zero at nzt + 1 for nested domains 682 !-- Instead it is set to the values of nzt (see routine vnest_boundary_conds 683 !-- or pmci_interp_tril_t) BEFORE calling the pressure solver. To avoid jumps 684 !-- while plotting profiles w at the top has to be set to the values in the 685 !-- height nzt after above modifications. Hint: w level nzt+1 does not impact 686 !-- results. 687 IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN 692 !-- The vertical velocity is not set to zero at nzt + 1 for nested domains. Instead it is set to the 693 !-- values of nzt (see routine vnest_boundary_conds or pmci_interp_tril_t) BEFORE calling the 694 !-- pressure solver. To avoid jumps while plotting profiles, w at the top has to be set to the 695 !-- values in height nzt after above modifications. Hint: w level nzt+1 does not impact results. 696 IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN 688 697 w(nzt+1,:,:) = w(nzt,:,:) 689 698 ENDIF … … 691 700 ! 692 701 !-- Sum up the volume flow through the right and north boundary 693 IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. & 694 nxr == nx ) THEN 702 IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. nxr == nx ) THEN 695 703 696 704 !$OMP PARALLEL PRIVATE (j,k) … … 699 707 !$OMP CRITICAL 700 708 DO k = nzb+1, nzt 701 volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) & 702 * MERGE( 1.0_wp, 0.0_wp, & 703 BTEST( wall_flags_total_0(k,j,nxr), 1 )& 704 ) 709 volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) & 710 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr), 1 ) ) 705 711 ENDDO 706 712 !$OMP END CRITICAL … … 710 716 ENDIF 711 717 712 IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. & 713 nyn == ny ) THEN 718 IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. nyn == ny ) THEN 714 719 715 720 !$OMP PARALLEL PRIVATE (i,k) … … 718 723 !$OMP CRITICAL 719 724 DO k = nzb+1, nzt 720 volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) & 721 * MERGE( 1.0_wp, 0.0_wp, & 722 BTEST( wall_flags_total_0(k,nyn,i), 2 )& 723 ) 725 volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) & 726 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn,i), 2 ) ) 724 727 ENDDO 725 728 !$OMP END CRITICAL … … 728 731 729 732 ENDIF 730 733 731 734 ! 732 735 !-- Conserve the volume flow 733 736 IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN 734 737 735 #if defined( __parallel ) 738 #if defined( __parallel ) 736 739 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 737 CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, & 738 MPI_SUM, comm2d, ierr ) 740 CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, MPI_SUM, comm2d, ierr ) 739 741 #else 740 volume_flow = volume_flow_l 741 #endif 742 743 volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) /&744 volume_flow_area(1:2)742 volume_flow = volume_flow_l 743 #endif 744 745 volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) & 746 / volume_flow_area(1:2) 745 747 746 748 !$OMP PARALLEL PRIVATE (i,j,k) … … 749 751 DO j = nys, nyn 750 752 DO k = nzb+1, nzt 751 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) & 752 * MERGE( 1.0_wp, 0.0_wp, & 753 BTEST( wall_flags_total_0(k,j,i), 1 ) & 754 ) 755 ENDDO 756 DO k = nzb+1, nzt 757 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) & 758 * MERGE( 1.0_wp, 0.0_wp, & 759 BTEST( wall_flags_total_0(k,j,i), 2 ) & 760 ) 753 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) & 754 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 755 ENDDO 756 DO k = nzb+1, nzt 757 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) & 758 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 761 759 ENDDO 762 760 ENDDO … … 774 772 775 773 ! 776 !-- Compute the divergence of the corrected velocity field ,777 !-- A possible PE-sum is computed in flow_statistics. Carry out computation 778 !-- only at lastRunge-Kutta step.779 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &774 !-- Compute the divergence of the corrected velocity field. 775 !-- A possible PE-sum is computed in flow_statistics. Carry out computation only at last 776 !-- Runge-Kutta step. 777 IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. & 780 778 intermediate_timestep_count == 0 ) THEN 781 779 CALL cpu_log( log_point_s(1), 'divergence', 'start' ) … … 783 781 784 782 ! 785 !-- d must be reset to zero because it can contain nonzero values below the 786 !-- topography 783 !-- d must be reset to zero because it can contain nonzero values below the topography 787 784 IF ( topography /= 'flat' ) d = 0.0_wp 788 785 … … 796 793 DO j = nys, nyn 797 794 DO k = nzb+1, nzt 798 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 799 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 800 ( w(k,j,i) * rho_air_zw(k) - & 801 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 802 ) * MERGE( 1.0_wp, 0.0_wp, & 803 BTEST( wall_flags_total_0(k,j,i), 0 ) & 804 ) 795 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 796 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 797 ( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) & 798 * ddzw(k) ) & 799 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 805 800 ENDDO 806 801 DO k = nzb+1, nzt … … 819 814 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 820 815 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 821 ( w(k,j,i) * rho_air_zw(k) - & 822 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 823 ) * MERGE( 1.0_wp, 0.0_wp, & 824 BTEST( wall_flags_total_0(k,j,i), 0 ) & 825 ) 816 ( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) & 817 * ddzw(k) ) & 818 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 826 819 ENDDO 827 820 ENDDO … … 846 839 847 840 ! 848 !-- For completeness, set the divergence sum of all statistic regions to those 849 !-- of the totaldomain841 !-- For completeness, set the divergence sum of all statistic regions to those of the total 842 !-- domain 850 843 sums_divnew_l(0:statistic_regions) = localsum 851 844
Note: See TracChangeset
for help on using the changeset viewer.
