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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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