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

Last change on this file since 368 was 257, checked in by heinze, 16 years ago

Output of messages replaced by message handling routine

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