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

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

parameter use_prior_plot1d_parameters removed; little reformatting

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