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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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