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

Last change on this file since 2106 was 2101, checked in by suehring, 8 years ago

last commit documented

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