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

Last change on this file since 622 was 622, checked in by raasch, 14 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


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