source: palm/trunk/SOURCE/poismg_fast_mod.f90 @ 1928

Last change on this file since 1928 was 1905, checked in by suehring, 9 years ago

last commit documented

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