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

Last change on this file since 2051 was 2038, checked in by knoop, 8 years ago

last commit documented

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