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

Last change on this file since 392 was 392, checked in by raasch, 12 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

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