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

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

preliminary updates for implementing buildings in poismg

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