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

Last change on this file since 778 was 778, checked in by fricke, 10 years ago

Modifications of the multigrid pressure solver

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