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

Last change on this file since 852 was 779, checked in by fricke, 12 years ago

last commit documented

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