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

Last change on this file since 2074 was 2074, checked in by raasch, 5 years ago

last commit documented

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