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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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