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

Last change on this file since 2001 was 2001, checked in by knoop, 8 years ago

last commit documented

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