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

Last change on this file since 3768 was 3725, checked in by raasch, 2 years ago

modifications to avoid compiler warnings about unused variables, temperton-fft: GOTO statements replaced, file re-formatted corresponding to coding standards, ssh-calls for compilations on remote systems modified to avoid output of login messages on specific systems changed again (palmbuild, reverted as before r3549), error messages for failed restarts extended (palmrun)

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