source: palm/trunk/SOURCE/poismg_mod.f90

Last change on this file was 4832, checked in by raasch, 9 months ago

bugfix in redblack algorithm: lower i,j indices need to start alternatively with even or odd value on the coarsest grid level, if the subdomain has an uneven number of gridpoints along x/y; some debug output flushed

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