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

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

last commit documented

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