source: palm/tags/release-3.8/SOURCE/poismg.f90 @ 716

Last change on this file since 716 was 708, checked in by raasch, 10 years ago

last commit documented

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