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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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