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

Last change on this file since 623 was 623, checked in by raasch, 11 years ago

last commit documented

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