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

Last change on this file since 198 was 198, checked in by raasch, 13 years ago

file headers updated for the next release 3.5

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