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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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