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

Last change on this file since 4401 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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