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

Last change on this file since 667 was 667, checked in by suehring, 11 years ago

summary:


Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90


Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90


Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90


Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90


When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90


Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90


Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90


Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90


Steering parameter dissipation_control added in inipar.

Makefile


Module advec_ws added.

Modules


Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90


New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90


Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90


Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90


Output of total array size was adapted to nbgp.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90


Calls of exchange_horiz are modified.

diffusion_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90


Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90


Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90


determination of target_id's moved to init_pegrid

init_pt_anomaly.f90


Call of exchange_horiz are modified.

init_rankine.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90


Output of advection scheme.

poismg.f90


Calls of exchange_horiz are modified.

prandtl_fluxes.f90


Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90


+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90


Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90


Allocation with nbgp.

wall_fluxes.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90


Array bounds and nx, ny adapted with nbgp.

sor.f90


bug removed in declaration of ddzw(), nz replaced by nzt+1

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