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

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

unused variables remove from several routines

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