source: palm/trunk/SOURCE/poismg.f90 @ 1056

Last change on this file since 1056 was 1056, checked in by raasch, 10 years ago

bugfix for mg-solver

  • Property svn:keywords set to Id
File size: 56.3 KB
Line 
1 SUBROUTINE poismg( r )
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Attention: Loop unrolling and cache optimization in SOR-Red/Black method
21!            still does not give the expected speedup! Further work required.
22!
23! Current revisions:
24! -----------------
25! Bugfix: all ghost points have to be used for allocating p3
26! arrays p2, f2, and f2_l changed from allocatable to automatic
27!
28! Former revisions:
29! -----------------
30! $Id: poismg.f90 1056 2012-11-16 15:28:04Z raasch $
31!
32! 1036 2012-10-22 13:43:42Z raasch
33! code put under GPL (PALM 3.9)
34!
35! 996 2012-09-07 10:41:47Z raasch
36! little reformatting
37!
38! 978 2012-08-09 08:28:32Z fricke
39! bc_lr/ns_dirneu/neudir added
40!
41! 880 2012-04-13 06:28:59Z raasch
42! Bugfix: preprocessor statements for parallel execution added
43!
44! 778 2011-11-07 14:18:25Z fricke
45! Allocation of p3 changes when multigrid is used and the collected field on PE0
46! has more grid points than the subdomain of an PE.
47!
48! 707 2011-03-29 11:39:40Z raasch
49! p_loc is used instead of p in the main routine (poismg).
50! On coarse grid levels, gathered data are identically processed on all PEs
51! (before, on PE0 only), so that the subsequent scattering of data is not
52! neccessary any more.
53! bc_lr/ns replaced by bc_lr/ns_cyc/dirrad/raddir
54! Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
55! resid and restrict. They were missed before which may have led to
56! unpredictable results.
57!
58! 667 2010-12-23 12:06:00Z suehring/gryschka
59! Calls of exchange_horiz are modified.
60!
61! 622 2010-12-10 08:08:13Z raasch
62! optional barriers included in order to speed up collective operations
63!
64! 257 2009-03-11 15:17:42Z heinze
65! Output of messages replaced by message handling routine.
66!
67! 181 2008-07-30 07:07:47Z raasch
68! Bugfix: grid_level+1 has to be used in restrict for flags-array
69!
70! 114 2007-10-10 00:03:15Z raasch
71! Boundary conditions at walls are implicitly set using flag arrays. Only
72! Neumann BC is allowed. Upper walls are still not realized.
73! Bottom and top BCs for array f_mg in restrict removed because boundary
74! values are not needed (right hand side of SOR iteration).
75!
76! 75 2007-03-22 09:54:05Z raasch
77! 2nd+3rd argument removed from exchange horiz
78!
79! RCS Log replace by Id keyword, revision history cleaned up
80!
81! Revision 1.6  2005/03/26 20:55:54  raasch
82! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
83! routine prolong simplified (one call of exchange_horiz spared)
84!
85! Revision 1.1  2001/07/20 13:10:51  raasch
86! Initial revision
87!
88!
89! Description:
90! ------------
91! Solves the Poisson equation for the perturbation pressure with a multigrid
92! V- or W-Cycle scheme.
93!
94! This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
95! September 2000 - July 2001.
96!------------------------------------------------------------------------------!
97
98    USE arrays_3d
99    USE control_parameters
100    USE cpulog   
101    USE grid_variables
102    USE indices
103    USE interfaces
104    USE pegrid
105
106    IMPLICIT NONE
107
108    REAL    ::  maxerror, maximum_mgcycles, residual_norm
109
110    REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r
111
112    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  p3
113
114
115    CALL cpu_log( log_point_s(29), 'poismg', 'start' )
116!
117!-- Initialize arrays and variables used in this subroutine
118
119!-- If the number of grid points of the gathered grid, which is collected
120!-- on PE0, is larger than the number of grid points of an PE, than array
121!-- p3 will be enlarged.
122    IF ( gathered_size > subdomain_size )  THEN
123       ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(               &
124                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
125                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
126                    mg_switch_to_pe0_level)+1) )
127    ELSE
128       ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
129    ENDIF
130
131    p3 = 0.0
132 
133!
134!-- Ghost boundaries have to be added to divergence array.
135!-- Exchange routine needs to know the grid level!
136    grid_level = maximum_grid_level
137    CALL exchange_horiz( d, 1)
138    d(nzb,:,:) = d(nzb+1,:,:)
139
140!
141!-- Initiation of the multigrid scheme. Does n cycles until the
142!-- residual is smaller than the given limit. The accuracy of the solution
143!-- of the poisson equation will increase with the number of cycles.
144!-- If the number of cycles is preset by the user, this number will be
145!-- carried out regardless of the accuracy.
146    grid_level_count =  0
147    mgcycles         =  0
148    IF ( mg_cycles == -1 )  THEN
149       maximum_mgcycles = 0
150       residual_norm    = 1.0 
151    ELSE
152       maximum_mgcycles = mg_cycles
153       residual_norm    = 0.0
154    ENDIF
155
156    DO WHILE ( residual_norm > residual_limit  .OR. &
157               mgcycles < maximum_mgcycles )
158 
159       CALL next_mg_level( d, p_loc, p3, r)
160
161!
162!--    Calculate the residual if the user has not preset the number of
163!--    cycles to be performed
164       IF ( maximum_mgcycles == 0 )  THEN
165          CALL resid( d, p_loc, r )
166          maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
167
168#if defined( __parallel )
169          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
170             CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
171                              comm2d, ierr)
172#else
173             residual_norm = maxerror
174#endif
175          residual_norm = SQRT( residual_norm )
176       ENDIF
177
178       mgcycles = mgcycles + 1
179
180!
181!--    If the user has not limited the number of cycles, stop the run in case
182!--    of insufficient convergence
183       IF ( mgcycles > 1000  .AND.  mg_cycles == -1 )  THEN
184          message_string = 'no sufficient convergence within 1000 cycles'
185          CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 )
186       ENDIF
187
188    ENDDO
189
190    DEALLOCATE( p3 )
191
192!
193!-- Unset the grid level. Variable is used to determine the MPI datatypes for
194!-- ghost point exchange
195    grid_level = 0 
196
197    CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
198
199 END SUBROUTINE poismg
200
201
202
203 SUBROUTINE resid( f_mg, p_mg, r )
204
205!------------------------------------------------------------------------------!
206! Description:
207! ------------
208! Computes the residual of the perturbation pressure.
209!------------------------------------------------------------------------------!
210
211    USE arrays_3d
212    USE control_parameters
213    USE grid_variables
214    USE indices
215    USE pegrid
216
217    IMPLICIT NONE
218
219    INTEGER ::  i, j, k, l
220
221    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                  &
222                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                 &
223                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg, p_mg, r
224
225!
226!-- Calculate the residual
227    l = grid_level
228
229!
230!-- Choose flag array of this level
231    SELECT CASE ( l )
232       CASE ( 1 )
233          flags => wall_flags_1
234       CASE ( 2 )
235          flags => wall_flags_2
236       CASE ( 3 )
237          flags => wall_flags_3
238       CASE ( 4 )
239          flags => wall_flags_4
240       CASE ( 5 )
241          flags => wall_flags_5
242       CASE ( 6 )
243          flags => wall_flags_6
244       CASE ( 7 )
245          flags => wall_flags_7
246       CASE ( 8 )
247          flags => wall_flags_8
248       CASE ( 9 )
249          flags => wall_flags_9
250       CASE ( 10 )
251          flags => wall_flags_10
252    END SELECT
253
254!$OMP PARALLEL PRIVATE (i,j,k)
255!$OMP DO
256    DO  i = nxl_mg(l), nxr_mg(l)
257       DO  j = nys_mg(l), nyn_mg(l) 
258          DO  k = nzb+1, nzt_mg(l)
259             r(k,j,i) = f_mg(k,j,i)                                         &
260                        - ddx2_mg(l) *                                      &
261                            ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
262                                          ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
263                              p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
264                                          ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
265                        - ddy2_mg(l) *                                      &
266                            ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
267                                          ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
268                              p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
269                                          ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
270                        - f2_mg(k,l) * p_mg(k+1,j,i)                        &
271                        - f3_mg(k,l) *                                      &
272                            ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
273                                          ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
274                        + f1_mg(k,l) * p_mg(k,j,i)
275!
276!--          Residual within topography should be zero
277             r(k,j,i) = r(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) )
278          ENDDO
279       ENDDO
280    ENDDO
281!$OMP END PARALLEL
282
283!
284!-- Horizontal boundary conditions
285    CALL exchange_horiz( r, 1)
286
287    IF ( .NOT. bc_lr_cyc )  THEN
288       IF ( inflow_l .OR. outflow_l )  r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
289       IF ( inflow_r .OR. outflow_r )  r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
290    ENDIF
291
292    IF ( .NOT. bc_ns_cyc )  THEN
293       IF ( inflow_n .OR. outflow_n )  r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
294       IF ( inflow_s .OR. outflow_s )  r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
295    ENDIF
296
297!
298!-- Boundary conditions at bottom and top of the domain.
299!-- These points are not handled by the above loop. Points may be within
300!-- buildings, but that doesn't matter.
301    IF ( ibc_p_b == 1 )  THEN
302       r(nzb,:,: ) = r(nzb+1,:,:)
303    ELSE
304       r(nzb,:,: ) = 0.0
305    ENDIF
306
307    IF ( ibc_p_t == 1 )  THEN
308       r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
309    ELSE
310       r(nzt_mg(l)+1,:,: ) = 0.0
311    ENDIF
312
313
314 END SUBROUTINE resid
315
316
317
318 SUBROUTINE restrict( f_mg, r )
319
320!------------------------------------------------------------------------------!
321! Description:
322! ------------
323! Interpolates the residual on the next coarser grid with "full weighting"
324! scheme
325!------------------------------------------------------------------------------!
326
327    USE control_parameters
328    USE grid_variables
329    USE indices
330    USE pegrid
331
332    IMPLICIT NONE
333
334    INTEGER ::  i, ic, j, jc, k, kc, l
335
336    REAL ::  rkjim, rkjip, rkjmi, rkjmim, rkjmip, rkjpi, rkjpim, rkjpip,       &
337             rkmji, rkmjim, rkmjip, rkmjmi, rkmjmim, rkmjmip, rkmjpi, rkmjpim, &
338             rkmjpip
339
340    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
341                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
342                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg
343
344    REAL, DIMENSION(nzb:nzt_mg(grid_level+1)+1,                          &
345                    nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,       &
346                    nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r
347
348!
349!-- Interpolate the residual
350    l = grid_level
351
352!
353!-- Choose flag array of the upper level
354    SELECT CASE ( l+1 )
355       CASE ( 1 )
356          flags => wall_flags_1
357       CASE ( 2 )
358          flags => wall_flags_2
359       CASE ( 3 )
360          flags => wall_flags_3
361       CASE ( 4 )
362          flags => wall_flags_4
363       CASE ( 5 )
364          flags => wall_flags_5
365       CASE ( 6 )
366          flags => wall_flags_6
367       CASE ( 7 )
368          flags => wall_flags_7
369       CASE ( 8 )
370          flags => wall_flags_8
371       CASE ( 9 )
372          flags => wall_flags_9
373       CASE ( 10 )
374          flags => wall_flags_10
375    END SELECT
376   
377!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
378!$OMP DO
379    DO  ic = nxl_mg(l), nxr_mg(l)   
380       i = 2*ic
381       DO  jc = nys_mg(l), nyn_mg(l) 
382          j = 2*jc
383          DO  kc = nzb+1, nzt_mg(l)
384             k = 2*kc-1
385!
386!--          Use implicit Neumann BCs if the respective gridpoint is inside
387!--          the building
388             rkjim   = r(k,j,i-1)     + IBITS( flags(k,j,i-1), 6, 1 ) *     &
389                                        ( r(k,j,i) - r(k,j,i-1) )
390             rkjip   = r(k,j,i+1)     + IBITS( flags(k,j,i+1), 6, 1 ) *     &
391                                        ( r(k,j,i) - r(k,j,i+1) )
392             rkjpi   = r(k,j+1,i)     + IBITS( flags(k,j+1,i), 6, 1 ) *     &
393                                        ( r(k,j,i) - r(k,j+1,i) )
394             rkjmi   = r(k,j-1,i)     + IBITS( flags(k,j-1,i), 6, 1 ) *     &
395                                        ( r(k,j,i) - r(k,j-1,i) )
396             rkjmim  = r(k,j-1,i-1)   + IBITS( flags(k,j-1,i-1), 6, 1 ) *   &
397                                        ( r(k,j,i) - r(k,j-1,i-1) )
398             rkjpim  = r(k,j+1,i-1)   + IBITS( flags(k,j+1,i-1), 6, 1 ) *   &
399                                        ( r(k,j,i) - r(k,j+1,i-1) )
400             rkjmip  = r(k,j-1,i+1)   + IBITS( flags(k,j-1,i+1), 6, 1 ) *   &
401                                        ( r(k,j,i) - r(k,j-1,i+1) )
402             rkjpip  = r(k,j+1,i+1)   + IBITS( flags(k,j+1,i+1), 6, 1 ) *   &
403                                        ( r(k,j,i) - r(k,j+1,i+1) )
404             rkmji   = r(k-1,j,i)     + IBITS( flags(k-1,j,i), 6, 1 ) *     &
405                                        ( r(k,j,i) - r(k-1,j,i) )
406             rkmjim  = r(k-1,j,i-1)   + IBITS( flags(k-1,j,i-1), 6, 1 ) *   &
407                                        ( r(k,j,i) - r(k-1,j,i-1) )
408             rkmjip  = r(k-1,j,i+1)   + IBITS( flags(k-1,j,i+1), 6, 1 ) *   &
409                                        ( r(k,j,i) - r(k-1,j,i+1) )
410             rkmjpi  = r(k-1,j+1,i)   + IBITS( flags(k-1,j+1,i), 6, 1 ) *   &
411                                        ( r(k,j,i) - r(k-1,j+1,i) )
412             rkmjmi  = r(k-1,j-1,i)   + IBITS( flags(k-1,j-1,i), 6, 1 ) *   &
413                                        ( r(k,j,i) - r(k-1,j-1,i) )
414             rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * &
415                                        ( r(k,j,i) - r(k-1,j-1,i-1) )
416             rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * &
417                                        ( r(k,j,i) - r(k-1,j+1,i-1) )
418             rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * &
419                                        ( r(k,j,i) - r(k-1,j-1,i+1) )
420             rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * &
421                                        ( r(k,j,i) - r(k-1,j+1,i+1) )
422
423             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
424                              8.0 * r(k,j,i)                            &
425                            + 4.0 * ( rkjim   + rkjip   +               &
426                                      rkjpi   + rkjmi   )               &
427                            + 2.0 * ( rkjmim  + rkjpim  +               &
428                                      rkjmip  + rkjpip  )               &
429                            + 4.0 * rkmji                               &
430                            + 2.0 * ( rkmjim  + rkmjim  +               &
431                                      rkmjpi  + rkmjmi  )               &
432                            +       ( rkmjmim + rkmjpim +               &
433                                      rkmjmip + rkmjpip )               &
434                            + 4.0 * r(k+1,j,i)                          &
435                            + 2.0 * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
436                                      r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
437                            +       ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
438                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
439                                           ) 
440
441!             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
442!                              8.0 * r(k,j,i)                            &
443!                            + 4.0 * ( r(k,j,i-1)     + r(k,j,i+1)     + &
444!                                      r(k,j+1,i)     + r(k,j-1,i)     ) &
445!                            + 2.0 * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
446!                                      r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
447!                            + 4.0 * r(k-1,j,i)                          &
448!                            + 2.0 * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
449!                                      r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
450!                            +       ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
451!                                      r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
452!                            + 4.0 * r(k+1,j,i)                          &
453!                            + 2.0 * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
454!                                      r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
455!                            +       ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
456!                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
457!                                           )
458          ENDDO
459       ENDDO
460    ENDDO
461!$OMP END PARALLEL
462
463!
464!-- Horizontal boundary conditions
465    CALL exchange_horiz( f_mg, 1)
466
467    IF ( .NOT. bc_lr_cyc )  THEN
468       IF (inflow_l .OR. outflow_l)  f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
469       IF (inflow_r .OR. outflow_r)  f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
470    ENDIF
471
472    IF ( .NOT. bc_ns_cyc )  THEN
473       IF (inflow_n .OR. outflow_n)  f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
474       IF (inflow_s .OR. outflow_s)  f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
475    ENDIF
476
477!
478!-- Boundary conditions at bottom and top of the domain.
479!-- These points are not handled by the above loop. Points may be within
480!-- buildings, but that doesn't matter.
481    IF ( ibc_p_b == 1 )  THEN
482       f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
483    ELSE
484       f_mg(nzb,:,: ) = 0.0
485    ENDIF
486
487    IF ( ibc_p_t == 1 )  THEN
488       f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
489    ELSE
490       f_mg(nzt_mg(l)+1,:,: ) = 0.0
491    ENDIF
492
493
494END SUBROUTINE restrict
495
496
497
498 SUBROUTINE prolong( p, temp )
499
500!------------------------------------------------------------------------------!
501! Description:
502! ------------
503! Interpolates the correction of the perturbation pressure
504! to the next finer grid.
505!------------------------------------------------------------------------------!
506
507    USE control_parameters
508    USE pegrid
509    USE indices
510
511    IMPLICIT NONE
512
513    INTEGER ::  i, j, k, l
514
515    REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1,                           &
516                    nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,        &
517                    nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p
518
519    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
520                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
521                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  temp
522
523
524!
525!-- First, store elements of the coarser grid on the next finer grid
526    l = grid_level
527
528!$OMP PARALLEL PRIVATE (i,j,k)
529!$OMP DO
530    DO  i = nxl_mg(l-1), nxr_mg(l-1)
531       DO  j = nys_mg(l-1), nyn_mg(l-1)
532!CDIR NODEP
533          DO  k = nzb+1, nzt_mg(l-1)
534!
535!--          Points of the coarse grid are directly stored on the next finer
536!--          grid
537             temp(2*k-1,2*j,2*i) = p(k,j,i) 
538!
539!--          Points between two coarse-grid points
540             temp(2*k-1,2*j,2*i+1) = 0.5 * ( p(k,j,i) + p(k,j,i+1) )
541             temp(2*k-1,2*j+1,2*i) = 0.5 * ( p(k,j,i) + p(k,j+1,i) )
542             temp(2*k,2*j,2*i)     = 0.5 * ( p(k,j,i) + p(k+1,j,i) )
543!
544!--          Points in the center of the planes stretched by four points
545!--          of the coarse grid cube
546             temp(2*k-1,2*j+1,2*i+1) = 0.25 * ( p(k,j,i)   + p(k,j,i+1) + &
547                                                p(k,j+1,i) + p(k,j+1,i+1) )
548             temp(2*k,2*j,2*i+1)     = 0.25 * ( p(k,j,i)   + p(k,j,i+1) + &
549                                                p(k+1,j,i) + p(k+1,j,i+1) )
550             temp(2*k,2*j+1,2*i)     = 0.25 * ( p(k,j,i)   + p(k,j+1,i) + &
551                                                p(k+1,j,i) + p(k+1,j+1,i) )
552!
553!--          Points in the middle of coarse grid cube
554             temp(2*k,2*j+1,2*i+1) = 0.125 * ( p(k,j,i)     + p(k,j,i+1)   + &
555                                               p(k,j+1,i)   + p(k,j+1,i+1) + &
556                                               p(k+1,j,i)   + p(k+1,j,i+1) + &
557                                               p(k+1,j+1,i) + p(k+1,j+1,i+1) )
558          ENDDO
559       ENDDO
560    ENDDO
561!$OMP END PARALLEL
562                         
563!
564!-- Horizontal boundary conditions
565    CALL exchange_horiz( temp, 1)
566
567    IF ( .NOT. bc_lr_cyc )  THEN
568       IF (inflow_l .OR. outflow_l)  temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
569       IF (inflow_r .OR. outflow_r)  temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
570    ENDIF
571
572    IF ( .NOT. bc_ns_cyc )  THEN
573       IF (inflow_n .OR. outflow_n)  temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
574       IF (inflow_s .OR. outflow_s)  temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
575    ENDIF
576
577!
578!-- Bottom and top boundary conditions
579    IF ( ibc_p_b == 1 )  THEN
580       temp(nzb,:,: ) = temp(nzb+1,:,:)
581    ELSE
582       temp(nzb,:,: ) = 0.0
583    ENDIF
584
585    IF ( ibc_p_t == 1 )  THEN
586       temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
587    ELSE
588       temp(nzt_mg(l)+1,:,: ) = 0.0
589    ENDIF
590
591 
592 END SUBROUTINE prolong
593
594
595 SUBROUTINE redblack( f_mg, p_mg )
596
597!------------------------------------------------------------------------------!
598! Description:
599! ------------
600! Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
601! 3D-Red-Black decomposition (GS-RB) is used.
602!------------------------------------------------------------------------------!
603
604    USE arrays_3d
605    USE control_parameters
606    USE cpulog
607    USE grid_variables
608    USE indices
609    USE interfaces
610    USE pegrid
611
612    IMPLICIT NONE
613
614    INTEGER :: colour, i, ic, j, jc, jj, k, l, n
615
616    LOGICAL :: unroll
617
618    REAL ::  wall_left, wall_north, wall_right, wall_south, wall_total, wall_top
619
620    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                 &
621                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                &
622                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg, p_mg
623
624    l = grid_level
625
626!
627!-- Choose flag array of this level
628    SELECT CASE ( l )
629       CASE ( 1 )
630          flags => wall_flags_1
631       CASE ( 2 )
632          flags => wall_flags_2
633       CASE ( 3 )
634          flags => wall_flags_3
635       CASE ( 4 )
636          flags => wall_flags_4
637       CASE ( 5 )
638          flags => wall_flags_5
639       CASE ( 6 )
640          flags => wall_flags_6
641       CASE ( 7 )
642          flags => wall_flags_7
643       CASE ( 8 )
644          flags => wall_flags_8
645       CASE ( 9 )
646          flags => wall_flags_9
647       CASE ( 10 )
648          flags => wall_flags_10
649    END SELECT
650
651    unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
652               MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
653
654    DO  n = 1, ngsrb
655       
656       DO  colour = 1, 2
657
658          IF ( .NOT. unroll )  THEN
659
660             CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
661
662!
663!--          Without unrolling of loops, no cache optimization
664             DO  i = nxl_mg(l), nxr_mg(l), 2
665                DO  j = nys_mg(l) + 2 - colour, nyn_mg(l), 2 
666                   DO  k = nzb+1, nzt_mg(l), 2
667!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
668!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
669!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
670!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
671!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
672!                                                       )
673
674                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
675                             ddx2_mg(l) *                                      &
676                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
677                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
678                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
679                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
680                           + ddy2_mg(l) *                                      &
681                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
682                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
683                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
684                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
685                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
686                           + f3_mg(k,l) *                                      &
687                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
688                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
689                           - f_mg(k,j,i)               )
690                   ENDDO
691                ENDDO
692             ENDDO
693   
694             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
695                DO  j = nys_mg(l) + (colour-1), nyn_mg(l), 2
696                   DO  k = nzb+1, nzt_mg(l), 2 
697                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
698                             ddx2_mg(l) *                                      &
699                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
700                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
701                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
702                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
703                           + ddy2_mg(l) *                                      &
704                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
705                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
706                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
707                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
708                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
709                           + f3_mg(k,l) *                                      &
710                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
711                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
712                           - f_mg(k,j,i)               )
713                   ENDDO
714                ENDDO
715             ENDDO
716 
717             DO  i = nxl_mg(l), nxr_mg(l), 2
718                DO  j = nys_mg(l) + (colour-1), nyn_mg(l), 2
719                   DO  k = nzb+2, nzt_mg(l), 2
720                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
721                             ddx2_mg(l) *                                      &
722                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
723                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
724                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
725                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
726                           + ddy2_mg(l) *                                      &
727                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
728                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
729                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
730                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
731                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
732                           + f3_mg(k,l) *                                      &
733                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
734                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
735                           - f_mg(k,j,i)               )
736                   ENDDO
737                ENDDO
738             ENDDO
739
740             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
741                DO  j = nys_mg(l) + 2 - colour, nyn_mg(l), 2
742                   DO  k = nzb+2, nzt_mg(l), 2
743                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
744                             ddx2_mg(l) *                                      &
745                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
746                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
747                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
748                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
749                           + ddy2_mg(l) *                                      &
750                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
751                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
752                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
753                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
754                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
755                           + f3_mg(k,l) *                                      &
756                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
757                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
758                           - f_mg(k,j,i)               )
759                   ENDDO
760                ENDDO
761             ENDDO
762             CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
763
764          ELSE
765
766!
767!--          Loop unrolling along y, only one i loop for better cache use
768             CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
769             DO  ic = nxl_mg(l), nxr_mg(l), 2
770                DO  jc = nys_mg(l), nyn_mg(l), 4
771                   i  = ic
772                   jj = jc+2-colour
773                   DO  k = nzb+1, nzt_mg(l), 2
774                      j = jj
775                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
776                             ddx2_mg(l) *                                      &
777                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
778                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
779                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
780                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
781                           + ddy2_mg(l) *                                      &
782                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
783                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
784                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
785                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
786                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
787                           + f3_mg(k,l) *                                      &
788                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
789                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
790                           - f_mg(k,j,i)               )
791                      j = jj+2
792                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
793                             ddx2_mg(l) *                                      &
794                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
795                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
796                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
797                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
798                           + ddy2_mg(l) *                                      &
799                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
800                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
801                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
802                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
803                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
804                           + f3_mg(k,l) *                                      &
805                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
806                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
807                           - f_mg(k,j,i)               )
808                   ENDDO
809   
810                   i  = ic+1
811                   jj = jc+colour-1
812                   DO  k = nzb+1, nzt_mg(l), 2 
813                      j =jj
814                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
815                             ddx2_mg(l) *                                      &
816                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
817                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
818                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
819                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
820                           + ddy2_mg(l) *                                      &
821                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
822                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
823                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
824                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
825                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
826                           + f3_mg(k,l) *                                      &
827                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
828                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
829                           - f_mg(k,j,i)               )
830                      j = jj+2
831                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
832                             ddx2_mg(l) *                                      &
833                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
834                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
835                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
836                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
837                           + ddy2_mg(l) *                                      &
838                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
839                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
840                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
841                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
842                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
843                           + f3_mg(k,l) *                                      &
844                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
845                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
846                           - f_mg(k,j,i)               )
847                   ENDDO
848
849                   i  = ic
850                   jj = jc+colour-1
851                   DO  k = nzb+2, nzt_mg(l), 2
852                      j =jj
853                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
854                             ddx2_mg(l) *                                      &
855                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
856                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
857                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
858                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
859                           + ddy2_mg(l) *                                      &
860                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
861                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
862                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
863                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
864                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
865                           + f3_mg(k,l) *                                      &
866                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
867                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
868                           - f_mg(k,j,i)               )
869                      j = jj+2
870                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
871                             ddx2_mg(l) *                                      &
872                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
873                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
874                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
875                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
876                           + ddy2_mg(l) *                                      &
877                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
878                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
879                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
880                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
881                           + f2_mg(k,l) * 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
888                   i  = ic+1
889                   jj = jc+2-colour
890                   DO  k = nzb+2, nzt_mg(l), 2
891                      j =jj
892                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
893                             ddx2_mg(l) *                                      &
894                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
895                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
896                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
897                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
898                           + ddy2_mg(l) *                                      &
899                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
900                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
901                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
902                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
903                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
904                           + f3_mg(k,l) *                                      &
905                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
906                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
907                           - f_mg(k,j,i)               )
908                      j = jj+2
909                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
910                             ddx2_mg(l) *                                      &
911                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
912                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
913                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
914                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
915                           + ddy2_mg(l) *                                      &
916                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
917                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
918                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
919                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
920                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
921                           + f3_mg(k,l) *                                      &
922                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
923                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
924                           - f_mg(k,j,i)               )
925                   ENDDO
926
927                ENDDO
928             ENDDO
929             CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
930
931          ENDIF
932
933!
934!--       Horizontal boundary conditions
935          CALL exchange_horiz( p_mg, 1 )
936
937          IF ( .NOT. bc_lr_cyc )  THEN
938             IF ( inflow_l .OR. outflow_l )  THEN
939                p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
940             ENDIF
941             IF ( inflow_r .OR. outflow_r )  THEN
942                p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
943             ENDIF
944          ENDIF
945
946          IF ( .NOT. bc_ns_cyc )  THEN
947             IF ( inflow_n .OR. outflow_n )  THEN
948                p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
949             ENDIF
950             IF ( inflow_s .OR. outflow_s )  THEN
951                p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
952             ENDIF
953          ENDIF
954
955!
956!--       Bottom and top boundary conditions
957          IF ( ibc_p_b == 1 )  THEN
958             p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
959          ELSE
960             p_mg(nzb,:,: ) = 0.0
961          ENDIF
962
963          IF ( ibc_p_t == 1 )  THEN
964             p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
965          ELSE
966             p_mg(nzt_mg(l)+1,:,: ) = 0.0
967          ENDIF
968
969       ENDDO
970
971    ENDDO
972
973!
974!-- Set pressure within topography and at the topography surfaces
975!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
976!$OMP DO
977    DO  i = nxl_mg(l), nxr_mg(l)
978       DO  j = nys_mg(l), nyn_mg(l) 
979          DO  k = nzb, nzt_mg(l)
980!
981!--          First, set pressure inside topography to zero
982             p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) )
983!
984!--          Second, determine if the gridpoint inside topography is adjacent
985!--          to a wall and set its value to a value given by the average of
986!--          those values obtained from Neumann boundary condition
987             wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
988             wall_right = IBITS( flags(k,j,i+1), 4, 1 )
989             wall_south = IBITS( flags(k,j-1,i), 3, 1 )
990             wall_north = IBITS( flags(k,j+1,i), 2, 1 )
991             wall_top   = IBITS( flags(k+1,j,i), 0, 1 )
992             wall_total = wall_left + wall_right + wall_south + wall_north + &
993                          wall_top
994
995             IF ( wall_total > 0.0 )  THEN
996                p_mg(k,j,i) = 1.0 / wall_total *                 &
997                                  ( wall_left  * p_mg(k,j,i-1) + &
998                                    wall_right * p_mg(k,j,i+1) + &
999                                    wall_south * p_mg(k,j-1,i) + &
1000                                    wall_north * p_mg(k,j+1,i) + &
1001                                    wall_top   * p_mg(k+1,j,i) )
1002             ENDIF
1003          ENDDO
1004       ENDDO
1005    ENDDO
1006!$OMP END PARALLEL
1007
1008!
1009!-- One more time horizontal boundary conditions
1010    CALL exchange_horiz( p_mg, 1)
1011
1012
1013 END SUBROUTINE redblack
1014
1015
1016
1017 SUBROUTINE mg_gather( f2, f2_sub )
1018
1019    USE control_parameters
1020    USE cpulog
1021    USE indices
1022    USE interfaces
1023    USE pegrid
1024
1025    IMPLICIT NONE
1026
1027    INTEGER ::  i, il, ir, j, jn, js, k, n, nwords, sender
1028
1029    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
1030                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
1031                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2, f2_l
1032
1033    REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1,                            &
1034                    mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,           &
1035                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub
1036
1037
1038#if defined( __parallel )
1039    CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
1040
1041    f2_l = 0.0
1042
1043!
1044!-- Store the local subdomain array on the total array
1045    js = mg_loc_ind(3,myid)
1046    IF ( south_border_pe )  js = js - 1
1047    jn = mg_loc_ind(4,myid)
1048    IF ( north_border_pe )  jn = jn + 1
1049    il = mg_loc_ind(1,myid)
1050    IF ( left_border_pe )   il = il - 1
1051    ir = mg_loc_ind(2,myid)
1052    IF ( right_border_pe )  ir = ir + 1
1053    DO  i = il, ir
1054       DO  j = js, jn
1055          DO  k = nzb, nzt_mg(grid_level)+1
1056             f2_l(k,j,i) = f2_sub(k,j,i)
1057          ENDDO
1058       ENDDO
1059    ENDDO
1060
1061!
1062!-- Find out the number of array elements of the total array
1063    nwords = SIZE( f2 )
1064
1065!
1066!-- Gather subdomain data from all PEs
1067    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1068    CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1069                        f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),   &
1070                        nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
1071
1072    CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
1073#endif
1074   
1075 END SUBROUTINE mg_gather
1076
1077
1078
1079 SUBROUTINE mg_scatter( p2, p2_sub )
1080!
1081!-- TODO: It may be possible to improve the speed of this routine by using
1082!-- non-blocking communication
1083
1084    USE control_parameters
1085    USE cpulog
1086    USE indices
1087    USE interfaces
1088    USE pegrid
1089
1090    IMPLICIT NONE
1091
1092    INTEGER ::  n, nwords, sender
1093
1094    REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1,                            &
1095                    nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,         &
1096                    nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2
1097
1098    REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1,                              &
1099                    mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,             &
1100                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub
1101
1102!
1103!-- Find out the number of array elements of the subdomain array
1104    nwords = SIZE( p2_sub )
1105
1106#if defined( __parallel )
1107    CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
1108
1109    p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1110                  mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1111
1112    CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
1113#endif
1114   
1115 END SUBROUTINE mg_scatter
1116
1117
1118
1119 RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
1120
1121!------------------------------------------------------------------------------!
1122! Description:
1123! ------------
1124! This is where the multigrid technique takes place. V- and W- Cycle are
1125! implemented and steered by the parameter "gamma". Parameter "nue" determines
1126! the convergence of the multigrid iterative solution. There are nue times
1127! RB-GS iterations. It should be set to "1" or "2", considering the time effort
1128! one would like to invest. Last choice shows a very good converging factor,
1129! but leads to an increase in computing time.
1130!------------------------------------------------------------------------------!
1131
1132    USE arrays_3d
1133    USE control_parameters
1134    USE grid_variables
1135    USE indices
1136    USE pegrid
1137
1138    IMPLICIT NONE
1139
1140    INTEGER ::  i, j, k, nxl_mg_save, nxr_mg_save, nyn_mg_save, nys_mg_save, &
1141                nzt_mg_save
1142
1143    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                  &
1144                 nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                    &
1145                 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg, p3, r
1146
1147    REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1,                                &
1148                    nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,             &
1149                    nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2, p2
1150
1151    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub, p2_sub
1152
1153!
1154!-- Restriction to the coarsest grid
1155 10 IF ( grid_level == 1 )  THEN
1156
1157!
1158!--    Solution on the coarsest grid. Double the number of Gauss-Seidel
1159!--    iterations in order to get a more accurate solution.
1160       ngsrb = 2 * ngsrb
1161
1162       CALL redblack( f_mg, p_mg )
1163
1164       ngsrb = ngsrb / 2
1165
1166
1167    ELSEIF ( grid_level /= 1 )  THEN
1168
1169       grid_level_count(grid_level) = grid_level_count(grid_level) + 1
1170
1171!
1172!--    Solution on the actual grid level
1173       CALL redblack( f_mg, p_mg )
1174
1175!
1176!--    Determination of the actual residual
1177       CALL resid( f_mg, p_mg, r )
1178
1179!
1180!--    Restriction of the residual (finer grid values!) to the next coarser
1181!--    grid. Therefore, the grid level has to be decremented now. nxl..nzt have
1182!--    to be set to the coarse grid values, because these variables are needed
1183!--    for the exchange of ghost points in routine exchange_horiz
1184       grid_level = grid_level - 1
1185       nxl = nxl_mg(grid_level)
1186       nys = nys_mg(grid_level)
1187       nxr = nxr_mg(grid_level)
1188       nyn = nyn_mg(grid_level)
1189       nzt = nzt_mg(grid_level)
1190
1191       IF ( grid_level == mg_switch_to_pe0_level )  THEN
1192
1193!
1194!--       From this level on, calculations are done on PE0 only.
1195!--       First, carry out restriction on the subdomain.
1196!--       Therefore, indices of the level have to be changed to subdomain values
1197!--       in between (otherwise, the restrict routine would expect
1198!--       the gathered array)
1199
1200          nxl_mg_save = nxl_mg(grid_level)
1201          nxr_mg_save = nxr_mg(grid_level)
1202          nys_mg_save = nys_mg(grid_level)
1203          nyn_mg_save = nyn_mg(grid_level)
1204          nzt_mg_save = nzt_mg(grid_level)
1205          nxl_mg(grid_level) = mg_loc_ind(1,myid)
1206          nxr_mg(grid_level) = mg_loc_ind(2,myid)
1207          nys_mg(grid_level) = mg_loc_ind(3,myid)
1208          nyn_mg(grid_level) = mg_loc_ind(4,myid)
1209          nzt_mg(grid_level) = mg_loc_ind(5,myid)
1210          nxl = mg_loc_ind(1,myid)
1211          nxr = mg_loc_ind(2,myid)
1212          nys = mg_loc_ind(3,myid)
1213          nyn = mg_loc_ind(4,myid)
1214          nzt = mg_loc_ind(5,myid)
1215
1216          ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,                    &
1217                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
1218                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
1219
1220          CALL restrict( f2_sub, r )
1221
1222!
1223!--       Restore the correct indices of this level
1224          nxl_mg(grid_level) = nxl_mg_save
1225          nxr_mg(grid_level) = nxr_mg_save
1226          nys_mg(grid_level) = nys_mg_save
1227          nyn_mg(grid_level) = nyn_mg_save
1228          nzt_mg(grid_level) = nzt_mg_save
1229          nxl = nxl_mg(grid_level)
1230          nxr = nxr_mg(grid_level)
1231          nys = nys_mg(grid_level)
1232          nyn = nyn_mg(grid_level)
1233          nzt = nzt_mg(grid_level)
1234!
1235!--       Gather all arrays from the subdomains on PE0
1236          CALL mg_gather( f2, f2_sub )
1237
1238!
1239!--       Set switch for routine exchange_horiz, that no ghostpoint exchange
1240!--       has to be carried out from now on
1241          mg_switch_to_pe0 = .TRUE.
1242
1243!
1244!--       In case of non-cyclic lateral boundary conditions, both in- and
1245!--       outflow conditions have to be used on all PEs after the switch,
1246!--       because then they have the total domain.
1247          IF ( bc_lr_dirrad  .OR.  bc_lr_dirneu )  THEN
1248             inflow_l  = .TRUE.
1249             inflow_r  = .FALSE.
1250             outflow_l = .FALSE.
1251             outflow_r = .TRUE.
1252          ELSEIF ( bc_lr_raddir  .OR.  bc_lr_neudir )  THEN
1253             inflow_l  = .FALSE.
1254             inflow_r  = .TRUE.
1255             outflow_l = .TRUE.
1256             outflow_r = .FALSE.
1257          ENDIF
1258
1259          IF ( bc_ns_dirrad  .OR.  bc_ns_dirneu )  THEN
1260             inflow_n  = .TRUE.
1261             inflow_s  = .FALSE.
1262             outflow_n = .FALSE.
1263             outflow_s = .TRUE.
1264          ELSEIF ( bc_ns_raddir  .OR.  bc_ns_neudir )  THEN
1265             inflow_n  = .FALSE.
1266             inflow_s  = .TRUE.
1267             outflow_n = .TRUE.
1268             outflow_s = .FALSE.
1269          ENDIF
1270
1271          DEALLOCATE( f2_sub )
1272
1273       ELSE
1274
1275          CALL restrict( f2, r )
1276
1277       ENDIF
1278
1279       p2 = 0.0
1280
1281!
1282!--    Repeat the same procedure till the coarsest grid is reached
1283       CALL next_mg_level( f2, p2, p3, r )
1284
1285    ENDIF
1286
1287!
1288!-- Now follows the prolongation
1289    IF ( grid_level >= 2 )  THEN
1290
1291!
1292!--    Prolongation of the new residual. The values are transferred
1293!--    from the coarse to the next finer grid.
1294       IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
1295
1296#if defined( __parallel )
1297!
1298!--       At this level, the new residual first has to be scattered from
1299!--       PE0 to the other PEs
1300          ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,             &
1301                    mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,   &
1302                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
1303
1304          CALL mg_scatter( p2, p2_sub )
1305
1306!
1307!--       Therefore, indices of the previous level have to be changed to
1308!--       subdomain values in between (otherwise, the prolong routine would
1309!--       expect the gathered array)
1310          nxl_mg_save = nxl_mg(grid_level-1)
1311          nxr_mg_save = nxr_mg(grid_level-1)
1312          nys_mg_save = nys_mg(grid_level-1)
1313          nyn_mg_save = nyn_mg(grid_level-1)
1314          nzt_mg_save = nzt_mg(grid_level-1)
1315          nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1316          nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1317          nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1318          nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1319          nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1320
1321!
1322!--       Set switch for routine exchange_horiz, that ghostpoint exchange
1323!--       has to be carried again out from now on
1324          mg_switch_to_pe0 = .FALSE.
1325
1326!
1327!--       For non-cyclic lateral boundary conditions, restore the
1328!--       in-/outflow conditions
1329          inflow_l  = .FALSE.;  inflow_r  = .FALSE.
1330          inflow_n  = .FALSE.;  inflow_s  = .FALSE.
1331          outflow_l = .FALSE.;  outflow_r = .FALSE.
1332          outflow_n = .FALSE.;  outflow_s = .FALSE.
1333
1334          IF ( pleft == MPI_PROC_NULL )  THEN
1335             IF ( bc_lr_dirrad  .OR.  bc_lr_dirneu )  THEN
1336                inflow_l  = .TRUE.
1337             ELSEIF ( bc_lr_raddir  .OR.  bc_lr_neudir )  THEN
1338                outflow_l = .TRUE.
1339             ENDIF
1340          ENDIF
1341
1342          IF ( pright == MPI_PROC_NULL )  THEN
1343             IF ( bc_lr_dirrad  .OR.  bc_lr_dirneu )  THEN
1344                outflow_r = .TRUE.
1345             ELSEIF ( bc_lr_raddir  .OR.  bc_lr_neudir )  THEN
1346                inflow_r  = .TRUE.
1347             ENDIF
1348          ENDIF
1349
1350          IF ( psouth == MPI_PROC_NULL )  THEN
1351             IF ( bc_ns_dirrad  .OR.  bc_ns_dirneu )  THEN
1352                outflow_s = .TRUE.
1353             ELSEIF ( bc_ns_raddir  .OR.  bc_ns_neudir )  THEN
1354                inflow_s  = .TRUE.
1355             ENDIF
1356          ENDIF
1357
1358          IF ( pnorth == MPI_PROC_NULL )  THEN
1359             IF ( bc_ns_dirrad  .OR.  bc_ns_dirneu )  THEN
1360                inflow_n  = .TRUE.
1361             ELSEIF ( bc_ns_raddir  .OR.  bc_ns_neudir )  THEN
1362                outflow_n = .TRUE.
1363             ENDIF
1364          ENDIF
1365
1366          CALL prolong( p2_sub, p3 )
1367
1368!
1369!--       Restore the correct indices of the previous level
1370          nxl_mg(grid_level-1) = nxl_mg_save
1371          nxr_mg(grid_level-1) = nxr_mg_save
1372          nys_mg(grid_level-1) = nys_mg_save
1373          nyn_mg(grid_level-1) = nyn_mg_save
1374          nzt_mg(grid_level-1) = nzt_mg_save
1375
1376          DEALLOCATE( p2_sub )
1377#endif
1378
1379       ELSE
1380
1381          CALL prolong( p2, p3 )
1382
1383       ENDIF
1384
1385!
1386!--    Computation of the new pressure correction. Therefore,
1387!--    values from prior grids are added up automatically stage by stage.
1388       DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1389          DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1390             DO  k = nzb, nzt_mg(grid_level)+1 
1391                p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1392             ENDDO
1393          ENDDO
1394       ENDDO
1395
1396!
1397!--    Relaxation of the new solution
1398       CALL redblack( f_mg, p_mg )
1399
1400    ENDIF
1401
1402
1403!
1404!-- The following few lines serve the steering of the multigrid scheme
1405    IF ( grid_level == maximum_grid_level )  THEN
1406
1407       GOTO 20
1408
1409    ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND. &
1410             grid_level_count(grid_level) /= gamma_mg )  THEN
1411
1412       GOTO 10
1413
1414    ENDIF
1415
1416!
1417!-- Reset counter for the next call of poismg
1418    grid_level_count(grid_level) = 0
1419
1420!
1421!-- Continue with the next finer level. nxl..nzt have to be
1422!-- set to the finer grid values, because these variables are needed for the
1423!-- exchange of ghost points in routine exchange_horiz
1424    grid_level = grid_level + 1
1425    nxl = nxl_mg(grid_level)
1426    nxr = nxr_mg(grid_level)
1427    nys = nys_mg(grid_level)
1428    nyn = nyn_mg(grid_level)
1429    nzt = nzt_mg(grid_level)
1430
1431 20 CONTINUE
1432
1433 END SUBROUTINE next_mg_level
Note: See TracBrowser for help on using the repository browser.