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

Last change on this file since 879 was 879, checked in by raasch, 9 years ago

Bugfix: preprocessor statements for parallel execution added

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