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

Last change on this file since 4429 was 4429, checked in by raasch, 3 years ago

serial (non-MPI) test case added, several bugfixes for the serial mode

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