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

Last change on this file since 3046 was 2939, checked in by suehring, 3 years ago

Set lateral boundary conditions for divergence in multigrid solver

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