source: palm/trunk/SOURCE/poismg_mod.f90 @ 2021

Last change on this file since 2021 was 2021, checked in by suehring, 5 years ago

Bugfix, restore flags for nest boundaries in multigrid solver; bugfix: setting Neumann boundary conditions for topography arrays in case of non-cyclic boundary conditions

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