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

Last change on this file since 4421 was 4414, checked in by suehring, 4 years ago

Remove deprecated topography arrays; Move basic initialization of numerics into an extra module interface

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