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

Last change on this file since 1988 was 1935, checked in by hellstea, 8 years ago

last commit documented

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