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

Last change on this file since 2716 was 2716, checked in by kanani, 3 years ago

Correction of "Former revisions" section

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