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

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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