Home

Context Navigation

poismg.f90 @ 402

Last change on this file since 402 was 392, checked in by raasch, 15 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

Rev	Line
[1]	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	!
[257]	8	! Current revisions:
[1]	9	! -----------------
[392]	10	!
[1]	11	!
	12	! Former revisions:
	13	! -----------------
[3]	14	! $Id: poismg.f90 392 2009-09-24 10:39:14Z maronga $
[77]	15	!
[392]	16	! 257 2009-03-11 15:17:42Z heinze
	17	! Output of messages replaced by message handling routine.
	18	!
[198]	19	! 181 2008-07-30 07:07:47Z raasch
	20	! Bugfix: grid_level+1 has to be used in restrict for flags-array
	21	!
[139]	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	!
[77]	28	! 75 2007-03-22 09:54:05Z raasch
	29	! 2nd+3rd argument removed from exchange horiz
	30	!
[3]	31	! RCS Log replace by Id keyword, revision history cleaned up
	32	!
[1]	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
[75]	77	CALL exchange_horiz( d )
[1]	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
[257]	122	message_string = 'no sufficient convergence within 1000 cycles'
	123	CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 )
[1]	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
[114]	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
[1]	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)
[114]	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) ) ) &
[1]	207	+ f1_mg(k,l) * p_mg(k,j,i)
[114]	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 ) )
[1]	211	ENDDO
	212	ENDDO
	213	ENDDO
	214	!$OMP END PARALLEL
	215
	216	!
	217	!-- Horizontal boundary conditions
[75]	218	CALL exchange_horiz( r )
[1]	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	!
[114]	231	!-- Top boundary condition
	232	!-- A Neumann boundary condition for r is implicitly set in routine restrict
[1]	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
[114]	262	REAL :: rkjim, rkjip, rkjmi, rkjmim, rkjmip, rkjpi, rkjpim, rkjpip, &
	263	rkmji, rkmjim, rkmjip, rkmjmi, rkmjmim, rkmjmip, rkmjpi, rkmjpim, &
	264	rkmjpip
	265
[1]	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
[114]	278	!
	279	!-- Choose flag array of the upper level
[181]	280	SELECT CASE ( l+1 )
[114]	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
[1]	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
[114]	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
[1]	349	f_mg(kc,jc,ic) = 1.0 / 64.0 * ( &
	350	8.0 * r(k,j,i) &
[114]	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 ) &
[1]	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	)
[114]	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	! )
[1]	384	ENDDO
	385	ENDDO
	386	ENDDO
	387	!$OMP END PARALLEL
	388
	389	!
	390	!-- Horizontal boundary conditions
[75]	391	CALL exchange_horiz( f_mg )
[1]	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
[114]	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
[1]	416
	417
	418	END 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(2k-1,2j,2*i) = p(k,j,i)
	462	!
	463	!-- Points between two coarse-grid points
	464	temp(2k-1,2j,2i+1) = 0.5 ( p(k,j,i) + p(k,j,i+1) )
	465	temp(2k-1,2j+1,2i) = 0.5 ( p(k,j,i) + p(k,j+1,i) )
	466	temp(2k,2j,2i) = 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(2k-1,2j+1,2i+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(2k,2j,2i+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(2k,2j+1,2i) = 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(2k,2j+1,2i+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
[75]	489	CALL exchange_horiz( temp )
[1]	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
[114]	542	REAL :: wall_left, wall_north, wall_right, wall_south, wall_total, wall_top
	543
[1]	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
[114]	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
[1]	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
[114]	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
[1]	598	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
[114]	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) )
[1]	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) * ( &
[114]	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) )
[1]	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) * ( &
[114]	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) )
[1]	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) * ( &
[114]	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) )
[1]	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) * ( &
[114]	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) )
[1]	715	j = jj+2
	716	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
[114]	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) )
[1]	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) * ( &
[114]	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) )
[1]	754	j = jj+2
	755	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
[114]	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) )
[1]	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) * ( &
[114]	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) )
[1]	793	j = jj+2
	794	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
[114]	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) )
[1]	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) * ( &
[114]	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) )
[1]	832	j = jj+2
	833	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
[114]	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) )
[1]	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
[75]	859	CALL exchange_horiz( p_mg )
[1]	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
[114]	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
[1]	918
[114]	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
[1]	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.

Context Navigation

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

Download in other formats: