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

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

last commit documented

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