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

Last change on this file since 4784 was 4649, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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