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

Last change on this file since 668 was 668, checked in by suehring, 13 years ago

last commit documented

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