Changeset 4649 for palm/trunk/SOURCE/poismg_mod.f90
- Timestamp:
- Aug 25, 2020 12:11:17 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/poismg_mod.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/poismg_mod.f90
r4457 r4649 1 1 !> @file poismg.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 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! use statement for exchange horiz added 28 ! 27 ! File re-formatted to follow the PALM coding standard 28 ! 29 ! 30 ! 4457 2020-03-11 14:20:43Z raasch 31 ! Use statement for exchange horiz added 32 ! 29 33 ! 4432 2020-02-28 07:43:21Z raasch 30 ! bugfix for previous revision (vector directive was changed by mistake)31 ! 34 ! Bugfix for previous revision (vector directive was changed by mistake) 35 ! 32 36 ! 4429 2020-02-27 15:24:30Z raasch 33 ! statement added to avoid compile error due to unused dummy argument34 ! bugfix: cpp-directives added for serial mode35 ! 37 ! Statement added to avoid compile error due to unused dummy argument 38 ! Bugfix: cpp-directives added for serial mode 39 ! 36 40 ! 4360 2020-01-07 11:25:50Z suehring 37 41 ! Corrected "Former revisions" section 38 ! 42 ! 39 43 ! 3725 2019-02-07 10:11:02Z raasch 40 ! unused subroutine removed41 ! 44 ! Unused subroutine removed 45 ! 42 46 ! 3655 2019-01-07 16:51:22Z knoop 43 ! unnecessary check eliminated47 ! Unnecessary check eliminated 44 48 ! 45 49 ! Following optimisations have been made: 46 ! - vectorisation (for Intel-CPUs) of the red-black algorithm by resorting 47 ! array elements with even and odd indices 48 ! - explicit boundary conditions for building walls removed (solver is 49 ! running through the buildings 50 ! - reduced data transfer in case of ghost point exchange, because only 51 ! "red" or "black" data points need to be exchanged. This is not applied 52 ! for coarser grid levels, since for then the transfer time is latency bound 53 ! 54 ! 50 ! - Vectorisation (for Intel-CPUs) of the red-black algorithm by resorting array elements with even 51 ! and odd indices 52 ! - Explicit boundary conditions for building walls removed (solver is running through the 53 ! buildings) 54 ! - Reduced data transfer in case of ghost point exchange, because only "red" or "black" data points 55 ! need to be exchanged. This is not applied for coarser grid levels, since for then the transfer 56 ! time is latency bound 57 ! 58 ! 59 !--------------------------------------------------------------------------------------------------! 55 60 ! Description: 56 61 ! ------------ 57 !> Solves the Poisson equation for the perturbation pressure with a multigrid 58 !> V- or W-Cycle scheme. 62 !> Solves the Poisson equation for the perturbation pressure with a multigrid V- or W-Cycle scheme. 59 63 !> 60 64 !> This multigrid method was originally developed for PALM by Joerg Uhlenbrock, 61 !> September 2000 - July 2001. It has been optimised for speed by Klaus 62 !> Ketelsen in November 2014. 63 !> 64 !> @attention Loop unrolling and cache optimization in SOR-Red/Black method 65 !> still does not give the expected speedup! 65 !> September 2000 - July 2001. It has been optimised for speed by Klaus Ketelsen in November 2014. 66 !> 67 !> @attention Loop unrolling and cache optimization in SOR-Red/Black method still does not give the 68 ! expected speedup! 66 69 !> 67 70 !> @todo Further work required. 68 !------------------------------------------------------------------------------ !71 !--------------------------------------------------------------------------------------------------! 69 72 MODULE poismg_mod 70 71 USE control_parameters, & 72 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 73 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, & 74 bc_radiation_s, grid_level, nesting_offline 75 76 USE cpulog, & 77 ONLY: cpu_log, log_point_s 78 79 USE exchange_horiz_mod, & 73 74 USE control_parameters, & 75 ONLY: bc_dirichlet_l, & 76 bc_dirichlet_n, & 77 bc_dirichlet_r, & 78 bc_dirichlet_s, & 79 bc_radiation_l, & 80 bc_radiation_n, & 81 bc_radiation_r, & 82 bc_radiation_s, & 83 grid_level, & 84 nesting_offline 85 86 USE cpulog, & 87 ONLY: cpu_log, & 88 log_point_s 89 90 USE exchange_horiz_mod, & 80 91 ONLY: exchange_horiz 81 92 … … 87 98 88 99 INTEGER, SAVE :: ind_even_odd !< border index between even and odd k index 100 89 101 INTEGER, DIMENSION(:), SAVE, ALLOCATABLE :: even_odd_level !< stores ind_even_odd for all MG levels 90 102 … … 105 117 CONTAINS 106 118 107 !------------------------------------------------------------------------------ !119 !--------------------------------------------------------------------------------------------------! 108 120 ! Description: 109 121 ! ------------ 110 !> Solves the Poisson equation for the perturbation pressure with a multigrid 111 !> V- or W-Cycle scheme. 112 !------------------------------------------------------------------------------! 113 SUBROUTINE poismg( r ) 114 115 USE arrays_3d, & 116 ONLY: d, p_loc 117 118 USE control_parameters, & 119 ONLY: bc_lr_cyc, bc_ns_cyc, gathered_size, grid_level, & 120 grid_level_count, ibc_p_t, & 121 maximum_grid_level, message_string, mgcycles, mg_cycles, & 122 mg_switch_to_pe0_level, residual_limit, subdomain_size 123 124 USE cpulog, & 125 ONLY: cpu_log, log_point_s 126 127 USE indices, & 128 ONLY: nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn,& 129 nyng, nyn_mg, nzb, nzt, nzt_mg 130 131 IMPLICIT NONE 132 133 REAL(wp) :: maxerror !< 134 REAL(wp) :: maximum_mgcycles !< 135 REAL(wp) :: residual_norm !< 136 137 REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !< 138 139 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !< 140 141 142 CALL cpu_log( log_point_s(29), 'poismg', 'start' ) 143 ! 144 !-- Initialize arrays and variables used in this subroutine 145 146 !-- If the number of grid points of the gathered grid, which is collected 147 !-- on PE0, is larger than the number of grid points of an PE, than array 148 !-- p3 will be enlarged. 149 IF ( gathered_size > subdomain_size ) THEN 150 ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( & 151 mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,& 152 nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( & 153 mg_switch_to_pe0_level)+1) ) 154 ELSE 155 ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 122 !> Solves the Poisson equation for the perturbation pressure with a multigrid V- or W-Cycle scheme. 123 !--------------------------------------------------------------------------------------------------! 124 SUBROUTINE poismg( r ) 125 126 USE arrays_3d, & 127 ONLY: d, & 128 p_loc 129 130 USE control_parameters, & 131 ONLY: bc_lr_cyc, & 132 bc_ns_cyc, & 133 gathered_size, & 134 grid_level, & 135 grid_level_count, & 136 ibc_p_t, & 137 maximum_grid_level, & 138 message_string, & 139 mgcycles, & 140 mg_cycles, & 141 mg_switch_to_pe0_level, & 142 residual_limit, & 143 subdomain_size 144 145 USE cpulog, & 146 ONLY: cpu_log, & 147 log_point_s 148 149 USE indices, & 150 ONLY: nxl, & 151 nxlg, & 152 nxl_mg, & 153 nxr, & 154 nxrg, & 155 nxr_mg, & 156 nys, & 157 nysg, & 158 nys_mg, & 159 nyn, & 160 nyng, & 161 nyn_mg, & 162 nzb, & 163 nzt, & 164 nzt_mg 165 166 IMPLICIT NONE 167 168 REAL(wp) :: maxerror !< 169 REAL(wp) :: maximum_mgcycles !< 170 REAL(wp) :: residual_norm !< 171 172 REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !< 173 174 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !< 175 176 177 CALL cpu_log( log_point_s(29), 'poismg', 'start' ) 178 ! 179 !-- Initialize arrays and variables used in this subroutine 180 181 !-- If the number of grid points of the gathered grid, which is collected on PE0, is larger than the 182 !-- number of grid points of a PE, than array p3 will be enlarged. 183 IF ( gathered_size > subdomain_size ) THEN 184 ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( & 185 mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1, & 186 nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(mg_switch_to_pe0_level)+1) ) 187 ELSE 188 ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 189 ENDIF 190 191 p3 = 0.0_wp 192 193 194 ! 195 !-- Ghost boundaries have to be added to divergence array. 196 !-- Exchange routine needs to know the grid level! 197 grid_level = maximum_grid_level 198 CALL exchange_horiz( d, 1) 199 ! 200 !-- Set bottom and top boundary conditions 201 d(nzb,:,:) = d(nzb+1,:,:) 202 IF ( ibc_p_t == 1 ) d(nzt+1,:,: ) = d(nzt,:,:) 203 ! 204 !-- Set lateral boundary conditions in non-cyclic case 205 IF ( .NOT. bc_lr_cyc ) THEN 206 IF ( bc_dirichlet_l .OR. bc_radiation_l ) d(:,:,nxl-1) = d(:,:,nxl) 207 IF ( bc_dirichlet_r .OR. bc_radiation_r ) d(:,:,nxr+1) = d(:,:,nxr) 208 ENDIF 209 IF ( .NOT. bc_ns_cyc ) THEN 210 IF ( bc_dirichlet_n .OR. bc_radiation_n ) d(:,nyn+1,:) = d(:,nyn,:) 211 IF ( bc_dirichlet_s .OR. bc_radiation_s ) d(:,nys-1,:) = d(:,nys,:) 212 ENDIF 213 ! 214 !-- Initiation of the multigrid scheme. Does n cycles until the residual is smaller than the given 215 !-- limit. The accuracy of the solution of the poisson equation will increase with the number of 216 !-- cycles. If the number of cycles is preset by the user, this number will be carried out 217 !-- regardless of the accuracy. 218 grid_level_count = 0 219 mgcycles = 0 220 IF ( mg_cycles == -1 ) THEN 221 maximum_mgcycles = 0 222 residual_norm = 1.0_wp 223 ELSE 224 maximum_mgcycles = mg_cycles 225 residual_norm = 0.0_wp 226 ENDIF 227 228 ! 229 !-- Initial settings for sorting k-dimension from sequential order (alternate even/odd) into blocks 230 !-- of even and odd or vice versa 231 CALL init_even_odd_blocks 232 233 ! 234 !-- Sort input arrays in even/odd blocks along k-dimension 235 CALL sort_k_to_even_odd_blocks( d, grid_level ) 236 CALL sort_k_to_even_odd_blocks( p_loc, grid_level ) 237 238 ! 239 !-- The complete multigrid cycles are running in block mode, i.e. over seperate data blocks of even 240 !-- and odd indices 241 DO WHILE ( residual_norm > residual_limit .OR. mgcycles < maximum_mgcycles ) 242 243 CALL next_mg_level( d, p_loc, p3, r) 244 245 ! 246 !-- Calculate the residual if the user has not preset the number of cycles to be performed 247 IF ( maximum_mgcycles == 0 ) THEN 248 CALL resid( d, p_loc, r ) 249 maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 ) 250 251 #if defined( __parallel ) 252 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 253 CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, comm2d, ierr) 254 #else 255 residual_norm = maxerror 256 #endif 257 residual_norm = SQRT( residual_norm ) 156 258 ENDIF 157 259 158 p3 = 0.0_wp 159 160 161 ! 162 !-- Ghost boundaries have to be added to divergence array. 163 !-- Exchange routine needs to know the grid level! 164 grid_level = maximum_grid_level 165 CALL exchange_horiz( d, 1) 166 ! 167 !-- Set bottom and top boundary conditions 168 d(nzb,:,:) = d(nzb+1,:,:) 169 IF ( ibc_p_t == 1 ) d(nzt+1,:,: ) = d(nzt,:,:) 170 ! 171 !-- Set lateral boundary conditions in non-cyclic case 172 IF ( .NOT. bc_lr_cyc ) THEN 173 IF ( bc_dirichlet_l .OR. bc_radiation_l ) & 174 d(:,:,nxl-1) = d(:,:,nxl) 175 IF ( bc_dirichlet_r .OR. bc_radiation_r ) & 176 d(:,:,nxr+1) = d(:,:,nxr) 260 mgcycles = mgcycles + 1 261 262 ! 263 !-- If the user has not limited the number of cycles, stop the run in case of insufficient 264 !-- convergence 265 IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN 266 message_string = 'no sufficient convergence within 1000 cycles' 267 CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 ) 177 268 ENDIF 178 IF ( .NOT. bc_ns_cyc ) THEN 179 IF ( bc_dirichlet_n .OR. bc_radiation_n ) & 180 d(:,nyn+1,:) = d(:,nyn,:) 181 IF ( bc_dirichlet_s .OR. bc_radiation_s ) & 182 d(:,nys-1,:) = d(:,nys,:) 183 ENDIF 184 ! 185 !-- Initiation of the multigrid scheme. Does n cycles until the 186 !-- residual is smaller than the given limit. The accuracy of the solution 187 !-- of the poisson equation will increase with the number of cycles. 188 !-- If the number of cycles is preset by the user, this number will be 189 !-- carried out regardless of the accuracy. 190 grid_level_count = 0 191 mgcycles = 0 192 IF ( mg_cycles == -1 ) THEN 193 maximum_mgcycles = 0 194 residual_norm = 1.0_wp 195 ELSE 196 maximum_mgcycles = mg_cycles 197 residual_norm = 0.0_wp 198 ENDIF 199 200 ! 201 !-- Initial settings for sorting k-dimension from sequential order (alternate 202 !-- even/odd) into blocks of even and odd or vice versa 203 CALL init_even_odd_blocks 204 205 ! 206 !-- Sort input arrays in even/odd blocks along k-dimension 207 CALL sort_k_to_even_odd_blocks( d, grid_level ) 208 CALL sort_k_to_even_odd_blocks( p_loc, grid_level ) 209 210 ! 211 !-- The complete multigrid cycles are running in block mode, i.e. over 212 !-- seperate data blocks of even and odd indices 213 DO WHILE ( residual_norm > residual_limit .OR. & 214 mgcycles < maximum_mgcycles ) 215 216 CALL next_mg_level( d, p_loc, p3, r) 217 218 ! 219 !-- Calculate the residual if the user has not preset the number of 220 !-- cycles to be performed 221 IF ( maximum_mgcycles == 0 ) THEN 222 CALL resid( d, p_loc, r ) 223 maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 ) 224 225 #if defined( __parallel ) 226 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 227 CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, & 228 MPI_SUM, comm2d, ierr) 229 #else 230 residual_norm = maxerror 231 #endif 232 residual_norm = SQRT( residual_norm ) 233 ENDIF 234 235 mgcycles = mgcycles + 1 236 237 ! 238 !-- If the user has not limited the number of cycles, stop the run in case 239 !-- of insufficient convergence 240 IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN 241 message_string = 'no sufficient convergence within 1000 cycles' 242 CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 ) 243 ENDIF 244 245 ENDDO 246 247 DEALLOCATE( p3 ) 248 ! 249 !-- Result has to be sorted back from even/odd blocks to sequential order 250 CALL sort_k_to_sequential( p_loc ) 251 ! 252 !-- Unset the grid level. Variable is used to determine the MPI datatypes for 253 !-- ghost point exchange 254 grid_level = 0 255 256 CALL cpu_log( log_point_s(29), 'poismg', 'stop' ) 257 258 END SUBROUTINE poismg 259 260 261 !------------------------------------------------------------------------------! 269 270 ENDDO 271 272 DEALLOCATE( p3 ) 273 ! 274 !-- Result has to be sorted back from even/odd blocks to sequential order 275 CALL sort_k_to_sequential( p_loc ) 276 ! 277 !-- Unset the grid level. Variable is used to determine the MPI datatypes for ghost point exchange 278 grid_level = 0 279 280 CALL cpu_log( log_point_s(29), 'poismg', 'stop' ) 281 282 END SUBROUTINE poismg 283 284 285 !--------------------------------------------------------------------------------------------------! 262 286 ! Description: 263 287 ! ------------ 264 288 !> Computes the residual of the perturbation pressure. 265 !------------------------------------------------------------------------------! 266 SUBROUTINE resid( f_mg, p_mg, r ) 267 268 269 USE arrays_3d, & 270 ONLY: rho_air_mg 271 272 USE control_parameters, & 273 ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t 274 USE grid_variables, & 275 ONLY: ddx2_mg, ddy2_mg 276 277 USE indices, & 278 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 279 280 IMPLICIT NONE 281 282 INTEGER(iwp) :: i !< index variable along x 283 INTEGER(iwp) :: j !< index variable along y 284 INTEGER(iwp) :: k !< index variable along z 285 INTEGER(iwp) :: l !< index indicating grid level 286 INTEGER(iwp) :: km1 !< index variable along z dimension (k-1) 287 INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1) 288 289 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 290 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 291 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< velocity divergence 292 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 293 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 294 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< perturbation pressure 295 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 296 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 297 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !< residuum of perturbation pressure 298 299 ! 300 !-- Calculate the residual 301 l = grid_level 302 303 CALL cpu_log( log_point_s(53), 'resid', 'start' ) 304 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 305 !$OMP DO 306 DO i = nxl_mg(l), nxr_mg(l) 307 DO j = nys_mg(l), nyn_mg(l) 308 !DIR$ IVDEP 309 DO k = ind_even_odd+1, nzt_mg(l) 310 km1 = k-ind_even_odd-1 311 kp1 = k-ind_even_odd 312 r(k,j,i) = f_mg(k,j,i) & 313 - rho_air_mg(k,l) * ddx2_mg(l) * & 314 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 315 - rho_air_mg(k,l) * ddy2_mg(l) * & 316 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 317 - f2_mg_b(k,l) * p_mg(kp1,j,i) & 318 - f3_mg_b(k,l) * p_mg(km1,j,i) & 319 + f1_mg_b(k,l) * p_mg(k,j,i) 320 ENDDO 289 !--------------------------------------------------------------------------------------------------! 290 SUBROUTINE resid( f_mg, p_mg, r ) 291 292 293 USE arrays_3d, & 294 ONLY: rho_air_mg 295 296 USE control_parameters, & 297 ONLY: bc_lr_cyc, & 298 bc_ns_cyc, & 299 ibc_p_b, & 300 ibc_p_t 301 302 USE grid_variables, & 303 ONLY: ddx2_mg, & 304 ddy2_mg 305 306 USE indices, & 307 ONLY: nxl_mg, & 308 nxr_mg, & 309 nys_mg, & 310 nyn_mg, & 311 nzb, & 312 nzt_mg 313 314 IMPLICIT NONE 315 316 INTEGER(iwp) :: i !< index variable along x 317 INTEGER(iwp) :: j !< index variable along y 318 INTEGER(iwp) :: k !< index variable along z 319 INTEGER(iwp) :: l !< index indicating grid level 320 INTEGER(iwp) :: km1 !< index variable along z dimension (k-1) 321 INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1) 322 323 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 324 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< velocity divergence 325 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 326 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< perturbation pressure 327 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 328 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !< residuum of perturbation pressure 329 330 ! 331 !-- Calculate the residual 332 l = grid_level 333 334 CALL cpu_log( log_point_s(53), 'resid', 'start' ) 335 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 336 !$OMP DO 337 DO i = nxl_mg(l), nxr_mg(l) 338 DO j = nys_mg(l), nyn_mg(l) 321 339 !DIR$ IVDEP 322 DO k = nzb+1, ind_even_odd 323 km1 = k+ind_even_odd 324 kp1 = k+ind_even_odd+1 325 r(k,j,i) = f_mg(k,j,i) & 326 - rho_air_mg(k,l) * ddx2_mg(l) * & 327 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 328 - rho_air_mg(k,l) * ddy2_mg(l) * & 329 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 330 - f2_mg_b(k,l) * p_mg(kp1,j,i) & 331 - f3_mg_b(k,l) * p_mg(km1,j,i) & 332 + f1_mg_b(k,l) * p_mg(k,j,i) 333 ENDDO 340 DO k = ind_even_odd+1, nzt_mg(l) 341 km1 = k-ind_even_odd-1 342 kp1 = k-ind_even_odd 343 r(k,j,i) = f_mg(k,j,i) - rho_air_mg(k,l) * ddx2_mg(l) & 344 * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 345 - rho_air_mg(k,l) * ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 346 - f2_mg_b(k,l) * p_mg(kp1,j,i) - f3_mg_b(k,l) * p_mg(km1,j,i) & 347 + f1_mg_b(k,l) * p_mg(k,j,i) 348 ENDDO 349 !DIR$ IVDEP 350 DO k = nzb+1, ind_even_odd 351 km1 = k+ind_even_odd 352 kp1 = k+ind_even_odd+1 353 r(k,j,i) = f_mg(k,j,i) - rho_air_mg(k,l) * ddx2_mg(l) & 354 * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) - rho_air_mg(k,l) * ddy2_mg(l) & 355 * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) - f2_mg_b(k,l) * p_mg(kp1,j,i) & 356 - f3_mg_b(k,l) * p_mg(km1,j,i) + f1_mg_b(k,l) * p_mg(k,j,i) 334 357 ENDDO 335 358 ENDDO 336 !$OMP END PARALLEL 337 ! 338 !-- Horizontal boundary conditions 339 CALL exchange_horiz( r, 1) 340 341 IF ( .NOT. bc_lr_cyc ) THEN 342 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 343 r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l)) 344 ENDIF 345 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 346 r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l)) 347 ENDIF 359 ENDDO 360 !$OMP END PARALLEL 361 ! 362 !-- Horizontal boundary conditions 363 CALL exchange_horiz( r, 1) 364 365 IF ( .NOT. bc_lr_cyc ) THEN 366 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 367 r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l)) 348 368 ENDIF 349 350 IF ( .NOT. bc_ns_cyc ) THEN 351 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 352 r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:) 353 ENDIF 354 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 355 r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:) 356 ENDIF 369 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 370 r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l)) 357 371 ENDIF 358 359 ! 360 !-- Boundary conditions at bottom and top of the domain. Points may be within 361 !-- buildings, but that doesn't matter. 362 IF ( ibc_p_b == 1 ) THEN 363 ! 364 !-- equivalent to r(nzb,:,: ) = r(nzb+1,:,:) 365 r(nzb,:,: ) = r(ind_even_odd+1,:,:) 366 ELSE 367 r(nzb,:,: ) = 0.0_wp 372 ENDIF 373 374 IF ( .NOT. bc_ns_cyc ) THEN 375 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 376 r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:) 368 377 ENDIF 369 370 IF ( ibc_p_t == 1 ) THEN 371 ! 372 !-- equivalent to r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:) 373 r(nzt_mg(l)+1,:,: ) = r(ind_even_odd,:,:) 374 ELSE 375 r(nzt_mg(l)+1,:,: ) = 0.0_wp 378 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 379 r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:) 376 380 ENDIF 377 378 CALL cpu_log( log_point_s(53), 'resid', 'stop' ) 379 380 END SUBROUTINE resid 381 382 383 !------------------------------------------------------------------------------! 381 ENDIF 382 383 ! 384 !-- Boundary conditions at bottom and top of the domain. Points may be within buildings, but that 385 !-- doesn't matter. 386 IF ( ibc_p_b == 1 ) THEN 387 ! 388 !-- Equivalent to r(nzb,:,: ) = r(nzb+1,:,:) 389 r(nzb,:,: ) = r(ind_even_odd+1,:,:) 390 ELSE 391 r(nzb,:,: ) = 0.0_wp 392 ENDIF 393 394 IF ( ibc_p_t == 1 ) THEN 395 ! 396 !-- Equivalent to r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:) 397 r(nzt_mg(l)+1,:,: ) = r(ind_even_odd,:,:) 398 ELSE 399 r(nzt_mg(l)+1,:,: ) = 0.0_wp 400 ENDIF 401 402 CALL cpu_log( log_point_s(53), 'resid', 'stop' ) 403 404 END SUBROUTINE resid 405 406 407 !--------------------------------------------------------------------------------------------------! 384 408 ! Description: 385 409 ! ------------ 386 !> Interpolates the residual on the next coarser grid with "full weighting" 387 !> scheme 388 !------------------------------------------------------------------------------! 389 SUBROUTINE restrict( f_mg, r ) 390 391 392 USE control_parameters, & 393 ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t 394 395 USE indices, & 396 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 397 398 IMPLICIT NONE 399 400 INTEGER(iwp) :: i !< index variable along x on finer grid 401 INTEGER(iwp) :: ic !< index variable along x on coarser grid 402 INTEGER(iwp) :: j !< index variable along y on finer grid 403 INTEGER(iwp) :: jc !< index variable along y on coarser grid 404 INTEGER(iwp) :: k !< index variable along z on finer grid 405 INTEGER(iwp) :: kc !< index variable along z on coarser grid 406 INTEGER(iwp) :: l !< index indicating finer grid level 407 INTEGER(iwp) :: km1 !< index variable along z dimension (k-1 on finer level) 408 INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1 on finer level) 409 410 411 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 412 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 413 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 414 f_mg !< Residual on coarser grid level 415 416 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1, & 417 nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, & 418 nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: & 419 r !< Residual on finer grid level 420 421 ! 422 !-- Interpolate the residual 423 l = grid_level 424 425 CALL cpu_log( log_point_s(54), 'restrict', 'start' ) 426 ! 427 !-- No wall treatment 428 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1) 429 !$OMP DO SCHEDULE( STATIC ) 430 DO ic = nxl_mg(l), nxr_mg(l) 431 i = 2*ic 432 DO jc = nys_mg(l), nyn_mg(l) 433 ! 434 !-- Calculation for the first point along k 435 j = 2*jc 436 ! 437 !-- Calculation for the other points along k 438 !DIR$ IVDEP 439 DO k = ind_even_odd+1, nzt_mg(l+1) ! Fine grid at this point 440 km1 = k-ind_even_odd-1 441 kp1 = k-ind_even_odd 442 kc = k-ind_even_odd ! Coarse grid index 443 444 f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( & 445 8.0_wp * r(k,j,i) & 446 + 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + & 447 r(k,j+1,i) + r(k,j-1,i) ) & 448 + 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + & 449 r(k,j-1,i+1) + r(k,j+1,i+1) ) & 450 + 4.0_wp * r(km1,j,i) & 451 + 2.0_wp * ( r(km1,j,i-1) + r(km1,j,i+1) + & 452 r(km1,j+1,i) + r(km1,j-1,i) ) & 453 + ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + & 454 r(km1,j-1,i+1) + r(km1,j+1,i+1) ) & 455 + 4.0_wp * r(kp1,j,i) & 456 + 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + & 457 r(kp1,j+1,i) + r(kp1,j-1,i) ) & 458 + ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + & 459 r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) & 460 ) 461 ENDDO 410 !> Interpolates the residual on the next coarser grid with "full weighting" scheme 411 !--------------------------------------------------------------------------------------------------! 412 SUBROUTINE restrict( f_mg, r ) 413 414 415 USE control_parameters, & 416 ONLY: bc_lr_cyc, & 417 bc_ns_cyc, & 418 ibc_p_b, & 419 ibc_p_t 420 421 USE indices, & 422 ONLY: nxl_mg, & 423 nxr_mg, & 424 nys_mg, & 425 nyn_mg, & 426 nzb, & 427 nzt_mg 428 429 IMPLICIT NONE 430 431 INTEGER(iwp) :: i !< index variable along x on finer grid 432 INTEGER(iwp) :: ic !< index variable along x on coarser grid 433 INTEGER(iwp) :: j !< index variable along y on finer grid 434 INTEGER(iwp) :: jc !< index variable along y on coarser grid 435 INTEGER(iwp) :: k !< index variable along z on finer grid 436 INTEGER(iwp) :: kc !< index variable along z on coarser grid 437 INTEGER(iwp) :: l !< index indicating finer grid level 438 INTEGER(iwp) :: km1 !< index variable along z dimension (k-1 on finer level) 439 INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1 on finer level) 440 441 442 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 443 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< Residual on coarser grid level 444 445 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, & 446 nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r !< Residual on finer grid level 447 448 ! 449 !-- Interpolate the residual 450 l = grid_level 451 452 CALL cpu_log( log_point_s(54), 'restrict', 'start' ) 453 ! 454 !-- No wall treatment 455 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1) 456 !$OMP DO SCHEDULE( STATIC ) 457 DO ic = nxl_mg(l), nxr_mg(l) 458 i = 2 * ic 459 DO jc = nys_mg(l), nyn_mg(l) 460 ! 461 !-- Calculation for the first point along k 462 j = 2 * jc 463 ! 464 !-- Calculation for the other points along k 465 !DIR$ IVDEP 466 DO k = ind_even_odd+1, nzt_mg(l+1) ! Fine grid at this point 467 km1 = k-ind_even_odd-1 468 kp1 = k-ind_even_odd 469 kc = k-ind_even_odd ! Coarse grid index 470 471 f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp & 472 * ( 8.0_wp * r(k,j,i) + 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) & 473 + r(k,j+1,i) + r(k,j-1,i) ) & 474 + 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) & 475 + r(k,j-1,i+1) + r(k,j+1,i+1) ) & 476 + 4.0_wp * r(km1,j,i) & 477 + 2.0_wp * ( r(km1,j,i-1) + r(km1,j,i+1) & 478 + r(km1,j+1,i) + r(km1,j-1,i) ) & 479 + ( r(km1,j-1,i-1) + r(km1,j+1,i-1) & 480 + r(km1,j-1,i+1) + r(km1,j+1,i+1) ) & 481 + 4.0_wp * r(kp1,j,i) & 482 + 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) & 483 + r(kp1,j+1,i) + r(kp1,j-1,i) ) & 484 + ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) & 485 + r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) & 486 ) 462 487 ENDDO 463 488 ENDDO 464 !$OMP ENDDO 465 !$OMP END PARALLEL 466 467 ! 468 !-- Ghost point exchange 469 CALL exchange_horiz( f_mg, 1) 470 ! 471 !-- Horizontal boundary conditions 472 IF ( .NOT. bc_lr_cyc ) THEN 473 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 474 f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l)) 475 ENDIF 476 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 477 f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l)) 478 ENDIF 489 ENDDO 490 !$OMP ENDDO 491 !$OMP END PARALLEL 492 493 ! 494 !-- Ghost point exchange 495 CALL exchange_horiz( f_mg, 1) 496 ! 497 !-- Horizontal boundary conditions 498 IF ( .NOT. bc_lr_cyc ) THEN 499 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 500 f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l)) 479 501 ENDIF 480 481 IF ( .NOT. bc_ns_cyc ) THEN 482 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 483 f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:) 484 ENDIF 485 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 486 f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:) 487 ENDIF 502 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 503 f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l)) 488 504 ENDIF 489 490 ! 491 !-- Boundary conditions at bottom and top of the domain. 492 !-- These points are not handled by the above loop. Points may be within 493 !-- buildings, but that doesn't matter. Remark: f_mg is ordered sequentielly 494 !-- after interpolation on coarse grid (is ordered in odd-even blocks further 495 !-- below). 496 IF ( ibc_p_b == 1 ) THEN 497 f_mg(nzb,:,: ) = f_mg(nzb+1,:,:) 498 ELSE 499 f_mg(nzb,:,: ) = 0.0_wp 505 ENDIF 506 507 IF ( .NOT. bc_ns_cyc ) THEN 508 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 509 f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:) 500 510 ENDIF 501 502 IF ( ibc_p_t == 1 ) THEN 503 f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:) 504 ELSE 505 f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp 511 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 512 f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:) 506 513 ENDIF 507 508 CALL cpu_log( log_point_s(54), 'restrict', 'stop' ) 509 ! 510 !-- Since residual is in sequential order after interpolation, an additional 511 !-- sorting in odd-even blocks along z dimension is required at this point. 512 CALL sort_k_to_even_odd_blocks( f_mg , l) 513 514 END SUBROUTINE restrict 515 516 517 !------------------------------------------------------------------------------! 514 ENDIF 515 516 ! 517 !-- Boundary conditions at bottom and top of the domain. These points are not handled by the above 518 !-- loop. Points may be within buildings, but that doesn't matter. Remark: f_mg is ordered 519 !-- sequentielly after interpolation on coarse grid (is ordered in odd-even blocks further below). 520 IF ( ibc_p_b == 1 ) THEN 521 f_mg(nzb,:,: ) = f_mg(nzb+1,:,:) 522 ELSE 523 f_mg(nzb,:,: ) = 0.0_wp 524 ENDIF 525 526 IF ( ibc_p_t == 1 ) THEN 527 f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:) 528 ELSE 529 f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp 530 ENDIF 531 532 CALL cpu_log( log_point_s(54), 'restrict', 'stop' ) 533 ! 534 !-- Since residual is in sequential order after interpolation, an additional sorting in odd-even 535 !-- blocks along z dimension is required at this point. 536 CALL sort_k_to_even_odd_blocks( f_mg , l) 537 538 END SUBROUTINE restrict 539 540 541 !--------------------------------------------------------------------------------------------------! 518 542 ! Description: 519 543 ! ------------ 520 !> Interpolates the correction of the perturbation pressure 521 !> to the next finer grid. 522 !------------------------------------------------------------------------------! 523 SUBROUTINE prolong( p, temp ) 524 525 526 USE control_parameters, & 527 ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t 528 USE indices, & 529 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 530 531 IMPLICIT NONE 532 533 INTEGER(iwp) :: i !< index variable along x on coarser grid level 534 INTEGER(iwp) :: j !< index variable along y on coarser grid level 535 INTEGER(iwp) :: k !< index variable along z on coarser grid level 536 INTEGER(iwp) :: l !< index indicating finer grid level 537 INTEGER(iwp) :: kp1 !< index variable along z 538 INTEGER(iwp) :: ke !< index for prolog even 539 INTEGER(iwp) :: ko !< index for prolog odd 540 541 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, & 542 nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 543 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: & 544 p !< perturbation pressure on coarser grid level 545 546 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 547 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 548 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 549 temp !< perturbation pressure on finer grid level 550 551 552 CALL cpu_log( log_point_s(55), 'prolong', 'start' ) 553 554 ! 555 !-- First, store elements of the coarser grid on the next finer grid 556 l = grid_level 557 ind_even_odd = even_odd_level(grid_level-1) 558 559 !$OMP PARALLEL PRIVATE (i,j,k,kp1,ke,ko) 560 !$OMP DO 561 DO i = nxl_mg(l-1), nxr_mg(l-1) 562 DO j = nys_mg(l-1), nyn_mg(l-1) 563 564 !DIR$ IVDEP 565 DO k = ind_even_odd+1, nzt_mg(l-1) 566 kp1 = k - ind_even_odd 567 ke = 2 * ( k-ind_even_odd - 1 ) + 1 568 ko = 2 * k - 1 569 ! 570 !-- Points of the coarse grid are directly stored on the next finer 571 !-- grid 572 temp(ko,2*j,2*i) = p(k,j,i) 573 ! 574 !-- Points between two coarse-grid points 575 temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) ) 576 temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) ) 577 temp(ke,2*j,2*i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) ) 578 ! 579 !-- Points in the center of the planes stretched by four points 580 !-- of the coarse grid cube 581 temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 582 p(k,j+1,i) + p(k,j+1,i+1) ) 583 temp(ke,2*j,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 584 p(kp1,j,i) + p(kp1,j,i+1) ) 585 temp(ke,2*j+1,2*i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + & 586 p(kp1,j,i) + p(kp1,j+1,i) ) 587 ! 588 !-- Points in the middle of coarse grid cube 589 temp(ke,2*j+1,2*i+1) = 0.125_wp * & 590 ( p(k,j,i) + p(k,j,i+1) + & 591 p(k,j+1,i) + p(k,j+1,i+1) + & 592 p(kp1,j,i) + p(kp1,j,i+1) + & 593 p(kp1,j+1,i) + p(kp1,j+1,i+1) ) 594 595 ENDDO 596 597 !DIR$ IVDEP 598 DO k = nzb+1, ind_even_odd 599 kp1 = k + ind_even_odd + 1 600 ke = 2 * k 601 ko = 2 * ( k + ind_even_odd ) 602 ! 603 !-- Points of the coarse grid are directly stored on the next finer 604 !-- grid 605 temp(ko,2*j,2*i) = p(k,j,i) 606 ! 607 !-- Points between two coarse-grid points 608 temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) ) 609 temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) ) 610 temp(ke,2*j,2*i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) ) 611 ! 612 !-- Points in the center of the planes stretched by four points 613 !-- of the coarse grid cube 614 temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 615 p(k,j+1,i) + p(k,j+1,i+1) ) 616 temp(ke,2*j,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 617 p(kp1,j,i) + p(kp1,j,i+1) ) 618 temp(ke,2*j+1,2*i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + & 619 p(kp1,j,i) + p(kp1,j+1,i) ) 620 ! 621 !-- Points in the middle of coarse grid cube 622 temp(ke,2*j+1,2*i+1) = 0.125_wp * & 623 ( p(k,j,i) + p(k,j,i+1) + & 624 p(k,j+1,i) + p(k,j+1,i+1) + & 625 p(kp1,j,i) + p(kp1,j,i+1) + & 626 p(kp1,j+1,i) + p(kp1,j+1,i+1) ) 627 628 ENDDO 629 630 ENDDO 544 !> Interpolates the correction of the perturbation pressure to the next finer grid. 545 !--------------------------------------------------------------------------------------------------! 546 SUBROUTINE prolong( p, temp ) 547 548 549 USE control_parameters, & 550 ONLY: bc_lr_cyc, & 551 bc_ns_cyc, & 552 ibc_p_b, & 553 ibc_p_t 554 USE indices, & 555 ONLY: nxl_mg, & 556 nxr_mg, & 557 nys_mg, & 558 nyn_mg, & 559 nzb, & 560 nzt_mg 561 562 IMPLICIT NONE 563 564 INTEGER(iwp) :: i !< index variable along x on coarser grid level 565 INTEGER(iwp) :: j !< index variable along y on coarser grid level 566 INTEGER(iwp) :: k !< index variable along z on coarser grid level 567 INTEGER(iwp) :: l !< index indicating finer grid level 568 INTEGER(iwp) :: kp1 !< index variable along z 569 INTEGER(iwp) :: ke !< index for prolog even 570 INTEGER(iwp) :: ko !< index for prolog odd 571 572 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 573 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p !< perturbation pressure on coarser grid level 574 575 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 576 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp !< perturbation pressure on finer grid level 577 578 579 CALL cpu_log( log_point_s(55), 'prolong', 'start' ) 580 581 ! 582 !-- First, store elements of the coarser grid on the next finer grid 583 l = grid_level 584 ind_even_odd = even_odd_level(grid_level-1) 585 586 !$OMP PARALLEL PRIVATE (i,j,k,kp1,ke,ko) 587 !$OMP DO 588 DO i = nxl_mg(l-1), nxr_mg(l-1) 589 DO j = nys_mg(l-1), nyn_mg(l-1) 590 591 !DIR$ IVDEP 592 DO k = ind_even_odd+1, nzt_mg(l-1) 593 kp1 = k - ind_even_odd 594 ke = 2 * ( k-ind_even_odd - 1 ) + 1 595 ko = 2 * k - 1 596 ! 597 !-- Points of the coarse grid are directly stored on the next finer grid 598 temp(ko,2*j,2*i) = p(k,j,i) 599 ! 600 !-- Points between two coarse-grid points 601 temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) ) 602 temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) ) 603 temp(ke,2*j,2*i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) ) 604 ! 605 !-- Points in the center of the planes stretched by four points of the coarse grid cube 606 temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 607 p(k,j+1,i) + p(k,j+1,i+1) ) 608 temp(ke,2*j,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 609 p(kp1,j,i) + p(kp1,j,i+1) ) 610 temp(ke,2*j+1,2*i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + & 611 p(kp1,j,i) + p(kp1,j+1,i) ) 612 ! 613 !-- Points in the middle of coarse grid cube 614 temp(ke,2*j+1,2*i+1) = 0.125_wp * ( p(k,j,i) + p(k,j,i+1) + p(k,j+1,i) & 615 + p(k,j+1,i+1) + p(kp1,j,i) + p(kp1,j,i+1) & 616 + p(kp1,j+1,i) + p(kp1,j+1,i+1) ) 617 618 ENDDO 619 620 !DIR$ IVDEP 621 DO k = nzb+1, ind_even_odd 622 kp1 = k + ind_even_odd + 1 623 ke = 2 * k 624 ko = 2 * ( k + ind_even_odd ) 625 ! 626 !-- Points of the coarse grid are directly stored on the next finer grid 627 temp(ko,2*j,2*i) = p(k,j,i) 628 ! 629 !-- Points between two coarse-grid points 630 temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) ) 631 temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) ) 632 temp(ke,2*j,2*i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) ) 633 ! 634 !-- Points in the center of the planes stretched by four points 635 !-- of the coarse grid cube 636 temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 637 p(k,j+1,i) + p(k,j+1,i+1) ) 638 temp(ke,2*j,2*i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + & 639 p(kp1,j,i) + p(kp1,j,i+1) ) 640 temp(ke,2*j+1,2*i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + & 641 p(kp1,j,i) + p(kp1,j+1,i) ) 642 ! 643 !-- Points in the middle of coarse grid cube 644 temp(ke,2*j+1,2*i+1) = 0.125_wp * ( p(k,j,i) + p(k,j,i+1) + p(k,j+1,i) & 645 + p(k,j+1,i+1) + p(kp1,j,i) + p(kp1,j,i+1) & 646 + p(kp1,j+1,i) + p(kp1,j+1,i+1) ) 647 648 ENDDO 649 631 650 ENDDO 632 !$OMP END PARALLEL 633 634 ind_even_odd = even_odd_level(grid_level) 635 ! 636 !-- Horizontal boundary conditions 637 CALL exchange_horiz( temp, 1) 638 639 IF ( .NOT. bc_lr_cyc ) THEN 640 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 641 temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l)) 642 ENDIF 643 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 644 temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l)) 645 ENDIF 651 ENDDO 652 !$OMP END PARALLEL 653 654 ind_even_odd = even_odd_level(grid_level) 655 ! 656 !-- Horizontal boundary conditions 657 CALL exchange_horiz( temp, 1) 658 659 IF ( .NOT. bc_lr_cyc ) THEN 660 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 661 temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l)) 646 662 ENDIF 647 648 IF ( .NOT. bc_ns_cyc ) THEN 649 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 650 temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:) 651 ENDIF 652 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 653 temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:) 654 ENDIF 663 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 664 temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l)) 655 665 ENDIF 656 657 ! 658 !-- Bottom and top boundary conditions 659 IF ( ibc_p_b == 1 ) THEN 660 ! 661 !-- equivalent to temp(nzb,:,: ) = temp(nzb+1,:,:) 662 temp(nzb,:,: ) = temp(ind_even_odd+1,:,:) 663 ELSE 664 temp(nzb,:,: ) = 0.0_wp 666 ENDIF 667 668 IF ( .NOT. bc_ns_cyc ) THEN 669 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 670 temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:) 665 671 ENDIF 666 667 IF ( ibc_p_t == 1 ) THEN 668 ! 669 !-- equivalent to temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:) 670 temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:) 671 ELSE 672 temp(nzt_mg(l)+1,:,: ) = 0.0_wp 672 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 673 temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:) 673 674 ENDIF 674 675 CALL cpu_log( log_point_s(55), 'prolong', 'stop' ) 676 677 END SUBROUTINE prolong 678 679 680 !------------------------------------------------------------------------------! 675 ENDIF 676 677 ! 678 !-- Bottom and top boundary conditions 679 IF ( ibc_p_b == 1 ) THEN 680 ! 681 !-- Equivalent to temp(nzb,:,: ) = temp(nzb+1,:,:) 682 temp(nzb,:,: ) = temp(ind_even_odd+1,:,:) 683 ELSE 684 temp(nzb,:,: ) = 0.0_wp 685 ENDIF 686 687 IF ( ibc_p_t == 1 ) THEN 688 ! 689 !-- Equivalent to temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:) 690 temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:) 691 ELSE 692 temp(nzt_mg(l)+1,:,: ) = 0.0_wp 693 ENDIF 694 695 CALL cpu_log( log_point_s(55), 'prolong', 'stop' ) 696 697 END SUBROUTINE prolong 698 699 700 !--------------------------------------------------------------------------------------------------! 681 701 ! Description: 682 702 ! ------------ 683 !> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with 684 !> 3D-Red-Black decomposition (GS-RB) is used. 685 !------------------------------------------------------------------------------! 686 SUBROUTINE redblack( f_mg, p_mg ) 687 688 689 USE arrays_3d, & 690 ONLY: rho_air_mg 691 692 USE control_parameters, & 693 ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t, ngsrb 694 695 USE grid_variables, & 696 ONLY: ddx2_mg, ddy2_mg 697 698 USE indices, & 699 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 700 701 IMPLICIT NONE 702 703 INTEGER(iwp) :: color !< grid point color, either red or black 704 INTEGER(iwp) :: i !< index variable along x 705 INTEGER(iwp) :: ic !< index variable along x 706 INTEGER(iwp) :: j !< index variable along y 707 INTEGER(iwp) :: jc !< index variable along y 708 INTEGER(iwp) :: jj !< index variable along y 709 INTEGER(iwp) :: k !< index variable along z 710 INTEGER(iwp) :: l !< grid level 711 INTEGER(iwp) :: n !< loop variable GauÃ-Seidel iterations 712 INTEGER(iwp) :: km1 !< index variable (k-1) 713 INTEGER(iwp) :: kp1 !< index variable (k+1) 714 715 LOGICAL :: unroll !< flag indicating whether loop unrolling is possible 716 717 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 718 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 719 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 720 f_mg !< residual of perturbation pressure 721 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 722 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 723 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 724 p_mg !< perturbation pressure 725 726 l = grid_level 727 728 unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. & 729 MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 ) 730 731 DO n = 1, ngsrb 732 733 DO color = 1, 2 734 735 IF ( .NOT. unroll ) THEN 736 737 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' ) 738 ! 739 !-- Without unrolling of loops, no cache optimization 740 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 741 !$OMP DO 742 DO i = nxl_mg(l), nxr_mg(l), 2 743 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 744 !DIR$ IVDEP 745 DO k = ind_even_odd+1, nzt_mg(l) 746 km1 = k-ind_even_odd-1 747 kp1 = k-ind_even_odd 748 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 749 rho_air_mg(k,l) * ddx2_mg(l) * & 750 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 751 + rho_air_mg(k,l) * ddy2_mg(l) * & 752 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 753 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 754 + f3_mg_b(k,l) * p_mg(km1,j,i) & 755 - f_mg(k,j,i) ) 756 ENDDO 703 !> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with 3D-Red-Black 704 !> decomposition (GS-RB) is used. 705 !--------------------------------------------------------------------------------------------------! 706 SUBROUTINE redblack( f_mg, p_mg ) 707 708 709 USE arrays_3d, & 710 ONLY: rho_air_mg 711 712 USE control_parameters, & 713 ONLY: bc_lr_cyc, & 714 bc_ns_cyc, & 715 ibc_p_b, & 716 ibc_p_t, & 717 ngsrb 718 719 USE grid_variables, & 720 ONLY: ddx2_mg, & 721 ddy2_mg 722 723 USE indices, & 724 ONLY: nxl_mg, & 725 nxr_mg, & 726 nys_mg, & 727 nyn_mg, & 728 nzb, & 729 nzt_mg 730 731 IMPLICIT NONE 732 733 INTEGER(iwp) :: color !< grid point color, either red or black 734 INTEGER(iwp) :: i !< index variable along x 735 INTEGER(iwp) :: ic !< index variable along x 736 INTEGER(iwp) :: j !< index variable along y 737 INTEGER(iwp) :: jc !< index variable along y 738 INTEGER(iwp) :: jj !< index variable along y 739 INTEGER(iwp) :: k !< index variable along z 740 INTEGER(iwp) :: km1 !< index variable (k-1) 741 INTEGER(iwp) :: kp1 !< index variable (k+1) 742 INTEGER(iwp) :: l !< grid level 743 INTEGER(iwp) :: n !< loop variable GauÃ-Seidel iterations 744 745 746 LOGICAL :: unroll !< flag indicating whether loop unrolling is possible 747 748 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 749 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< residual of perturbation pressure 750 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 751 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< perturbation pressure 752 753 l = grid_level 754 755 unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 ) 756 757 DO n = 1, ngsrb 758 759 DO color = 1, 2 760 761 IF ( .NOT. unroll ) THEN 762 763 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' ) 764 ! 765 !-- Without unrolling of loops, no cache optimization 766 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 767 !$OMP DO 768 DO i = nxl_mg(l), nxr_mg(l), 2 769 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 770 !DIR$ IVDEP 771 DO k = ind_even_odd+1, nzt_mg(l) 772 km1 = k-ind_even_odd-1 773 kp1 = k-ind_even_odd 774 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) & 775 * ( rho_air_mg(k,l) * ddx2_mg(l) * & 776 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 777 + rho_air_mg(k,l) * ddy2_mg(l) * & 778 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 779 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 780 + f3_mg_b(k,l) * p_mg(km1,j,i) - f_mg(k,j,i) & 781 ) 757 782 ENDDO 758 783 ENDDO 759 760 !$OMP DO 761 DO i = nxl_mg(l)+1, nxr_mg(l), 2762 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2763 !DIR$ IVDEP764 DO k = ind_even_odd+1, nzt_mg(l)765 km1 = k-ind_even_odd-1766 kp1 = k-ind_even_odd767 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &768 rho_air_mg(k,l) * ddx2_mg(l) *&769 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )&770 + rho_air_mg(k,l) * ddy2_mg(l) *&771 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )&772 + f2_mg_b(k,l) * p_mg(kp1,j,i)&773 + f3_mg_b(k,l) * p_mg(km1,j,i)&774 - f_mg(k,j,i) )775 ENDDO784 ENDDO 785 786 !$OMP DO 787 DO i = nxl_mg(l)+1, nxr_mg(l), 2 788 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 789 !DIR$ IVDEP 790 DO k = ind_even_odd+1, nzt_mg(l) 791 km1 = k-ind_even_odd-1 792 kp1 = k-ind_even_odd 793 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 794 rho_air_mg(k,l) * ddx2_mg(l) * & 795 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 796 + rho_air_mg(k,l) * ddy2_mg(l) * & 797 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 798 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 799 + f3_mg_b(k,l) * p_mg(km1,j,i) & 800 - f_mg(k,j,i) ) 776 801 ENDDO 777 802 ENDDO 778 779 !$OMP DO 780 DO i = nxl_mg(l), nxr_mg(l), 2781 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2782 !DIR$ IVDEP783 DO k = nzb+1, ind_even_odd784 km1 = k+ind_even_odd785 kp1 = k+ind_even_odd+1786 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &787 rho_air_mg(k,l) * ddx2_mg(l) *&788 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )&789 + rho_air_mg(k,l) * ddy2_mg(l) *&790 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )&791 + f2_mg_b(k,l) * p_mg(kp1,j,i)&792 + f3_mg_b(k,l) * p_mg(km1,j,i)&793 - f_mg(k,j,i) )794 ENDDO803 ENDDO 804 805 !$OMP DO 806 DO i = nxl_mg(l), nxr_mg(l), 2 807 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 808 !DIR$ IVDEP 809 DO k = nzb+1, ind_even_odd 810 km1 = k+ind_even_odd 811 kp1 = k+ind_even_odd+1 812 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 813 rho_air_mg(k,l) * ddx2_mg(l) * & 814 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 815 + rho_air_mg(k,l) * ddy2_mg(l) * & 816 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 817 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 818 + f3_mg_b(k,l) * p_mg(km1,j,i) & 819 - f_mg(k,j,i) ) 795 820 ENDDO 796 821 ENDDO 797 798 !$OMP DO 799 DO i = nxl_mg(l)+1, nxr_mg(l), 2800 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2801 !DIR$ IVDEP802 DO k = nzb+1, ind_even_odd803 km1 = k+ind_even_odd804 kp1 = k+ind_even_odd+1805 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &806 rho_air_mg(k,l) * ddx2_mg(l) *&807 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )&808 + rho_air_mg(k,l) * ddy2_mg(l) *&809 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )&810 + f2_mg_b(k,l) * p_mg(kp1,j,i)&811 + f3_mg_b(k,l) * p_mg(km1,j,i)&812 - f_mg(k,j,i) )813 ENDDO822 ENDDO 823 824 !$OMP DO 825 DO i = nxl_mg(l)+1, nxr_mg(l), 2 826 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 827 !DIR$ IVDEP 828 DO k = nzb+1, ind_even_odd 829 km1 = k+ind_even_odd 830 kp1 = k+ind_even_odd+1 831 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 832 rho_air_mg(k,l) * ddx2_mg(l) * & 833 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 834 + rho_air_mg(k,l) * ddy2_mg(l) * & 835 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 836 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 837 + f3_mg_b(k,l) * p_mg(km1,j,i) & 838 - f_mg(k,j,i) ) 814 839 ENDDO 815 840 ENDDO 816 !$OMP END PARALLEL 817 818 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' ) 819 820 ELSE 821 ! 822 !-- Loop unrolling along y, only one i loop for better cache use 823 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' ) 824 825 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj) 826 !$OMP DO 827 DO ic = nxl_mg(l), nxr_mg(l), 2 828 DO jc = nys_mg(l), nyn_mg(l), 4 829 i = ic 830 jj = jc+2-color 831 !DIR$ IVDEP 832 DO k = ind_even_odd+1, nzt_mg(l) 833 km1 = k-ind_even_odd-1 834 kp1 = k-ind_even_odd 835 j = jj 836 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 837 rho_air_mg(k,l) * ddx2_mg(l) * & 838 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 839 + rho_air_mg(k,l) * ddy2_mg(l) * & 840 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 841 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 842 + f3_mg_b(k,l) * p_mg(km1,j,i) & 843 - f_mg(k,j,i) ) 844 j = jj+2 845 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 846 rho_air_mg(k,l) * ddx2_mg(l) * & 847 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 848 + rho_air_mg(k,l) * ddy2_mg(l) * & 849 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 850 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 851 + f3_mg_b(k,l) * p_mg(km1,j,i) & 852 - f_mg(k,j,i) ) 853 ENDDO 854 855 i = ic+1 856 jj = jc+color-1 857 !DIR$ IVDEP 858 DO k = ind_even_odd+1, nzt_mg(l) 859 km1 = k-ind_even_odd-1 860 kp1 = k-ind_even_odd 861 j = jj 862 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 863 rho_air_mg(k,l) * ddx2_mg(l) * & 864 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 865 + rho_air_mg(k,l) * ddy2_mg(l) * & 866 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 867 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 868 + f3_mg_b(k,l) * p_mg(km1,j,i) & 869 - f_mg(k,j,i) ) 870 j = jj+2 871 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 872 rho_air_mg(k,l) * ddx2_mg(l) * & 873 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 874 + rho_air_mg(k,l) * ddy2_mg(l) * & 875 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 876 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 877 + f3_mg_b(k,l) * p_mg(km1,j,i) & 878 - f_mg(k,j,i) ) 879 ENDDO 880 881 i = ic 882 jj = jc+color-1 883 !DIR$ IVDEP 884 DO k = nzb+1, ind_even_odd 885 km1 = k+ind_even_odd 886 kp1 = k+ind_even_odd+1 887 j = jj 888 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 889 rho_air_mg(k,l) * ddx2_mg(l) * & 890 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 891 + rho_air_mg(k,l) * ddy2_mg(l) * & 892 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 893 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 894 + f3_mg_b(k,l) * p_mg(km1,j,i) & 895 - f_mg(k,j,i) ) 896 j = jj+2 897 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 898 rho_air_mg(k,l) * ddx2_mg(l) * & 899 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 900 + rho_air_mg(k,l) * ddy2_mg(l) * & 901 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 902 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 903 + f3_mg_b(k,l) * p_mg(km1,j,i) & 904 - f_mg(k,j,i) ) 905 ENDDO 906 907 i = ic+1 908 jj = jc+2-color 909 !DIR$ IVDEP 910 DO k = nzb+1, ind_even_odd 911 km1 = k+ind_even_odd 912 kp1 = k+ind_even_odd+1 913 j = jj 914 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 915 rho_air_mg(k,l) * ddx2_mg(l) * & 916 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 917 + rho_air_mg(k,l) * ddy2_mg(l) * & 918 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 919 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 920 + f3_mg_b(k,l) * p_mg(km1,j,i) & 921 - f_mg(k,j,i) ) 922 j = jj+2 923 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 924 rho_air_mg(k,l) * ddx2_mg(l) * & 925 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 926 + rho_air_mg(k,l) * ddy2_mg(l) * & 927 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 928 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 929 + f3_mg_b(k,l) * p_mg(km1,j,i) & 930 - f_mg(k,j,i) ) 931 ENDDO 932 841 ENDDO 842 !$OMP END PARALLEL 843 844 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' ) 845 846 ELSE 847 ! 848 !-- Loop unrolling along y, only one i loop for better cache use 849 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' ) 850 851 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj) 852 !$OMP DO 853 DO ic = nxl_mg(l), nxr_mg(l), 2 854 DO jc = nys_mg(l), nyn_mg(l), 4 855 i = ic 856 jj = jc+2-color 857 !DIR$ IVDEP 858 DO k = ind_even_odd+1, nzt_mg(l) 859 km1 = k-ind_even_odd-1 860 kp1 = k-ind_even_odd 861 j = jj 862 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 863 rho_air_mg(k,l) * ddx2_mg(l) * & 864 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 865 + rho_air_mg(k,l) * ddy2_mg(l) * & 866 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 867 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 868 + f3_mg_b(k,l) * p_mg(km1,j,i) & 869 - f_mg(k,j,i) ) 870 871 j = jj+2 872 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 873 rho_air_mg(k,l) * ddx2_mg(l) * & 874 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 875 + rho_air_mg(k,l) * ddy2_mg(l) * & 876 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 877 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 878 + f3_mg_b(k,l) * p_mg(km1,j,i) & 879 - f_mg(k,j,i) ) 933 880 ENDDO 934 ENDDO 935 !$OMP END PARALLEL 936 937 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' ) 938 939 ENDIF 940 941 ! 942 !-- Horizontal boundary conditions 943 CALL special_exchange_horiz( p_mg, color ) 944 945 IF ( .NOT. bc_lr_cyc ) THEN 946 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 947 p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l)) 948 ENDIF 949 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 950 p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l)) 951 ENDIF 952 ENDIF 953 954 IF ( .NOT. bc_ns_cyc ) THEN 955 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 956 p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:) 957 ENDIF 958 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 959 p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:) 960 ENDIF 961 ENDIF 962 963 ! 964 !-- Bottom and top boundary conditions 965 IF ( ibc_p_b == 1 ) THEN 966 ! 967 !-- equivalent to p_mg(nzb,:,: ) = p_mg(nzb+1,:,:) 968 p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:) 969 ELSE 970 p_mg(nzb,:,: ) = 0.0_wp 971 ENDIF 972 973 IF ( ibc_p_t == 1 ) THEN 974 ! 975 !-- equivalent to p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:) 976 p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:) 977 ELSE 978 p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp 979 ENDIF 980 981 ENDDO 881 882 i = ic+1 883 jj = jc+color-1 884 !DIR$ IVDEP 885 DO k = ind_even_odd+1, nzt_mg(l) 886 km1 = k-ind_even_odd-1 887 kp1 = k-ind_even_odd 888 j = jj 889 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 890 rho_air_mg(k,l) * ddx2_mg(l) * & 891 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 892 + rho_air_mg(k,l) * ddy2_mg(l) * & 893 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 894 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 895 + f3_mg_b(k,l) * p_mg(km1,j,i) & 896 - f_mg(k,j,i) ) 897 898 j = jj+2 899 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 900 rho_air_mg(k,l) * ddx2_mg(l) * & 901 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 902 + rho_air_mg(k,l) * ddy2_mg(l) * & 903 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 904 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 905 + f3_mg_b(k,l) * p_mg(km1,j,i) & 906 - f_mg(k,j,i) ) 907 ENDDO 908 909 i = ic 910 jj = jc+color-1 911 !DIR$ IVDEP 912 DO k = nzb+1, ind_even_odd 913 km1 = k+ind_even_odd 914 kp1 = k+ind_even_odd+1 915 j = jj 916 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 917 rho_air_mg(k,l) * ddx2_mg(l) * & 918 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 919 + rho_air_mg(k,l) * ddy2_mg(l) * & 920 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 921 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 922 + f3_mg_b(k,l) * p_mg(km1,j,i) & 923 - f_mg(k,j,i) ) 924 925 j = jj+2 926 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 927 rho_air_mg(k,l) * ddx2_mg(l) * & 928 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 929 + rho_air_mg(k,l) * ddy2_mg(l) * & 930 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 931 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 932 + f3_mg_b(k,l) * p_mg(km1,j,i) & 933 - f_mg(k,j,i) ) 934 ENDDO 935 936 i = ic+1 937 jj = jc+2-color 938 !DIR$ IVDEP 939 DO k = nzb+1, ind_even_odd 940 km1 = k+ind_even_odd 941 kp1 = k+ind_even_odd+1 942 j = jj 943 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 944 rho_air_mg(k,l) * ddx2_mg(l) * & 945 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 946 + rho_air_mg(k,l) * ddy2_mg(l) * & 947 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 948 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 949 + f3_mg_b(k,l) * p_mg(km1,j,i) & 950 - f_mg(k,j,i) ) 951 952 j = jj+2 953 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 954 rho_air_mg(k,l) * ddx2_mg(l) * & 955 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 956 + rho_air_mg(k,l) * ddy2_mg(l) * & 957 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 958 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 959 + f3_mg_b(k,l) * p_mg(km1,j,i) & 960 - f_mg(k,j,i) ) 961 ENDDO 962 963 ENDDO 964 ENDDO 965 !$OMP END PARALLEL 966 967 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' ) 968 969 ENDIF 970 971 ! 972 !-- Horizontal boundary conditions 973 CALL special_exchange_horiz( p_mg, color ) 974 975 IF ( .NOT. bc_lr_cyc ) THEN 976 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 977 p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l)) 978 ENDIF 979 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 980 p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l)) 981 ENDIF 982 ENDIF 983 984 IF ( .NOT. bc_ns_cyc ) THEN 985 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 986 p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:) 987 ENDIF 988 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 989 p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:) 990 ENDIF 991 ENDIF 992 993 ! 994 !-- Bottom and top boundary conditions 995 IF ( ibc_p_b == 1 ) THEN 996 ! 997 !-- Equivalent to p_mg(nzb,:,: ) = p_mg(nzb+1,:,:) 998 p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:) 999 ELSE 1000 p_mg(nzb,:,: ) = 0.0_wp 1001 ENDIF 1002 1003 IF ( ibc_p_t == 1 ) THEN 1004 ! 1005 !-- Equivalent to p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:) 1006 p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:) 1007 ELSE 1008 p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp 1009 ENDIF 982 1010 983 1011 ENDDO 984 1012 985 END SUBROUTINE redblack 986 987 988 !------------------------------------------------------------------------------! 1013 ENDDO 1014 1015 END SUBROUTINE redblack 1016 1017 1018 !--------------------------------------------------------------------------------------------------! 989 1019 ! Description: 990 1020 ! ------------ 991 !> Sort k-Dimension from sequential into blocks of even and odd. 992 !> This is required to vectorize the red-black subroutine. 993 !> Version for 3D-REAL arrays 994 !------------------------------------------------------------------------------! 995 SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel ) 996 997 998 USE indices, & 999 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1000 1001 IMPLICIT NONE 1002 1003 INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1004 1005 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1, & 1006 nys_mg(glevel)-1:nyn_mg(glevel)+1, & 1007 nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: & 1008 p_mg !< array to be sorted 1009 ! 1010 !-- Local variables 1011 INTEGER(iwp) :: i !< index variable along x 1012 INTEGER(iwp) :: j !< index variable along y 1013 INTEGER(iwp) :: k !< index variable along z 1014 INTEGER(iwp) :: l !< grid level 1015 INTEGER(iwp) :: ind !< index variable along z 1016 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< odd-even sorted temporary array 1017 1018 1019 CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' ) 1020 1021 l = glevel 1022 ind_even_odd = even_odd_level(l) 1023 1024 !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp) 1025 !$OMP DO 1026 DO i = nxl_mg(l)-1, nxr_mg(l)+1 1027 DO j = nys_mg(l)-1, nyn_mg(l)+1 1028 1029 ! 1030 !-- Sort the data with even k index 1031 ind = nzb-1 1032 DO k = nzb, nzt_mg(l), 2 1033 ind = ind + 1 1034 tmp(ind) = p_mg(k,j,i) 1035 ENDDO 1036 ! 1037 !-- Sort the data with odd k index 1038 DO k = nzb+1, nzt_mg(l)+1, 2 1039 ind = ind + 1 1040 tmp(ind) = p_mg(k,j,i) 1041 ENDDO 1042 1043 p_mg(:,j,i) = tmp 1044 1045 ENDDO 1021 !> Sort k-dimension from sequential into blocks of even and odd. This is required to vectorize the 1022 !> red-black subroutine. Version for 3D-REAL arrays 1023 !--------------------------------------------------------------------------------------------------! 1024 SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel ) 1025 1026 1027 USE indices, & 1028 ONLY: nxl_mg, & 1029 nxr_mg, & 1030 nys_mg, & 1031 nyn_mg, & 1032 nzb, & 1033 nzt_mg 1034 1035 IMPLICIT NONE 1036 1037 INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1038 1039 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1,nys_mg(glevel)-1:nyn_mg(glevel)+1, & 1040 nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: p_mg !< array to be sorted 1041 ! 1042 !-- Local variables 1043 INTEGER(iwp) :: i !< index variable along x 1044 INTEGER(iwp) :: ind !< index variable along z 1045 INTEGER(iwp) :: j !< index variable along y 1046 INTEGER(iwp) :: k !< index variable along z 1047 INTEGER(iwp) :: l !< grid level 1048 1049 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< odd-even sorted temporary array 1050 1051 1052 CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' ) 1053 1054 l = glevel 1055 ind_even_odd = even_odd_level(l) 1056 1057 !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp) 1058 !$OMP DO 1059 DO i = nxl_mg(l)-1, nxr_mg(l)+1 1060 DO j = nys_mg(l)-1, nyn_mg(l)+1 1061 1062 ! 1063 !-- Sort the data with even k index 1064 ind = nzb-1 1065 DO k = nzb, nzt_mg(l), 2 1066 ind = ind + 1 1067 tmp(ind) = p_mg(k,j,i) 1068 ENDDO 1069 ! 1070 !-- Sort the data with odd k index 1071 DO k = nzb+1, nzt_mg(l)+1, 2 1072 ind = ind + 1 1073 tmp(ind) = p_mg(k,j,i) 1074 ENDDO 1075 1076 p_mg(:,j,i) = tmp 1077 1046 1078 ENDDO 1047 !$OMP END PARALLEL 1048 1049 CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' ) 1050 1051 END SUBROUTINE sort_k_to_even_odd_blocks 1052 1053 1054 !------------------------------------------------------------------------------! 1079 ENDDO 1080 !$OMP END PARALLEL 1081 1082 CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' ) 1083 1084 END SUBROUTINE sort_k_to_even_odd_blocks 1085 1086 1087 !--------------------------------------------------------------------------------------------------! 1055 1088 ! Description: 1056 1089 ! ------------ 1057 !> Sort k- Dimension from sequential into blocks of even and odd.1058 !> This is required to vectorize the red-black subroutine.1059 ! > Version for 1D-REAL arrays1060 !------------------------------------------------------------------------------! 1061 SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel ) 1062 1063 1064 USE indices,&1065 ONLY: nzb,nzt_mg1066 1067 1068 1069 1070 1071 1072 1073 1074 ! 1075 !-- 1076 INTEGER(iwp) :: ind!< index variable along z1077 INTEGER(iwp) :: k !< index variable along z1078 1079 1080 1081 ! 1082 !-- 1083 1084 1085 1086 1087 1088 1089 ! 1090 !-- 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 !------------------------------------------------------------------------------ !1090 !> Sort k-dimension from sequential into blocks of even and odd. This is required to vectorize the 1091 !> red-black subroutine. Version for 1D-REAL arrays 1092 !--------------------------------------------------------------------------------------------------! 1093 SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel ) 1094 1095 1096 USE indices, & 1097 ONLY: nzb, & 1098 nzt_mg 1099 1100 IMPLICIT NONE 1101 1102 INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1103 1104 REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) :: f_mg !< 1D input array 1105 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: f_mg_b !< 1D output array 1106 1107 ! 1108 !-- Local variables 1109 INTEGER(iwp) :: ind !< index variable along z 1110 INTEGER(iwp) :: k !< index variable along z 1111 1112 1113 ind = nzb - 1 1114 ! 1115 !-- Sort the data with even k index 1116 DO k = nzb, nzt_mg(glevel), 2 1117 ind = ind + 1 1118 IF ( k >= nzb+1 .AND. k <= nzt_mg(glevel) ) THEN 1119 f_mg_b(ind) = f_mg(k) 1120 ENDIF 1121 ENDDO 1122 ! 1123 !-- Sort the data with odd k index 1124 DO k = nzb+1, nzt_mg(glevel)+1, 2 1125 ind = ind + 1 1126 IF( k >= nzb+1 .AND. k <= nzt_mg(glevel) ) THEN 1127 f_mg_b(ind) = f_mg(k) 1128 ENDIF 1129 ENDDO 1130 1131 END SUBROUTINE sort_k_to_even_odd_blocks_1d 1132 1133 1134 !--------------------------------------------------------------------------------------------------! 1102 1135 ! Description: 1103 1136 ! ------------ 1104 !> Sort k-Dimension from sequential into blocks of even and odd. 1105 !> This is required to vectorize the red-black subroutine. 1106 !> Version for 2D-INTEGER arrays 1107 !------------------------------------------------------------------------------! 1108 ! SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel ) 1109 ! 1110 ! 1111 ! USE indices, & 1112 ! ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1113 ! 1114 ! IMPLICIT NONE 1115 ! 1116 ! INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1117 ! 1118 ! INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1, & 1119 ! nys_mg(glevel)-1:nyn_mg(glevel)+1, & 1120 ! nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: & 1121 ! i_mg !< array to be sorted 1122 ! 1123 !-- Local variables 1124 ! INTEGER(iwp) :: i !< index variabel along x 1125 ! INTEGER(iwp) :: j !< index variable along y 1126 ! INTEGER(iwp) :: k !< index variable along z 1127 ! INTEGER(iwp) :: l !< grid level 1128 ! INTEGER(iwp) :: ind !< index variable along z 1129 ! INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< temporary odd-even sorted array 1130 ! 1131 ! 1132 ! CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' ) 1133 ! 1134 ! l = glevel 1135 ! ind_even_odd = even_odd_level(l) 1136 ! 1137 ! DO i = nxl_mg(l)-1, nxr_mg(l)+1 1138 ! DO j = nys_mg(l)-1, nyn_mg(l)+1 1139 ! 1140 ! 1141 !-- Sort the data with even k index 1142 ! ind = nzb-1 1143 ! DO k = nzb, nzt_mg(l), 2 1144 ! ind = ind + 1 1145 ! tmp(ind) = i_mg(k,j,i) 1146 ! ENDDO 1147 ! 1148 !-- Sort the data with odd k index 1149 ! DO k = nzb+1, nzt_mg(l)+1, 2 1150 ! ind = ind + 1 1151 ! tmp(ind) = i_mg(k,j,i) 1152 ! ENDDO 1153 ! 1154 ! i_mg(:,j,i) = tmp 1155 ! 1137 !> Sort k-dimension from sequential into blocks of even and odd. This is required to vectorize the 1138 !> red-black subroutine. Version for 2D-INTEGER arrays 1139 !--------------------------------------------------------------------------------------------------! 1140 ! SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel ) 1141 ! 1142 ! 1143 ! USE indices, & 1144 ! ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1145 ! 1146 ! IMPLICIT NONE 1147 ! 1148 ! INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1149 ! 1150 ! INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1, & 1151 ! nys_mg(glevel)-1:nyn_mg(glevel)+1, & 1152 ! nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: & 1153 ! i_mg !< array to be sorted 1154 !! 1155 !!-- Local variables 1156 ! INTEGER(iwp) :: i !< index variabel along x 1157 ! INTEGER(iwp) :: j !< index variable along y 1158 ! INTEGER(iwp) :: k !< index variable along z 1159 ! INTEGER(iwp) :: l !< grid level 1160 ! INTEGER(iwp) :: ind !< index variable along z 1161 ! INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< temporary odd-even sorted array 1162 ! 1163 ! 1164 ! CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' ) 1165 ! 1166 ! l = glevel 1167 ! ind_even_odd = even_odd_level(l) 1168 ! 1169 ! DO i = nxl_mg(l)-1, nxr_mg(l)+1 1170 ! DO j = nys_mg(l)-1, nyn_mg(l)+1 1171 ! 1172 ! 1173 !-- Sort the data with even k index 1174 ! ind = nzb-1 1175 ! DO k = nzb, nzt_mg(l), 2 1176 ! ind = ind + 1 1177 ! tmp(ind) = i_mg(k,j,i) 1156 1178 ! ENDDO 1179 ! 1180 !-- Sort the data with odd k index 1181 ! DO k = nzb+1, nzt_mg(l)+1, 2 1182 ! ind = ind + 1 1183 ! tmp(ind) = i_mg(k,j,i) 1184 ! ENDDO 1185 ! 1186 ! i_mg(:,j,i) = tmp 1187 ! 1157 1188 ! ENDDO 1158 ! 1159 ! CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' ) 1160 ! 1161 ! END SUBROUTINE sort_k_to_even_odd_blocks_int 1162 1163 1164 !------------------------------------------------------------------------------! 1189 ! ENDDO 1190 ! 1191 ! CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' ) 1192 ! 1193 ! END SUBROUTINE sort_k_to_even_odd_blocks_int 1194 1195 1196 !--------------------------------------------------------------------------------------------------! 1165 1197 ! Description: 1166 1198 ! ------------ 1167 1199 !> Sort k-dimension from blocks of even and odd into sequential 1168 !------------------------------------------------------------------------------! 1169 SUBROUTINE sort_k_to_sequential( p_mg ) 1170 1171 1172 USE control_parameters, & 1173 ONLY: grid_level 1174 1175 USE indices, & 1176 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1177 1178 IMPLICIT NONE 1179 1180 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1181 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1182 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 1183 p_mg !< array to be sorted 1184 ! 1185 !-- Local variables 1186 INTEGER(iwp) :: i !< index variable along x 1187 INTEGER(iwp) :: j !< index variable along y 1188 INTEGER(iwp) :: k !< index variable along z 1189 INTEGER(iwp) :: l !< grid level 1190 INTEGER(iwp) :: ind !< index variable along z 1191 1192 REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) :: tmp 1193 1194 1195 l = grid_level 1196 1197 !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp) 1198 !$OMP DO 1199 DO i = nxl_mg(l)-1, nxr_mg(l)+1 1200 DO j = nys_mg(l)-1, nyn_mg(l)+1 1201 1202 ind = nzb - 1 1203 tmp = p_mg(:,j,i) 1204 DO k = nzb, nzt_mg(l), 2 1205 ind = ind + 1 1206 p_mg(k,j,i) = tmp(ind) 1207 ENDDO 1208 1209 DO k = nzb+1, nzt_mg(l)+1, 2 1210 ind = ind + 1 1211 p_mg(k,j,i) = tmp(ind) 1212 ENDDO 1200 !--------------------------------------------------------------------------------------------------! 1201 SUBROUTINE sort_k_to_sequential( p_mg ) 1202 1203 1204 USE control_parameters, & 1205 ONLY: grid_level 1206 1207 USE indices, & 1208 ONLY: nxl_mg, & 1209 nxr_mg, & 1210 nys_mg, & 1211 nyn_mg, & 1212 nzb, & 1213 nzt_mg 1214 1215 IMPLICIT NONE 1216 1217 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1218 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< array to be sorted 1219 ! 1220 !-- Local variables 1221 INTEGER(iwp) :: i !< index variable along x 1222 INTEGER(iwp) :: j !< index variable along y 1223 INTEGER(iwp) :: k !< index variable along z 1224 INTEGER(iwp) :: l !< grid level 1225 INTEGER(iwp) :: ind !< index variable along z 1226 1227 REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) :: tmp !< 1228 1229 1230 l = grid_level 1231 1232 !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp) 1233 !$OMP DO 1234 DO i = nxl_mg(l)-1, nxr_mg(l)+1 1235 DO j = nys_mg(l)-1, nyn_mg(l)+1 1236 1237 ind = nzb - 1 1238 tmp = p_mg(:,j,i) 1239 DO k = nzb, nzt_mg(l), 2 1240 ind = ind + 1 1241 p_mg(k,j,i) = tmp(ind) 1242 ENDDO 1243 1244 DO k = nzb+1, nzt_mg(l)+1, 2 1245 ind = ind + 1 1246 p_mg(k,j,i) = tmp(ind) 1213 1247 ENDDO 1214 1248 ENDDO 1215 !$OMP END PARALLEL 1216 1217 END SUBROUTINE sort_k_to_sequential 1218 1219 1220 !------------------------------------------------------------------------------! 1249 ENDDO 1250 !$OMP END PARALLEL 1251 1252 END SUBROUTINE sort_k_to_sequential 1253 1254 1255 !--------------------------------------------------------------------------------------------------! 1221 1256 ! Description: 1222 1257 ! ------------ 1223 1258 !> Gather subdomain data from all PEs. 1224 !------------------------------------------------------------------------------ !1259 !--------------------------------------------------------------------------------------------------! 1225 1260 #if defined( __parallel ) 1226 SUBROUTINE mg_gather( f2, f2_sub ) 1227 1228 USE control_parameters, & 1229 ONLY: grid_level 1230 1231 USE cpulog, & 1232 ONLY: cpu_log, log_point_s 1233 1234 USE indices, & 1235 ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1236 1237 IMPLICIT NONE 1238 1239 INTEGER(iwp) :: i !< 1240 INTEGER(iwp) :: il !< 1241 INTEGER(iwp) :: ir !< 1242 INTEGER(iwp) :: j !< 1243 INTEGER(iwp) :: jn !< 1244 INTEGER(iwp) :: js !< 1245 INTEGER(iwp) :: k !< 1246 INTEGER(iwp) :: nwords !< 1247 1248 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1249 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1250 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !< 1251 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1252 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1253 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !< 1254 1255 REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, & 1256 mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1257 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !< 1258 1259 1260 CALL cpu_log( log_point_s(34), 'mg_gather', 'start' ) 1261 1262 f2_l = 0.0_wp 1263 1264 ! 1265 !-- Store the local subdomain array on the total array 1266 js = mg_loc_ind(3,myid) 1267 IF ( south_border_pe ) js = js - 1 1268 jn = mg_loc_ind(4,myid) 1269 IF ( north_border_pe ) jn = jn + 1 1270 il = mg_loc_ind(1,myid) 1271 IF ( left_border_pe ) il = il - 1 1272 ir = mg_loc_ind(2,myid) 1273 IF ( right_border_pe ) ir = ir + 1 1274 DO i = il, ir 1275 DO j = js, jn 1276 DO k = nzb, nzt_mg(grid_level)+1 1277 f2_l(k,j,i) = f2_sub(k,j,i) 1278 ENDDO 1261 SUBROUTINE mg_gather( f2, f2_sub ) 1262 1263 USE control_parameters, & 1264 ONLY: grid_level 1265 1266 USE cpulog, & 1267 ONLY: cpu_log, & 1268 log_point_s 1269 1270 USE indices, & 1271 ONLY: mg_loc_ind, & 1272 nxl_mg, & 1273 nxr_mg, & 1274 nys_mg, & 1275 nyn_mg, & 1276 nzb, & 1277 nzt_mg 1278 1279 IMPLICIT NONE 1280 1281 INTEGER(iwp) :: i !< 1282 INTEGER(iwp) :: il !< 1283 INTEGER(iwp) :: ir !< 1284 INTEGER(iwp) :: j !< 1285 INTEGER(iwp) :: jn !< 1286 INTEGER(iwp) :: js !< 1287 INTEGER(iwp) :: k !< 1288 INTEGER(iwp) :: nwords !< 1289 1290 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1291 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !< 1292 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1293 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !< 1294 1295 REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1296 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !< 1297 1298 1299 CALL cpu_log( log_point_s(34), 'mg_gather', 'start' ) 1300 1301 f2_l = 0.0_wp 1302 1303 ! 1304 !-- Store the local subdomain array on the total array 1305 js = mg_loc_ind(3,myid) 1306 IF ( south_border_pe ) js = js - 1 1307 jn = mg_loc_ind(4,myid) 1308 IF ( north_border_pe ) jn = jn + 1 1309 il = mg_loc_ind(1,myid) 1310 IF ( left_border_pe ) il = il - 1 1311 ir = mg_loc_ind(2,myid) 1312 IF ( right_border_pe ) ir = ir + 1 1313 DO i = il, ir 1314 DO j = js, jn 1315 DO k = nzb, nzt_mg(grid_level)+1 1316 f2_l(k,j,i) = f2_sub(k,j,i) 1279 1317 ENDDO 1280 1318 ENDDO 1281 1282 ! 1283 !-- Find out the number of array elements of the total array 1284 nwords = SIZE( f2 ) 1285 1286 ! 1287 !-- Gather subdomain data from all PEs 1288 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1289 CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), & 1290 f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), & 1291 nwords, MPI_REAL, MPI_SUM, comm2d, ierr ) 1292 1293 CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' ) 1294 1295 END SUBROUTINE mg_gather 1319 ENDDO 1320 1321 ! 1322 !-- Find out the number of array elements of the total array 1323 nwords = SIZE( f2 ) 1324 1325 ! 1326 !-- Gather subdomain data from all PEs 1327 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1328 CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), & 1329 f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), nwords, MPI_REAL, & 1330 MPI_SUM, comm2d, ierr ) 1331 1332 CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' ) 1333 1334 END SUBROUTINE mg_gather 1296 1335 #endif 1297 1336 1298 1337 1299 !------------------------------------------------------------------------------ !1338 !--------------------------------------------------------------------------------------------------! 1300 1339 ! Description: 1301 1340 ! ------------ 1302 !> @todo It might be possible to improve the speed of this routine by using 1303 !> non-blockingcommunication1304 !------------------------------------------------------------------------------ !1341 !> @todo It might be possible to improve the speed of this routine by using non-blocking 1342 !> communication 1343 !--------------------------------------------------------------------------------------------------! 1305 1344 #if defined( __parallel ) 1306 SUBROUTINE mg_scatter( p2, p2_sub ) 1307 1308 USE control_parameters, & 1309 ONLY: grid_level 1310 1311 USE cpulog, & 1312 ONLY: cpu_log, log_point_s 1313 1314 USE indices, & 1315 ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1316 1317 IMPLICIT NONE 1318 1319 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, & 1320 nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 1321 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !< 1322 1323 REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, & 1324 mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1325 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !< 1326 1327 1328 CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' ) 1329 1330 p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1331 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) 1332 1333 CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' ) 1334 1335 END SUBROUTINE mg_scatter 1345 SUBROUTINE mg_scatter( p2, p2_sub ) 1346 1347 USE control_parameters, & 1348 ONLY: grid_level 1349 1350 USE cpulog, & 1351 ONLY: cpu_log, & 1352 log_point_s 1353 1354 USE indices, & 1355 ONLY: mg_loc_ind, & 1356 nxl_mg, & 1357 nxr_mg, & 1358 nys_mg, & 1359 nyn_mg, & 1360 nzb, & 1361 nzt_mg 1362 1363 IMPLICIT NONE 1364 1365 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 1366 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !< 1367 1368 REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1369 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !< 1370 1371 1372 CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' ) 1373 1374 p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1375 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) 1376 1377 CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' ) 1378 1379 END SUBROUTINE mg_scatter 1336 1380 #endif 1337 1381 1338 !------------------------------------------------------------------------------ !1382 !--------------------------------------------------------------------------------------------------! 1339 1383 ! Description: 1340 1384 ! ------------ 1341 !> This is where the multigrid technique takes place. V- and W- Cycle are 1342 !> implemented and steered by the parameter "gamma". Parameter "nue" determines 1343 !> the convergence of the multigrid iterative solution. There are nue times 1344 !> RB-GS iterations. It should be set to "1" or "2", considering the time effort 1345 !> one would like to invest. Last choice shows a very good converging factor, 1346 !> but leads to an increase in computing time. 1347 !------------------------------------------------------------------------------! 1348 RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r ) 1349 1350 USE control_parameters, & 1351 ONLY: bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir, & 1352 child_domain, gamma_mg, grid_level_count, maximum_grid_level,& 1353 mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb 1354 1355 USE indices, & 1356 ONLY: mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn, & 1357 nyn_mg, nzb, nzt, nzt_mg 1358 1359 IMPLICIT NONE 1360 1361 INTEGER(iwp) :: i !< index variable along x 1362 INTEGER(iwp) :: j !< index variable along y 1363 INTEGER(iwp) :: k !< index variable along z 1364 INTEGER(iwp) :: nxl_mg_save !< 1365 INTEGER(iwp) :: nxr_mg_save !< 1366 INTEGER(iwp) :: nyn_mg_save !< 1367 INTEGER(iwp) :: nys_mg_save !< 1368 INTEGER(iwp) :: nzt_mg_save !< 1369 1370 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1371 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1372 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< 1373 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1374 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1375 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< 1376 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1377 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1378 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !< 1379 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1380 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1381 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !< 1382 1383 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, & 1384 nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 1385 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !< 1386 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, & 1387 nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 1388 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !< 1389 1390 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !< 1385 !> This is where the multigrid technique takes place. V- and W- Cycle are implemented and steered by 1386 !> the parameter "gamma". Parameter "nue" determines the convergence of the multigrid iterative 1387 !> solution. There are nue times RB-GS iterations. It should be set to "1" or "2", considering the 1388 !> time effort one would like to invest. Last choice shows a very good converging factor, but leads 1389 !> to an increase in computing time. 1390 !--------------------------------------------------------------------------------------------------! 1391 RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r ) 1392 1393 USE control_parameters, & 1394 ONLY: bc_lr_dirrad, & 1395 bc_lr_raddir, & 1396 bc_ns_dirrad, & 1397 bc_ns_raddir, & 1398 child_domain, & 1399 gamma_mg, & 1400 grid_level_count, & 1401 maximum_grid_level, & 1402 mg_switch_to_pe0_level, & 1403 mg_switch_to_pe0, & 1404 ngsrb 1405 1406 USE indices, & 1407 ONLY: mg_loc_ind, & 1408 nxl, & 1409 nxl_mg, & 1410 nxr, & 1411 nxr_mg, & 1412 nys, & 1413 nys_mg, & 1414 nyn, & 1415 nyn_mg, & 1416 nzb, & 1417 nzt, & 1418 nzt_mg 1419 1420 IMPLICIT NONE 1421 1422 INTEGER(iwp) :: i !< index variable along x 1423 INTEGER(iwp) :: j !< index variable along y 1424 INTEGER(iwp) :: k !< index variable along z 1425 INTEGER(iwp) :: nxl_mg_save !< 1426 INTEGER(iwp) :: nxr_mg_save !< 1427 INTEGER(iwp) :: nyn_mg_save !< 1428 INTEGER(iwp) :: nys_mg_save !< 1429 INTEGER(iwp) :: nzt_mg_save !< 1430 1431 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1432 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< 1433 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1434 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< 1435 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1436 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !< 1437 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1438 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !< 1439 1440 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 1441 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !< 1442 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 1443 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !< 1444 1445 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !< 1391 1446 1392 1447 #if defined( __parallel ) 1393 1448 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p2_sub !< 1394 1449 #endif 1395 1450 1396 1451 ! 1397 !-- Restriction to the coarsest grid 1398 10 IF ( grid_level == 1 ) THEN 1399 1400 ! 1401 !-- Solution on the coarsest grid. Double the number of Gauss-Seidel 1402 !-- iterations in order to get a more accurate solution. 1403 ngsrb = 2 * ngsrb 1404 1405 ind_even_odd = even_odd_level(grid_level) 1406 1407 CALL redblack( f_mg, p_mg ) 1408 1409 ngsrb = ngsrb / 2 1410 1411 1412 ELSEIF ( grid_level /= 1 ) THEN 1413 1414 grid_level_count(grid_level) = grid_level_count(grid_level) + 1 1415 1416 ! 1417 !-- Solution on the actual grid level 1418 ind_even_odd = even_odd_level(grid_level) 1419 1420 CALL redblack( f_mg, p_mg ) 1421 1422 ! 1423 !-- Determination of the actual residual 1424 CALL resid( f_mg, p_mg, r ) 1425 1426 !-- Restriction of the residual (finer grid values!) to the next coarser 1427 !-- grid. Therefore, the grid level has to be decremented now. nxl..nzt have 1428 !-- to be set to the coarse grid values, because these variables are needed 1429 !-- for the exchange of ghost points in routine exchange_horiz 1430 grid_level = grid_level - 1 1431 1452 !-- Restriction to the coarsest grid 1453 10 IF ( grid_level == 1 ) THEN 1454 1455 ! 1456 !-- Solution on the coarsest grid. Double the number of Gauss-Seidel iterations in order to get a 1457 !-- more accurate solution. 1458 ngsrb = 2 * ngsrb 1459 1460 ind_even_odd = even_odd_level(grid_level) 1461 1462 CALL redblack( f_mg, p_mg ) 1463 1464 ngsrb = ngsrb / 2 1465 1466 1467 ELSEIF ( grid_level /= 1 ) THEN 1468 1469 grid_level_count(grid_level) = grid_level_count(grid_level) + 1 1470 1471 ! 1472 !-- Solution on the actual grid level 1473 ind_even_odd = even_odd_level(grid_level) 1474 1475 CALL redblack( f_mg, p_mg ) 1476 1477 ! 1478 !-- Determination of the actual residual 1479 CALL resid( f_mg, p_mg, r ) 1480 1481 !-- Restriction of the residual (finer grid values!) to the next coarser grid. Therefore, the 1482 !-- grid level has to be decremented now. nxl..nzt have to be set to the coarse grid values, 1483 !-- because these variables are needed for the exchange of ghost points in routine exchange_horiz 1484 grid_level = grid_level - 1 1485 1486 nxl = nxl_mg(grid_level) 1487 nys = nys_mg(grid_level) 1488 nxr = nxr_mg(grid_level) 1489 nyn = nyn_mg(grid_level) 1490 nzt = nzt_mg(grid_level) 1491 1492 IF ( grid_level == mg_switch_to_pe0_level ) THEN 1493 1494 ! 1495 !-- From this level on, calculations are done on PE0 only. First, carry out restriction on the 1496 !-- subdomain. Therefore, indices of the level have to be changed to subdomain values in 1497 !-- between (otherwise, the restrict routine would expect the gathered array). 1498 1499 nxl_mg_save = nxl_mg(grid_level) 1500 nxr_mg_save = nxr_mg(grid_level) 1501 nys_mg_save = nys_mg(grid_level) 1502 nyn_mg_save = nyn_mg(grid_level) 1503 nzt_mg_save = nzt_mg(grid_level) 1504 nxl_mg(grid_level) = mg_loc_ind(1,myid) 1505 nxr_mg(grid_level) = mg_loc_ind(2,myid) 1506 nys_mg(grid_level) = mg_loc_ind(3,myid) 1507 nyn_mg(grid_level) = mg_loc_ind(4,myid) 1508 nzt_mg(grid_level) = mg_loc_ind(5,myid) 1509 nxl = mg_loc_ind(1,myid) 1510 nxr = mg_loc_ind(2,myid) 1511 nys = mg_loc_ind(3,myid) 1512 nyn = mg_loc_ind(4,myid) 1513 nzt = mg_loc_ind(5,myid) 1514 1515 ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1516 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ) 1517 1518 CALL restrict( f2_sub, r ) 1519 1520 ! 1521 !-- Restore the correct indices of this level 1522 nxl_mg(grid_level) = nxl_mg_save 1523 nxr_mg(grid_level) = nxr_mg_save 1524 nys_mg(grid_level) = nys_mg_save 1525 nyn_mg(grid_level) = nyn_mg_save 1526 nzt_mg(grid_level) = nzt_mg_save 1432 1527 nxl = nxl_mg(grid_level) 1528 nxr = nxr_mg(grid_level) 1433 1529 nys = nys_mg(grid_level) 1434 nxr = nxr_mg(grid_level)1435 1530 nyn = nyn_mg(grid_level) 1436 1531 nzt = nzt_mg(grid_level) 1437 1438 IF ( grid_level == mg_switch_to_pe0_level ) THEN 1439 1440 ! 1441 !-- From this level on, calculations are done on PE0 only. 1442 !-- First, carry out restriction on the subdomain. 1443 !-- Therefore, indices of the level have to be changed to subdomain 1444 !-- values in between (otherwise, the restrict routine would expect 1445 !-- the gathered array) 1446 1447 nxl_mg_save = nxl_mg(grid_level) 1448 nxr_mg_save = nxr_mg(grid_level) 1449 nys_mg_save = nys_mg(grid_level) 1450 nyn_mg_save = nyn_mg(grid_level) 1451 nzt_mg_save = nzt_mg(grid_level) 1452 nxl_mg(grid_level) = mg_loc_ind(1,myid) 1453 nxr_mg(grid_level) = mg_loc_ind(2,myid) 1454 nys_mg(grid_level) = mg_loc_ind(3,myid) 1455 nyn_mg(grid_level) = mg_loc_ind(4,myid) 1456 nzt_mg(grid_level) = mg_loc_ind(5,myid) 1457 nxl = mg_loc_ind(1,myid) 1458 nxr = mg_loc_ind(2,myid) 1459 nys = mg_loc_ind(3,myid) 1460 nyn = mg_loc_ind(4,myid) 1461 nzt = mg_loc_ind(5,myid) 1462 1463 ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1, & 1464 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1465 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ) 1466 1467 CALL restrict( f2_sub, r ) 1468 1469 ! 1470 !-- Restore the correct indices of this level 1471 nxl_mg(grid_level) = nxl_mg_save 1472 nxr_mg(grid_level) = nxr_mg_save 1473 nys_mg(grid_level) = nys_mg_save 1474 nyn_mg(grid_level) = nyn_mg_save 1475 nzt_mg(grid_level) = nzt_mg_save 1476 nxl = nxl_mg(grid_level) 1477 nxr = nxr_mg(grid_level) 1478 nys = nys_mg(grid_level) 1479 nyn = nyn_mg(grid_level) 1480 nzt = nzt_mg(grid_level) 1481 ! 1482 !-- Gather all arrays from the subdomains on PE0 1532 ! 1533 !-- Gather all arrays from the subdomains on PE0 1483 1534 #if defined( __parallel ) 1484 1535 CALL mg_gather( f2, f2_sub ) 1485 1536 #endif 1486 1537 1487 1538 ! 1488 !-- Set switch for routine exchange_horiz, that no ghostpoint exchange 1489 !-- has to be carried out from now on 1490 mg_switch_to_pe0 = .TRUE. 1491 1492 ! 1493 !-- In case of non-cyclic lateral boundary conditions, both in- and 1494 !-- outflow conditions have to be used on all PEs after the switch, 1495 !-- because then they have the total domain. 1539 !-- Set switch for routine exchange_horiz, that no ghostpoint exchange has to be carried out 1540 !-- from now on 1541 mg_switch_to_pe0 = .TRUE. 1542 1543 ! 1544 !-- In case of non-cyclic lateral boundary conditions, both in- and outflow conditions have to 1545 !-- be used on all PEs after the switch, because then they have the total domain. 1546 IF ( bc_lr_dirrad ) THEN 1547 bc_dirichlet_l = .TRUE. 1548 bc_dirichlet_r = .FALSE. 1549 bc_radiation_l = .FALSE. 1550 bc_radiation_r = .TRUE. 1551 ELSEIF ( bc_lr_raddir ) THEN 1552 bc_dirichlet_l = .FALSE. 1553 bc_dirichlet_r = .TRUE. 1554 bc_radiation_l = .TRUE. 1555 bc_radiation_r = .FALSE. 1556 ELSEIF ( child_domain .OR. nesting_offline ) THEN 1557 bc_dirichlet_l = .TRUE. 1558 bc_dirichlet_r = .TRUE. 1559 ENDIF 1560 1561 IF ( bc_ns_dirrad ) THEN 1562 bc_dirichlet_n = .TRUE. 1563 bc_dirichlet_s = .FALSE. 1564 bc_radiation_n = .FALSE. 1565 bc_radiation_s = .TRUE. 1566 ELSEIF ( bc_ns_raddir ) THEN 1567 bc_dirichlet_n = .FALSE. 1568 bc_dirichlet_s = .TRUE. 1569 bc_radiation_n = .TRUE. 1570 bc_radiation_s = .FALSE. 1571 ELSEIF ( child_domain .OR. nesting_offline) THEN 1572 bc_dirichlet_s = .TRUE. 1573 bc_dirichlet_n = .TRUE. 1574 ENDIF 1575 1576 DEALLOCATE( f2_sub ) 1577 1578 ELSE 1579 1580 CALL restrict( f2, r ) 1581 1582 ind_even_odd = even_odd_level(grid_level) ! Must be after restrict 1583 1584 ENDIF 1585 1586 p2 = 0.0_wp 1587 1588 ! 1589 !-- Repeat the same procedure untill the coarsest grid is reached 1590 CALL next_mg_level( f2, p2, p3, r ) 1591 1592 ENDIF 1593 1594 ! 1595 !-- Now follows the prolongation 1596 IF ( grid_level >= 2 ) THEN 1597 1598 ! 1599 !-- Prolongation of the new residual. The values are transferred from the coarse to the next 1600 !-- finer grid. 1601 IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN 1602 1603 #if defined( __parallel ) 1604 ! 1605 !-- At this level, the new residual first has to be scattered from PE0 to the other PEs 1606 ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1607 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ) 1608 1609 CALL mg_scatter( p2, p2_sub ) 1610 1611 ! 1612 !-- Therefore, indices of the previous level have to be changed to subdomain values in between 1613 !-- (otherwise, the prolong routine would expect the gathered array). 1614 nxl_mg_save = nxl_mg(grid_level-1) 1615 nxr_mg_save = nxr_mg(grid_level-1) 1616 nys_mg_save = nys_mg(grid_level-1) 1617 nyn_mg_save = nyn_mg(grid_level-1) 1618 nzt_mg_save = nzt_mg(grid_level-1) 1619 nxl_mg(grid_level-1) = mg_loc_ind(1,myid) 1620 nxr_mg(grid_level-1) = mg_loc_ind(2,myid) 1621 nys_mg(grid_level-1) = mg_loc_ind(3,myid) 1622 nyn_mg(grid_level-1) = mg_loc_ind(4,myid) 1623 nzt_mg(grid_level-1) = mg_loc_ind(5,myid) 1624 1625 ! 1626 !-- Set switch for routine exchange_horiz, that ghost point exchange has to be carried out 1627 !-- again from now on 1628 mg_switch_to_pe0 = .FALSE. 1629 1630 ! 1631 !-- For non-cyclic lateral boundary conditions and in case of nesting, restore the in-/outflow 1632 !-- conditions. 1633 bc_dirichlet_l = .FALSE.; bc_dirichlet_r = .FALSE. 1634 bc_dirichlet_n = .FALSE.; bc_dirichlet_s = .FALSE. 1635 bc_radiation_l = .FALSE.; bc_radiation_r = .FALSE. 1636 bc_radiation_n = .FALSE.; bc_radiation_s = .FALSE. 1637 1638 IF ( pleft == MPI_PROC_NULL ) THEN 1639 IF ( bc_lr_dirrad .OR. child_domain .OR. nesting_offline ) THEN 1640 bc_dirichlet_l = .TRUE. 1641 ELSEIF ( bc_lr_raddir ) THEN 1642 bc_radiation_l = .TRUE. 1643 ENDIF 1644 ENDIF 1645 1646 IF ( pright == MPI_PROC_NULL ) THEN 1496 1647 IF ( bc_lr_dirrad ) THEN 1497 bc_dirichlet_l = .TRUE.1498 bc_dirichlet_r = .FALSE.1499 bc_radiation_l = .FALSE.1500 1648 bc_radiation_r = .TRUE. 1501 ELSEIF ( bc_lr_raddir ) THEN 1502 bc_dirichlet_l = .FALSE. 1503 bc_dirichlet_r = .TRUE. 1504 bc_radiation_l = .TRUE. 1505 bc_radiation_r = .FALSE. 1506 ELSEIF ( child_domain .OR. nesting_offline ) THEN 1507 bc_dirichlet_l = .TRUE. 1649 ELSEIF ( bc_lr_raddir .OR. child_domain .OR. nesting_offline ) THEN 1508 1650 bc_dirichlet_r = .TRUE. 1509 1651 ENDIF 1510 1652 ENDIF 1653 1654 IF ( psouth == MPI_PROC_NULL ) THEN 1511 1655 IF ( bc_ns_dirrad ) THEN 1512 bc_dirichlet_n = .TRUE.1513 bc_dirichlet_s = .FALSE.1514 bc_radiation_n = .FALSE.1515 1656 bc_radiation_s = .TRUE. 1657 ELSEIF ( bc_ns_raddir .OR. child_domain .OR. nesting_offline ) THEN 1658 bc_dirichlet_s = .TRUE. 1659 ENDIF 1660 ENDIF 1661 1662 IF ( pnorth == MPI_PROC_NULL ) THEN 1663 IF ( bc_ns_dirrad .OR. child_domain .OR. nesting_offline ) THEN 1664 bc_dirichlet_n = .TRUE. 1516 1665 ELSEIF ( bc_ns_raddir ) THEN 1517 bc_dirichlet_n = .FALSE.1518 bc_dirichlet_s = .TRUE.1519 1666 bc_radiation_n = .TRUE. 1520 bc_radiation_s = .FALSE. 1521 ELSEIF ( child_domain .OR. nesting_offline) THEN 1522 bc_dirichlet_s = .TRUE. 1523 bc_dirichlet_n = .TRUE. 1524 ENDIF 1525 1526 DEALLOCATE( f2_sub ) 1527 1528 ELSE 1529 1530 CALL restrict( f2, r ) 1531 1532 ind_even_odd = even_odd_level(grid_level) ! must be after restrict 1533 1667 ENDIF 1534 1668 ENDIF 1535 1669 1536 p2 = 0.0_wp 1537 1538 ! 1539 !-- Repeat the same procedure till the coarsest grid is reached 1540 CALL next_mg_level( f2, p2, p3, r ) 1670 CALL prolong( p2_sub, p3 ) 1671 1672 ! 1673 !-- Restore the correct indices of the previous level 1674 nxl_mg(grid_level-1) = nxl_mg_save 1675 nxr_mg(grid_level-1) = nxr_mg_save 1676 nys_mg(grid_level-1) = nys_mg_save 1677 nyn_mg(grid_level-1) = nyn_mg_save 1678 nzt_mg(grid_level-1) = nzt_mg_save 1679 1680 DEALLOCATE( p2_sub ) 1681 #endif 1682 1683 ELSE 1684 1685 CALL prolong( p2, p3 ) 1541 1686 1542 1687 ENDIF 1543 1688 1544 1689 ! 1545 !-- Now follows the prolongation 1546 IF ( grid_level >= 2 ) THEN 1547 1548 ! 1549 !-- Prolongation of the new residual. The values are transferred 1550 !-- from the coarse to the next finer grid. 1551 IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN 1690 !-- Computation of the new pressure correction. Therefore, values from prior grids are added up 1691 !-- automatically stage by stage. 1692 DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1 1693 DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1 1694 DO k = nzb, nzt_mg(grid_level)+1 1695 p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i) 1696 ENDDO 1697 ENDDO 1698 ENDDO 1699 1700 ! 1701 !-- Relaxation of the new solution 1702 CALL redblack( f_mg, p_mg ) 1703 1704 ENDIF 1705 1706 1707 ! 1708 !-- The following few lines serve the steering of the multigrid scheme 1709 IF ( grid_level == maximum_grid_level ) THEN 1710 1711 GOTO 20 1712 1713 ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. & 1714 grid_level_count(grid_level) /= gamma_mg ) THEN 1715 1716 GOTO 10 1717 1718 ENDIF 1719 1720 ! 1721 !-- Reset counter for the next call of poismg 1722 grid_level_count(grid_level) = 0 1723 1724 ! 1725 !-- Continue with the next finer level. nxl..nzt have to be set to the finer grid values, because 1726 !-- these variables are needed for the exchange of ghost points in routine exchange_horiz. 1727 grid_level = grid_level + 1 1728 ind_even_odd = even_odd_level(grid_level) 1729 1730 nxl = nxl_mg(grid_level) 1731 nxr = nxr_mg(grid_level) 1732 nys = nys_mg(grid_level) 1733 nyn = nyn_mg(grid_level) 1734 nzt = nzt_mg(grid_level) 1735 1736 20 CONTINUE 1737 1738 END SUBROUTINE next_mg_level 1739 1740 1741 !--------------------------------------------------------------------------------------------------! 1742 ! Description: 1743 ! ------------ 1744 !> Initial settings for sorting k-dimension from sequential order (alternate even/odd) into blocks 1745 !> of even and odd or vice versa 1746 !--------------------------------------------------------------------------------------------------! 1747 SUBROUTINE init_even_odd_blocks 1748 1749 1750 USE arrays_3d, & 1751 ONLY: f1_mg, & 1752 f2_mg, & 1753 f3_mg 1754 1755 USE control_parameters, & 1756 ONLY: grid_level, & 1757 maximum_grid_level 1758 1759 USE indices, & 1760 ONLY: nzb, & 1761 nzt, & 1762 nzt_mg 1763 1764 USE indices, & 1765 ONLY: nzb, & 1766 nzt_mg 1767 1768 IMPLICIT NONE 1769 ! 1770 !-- Local variables 1771 INTEGER(iwp) :: i !< 1772 INTEGER(iwp) :: l !< 1773 1774 LOGICAL, SAVE :: lfirst = .TRUE. !< 1775 1776 1777 IF ( .NOT. lfirst ) RETURN 1778 1779 ALLOCATE( even_odd_level(maximum_grid_level) ) 1780 1781 ALLOCATE( f1_mg_b(nzb:nzt+1,maximum_grid_level), f2_mg_b(nzb:nzt+1,maximum_grid_level), & 1782 f3_mg_b(nzb:nzt+1,maximum_grid_level) ) 1783 1784 ! 1785 !-- Set border index between the even and odd block 1786 DO i = maximum_grid_level, 1, -1 1787 even_odd_level(i) = nzt_mg(i) / 2 1788 ENDDO 1789 1790 ! 1791 !-- Sort grid coefficients used in red/black scheme and for calculating the residual to block 1792 !-- (even/odd) structure 1793 DO l = maximum_grid_level, 1 , -1 1794 CALL sort_k_to_even_odd_blocks( f1_mg(nzb+1:nzt_mg(grid_level),l), & 1795 f1_mg_b(nzb:nzt_mg(grid_level)+1,l), l ) 1796 CALL sort_k_to_even_odd_blocks( f2_mg(nzb+1:nzt_mg(grid_level),l), & 1797 f2_mg_b(nzb:nzt_mg(grid_level)+1,l), l ) 1798 CALL sort_k_to_even_odd_blocks( f3_mg(nzb+1:nzt_mg(grid_level),l), & 1799 f3_mg_b(nzb:nzt_mg(grid_level)+1,l), l ) 1800 ENDDO 1801 1802 lfirst = .FALSE. 1803 1804 END SUBROUTINE init_even_odd_blocks 1805 1806 1807 !--------------------------------------------------------------------------------------------------! 1808 ! Description: 1809 ! ------------ 1810 !> Special exchange_horiz subroutine for use in redblack. Transfers only "red" or "black" data 1811 !> points. 1812 !--------------------------------------------------------------------------------------------------! 1813 SUBROUTINE special_exchange_horiz( p_mg, color ) 1814 1815 1816 USE control_parameters, & 1817 ONLY: grid_level 1552 1818 1553 1819 #if defined( __parallel ) 1554 ! 1555 !-- At this level, the new residual first has to be scattered from 1556 !-- PE0 to the other PEs 1557 ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1, & 1558 mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, & 1559 mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ) 1560 1561 CALL mg_scatter( p2, p2_sub ) 1562 1563 ! 1564 !-- Therefore, indices of the previous level have to be changed to 1565 !-- subdomain values in between (otherwise, the prolong routine would 1566 !-- expect the gathered array) 1567 nxl_mg_save = nxl_mg(grid_level-1) 1568 nxr_mg_save = nxr_mg(grid_level-1) 1569 nys_mg_save = nys_mg(grid_level-1) 1570 nyn_mg_save = nyn_mg(grid_level-1) 1571 nzt_mg_save = nzt_mg(grid_level-1) 1572 nxl_mg(grid_level-1) = mg_loc_ind(1,myid) 1573 nxr_mg(grid_level-1) = mg_loc_ind(2,myid) 1574 nys_mg(grid_level-1) = mg_loc_ind(3,myid) 1575 nyn_mg(grid_level-1) = mg_loc_ind(4,myid) 1576 nzt_mg(grid_level-1) = mg_loc_ind(5,myid) 1577 1578 ! 1579 !-- Set switch for routine exchange_horiz, that ghostpoint exchange 1580 !-- has to be carried again out from now on 1581 mg_switch_to_pe0 = .FALSE. 1582 1583 ! 1584 !-- For non-cyclic lateral boundary conditions and in case of nesting, 1585 !-- restore the in-/outflow conditions. 1586 bc_dirichlet_l = .FALSE.; bc_dirichlet_r = .FALSE. 1587 bc_dirichlet_n = .FALSE.; bc_dirichlet_s = .FALSE. 1588 bc_radiation_l = .FALSE.; bc_radiation_r = .FALSE. 1589 bc_radiation_n = .FALSE.; bc_radiation_s = .FALSE. 1590 1591 IF ( pleft == MPI_PROC_NULL ) THEN 1592 IF ( bc_lr_dirrad .OR. child_domain .OR. nesting_offline ) & 1593 THEN 1594 bc_dirichlet_l = .TRUE. 1595 ELSEIF ( bc_lr_raddir ) THEN 1596 bc_radiation_l = .TRUE. 1597 ENDIF 1598 ENDIF 1599 1600 IF ( pright == MPI_PROC_NULL ) THEN 1601 IF ( bc_lr_dirrad ) THEN 1602 bc_radiation_r = .TRUE. 1603 ELSEIF ( bc_lr_raddir .OR. child_domain .OR. & 1604 nesting_offline ) THEN 1605 bc_dirichlet_r = .TRUE. 1606 ENDIF 1607 ENDIF 1608 1609 IF ( psouth == MPI_PROC_NULL ) THEN 1610 IF ( bc_ns_dirrad ) THEN 1611 bc_radiation_s = .TRUE. 1612 ELSEIF ( bc_ns_raddir .OR. child_domain .OR. & 1613 nesting_offline ) THEN 1614 bc_dirichlet_s = .TRUE. 1615 ENDIF 1616 ENDIF 1617 1618 IF ( pnorth == MPI_PROC_NULL ) THEN 1619 IF ( bc_ns_dirrad .OR. child_domain .OR. nesting_offline ) & 1620 THEN 1621 bc_dirichlet_n = .TRUE. 1622 ELSEIF ( bc_ns_raddir ) THEN 1623 bc_radiation_n = .TRUE. 1624 ENDIF 1625 ENDIF 1626 1627 CALL prolong( p2_sub, p3 ) 1628 1629 ! 1630 !-- Restore the correct indices of the previous level 1631 nxl_mg(grid_level-1) = nxl_mg_save 1632 nxr_mg(grid_level-1) = nxr_mg_save 1633 nys_mg(grid_level-1) = nys_mg_save 1634 nyn_mg(grid_level-1) = nyn_mg_save 1635 nzt_mg(grid_level-1) = nzt_mg_save 1636 1637 DEALLOCATE( p2_sub ) 1820 USE control_parameters, & 1821 ONLY: mg_switch_to_pe0_level, & 1822 synchronous_exchange 1638 1823 #endif 1639 1824 1640 ELSE 1641 1642 CALL prolong( p2, p3 ) 1643 1644 ENDIF 1645 1646 ! 1647 !-- Computation of the new pressure correction. Therefore, 1648 !-- values from prior grids are added up automatically stage by stage. 1649 DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1 1650 DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1 1651 DO k = nzb, nzt_mg(grid_level)+1 1652 p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i) 1653 ENDDO 1654 ENDDO 1655 ENDDO 1656 1657 ! 1658 !-- Relaxation of the new solution 1659 CALL redblack( f_mg, p_mg ) 1660 1661 ENDIF 1662 1663 1664 ! 1665 !-- The following few lines serve the steering of the multigrid scheme 1666 IF ( grid_level == maximum_grid_level ) THEN 1667 1668 GOTO 20 1669 1670 ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. & 1671 grid_level_count(grid_level) /= gamma_mg ) THEN 1672 1673 GOTO 10 1674 1675 ENDIF 1676 1677 ! 1678 !-- Reset counter for the next call of poismg 1679 grid_level_count(grid_level) = 0 1680 1681 ! 1682 !-- Continue with the next finer level. nxl..nzt have to be 1683 !-- set to the finer grid values, because these variables are needed for the 1684 !-- exchange of ghost points in routine exchange_horiz 1685 grid_level = grid_level + 1 1686 ind_even_odd = even_odd_level(grid_level) 1825 USE indices, & 1826 ONLY: nxl_mg, & 1827 nxr_mg, & 1828 nys_mg, & 1829 nyn_mg, & 1830 nzb, & 1831 nzt_mg 1832 1833 #if defined( __parallel ) 1834 USE indices, & 1835 ONLY: nxl, & 1836 nxr, & 1837 nys, & 1838 nyn, & 1839 nzt 1840 #endif 1841 1842 IMPLICIT NONE 1843 1844 INTEGER(iwp), INTENT(IN) :: color !< flag for grid point type (red or black) 1845 1846 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1847 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< treated array 1848 1849 1850 #if defined ( __parallel ) 1851 ! 1852 !-- Local variables 1853 INTEGER(iwp) :: i !< index variable along x 1854 INTEGER(iwp) :: i1 !< index variable along x on coarse level 1855 INTEGER(iwp) :: i2 !< index variable along x on coarse level 1856 1857 INTEGER(iwp) :: j !< index variable along y 1858 INTEGER(iwp) :: j1 !< index variable along y on coarse level 1859 INTEGER(iwp) :: j2 !< index variable along y on coarse level 1860 INTEGER(iwp) :: k !< index variable along z 1861 INTEGER(iwp) :: l !< short for grid level 1862 INTEGER(iwp) :: jys !< index for lower local PE boundary along y 1863 INTEGER(iwp) :: jyn !< index for upper local PE boundary along y 1864 INTEGER(iwp) :: ixl !< index for lower local PE boundary along x 1865 INTEGER(iwp) :: ixr !< index for upper local PE boundary along x 1866 1867 LOGICAL :: synchronous_exchange_save !< dummy to reset synchronous_exchange to prescribed value 1868 1869 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: temp !< temporary array on next coarser grid level 1870 1871 synchronous_exchange_save = synchronous_exchange 1872 synchronous_exchange = .FALSE. 1873 1874 l = grid_level 1875 1876 ind_even_odd = even_odd_level(grid_level) 1877 1878 ! 1879 !-- Restricted transfer only on finer levels with enough data. Restricted transfer is not possible 1880 !-- for levels smaller or equal to 'switch to PE0 levels', since array bounds do not fit. 1881 !-- Moreover, it is not possible for the coarsest grid level, since the dimensions of temp are not 1882 !-- defined. For such cases, normal exchange_horiz is called. 1883 IF ( l > 1 .AND. l > mg_switch_to_pe0_level + 1 .AND. & 1884 ( ngp_xz(grid_level) >= 900 .OR. ngp_yz(grid_level) >= 900 ) ) THEN 1885 1886 jys = nys_mg(grid_level-1) 1887 jyn = nyn_mg(grid_level-1) 1888 ixl = nxl_mg(grid_level-1) 1889 ixr = nxr_mg(grid_level-1) 1890 ALLOCATE( temp(nzb:nzt_mg(l-1)+1,jys-1:jyn+1,ixl-1:ixr+1) ) 1891 ! 1892 !-- Handling the even k Values 1893 !-- Collecting data for the north - south exchange 1894 !-- Since only every second value has to be transfered, data are stored on the next coarser grid 1895 !-- level, because the arrays on that level have just the required size 1896 i1 = nxl_mg(grid_level-1) 1897 i2 = nxl_mg(grid_level-1) 1898 1899 DO i = nxl_mg(l), nxr_mg(l), 2 1900 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1901 1902 IF ( j == nys_mg(l) ) THEN 1903 !DIR$ IVDEP 1904 DO k = ind_even_odd+1, nzt_mg(l) 1905 temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i) 1906 ENDDO 1907 i1 = i1 + 1 1908 1909 ENDIF 1910 1911 IF ( j == nyn_mg(l) ) THEN 1912 !DIR$ IVDEP 1913 DO k = ind_even_odd+1, nzt_mg(l) 1914 temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i) 1915 ENDDO 1916 i2 = i2 + 1 1917 1918 ENDIF 1919 1920 ENDDO 1921 ENDDO 1922 1923 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1924 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 1925 1926 IF ( j == nys_mg(l) ) THEN 1927 !DIR$ IVDEP 1928 DO k = ind_even_odd+1, nzt_mg(l) 1929 temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i) 1930 ENDDO 1931 i1 = i1 + 1 1932 1933 ENDIF 1934 1935 IF ( j == nyn_mg(l) ) THEN 1936 !DIR$ IVDEP 1937 DO k = ind_even_odd+1, nzt_mg(l) 1938 temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i) 1939 ENDDO 1940 i2 = i2 + 1 1941 1942 ENDIF 1943 1944 ENDDO 1945 ENDDO 1946 1947 grid_level = grid_level-1 1687 1948 1688 1949 nxl = nxl_mg(grid_level) 1950 nys = nys_mg(grid_level) 1689 1951 nxr = nxr_mg(grid_level) 1690 nys = nys_mg(grid_level)1691 1952 nyn = nyn_mg(grid_level) 1692 1953 nzt = nzt_mg(grid_level) 1693 1954 1694 20 CONTINUE 1695 1696 END SUBROUTINE next_mg_level 1697 1698 1699 !------------------------------------------------------------------------------! 1700 ! Description: 1701 ! ------------ 1702 !> Initial settings for sorting k-dimension from sequential order (alternate 1703 !> even/odd) into blocks of even and odd or vice versa 1704 !------------------------------------------------------------------------------! 1705 SUBROUTINE init_even_odd_blocks 1706 1707 1708 USE arrays_3d, & 1709 ONLY: f1_mg, f2_mg, f3_mg 1710 1711 USE control_parameters, & 1712 ONLY: grid_level, maximum_grid_level 1713 1714 USE indices, & 1715 ONLY: nzb, nzt, nzt_mg 1716 1717 USE indices, & 1718 ONLY: nzb, nzt_mg 1719 1720 IMPLICIT NONE 1721 ! 1722 !-- Local variables 1723 INTEGER(iwp) :: i !< 1724 INTEGER(iwp) :: l !< 1725 1726 LOGICAL, SAVE :: lfirst = .TRUE. 1727 1728 1729 IF ( .NOT. lfirst ) RETURN 1730 1731 ALLOCATE( even_odd_level(maximum_grid_level) ) 1732 1733 ALLOCATE( f1_mg_b(nzb:nzt+1,maximum_grid_level), & 1734 f2_mg_b(nzb:nzt+1,maximum_grid_level), & 1735 f3_mg_b(nzb:nzt+1,maximum_grid_level) ) 1736 1737 ! 1738 !-- Set border index between the even and odd block 1739 DO i = maximum_grid_level, 1, -1 1740 even_odd_level(i) = nzt_mg(i) / 2 1955 send_receive = 'ns' 1956 CALL exchange_horiz( temp, 1 ) 1957 1958 grid_level = grid_level+1 1959 1960 i1 = nxl_mg(grid_level-1) 1961 i2 = nxl_mg(grid_level-1) 1962 1963 DO i = nxl_mg(l), nxr_mg(l), 2 1964 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1965 1966 IF ( j == nys_mg(l) ) THEN 1967 !DIR$ IVDEP 1968 DO k = ind_even_odd+1, nzt_mg(l) 1969 p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1) 1970 ENDDO 1971 i1 = i1 + 1 1972 1973 ENDIF 1974 1975 IF ( j == nyn_mg(l) ) THEN 1976 !DIR$ IVDEP 1977 DO k = ind_even_odd+1, nzt_mg(l) 1978 p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2) 1979 ENDDO 1980 i2 = i2 + 1 1981 1982 ENDIF 1983 1984 ENDDO 1741 1985 ENDDO 1742 1986 1743 ! 1744 !-- Sort grid coefficients used in red/black scheme and for calculating the 1745 !-- residual to block (even/odd) structure 1746 DO l = maximum_grid_level, 1 , -1 1747 CALL sort_k_to_even_odd_blocks( f1_mg(nzb+1:nzt_mg(grid_level),l), & 1748 f1_mg_b(nzb:nzt_mg(grid_level)+1,l), & 1749 l ) 1750 CALL sort_k_to_even_odd_blocks( f2_mg(nzb+1:nzt_mg(grid_level),l), & 1751 f2_mg_b(nzb:nzt_mg(grid_level)+1,l), & 1752 l ) 1753 CALL sort_k_to_even_odd_blocks( f3_mg(nzb+1:nzt_mg(grid_level),l), & 1754 f3_mg_b(nzb:nzt_mg(grid_level)+1,l), & 1755 l ) 1987 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1988 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 1989 1990 IF ( j == nys_mg(l) ) THEN 1991 !DIR$ IVDEP 1992 DO k = ind_even_odd+1, nzt_mg(l) 1993 p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1) 1994 ENDDO 1995 i1 = i1 + 1 1996 1997 ENDIF 1998 1999 IF ( j == nyn_mg(l) ) THEN 2000 !DIR$ IVDEP 2001 DO k = ind_even_odd+1, nzt_mg(l) 2002 p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2) 2003 ENDDO 2004 i2 = i2 + 1 2005 2006 ENDIF 2007 2008 ENDDO 1756 2009 ENDDO 1757 2010 1758 lfirst = .FALSE. 1759 1760 END SUBROUTINE init_even_odd_blocks 1761 1762 1763 !------------------------------------------------------------------------------! 1764 ! Description: 1765 ! ------------ 1766 !> Special exchange_horiz subroutine for use in redblack. Transfers only 1767 !> "red" or "black" data points. 1768 !------------------------------------------------------------------------------! 1769 SUBROUTINE special_exchange_horiz( p_mg, color ) 1770 1771 1772 USE control_parameters, & 1773 ONLY: grid_level 1774 1775 #if defined( __parallel ) 1776 USE control_parameters, & 1777 ONLY: mg_switch_to_pe0_level, synchronous_exchange 2011 ! 2012 !-- Collecting data for the left - right exchange. 2013 !-- Since only every second value has to be transfered, data are stored on the next coarser grid 2014 !-- level, because the arrays on that level have just the required size. 2015 j1 = nys_mg(grid_level-1) 2016 j2 = nys_mg(grid_level-1) 2017 2018 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2019 DO i = nxl_mg(l), nxr_mg(l), 2 2020 2021 IF ( i == nxl_mg(l) ) THEN 2022 !DIR$ IVDEP 2023 DO k = ind_even_odd+1, nzt_mg(l) 2024 temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i) 2025 ENDDO 2026 j1 = j1 + 1 2027 2028 ENDIF 2029 2030 IF ( i == nxr_mg(l) ) THEN 2031 !DIR$ IVDEP 2032 DO k = ind_even_odd+1, nzt_mg(l) 2033 temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i) 2034 ENDDO 2035 j2 = j2 + 1 2036 2037 ENDIF 2038 2039 ENDDO 2040 ENDDO 2041 2042 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2043 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2044 2045 IF ( i == nxl_mg(l) ) THEN 2046 !DIR$ IVDEP 2047 DO k = ind_even_odd+1, nzt_mg(l) 2048 temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i) 2049 ENDDO 2050 j1 = j1 + 1 2051 2052 ENDIF 2053 2054 IF ( i == nxr_mg(l) ) THEN 2055 !DIR$ IVDEP 2056 DO k = ind_even_odd+1, nzt_mg(l) 2057 temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i) 2058 ENDDO 2059 j2 = j2 + 1 2060 2061 ENDIF 2062 2063 ENDDO 2064 ENDDO 2065 2066 grid_level = grid_level-1 2067 send_receive = 'lr' 2068 2069 CALL exchange_horiz( temp, 1 ) 2070 2071 grid_level = grid_level+1 2072 2073 j1 = nys_mg(grid_level-1) 2074 j2 = nys_mg(grid_level-1) 2075 2076 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2077 DO i = nxl_mg(l), nxr_mg(l), 2 2078 2079 IF ( i == nxl_mg(l) ) THEN 2080 !DIR$ IVDEP 2081 DO k = ind_even_odd+1, nzt_mg(l) 2082 p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1) 2083 ENDDO 2084 j1 = j1 + 1 2085 2086 ENDIF 2087 2088 IF ( i == nxr_mg(l) ) THEN 2089 !DIR$ IVDEP 2090 DO k = ind_even_odd+1, nzt_mg(l) 2091 p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1) 2092 ENDDO 2093 j2 = j2 + 1 2094 2095 ENDIF 2096 2097 ENDDO 2098 ENDDO 2099 2100 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2101 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2102 2103 IF ( i == nxl_mg(l) ) THEN 2104 !DIR$ IVDEP 2105 DO k = ind_even_odd+1, nzt_mg(l) 2106 p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1) 2107 ENDDO 2108 j1 = j1 + 1 2109 2110 ENDIF 2111 2112 IF ( i == nxr_mg(l) ) THEN 2113 !DIR$ IVDEP 2114 DO k = ind_even_odd+1, nzt_mg(l) 2115 p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1) 2116 ENDDO 2117 j2 = j2 + 1 2118 2119 ENDIF 2120 2121 ENDDO 2122 ENDDO 2123 2124 ! 2125 !-- Now handling the even k values 2126 !-- Collecting data for the north - south exchange. 2127 !-- Since only every second value has to be transfered, data are stored on the next coarser grid 2128 !-- level, because the arrays on that level have just the required size. 2129 i1 = nxl_mg(grid_level-1) 2130 i2 = nxl_mg(grid_level-1) 2131 2132 DO i = nxl_mg(l), nxr_mg(l), 2 2133 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2134 2135 IF ( j == nys_mg(l) ) THEN 2136 !DIR$ IVDEP 2137 DO k = nzb+1, ind_even_odd 2138 temp(k,jys,i1) = p_mg(k,j,i) 2139 ENDDO 2140 i1 = i1 + 1 2141 2142 ENDIF 2143 2144 IF ( j == nyn_mg(l) ) THEN 2145 !DIR$ IVDEP 2146 DO k = nzb+1, ind_even_odd 2147 temp(k,jyn,i2) = p_mg(k,j,i) 2148 ENDDO 2149 i2 = i2 + 1 2150 2151 ENDIF 2152 2153 ENDDO 2154 ENDDO 2155 2156 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2157 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2158 2159 IF ( j == nys_mg(l) ) THEN 2160 !DIR$ IVDEP 2161 DO k = nzb+1, ind_even_odd 2162 temp(k,jys,i1) = p_mg(k,j,i) 2163 ENDDO 2164 i1 = i1 + 1 2165 2166 ENDIF 2167 2168 IF ( j == nyn_mg(l) ) THEN 2169 !DIR$ IVDEP 2170 DO k = nzb+1, ind_even_odd 2171 temp(k,jyn,i2) = p_mg(k,j,i) 2172 ENDDO 2173 i2 = i2 + 1 2174 2175 ENDIF 2176 2177 ENDDO 2178 ENDDO 2179 2180 grid_level = grid_level-1 2181 2182 send_receive = 'ns' 2183 CALL exchange_horiz( temp, 1 ) 2184 2185 grid_level = grid_level+1 2186 2187 i1 = nxl_mg(grid_level-1) 2188 i2 = nxl_mg(grid_level-1) 2189 2190 DO i = nxl_mg(l), nxr_mg(l), 2 2191 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2192 2193 IF ( j == nys_mg(l) ) THEN 2194 !DIR$ IVDEP 2195 DO k = nzb+1, ind_even_odd 2196 p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1) 2197 ENDDO 2198 i1 = i1 + 1 2199 2200 ENDIF 2201 2202 IF ( j == nyn_mg(l) ) THEN 2203 !DIR$ IVDEP 2204 DO k = nzb+1, ind_even_odd 2205 p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2) 2206 ENDDO 2207 i2 = i2 + 1 2208 2209 ENDIF 2210 2211 ENDDO 2212 ENDDO 2213 2214 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2215 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2216 2217 IF ( j == nys_mg(l) ) THEN 2218 !DIR$ IVDEP 2219 DO k = nzb+1, ind_even_odd 2220 p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1) 2221 ENDDO 2222 i1 = i1 + 1 2223 2224 ENDIF 2225 2226 IF ( j == nyn_mg(l) ) THEN 2227 !DIR$ IVDEP 2228 DO k = nzb+1, ind_even_odd 2229 p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2) 2230 ENDDO 2231 i2 = i2 + 1 2232 2233 ENDIF 2234 2235 ENDDO 2236 ENDDO 2237 2238 j1 = nys_mg(grid_level-1) 2239 j2 = nys_mg(grid_level-1) 2240 2241 DO i = nxl_mg(l), nxr_mg(l), 2 2242 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2243 2244 IF ( i == nxl_mg(l) ) THEN 2245 !DIR$ IVDEP 2246 DO k = nzb+1, ind_even_odd 2247 temp(k,j1,ixl) = p_mg(k,j,i) 2248 ENDDO 2249 j1 = j1 + 1 2250 2251 ENDIF 2252 2253 IF ( i == nxr_mg(l) ) THEN 2254 !DIR$ IVDEP 2255 DO k = nzb+1, ind_even_odd 2256 temp(k,j2,ixr) = p_mg(k,j,i) 2257 ENDDO 2258 j2 = j2 + 1 2259 2260 ENDIF 2261 2262 ENDDO 2263 ENDDO 2264 2265 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2266 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2267 2268 IF ( i == nxl_mg(l) ) THEN 2269 !DIR$ IVDEP 2270 DO k = nzb+1, ind_even_odd 2271 temp(k,j1,ixl) = p_mg(k,j,i) 2272 ENDDO 2273 j1 = j1 + 1 2274 2275 ENDIF 2276 2277 IF ( i == nxr_mg(l) ) THEN 2278 !DIR$ IVDEP 2279 DO k = nzb+1, ind_even_odd 2280 temp(k,j2,ixr) = p_mg(k,j,i) 2281 ENDDO 2282 j2 = j2 + 1 2283 2284 ENDIF 2285 2286 ENDDO 2287 ENDDO 2288 2289 grid_level = grid_level-1 2290 2291 send_receive = 'lr' 2292 CALL exchange_horiz( temp, 1 ) 2293 2294 grid_level = grid_level+1 2295 2296 nxl = nxl_mg(grid_level) 2297 nys = nys_mg(grid_level) 2298 nxr = nxr_mg(grid_level) 2299 nyn = nyn_mg(grid_level) 2300 nzt = nzt_mg(grid_level) 2301 2302 j1 = nys_mg(grid_level-1) 2303 j2 = nys_mg(grid_level-1) 2304 2305 DO i = nxl_mg(l), nxr_mg(l), 2 2306 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2307 2308 IF ( i == nxl_mg(l) ) THEN 2309 !DIR$ IVDEP 2310 DO k = nzb+1, ind_even_odd 2311 p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1) 2312 ENDDO 2313 j1 = j1 + 1 2314 2315 ENDIF 2316 2317 IF ( i == nxr_mg(l) ) THEN 2318 !DIR$ IVDEP 2319 DO k = nzb+1, ind_even_odd 2320 p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1) 2321 ENDDO 2322 j2 = j2 + 1 2323 2324 ENDIF 2325 2326 ENDDO 2327 ENDDO 2328 2329 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2330 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2331 2332 IF ( i == nxl_mg(l) ) THEN 2333 !DIR$ IVDEP 2334 DO k = nzb+1, ind_even_odd 2335 p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1) 2336 ENDDO 2337 j1 = j1 + 1 2338 2339 ENDIF 2340 2341 IF ( i == nxr_mg(l) ) THEN 2342 !DIR$ IVDEP 2343 DO k = nzb+1, ind_even_odd 2344 p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1) 2345 ENDDO 2346 j2 = j2 + 1 2347 2348 ENDIF 2349 2350 ENDDO 2351 ENDDO 2352 2353 DEALLOCATE( temp ) 2354 2355 ELSE 2356 2357 ! 2358 !-- Standard horizontal ghost boundary exchange for small coarse grid levels, where the transfer 2359 !-- time is latency bound 2360 CALL exchange_horiz( p_mg, 1 ) 2361 2362 ENDIF 2363 2364 ! 2365 !-- Reset values to default PALM setup 2366 synchronous_exchange = synchronous_exchange_save 2367 send_receive = 'al' 2368 #else 2369 2370 ! 2371 !-- Next line is to avoid compiler error due to unused dummy argument 2372 IF ( color == 1234567 ) RETURN 2373 ! 2374 !-- Standard horizontal ghost boundary exchange for small coarse grid levels, where the transfer 2375 !-- time is latency bound. 2376 CALL exchange_horiz( p_mg, 1 ) 1778 2377 #endif 1779 2378 1780 USE indices, & 1781 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1782 1783 #if defined( __parallel ) 1784 USE indices, & 1785 ONLY: nxl, nxr, nys, nyn, nzt 1786 #endif 1787 1788 IMPLICIT NONE 1789 1790 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1791 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1792 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 1793 p_mg !< treated array 1794 1795 INTEGER(iwp), INTENT(IN) :: color !< flag for grid point type (red or black) 1796 1797 #if defined ( __parallel ) 1798 ! 1799 !-- Local variables 1800 INTEGER(iwp) :: i !< index variable along x 1801 INTEGER(iwp) :: i1 !< index variable along x on coarse level 1802 INTEGER(iwp) :: i2 !< index variable along x on coarse level 1803 1804 INTEGER(iwp) :: j !< index variable along y 1805 INTEGER(iwp) :: j1 !< index variable along y on coarse level 1806 INTEGER(iwp) :: j2 !< index variable along y on coarse level 1807 INTEGER(iwp) :: k !< index variable along z 1808 INTEGER(iwp) :: l !< short for grid level 1809 INTEGER(iwp) :: jys !< index for lower local PE boundary along y 1810 INTEGER(iwp) :: jyn !< index for upper local PE boundary along y 1811 INTEGER(iwp) :: ixl !< index for lower local PE boundary along x 1812 INTEGER(iwp) :: ixr !< index for upper local PE boundary along x 1813 1814 LOGICAL :: synchronous_exchange_save !< dummy to reset synchronous_exchange to prescribed value 1815 1816 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: temp !< temporary array on next coarser grid level 1817 1818 synchronous_exchange_save = synchronous_exchange 1819 synchronous_exchange = .FALSE. 1820 1821 l = grid_level 1822 1823 ind_even_odd = even_odd_level(grid_level) 1824 1825 ! 1826 !-- Restricted transfer only on finer levels with enough data. 1827 !-- Restricted transfer is not possible for levels smaller or equal to 1828 !-- 'switch to PE0 levels', since array bounds does not fit. Moreover, 1829 !-- it is not possible for the coarsest grid level, since the dimensions 1830 !-- of temp are not defined. For such cases, normal exchange_horiz is called. 1831 IF ( l > 1 .AND. l > mg_switch_to_pe0_level + 1 .AND. & 1832 ( ngp_xz(grid_level) >= 900 .OR. ngp_yz(grid_level) >= 900 ) ) THEN 1833 1834 jys = nys_mg(grid_level-1) 1835 jyn = nyn_mg(grid_level-1) 1836 ixl = nxl_mg(grid_level-1) 1837 ixr = nxr_mg(grid_level-1) 1838 ALLOCATE( temp(nzb:nzt_mg(l-1)+1,jys-1:jyn+1,ixl-1:ixr+1) ) 1839 ! 1840 !-- Handling the even k Values 1841 !-- Collecting data for the north - south exchange 1842 !-- Since only every second value has to be transfered, data are stored 1843 !-- on the next coarser grid level, because the arrays on that level 1844 !-- have just the required size 1845 i1 = nxl_mg(grid_level-1) 1846 i2 = nxl_mg(grid_level-1) 1847 1848 DO i = nxl_mg(l), nxr_mg(l), 2 1849 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1850 1851 IF ( j == nys_mg(l) ) THEN 1852 !DIR$ IVDEP 1853 DO k = ind_even_odd+1, nzt_mg(l) 1854 temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i) 1855 ENDDO 1856 i1 = i1 + 1 1857 1858 ENDIF 1859 1860 IF ( j == nyn_mg(l) ) THEN 1861 !DIR$ IVDEP 1862 DO k = ind_even_odd+1, nzt_mg(l) 1863 temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i) 1864 ENDDO 1865 i2 = i2 + 1 1866 1867 ENDIF 1868 1869 ENDDO 1870 ENDDO 1871 1872 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1873 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 1874 1875 IF ( j == nys_mg(l) ) THEN 1876 !DIR$ IVDEP 1877 DO k = ind_even_odd+1, nzt_mg(l) 1878 temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i) 1879 ENDDO 1880 i1 = i1 + 1 1881 1882 ENDIF 1883 1884 IF ( j == nyn_mg(l) ) THEN 1885 !DIR$ IVDEP 1886 DO k = ind_even_odd+1, nzt_mg(l) 1887 temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i) 1888 ENDDO 1889 i2 = i2 + 1 1890 1891 ENDIF 1892 1893 ENDDO 1894 ENDDO 1895 1896 grid_level = grid_level-1 1897 1898 nxl = nxl_mg(grid_level) 1899 nys = nys_mg(grid_level) 1900 nxr = nxr_mg(grid_level) 1901 nyn = nyn_mg(grid_level) 1902 nzt = nzt_mg(grid_level) 1903 1904 send_receive = 'ns' 1905 CALL exchange_horiz( temp, 1 ) 1906 1907 grid_level = grid_level+1 1908 1909 i1 = nxl_mg(grid_level-1) 1910 i2 = nxl_mg(grid_level-1) 1911 1912 DO i = nxl_mg(l), nxr_mg(l), 2 1913 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1914 1915 IF ( j == nys_mg(l) ) THEN 1916 !DIR$ IVDEP 1917 DO k = ind_even_odd+1, nzt_mg(l) 1918 p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1) 1919 ENDDO 1920 i1 = i1 + 1 1921 1922 ENDIF 1923 1924 IF ( j == nyn_mg(l) ) THEN 1925 !DIR$ IVDEP 1926 DO k = ind_even_odd+1, nzt_mg(l) 1927 p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2) 1928 ENDDO 1929 i2 = i2 + 1 1930 1931 ENDIF 1932 1933 ENDDO 1934 ENDDO 1935 1936 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1937 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 1938 1939 IF ( j == nys_mg(l) ) THEN 1940 !DIR$ IVDEP 1941 DO k = ind_even_odd+1, nzt_mg(l) 1942 p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1) 1943 ENDDO 1944 i1 = i1 + 1 1945 1946 ENDIF 1947 1948 IF ( j == nyn_mg(l) ) THEN 1949 !DIR$ IVDEP 1950 DO k = ind_even_odd+1, nzt_mg(l) 1951 p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2) 1952 ENDDO 1953 i2 = i2 + 1 1954 1955 ENDIF 1956 1957 ENDDO 1958 ENDDO 1959 1960 ! 1961 !-- Collecting data for the left - right exchange 1962 !-- Since only every second value has to be transfered, data are stored 1963 !-- on the next coarser grid level, because the arrays on that level 1964 !-- have just the required size 1965 j1 = nys_mg(grid_level-1) 1966 j2 = nys_mg(grid_level-1) 1967 1968 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1969 DO i = nxl_mg(l), nxr_mg(l), 2 1970 1971 IF ( i == nxl_mg(l) ) THEN 1972 !DIR$ IVDEP 1973 DO k = ind_even_odd+1, nzt_mg(l) 1974 temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i) 1975 ENDDO 1976 j1 = j1 + 1 1977 1978 ENDIF 1979 1980 IF ( i == nxr_mg(l) ) THEN 1981 !DIR$ IVDEP 1982 DO k = ind_even_odd+1, nzt_mg(l) 1983 temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i) 1984 ENDDO 1985 j2 = j2 + 1 1986 1987 ENDIF 1988 1989 ENDDO 1990 ENDDO 1991 1992 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 1993 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1994 1995 IF ( i == nxl_mg(l) ) THEN 1996 !DIR$ IVDEP 1997 DO k = ind_even_odd+1, nzt_mg(l) 1998 temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i) 1999 ENDDO 2000 j1 = j1 + 1 2001 2002 ENDIF 2003 2004 IF ( i == nxr_mg(l) ) THEN 2005 !DIR$ IVDEP 2006 DO k = ind_even_odd+1, nzt_mg(l) 2007 temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i) 2008 ENDDO 2009 j2 = j2 + 1 2010 2011 ENDIF 2012 2013 ENDDO 2014 ENDDO 2015 2016 grid_level = grid_level-1 2017 send_receive = 'lr' 2018 2019 CALL exchange_horiz( temp, 1 ) 2020 2021 grid_level = grid_level+1 2022 2023 j1 = nys_mg(grid_level-1) 2024 j2 = nys_mg(grid_level-1) 2025 2026 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2027 DO i = nxl_mg(l), nxr_mg(l), 2 2028 2029 IF ( i == nxl_mg(l) ) THEN 2030 !DIR$ IVDEP 2031 DO k = ind_even_odd+1, nzt_mg(l) 2032 p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1) 2033 ENDDO 2034 j1 = j1 + 1 2035 2036 ENDIF 2037 2038 IF ( i == nxr_mg(l) ) THEN 2039 !DIR$ IVDEP 2040 DO k = ind_even_odd+1, nzt_mg(l) 2041 p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1) 2042 ENDDO 2043 j2 = j2 + 1 2044 2045 ENDIF 2046 2047 ENDDO 2048 ENDDO 2049 2050 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2051 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2052 2053 IF ( i == nxl_mg(l) ) THEN 2054 !DIR$ IVDEP 2055 DO k = ind_even_odd+1, nzt_mg(l) 2056 p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1) 2057 ENDDO 2058 j1 = j1 + 1 2059 2060 ENDIF 2061 2062 IF ( i == nxr_mg(l) ) THEN 2063 !DIR$ IVDEP 2064 DO k = ind_even_odd+1, nzt_mg(l) 2065 p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1) 2066 ENDDO 2067 j2 = j2 + 1 2068 2069 ENDIF 2070 2071 ENDDO 2072 ENDDO 2073 2074 ! 2075 !-- Now handling the even k values 2076 !-- Collecting data for the north - south exchange 2077 !-- Since only every second value has to be transfered, data are stored 2078 !-- on the next coarser grid level, because the arrays on that level 2079 !-- have just the required size 2080 i1 = nxl_mg(grid_level-1) 2081 i2 = nxl_mg(grid_level-1) 2082 2083 DO i = nxl_mg(l), nxr_mg(l), 2 2084 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2085 2086 IF ( j == nys_mg(l) ) THEN 2087 !DIR$ IVDEP 2088 DO k = nzb+1, ind_even_odd 2089 temp(k,jys,i1) = p_mg(k,j,i) 2090 ENDDO 2091 i1 = i1 + 1 2092 2093 ENDIF 2094 2095 IF ( j == nyn_mg(l) ) THEN 2096 !DIR$ IVDEP 2097 DO k = nzb+1, ind_even_odd 2098 temp(k,jyn,i2) = p_mg(k,j,i) 2099 ENDDO 2100 i2 = i2 + 1 2101 2102 ENDIF 2103 2104 ENDDO 2105 ENDDO 2106 2107 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2108 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2109 2110 IF ( j == nys_mg(l) ) THEN 2111 !DIR$ IVDEP 2112 DO k = nzb+1, ind_even_odd 2113 temp(k,jys,i1) = p_mg(k,j,i) 2114 ENDDO 2115 i1 = i1 + 1 2116 2117 ENDIF 2118 2119 IF ( j == nyn_mg(l) ) THEN 2120 !DIR$ IVDEP 2121 DO k = nzb+1, ind_even_odd 2122 temp(k,jyn,i2) = p_mg(k,j,i) 2123 ENDDO 2124 i2 = i2 + 1 2125 2126 ENDIF 2127 2128 ENDDO 2129 ENDDO 2130 2131 grid_level = grid_level-1 2132 2133 send_receive = 'ns' 2134 CALL exchange_horiz( temp, 1 ) 2135 2136 grid_level = grid_level+1 2137 2138 i1 = nxl_mg(grid_level-1) 2139 i2 = nxl_mg(grid_level-1) 2140 2141 DO i = nxl_mg(l), nxr_mg(l), 2 2142 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2143 2144 IF ( j == nys_mg(l) ) THEN 2145 !DIR$ IVDEP 2146 DO k = nzb+1, ind_even_odd 2147 p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1) 2148 ENDDO 2149 i1 = i1 + 1 2150 2151 ENDIF 2152 2153 IF ( j == nyn_mg(l) ) THEN 2154 !DIR$ IVDEP 2155 DO k = nzb+1, ind_even_odd 2156 p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2) 2157 ENDDO 2158 i2 = i2 + 1 2159 2160 ENDIF 2161 2162 ENDDO 2163 ENDDO 2164 2165 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2166 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2167 2168 IF ( j == nys_mg(l) ) THEN 2169 !DIR$ IVDEP 2170 DO k = nzb+1, ind_even_odd 2171 p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1) 2172 ENDDO 2173 i1 = i1 + 1 2174 2175 ENDIF 2176 2177 IF ( j == nyn_mg(l) ) THEN 2178 !DIR$ IVDEP 2179 DO k = nzb+1, ind_even_odd 2180 p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2) 2181 ENDDO 2182 i2 = i2 + 1 2183 2184 ENDIF 2185 2186 ENDDO 2187 ENDDO 2188 2189 j1 = nys_mg(grid_level-1) 2190 j2 = nys_mg(grid_level-1) 2191 2192 DO i = nxl_mg(l), nxr_mg(l), 2 2193 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2194 2195 IF ( i == nxl_mg(l) ) THEN 2196 !DIR$ IVDEP 2197 DO k = nzb+1, ind_even_odd 2198 temp(k,j1,ixl) = p_mg(k,j,i) 2199 ENDDO 2200 j1 = j1 + 1 2201 2202 ENDIF 2203 2204 IF ( i == nxr_mg(l) ) THEN 2205 !DIR$ IVDEP 2206 DO k = nzb+1, ind_even_odd 2207 temp(k,j2,ixr) = p_mg(k,j,i) 2208 ENDDO 2209 j2 = j2 + 1 2210 2211 ENDIF 2212 2213 ENDDO 2214 ENDDO 2215 2216 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2217 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2218 2219 IF ( i == nxl_mg(l) ) THEN 2220 !DIR$ IVDEP 2221 DO k = nzb+1, ind_even_odd 2222 temp(k,j1,ixl) = p_mg(k,j,i) 2223 ENDDO 2224 j1 = j1 + 1 2225 2226 ENDIF 2227 2228 IF ( i == nxr_mg(l) ) THEN 2229 !DIR$ IVDEP 2230 DO k = nzb+1, ind_even_odd 2231 temp(k,j2,ixr) = p_mg(k,j,i) 2232 ENDDO 2233 j2 = j2 + 1 2234 2235 ENDIF 2236 2237 ENDDO 2238 ENDDO 2239 2240 grid_level = grid_level-1 2241 2242 send_receive = 'lr' 2243 CALL exchange_horiz( temp, 1 ) 2244 2245 grid_level = grid_level+1 2246 2247 nxl = nxl_mg(grid_level) 2248 nys = nys_mg(grid_level) 2249 nxr = nxr_mg(grid_level) 2250 nyn = nyn_mg(grid_level) 2251 nzt = nzt_mg(grid_level) 2252 2253 j1 = nys_mg(grid_level-1) 2254 j2 = nys_mg(grid_level-1) 2255 2256 DO i = nxl_mg(l), nxr_mg(l), 2 2257 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 2258 2259 IF ( i == nxl_mg(l) ) THEN 2260 !DIR$ IVDEP 2261 DO k = nzb+1, ind_even_odd 2262 p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1) 2263 ENDDO 2264 j1 = j1 + 1 2265 2266 ENDIF 2267 2268 IF ( i == nxr_mg(l) ) THEN 2269 !DIR$ IVDEP 2270 DO k = nzb+1, ind_even_odd 2271 p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1) 2272 ENDDO 2273 j2 = j2 + 1 2274 2275 ENDIF 2276 2277 ENDDO 2278 ENDDO 2279 2280 DO i = nxl_mg(l)+1, nxr_mg(l), 2 2281 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 2282 2283 IF ( i == nxl_mg(l) ) THEN 2284 !DIR$ IVDEP 2285 DO k = nzb+1, ind_even_odd 2286 p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1) 2287 ENDDO 2288 j1 = j1 + 1 2289 2290 ENDIF 2291 2292 IF ( i == nxr_mg(l) ) THEN 2293 !DIR$ IVDEP 2294 DO k = nzb+1, ind_even_odd 2295 p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1) 2296 ENDDO 2297 j2 = j2 + 1 2298 2299 ENDIF 2300 2301 ENDDO 2302 ENDDO 2303 2304 DEALLOCATE( temp ) 2305 2306 ELSE 2307 2308 ! 2309 !-- Standard horizontal ghost boundary exchange for small coarse grid 2310 !-- levels, where the transfer time is latency bound 2311 CALL exchange_horiz( p_mg, 1 ) 2312 2313 ENDIF 2314 2315 ! 2316 !-- Reset values to default PALM setup 2317 synchronous_exchange = synchronous_exchange_save 2318 send_receive = 'al' 2319 #else 2320 2321 ! 2322 !-- Next line is to avoid compile error due to unused dummy argument 2323 IF ( color == 1234567 ) RETURN 2324 ! 2325 !-- Standard horizontal ghost boundary exchange for small coarse grid 2326 !-- levels, where the transfer time is latency bound 2327 CALL exchange_horiz( p_mg, 1 ) 2328 #endif 2329 2330 END SUBROUTINE special_exchange_horiz 2379 END SUBROUTINE special_exchange_horiz 2331 2380 2332 2381 END MODULE poismg_mod
Note: See TracChangeset
for help on using the changeset viewer.
