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

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

Bugfix: preprocessor statements for parallel execution added

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