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

Last change on this file since 1318 was 1318, checked in by raasch, 8 years ago

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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