source: palm/tags/release-6.0/SOURCE/poismg_mod.f90 @ 3984

Last change on this file since 3984 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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