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

Last change on this file since 4456 was 4432, checked in by raasch, 4 years ago

bugfix for previous revision (vector directive was changed by mistake)

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