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

Last change on this file since 2564 was 2298, checked in by raasch, 7 years ago

write_binary is of type LOGICAL now, MPI2-related code removed, obsolete variables removed, sendrecv_in_background related parts removed, missing variable descriptions added

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