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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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