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

Last change on this file since 881 was 881, checked in by raasch, 12 years ago

last commit documented

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