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

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

last commit documented

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