source: palm/trunk/SOURCE/poismg_noopt_mod.f90 @ 3083

Last change on this file since 3083 was 2939, checked in by suehring, 6 years ago

Set lateral boundary conditions for divergence in multigrid solver

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