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

Last change on this file since 181 was 181, checked in by raasch, 16 years ago

bugfixes + adjustments for SGI ICE system

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