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

Last change on this file since 1093 was 1093, checked in by raasch, 9 years ago

last commit documented

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