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

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