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

Last change on this file since 1904 was 1904, checked in by suehring, 8 years ago

Bugfix: enable special_exchange_horiz only for finer grid levels

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