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

Last change on this file since 4628 was 4457, checked in by raasch, 5 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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