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

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