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

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

code has been put under the GNU General Public License (v3)

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