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

Last change on this file since 707 was 707, checked in by raasch, 14 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

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