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

Last change on this file since 1160 was 1160, checked in by fricke, 8 years ago

last commit documented

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