Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |