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

Last change on this file since 2699 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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