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

Last change on this file since 2084 was 2084, checked in by knoop, 5 years ago

Bugfix: Multigrid now usable with anelastic approximation

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