Home

Context Navigation

source: palm/trunk/SOURCE/poismg_mod.f90 @ 3182

Last change on this file since 3182 was 3182, checked in by suehring, 7 years ago
New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed
Property svn:keywords set to `Id`
File size: 87.5 KB

Rev	Line
[1931]	1	!> @file poismg.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1575]	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1575]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[2718]	17	! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1575]	19	!
	20	! Current revisions:
	21	! -----------------
[3182]	22	! Rename variables in mesoscale-offline nesting mode
[1610]	23	!
[1576]	24	! Former revisions:
	25	! -----------------
	26	! $Id: poismg_mod.f90 3182 2018-07-27 13:36:03Z suehring $
[2939]	27	! Set lateral boundary conditions for divergence
[3182]	28	!
[2939]	29	! 2937 2018-03-27 14:58:33Z suehring
[2716]	30	! Corrected "Former revisions" section
[2696]	31	!
[2716]	32	! 2696 2017-12-14 17:12:51Z kanani
	33	! Change in file header (GPL part)
	34	! Large-scale forcing implemented (MS)
	35	!
[2696]	36	! 2298 2017-06-29 09:28:18Z raasch
[2298]	37	! sendrecv_in_background related parts removed
	38	!
	39	! 2232 2017-05-30 17:47:52Z suehring
[1576]	40	!
[2074]	41	! 2073 2016-11-30 14:34:05Z raasch
	42	! change of openmp directives in restrict
	43	!
[2038]	44	! 2037 2016-10-26 11:15:40Z knoop
	45	! Anelastic approximation implemented (stll error in optimized multigrid)
	46	!
[2022]	47	! 2021 2016-10-07 14:08:57Z suehring
	48	! Bugfix: restore nest_bound_(l/r/s/n) in case of mg_switch_to_pe0
	49	!
[2001]	50	! 2000 2016-08-20 18:09:15Z knoop
	51	! Forced header and separation lines into 80 columns
	52	!
[1935]	53	! 1934 2016-06-13 09:46:57Z hellstea
	54	! Rename subroutines and cpu-measure log points to indicate default version
	55	!
[1905]	56	! 1904 2016-05-11 13:06:12Z suehring
	57	! Bugfix: enable special_exchange_horiz only for finer grid levels.
	58	! Some formatting adjustments and variable descriptions.
	59	!
[1899]	60	! 1898 2016-05-03 11:27:17Z suehring
	61	! Bugfix: bottom and top boundary condition in resid_fast
	62	! Bugfix: restriction at nzb+1
	63	! formatting adjustments, variable descriptions added in some declaration blocks
	64	!
[1851]	65	! 1850 2016-04-08 13:29:27Z maronga
	66	! Module renamed
	67	!
	68	!
[1763]	69	! 1762 2016-02-25 12:31:13Z hellstea
	70	! Introduction of nested domain feature
	71	!
[1683]	72	! 1682 2015-10-07 23:56:08Z knoop
	73	! Code annotations made doxygen readable
	74	!
[1610]	75	! 1609 2015-07-03 15:37:58Z maronga
	76	! Bugfix: allow compilation without __parallel.
	77	!
[1576]	78	! 1575 2015-03-27 09:56:27Z raasch
[1575]	79	! Initial revision.
	80	! Routine re-written and optimised based on poismg.
	81	!
	82	! Following optimisations have been made:
	83	! - vectorisation (for Intel-CPUs) of the red-black algorithm by resorting
	84	! array elements with even and odd indices
	85	! - explicit boundary conditions for building walls removed (solver is
	86	! running through the buildings
	87	! - reduced data transfer in case of ghost point exchange, because only
	88	! "red" or "black" data points need to be exchanged. This is not applied
	89	! for coarser grid levels, since for then the transfer time is latency bound
	90	!
	91	!
	92	! Description:
	93	! ------------
[1682]	94	!> Solves the Poisson equation for the perturbation pressure with a multigrid
	95	!> V- or W-Cycle scheme.
	96	!>
	97	!> This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
	98	!> September 2000 - July 2001. It has been optimised for speed by Klaus
	99	!> Ketelsen in November 2014.
	100	!>
	101	!> @attention Loop unrolling and cache optimization in SOR-Red/Black method
	102	!> still does not give the expected speedup!
	103	!>
	104	!> @todo Further work required.
[1575]	105	!------------------------------------------------------------------------------!
[1682]	106	MODULE poismg_mod
	107
[2696]	108	USE control_parameters, &
[3182]	109	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
	110	bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
	111	bc_radiation_s, grid_level, nesting_offline
[1575]	112
	113	USE cpulog, &
	114	ONLY: cpu_log, log_point_s
[1609]	115
[1575]	116	USE kinds
[1609]	117
[1575]	118	USE pegrid
	119
	120	PRIVATE
	121
[1682]	122	INTEGER, SAVE :: ind_even_odd !< border index between even and odd k index
	123	INTEGER, DIMENSION(:), SAVE, ALLOCATABLE :: even_odd_level !< stores ind_even_odd for all MG levels
[1575]	124
[1682]	125	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: f1_mg_b, f2_mg_b, f3_mg_b !< blocked version of f1_mg ...
[1575]	126
[1931]	127	INTERFACE poismg
	128	MODULE PROCEDURE poismg
	129	END INTERFACE poismg
[1575]	130
	131	INTERFACE sort_k_to_even_odd_blocks
	132	MODULE PROCEDURE sort_k_to_even_odd_blocks
	133	MODULE PROCEDURE sort_k_to_even_odd_blocks_int
	134	MODULE PROCEDURE sort_k_to_even_odd_blocks_1d
	135	END INTERFACE sort_k_to_even_odd_blocks
	136
[1931]	137	PUBLIC poismg
[1575]	138
	139	CONTAINS
	140
[1682]	141	!------------------------------------------------------------------------------!
	142	! Description:
	143	! ------------
	144	!> Solves the Poisson equation for the perturbation pressure with a multigrid
	145	!> V- or W-Cycle scheme.
	146	!------------------------------------------------------------------------------!
[1931]	147	SUBROUTINE poismg( r )
[1575]	148
	149	USE arrays_3d, &
	150	ONLY: d, p_loc
	151
	152	USE control_parameters, &
[2939]	153	ONLY: bc_lr_cyc, bc_ns_cyc, gathered_size, grid_level, &
	154	grid_level_count, ibc_p_t, &
[1575]	155	maximum_grid_level, message_string, mgcycles, mg_cycles, &
	156	mg_switch_to_pe0_level, residual_limit, subdomain_size
	157
	158	USE cpulog, &
	159	ONLY: cpu_log, log_point_s
	160
	161	USE indices, &
	162	ONLY: nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn,&
	163	nyng, nyn_mg, nzb, nzt, nzt_mg
	164
	165	IMPLICIT NONE
	166
[1682]	167	REAL(wp) :: maxerror !<
	168	REAL(wp) :: maximum_mgcycles !<
	169	REAL(wp) :: residual_norm !<
[1575]	170
[1682]	171	REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !<
[1575]	172
[1682]	173	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !<
[1575]	174
	175
[1931]	176	CALL cpu_log( log_point_s(29), 'poismg', 'start' )
[1575]	177	!
	178	!-- Initialize arrays and variables used in this subroutine
	179
	180	!-- If the number of grid points of the gathered grid, which is collected
	181	!-- on PE0, is larger than the number of grid points of an PE, than array
	182	!-- p3 will be enlarged.
	183	IF ( gathered_size > subdomain_size ) THEN
	184	ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
	185	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
	186	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
	187	mg_switch_to_pe0_level)+1) )
	188	ELSE
	189	ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	190	ENDIF
	191
	192	p3 = 0.0_wp
	193
	194
	195	!
	196	!-- Ghost boundaries have to be added to divergence array.
	197	!-- Exchange routine needs to know the grid level!
	198	grid_level = maximum_grid_level
	199	CALL exchange_horiz( d, 1)
[2939]	200	!
	201	!-- Set bottom and top boundary conditions
[1575]	202	d(nzb,:,:) = d(nzb+1,:,:)
[2939]	203	IF ( ibc_p_t == 1 ) d(nzt+1,:,: ) = d(nzt,:,:)
[1575]	204	!
[2939]	205	!-- Set lateral boundary conditions in non-cyclic case
	206	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	207	IF ( bc_dirichlet_l .OR. bc_radiation_l ) &
[2939]	208	d(:,:,nxl-1) = d(:,:,nxl)
[3182]	209	IF ( bc_dirichlet_r .OR. bc_radiation_r ) &
[2939]	210	d(:,:,nxr+1) = d(:,:,nxr)
	211	ENDIF
	212	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	213	IF ( bc_dirichlet_n .OR. bc_radiation_n ) &
[2939]	214	d(:,nyn+1,:) = d(:,nyn,:)
[3182]	215	IF ( bc_dirichlet_s .OR. bc_radiation_s ) &
[2939]	216	d(:,nys-1,:) = d(:,nys,:)
	217	ENDIF
	218	!
[1575]	219	!-- Initiation of the multigrid scheme. Does n cycles until the
	220	!-- residual is smaller than the given limit. The accuracy of the solution
	221	!-- of the poisson equation will increase with the number of cycles.
	222	!-- If the number of cycles is preset by the user, this number will be
	223	!-- carried out regardless of the accuracy.
	224	grid_level_count = 0
	225	mgcycles = 0
	226	IF ( mg_cycles == -1 ) THEN
	227	maximum_mgcycles = 0
	228	residual_norm = 1.0_wp
	229	ELSE
	230	maximum_mgcycles = mg_cycles
	231	residual_norm = 0.0_wp
	232	ENDIF
	233
	234	!
	235	!-- Initial settings for sorting k-dimension from sequential order (alternate
	236	!-- even/odd) into blocks of even and odd or vice versa
	237	CALL init_even_odd_blocks
	238
	239	!
	240	!-- Sort input arrays in even/odd blocks along k-dimension
	241	CALL sort_k_to_even_odd_blocks( d, grid_level )
	242	CALL sort_k_to_even_odd_blocks( p_loc, grid_level )
	243
	244	!
	245	!-- The complete multigrid cycles are running in block mode, i.e. over
	246	!-- seperate data blocks of even and odd indices
	247	DO WHILE ( residual_norm > residual_limit .OR. &
	248	mgcycles < maximum_mgcycles )
	249
[1931]	250	CALL next_mg_level( d, p_loc, p3, r)
[1575]	251
	252	!
	253	!-- Calculate the residual if the user has not preset the number of
	254	!-- cycles to be performed
	255	IF ( maximum_mgcycles == 0 ) THEN
[1931]	256	CALL resid( d, p_loc, r )
[1575]	257	maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
	258
	259	#if defined( __parallel )
	260	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	261	CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, &
	262	MPI_SUM, comm2d, ierr)
	263	#else
	264	residual_norm = maxerror
	265	#endif
	266	residual_norm = SQRT( residual_norm )
	267	ENDIF
	268
	269	mgcycles = mgcycles + 1
	270
	271	!
	272	!-- If the user has not limited the number of cycles, stop the run in case
	273	!-- of insufficient convergence
	274	IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN
	275	message_string = 'no sufficient convergence within 1000 cycles'
[1931]	276	CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 )
[1575]	277	ENDIF
	278
	279	ENDDO
	280
	281	DEALLOCATE( p3 )
	282	!
	283	!-- Result has to be sorted back from even/odd blocks to sequential order
	284	CALL sort_k_to_sequential( p_loc )
	285	!
	286	!-- Unset the grid level. Variable is used to determine the MPI datatypes for
	287	!-- ghost point exchange
	288	grid_level = 0
	289
[1931]	290	CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
[1575]	291
[1931]	292	END SUBROUTINE poismg
[1575]	293
	294
	295	!------------------------------------------------------------------------------!
	296	! Description:
	297	! ------------
[1682]	298	!> Computes the residual of the perturbation pressure.
[1575]	299	!------------------------------------------------------------------------------!
[1931]	300	SUBROUTINE resid( f_mg, p_mg, r )
[1575]	301
[1682]	302
[1575]	303	USE arrays_3d, &
[2037]	304	ONLY: f1_mg, f2_mg, f3_mg, rho_air_mg
[1575]	305
	306	USE control_parameters, &
[2696]	307	ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
[1575]	308	USE grid_variables, &
	309	ONLY: ddx2_mg, ddy2_mg
	310
	311	USE indices, &
[1898]	312	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]	313
	314	IMPLICIT NONE
	315
[1898]	316	INTEGER(iwp) :: i !< index variable along x
	317	INTEGER(iwp) :: j !< index variable along y
	318	INTEGER(iwp) :: k !< index variable along z
	319	INTEGER(iwp) :: l !< index indicating grid level
	320	INTEGER(iwp) :: km1 !< index variable along z dimension (k-1)
	321	INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1)
[1575]	322
	323	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	324	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	325	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< velocity divergence
[1575]	326	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	327	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	328	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< perturbation pressure
[1575]	329	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	330	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	331	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !< residuum of perturbation pressure
[1575]	332
	333	!
	334	!-- Calculate the residual
	335	l = grid_level
	336
	337	CALL cpu_log( log_point_s(53), 'resid', 'start' )
[1898]	338	!$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
	339	!$OMP DO
	340	DO i = nxl_mg(l), nxr_mg(l)
	341	DO j = nys_mg(l), nyn_mg(l)
[1575]	342	!DIR$ IVDEP
[1898]	343	DO k = ind_even_odd+1, nzt_mg(l)
	344	km1 = k-ind_even_odd-1
	345	kp1 = k-ind_even_odd
	346	r(k,j,i) = f_mg(k,j,i) &
[2232]	347	- rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	348	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	349	- rho_air_mg(k,l) * ddy2_mg(l) * &
[1898]	350	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	351	- f2_mg_b(k,l) * p_mg(kp1,j,i) &
	352	- f3_mg_b(k,l) * p_mg(km1,j,i) &
[1575]	353	+ f1_mg_b(k,l) * p_mg(k,j,i)
[1898]	354	ENDDO
	355	!DIR$ IVDEP
	356	DO k = nzb+1, ind_even_odd
	357	km1 = k+ind_even_odd
	358	kp1 = k+ind_even_odd+1
	359	r(k,j,i) = f_mg(k,j,i) &
[2232]	360	- rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	361	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	362	- rho_air_mg(k,l) * ddy2_mg(l) * &
[1898]	363	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	364	- f2_mg_b(k,l) * p_mg(kp1,j,i) &
	365	- f3_mg_b(k,l) * p_mg(km1,j,i) &
[1575]	366	+ f1_mg_b(k,l) * p_mg(k,j,i)
	367	ENDDO
	368	ENDDO
[1898]	369	ENDDO
	370	!$OMP END PARALLEL
[1575]	371	!
	372	!-- Horizontal boundary conditions
	373	CALL exchange_horiz( r, 1)
	374
	375	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	376	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[1762]	377	r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
	378	ENDIF
[3182]	379	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
[1762]	380	r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
	381	ENDIF
[1575]	382	ENDIF
	383
	384	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	385	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
[1762]	386	r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
	387	ENDIF
[3182]	388	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[1762]	389	r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
	390	ENDIF
[1575]	391	ENDIF
	392
	393	!
[1898]	394	!-- Boundary conditions at bottom and top of the domain. Points may be within
[1575]	395	!-- buildings, but that doesn't matter.
	396	IF ( ibc_p_b == 1 ) THEN
[1898]	397	!
	398	!-- equivalent to r(nzb,:,: ) = r(nzb+1,:,:)
	399	r(nzb,:,: ) = r(ind_even_odd+1,:,:)
[1575]	400	ELSE
	401	r(nzb,:,: ) = 0.0_wp
	402	ENDIF
	403
	404	IF ( ibc_p_t == 1 ) THEN
[1898]	405	!
	406	!-- equivalent to r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
	407	r(nzt_mg(l)+1,:,: ) = r(ind_even_odd,:,:)
[1575]	408	ELSE
	409	r(nzt_mg(l)+1,:,: ) = 0.0_wp
	410	ENDIF
	411
	412	CALL cpu_log( log_point_s(53), 'resid', 'stop' )
	413
[1931]	414	END SUBROUTINE resid
[1575]	415
	416
	417	!------------------------------------------------------------------------------!
	418	! Description:
	419	! ------------
[1682]	420	!> Interpolates the residual on the next coarser grid with "full weighting"
	421	!> scheme
[1575]	422	!------------------------------------------------------------------------------!
[1931]	423	SUBROUTINE restrict( f_mg, r )
[1575]	424
[1682]	425
[1575]	426	USE control_parameters, &
[2696]	427	ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
[1575]	428
	429	USE indices, &
[1898]	430	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]	431
	432	IMPLICIT NONE
	433
[1898]	434	INTEGER(iwp) :: i !< index variable along x on finer grid
	435	INTEGER(iwp) :: ic !< index variable along x on coarser grid
	436	INTEGER(iwp) :: j !< index variable along y on finer grid
	437	INTEGER(iwp) :: jc !< index variable along y on coarser grid
	438	INTEGER(iwp) :: k !< index variable along z on finer grid
	439	INTEGER(iwp) :: kc !< index variable along z on coarser grid
	440	INTEGER(iwp) :: l !< index indicating finer grid level
	441	INTEGER(iwp) :: km1 !< index variable along z dimension (k-1 on finer level)
	442	INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1 on finer level)
[1575]	443
	444
	445	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	446	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	447	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: &
	448	f_mg !< Residual on coarser grid level
[1575]	449
	450	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1, &
	451	nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, &
[1898]	452	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: &
	453	r !< Residual on finer grid level
[1575]	454
	455	!
	456	!-- Interpolate the residual
	457	l = grid_level
	458
	459	CALL cpu_log( log_point_s(54), 'restrict', 'start' )
	460	!
[1898]	461	!-- No wall treatment
	462	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1)
[2073]	463	!$OMP DO SCHEDULE( STATIC )
[1898]	464	DO ic = nxl_mg(l), nxr_mg(l)
	465	i = 2*ic
	466	DO jc = nys_mg(l), nyn_mg(l)
[1575]	467	!
[1898]	468	!-- Calculation for the first point along k
	469	j = 2*jc
	470	!
	471	!-- Calculation for the other points along k
	472	!DIR$ IVDEP
	473	DO k = ind_even_odd+1, nzt_mg(l+1) ! Fine grid at this point
	474	km1 = k-ind_even_odd-1
[1575]	475	kp1 = k-ind_even_odd
[1898]	476	kc = k-ind_even_odd ! Coarse grid index
[1575]	477
[1898]	478	f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( &
	479	8.0_wp * r(k,j,i) &
	480	+ 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + &
	481	r(k,j+1,i) + r(k,j-1,i) ) &
	482	+ 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
	483	r(k,j-1,i+1) + r(k,j+1,i+1) ) &
	484	+ 4.0_wp * r(km1,j,i) &
	485	+ 2.0_wp * ( r(km1,j,i-1) + r(km1,j,i+1) + &
	486	r(km1,j+1,i) + r(km1,j-1,i) ) &
	487	+ ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + &
	488	r(km1,j-1,i+1) + r(km1,j+1,i+1) ) &
	489	+ 4.0_wp * r(kp1,j,i) &
	490	+ 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + &
	491	r(kp1,j+1,i) + r(kp1,j-1,i) ) &
	492	+ ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
	493	r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
	494	)
[1575]	495	ENDDO
	496	ENDDO
[1898]	497	ENDDO
[2073]	498	!$OMP ENDDO
[1898]	499	!$OMP END PARALLEL
[1575]	500
	501	!
[1898]	502	!-- Ghost point exchange
	503	CALL exchange_horiz( f_mg, 1)
[1575]	504	!
	505	!-- Horizontal boundary conditions
	506	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	507	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[1762]	508	f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
	509	ENDIF
[3182]	510	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
[1762]	511	f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
	512	ENDIF
[1575]	513	ENDIF
	514
	515	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	516	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
[1762]	517	f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
	518	ENDIF
[3182]	519	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[1762]	520	f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
	521	ENDIF
[1575]	522	ENDIF
	523
	524	!
	525	!-- Boundary conditions at bottom and top of the domain.
	526	!-- These points are not handled by the above loop. Points may be within
[1898]	527	!-- buildings, but that doesn't matter. Remark: f_mg is ordered sequentielly
	528	!-- after interpolation on coarse grid (is ordered in odd-even blocks further
	529	!-- below).
[1575]	530	IF ( ibc_p_b == 1 ) THEN
	531	f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
	532	ELSE
	533	f_mg(nzb,:,: ) = 0.0_wp
	534	ENDIF
	535
	536	IF ( ibc_p_t == 1 ) THEN
	537	f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
	538	ELSE
	539	f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
	540	ENDIF
	541
	542	CALL cpu_log( log_point_s(54), 'restrict', 'stop' )
	543	!
[1898]	544	!-- Since residual is in sequential order after interpolation, an additional
	545	!-- sorting in odd-even blocks along z dimension is required at this point.
[1575]	546	CALL sort_k_to_even_odd_blocks( f_mg , l)
	547
[1931]	548	END SUBROUTINE restrict
[1575]	549
	550
	551	!------------------------------------------------------------------------------!
	552	! Description:
	553	! ------------
[1682]	554	!> Interpolates the correction of the perturbation pressure
	555	!> to the next finer grid.
[1575]	556	!------------------------------------------------------------------------------!
[1931]	557	SUBROUTINE prolong( p, temp )
[1575]	558
[1682]	559
[1575]	560	USE control_parameters, &
[2696]	561	ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
[1575]	562	USE indices, &
	563	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	564
	565	IMPLICIT NONE
	566
[1898]	567	INTEGER(iwp) :: i !< index variable along x on coarser grid level
	568	INTEGER(iwp) :: j !< index variable along y on coarser grid level
	569	INTEGER(iwp) :: k !< index variable along z on coarser grid level
	570	INTEGER(iwp) :: l !< index indicating finer grid level
	571	INTEGER(iwp) :: kp1 !< index variable along z
	572	INTEGER(iwp) :: ke !< index for prolog even
	573	INTEGER(iwp) :: ko !< index for prolog odd
[1575]	574
	575	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	576	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1898]	577	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: &
	578	p !< perturbation pressure on coarser grid level
[1575]	579
	580	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	581	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	582	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: &
	583	temp !< perturbation pressure on finer grid level
[1575]	584
	585
	586	CALL cpu_log( log_point_s(55), 'prolong', 'start' )
	587
	588	!
	589	!-- First, store elements of the coarser grid on the next finer grid
	590	l = grid_level
	591	ind_even_odd = even_odd_level(grid_level-1)
	592
	593	!$OMP PARALLEL PRIVATE (i,j,k,kp1,ke,ko)
	594	!$OMP DO
	595	DO i = nxl_mg(l-1), nxr_mg(l-1)
	596	DO j = nys_mg(l-1), nyn_mg(l-1)
	597
	598	!DIR$ IVDEP
	599	DO k = ind_even_odd+1, nzt_mg(l-1)
	600	kp1 = k - ind_even_odd
	601	ke = 2 * ( k-ind_even_odd - 1 ) + 1
	602	ko = 2 * k - 1
	603	!
	604	!-- Points of the coarse grid are directly stored on the next finer
	605	!-- grid
	606	temp(ko,2j,2i) = p(k,j,i)
	607	!
	608	!-- Points between two coarse-grid points
	609	temp(ko,2j,2i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
	610	temp(ko,2j+1,2i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
	611	temp(ke,2j,2i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
	612	!
	613	!-- Points in the center of the planes stretched by four points
	614	!-- of the coarse grid cube
	615	temp(ko,2j+1,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
	616	p(k,j+1,i) + p(k,j+1,i+1) )
	617	temp(ke,2j,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
	618	p(kp1,j,i) + p(kp1,j,i+1) )
	619	temp(ke,2j+1,2i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + &
	620	p(kp1,j,i) + p(kp1,j+1,i) )
	621	!
	622	!-- Points in the middle of coarse grid cube
	623	temp(ke,2j+1,2i+1) = 0.125_wp * &
	624	( p(k,j,i) + p(k,j,i+1) + &
	625	p(k,j+1,i) + p(k,j+1,i+1) + &
	626	p(kp1,j,i) + p(kp1,j,i+1) + &
	627	p(kp1,j+1,i) + p(kp1,j+1,i+1) )
[1898]	628
[1575]	629	ENDDO
	630
	631	!DIR$ IVDEP
	632	DO k = nzb+1, ind_even_odd
	633	kp1 = k + ind_even_odd + 1
	634	ke = 2 * k
	635	ko = 2 * ( k + ind_even_odd )
	636	!
	637	!-- Points of the coarse grid are directly stored on the next finer
	638	!-- grid
	639	temp(ko,2j,2i) = p(k,j,i)
	640	!
	641	!-- Points between two coarse-grid points
	642	temp(ko,2j,2i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
	643	temp(ko,2j+1,2i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
	644	temp(ke,2j,2i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
	645	!
	646	!-- Points in the center of the planes stretched by four points
	647	!-- of the coarse grid cube
	648	temp(ko,2j+1,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
	649	p(k,j+1,i) + p(k,j+1,i+1) )
	650	temp(ke,2j,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
	651	p(kp1,j,i) + p(kp1,j,i+1) )
	652	temp(ke,2j+1,2i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + &
	653	p(kp1,j,i) + p(kp1,j+1,i) )
	654	!
	655	!-- Points in the middle of coarse grid cube
	656	temp(ke,2j+1,2i+1) = 0.125_wp * &
	657	( p(k,j,i) + p(k,j,i+1) + &
	658	p(k,j+1,i) + p(k,j+1,i+1) + &
	659	p(kp1,j,i) + p(kp1,j,i+1) + &
	660	p(kp1,j+1,i) + p(kp1,j+1,i+1) )
[1898]	661
[1575]	662	ENDDO
	663
	664	ENDDO
	665	ENDDO
	666	!$OMP END PARALLEL
	667
	668	ind_even_odd = even_odd_level(grid_level)
	669	!
	670	!-- Horizontal boundary conditions
	671	CALL exchange_horiz( temp, 1)
	672
	673	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	674	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[1762]	675	temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
	676	ENDIF
[3182]	677	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
[1762]	678	temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
	679	ENDIF
[1575]	680	ENDIF
	681
	682	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	683	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
[1762]	684	temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
	685	ENDIF
[3182]	686	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[1762]	687	temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
	688	ENDIF
[1575]	689	ENDIF
	690
	691	!
	692	!-- Bottom and top boundary conditions
	693	IF ( ibc_p_b == 1 ) THEN
[1898]	694	!
	695	!-- equivalent to temp(nzb,:,: ) = temp(nzb+1,:,:)
[1575]	696	temp(nzb,:,: ) = temp(ind_even_odd+1,:,:)
	697	ELSE
	698	temp(nzb,:,: ) = 0.0_wp
	699	ENDIF
	700
	701	IF ( ibc_p_t == 1 ) THEN
[1898]	702	!
	703	!-- equivalent to temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
[1575]	704	temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:)
	705	ELSE
	706	temp(nzt_mg(l)+1,:,: ) = 0.0_wp
	707	ENDIF
	708
	709	CALL cpu_log( log_point_s(55), 'prolong', 'stop' )
	710
[1931]	711	END SUBROUTINE prolong
[1575]	712
	713
	714	!------------------------------------------------------------------------------!
	715	! Description:
	716	! ------------
[1682]	717	!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
	718	!> 3D-Red-Black decomposition (GS-RB) is used.
[1575]	719	!------------------------------------------------------------------------------!
[1931]	720	SUBROUTINE redblack( f_mg, p_mg )
[1575]	721
[1682]	722
[1575]	723	USE arrays_3d, &
[2037]	724	ONLY: f1_mg, f2_mg, f3_mg, rho_air_mg
[1575]	725
	726	USE control_parameters, &
[2696]	727	ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t, ngsrb
[1575]	728
	729	USE grid_variables, &
	730	ONLY: ddx2_mg, ddy2_mg
	731
	732	USE indices, &
[1898]	733	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]	734
	735	IMPLICIT NONE
	736
[1898]	737	INTEGER(iwp) :: color !< grid point color, either red or black
	738	INTEGER(iwp) :: i !< index variable along x
	739	INTEGER(iwp) :: ic !< index variable along x
	740	INTEGER(iwp) :: j !< index variable along y
	741	INTEGER(iwp) :: jc !< index variable along y
	742	INTEGER(iwp) :: jj !< index variable along y
	743	INTEGER(iwp) :: k !< index variable along z
	744	INTEGER(iwp) :: l !< grid level
	745	INTEGER(iwp) :: n !< loop variable GauÃ-Seidel iterations
	746	INTEGER(iwp) :: km1 !< index variable (k-1)
	747	INTEGER(iwp) :: kp1 !< index variable (k+1)
[1575]	748
[1898]	749	LOGICAL :: unroll !< flag indicating whether loop unrolling is possible
[1575]	750
	751	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	752	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	753	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: &
	754	f_mg !< residual of perturbation pressure
[1575]	755	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	756	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	757	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: &
	758	p_mg !< perturbation pressure
[1575]	759
	760	l = grid_level
	761
[1898]	762	unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. &
[1575]	763	MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
	764
	765	DO n = 1, ngsrb
	766
	767	DO color = 1, 2
	768
[1898]	769	IF ( .NOT. unroll ) THEN
[1575]	770
[1898]	771	CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' )
[1575]	772	!
[1898]	773	!-- Without unrolling of loops, no cache optimization
	774	!$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
	775	!$OMP DO
	776	DO i = nxl_mg(l), nxr_mg(l), 2
	777	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	778	!DIR$ IVDEP
	779	DO k = ind_even_odd+1, nzt_mg(l)
	780	km1 = k-ind_even_odd-1
	781	kp1 = k-ind_even_odd
	782	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	783	rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	784	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	785	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1898]	786	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
[1575]	787	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	788	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	789	- f_mg(k,j,i) )
[1575]	790	ENDDO
	791	ENDDO
[1898]	792	ENDDO
[1575]	793
[1898]	794	!$OMP DO
	795	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	796	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	797	!DIR$ IVDEP
	798	DO k = ind_even_odd+1, nzt_mg(l)
	799	km1 = k-ind_even_odd-1
	800	kp1 = k-ind_even_odd
	801	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	802	rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	803	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	804	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	805	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	806	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	807	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	808	- f_mg(k,j,i) )
[1575]	809	ENDDO
	810	ENDDO
[1898]	811	ENDDO
[1575]	812
[1898]	813	!$OMP DO
	814	DO i = nxl_mg(l), nxr_mg(l), 2
	815	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	816	!DIR$ IVDEP
	817	DO k = nzb+1, ind_even_odd
	818	km1 = k+ind_even_odd
	819	kp1 = k+ind_even_odd+1
	820	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	821	rho_air_mg(k,l) * ddx2_mg(l) * &
[1575]	822	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	823	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	824	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	825	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	826	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	827	- f_mg(k,j,i) )
[1575]	828	ENDDO
	829	ENDDO
[1898]	830	ENDDO
[1575]	831
[1898]	832	!$OMP DO
	833	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	834	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	835	!DIR$ IVDEP
	836	DO k = nzb+1, ind_even_odd
	837	km1 = k+ind_even_odd
	838	kp1 = k+ind_even_odd+1
	839	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	840	rho_air_mg(k,l) * ddx2_mg(l) * &
[1575]	841	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	842	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	843	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	844	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	845	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	846	- f_mg(k,j,i) )
[1575]	847	ENDDO
	848	ENDDO
[1898]	849	ENDDO
	850	!$OMP END PARALLEL
[1575]	851
[1898]	852	CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' )
[1575]	853
[1898]	854	ELSE
[1575]	855	!
[1898]	856	!-- Loop unrolling along y, only one i loop for better cache use
	857	CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' )
[1575]	858
[1898]	859	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
	860	!$OMP DO
	861	DO ic = nxl_mg(l), nxr_mg(l), 2
	862	DO jc = nys_mg(l), nyn_mg(l), 4
	863	i = ic
	864	jj = jc+2-color
	865	!DIR$ IVDEP
	866	DO k = ind_even_odd+1, nzt_mg(l)
	867	km1 = k-ind_even_odd-1
	868	kp1 = k-ind_even_odd
	869	j = jj
	870	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	871	rho_air_mg(k,l) * ddx2_mg(l) * &
[1575]	872	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	873	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	874	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	875	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	876	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	877	- f_mg(k,j,i) )
	878	j = jj+2
	879	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	880	rho_air_mg(k,l) * ddx2_mg(l) * &
[1575]	881	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	882	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	883	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	884	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	885	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	886	- f_mg(k,j,i) )
	887	ENDDO
[1575]	888
[1898]	889	i = ic+1
	890	jj = jc+color-1
	891	!DIR$ IVDEP
	892	DO k = ind_even_odd+1, nzt_mg(l)
	893	km1 = k-ind_even_odd-1
	894	kp1 = k-ind_even_odd
	895	j = jj
	896	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	897	rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	898	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	899	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	900	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	901	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	902	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	903	- f_mg(k,j,i) )
	904	j = jj+2
	905	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	906	rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	907	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	908	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1898]	909	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
[1575]	910	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	911	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	912	- f_mg(k,j,i) )
	913	ENDDO
[1575]	914
[1898]	915	i = ic
	916	jj = jc+color-1
	917	!DIR$ IVDEP
	918	DO k = nzb+1, ind_even_odd
	919	km1 = k+ind_even_odd
	920	kp1 = k+ind_even_odd+1
	921	j = jj
	922	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	923	rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	924	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	925	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1898]	926	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
[1575]	927	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	928	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	929	- f_mg(k,j,i) )
	930	j = jj+2
	931	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	932	rho_air_mg(k,l) * ddx2_mg(l) * &
[1575]	933	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	934	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	935	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	936	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	937	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	938	- f_mg(k,j,i) )
	939	ENDDO
[1575]	940
[1898]	941	i = ic+1
	942	jj = jc+2-color
	943	!DIR$ IVDEP
	944	DO k = nzb+1, ind_even_odd
	945	km1 = k+ind_even_odd
	946	kp1 = k+ind_even_odd+1
	947	j = jj
	948	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	949	rho_air_mg(k,l) * ddx2_mg(l) * &
[1898]	950	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	951	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1575]	952	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	953	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	954	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	955	- f_mg(k,j,i) )
	956	j = jj+2
	957	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
[2232]	958	rho_air_mg(k,l) * ddx2_mg(l) * &
[1575]	959	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
[2232]	960	+ rho_air_mg(k,l) * ddy2_mg(l) * &
[1898]	961	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
[1575]	962	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
	963	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
[1898]	964	- f_mg(k,j,i) )
[1575]	965	ENDDO
	966
	967	ENDDO
[1898]	968	ENDDO
	969	!$OMP END PARALLEL
[1575]	970
[1898]	971	CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' )
[1575]	972
	973	ENDIF
	974
	975	!
	976	!-- Horizontal boundary conditions
	977	CALL special_exchange_horiz( p_mg, color )
	978
[1762]	979	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	980	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[1575]	981	p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
	982	ENDIF
[3182]	983	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
[1575]	984	p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
	985	ENDIF
	986	ENDIF
	987
	988	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	989	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
[1575]	990	p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
	991	ENDIF
[3182]	992	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[1575]	993	p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
	994	ENDIF
	995	ENDIF
	996
	997	!
	998	!-- Bottom and top boundary conditions
	999	IF ( ibc_p_b == 1 ) THEN
[1898]	1000	!
	1001	!-- equivalent to p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
[1575]	1002	p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:)
	1003	ELSE
	1004	p_mg(nzb,:,: ) = 0.0_wp
	1005	ENDIF
	1006
	1007	IF ( ibc_p_t == 1 ) THEN
[1898]	1008	!
	1009	!-- equivalent to p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
[1575]	1010	p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:)
	1011	ELSE
	1012	p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
	1013	ENDIF
	1014
	1015	ENDDO
[1898]	1016
[1575]	1017	ENDDO
	1018
[1931]	1019	END SUBROUTINE redblack
[1575]	1020
	1021
	1022	!------------------------------------------------------------------------------!
	1023	! Description:
	1024	! ------------
[1682]	1025	!> Sort k-Dimension from sequential into blocks of even and odd.
	1026	!> This is required to vectorize the red-black subroutine.
	1027	!> Version for 3D-REAL arrays
[1575]	1028	!------------------------------------------------------------------------------!
[1682]	1029	SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel )
[1575]	1030
[1682]	1031
[1575]	1032	USE control_parameters, &
	1033	ONLY: grid_level
	1034
	1035	USE indices, &
	1036	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	1037
	1038	IMPLICIT NONE
	1039
[1898]	1040	INTEGER(iwp), INTENT(IN) :: glevel !< grid level
[1575]	1041
	1042	REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1, &
	1043	nys_mg(glevel)-1:nyn_mg(glevel)+1, &
[1898]	1044	nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: &
	1045	p_mg !< array to be sorted
[1575]	1046	!
	1047	!-- Local variables
[1898]	1048	INTEGER(iwp) :: i !< index variable along x
	1049	INTEGER(iwp) :: j !< index variable along y
	1050	INTEGER(iwp) :: k !< index variable along z
	1051	INTEGER(iwp) :: l !< grid level
	1052	INTEGER(iwp) :: ind !< index variable along z
	1053	REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< odd-even sorted temporary array
[1575]	1054
	1055
	1056	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
	1057
	1058	l = glevel
	1059	ind_even_odd = even_odd_level(l)
	1060
	1061	!$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
	1062	!$OMP DO
	1063	DO i = nxl_mg(l)-1, nxr_mg(l)+1
	1064	DO j = nys_mg(l)-1, nyn_mg(l)+1
	1065
	1066	!
	1067	!-- Sort the data with even k index
	1068	ind = nzb-1
	1069	DO k = nzb, nzt_mg(l), 2
	1070	ind = ind + 1
	1071	tmp(ind) = p_mg(k,j,i)
	1072	ENDDO
	1073	!
	1074	!-- Sort the data with odd k index
	1075	DO k = nzb+1, nzt_mg(l)+1, 2
	1076	ind = ind + 1
	1077	tmp(ind) = p_mg(k,j,i)
	1078	ENDDO
	1079
	1080	p_mg(:,j,i) = tmp
	1081
	1082	ENDDO
	1083	ENDDO
	1084	!$OMP END PARALLEL
	1085
	1086	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
	1087
	1088	END SUBROUTINE sort_k_to_even_odd_blocks
	1089
	1090
	1091	!------------------------------------------------------------------------------!
	1092	! Description:
	1093	! ------------
[1682]	1094	!> Sort k-Dimension from sequential into blocks of even and odd.
	1095	!> This is required to vectorize the red-black subroutine.
	1096	!> Version for 1D-REAL arrays
[1575]	1097	!------------------------------------------------------------------------------!
[1682]	1098	SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel )
[1575]	1099
[1682]	1100
[1575]	1101	USE indices, &
	1102	ONLY: nzb, nzt_mg
	1103
	1104	IMPLICIT NONE
	1105
[1898]	1106	INTEGER(iwp), INTENT(IN) :: glevel !< grid level
[1575]	1107
[1898]	1108	REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) :: f_mg !< 1D input array
	1109	REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: f_mg_b !< 1D output array
[1575]	1110
	1111	!
	1112	!-- Local variables
[1898]	1113	INTEGER(iwp) :: ind !< index variable along z
	1114	INTEGER(iwp) :: k !< index variable along z
[1575]	1115
	1116
	1117	ind = nzb - 1
	1118	!
	1119	!-- Sort the data with even k index
	1120	DO k = nzb, nzt_mg(glevel), 2
	1121	ind = ind + 1
	1122	IF ( k >= nzb+1 .AND. k <= nzt_mg(glevel) ) THEN
	1123	f_mg_b(ind) = f_mg(k)
	1124	ENDIF
	1125	ENDDO
	1126	!
	1127	!-- Sort the data with odd k index
	1128	DO k = nzb+1, nzt_mg(glevel)+1, 2
	1129	ind = ind + 1
	1130	IF( k >= nzb+1 .AND. k <= nzt_mg(glevel) ) THEN
	1131	f_mg_b(ind) = f_mg(k)
	1132	ENDIF
	1133	ENDDO
	1134
	1135	END SUBROUTINE sort_k_to_even_odd_blocks_1d
	1136
	1137
	1138	!------------------------------------------------------------------------------!
	1139	! Description:
	1140	! ------------
[1682]	1141	!> Sort k-Dimension from sequential into blocks of even and odd.
	1142	!> This is required to vectorize the red-black subroutine.
	1143	!> Version for 2D-INTEGER arrays
[1575]	1144	!------------------------------------------------------------------------------!
[1682]	1145	SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
[1575]	1146
[1682]	1147
[1575]	1148	USE control_parameters, &
	1149	ONLY: grid_level
	1150
	1151	USE indices, &
	1152	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	1153
	1154	IMPLICIT NONE
	1155
[1898]	1156	INTEGER(iwp), INTENT(IN) :: glevel !< grid level
[1575]	1157
	1158	INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1, &
	1159	nys_mg(glevel)-1:nyn_mg(glevel)+1, &
[1898]	1160	nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: &
	1161	i_mg !< array to be sorted
[1575]	1162	!
	1163	!-- Local variables
[1898]	1164	INTEGER(iwp) :: i !< index variabel along x
	1165	INTEGER(iwp) :: j !< index variable along y
	1166	INTEGER(iwp) :: k !< index variable along z
	1167	INTEGER(iwp) :: l !< grid level
	1168	INTEGER(iwp) :: ind !< index variable along z
	1169	INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< temporary odd-even sorted array
[1575]	1170
	1171
	1172	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
	1173
	1174	l = glevel
	1175	ind_even_odd = even_odd_level(l)
	1176
	1177	DO i = nxl_mg(l)-1, nxr_mg(l)+1
	1178	DO j = nys_mg(l)-1, nyn_mg(l)+1
	1179
	1180	!
	1181	!-- Sort the data with even k index
	1182	ind = nzb-1
	1183	DO k = nzb, nzt_mg(l), 2
	1184	ind = ind + 1
	1185	tmp(ind) = i_mg(k,j,i)
	1186	ENDDO
	1187	!
	1188	!++ ATTENTION: Check reason for this error. Remove it or replace WRITE
	1189	!++ by PALM message
[1609]	1190	#if defined ( __parallel )
[1575]	1191	IF ( ind /= ind_even_odd ) THEN
	1192	WRITE (0,*) 'ERROR ==> illegal ind_even_odd ',ind,ind_even_odd,l
	1193	CALL MPI_ABORT(MPI_COMM_WORLD,i,j)
	1194	ENDIF
[1609]	1195	#endif
[1575]	1196	!
	1197	!-- Sort the data with odd k index
	1198	DO k = nzb+1, nzt_mg(l)+1, 2
	1199	ind = ind + 1
	1200	tmp(ind) = i_mg(k,j,i)
	1201	ENDDO
	1202
	1203	i_mg(:,j,i) = tmp
	1204
	1205	ENDDO
	1206	ENDDO
	1207
	1208	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
	1209
	1210	END SUBROUTINE sort_k_to_even_odd_blocks_int
	1211
	1212
	1213	!------------------------------------------------------------------------------!
	1214	! Description:
	1215	! ------------
[1682]	1216	!> Sort k-dimension from blocks of even and odd into sequential
[1575]	1217	!------------------------------------------------------------------------------!
[1682]	1218	SUBROUTINE sort_k_to_sequential( p_mg )
[1575]	1219
[1682]	1220
[1575]	1221	USE control_parameters, &
	1222	ONLY: grid_level
	1223
	1224	USE indices, &
	1225	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	1226
	1227	IMPLICIT NONE
	1228
	1229	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1230	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1898]	1231	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: &
	1232	p_mg !< array to be sorted
[1575]	1233	!
	1234	!-- Local variables
[1898]	1235	INTEGER(iwp) :: i !< index variable along x
	1236	INTEGER(iwp) :: j !< index variable along y
	1237	INTEGER(iwp) :: k !< index variable along z
	1238	INTEGER(iwp) :: l !< grid level
	1239	INTEGER(iwp) :: ind !< index variable along z
[1575]	1240
	1241	REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) :: tmp
	1242
	1243
	1244	l = grid_level
	1245
	1246	!$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
	1247	!$OMP DO
	1248	DO i = nxl_mg(l)-1, nxr_mg(l)+1
	1249	DO j = nys_mg(l)-1, nyn_mg(l)+1
	1250
	1251	ind = nzb - 1
	1252	tmp = p_mg(:,j,i)
	1253	DO k = nzb, nzt_mg(l), 2
	1254	ind = ind + 1
	1255	p_mg(k,j,i) = tmp(ind)
	1256	ENDDO
	1257
	1258	DO k = nzb+1, nzt_mg(l)+1, 2
	1259	ind = ind + 1
	1260	p_mg(k,j,i) = tmp(ind)
	1261	ENDDO
	1262	ENDDO
	1263	ENDDO
	1264	!$OMP END PARALLEL
	1265
	1266	END SUBROUTINE sort_k_to_sequential
	1267
	1268
[1682]	1269	!------------------------------------------------------------------------------!
	1270	! Description:
	1271	! ------------
	1272	!> Gather subdomain data from all PEs.
	1273	!------------------------------------------------------------------------------!
[1931]	1274	SUBROUTINE mg_gather( f2, f2_sub )
[1575]	1275
	1276	USE control_parameters, &
	1277	ONLY: grid_level
	1278
	1279	USE cpulog, &
	1280	ONLY: cpu_log, log_point_s
	1281
	1282	USE indices, &
	1283	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	1284
	1285	IMPLICIT NONE
	1286
[1682]	1287	INTEGER(iwp) :: i !<
	1288	INTEGER(iwp) :: il !<
	1289	INTEGER(iwp) :: ir !<
	1290	INTEGER(iwp) :: j !<
	1291	INTEGER(iwp) :: jn !<
	1292	INTEGER(iwp) :: js !<
	1293	INTEGER(iwp) :: k !<
	1294	INTEGER(iwp) :: nwords !<
[1575]	1295
	1296	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1297	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1298	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !<
[1575]	1299	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1300	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1301	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !<
[1575]	1302
	1303	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
	1304	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
[1682]	1305	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !<
[1575]	1306
	1307
	1308	#if defined( __parallel )
	1309	CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
	1310
	1311	f2_l = 0.0_wp
	1312
	1313	!
	1314	!-- Store the local subdomain array on the total array
	1315	js = mg_loc_ind(3,myid)
	1316	IF ( south_border_pe ) js = js - 1
	1317	jn = mg_loc_ind(4,myid)
	1318	IF ( north_border_pe ) jn = jn + 1
	1319	il = mg_loc_ind(1,myid)
	1320	IF ( left_border_pe ) il = il - 1
	1321	ir = mg_loc_ind(2,myid)
	1322	IF ( right_border_pe ) ir = ir + 1
	1323	DO i = il, ir
	1324	DO j = js, jn
	1325	DO k = nzb, nzt_mg(grid_level)+1
	1326	f2_l(k,j,i) = f2_sub(k,j,i)
	1327	ENDDO
	1328	ENDDO
	1329	ENDDO
	1330
	1331	!
	1332	!-- Find out the number of array elements of the total array
	1333	nwords = SIZE( f2 )
	1334
	1335	!
	1336	!-- Gather subdomain data from all PEs
	1337	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1338	CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
	1339	f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
	1340	nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
	1341
	1342	CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
	1343	#endif
	1344
[1931]	1345	END SUBROUTINE mg_gather
[1575]	1346
	1347
	1348
[1682]	1349	!------------------------------------------------------------------------------!
	1350	! Description:
	1351	! ------------
	1352	!> @todo It might be possible to improve the speed of this routine by using
	1353	!> non-blocking communication
	1354	!------------------------------------------------------------------------------!
[1931]	1355	SUBROUTINE mg_scatter( p2, p2_sub )
[1575]	1356
	1357	USE control_parameters, &
	1358	ONLY: grid_level
	1359
	1360	USE cpulog, &
	1361	ONLY: cpu_log, log_point_s
	1362
	1363	USE indices, &
	1364	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	1365
	1366	IMPLICIT NONE
	1367
[1682]	1368	INTEGER(iwp) :: nwords !<
[1575]	1369
	1370	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	1371	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	1372	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[1575]	1373
	1374	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
	1375	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
[1682]	1376	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !<
[1575]	1377
	1378	!
	1379	!-- Find out the number of array elements of the subdomain array
	1380	nwords = SIZE( p2_sub )
	1381
	1382	#if defined( __parallel )
	1383	CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
	1384
	1385	p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	1386	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
	1387
	1388	CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
	1389	#endif
	1390
[1931]	1391	END SUBROUTINE mg_scatter
[1575]	1392
	1393
	1394	!------------------------------------------------------------------------------!
	1395	! Description:
	1396	! ------------
[1682]	1397	!> This is where the multigrid technique takes place. V- and W- Cycle are
	1398	!> implemented and steered by the parameter "gamma". Parameter "nue" determines
	1399	!> the convergence of the multigrid iterative solution. There are nue times
	1400	!> RB-GS iterations. It should be set to "1" or "2", considering the time effort
	1401	!> one would like to invest. Last choice shows a very good converging factor,
	1402	!> but leads to an increase in computing time.
[1575]	1403	!------------------------------------------------------------------------------!
[1931]	1404	RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
[1575]	1405
	1406	USE control_parameters, &
	1407	ONLY: bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir, &
[3182]	1408	child_domain, gamma_mg, grid_level_count, ibc_p_b, ibc_p_t, &
[2696]	1409	maximum_grid_level, mg_switch_to_pe0_level, &
[3182]	1410	mg_switch_to_pe0, ngsrb
[1575]	1411
	1412	USE indices, &
	1413	ONLY: mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn, &
	1414	nyn_mg, nzb, nzt, nzt_mg
	1415
	1416	IMPLICIT NONE
	1417
[1898]	1418	INTEGER(iwp) :: i !< index variable along x
	1419	INTEGER(iwp) :: j !< index variable along y
	1420	INTEGER(iwp) :: k !< index variable along z
[1682]	1421	INTEGER(iwp) :: nxl_mg_save !<
	1422	INTEGER(iwp) :: nxr_mg_save !<
	1423	INTEGER(iwp) :: nyn_mg_save !<
	1424	INTEGER(iwp) :: nys_mg_save !<
	1425	INTEGER(iwp) :: nzt_mg_save !<
[1575]	1426
	1427	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1428	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1429	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1575]	1430	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1431	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1432	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1575]	1433	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1434	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1435	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !<
[1575]	1436	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1437	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1438	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
[1575]	1439
	1440	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	1441	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	1442	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !<
[1575]	1443	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	1444	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	1445	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[1575]	1446
[1682]	1447	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !<
	1448	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p2_sub !<
[1575]	1449
	1450	!
	1451	!-- Restriction to the coarsest grid
	1452	10 IF ( grid_level == 1 ) THEN
	1453
	1454	!
	1455	!-- Solution on the coarsest grid. Double the number of Gauss-Seidel
	1456	!-- iterations in order to get a more accurate solution.
	1457	ngsrb = 2 * ngsrb
	1458
	1459	ind_even_odd = even_odd_level(grid_level)
	1460
[1931]	1461	CALL redblack( f_mg, p_mg )
[1575]	1462
	1463	ngsrb = ngsrb / 2
	1464
	1465
	1466	ELSEIF ( grid_level /= 1 ) THEN
	1467
	1468	grid_level_count(grid_level) = grid_level_count(grid_level) + 1
	1469
	1470	!
	1471	!-- Solution on the actual grid level
	1472	ind_even_odd = even_odd_level(grid_level)
	1473
[1931]	1474	CALL redblack( f_mg, p_mg )
[1575]	1475
	1476	!
	1477	!-- Determination of the actual residual
[1931]	1478	CALL resid( f_mg, p_mg, r )
[1575]	1479
	1480	!-- Restriction of the residual (finer grid values!) to the next coarser
	1481	!-- grid. Therefore, the grid level has to be decremented now. nxl..nzt have
	1482	!-- to be set to the coarse grid values, because these variables are needed
	1483	!-- for the exchange of ghost points in routine exchange_horiz
	1484	grid_level = grid_level - 1
	1485
	1486	nxl = nxl_mg(grid_level)
	1487	nys = nys_mg(grid_level)
	1488	nxr = nxr_mg(grid_level)
	1489	nyn = nyn_mg(grid_level)
	1490	nzt = nzt_mg(grid_level)
	1491
	1492	IF ( grid_level == mg_switch_to_pe0_level ) THEN
	1493
	1494	!
	1495	!-- From this level on, calculations are done on PE0 only.
	1496	!-- First, carry out restriction on the subdomain.
	1497	!-- Therefore, indices of the level have to be changed to subdomain
	1498	!-- values in between (otherwise, the restrict routine would expect
	1499	!-- the gathered array)
	1500
	1501	nxl_mg_save = nxl_mg(grid_level)
	1502	nxr_mg_save = nxr_mg(grid_level)
	1503	nys_mg_save = nys_mg(grid_level)
	1504	nyn_mg_save = nyn_mg(grid_level)
	1505	nzt_mg_save = nzt_mg(grid_level)
	1506	nxl_mg(grid_level) = mg_loc_ind(1,myid)
	1507	nxr_mg(grid_level) = mg_loc_ind(2,myid)
	1508	nys_mg(grid_level) = mg_loc_ind(3,myid)
	1509	nyn_mg(grid_level) = mg_loc_ind(4,myid)
	1510	nzt_mg(grid_level) = mg_loc_ind(5,myid)
	1511	nxl = mg_loc_ind(1,myid)
	1512	nxr = mg_loc_ind(2,myid)
	1513	nys = mg_loc_ind(3,myid)
	1514	nyn = mg_loc_ind(4,myid)
	1515	nzt = mg_loc_ind(5,myid)
	1516
	1517	ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1, &
	1518	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1519	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
	1520
[1931]	1521	CALL restrict( f2_sub, r )
[1575]	1522
	1523	!
	1524	!-- Restore the correct indices of this level
	1525	nxl_mg(grid_level) = nxl_mg_save
	1526	nxr_mg(grid_level) = nxr_mg_save
	1527	nys_mg(grid_level) = nys_mg_save
	1528	nyn_mg(grid_level) = nyn_mg_save
	1529	nzt_mg(grid_level) = nzt_mg_save
	1530	nxl = nxl_mg(grid_level)
	1531	nxr = nxr_mg(grid_level)
	1532	nys = nys_mg(grid_level)
	1533	nyn = nyn_mg(grid_level)
	1534	nzt = nzt_mg(grid_level)
	1535	!
	1536	!-- Gather all arrays from the subdomains on PE0
[1931]	1537	CALL mg_gather( f2, f2_sub )
[1575]	1538
	1539	!
	1540	!-- Set switch for routine exchange_horiz, that no ghostpoint exchange
	1541	!-- has to be carried out from now on
	1542	mg_switch_to_pe0 = .TRUE.
	1543
	1544	!
	1545	!-- In case of non-cyclic lateral boundary conditions, both in- and
	1546	!-- outflow conditions have to be used on all PEs after the switch,
	1547	!-- because then they have the total domain.
	1548	IF ( bc_lr_dirrad ) THEN
[3182]	1549	bc_dirichlet_l = .TRUE.
	1550	bc_dirichlet_r = .FALSE.
	1551	bc_radiation_l = .FALSE.
	1552	bc_radiation_r = .TRUE.
[1575]	1553	ELSEIF ( bc_lr_raddir ) THEN
[3182]	1554	bc_dirichlet_l = .FALSE.
	1555	bc_dirichlet_r = .TRUE.
	1556	bc_radiation_l = .TRUE.
	1557	bc_radiation_r = .FALSE.
	1558	ELSEIF ( child_domain .OR. nesting_offline ) THEN
	1559	bc_dirichlet_l = .TRUE.
	1560	bc_dirichlet_r = .TRUE.
[1575]	1561	ENDIF
	1562
	1563	IF ( bc_ns_dirrad ) THEN
[3182]	1564	bc_dirichlet_n = .TRUE.
	1565	bc_dirichlet_s = .FALSE.
	1566	bc_radiation_n = .FALSE.
	1567	bc_radiation_s = .TRUE.
[1575]	1568	ELSEIF ( bc_ns_raddir ) THEN
[3182]	1569	bc_dirichlet_n = .FALSE.
	1570	bc_dirichlet_s = .TRUE.
	1571	bc_radiation_n = .TRUE.
	1572	bc_radiation_s = .FALSE.
	1573	ELSEIF ( child_domain .OR. nesting_offline) THEN
	1574	bc_dirichlet_s = .TRUE.
	1575	bc_dirichlet_n = .TRUE.
[1575]	1576	ENDIF
	1577
	1578	DEALLOCATE( f2_sub )
	1579
	1580	ELSE
	1581
[1931]	1582	CALL restrict( f2, r )
[1575]	1583
	1584	ind_even_odd = even_odd_level(grid_level) ! must be after restrict
	1585
	1586	ENDIF
	1587
	1588	p2 = 0.0_wp
	1589
	1590	!
	1591	!-- Repeat the same procedure till the coarsest grid is reached
[1931]	1592	CALL next_mg_level( f2, p2, p3, r )
[1575]	1593
	1594	ENDIF
	1595
	1596	!
	1597	!-- Now follows the prolongation
	1598	IF ( grid_level >= 2 ) THEN
	1599
	1600	!
	1601	!-- Prolongation of the new residual. The values are transferred
	1602	!-- from the coarse to the next finer grid.
	1603	IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN
	1604
	1605	#if defined( __parallel )
	1606	!
	1607	!-- At this level, the new residual first has to be scattered from
	1608	!-- PE0 to the other PEs
	1609	ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1, &
	1610	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	1611	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
	1612
[1931]	1613	CALL mg_scatter( p2, p2_sub )
[1575]	1614
	1615	!
	1616	!-- Therefore, indices of the previous level have to be changed to
	1617	!-- subdomain values in between (otherwise, the prolong routine would
	1618	!-- expect the gathered array)
	1619	nxl_mg_save = nxl_mg(grid_level-1)
	1620	nxr_mg_save = nxr_mg(grid_level-1)
	1621	nys_mg_save = nys_mg(grid_level-1)
	1622	nyn_mg_save = nyn_mg(grid_level-1)
	1623	nzt_mg_save = nzt_mg(grid_level-1)
	1624	nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
	1625	nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
	1626	nys_mg(grid_level-1) = mg_loc_ind(3,myid)
	1627	nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
	1628	nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
	1629
	1630	!
	1631	!-- Set switch for routine exchange_horiz, that ghostpoint exchange
	1632	!-- has to be carried again out from now on
	1633	mg_switch_to_pe0 = .FALSE.
	1634
	1635	!
[2021]	1636	!-- For non-cyclic lateral boundary conditions and in case of nesting,
	1637	!-- restore the in-/outflow conditions.
[3182]	1638	bc_dirichlet_l = .FALSE.; bc_dirichlet_r = .FALSE.
	1639	bc_dirichlet_n = .FALSE.; bc_dirichlet_s = .FALSE.
	1640	bc_radiation_l = .FALSE.; bc_radiation_r = .FALSE.
	1641	bc_radiation_n = .FALSE.; bc_radiation_s = .FALSE.
[1575]	1642
	1643	IF ( pleft == MPI_PROC_NULL ) THEN
[3182]	1644	IF ( bc_lr_dirrad .OR. child_domain .OR. nesting_offline ) &
	1645	THEN
	1646	bc_dirichlet_l = .TRUE.
[1575]	1647	ELSEIF ( bc_lr_raddir ) THEN
[3182]	1648	bc_radiation_l = .TRUE.
[1575]	1649	ENDIF
	1650	ENDIF
	1651
	1652	IF ( pright == MPI_PROC_NULL ) THEN
	1653	IF ( bc_lr_dirrad ) THEN
[3182]	1654	bc_radiation_r = .TRUE.
	1655	ELSEIF ( bc_lr_raddir .OR. child_domain .OR. &
	1656	nesting_offline ) THEN
	1657	bc_dirichlet_r = .TRUE.
[1575]	1658	ENDIF
	1659	ENDIF
	1660
	1661	IF ( psouth == MPI_PROC_NULL ) THEN
	1662	IF ( bc_ns_dirrad ) THEN
[3182]	1663	bc_radiation_s = .TRUE.
	1664	ELSEIF ( bc_ns_raddir .OR. child_domain .OR. &
	1665	nesting_offline ) THEN
	1666	bc_dirichlet_s = .TRUE.
[1575]	1667	ENDIF
	1668	ENDIF
	1669
	1670	IF ( pnorth == MPI_PROC_NULL ) THEN
[3182]	1671	IF ( bc_ns_dirrad .OR. child_domain .OR. nesting_offline ) &
	1672	THEN
	1673	bc_dirichlet_n = .TRUE.
[1575]	1674	ELSEIF ( bc_ns_raddir ) THEN
[3182]	1675	bc_radiation_n = .TRUE.
[1575]	1676	ENDIF
	1677	ENDIF
	1678
[1931]	1679	CALL prolong( p2_sub, p3 )
[1575]	1680
	1681	!
	1682	!-- Restore the correct indices of the previous level
	1683	nxl_mg(grid_level-1) = nxl_mg_save
	1684	nxr_mg(grid_level-1) = nxr_mg_save
	1685	nys_mg(grid_level-1) = nys_mg_save
	1686	nyn_mg(grid_level-1) = nyn_mg_save
	1687	nzt_mg(grid_level-1) = nzt_mg_save
	1688
	1689	DEALLOCATE( p2_sub )
	1690	#endif
	1691
	1692	ELSE
	1693
[1931]	1694	CALL prolong( p2, p3 )
[1575]	1695
	1696	ENDIF
	1697
	1698	!
	1699	!-- Computation of the new pressure correction. Therefore,
	1700	!-- values from prior grids are added up automatically stage by stage.
	1701	DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
	1702	DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
	1703	DO k = nzb, nzt_mg(grid_level)+1
	1704	p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
	1705	ENDDO
	1706	ENDDO
	1707	ENDDO
	1708
	1709	!
	1710	!-- Relaxation of the new solution
[1931]	1711	CALL redblack( f_mg, p_mg )
[1575]	1712
	1713	ENDIF
	1714
	1715
	1716	!
	1717	!-- The following few lines serve the steering of the multigrid scheme
	1718	IF ( grid_level == maximum_grid_level ) THEN
	1719
	1720	GOTO 20
	1721
	1722	ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. &
	1723	grid_level_count(grid_level) /= gamma_mg ) THEN
	1724
	1725	GOTO 10
	1726
	1727	ENDIF
	1728
	1729	!
[1931]	1730	!-- Reset counter for the next call of poismg
[1575]	1731	grid_level_count(grid_level) = 0
	1732
	1733	!
	1734	!-- Continue with the next finer level. nxl..nzt have to be
	1735	!-- set to the finer grid values, because these variables are needed for the
	1736	!-- exchange of ghost points in routine exchange_horiz
	1737	grid_level = grid_level + 1
	1738	ind_even_odd = even_odd_level(grid_level)
	1739
	1740	nxl = nxl_mg(grid_level)
	1741	nxr = nxr_mg(grid_level)
	1742	nys = nys_mg(grid_level)
	1743	nyn = nyn_mg(grid_level)
	1744	nzt = nzt_mg(grid_level)
	1745
	1746	20 CONTINUE
	1747
[1931]	1748	END SUBROUTINE next_mg_level
[1575]	1749
	1750
	1751	!------------------------------------------------------------------------------!
	1752	! Description:
	1753	! ------------
[1682]	1754	!> Initial settings for sorting k-dimension from sequential order (alternate
	1755	!> even/odd) into blocks of even and odd or vice versa
[1575]	1756	!------------------------------------------------------------------------------!
[1682]	1757	SUBROUTINE init_even_odd_blocks
[1575]	1758
[1682]	1759
[1575]	1760	USE arrays_3d, &
[2232]	1761	ONLY: f1_mg, f2_mg, f3_mg
[1575]	1762
	1763	USE control_parameters, &
[1898]	1764	ONLY: grid_level, maximum_grid_level
[1575]	1765
	1766	USE indices, &
	1767	ONLY: nzb, nzt, nzt_mg
	1768
	1769	USE indices, &
[1898]	1770	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]	1771
	1772	IMPLICIT NONE
	1773	!
	1774	!-- Local variables
[1898]	1775	INTEGER(iwp) :: i !<
[1682]	1776	INTEGER(iwp) :: l !<
[1575]	1777
	1778	LOGICAL, SAVE :: lfirst = .TRUE.
	1779
	1780
	1781	IF ( .NOT. lfirst ) RETURN
	1782
	1783	ALLOCATE( even_odd_level(maximum_grid_level) )
	1784
	1785	ALLOCATE( f1_mg_b(nzb:nzt+1,maximum_grid_level), &
	1786	f2_mg_b(nzb:nzt+1,maximum_grid_level), &
[2232]	1787	f3_mg_b(nzb:nzt+1,maximum_grid_level) )
[1575]	1788
	1789	!
	1790	!-- Set border index between the even and odd block
	1791	DO i = maximum_grid_level, 1, -1
	1792	even_odd_level(i) = nzt_mg(i) / 2
	1793	ENDDO
	1794
	1795	!
	1796	!-- Sort grid coefficients used in red/black scheme and for calculating the
	1797	!-- residual to block (even/odd) structure
	1798	DO l = maximum_grid_level, 1 , -1
	1799	CALL sort_k_to_even_odd_blocks( f1_mg(nzb+1:nzt_mg(grid_level),l), &
	1800	f1_mg_b(nzb:nzt_mg(grid_level)+1,l), &
	1801	l )
	1802	CALL sort_k_to_even_odd_blocks( f2_mg(nzb+1:nzt_mg(grid_level),l), &
	1803	f2_mg_b(nzb:nzt_mg(grid_level)+1,l), &
	1804	l )
	1805	CALL sort_k_to_even_odd_blocks( f3_mg(nzb+1:nzt_mg(grid_level),l), &
	1806	f3_mg_b(nzb:nzt_mg(grid_level)+1,l), &
	1807	l )
	1808	ENDDO
	1809
	1810	lfirst = .FALSE.
	1811
	1812	END SUBROUTINE init_even_odd_blocks
	1813
	1814
	1815	!------------------------------------------------------------------------------!
	1816	! Description:
	1817	! ------------
[1682]	1818	!> Special exchange_horiz subroutine for use in redblack. Transfers only
	1819	!> "red" or "black" data points.
[1575]	1820	!------------------------------------------------------------------------------!
[1682]	1821	SUBROUTINE special_exchange_horiz ( p_mg, color )
[1575]	1822
[1682]	1823
[1575]	1824	USE control_parameters, &
	1825	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, &
[2696]	1826	maximum_grid_level, &
	1827	mg_switch_to_pe0_level, synchronous_exchange
[1575]	1828
	1829	USE indices, &
	1830	ONLY: mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn, &
	1831	nyn_mg, nzb, nzt, nzt_mg
	1832
	1833	IMPLICIT NONE
	1834
	1835	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1836	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1904]	1837	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: &
	1838	p_mg !< treated array
[1575]	1839
[1904]	1840	INTEGER(iwp), intent(IN) :: color !< flag for grid point type (red or black)
[1575]	1841	!
	1842	!-- Local variables
[1904]	1843	INTEGER(iwp) :: i !< index variable along x
	1844	INTEGER(iwp) :: i1 !< index variable along x on coarse level
	1845	INTEGER(iwp) :: i2 !< index variable along x on coarse level
[1575]	1846
[1904]	1847	INTEGER(iwp) :: j !< index variable along y
	1848	INTEGER(iwp) :: j1 !< index variable along y on coarse level
	1849	INTEGER(iwp) :: j2 !< index variable along y on coarse level
	1850	INTEGER(iwp) :: k !< index variable along z
	1851	INTEGER(iwp) :: l !< short for grid level
	1852	INTEGER(iwp) :: jys !< index for lower local PE boundary along y
	1853	INTEGER(iwp) :: jyn !< index for upper local PE boundary along y
	1854	INTEGER(iwp) :: ixl !< index for lower local PE boundary along x
	1855	INTEGER(iwp) :: ixr !< index for upper local PE boundary along x
[1575]	1856
[1904]	1857	LOGICAL :: synchronous_exchange_save !< dummy to reset synchronous_exchange to prescribed value
	1858
	1859	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: temp !< temporary array on next coarser grid level
	1860
[1609]	1861	#if defined ( __parallel )
[1575]	1862	synchronous_exchange_save = synchronous_exchange
	1863	synchronous_exchange = .FALSE.
	1864
	1865	l = grid_level
	1866
	1867	ind_even_odd = even_odd_level(grid_level)
	1868
[1904]	1869	!
	1870	!-- Restricted transfer only on finer levels with enough data.
	1871	!-- Restricted transfer is not possible for levels smaller or equal to
	1872	!-- 'switch to PE0 levels', since array bounds does not fit. Moreover,
	1873	!-- it is not possible for the coarsest grid level, since the dimensions
	1874	!-- of temp are not defined. For such cases, normal exchange_horiz is called.
	1875	IF ( l > 1 .AND. l > mg_switch_to_pe0_level + 1 .AND. &
	1876	( ngp_xz(grid_level) >= 900 .OR. ngp_yz(grid_level) >= 900 ) ) THEN
[1575]	1877
[1904]	1878	jys = nys_mg(grid_level-1)
	1879	jyn = nyn_mg(grid_level-1)
	1880	ixl = nxl_mg(grid_level-1)
	1881	ixr = nxr_mg(grid_level-1)
	1882	ALLOCATE( temp(nzb:nzt_mg(l-1)+1,jys-1:jyn+1,ixl-1:ixr+1) )
[1575]	1883	!
	1884	!-- Handling the even k Values
	1885	!-- Collecting data for the north - south exchange
	1886	!-- Since only every second value has to be transfered, data are stored
	1887	!-- on the next coarser grid level, because the arrays on that level
	1888	!-- have just the required size
	1889	i1 = nxl_mg(grid_level-1)
	1890	i2 = nxl_mg(grid_level-1)
	1891
	1892	DO i = nxl_mg(l), nxr_mg(l), 2
	1893	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	1894
	1895	IF ( j == nys_mg(l) ) THEN
	1896	!DIR$ IVDEP
	1897	DO k = ind_even_odd+1, nzt_mg(l)
	1898	temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
	1899	ENDDO
	1900	i1 = i1 + 1
	1901
	1902	ENDIF
	1903
	1904	IF ( j == nyn_mg(l) ) THEN
	1905	!DIR$ IVDEP
	1906	DO k = ind_even_odd+1, nzt_mg(l)
	1907	temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
	1908	ENDDO
	1909	i2 = i2 + 1
	1910
	1911	ENDIF
	1912
	1913	ENDDO
	1914	ENDDO
	1915
	1916	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	1917	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	1918
	1919	IF ( j == nys_mg(l) ) THEN
	1920	!DIR$ IVDEP
	1921	DO k = ind_even_odd+1, nzt_mg(l)
	1922	temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
	1923	ENDDO
	1924	i1 = i1 + 1
	1925
	1926	ENDIF
	1927
	1928	IF ( j == nyn_mg(l) ) THEN
	1929	!DIR$ IVDEP
	1930	DO k = ind_even_odd+1, nzt_mg(l)
	1931	temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
	1932	ENDDO
	1933	i2 = i2 + 1
	1934
	1935	ENDIF
	1936
	1937	ENDDO
	1938	ENDDO
	1939
	1940	grid_level = grid_level-1
	1941
	1942	nxl = nxl_mg(grid_level)
	1943	nys = nys_mg(grid_level)
	1944	nxr = nxr_mg(grid_level)
	1945	nyn = nyn_mg(grid_level)
	1946	nzt = nzt_mg(grid_level)
	1947
	1948	send_receive = 'ns'
	1949	CALL exchange_horiz( temp, 1 )
	1950
	1951	grid_level = grid_level+1
	1952
	1953	i1 = nxl_mg(grid_level-1)
	1954	i2 = nxl_mg(grid_level-1)
	1955
	1956	DO i = nxl_mg(l), nxr_mg(l), 2
	1957	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	1958
	1959	IF ( j == nys_mg(l) ) THEN
	1960	!DIR$ IVDEP
	1961	DO k = ind_even_odd+1, nzt_mg(l)
	1962	p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
	1963	ENDDO
	1964	i1 = i1 + 1
	1965
	1966	ENDIF
	1967
	1968	IF ( j == nyn_mg(l) ) THEN
	1969	!DIR$ IVDEP
	1970	DO k = ind_even_odd+1, nzt_mg(l)
	1971	p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
	1972	ENDDO
	1973	i2 = i2 + 1
	1974
	1975	ENDIF
	1976
	1977	ENDDO
	1978	ENDDO
	1979
	1980	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	1981	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	1982
	1983	IF ( j == nys_mg(l) ) THEN
	1984	!DIR$ IVDEP
	1985	DO k = ind_even_odd+1, nzt_mg(l)
	1986	p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
	1987	ENDDO
	1988	i1 = i1 + 1
	1989
	1990	ENDIF
	1991
	1992	IF ( j == nyn_mg(l) ) THEN
	1993	!DIR$ IVDEP
	1994	DO k = ind_even_odd+1, nzt_mg(l)
	1995	p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
	1996	ENDDO
	1997	i2 = i2 + 1
	1998
	1999	ENDIF
	2000
	2001	ENDDO
	2002	ENDDO
	2003
	2004	!
	2005	!-- Collecting data for the left - right exchange
	2006	!-- Since only every second value has to be transfered, data are stored
	2007	!-- on the next coarser grid level, because the arrays on that level
	2008	!-- have just the required size
	2009	j1 = nys_mg(grid_level-1)
	2010	j2 = nys_mg(grid_level-1)
	2011
	2012	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	2013	DO i = nxl_mg(l), nxr_mg(l), 2
	2014
	2015	IF ( i == nxl_mg(l) ) THEN
	2016	!DIR$ IVDEP
	2017	DO k = ind_even_odd+1, nzt_mg(l)
	2018	temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
	2019	ENDDO
	2020	j1 = j1 + 1
	2021
	2022	ENDIF
	2023
	2024	IF ( i == nxr_mg(l) ) THEN
	2025	!DIR$ IVDEP
	2026	DO k = ind_even_odd+1, nzt_mg(l)
	2027	temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
	2028	ENDDO
	2029	j2 = j2 + 1
	2030
	2031	ENDIF
	2032
	2033	ENDDO
	2034	ENDDO
	2035
	2036	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	2037	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	2038
	2039	IF ( i == nxl_mg(l) ) THEN
	2040	!DIR$ IVDEP
	2041	DO k = ind_even_odd+1, nzt_mg(l)
	2042	temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
	2043	ENDDO
	2044	j1 = j1 + 1
	2045
	2046	ENDIF
	2047
	2048	IF ( i == nxr_mg(l) ) THEN
	2049	!DIR$ IVDEP
	2050	DO k = ind_even_odd+1, nzt_mg(l)
	2051	temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
	2052	ENDDO
	2053	j2 = j2 + 1
	2054
	2055	ENDIF
	2056
	2057	ENDDO
	2058	ENDDO
	2059
	2060	grid_level = grid_level-1
	2061	send_receive = 'lr'
	2062
	2063	CALL exchange_horiz( temp, 1 )
	2064
	2065	grid_level = grid_level+1
	2066
	2067	j1 = nys_mg(grid_level-1)
	2068	j2 = nys_mg(grid_level-1)
	2069
	2070	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	2071	DO i = nxl_mg(l), nxr_mg(l), 2
	2072
	2073	IF ( i == nxl_mg(l) ) THEN
	2074	!DIR$ IVDEP
	2075	DO k = ind_even_odd+1, nzt_mg(l)
	2076	p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1)
	2077	ENDDO
	2078	j1 = j1 + 1
	2079
	2080	ENDIF
	2081
	2082	IF ( i == nxr_mg(l) ) THEN
	2083	!DIR$ IVDEP
	2084	DO k = ind_even_odd+1, nzt_mg(l)
	2085	p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1)
	2086	ENDDO
	2087	j2 = j2 + 1
	2088
	2089	ENDIF
	2090
	2091	ENDDO
	2092	ENDDO
	2093
	2094	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	2095	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	2096
	2097	IF ( i == nxl_mg(l) ) THEN
	2098	!DIR$ IVDEP
	2099	DO k = ind_even_odd+1, nzt_mg(l)
	2100	p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1)
	2101	ENDDO
	2102	j1 = j1 + 1
	2103
	2104	ENDIF
	2105
	2106	IF ( i == nxr_mg(l) ) THEN
	2107	!DIR$ IVDEP
	2108	DO k = ind_even_odd+1, nzt_mg(l)
	2109	p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1)
	2110	ENDDO
	2111	j2 = j2 + 1
	2112
	2113	ENDIF
	2114
	2115	ENDDO
	2116	ENDDO
	2117
	2118	!
	2119	!-- Now handling the even k values
	2120	!-- Collecting data for the north - south exchange
	2121	!-- Since only every second value has to be transfered, data are stored
	2122	!-- on the next coarser grid level, because the arrays on that level
	2123	!-- have just the required size
	2124	i1 = nxl_mg(grid_level-1)
	2125	i2 = nxl_mg(grid_level-1)
	2126
	2127	DO i = nxl_mg(l), nxr_mg(l), 2
	2128	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	2129
	2130	IF ( j == nys_mg(l) ) THEN
	2131	!DIR$ IVDEP
	2132	DO k = nzb+1, ind_even_odd
	2133	temp(k,jys,i1) = p_mg(k,j,i)
	2134	ENDDO
	2135	i1 = i1 + 1
	2136
	2137	ENDIF
	2138
	2139	IF ( j == nyn_mg(l) ) THEN
	2140	!DIR$ IVDEP
	2141	DO k = nzb+1, ind_even_odd
	2142	temp(k,jyn,i2) = p_mg(k,j,i)
	2143	ENDDO
	2144	i2 = i2 + 1
	2145
	2146	ENDIF
	2147
	2148	ENDDO
	2149	ENDDO
	2150
	2151	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	2152	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	2153
	2154	IF ( j == nys_mg(l) ) THEN
	2155	!DIR$ IVDEP
	2156	DO k = nzb+1, ind_even_odd
	2157	temp(k,jys,i1) = p_mg(k,j,i)
	2158	ENDDO
	2159	i1 = i1 + 1
	2160
	2161	ENDIF
	2162
	2163	IF ( j == nyn_mg(l) ) THEN
	2164	!DIR$ IVDEP
	2165	DO k = nzb+1, ind_even_odd
	2166	temp(k,jyn,i2) = p_mg(k,j,i)
	2167	ENDDO
	2168	i2 = i2 + 1
	2169
	2170	ENDIF
	2171
	2172	ENDDO
	2173	ENDDO
	2174
	2175	grid_level = grid_level-1
	2176
	2177	send_receive = 'ns'
	2178	CALL exchange_horiz( temp, 1 )
	2179
	2180	grid_level = grid_level+1
	2181
	2182	i1 = nxl_mg(grid_level-1)
	2183	i2 = nxl_mg(grid_level-1)
	2184
	2185	DO i = nxl_mg(l), nxr_mg(l), 2
	2186	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	2187
	2188	IF ( j == nys_mg(l) ) THEN
	2189	!DIR$ IVDEP
	2190	DO k = nzb+1, ind_even_odd
	2191	p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
	2192	ENDDO
	2193	i1 = i1 + 1
	2194
	2195	ENDIF
	2196
	2197	IF ( j == nyn_mg(l) ) THEN
	2198	!DIR$ IVDEP
	2199	DO k = nzb+1, ind_even_odd
	2200	p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
	2201	ENDDO
	2202	i2 = i2 + 1
	2203
	2204	ENDIF
	2205
	2206	ENDDO
	2207	ENDDO
	2208
	2209	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	2210	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	2211
	2212	IF ( j == nys_mg(l) ) THEN
	2213	!DIR$ IVDEP
	2214	DO k = nzb+1, ind_even_odd
	2215	p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
	2216	ENDDO
	2217	i1 = i1 + 1
	2218
	2219	ENDIF
	2220
	2221	IF ( j == nyn_mg(l) ) THEN
	2222	!DIR$ IVDEP
	2223	DO k = nzb+1, ind_even_odd
	2224	p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
	2225	ENDDO
	2226	i2 = i2 + 1
	2227
	2228	ENDIF
	2229
	2230	ENDDO
	2231	ENDDO
	2232
	2233	j1 = nys_mg(grid_level-1)
	2234	j2 = nys_mg(grid_level-1)
	2235
	2236	DO i = nxl_mg(l), nxr_mg(l), 2
	2237	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	2238
	2239	IF ( i == nxl_mg(l) ) THEN
	2240	!DIR$ IVDEP
	2241	DO k = nzb+1, ind_even_odd
	2242	temp(k,j1,ixl) = p_mg(k,j,i)
	2243	ENDDO
	2244	j1 = j1 + 1
	2245
	2246	ENDIF
	2247
	2248	IF ( i == nxr_mg(l) ) THEN
	2249	!DIR$ IVDEP
	2250	DO k = nzb+1, ind_even_odd
	2251	temp(k,j2,ixr) = p_mg(k,j,i)
	2252	ENDDO
	2253	j2 = j2 + 1
	2254
	2255	ENDIF
	2256
	2257	ENDDO
	2258	ENDDO
	2259
	2260	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	2261	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	2262
	2263	IF ( i == nxl_mg(l) ) THEN
	2264	!DIR$ IVDEP
	2265	DO k = nzb+1, ind_even_odd
	2266	temp(k,j1,ixl) = p_mg(k,j,i)
	2267	ENDDO
	2268	j1 = j1 + 1
	2269
	2270	ENDIF
	2271
	2272	IF ( i == nxr_mg(l) ) THEN
	2273	!DIR$ IVDEP
	2274	DO k = nzb+1, ind_even_odd
	2275	temp(k,j2,ixr) = p_mg(k,j,i)
	2276	ENDDO
	2277	j2 = j2 + 1
	2278
	2279	ENDIF
	2280
	2281	ENDDO
	2282	ENDDO
	2283
	2284	grid_level = grid_level-1
	2285
	2286	send_receive = 'lr'
	2287	CALL exchange_horiz( temp, 1 )
	2288
	2289	grid_level = grid_level+1
	2290
	2291	nxl = nxl_mg(grid_level)
	2292	nys = nys_mg(grid_level)
	2293	nxr = nxr_mg(grid_level)
	2294	nyn = nyn_mg(grid_level)
	2295	nzt = nzt_mg(grid_level)
	2296
	2297	j1 = nys_mg(grid_level-1)
	2298	j2 = nys_mg(grid_level-1)
	2299
	2300	DO i = nxl_mg(l), nxr_mg(l), 2
	2301	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	2302
	2303	IF ( i == nxl_mg(l) ) THEN
	2304	!DIR$ IVDEP
	2305	DO k = nzb+1, ind_even_odd
	2306	p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1)
	2307	ENDDO
	2308	j1 = j1 + 1
	2309
	2310	ENDIF
	2311
	2312	IF ( i == nxr_mg(l) ) THEN
	2313	!DIR$ IVDEP
	2314	DO k = nzb+1, ind_even_odd
	2315	p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1)
	2316	ENDDO
	2317	j2 = j2 + 1
	2318
	2319	ENDIF
	2320
	2321	ENDDO
	2322	ENDDO
	2323
	2324	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	2325	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	2326
	2327	IF ( i == nxl_mg(l) ) THEN
	2328	!DIR$ IVDEP
	2329	DO k = nzb+1, ind_even_odd
	2330	p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1)
	2331	ENDDO
	2332	j1 = j1 + 1
	2333
	2334	ENDIF
	2335
	2336	IF ( i == nxr_mg(l) ) THEN
	2337	!DIR$ IVDEP
	2338	DO k = nzb+1, ind_even_odd
	2339	p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1)
	2340	ENDDO
	2341	j2 = j2 + 1
	2342
	2343	ENDIF
	2344
	2345	ENDDO
	2346	ENDDO
	2347
[1904]	2348	DEALLOCATE( temp )
	2349
[1575]	2350	ELSE
[1609]	2351
[1575]	2352	!
	2353	!-- Standard horizontal ghost boundary exchange for small coarse grid
	2354	!-- levels, where the transfer time is latency bound
	2355	CALL exchange_horiz( p_mg, 1 )
	2356
	2357	ENDIF
	2358
	2359	!
	2360	!-- Reset values to default PALM setup
	2361	synchronous_exchange = synchronous_exchange_save
	2362	send_receive = 'al'
[1609]	2363	#else
[1575]	2364
[1609]	2365	!
	2366	!-- Standard horizontal ghost boundary exchange for small coarse grid
	2367	!-- levels, where the transfer time is latency bound
	2368	CALL exchange_horiz( p_mg, 1 )
	2369	#endif
	2370
[1575]	2371	END SUBROUTINE special_exchange_horiz
	2372
	2373	END MODULE poismg_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |