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

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

last commit documented

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