source: palm/trunk/SOURCE/poismg_noopt_mod.f90

Last change on this file was 4832, checked in by raasch, 5 months ago

bugfix in redblack algorithm: lower i,j indices need to start alternatively with even or odd value on the coarsest grid level, if the subdomain has an uneven number of gridpoints along x/y; some debug output flushed

  • Property svn:keywords set to Id
File size: 96.6 KB
Line 
1!> @file poismg_noopt_mod.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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: poismg_noopt_mod.f90 4832 2021-01-07 07:15:12Z banzhafs $
27! bugfix in redblack algorithm: lower i,j indices need to start alternatively with even or odd
28! value on the coarsest grid level, if the subdomain has an uneven number of gridpoints along x/y
29!
30! 4828 2021-01-05 11:21:41Z Giersch
31! File re-formatted to follow the PALM coding standard
32!
33! 4457 2020-03-11 14:20:43Z raasch
34! Use statement for exchange horiz added
35!
36! 4429 2020-02-27 15:24:30Z raasch
37! Bugfix: cpp-directives added for serial mode
38!
39! 4414 2020-02-19 20:16:04Z suehring
40! Remove double-declared use only construct.
41!
42! 4360 2020-01-07 11:25:50Z suehring
43! Introduction of wall_flags_total_0, which currently sets bits based on static
44! topography information used in wall_flags_static_0
45!
46! 4329 2019-12-10 15:46:36Z motisi
47! Renamed wall_flags_0 to wall_flags_static_0
48!
49! 4182 2019-08-22 15:20:23Z scharf
50! Corrected "Former revisions" section
51!
52! 3655 2019-01-07 16:51:22Z knoop
53! Unused variables removed
54!
55! Revision 1.1  2001/07/20 13:10:51  raasch
56! Initial revision
57!
58!--------------------------------------------------------------------------------------------------!
59! Description:
60! ------------
61!> Solves the Poisson equation for the perturbation pressure with a multigrid V- or W-Cycle scheme.
62!>
63!> This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
64!> September 2000 - July 2001.
65!>
66!> @attention Loop unrolling and cache optimization in SOR-Red/Black method still does not give the
67!>            expected speedup!
68!>
69!> @todo Further work required.
70!> @todo Formatting adjustments required (indention after modularization)
71!--------------------------------------------------------------------------------------------------!
72 MODULE poismg_noopt_mod
73
74    USE control_parameters,                                                                        &
75        ONLY:  bc_dirichlet_l,                                                                     &
76               bc_dirichlet_n,                                                                     &
77               bc_dirichlet_r,                                                                     &
78               bc_dirichlet_s,                                                                     &
79               bc_radiation_l,                                                                     &
80               bc_radiation_n,                                                                     &
81               bc_radiation_r,                                                                     &
82               bc_radiation_s,                                                                     &
83               child_domain,                                                                       &
84               grid_level,                                                                         &
85               nesting_offline
86
87    USE cpulog,                                                                                    &
88        ONLY:  cpu_log,                                                                            &
89               log_point_s
90
91    USE exchange_horiz_mod,                                                                        &
92        ONLY:  exchange_horiz,                                                                     &
93               exchange_horiz_int
94
95    USE kinds
96
97    USE pegrid
98
99    PRIVATE
100
101    INTERFACE poismg_noopt
102       MODULE PROCEDURE poismg_noopt
103    END INTERFACE poismg_noopt
104
105    INTERFACE poismg_noopt_init
106       MODULE PROCEDURE poismg_noopt_init
107    END INTERFACE poismg_noopt_init
108
109    PUBLIC poismg_noopt, poismg_noopt_init
110
111 CONTAINS
112
113    SUBROUTINE poismg_noopt( r )
114
115       USE arrays_3d,                                                                              &
116           ONLY:  d,                                                                               &
117                  p_loc
118
119       USE control_parameters,                                                                     &
120           ONLY:  bc_lr_cyc,                                                                       &
121                  bc_ns_cyc,                                                                       &
122                  gathered_size,                                                                   &
123                  grid_level,                                                                      &
124                  grid_level_count,                                                                &
125                  ibc_p_t,                                                                         &
126                  maximum_grid_level,                                                              &
127                  message_string,                                                                  &
128                  mgcycles,                                                                        &
129                  mg_cycles,                                                                       &
130                  mg_switch_to_pe0_level,                                                          &
131                  residual_limit,                                                                  &
132                  subdomain_size
133
134       USE cpulog,                                                                                 &
135           ONLY:  cpu_log,                                                                         &
136                  log_point_s
137
138       USE indices,                                                                                &
139           ONLY:  nxl,                                                                             &
140                  nxlg,                                                                            &
141                  nxl_mg,                                                                          &
142                  nxr,                                                                             &
143                  nxrg,                                                                            &
144                  nxr_mg,                                                                          &
145                  nys,                                                                             &
146                  nysg,                                                                            &
147                  nys_mg,                                                                          &
148                  nyn,                                                                             &
149                  nyng,                                                                            &
150                  nyn_mg,                                                                          &
151                  nzb,                                                                             &
152                  nzt,                                                                             &
153                  nzt_mg
154
155       USE kinds
156
157       USE pegrid
158
159       IMPLICIT NONE
160
161       REAL(wp) ::  maxerror          !<
162       REAL(wp) ::  maximum_mgcycles  !<
163       REAL(wp) ::  residual_norm     !<
164
165       REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r  !<
166
167       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p3  !<
168
169
170       CALL cpu_log( log_point_s(29), 'poismg_noopt', 'start' )
171!
172!--    Initialize arrays and variables used in this subroutine
173
174!--    If the number of grid points of the gathered grid, which is collected on PE0, is larger than
175!--    the number of grid points of an PE, than array p3 will be enlarged.
176       IF ( gathered_size > subdomain_size )  THEN
177          ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(                                &
178                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,                    &
179                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(mg_switch_to_pe0_level)+1) )
180       ELSE
181          ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
182       ENDIF
183
184       p3 = 0.0_wp
185
186!
187!--    Ghost boundaries have to be added to divergence array.
188!--    Exchange routine needs to know the grid level!
189       grid_level = maximum_grid_level
190       CALL exchange_horiz( d, 1)
191!
192!--    Set bottom and top boundary conditions
193       d(nzb,:,:) = d(nzb+1,:,:)
194       IF ( ibc_p_t == 1 )  d(nzt+1,:,: ) = d(nzt,:,:)
195!
196!--    Set lateral boundary conditions in non-cyclic case
197       IF ( .NOT. bc_lr_cyc )  THEN
198          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  d(:,:,nxl-1) = d(:,:,nxl)
199          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  d(:,:,nxr+1) = d(:,:,nxr)
200       ENDIF
201       IF ( .NOT. bc_ns_cyc )  THEN
202          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  d(:,nyn+1,:) = d(:,nyn,:)
203          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  d(:,nys-1,:) = d(:,nys,:)
204       ENDIF
205
206!
207!--    Initiation of the multigrid scheme. Does n cycles until the residual is smaller than the
208!--    given limit. The accuracy of the solution of the poisson equation will increase with the
209!--    number of cycles. If the number of cycles is preset by the user, this number will be carried
210!--    out regardless of the accuracy.
211       grid_level_count =  0
212       mgcycles         =  0
213       IF ( mg_cycles == -1 )  THEN
214          maximum_mgcycles = 0
215          residual_norm    = 1.0_wp
216       ELSE
217          maximum_mgcycles = mg_cycles
218          residual_norm    = 0.0_wp
219       ENDIF
220
221       DO WHILE ( residual_norm > residual_limit  .OR.  mgcycles < maximum_mgcycles )
222
223          CALL next_mg_level_noopt( d, p_loc, p3, r)
224
225!
226!--       Calculate the residual if the user has not preset the number of cycles to be performed
227          IF ( maximum_mgcycles == 0 )  THEN
228             CALL resid_noopt( d, p_loc, r )
229             maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
230
231#if defined( __parallel )
232             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
233                CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
234#else
235                residual_norm = maxerror
236#endif
237             residual_norm = SQRT( residual_norm )
238          ENDIF
239
240          mgcycles = mgcycles + 1
241
242!
243!--       If the user has not limited the number of cycles, stop the run in case of insufficient
244!--       convergence
245          IF ( mgcycles > 1000  .AND.  mg_cycles == -1 )  THEN
246             message_string = 'no sufficient convergence within 1000 cycles'
247             CALL message( 'poismg_noopt', 'PA0283', 1, 2, 0, 6, 0 )
248          ENDIF
249
250       ENDDO
251
252       DEALLOCATE( p3 )
253
254!
255!--    Unset the grid level. Variable is used to determine the MPI datatypes for ghost point
256!--    exchange
257       grid_level = 0
258
259       CALL cpu_log( log_point_s(29), 'poismg_noopt', 'stop' )
260
261    END SUBROUTINE poismg_noopt
262
263
264!--------------------------------------------------------------------------------------------------!
265! Description:
266! ------------
267!> Computes the residual of the perturbation pressure.
268!--------------------------------------------------------------------------------------------------!
269 SUBROUTINE resid_noopt( f_mg, p_mg, r )
270
271
272    USE arrays_3d,                                                                                 &
273        ONLY:  f1_mg,                                                                              &
274               f2_mg,                                                                              &
275               f3_mg,                                                                              &
276               rho_air_mg
277
278    USE control_parameters,                                                                        &
279        ONLY:  bc_lr_cyc,                                                                          &
280               bc_ns_cyc,                                                                          &
281               ibc_p_b,                                                                            &
282               ibc_p_t
283
284    USE grid_variables,                                                                            &
285        ONLY:  ddx2_mg,                                                                            &
286               ddy2_mg
287
288    USE indices,                                                                                   &
289        ONLY:  flags,                                                                              &
290               nxl_mg,                                                                             &
291               nxr_mg,                                                                             &
292               nys_mg,                                                                             &
293               nyn_mg,                                                                             &
294               nzb,                                                                                &
295               nzt_mg,                                                                             &
296               wall_flags_1,                                                                       &
297               wall_flags_2,                                                                       &
298               wall_flags_3,                                                                       &
299               wall_flags_4,                                                                       &
300               wall_flags_5,                                                                       &
301               wall_flags_6,                                                                       &
302               wall_flags_7,                                                                       &
303               wall_flags_8,                                                                       &
304               wall_flags_9,                                                                       &
305               wall_flags_10
306
307
308
309
310    USE kinds
311
312    IMPLICIT NONE
313
314    INTEGER(iwp) ::  i  !<
315    INTEGER(iwp) ::  j  !<
316    INTEGER(iwp) ::  k  !<
317    INTEGER(iwp) ::  l  !<
318
319    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
320                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
321    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
322                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
323    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
324                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !<
325
326!
327!-- Calculate the residual
328    l = grid_level
329
330!
331!-- Choose flag array of this level
332    SELECT CASE ( l )
333       CASE ( 1 )
334          flags => wall_flags_1
335       CASE ( 2 )
336          flags => wall_flags_2
337       CASE ( 3 )
338          flags => wall_flags_3
339       CASE ( 4 )
340          flags => wall_flags_4
341       CASE ( 5 )
342          flags => wall_flags_5
343       CASE ( 6 )
344          flags => wall_flags_6
345       CASE ( 7 )
346          flags => wall_flags_7
347       CASE ( 8 )
348          flags => wall_flags_8
349       CASE ( 9 )
350          flags => wall_flags_9
351       CASE ( 10 )
352          flags => wall_flags_10
353    END SELECT
354
355!$OMP PARALLEL PRIVATE (i,j,k)
356!$OMP DO
357    DO  i = nxl_mg(l), nxr_mg(l)
358       DO  j = nys_mg(l), nyn_mg(l)
359          DO  k = nzb+1, nzt_mg(l)
360             r(k,j,i) = f_mg(k,j,i) - rho_air_mg(k,l) * ddx2_mg(l) *                               &
361                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *                            &
362                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +                            &
363                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *                            &
364                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )                            &
365                        - rho_air_mg(k,l) * ddy2_mg(l) *                                           &
366                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *                            &
367                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +                            &
368                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *                            &
369                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )                            &
370                        - f2_mg(k,l) *                                                             &
371                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *                            &
372                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )                            &
373                        - f3_mg(k,l) *                                                             &
374                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *                            &
375                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )                            &
376                        + f1_mg(k,l) * p_mg(k,j,i)
377!
378!--          Residual within topography should be zero
379             r(k,j,i) = r(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
380          ENDDO
381       ENDDO
382    ENDDO
383!$OMP END PARALLEL
384
385!
386!-- Horizontal boundary conditions
387    CALL exchange_horiz( r, 1)
388
389    IF ( .NOT. bc_lr_cyc )  THEN
390       IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
391          r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
392       ENDIF
393       IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
394          r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
395       ENDIF
396    ENDIF
397
398    IF ( .NOT. bc_ns_cyc )  THEN
399       IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
400          r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
401       ENDIF
402       IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
403          r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
404       ENDIF
405    ENDIF
406
407!
408!-- Boundary conditions at bottom and top of the domain.
409!-- These points are not handled by the above loop. Points may be within buildings, but that
410!-- doesn't matter.
411    IF ( ibc_p_b == 1 )  THEN
412       r(nzb,:,: ) = r(nzb+1,:,:)
413    ELSE
414       r(nzb,:,: ) = 0.0_wp
415    ENDIF
416
417    IF ( ibc_p_t == 1 )  THEN
418       r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
419    ELSE
420       r(nzt_mg(l)+1,:,: ) = 0.0_wp
421    ENDIF
422
423
424 END SUBROUTINE resid_noopt
425
426
427!--------------------------------------------------------------------------------------------------!
428! Description:
429! ------------
430!> Interpolates the residual on the next coarser grid with "full weighting" scheme.
431!--------------------------------------------------------------------------------------------------!
432 SUBROUTINE restrict_noopt( f_mg, r )
433
434
435    USE control_parameters,                                                                        &
436        ONLY:  bc_lr_cyc,                                                                          &
437               bc_ns_cyc,                                                                          &
438               grid_level,                                                                         &
439               ibc_p_b,                                                                            &
440               ibc_p_t
441
442    USE indices,                                                                                   &
443        ONLY:  flags,                                                                              &
444               nxl_mg,                                                                             &
445               nxr_mg,                                                                             &
446               nys_mg,                                                                             &
447               nyn_mg,                                                                             &
448               nzb,                                                                                &
449               nzt_mg,                                                                             &
450               wall_flags_1,                                                                       &
451               wall_flags_2,                                                                       &
452               wall_flags_3,                                                                       &
453               wall_flags_4,                                                                       &
454               wall_flags_5,                                                                       &
455               wall_flags_6,                                                                       &
456               wall_flags_7,                                                                       &
457               wall_flags_8,                                                                       &
458               wall_flags_9,                                                                       &
459               wall_flags_10
460
461
462    USE kinds
463
464    IMPLICIT NONE
465
466    INTEGER(iwp) ::  i   !<
467    INTEGER(iwp) ::  ic  !<
468    INTEGER(iwp) ::  j   !<
469    INTEGER(iwp) ::  jc  !<
470    INTEGER(iwp) ::  k   !<
471    INTEGER(iwp) ::  kc  !<
472    INTEGER(iwp) ::  l   !<
473
474    REAL(wp) ::  rkjim    !<
475    REAL(wp) ::  rkjip    !<
476    REAL(wp) ::  rkjmi    !<
477    REAL(wp) ::  rkjmim   !<
478    REAL(wp) ::  rkjmip   !<
479    REAL(wp) ::  rkjpi    !<
480    REAL(wp) ::  rkjpim   !<
481    REAL(wp) ::  rkjpip   !<
482    REAL(wp) ::  rkmji    !<
483    REAL(wp) ::  rkmjim   !<
484    REAL(wp) ::  rkmjip   !<
485    REAL(wp) ::  rkmjmi   !<
486    REAL(wp) ::  rkmjmim  !<
487    REAL(wp) ::  rkmjmip  !<
488    REAL(wp) ::  rkmjpi   !<
489    REAL(wp) ::  rkmjpim  !<
490    REAL(wp) ::  rkmjpip  !<
491
492    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
493                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1)    ::  f_mg  !<
494
495    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,  &
496                        nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r    !<
497
498!
499!-- Interpolate the residual
500    l = grid_level
501
502!
503!-- Choose flag array of the upper level
504    SELECT CASE ( l+1 )
505       CASE ( 1 )
506          flags => wall_flags_1
507       CASE ( 2 )
508          flags => wall_flags_2
509       CASE ( 3 )
510          flags => wall_flags_3
511       CASE ( 4 )
512          flags => wall_flags_4
513       CASE ( 5 )
514          flags => wall_flags_5
515       CASE ( 6 )
516          flags => wall_flags_6
517       CASE ( 7 )
518          flags => wall_flags_7
519       CASE ( 8 )
520          flags => wall_flags_8
521       CASE ( 9 )
522          flags => wall_flags_9
523       CASE ( 10 )
524          flags => wall_flags_10
525    END SELECT
526
527!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc, rkjim,rkjip,rkjpi,rkjmi,rkjmim,rkjpim, &
528!$OMP rkjmip, rkjpip,rkmji,rkmjim,rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim,rkmjmip,&
529!$OMP rkmjpip          )
530!$OMP DO
531    DO  ic = nxl_mg(l), nxr_mg(l)
532       i = 2*ic
533       DO  jc = nys_mg(l), nyn_mg(l)
534          j = 2*jc
535          DO  kc = nzb+1, nzt_mg(l)
536             k = 2*kc-1
537!
538!--          Use implicit Neumann BCs if the respective gridpoint is inside the building
539             rkjim   = r(k,j,i-1) + IBITS( flags(k,j,i-1), 6, 1 ) * ( r(k,j,i) - r(k,j,i-1) )
540             rkjip   = r(k,j,i+1) + IBITS( flags(k,j,i+1), 6, 1 ) * ( r(k,j,i) - r(k,j,i+1) )
541             rkjpi   = r(k,j+1,i) + IBITS( flags(k,j+1,i), 6, 1 ) * ( r(k,j,i) - r(k,j+1,i) )
542             rkjmi   = r(k,j-1,i) + IBITS( flags(k,j-1,i), 6, 1 ) * ( r(k,j,i) - r(k,j-1,i) )
543             rkjmim  = r(k,j-1,i-1) + IBITS( flags(k,j-1,i-1), 6, 1 ) * ( r(k,j,i) - r(k,j-1,i-1) )
544             rkjpim  = r(k,j+1,i-1) + IBITS( flags(k,j+1,i-1), 6, 1 ) * ( r(k,j,i) - r(k,j+1,i-1) )
545             rkjmip  = r(k,j-1,i+1) + IBITS( flags(k,j-1,i+1), 6, 1 ) * ( r(k,j,i) - r(k,j-1,i+1) )
546             rkjpip  = r(k,j+1,i+1) + IBITS( flags(k,j+1,i+1), 6, 1 ) * ( r(k,j,i) - r(k,j+1,i+1) )
547             rkmji   = r(k-1,j,i) + IBITS( flags(k-1,j,i), 6, 1 ) * ( r(k,j,i) - r(k-1,j,i) )
548             rkmjim  = r(k-1,j,i-1) + IBITS( flags(k-1,j,i-1), 6, 1 ) * ( r(k,j,i) - r(k-1,j,i-1) )
549             rkmjip  = r(k-1,j,i+1) + IBITS( flags(k-1,j,i+1), 6, 1 ) * ( r(k,j,i) - r(k-1,j,i+1) )
550             rkmjpi  = r(k-1,j+1,i) + IBITS( flags(k-1,j+1,i), 6, 1 ) * ( r(k,j,i) - r(k-1,j+1,i) )
551             rkmjmi  = r(k-1,j-1,i) + IBITS( flags(k-1,j-1,i), 6, 1 ) * ( r(k,j,i) - r(k-1,j-1,i) )
552             rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 )                          &
553                       * ( r(k,j,i) - r(k-1,j-1,i-1) )
554             rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 )                          &
555                       * ( r(k,j,i) - r(k-1,j+1,i-1) )
556             rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 )                          &
557                       * ( r(k,j,i) - r(k-1,j-1,i+1) )
558             rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 )                          &
559                       * ( r(k,j,i) - r(k-1,j+1,i+1) )
560
561             f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp                                                     &
562                              * ( 8.0_wp * r(k,j,i)                                                &
563                                  + 4.0_wp * ( rkjim  + rkjip  + rkjpi   + rkjmi  )                &
564                                  + 2.0_wp * ( rkjmim + rkjpim + rkjmip  + rkjpip )                &
565                                  + 4.0_wp * rkmji                                                 &
566                                  + 2.0_wp * ( rkmjim  + rkmjim  + rkmjpi  + rkmjmi  )             &
567                                  +          ( rkmjmim + rkmjpim + rkmjmip + rkmjpip )             &
568                                  + 4.0_wp * r(k+1,j,i)                                            &
569                                  + 2.0_wp * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   +                   &
570                                               r(k+1,j+1,i)   + r(k+1,j-1,i)   )                   &
571                                  +          ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) +                   &
572                                               r(k+1,j-1,i+1) + r(k+1,j+1,i+1) )                   &
573                                )
574
575!          f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                                                    &
576!                           8.0_wp * r(k,j,i)                                                       &
577!                         + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     +                            &
578!                                      r(k,j+1,i)     + r(k,j-1,i)     )                            &
579!                         + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   +                            &
580!                                      r(k,j-1,i+1)   + r(k,j+1,i+1)   )                            &
581!                         + 4.0_wp * r(k-1,j,i)                                                     &
582!                         + 2.0_wp * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   +                            &
583!                                      r(k-1,j+1,i)   + r(k-1,j-1,i)   )                            &
584!                         +          ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) +                            &
585!                                      r(k-1,j-1,i+1) + r(k-1,j+1,i+1) )                            &
586!                         + 4.0_wp * r(k+1,j,i)                                                     &
587!                         + 2.0_wp * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   +                            &
588!                                      r(k+1,j+1,i)   + r(k+1,j-1,i)   )                            &
589!                         +          ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) +                            &
590!                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) )                            &
591!                                             )
592          ENDDO
593       ENDDO
594    ENDDO
595!$OMP END PARALLEL
596
597!
598!-- Horizontal boundary conditions
599    CALL exchange_horiz( f_mg, 1)
600
601    IF ( .NOT. bc_lr_cyc )  THEN
602       IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
603          f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
604       ENDIF
605       IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
606          f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
607       ENDIF
608    ENDIF
609
610    IF ( .NOT. bc_ns_cyc )  THEN
611       IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
612          f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
613       ENDIF
614       IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
615          f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
616       ENDIF
617    ENDIF
618
619!
620!-- Boundary conditions at bottom and top of the domain.
621!-- These points are not handled by the above loop. Points may be within buildings, but that
622!-- doesn't matter.
623    IF ( ibc_p_b == 1 )  THEN
624       f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
625    ELSE
626       f_mg(nzb,:,: ) = 0.0_wp
627    ENDIF
628
629    IF ( ibc_p_t == 1 )  THEN
630       f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
631    ELSE
632       f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
633    ENDIF
634
635
636 END SUBROUTINE restrict_noopt
637
638
639!--------------------------------------------------------------------------------------------------!
640! Description:
641! ------------
642!> Interpolates the correction of the perturbation pressure to the next finer grid.
643!--------------------------------------------------------------------------------------------------!
644 SUBROUTINE prolong_noopt( p, temp )
645
646
647    USE control_parameters,                                                                        &
648        ONLY:  bc_lr_cyc,                                                                          &
649               bc_ns_cyc,                                                                          &
650               ibc_p_b,                                                                            &
651               ibc_p_t
652    USE indices,                                                                                   &
653        ONLY:  nxl_mg,                                                                             &
654               nxr_mg,                                                                             &
655               nys_mg,                                                                             &
656               nyn_mg,                                                                             &
657               nzb,                                                                                &
658               nzt_mg
659
660    USE kinds
661
662    IMPLICIT NONE
663
664    INTEGER(iwp) ::  i  !<
665    INTEGER(iwp) ::  j  !<
666    INTEGER(iwp) ::  k  !<
667    INTEGER(iwp) ::  l  !<
668
669    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,  &
670                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p  !<
671
672    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
673                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1)      ::  temp  !<
674
675
676!
677!-- First, store elements of the coarser grid on the next finer grid
678    l = grid_level
679
680!$OMP PARALLEL PRIVATE (i,j,k)
681!$OMP DO
682    DO  i = nxl_mg(l-1), nxr_mg(l-1)
683       DO  j = nys_mg(l-1), nyn_mg(l-1)
684!CDIR NODEP
685          DO  k = nzb+1, nzt_mg(l-1)
686!
687!--          Points of the coarse grid are directly stored on the next finer grid
688             temp(2*k-1,2*j,2*i) = p(k,j,i)
689!
690!--          Points between two coarse-grid points
691             temp(2*k-1,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
692             temp(2*k-1,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
693             temp(2*k,2*j,2*i)     = 0.5_wp * ( p(k,j,i) + p(k+1,j,i) )
694!
695!--          Points in the center of the planes stretched by four points of the coarse grid cube
696             temp(2*k-1,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +                       &
697                                                   p(k,j+1,i) + p(k,j+1,i+1) )
698             temp(2*k,2*j,2*i+1)     = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +                       &
699                                                   p(k+1,j,i) + p(k+1,j,i+1) )
700             temp(2*k,2*j+1,2*i)     = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) +                       &
701                                                   p(k+1,j,i) + p(k+1,j+1,i) )
702!
703!--          Points in the middle of coarse grid cube
704             temp(2*k,2*j+1,2*i+1) = 0.125_wp * ( p(k,j,i)     + p(k,j,i+1)   + p(k,j+1,i)   +     &
705                                                  p(k,j+1,i+1) + p(k+1,j,i)   + p(k+1,j,i+1) +     &
706                                                  p(k+1,j+1,i) + p(k+1,j+1,i+1) )
707          ENDDO
708       ENDDO
709    ENDDO
710!$OMP END PARALLEL
711
712!
713!-- Horizontal boundary conditions
714    CALL exchange_horiz( temp, 1)
715
716    IF ( .NOT. bc_lr_cyc )  THEN
717       IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
718          temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
719       ENDIF
720       IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
721          temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
722       ENDIF
723    ENDIF
724
725    IF ( .NOT. bc_ns_cyc )  THEN
726       IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
727          temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
728       ENDIF
729       IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
730          temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
731       ENDIF
732    ENDIF
733
734!
735!-- Bottom and top boundary conditions
736    IF ( ibc_p_b == 1 )  THEN
737       temp(nzb,:,: ) = temp(nzb+1,:,:)
738    ELSE
739       temp(nzb,:,: ) = 0.0_wp
740    ENDIF
741
742    IF ( ibc_p_t == 1 )  THEN
743       temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
744    ELSE
745       temp(nzt_mg(l)+1,:,: ) = 0.0_wp
746    ENDIF
747
748
749 END SUBROUTINE prolong_noopt
750
751
752!--------------------------------------------------------------------------------------------------!
753! Description:
754! ------------
755!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with 3D-Red-Black
756!> decomposition (GS-RB) is used.
757!--------------------------------------------------------------------------------------------------!
758 SUBROUTINE redblack_noopt( f_mg, p_mg )
759
760
761    USE arrays_3d,                                                                                 &
762        ONLY:  f1_mg,                                                                              &
763               f2_mg,                                                                              &
764               f3_mg,                                                                              &
765               rho_air_mg
766
767    USE control_parameters,                                                                        &
768        ONLY:  bc_lr_cyc,                                                                          &
769               bc_ns_cyc,                                                                          &
770               ibc_p_b,                                                                            &
771               ibc_p_t,                                                                            &
772               ngsrb
773
774    USE cpulog,                                                                                    &
775        ONLY:  cpu_log,                                                                            &
776               log_point_s
777
778    USE grid_variables,                                                                            &
779        ONLY:  ddx2_mg,                                                                            &
780               ddy2_mg
781
782    USE indices,                                                                                   &
783        ONLY:  flags,                                                                              &
784               nxl_mg,                                                                             &
785               nxr_mg,                                                                             &
786               nys_mg,                                                                             &
787               nyn_mg,                                                                             &
788               nzb,                                                                                &
789               nzt_mg,                                                                             &
790               wall_flags_1,                                                                       &
791               wall_flags_2,                                                                       &
792               wall_flags_3,                                                                       &
793               wall_flags_4,                                                                       &
794               wall_flags_5,                                                                       &
795               wall_flags_6,                                                                       &
796               wall_flags_7,                                                                       &
797               wall_flags_8,                                                                       &
798               wall_flags_9,                                                                       &
799               wall_flags_10
800
801
802    USE kinds
803
804    IMPLICIT NONE
805
806    INTEGER(iwp) ::  color  !<
807    INTEGER(iwp) ::  i      !<
808    INTEGER(iwp) ::  ic     !<
809    INTEGER(iwp) ::  j      !<
810    INTEGER(iwp) ::  jc     !<
811    INTEGER(iwp) ::  jj     !<
812    INTEGER(iwp) ::  k      !<
813    INTEGER(iwp) ::  l      !<
814    INTEGER(iwp) ::  n      !<
815    INTEGER(iwp) ::  save_nxl_mg  !< to save nxl_mg on coarsest level 1
816    INTEGER(iwp) ::  save_nys_mg  !< to save nys_mg on coarsest level 1
817
818    LOGICAL ::  adjust_lower_i_index  !< adjust lower limit of i loop in case of odd number of grid points
819    LOGICAL ::  adjust_lower_j_index  !< adjust lower limit of j loop in case of odd number of grid points
820    LOGICAL ::  unroll  !<
821
822    REAL(wp) ::  wall_left   !<
823    REAL(wp) ::  wall_north  !<
824    REAL(wp) ::  wall_right  !<
825    REAL(wp) ::  wall_south  !<
826    REAL(wp) ::  wall_total  !<
827    REAL(wp) ::  wall_top    !<
828
829    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
830                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
831    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
832                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
833
834    l = grid_level
835
836!
837!-- Choose flag array of this level
838    SELECT CASE ( l )
839       CASE ( 1 )
840          flags => wall_flags_1
841       CASE ( 2 )
842          flags => wall_flags_2
843       CASE ( 3 )
844          flags => wall_flags_3
845       CASE ( 4 )
846          flags => wall_flags_4
847       CASE ( 5 )
848          flags => wall_flags_5
849       CASE ( 6 )
850          flags => wall_flags_6
851       CASE ( 7 )
852          flags => wall_flags_7
853       CASE ( 8 )
854          flags => wall_flags_8
855       CASE ( 9 )
856          flags => wall_flags_9
857       CASE ( 10 )
858          flags => wall_flags_10
859    END SELECT
860
861    unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND.                                         &
862               MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
863
864!
865!-- The red/black decomposition requires that on the lower i,j indices need to start alternatively with an
866!-- even or odd value on the coarsest grid level, depending on the core-id, and if the subdomain has an
867!-- uneven number of gridpoints along x/y. Set the respective steering switches here.
868    IF ( l == 1  .AND.  MOD( myidx, 2 ) /= 0  .AND.  MOD( nxl_mg(l) - nxr_mg(l), 2 ) == 0 )  THEN
869       adjust_lower_i_index = .TRUE.
870       save_nxl_mg = nxl_mg(1)
871    ELSE
872       adjust_lower_i_index = .FALSE.
873    ENDIF
874    IF ( l == 1  .AND.  MOD( myidy, 2 ) /= 0  .AND.  MOD( nyn_mg(l) - nys_mg(l), 2 ) == 0 )  THEN
875       adjust_lower_j_index = .TRUE.
876       save_nys_mg = nys_mg(l)
877    ELSE
878       adjust_lower_j_index = .FALSE.
879    ENDIF
880
881
882    DO  n = 1, ngsrb
883
884       DO  color = 1, 2
885
886          IF ( .NOT. unroll )  THEN
887
888             CALL cpu_log( log_point_s(36), 'redblack_no_unroll_noopt', 'start' )
889
890             IF ( adjust_lower_i_index )  THEN
891                nxl_mg(l) = save_nxl_mg + 1
892             ENDIF
893
894             IF ( adjust_lower_j_index )  THEN
895                IF ( color == 1 )  THEN
896                  nys_mg(l) = save_nys_mg - 1
897                ELSE
898                  nys_mg(l) = save_nys_mg + 1
899                ENDIF
900             ENDIF
901!
902!--          Without unrolling of loops, no cache optimization
903             DO  i = nxl_mg(l), nxr_mg(l), 2
904                DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
905                   DO  k = nzb+1, nzt_mg(l), 2
906!                   p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
907!                             ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
908!                           + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
909!                           + f2_mg(k,l) * p_mg(k+1,j,i)                     &
910!                           + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
911!                                                    )
912
913                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
914                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
915                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
916                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
917                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
918                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
919                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
920                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
921                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
922                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
923                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
924                                    + f2_mg(k,l) *                                                 &
925                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
926                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
927                                    + f3_mg(k,l) *                                                 &
928                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
929                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
930                                    - f_mg(k,j,i)                                                  &
931                                      )
932                   ENDDO
933                ENDDO
934             ENDDO
935
936             IF ( adjust_lower_i_index )  THEN
937                nxl_mg(l) = save_nxl_mg - 1
938             ENDIF
939
940             IF ( adjust_lower_j_index )  THEN
941                IF ( color == 1 )  THEN
942                  nys_mg(l) = save_nys_mg + 1
943                ELSE
944                  nys_mg(l) = save_nys_mg - 1
945                ENDIF
946             ENDIF
947
948             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
949                DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
950                   DO  k = nzb+1, nzt_mg(l), 2
951                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
952                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
953                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
954                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
955                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
956                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
957                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
958                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
959                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
960                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
961                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
962                                    + f2_mg(k,l) *                                                 &
963                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
964                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
965                                    + f3_mg(k,l) *                                                 &
966                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
967                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
968                                    - f_mg(k,j,i)                                                  &
969                                      )
970                   ENDDO
971                ENDDO
972             ENDDO
973
974             IF ( adjust_lower_i_index )  THEN
975                nxl_mg(l) = save_nxl_mg + 1
976             ENDIF
977
978             IF ( adjust_lower_j_index )  THEN
979                IF ( color == 1 )  THEN
980                  nys_mg(l) = save_nys_mg + 1
981                ELSE
982                  nys_mg(l) = save_nys_mg - 1
983                ENDIF
984             ENDIF
985
986             DO  i = nxl_mg(l), nxr_mg(l), 2
987                DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
988                   DO  k = nzb+2, nzt_mg(l), 2
989                     p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                             &
990                                   * ( rho_air_mg(k,l) * ddx2_mg(l) *                              &
991                                       ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *             &
992                                                     ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +             &
993                                         p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *             &
994                                                     ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )             &
995                                   + rho_air_mg(k,l) * ddy2_mg(l) *                                &
996                                       ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *             &
997                                                     ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +             &
998                                         p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *             &
999                                                     ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )             &
1000                                   + f2_mg(k,l) *                                                  &
1001                                       ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *             &
1002                                                     ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )             &
1003                                   + f3_mg(k,l) *                                                  &
1004                                       ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *             &
1005                                                     ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )             &
1006                                   - f_mg(k,j,i)                                                   &
1007                                     )
1008                   ENDDO
1009                ENDDO
1010             ENDDO
1011
1012             IF ( adjust_lower_i_index )  THEN
1013                nxl_mg(l) = save_nxl_mg - 1
1014             ENDIF
1015
1016             IF ( adjust_lower_j_index )  THEN
1017                IF ( color == 1 )  THEN
1018                  nys_mg(l) = save_nys_mg - 1
1019                ELSE
1020                  nys_mg(l) = save_nys_mg + 1
1021                ENDIF
1022             ENDIF
1023
1024             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
1025                DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1026                   DO  k = nzb+2, nzt_mg(l), 2
1027                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1028                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1029                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1030                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1031                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1032                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1033                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1034                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1035                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1036                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1037                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1038                                    + f2_mg(k,l) *                                                 &
1039                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1040                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1041                                    + f3_mg(k,l) *                                                 &
1042                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1043                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1044                                    - f_mg(k,j,i)                                                  &
1045                                      )
1046                   ENDDO
1047                ENDDO
1048             ENDDO
1049             CALL cpu_log( log_point_s(36), 'redblack_no_unroll_noopt', 'stop' )
1050
1051          ELSE
1052
1053!
1054!--          Loop unrolling along y, only one i loop for better cache use
1055             CALL cpu_log( log_point_s(38), 'redblack_unroll_noopt', 'start' )
1056             DO  ic = nxl_mg(l), nxr_mg(l), 2
1057                DO  jc = nys_mg(l), nyn_mg(l), 4
1058                   i  = ic
1059                   jj = jc+2-color
1060                   DO  k = nzb+1, nzt_mg(l), 2
1061                      j = jj
1062                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1063                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1064                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1065                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1066                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1067                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1068                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1069                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1070                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1071                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1072                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1073                                    + f2_mg(k,l) *                                                 &
1074                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1075                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1076                                    + f3_mg(k,l) *                                                 &
1077                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1078                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1079                                    - f_mg(k,j,i)                                                  &
1080                                      )
1081
1082                      j = jj+2
1083                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1084                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1085                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1086                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1087                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1088                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1089                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1090                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1091                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1092                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1093                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1094                                    + f2_mg(k,l) *                                                 &
1095                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1096                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1097                                    + f3_mg(k,l) *                                                 &
1098                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1099                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1100                                    - f_mg(k,j,i)                                                  &
1101                                      )
1102                   ENDDO
1103
1104                   i  = ic+1
1105                   jj = jc+color-1
1106                   DO  k = nzb+1, nzt_mg(l), 2
1107                      j =jj
1108                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1109                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1110                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1111                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1112                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1113                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1114                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1115                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1116                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1117                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1118                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1119                                    + f2_mg(k,l) *                                                 &
1120                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1121                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1122                                    + f3_mg(k,l) *                                                 &
1123                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1124                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1125                                    - f_mg(k,j,i)                                                  &
1126                                      )
1127
1128                      j = jj+2
1129                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1130                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1131                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1132                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1133                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1134                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1135                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1136                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1137                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1138                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1139                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1140                                    + f2_mg(k,l) *                                                 &
1141                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1142                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1143                                    + f3_mg(k,l) *                                                 &
1144                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1145                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1146                                    - f_mg(k,j,i)                                                  &
1147                                      )
1148                   ENDDO
1149
1150                   i  = ic
1151                   jj = jc+color-1
1152                   DO  k = nzb+2, nzt_mg(l), 2
1153                      j =jj
1154                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1155                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1156                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1157                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1158                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1159                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1160                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1161                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1162                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1163                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1164                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1165                                    + f2_mg(k,l) *                                                 &
1166                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1167                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1168                                    + f3_mg(k,l) *                                                 &
1169                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1170                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1171                                    - f_mg(k,j,i)                                                  &
1172                                      )
1173
1174                      j = jj+2
1175                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1176                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1177                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1178                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1179                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1180                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1181                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1182                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1183                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1184                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1185                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1186                                    + f2_mg(k,l) *                                                 &
1187                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1188                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1189                                    + f3_mg(k,l) *                                                 &
1190                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1191                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1192                                    - f_mg(k,j,i)                                                  &
1193                                      )
1194                   ENDDO
1195
1196                   i  = ic+1
1197                   jj = jc+2-color
1198                   DO  k = nzb+2, nzt_mg(l), 2
1199                      j =jj
1200                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1201                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1202                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1203                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1204                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1205                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1206                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1207                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1208                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1209                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1210                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1211                                    + f2_mg(k,l) *                                                 &
1212                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1213                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1214                                    + f3_mg(k,l) *                                                 &
1215                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1216                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1217                                    - f_mg(k,j,i)                                                  &
1218                                      )
1219
1220                      j = jj+2
1221                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l)                                            &
1222                                    * ( rho_air_mg(k,l) * ddx2_mg(l) *                             &
1223                                        ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) *            &
1224                                                      ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +            &
1225                                          p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *            &
1226                                                      ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )            &
1227                                    + rho_air_mg(k,l) * ddy2_mg(l) *                               &
1228                                        ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) *            &
1229                                                      ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +            &
1230                                          p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *            &
1231                                                      ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )            &
1232                                    + f2_mg(k,l) *                                                 &
1233                                        ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) *            &
1234                                                      ( p_mg(k,j,i) - p_mg(k+1,j,i) ) )            &
1235                                    + f3_mg(k,l) *                                                 &
1236                                        ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) *            &
1237                                                      ( p_mg(k,j,i) - p_mg(k-1,j,i) ) )            &
1238                                    - f_mg(k,j,i)                                                  &
1239                                      )
1240                   ENDDO
1241
1242                ENDDO
1243             ENDDO
1244             CALL cpu_log( log_point_s(38), 'redblack_unroll_noopt', 'stop' )
1245
1246          ENDIF
1247
1248!
1249!--       Horizontal boundary conditions
1250          CALL exchange_horiz( p_mg, 1 )
1251
1252          IF ( .NOT. bc_lr_cyc )  THEN
1253             IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
1254                p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
1255             ENDIF
1256             IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
1257                p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
1258             ENDIF
1259          ENDIF
1260
1261          IF ( .NOT. bc_ns_cyc )  THEN
1262             IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
1263                p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
1264             ENDIF
1265             IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
1266                p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
1267             ENDIF
1268          ENDIF
1269
1270!
1271!--       Bottom and top boundary conditions
1272          IF ( ibc_p_b == 1 )  THEN
1273             p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
1274          ELSE
1275             p_mg(nzb,:,: ) = 0.0_wp
1276          ENDIF
1277
1278          IF ( ibc_p_t == 1 )  THEN
1279             p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
1280          ELSE
1281             p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
1282          ENDIF
1283
1284       ENDDO
1285
1286    ENDDO
1287
1288!
1289!-- Reset lower index limits to their standard values (may happen on coarsest levels only)
1290    IF ( adjust_lower_i_index )  THEN
1291       nxl_mg(l) = save_nxl_mg
1292    ENDIF
1293
1294    IF ( adjust_lower_j_index )  THEN
1295       nys_mg(l) = save_nys_mg
1296    ENDIF
1297
1298!
1299!-- Set pressure within topography and at the topography surfaces
1300!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
1301!$OMP DO
1302    DO  i = nxl_mg(l), nxr_mg(l)
1303       DO  j = nys_mg(l), nyn_mg(l)
1304          DO  k = nzb, nzt_mg(l)
1305!
1306!--          First, set pressure inside topography to zero
1307             p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
1308!
1309!--          Second, determine if the grid point inside topography is adjacent to a wall and set its
1310!--          value to a value given by the average of those values obtained from Neumann boundary
1311!--          condition.
1312             wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
1313             wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1314             wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1315             wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1316             wall_top   = IBITS( flags(k+1,j,i), 0, 1 )
1317             wall_total = wall_left + wall_right + wall_south + wall_north + wall_top
1318
1319             IF ( wall_total > 0.0_wp )  THEN
1320                p_mg(k,j,i) = 1.0_wp / wall_total * ( wall_left  * p_mg(k,j,i-1) +                 &
1321                                                      wall_right * p_mg(k,j,i+1) +                 &
1322                                                      wall_south * p_mg(k,j-1,i) +                 &
1323                                                      wall_north * p_mg(k,j+1,i) +                 &
1324                                                      wall_top   * p_mg(k+1,j,i) )
1325             ENDIF
1326          ENDDO
1327       ENDDO
1328    ENDDO
1329!$OMP END PARALLEL
1330
1331!
1332!-- One more time horizontal boundary conditions
1333    CALL exchange_horiz( p_mg, 1)
1334
1335
1336 END SUBROUTINE redblack_noopt
1337
1338
1339
1340!--------------------------------------------------------------------------------------------------!
1341! Description:
1342! ------------
1343!> Gather subdomain data from all PEs.
1344!--------------------------------------------------------------------------------------------------!
1345#if defined( __parallel )
1346 SUBROUTINE mg_gather_noopt( f2, f2_sub )
1347
1348    USE cpulog,                                                                                    &
1349        ONLY:  cpu_log,                                                                            &
1350               log_point_s
1351
1352    USE indices,                                                                                   &
1353        ONLY:  mg_loc_ind,                                                                         &
1354               nxl_mg,                                                                             &
1355               nxr_mg,                                                                             &
1356               nys_mg,                                                                             &
1357               nyn_mg,                                                                             &
1358               nzb,                                                                                &
1359               nzt_mg
1360
1361    USE kinds
1362
1363    USE pegrid
1364
1365    IMPLICIT NONE
1366
1367    INTEGER(iwp) ::  i       !<
1368    INTEGER(iwp) ::  il      !<
1369    INTEGER(iwp) ::  ir      !<
1370    INTEGER(iwp) ::  j       !<
1371    INTEGER(iwp) ::  jn      !<
1372    INTEGER(iwp) ::  js      !<
1373    INTEGER(iwp) ::  k       !<
1374    INTEGER(iwp) ::  nwords  !<
1375
1376    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
1377                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2    !<
1378    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
1379                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2_l  !<
1380
1381    REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,        &
1382                        mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub  !<
1383
1384
1385    CALL cpu_log( log_point_s(34), 'mg_gather_noopt', 'start' )
1386
1387    f2_l = 0.0_wp
1388
1389!
1390!-- Store the local subdomain array on the total array
1391    js = mg_loc_ind(3,myid)
1392    IF ( south_border_pe )  js = js - 1
1393    jn = mg_loc_ind(4,myid)
1394    IF ( north_border_pe )  jn = jn + 1
1395    il = mg_loc_ind(1,myid)
1396    IF ( left_border_pe )   il = il - 1
1397    ir = mg_loc_ind(2,myid)
1398    IF ( right_border_pe )  ir = ir + 1
1399    DO  i = il, ir
1400       DO  j = js, jn
1401          DO  k = nzb, nzt_mg(grid_level)+1
1402             f2_l(k,j,i) = f2_sub(k,j,i)
1403          ENDDO
1404       ENDDO
1405    ENDDO
1406
1407!
1408!-- Find out the number of array elements of the total array
1409    nwords = SIZE( f2 )
1410
1411!
1412!-- Gather subdomain data from all PEs
1413    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1414    CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),                       &
1415                        f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), nwords, MPI_REAL,       &
1416                        MPI_SUM, comm2d, ierr )
1417
1418    CALL cpu_log( log_point_s(34), 'mg_gather_noopt', 'stop' )
1419
1420 END SUBROUTINE mg_gather_noopt
1421#endif
1422
1423
1424!--------------------------------------------------------------------------------------------------!
1425! Description:
1426! ------------
1427!> @todo It may be possible to improve the speed of this routine by using non-blocking communication
1428!--------------------------------------------------------------------------------------------------!
1429#if defined( __parallel )
1430 SUBROUTINE mg_scatter_noopt( p2, p2_sub )
1431
1432    USE cpulog,                                                                                    &
1433        ONLY:  cpu_log,                                                                            &
1434               log_point_s
1435
1436    USE indices,                                                                                   &
1437        ONLY:  mg_loc_ind,                                                                         &
1438               nxl_mg,                                                                             &
1439               nxr_mg,                                                                             &
1440               nys_mg,                                                                             &
1441               nyn_mg,                                                                             &
1442               nzb,                                                                                &
1443               nzt_mg
1444
1445    USE kinds
1446
1447    USE pegrid
1448
1449    IMPLICIT NONE
1450
1451    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,  &
1452                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
1453
1454    REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,        &
1455                        mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub  !<
1456
1457
1458    CALL cpu_log( log_point_s(35), 'mg_scatter_noopt', 'start' )
1459
1460    p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,                                       &
1461                  mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1462
1463    CALL cpu_log( log_point_s(35), 'mg_scatter_noopt', 'stop' )
1464
1465 END SUBROUTINE mg_scatter_noopt
1466#endif
1467
1468
1469!--------------------------------------------------------------------------------------------------!
1470! Description:
1471! ------------
1472!> This is where the multigrid technique takes place. V- and W- Cycle are implemented and steered by
1473!> the parameter "gamma". Parameter "nue" determines the convergence of the multigrid iterative
1474!> solution. There are nue times RB-GS iterations. It should be set to "1" or "2", considering the
1475!> time effort one would like to invest. Last choice shows a very good converging factor, but leads
1476!> to an increase in computing time.
1477!--------------------------------------------------------------------------------------------------!
1478 RECURSIVE SUBROUTINE next_mg_level_noopt( f_mg, p_mg, p3, r )
1479
1480    USE control_parameters,                                                                        &
1481        ONLY:  bc_lr_dirrad,                                                                       &
1482               bc_lr_raddir,                                                                       &
1483               bc_ns_dirrad,                                                                       &
1484               bc_ns_raddir,                                                                       &
1485               gamma_mg,                                                                           &
1486               grid_level_count,                                                                   &
1487               maximum_grid_level,                                                                 &
1488               mg_switch_to_pe0_level,                                                             &
1489               mg_switch_to_pe0,                                                                   &
1490               ngsrb
1491
1492
1493    USE indices,                                                                                   &
1494        ONLY:  mg_loc_ind,                                                                         &
1495               nxl,                                                                                &
1496               nxl_mg,                                                                             &
1497               nxr,                                                                                &
1498               nxr_mg,                                                                             &
1499               nys,                                                                                &
1500               nys_mg,                                                                             &
1501               nyn,                                                                                &
1502               nyn_mg,                                                                             &
1503               nzb,                                                                                &
1504               nzt,                                                                                &
1505               nzt_mg
1506
1507    USE kinds
1508
1509    USE pegrid
1510
1511    IMPLICIT NONE
1512
1513    INTEGER(iwp) ::  i            !<
1514    INTEGER(iwp) ::  j            !<
1515    INTEGER(iwp) ::  k            !<
1516    INTEGER(iwp) ::  nxl_mg_save  !<
1517    INTEGER(iwp) ::  nxr_mg_save  !<
1518    INTEGER(iwp) ::  nyn_mg_save  !<
1519    INTEGER(iwp) ::  nys_mg_save  !<
1520    INTEGER(iwp) ::  nzt_mg_save  !<
1521
1522    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
1523                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg  !<
1524    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
1525                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg  !<
1526    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
1527                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3    !<
1528    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,        &
1529                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r     !<
1530
1531    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,  &
1532                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2  !<
1533    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,  &
1534                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
1535
1536    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub  !<
1537
1538#if defined( __parallel )
1539    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p2_sub  !<
1540#endif
1541
1542!
1543!-- Restriction to the coarsest grid
1544 10 IF ( grid_level == 1 )  THEN
1545
1546!
1547!--    Solution on the coarsest grid. Double the number of Gauss-Seidel iterations in order to get a
1548!--    more accurate solution.
1549       ngsrb = 2 * ngsrb
1550
1551       CALL redblack_noopt( f_mg, p_mg )
1552
1553       ngsrb = ngsrb / 2
1554
1555
1556    ELSEIF ( grid_level /= 1 )  THEN
1557
1558       grid_level_count(grid_level) = grid_level_count(grid_level) + 1
1559
1560!
1561!--    Solution on the actual grid level
1562       CALL redblack_noopt( f_mg, p_mg )
1563
1564!
1565!--    Determination of the actual residual
1566       CALL resid_noopt( f_mg, p_mg, r )
1567
1568!
1569!--    Restriction of the residual (finer grid values!) to the next coarser grid. Therefore, the
1570!--    grid level has to be decremented now. nxl..nzt have to be set to the coarse grid values,
1571!--    because these variables are needed for the exchange of ghost points in routine exchange_horiz
1572       grid_level = grid_level - 1
1573       nxl = nxl_mg(grid_level)
1574       nys = nys_mg(grid_level)
1575       nxr = nxr_mg(grid_level)
1576       nyn = nyn_mg(grid_level)
1577       nzt = nzt_mg(grid_level)
1578
1579       IF ( grid_level == mg_switch_to_pe0_level )  THEN
1580
1581!
1582!--       From this level on, calculations are done on PE0 only. First, carry out restriction on the
1583!--       subdomain. Therefore, indices of the level have to be changed to subdomain values in
1584!--       between (otherwise, the restrict routine would expect the gathered array)
1585
1586          nxl_mg_save = nxl_mg(grid_level)
1587          nxr_mg_save = nxr_mg(grid_level)
1588          nys_mg_save = nys_mg(grid_level)
1589          nyn_mg_save = nyn_mg(grid_level)
1590          nzt_mg_save = nzt_mg(grid_level)
1591          nxl_mg(grid_level) = mg_loc_ind(1,myid)
1592          nxr_mg(grid_level) = mg_loc_ind(2,myid)
1593          nys_mg(grid_level) = mg_loc_ind(3,myid)
1594          nyn_mg(grid_level) = mg_loc_ind(4,myid)
1595          nzt_mg(grid_level) = mg_loc_ind(5,myid)
1596          nxl = mg_loc_ind(1,myid)
1597          nxr = mg_loc_ind(2,myid)
1598          nys = mg_loc_ind(3,myid)
1599          nyn = mg_loc_ind(4,myid)
1600          nzt = mg_loc_ind(5,myid)
1601
1602          ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1,     &
1603                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
1604
1605          CALL restrict_noopt( f2_sub, r )
1606
1607!
1608!--       Restore the correct indices of this level
1609          nxl_mg(grid_level) = nxl_mg_save
1610          nxr_mg(grid_level) = nxr_mg_save
1611          nys_mg(grid_level) = nys_mg_save
1612          nyn_mg(grid_level) = nyn_mg_save
1613          nzt_mg(grid_level) = nzt_mg_save
1614          nxl = nxl_mg(grid_level)
1615          nxr = nxr_mg(grid_level)
1616          nys = nys_mg(grid_level)
1617          nyn = nyn_mg(grid_level)
1618          nzt = nzt_mg(grid_level)
1619!
1620!--       Gather all arrays from the subdomains on PE0
1621#if defined( __parallel )
1622          CALL mg_gather_noopt( f2, f2_sub )
1623#endif
1624
1625!
1626!--       Set switch for routine exchange_horiz, that no ghostpoint exchange has to be carried out
1627!--       from now on
1628          mg_switch_to_pe0 = .TRUE.
1629
1630!
1631!--       In case of non-cyclic lateral boundary conditions, both in- and outflow conditions have to
1632!--       be used on all PEs after the switch, because then they have the total domain.
1633          IF ( bc_lr_dirrad )  THEN
1634             bc_dirichlet_l  = .TRUE.
1635             bc_dirichlet_r  = .FALSE.
1636             bc_radiation_l = .FALSE.
1637             bc_radiation_r = .TRUE.
1638          ELSEIF ( bc_lr_raddir )  THEN
1639             bc_dirichlet_l  = .FALSE.
1640             bc_dirichlet_r  = .TRUE.
1641             bc_radiation_l = .TRUE.
1642             bc_radiation_r = .FALSE.
1643          ELSEIF ( child_domain  .OR.  nesting_offline )  THEN
1644             bc_dirichlet_l = .TRUE.
1645             bc_dirichlet_r = .TRUE.
1646          ENDIF
1647
1648          IF ( bc_ns_dirrad )  THEN
1649             bc_dirichlet_n  = .TRUE.
1650             bc_dirichlet_s  = .FALSE.
1651             bc_radiation_n = .FALSE.
1652             bc_radiation_s = .TRUE.
1653          ELSEIF ( bc_ns_raddir )  THEN
1654             bc_dirichlet_n  = .FALSE.
1655             bc_dirichlet_s  = .TRUE.
1656             bc_radiation_n = .TRUE.
1657             bc_radiation_s = .FALSE.
1658          ELSEIF ( child_domain  .OR.  nesting_offline )  THEN
1659             bc_dirichlet_s = .TRUE.
1660             bc_dirichlet_n = .TRUE.
1661          ENDIF
1662
1663          DEALLOCATE( f2_sub )
1664
1665       ELSE
1666
1667          CALL restrict_noopt( f2, r )
1668
1669       ENDIF
1670
1671       p2 = 0.0_wp
1672
1673!
1674!--    Repeat the same procedure till the coarsest grid is reached
1675       CALL next_mg_level_noopt( f2, p2, p3, r )
1676
1677    ENDIF
1678
1679!
1680!-- Now follows the prolongation
1681    IF ( grid_level >= 2 )  THEN
1682
1683!
1684!--    Prolongation of the new residual. The values are transferred from the coarse to the next
1685!--    finer grid.
1686       IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
1687
1688#if defined( __parallel )
1689!
1690!--       At this level, the new residual first has to be scattered from PE0 to the other PEs
1691          ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,     &
1692                                                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
1693
1694          CALL mg_scatter_noopt( p2, p2_sub )
1695
1696!
1697!--       Therefore, indices of the previous level have to be changed to subdomain values in between
1698!--       (otherwise, the prolong routine would expect the gathered array)
1699          nxl_mg_save = nxl_mg(grid_level-1)
1700          nxr_mg_save = nxr_mg(grid_level-1)
1701          nys_mg_save = nys_mg(grid_level-1)
1702          nyn_mg_save = nyn_mg(grid_level-1)
1703          nzt_mg_save = nzt_mg(grid_level-1)
1704          nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1705          nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1706          nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1707          nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1708          nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1709
1710!
1711!--       Set switch for routine exchange_horiz, that ghostpoint exchange has to be carried again
1712!--       out from now on.
1713          mg_switch_to_pe0 = .FALSE.
1714
1715!
1716!--       For non-cyclic lateral boundary conditions and in case of nesting, restore the in-/outflow
1717!--       conditions.
1718          bc_dirichlet_l = .FALSE.;  bc_dirichlet_r = .FALSE.
1719          bc_dirichlet_n = .FALSE.;  bc_dirichlet_s = .FALSE.
1720          bc_radiation_l = .FALSE.;  bc_radiation_r = .FALSE.
1721          bc_radiation_n = .FALSE.;  bc_radiation_s = .FALSE.
1722
1723          IF ( pleft == MPI_PROC_NULL )  THEN
1724             IF ( bc_lr_dirrad  .OR.  child_domain  .OR.  nesting_offline )  THEN
1725                bc_dirichlet_l = .TRUE.
1726             ELSEIF ( bc_lr_raddir )  THEN
1727                bc_radiation_l = .TRUE.
1728             ENDIF
1729          ENDIF
1730
1731          IF ( pright == MPI_PROC_NULL )  THEN
1732             IF ( bc_lr_dirrad )  THEN
1733                bc_radiation_r = .TRUE.
1734             ELSEIF ( bc_lr_raddir  .OR.  child_domain  .OR.  nesting_offline )  THEN
1735                bc_dirichlet_r = .TRUE.
1736             ENDIF
1737          ENDIF
1738
1739          IF ( psouth == MPI_PROC_NULL )  THEN
1740             IF ( bc_ns_dirrad )  THEN
1741                bc_radiation_s = .TRUE.
1742             ELSEIF ( bc_ns_raddir  .OR.  child_domain  .OR.  nesting_offline )  THEN
1743                bc_dirichlet_s = .TRUE.
1744             ENDIF
1745          ENDIF
1746
1747          IF ( pnorth == MPI_PROC_NULL )  THEN
1748             IF ( bc_ns_dirrad  .OR.  child_domain  .OR.  nesting_offline )  THEN
1749                bc_dirichlet_n = .TRUE.
1750             ELSEIF ( bc_ns_raddir )  THEN
1751                bc_radiation_n = .TRUE.
1752             ENDIF
1753          ENDIF
1754
1755          CALL prolong_noopt( p2_sub, p3 )
1756
1757!
1758!--       Restore the correct indices of the previous level
1759          nxl_mg(grid_level-1) = nxl_mg_save
1760          nxr_mg(grid_level-1) = nxr_mg_save
1761          nys_mg(grid_level-1) = nys_mg_save
1762          nyn_mg(grid_level-1) = nyn_mg_save
1763          nzt_mg(grid_level-1) = nzt_mg_save
1764
1765          DEALLOCATE( p2_sub )
1766#endif
1767
1768       ELSE
1769
1770          CALL prolong_noopt( p2, p3 )
1771
1772       ENDIF
1773
1774!
1775!--    Computation of the new pressure correction. Therefore, values from prior grids are added up
1776!--    automatically stage by stage.
1777       DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1778          DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1779             DO  k = nzb, nzt_mg(grid_level)+1
1780                p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1781             ENDDO
1782          ENDDO
1783       ENDDO
1784
1785!
1786!--    Relaxation of the new solution
1787       CALL redblack_noopt( f_mg, p_mg )
1788
1789    ENDIF
1790
1791
1792!
1793!-- The following few lines serve the steering of the multigrid scheme
1794    IF ( grid_level == maximum_grid_level )  THEN
1795
1796       GOTO 20
1797
1798    ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND.                       &
1799             grid_level_count(grid_level) /= gamma_mg )  THEN
1800
1801       GOTO 10
1802
1803    ENDIF
1804
1805!
1806!-- Reset counter for the next call of poismg_noopt
1807    grid_level_count(grid_level) = 0
1808
1809!
1810!-- Continue with the next finer level. nxl..nzt have to be set to the finer grid values, because
1811!-- these variables are needed for the exchange of ghost points in routine exchange_horiz.
1812    grid_level = grid_level + 1
1813    nxl = nxl_mg(grid_level)
1814    nxr = nxr_mg(grid_level)
1815    nys = nys_mg(grid_level)
1816    nyn = nyn_mg(grid_level)
1817    nzt = nzt_mg(grid_level)
1818
1819 20 CONTINUE
1820
1821 END SUBROUTINE next_mg_level_noopt
1822
1823
1824
1825!--------------------------------------------------------------------------------------------------!
1826! Description:
1827! ------------
1828!> @TODO: Missing subroutine description
1829!--------------------------------------------------------------------------------------------------!
1830 SUBROUTINE poismg_noopt_init
1831
1832    USE control_parameters,                                                                        &
1833        ONLY:  bc_lr_cyc,                                                                          &
1834               bc_ns_cyc,                                                                          &
1835               masking_method,                                                                     &
1836               maximum_grid_level,                                                                 &
1837               psolver
1838
1839    USE indices,                                                                                   &
1840        ONLY:  flags,                                                                              &
1841               nxl_mg,                                                                             &
1842               nxr_mg,                                                                             &
1843               nyn_mg,                                                                             &
1844               nys_mg,                                                                             &
1845               nzb,                                                                                &
1846               nzt_mg,                                                                             &
1847               wall_flags_total_0,                                                                 &
1848               wall_flags_1,                                                                       &
1849               wall_flags_10,                                                                      &
1850               wall_flags_2,                                                                       &
1851               wall_flags_3,                                                                       &
1852               wall_flags_4,                                                                       &
1853               wall_flags_5,                                                                       &
1854               wall_flags_6,                                                                       &
1855               wall_flags_7,                                                                       &
1856               wall_flags_8,                                                                       &
1857               wall_flags_9
1858
1859    IMPLICIT NONE
1860
1861    INTEGER(iwp) ::  i      !< index variable along x
1862    INTEGER(iwp) ::  inc    !< incremental parameter for coarsening grid level
1863    INTEGER(iwp) ::  j      !< index variable along y
1864    INTEGER(iwp) ::  k      !< index variable along z
1865    INTEGER(iwp) ::  l      !< loop variable indication current grid level
1866    INTEGER(iwp) ::  nxl_l  !< index of left PE boundary for multigrid level
1867    INTEGER(iwp) ::  nxr_l  !< index of right PE boundary for multigrid level
1868    INTEGER(iwp) ::  nyn_l  !< index of north PE boundary for multigrid level
1869    INTEGER(iwp) ::  nys_l  !< index of south PE boundary for multigrid level
1870    INTEGER(iwp) ::  nzt_l  !< index of top PE boundary for multigrid level
1871
1872    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo_tmp  !<
1873
1874    IF ( psolver /= 'multigrid_noopt' )  RETURN
1875!
1876!-- Grid point increment of the current level.
1877    inc = 1
1878    DO  l = maximum_grid_level, 1 , -1
1879!
1880!--    Set grid_level as it is required for exchange_horiz_2d_int
1881       grid_level = l
1882
1883       nxl_l = nxl_mg(l)
1884       nxr_l = nxr_mg(l)
1885       nys_l = nys_mg(l)
1886       nyn_l = nyn_mg(l)
1887       nzt_l = nzt_mg(l)
1888!
1889!--    Assign the flag level to be calculated
1890       SELECT CASE ( l )
1891          CASE ( 1 )
1892             flags => wall_flags_1
1893          CASE ( 2 )
1894             flags => wall_flags_2
1895          CASE ( 3 )
1896             flags => wall_flags_3
1897          CASE ( 4 )
1898             flags => wall_flags_4
1899          CASE ( 5 )
1900             flags => wall_flags_5
1901          CASE ( 6 )
1902             flags => wall_flags_6
1903          CASE ( 7 )
1904             flags => wall_flags_7
1905          CASE ( 8 )
1906             flags => wall_flags_8
1907          CASE ( 9 )
1908             flags => wall_flags_9
1909          CASE ( 10 )
1910             flags => wall_flags_10
1911       END SELECT
1912
1913!
1914!--    Depending on the grid level, set the respective bits in case of neighbouring walls
1915!--    Bit 0:  wall to the bottom
1916!--    Bit 1:  wall to the top (not realized in remaining PALM code so far)
1917!--    Bit 2:  wall to the south
1918!--    Bit 3:  wall to the north
1919!--    Bit 4:  wall to the left
1920!--    Bit 5:  wall to the right
1921!--    Bit 6:  inside building
1922
1923       flags = 0
1924!
1925!--    In case of masking method, flags are not set and multigrid method works like FFT-solver
1926       IF ( .NOT. masking_method )  THEN
1927
1928!
1929!--       Allocate temporary array for topography heights on coarser grid level. Please note, 2
1930!--       ghost points are required, in order to calculate flags() on the interior ghost point.
1931          ALLOCATE( topo_tmp(nzb:nzt_l+1,nys_l-1:nyn_l+1,nxl_l-1:nxr_l+1) )
1932          topo_tmp = 0
1933
1934          DO  i = nxl_l, nxr_l
1935             DO  j = nys_l, nyn_l
1936                DO  k = nzb, nzt_l
1937                   topo_tmp(k,j,i) = wall_flags_total_0(k*inc,j*inc,i*inc)
1938                ENDDO
1939             ENDDO
1940          ENDDO
1941          topo_tmp(nzt_l+1,:,:) = topo_tmp(nzt_l,:,:)
1942!
1943!--       Exchange ghost points on respective multigrid level. 2 ghost points are required, in order
1944!--       to calculate flags on.
1945!--       nys_l-1 / nyn_l+1 / nxl_l-1 / nxr_l+1.
1946          CALL exchange_horiz_int( topo_tmp, nys_l, nyn_l, nxl_l, nxr_l, nzt_l, 1 )
1947!
1948!--       Set non-cyclic boundary conditions on respective multigrid level
1949          IF ( .NOT. bc_ns_cyc )  THEN
1950             IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
1951!                 topo_tmp(:,-2,:) = topo_tmp(:,0,:)
1952                topo_tmp(:,-1,:) = topo_tmp(:,0,:)
1953             ENDIF
1954             IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
1955!                 topo_tmp(:,nyn_l+2,:) = topo_tmp(:,nyn_l,:)
1956                topo_tmp(:,nyn_l+1,:) = topo_tmp(:,nyn_l,:)
1957             ENDIF
1958          ENDIF
1959          IF ( .NOT. bc_lr_cyc )  THEN
1960             IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
1961!                 topo_tmp(:,:,-2) = topo_tmp(:,:,0)
1962                topo_tmp(:,:,-1) = topo_tmp(:,:,0)
1963             ENDIF
1964             IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
1965!                 topo_tmp(:,:,nxr_l+2) = topo_tmp(:,:,nxr_l)
1966                topo_tmp(:,:,nxr_l+1) = topo_tmp(:,:,nxr_l)
1967             ENDIF
1968          ENDIF
1969
1970          DO  i = nxl_l, nxr_l
1971             DO  j = nys_l, nyn_l
1972                DO  k = nzb, nzt_l
1973!
1974!--                Inside/outside building (inside building does not need further tests for walls)
1975                   IF ( .NOT. BTEST( topo_tmp(k,j,i), 0 ) )  THEN
1976
1977                      flags(k,j,i) = IBSET( flags(k,j,i), 6 )
1978
1979                   ELSE
1980!
1981!--                   Bottom wall
1982                      IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) )  THEN
1983                         flags(k,j,i) = IBSET( flags(k,j,i), 0 )
1984                      ENDIF
1985!
1986!--                   South wall
1987                      IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) )  THEN
1988                         flags(k,j,i) = IBSET( flags(k,j,i), 2 )
1989                      ENDIF
1990!
1991!--                   North wall
1992                      IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) )  THEN
1993                         flags(k,j,i) = IBSET( flags(k,j,i), 3 )
1994                      ENDIF
1995!
1996!--                   Left wall
1997                      IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) )  THEN
1998                         flags(k,j,i) = IBSET( flags(k,j,i), 4 )
1999                      ENDIF
2000!
2001!--                   Right wall
2002                      IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) )  THEN
2003                         flags(k,j,i) = IBSET( flags(k,j,i), 5 )
2004                      ENDIF
2005!
2006!--                   Top wall
2007                      IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) )  THEN
2008                         flags(k,j,i) = IBSET( flags(k,j,i), 7 )
2009                      ENDIF
2010
2011                   ENDIF
2012
2013                ENDDO
2014             ENDDO
2015          ENDDO
2016          flags(nzt_l+1,:,:) = flags(nzt_l,:,:)
2017
2018          CALL exchange_horiz_int( flags, nys_l, nyn_l, nxl_l, nxr_l, nzt_l, 1 )
2019
2020          DEALLOCATE( topo_tmp )
2021
2022       ENDIF
2023
2024       inc = inc * 2
2025
2026    ENDDO
2027!
2028!-- Reset grid_level to "normal" grid
2029    grid_level = 0
2030
2031 END SUBROUTINE poismg_noopt_init
2032
2033 END MODULE poismg_noopt_mod
Note: See TracBrowser for help on using the repository browser.