Home

Context Navigation

source: palm/trunk/SOURCE/poismg_noopt_mod.f90 @ 4798

Last change on this file since 4798 was 4649, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 93.7 KB

Rev	Line
[2586]	1	!> @file poismg_noopt_mod.f90
[4649]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[2586]	4	!
[4649]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[2586]	8	!
[4649]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[2586]	12	!
[4649]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[2586]	15	!
[4360]	16	! Copyright 1997-2020 Leibniz Universitaet Hannover
[4649]	17	!--------------------------------------------------------------------------------------------------!
[2586]	18	!
[4649]	19	!
[2586]	20	! Current revisions:
	21	! -----------------
[4649]	22	!
	23	!
[2591]	24	! Former revisions:
	25	! -----------------
	26	! $Id: poismg_noopt_mod.f90 4649 2020-08-25 12:11:17Z schwenkel $
[4649]	27	! File re-formatted to follow the PALM coding standard
	28	!
	29	!
	30	! 4457 2020-03-11 14:20:43Z raasch
	31	! Use statement for exchange horiz added
	32	!
[4457]	33	! 4429 2020-02-27 15:24:30Z raasch
[4649]	34	! Bugfix: cpp-directives added for serial mode
	35	!
[4429]	36	! 4414 2020-02-19 20:16:04Z suehring
[4414]	37	! Remove double-declared use only construct.
[4649]	38	!
[4414]	39	! 4360 2020-01-07 11:25:50Z suehring
[4346]	40	! Introduction of wall_flags_total_0, which currently sets bits based on static
	41	! topography information used in wall_flags_static_0
[4649]	42	!
[4346]	43	! 4329 2019-12-10 15:46:36Z motisi
[4329]	44	! Renamed wall_flags_0 to wall_flags_static_0
[4649]	45	!
[4329]	46	! 4182 2019-08-22 15:20:23Z scharf
[4182]	47	! Corrected "Former revisions" section
[4649]	48	!
[4182]	49	! 3655 2019-01-07 16:51:22Z knoop
[4649]	50	! Unused variables removed
[2586]	51	!
[4182]	52	! Revision 1.1 2001/07/20 13:10:51 raasch
	53	! Initial revision
	54	!
[4649]	55	!--------------------------------------------------------------------------------------------------!
[2586]	56	! Description:
	57	! ------------
[4649]	58	!> Solves the Poisson equation for the perturbation pressure with a multigrid V- or W-Cycle scheme.
[2586]	59	!>
	60	!> This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
	61	!> September 2000 - July 2001.
	62	!>
[4649]	63	!> @attention Loop unrolling and cache optimization in SOR-Red/Black method still does not give the
	64	!> expected speedup!
	65	!>
[2586]	66	!> @todo Further work required.
[4649]	67	!> @todo Formatting adjustments required (indention after modularization)
	68	!--------------------------------------------------------------------------------------------------!
[2586]	69	MODULE poismg_noopt_mod
	70
[4649]	71	USE control_parameters, &
	72	ONLY: bc_dirichlet_l, &
	73	bc_dirichlet_n, &
	74	bc_dirichlet_r, &
	75	bc_dirichlet_s, &
	76	bc_radiation_l, &
	77	bc_radiation_n, &
	78	bc_radiation_r, &
	79	bc_radiation_s, &
	80	child_domain, &
	81	grid_level, &
	82	nesting_offline
[2586]	83
[4649]	84	USE cpulog, &
	85	ONLY: cpu_log, &
	86	log_point_s
[4457]	87
[4649]	88	USE exchange_horiz_mod, &
	89	ONLY: exchange_horiz, &
	90	exchange_horiz_int
	91
[2586]	92	USE kinds
	93
	94	USE pegrid
	95
	96	PRIVATE
	97
	98	INTERFACE poismg_noopt
	99	MODULE PROCEDURE poismg_noopt
	100	END INTERFACE poismg_noopt
	101
	102	INTERFACE poismg_noopt_init
	103	MODULE PROCEDURE poismg_noopt_init
	104	END INTERFACE poismg_noopt_init
	105
	106	PUBLIC poismg_noopt, poismg_noopt_init
	107
	108	CONTAINS
	109
	110	SUBROUTINE poismg_noopt( r )
	111
[4649]	112	USE arrays_3d, &
	113	ONLY: d, &
	114	p_loc
[2586]	115
[4649]	116	USE control_parameters, &
	117	ONLY: bc_lr_cyc, &
	118	bc_ns_cyc, &
	119	gathered_size, &
	120	grid_level, &
	121	grid_level_count, &
	122	ibc_p_t, &
	123	maximum_grid_level, &
	124	message_string, &
	125	mgcycles, &
	126	mg_cycles, &
	127	mg_switch_to_pe0_level, &
	128	residual_limit, &
	129	subdomain_size
[2586]	130
[4649]	131	USE cpulog, &
	132	ONLY: cpu_log, &
	133	log_point_s
[2586]	134
[4649]	135	USE indices, &
	136	ONLY: nxl, &
	137	nxlg, &
	138	nxl_mg, &
	139	nxr, &
	140	nxrg, &
	141	nxr_mg, &
	142	nys, &
	143	nysg, &
	144	nys_mg, &
	145	nyn, &
	146	nyng, &
	147	nyn_mg, &
	148	nzb, &
	149	nzt, &
	150	nzt_mg
[2586]	151
	152	USE kinds
	153
	154	USE pegrid
	155
	156	IMPLICIT NONE
	157
	158	REAL(wp) :: maxerror !<
	159	REAL(wp) :: maximum_mgcycles !<
	160	REAL(wp) :: residual_norm !<
	161
	162	REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !<
	163
	164	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !<
	165
	166
	167	CALL cpu_log( log_point_s(29), 'poismg_noopt', 'start' )
	168	!
	169	!-- Initialize arrays and variables used in this subroutine
	170
[4649]	171	!-- If the number of grid points of the gathered grid, which is collected on PE0, is larger than
	172	!-- the number of grid points of an PE, than array p3 will be enlarged.
	173	IF ( gathered_size > subdomain_size ) THEN
	174	ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
	175	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1, &
	176	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(mg_switch_to_pe0_level)+1) )
[2586]	177	ELSE
	178	ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	179	ENDIF
	180
	181	p3 = 0.0_wp
[4649]	182
[2586]	183	!
	184	!-- Ghost boundaries have to be added to divergence array.
	185	!-- Exchange routine needs to know the grid level!
	186	grid_level = maximum_grid_level
	187	CALL exchange_horiz( d, 1)
[2939]	188	!
	189	!-- Set bottom and top boundary conditions
[2586]	190	d(nzb,:,:) = d(nzb+1,:,:)
[2939]	191	IF ( ibc_p_t == 1 ) d(nzt+1,:,: ) = d(nzt,:,:)
	192	!
	193	!-- Set lateral boundary conditions in non-cyclic case
	194	IF ( .NOT. bc_lr_cyc ) THEN
[4649]	195	IF ( bc_dirichlet_l .OR. bc_radiation_l ) d(:,:,nxl-1) = d(:,:,nxl)
	196	IF ( bc_dirichlet_r .OR. bc_radiation_r ) d(:,:,nxr+1) = d(:,:,nxr)
[2939]	197	ENDIF
	198	IF ( .NOT. bc_ns_cyc ) THEN
[4649]	199	IF ( bc_dirichlet_n .OR. bc_radiation_n ) d(:,nyn+1,:) = d(:,nyn,:)
	200	IF ( bc_dirichlet_s .OR. bc_radiation_s ) d(:,nys-1,:) = d(:,nys,:)
[2939]	201	ENDIF
[2586]	202
	203	!
[4649]	204	!-- Initiation of the multigrid scheme. Does n cycles until the residual is smaller than the
	205	!-- given limit. The accuracy of the solution of the poisson equation will increase with the
	206	!-- number of cycles. If the number of cycles is preset by the user, this number will be carried
	207	!-- out regardless of the accuracy.
[2586]	208	grid_level_count = 0
	209	mgcycles = 0
	210	IF ( mg_cycles == -1 ) THEN
	211	maximum_mgcycles = 0
	212	residual_norm = 1.0_wp
	213	ELSE
	214	maximum_mgcycles = mg_cycles
	215	residual_norm = 0.0_wp
	216	ENDIF
	217
[4649]	218	DO WHILE ( residual_norm > residual_limit .OR. mgcycles < maximum_mgcycles )
	219
[2586]	220	CALL next_mg_level_noopt( d, p_loc, p3, r)
	221
	222	!
[4649]	223	!-- Calculate the residual if the user has not preset the number of cycles to be performed
[2586]	224	IF ( maximum_mgcycles == 0 ) THEN
	225	CALL resid_noopt( d, p_loc, r )
	226	maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
	227
	228	#if defined( __parallel )
	229	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4649]	230	CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
[2586]	231	#else
	232	residual_norm = maxerror
	233	#endif
	234	residual_norm = SQRT( residual_norm )
	235	ENDIF
	236
	237	mgcycles = mgcycles + 1
	238
	239	!
[4649]	240	!-- If the user has not limited the number of cycles, stop the run in case of insufficient
	241	!-- convergence
[2586]	242	IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN
	243	message_string = 'no sufficient convergence within 1000 cycles'
	244	CALL message( 'poismg_noopt', 'PA0283', 1, 2, 0, 6, 0 )
	245	ENDIF
	246
	247	ENDDO
	248
	249	DEALLOCATE( p3 )
	250
	251	!
[4649]	252	!-- Unset the grid level. Variable is used to determine the MPI datatypes for ghost point
	253	!-- exchange
	254	grid_level = 0
[2586]	255
	256	CALL cpu_log( log_point_s(29), 'poismg_noopt', 'stop' )
	257
	258	END SUBROUTINE poismg_noopt
	259
	260
[4649]	261	!--------------------------------------------------------------------------------------------------!
[2586]	262	! Description:
	263	! ------------
	264	!> Computes the residual of the perturbation pressure.
[4649]	265	!--------------------------------------------------------------------------------------------------!
	266	SUBROUTINE resid_noopt( f_mg, p_mg, r )
[2586]	267
	268
[4649]	269	USE arrays_3d, &
	270	ONLY: f1_mg, &
	271	f2_mg, &
	272	f3_mg, &
	273	rho_air_mg
[2586]	274
[4649]	275	USE control_parameters, &
	276	ONLY: bc_lr_cyc, &
	277	bc_ns_cyc, &
	278	ibc_p_b, &
	279	ibc_p_t
[2586]	280
[4649]	281	USE grid_variables, &
	282	ONLY: ddx2_mg, &
	283	ddy2_mg
[2586]	284
[4649]	285	USE indices, &
	286	ONLY: flags, &
	287	nxl_mg, &
	288	nxr_mg, &
	289	nys_mg, &
	290	nyn_mg, &
	291	nzb, &
	292	nzt_mg, &
	293	wall_flags_1, &
	294	wall_flags_2, &
	295	wall_flags_3, &
	296	wall_flags_4, &
	297	wall_flags_5, &
	298	wall_flags_6, &
	299	wall_flags_7, &
	300	wall_flags_8, &
	301	wall_flags_9, &
	302	wall_flags_10
[2586]	303
	304
	305
	306
[4649]	307	USE kinds
[2586]	308
[4649]	309	IMPLICIT NONE
	310
	311	INTEGER(iwp) :: i !<
	312	INTEGER(iwp) :: j !<
	313	INTEGER(iwp) :: k !<
	314	INTEGER(iwp) :: l !<
	315
	316	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	317	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
	318	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	319	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
	320	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	321	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
	322
[2586]	323	!
[4649]	324	!-- Calculate the residual
	325	l = grid_level
[2586]	326
	327	!
[4649]	328	!-- Choose flag array of this level
	329	SELECT CASE ( l )
	330	CASE ( 1 )
	331	flags => wall_flags_1
	332	CASE ( 2 )
	333	flags => wall_flags_2
	334	CASE ( 3 )
	335	flags => wall_flags_3
	336	CASE ( 4 )
	337	flags => wall_flags_4
	338	CASE ( 5 )
	339	flags => wall_flags_5
	340	CASE ( 6 )
	341	flags => wall_flags_6
	342	CASE ( 7 )
	343	flags => wall_flags_7
	344	CASE ( 8 )
	345	flags => wall_flags_8
	346	CASE ( 9 )
	347	flags => wall_flags_9
	348	CASE ( 10 )
	349	flags => wall_flags_10
	350	END SELECT
[2586]	351
	352	!$OMP PARALLEL PRIVATE (i,j,k)
	353	!$OMP DO
[4649]	354	DO i = nxl_mg(l), nxr_mg(l)
	355	DO j = nys_mg(l), nyn_mg(l)
	356	DO k = nzb+1, nzt_mg(l)
	357	r(k,j,i) = f_mg(k,j,i) - rho_air_mg(k,l) * ddx2_mg(l) * &
	358	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	359	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	360	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	361	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	362	- rho_air_mg(k,l) * ddy2_mg(l) * &
	363	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	364	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	365	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	366	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	367	- f2_mg(k,l) * &
	368	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	369	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	370	- f3_mg(k,l) * &
	371	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	372	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	373	+ f1_mg(k,l) * p_mg(k,j,i)
[2586]	374	!
[4649]	375	!-- Residual within topography should be zero
	376	r(k,j,i) = r(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
[2586]	377	ENDDO
	378	ENDDO
[4649]	379	ENDDO
[2586]	380	!$OMP END PARALLEL
	381
	382	!
[4649]	383	!-- Horizontal boundary conditions
	384	CALL exchange_horiz( r, 1)
[2586]	385
[4649]	386	IF ( .NOT. bc_lr_cyc ) THEN
	387	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
	388	r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
[2586]	389	ENDIF
[4649]	390	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
	391	r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
	392	ENDIF
	393	ENDIF
[2586]	394
[4649]	395	IF ( .NOT. bc_ns_cyc ) THEN
	396	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
	397	r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
[2586]	398	ENDIF
[4649]	399	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
	400	r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
	401	ENDIF
	402	ENDIF
[2586]	403
	404	!
[4649]	405	!-- Boundary conditions at bottom and top of the domain.
	406	!-- These points are not handled by the above loop. Points may be within buildings, but that
	407	!-- doesn't matter.
	408	IF ( ibc_p_b == 1 ) THEN
	409	r(nzb,:,: ) = r(nzb+1,:,:)
	410	ELSE
	411	r(nzb,:,: ) = 0.0_wp
	412	ENDIF
[2586]	413
[4649]	414	IF ( ibc_p_t == 1 ) THEN
	415	r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
	416	ELSE
	417	r(nzt_mg(l)+1,:,: ) = 0.0_wp
	418	ENDIF
[2586]	419
	420
[4649]	421	END SUBROUTINE resid_noopt
[2586]	422
	423
[4649]	424	!--------------------------------------------------------------------------------------------------!
[2586]	425	! Description:
	426	! ------------
[4649]	427	!> Interpolates the residual on the next coarser grid with "full weighting" scheme.
	428	!--------------------------------------------------------------------------------------------------!
	429	SUBROUTINE restrict_noopt( f_mg, r )
[2586]	430
	431
[4649]	432	USE control_parameters, &
	433	ONLY: bc_lr_cyc, &
	434	bc_ns_cyc, &
	435	grid_level, &
	436	ibc_p_b, &
	437	ibc_p_t
[2586]	438
[4649]	439	USE indices, &
	440	ONLY: flags, &
	441	nxl_mg, &
	442	nxr_mg, &
	443	nys_mg, &
	444	nyn_mg, &
	445	nzb, &
	446	nzt_mg, &
	447	wall_flags_1, &
	448	wall_flags_2, &
	449	wall_flags_3, &
	450	wall_flags_4, &
	451	wall_flags_5, &
	452	wall_flags_6, &
	453	wall_flags_7, &
	454	wall_flags_8, &
	455	wall_flags_9, &
	456	wall_flags_10
[2586]	457
	458
[4649]	459	USE kinds
[2586]	460
[4649]	461	IMPLICIT NONE
[2586]	462
[4649]	463	INTEGER(iwp) :: i !<
	464	INTEGER(iwp) :: ic !<
	465	INTEGER(iwp) :: j !<
	466	INTEGER(iwp) :: jc !<
	467	INTEGER(iwp) :: k !<
	468	INTEGER(iwp) :: kc !<
	469	INTEGER(iwp) :: l !<
[2586]	470
[4649]	471	REAL(wp) :: rkjim !<
	472	REAL(wp) :: rkjip !<
	473	REAL(wp) :: rkjmi !<
	474	REAL(wp) :: rkjmim !<
	475	REAL(wp) :: rkjmip !<
	476	REAL(wp) :: rkjpi !<
	477	REAL(wp) :: rkjpim !<
	478	REAL(wp) :: rkjpip !<
	479	REAL(wp) :: rkmji !<
	480	REAL(wp) :: rkmjim !<
	481	REAL(wp) :: rkmjip !<
	482	REAL(wp) :: rkmjmi !<
	483	REAL(wp) :: rkmjmim !<
	484	REAL(wp) :: rkmjmip !<
	485	REAL(wp) :: rkmjpi !<
	486	REAL(wp) :: rkmjpim !<
	487	REAL(wp) :: rkmjpip !<
[2586]	488
[4649]	489	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	490	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[2586]	491
[4649]	492	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, &
	493	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r !<
	494
[2586]	495	!
[4649]	496	!-- Interpolate the residual
	497	l = grid_level
[2586]	498
	499	!
[4649]	500	!-- Choose flag array of the upper level
	501	SELECT CASE ( l+1 )
	502	CASE ( 1 )
	503	flags => wall_flags_1
	504	CASE ( 2 )
	505	flags => wall_flags_2
	506	CASE ( 3 )
	507	flags => wall_flags_3
	508	CASE ( 4 )
	509	flags => wall_flags_4
	510	CASE ( 5 )
	511	flags => wall_flags_5
	512	CASE ( 6 )
	513	flags => wall_flags_6
	514	CASE ( 7 )
	515	flags => wall_flags_7
	516	CASE ( 8 )
	517	flags => wall_flags_8
	518	CASE ( 9 )
	519	flags => wall_flags_9
	520	CASE ( 10 )
	521	flags => wall_flags_10
	522	END SELECT
	523
[2586]	524	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc, rkjim,rkjip,rkjpi,rkjmi,rkjmim,rkjpim, &
	525	!$OMP rkjmip, rkjpip,rkmji,rkmjim,rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim,rkmjmip,&
	526	!$OMP rkmjpip )
	527	!$OMP DO
[4649]	528	DO ic = nxl_mg(l), nxr_mg(l)
	529	i = 2*ic
	530	DO jc = nys_mg(l), nyn_mg(l)
	531	j = 2*jc
	532	DO kc = nzb+1, nzt_mg(l)
	533	k = 2*kc-1
[2586]	534	!
[4649]	535	!-- Use implicit Neumann BCs if the respective gridpoint is inside the building
	536	rkjim = r(k,j,i-1) + IBITS( flags(k,j,i-1), 6, 1 ) * ( r(k,j,i) - r(k,j,i-1) )
	537	rkjip = r(k,j,i+1) + IBITS( flags(k,j,i+1), 6, 1 ) * ( r(k,j,i) - r(k,j,i+1) )
	538	rkjpi = r(k,j+1,i) + IBITS( flags(k,j+1,i), 6, 1 ) * ( r(k,j,i) - r(k,j+1,i) )
	539	rkjmi = r(k,j-1,i) + IBITS( flags(k,j-1,i), 6, 1 ) * ( r(k,j,i) - r(k,j-1,i) )
	540	rkjmim = r(k,j-1,i-1) + IBITS( flags(k,j-1,i-1), 6, 1 ) * ( r(k,j,i) - r(k,j-1,i-1) )
	541	rkjpim = r(k,j+1,i-1) + IBITS( flags(k,j+1,i-1), 6, 1 ) * ( r(k,j,i) - r(k,j+1,i-1) )
	542	rkjmip = r(k,j-1,i+1) + IBITS( flags(k,j-1,i+1), 6, 1 ) * ( r(k,j,i) - r(k,j-1,i+1) )
	543	rkjpip = r(k,j+1,i+1) + IBITS( flags(k,j+1,i+1), 6, 1 ) * ( r(k,j,i) - r(k,j+1,i+1) )
	544	rkmji = r(k-1,j,i) + IBITS( flags(k-1,j,i), 6, 1 ) * ( r(k,j,i) - r(k-1,j,i) )
	545	rkmjim = r(k-1,j,i-1) + IBITS( flags(k-1,j,i-1), 6, 1 ) * ( r(k,j,i) - r(k-1,j,i-1) )
	546	rkmjip = r(k-1,j,i+1) + IBITS( flags(k-1,j,i+1), 6, 1 ) * ( r(k,j,i) - r(k-1,j,i+1) )
	547	rkmjpi = r(k-1,j+1,i) + IBITS( flags(k-1,j+1,i), 6, 1 ) * ( r(k,j,i) - r(k-1,j+1,i) )
	548	rkmjmi = r(k-1,j-1,i) + IBITS( flags(k-1,j-1,i), 6, 1 ) * ( r(k,j,i) - r(k-1,j-1,i) )
	549	rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) &
	550	* ( r(k,j,i) - r(k-1,j-1,i-1) )
	551	rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) &
	552	* ( r(k,j,i) - r(k-1,j+1,i-1) )
	553	rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) &
	554	* ( r(k,j,i) - r(k-1,j-1,i+1) )
	555	rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) &
	556	* ( r(k,j,i) - r(k-1,j+1,i+1) )
[2586]	557
[4649]	558	f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp &
	559	* ( 8.0_wp * r(k,j,i) &
	560	+ 4.0_wp * ( rkjim + rkjip + rkjpi + rkjmi ) &
	561	+ 2.0_wp * ( rkjmim + rkjpim + rkjmip + rkjpip ) &
	562	+ 4.0_wp * rkmji &
	563	+ 2.0_wp * ( rkmjim + rkmjim + rkmjpi + rkmjmi ) &
	564	+ ( rkmjmim + rkmjpim + rkmjmip + rkmjpip ) &
	565	+ 4.0_wp * r(k+1,j,i) &
	566	+ 2.0_wp * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
	567	r(k+1,j+1,i) + r(k+1,j-1,i) ) &
	568	+ ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
	569	r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
	570	)
[2586]	571
[4649]	572	! f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( &
	573	! 8.0_wp * r(k,j,i) &
	574	! + 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + &
	575	! r(k,j+1,i) + r(k,j-1,i) ) &
	576	! + 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
	577	! r(k,j-1,i+1) + r(k,j+1,i+1) ) &
	578	! + 4.0_wp * r(k-1,j,i) &
	579	! + 2.0_wp * ( r(k-1,j,i-1) + r(k-1,j,i+1) + &
	580	! r(k-1,j+1,i) + r(k-1,j-1,i) ) &
	581	! + ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
	582	! r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
	583	! + 4.0_wp * r(k+1,j,i) &
	584	! + 2.0_wp * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
	585	! r(k+1,j+1,i) + r(k+1,j-1,i) ) &
	586	! + ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
	587	! r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
	588	! )
[2586]	589	ENDDO
	590	ENDDO
[4649]	591	ENDDO
[2586]	592	!$OMP END PARALLEL
	593
	594	!
[4649]	595	!-- Horizontal boundary conditions
	596	CALL exchange_horiz( f_mg, 1)
[2586]	597
[4649]	598	IF ( .NOT. bc_lr_cyc ) THEN
	599	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
	600	f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
[2586]	601	ENDIF
[4649]	602	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
	603	f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
	604	ENDIF
	605	ENDIF
[2586]	606
[4649]	607	IF ( .NOT. bc_ns_cyc ) THEN
	608	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
	609	f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
[2586]	610	ENDIF
[4649]	611	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
	612	f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
	613	ENDIF
	614	ENDIF
[2586]	615
	616	!
[4649]	617	!-- Boundary conditions at bottom and top of the domain.
	618	!-- These points are not handled by the above loop. Points may be within buildings, but that
	619	!-- doesn't matter.
	620	IF ( ibc_p_b == 1 ) THEN
	621	f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
	622	ELSE
	623	f_mg(nzb,:,: ) = 0.0_wp
	624	ENDIF
[2586]	625
[4649]	626	IF ( ibc_p_t == 1 ) THEN
	627	f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
	628	ELSE
	629	f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
	630	ENDIF
[2586]	631
	632
[4649]	633	END SUBROUTINE restrict_noopt
[2586]	634
	635
[4649]	636	!--------------------------------------------------------------------------------------------------!
[2586]	637	! Description:
	638	! ------------
[4649]	639	!> Interpolates the correction of the perturbation pressure to the next finer grid.
	640	!--------------------------------------------------------------------------------------------------!
[2586]	641	SUBROUTINE prolong_noopt( p, temp )
	642
	643
[4649]	644	USE control_parameters, &
	645	ONLY: bc_lr_cyc, &
	646	bc_ns_cyc, &
	647	ibc_p_b, &
	648	ibc_p_t
	649	USE indices, &
	650	ONLY: nxl_mg, &
	651	nxr_mg, &
	652	nys_mg, &
	653	nyn_mg, &
	654	nzb, &
	655	nzt_mg
[2586]	656
	657	USE kinds
	658
	659	IMPLICIT NONE
	660
	661	INTEGER(iwp) :: i !<
	662	INTEGER(iwp) :: j !<
	663	INTEGER(iwp) :: k !<
	664	INTEGER(iwp) :: l !<
	665
[4649]	666	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[2586]	667	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p !<
	668
[4649]	669	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	670	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp !<
[2586]	671
	672
	673	!
	674	!-- First, store elements of the coarser grid on the next finer grid
	675	l = grid_level
	676
	677	!$OMP PARALLEL PRIVATE (i,j,k)
	678	!$OMP DO
	679	DO i = nxl_mg(l-1), nxr_mg(l-1)
	680	DO j = nys_mg(l-1), nyn_mg(l-1)
	681	!CDIR NODEP
	682	DO k = nzb+1, nzt_mg(l-1)
	683	!
[4649]	684	!-- Points of the coarse grid are directly stored on the next finer grid
	685	temp(2k-1,2j,2*i) = p(k,j,i)
	686	!
[2586]	687	!-- Points between two coarse-grid points
	688	temp(2k-1,2j,2i+1) = 0.5_wp ( p(k,j,i) + p(k,j,i+1) )
	689	temp(2k-1,2j+1,2i) = 0.5_wp ( p(k,j,i) + p(k,j+1,i) )
	690	temp(2k,2j,2i) = 0.5_wp ( p(k,j,i) + p(k+1,j,i) )
	691	!
[4649]	692	!-- Points in the center of the planes stretched by four points of the coarse grid cube
	693	temp(2k-1,2j+1,2i+1) = 0.25_wp ( p(k,j,i) + p(k,j,i+1) + &
[2586]	694	p(k,j+1,i) + p(k,j+1,i+1) )
[4649]	695	temp(2k,2j,2i+1) = 0.25_wp ( p(k,j,i) + p(k,j,i+1) + &
[2586]	696	p(k+1,j,i) + p(k+1,j,i+1) )
[4649]	697	temp(2k,2j+1,2i) = 0.25_wp ( p(k,j,i) + p(k,j+1,i) + &
[2586]	698	p(k+1,j,i) + p(k+1,j+1,i) )
	699	!
	700	!-- Points in the middle of coarse grid cube
[4649]	701	temp(2k,2j+1,2i+1) = 0.125_wp ( p(k,j,i) + p(k,j,i+1) + p(k,j+1,i) + &
	702	p(k,j+1,i+1) + p(k+1,j,i) + p(k+1,j,i+1) + &
[2586]	703	p(k+1,j+1,i) + p(k+1,j+1,i+1) )
	704	ENDDO
	705	ENDDO
	706	ENDDO
[4649]	707	!$OMP END PARALLEL
	708
[2586]	709	!
	710	!-- Horizontal boundary conditions
	711	CALL exchange_horiz( temp, 1)
	712
	713	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	714	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[2586]	715	temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
	716	ENDIF
[3182]	717	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
[2586]	718	temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
	719	ENDIF
	720	ENDIF
	721
	722	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	723	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
[2586]	724	temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
	725	ENDIF
[3182]	726	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[2586]	727	temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
	728	ENDIF
	729	ENDIF
	730
	731	!
	732	!-- Bottom and top boundary conditions
	733	IF ( ibc_p_b == 1 ) THEN
	734	temp(nzb,:,: ) = temp(nzb+1,:,:)
	735	ELSE
	736	temp(nzb,:,: ) = 0.0_wp
	737	ENDIF
	738
	739	IF ( ibc_p_t == 1 ) THEN
	740	temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
	741	ELSE
	742	temp(nzt_mg(l)+1,:,: ) = 0.0_wp
	743	ENDIF
	744
[4649]	745
[2586]	746	END SUBROUTINE prolong_noopt
	747
	748
[4649]	749	!--------------------------------------------------------------------------------------------------!
[2586]	750	! Description:
	751	! ------------
[4649]	752	!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with 3D-Red-Black
	753	!> decomposition (GS-RB) is used.
	754	!--------------------------------------------------------------------------------------------------!
	755	SUBROUTINE redblack_noopt( f_mg, p_mg )
[2586]	756
	757
[4649]	758	USE arrays_3d, &
	759	ONLY: f1_mg, &
	760	f2_mg, &
	761	f3_mg, &
	762	rho_air_mg
[2586]	763
[4649]	764	USE control_parameters, &
	765	ONLY: bc_lr_cyc, &
	766	bc_ns_cyc, &
	767	ibc_p_b, &
	768	ibc_p_t, &
	769	ngsrb
[2586]	770
[4649]	771	USE cpulog, &
	772	ONLY: cpu_log, &
	773	log_point_s
[2586]	774
[4649]	775	USE grid_variables, &
	776	ONLY: ddx2_mg, &
	777	ddy2_mg
[2586]	778
[4649]	779	USE indices, &
	780	ONLY: flags, &
	781	nxl_mg, &
	782	nxr_mg, &
	783	nys_mg, &
	784	nyn_mg, &
	785	nzb, &
	786	nzt_mg, &
	787	wall_flags_1, &
	788	wall_flags_2, &
	789	wall_flags_3, &
	790	wall_flags_4, &
	791	wall_flags_5, &
	792	wall_flags_6, &
	793	wall_flags_7, &
	794	wall_flags_8, &
	795	wall_flags_9, &
	796	wall_flags_10
[2586]	797
	798
[4649]	799	USE kinds
[2586]	800
[4649]	801	IMPLICIT NONE
[2586]	802
[4649]	803	INTEGER(iwp) :: color !<
	804	INTEGER(iwp) :: i !<
	805	INTEGER(iwp) :: ic !<
	806	INTEGER(iwp) :: j !<
	807	INTEGER(iwp) :: jc !<
	808	INTEGER(iwp) :: jj !<
	809	INTEGER(iwp) :: k !<
	810	INTEGER(iwp) :: l !<
	811	INTEGER(iwp) :: n !<
[2586]	812
[4649]	813	LOGICAL :: unroll !<
[2586]	814
[4649]	815	REAL(wp) :: wall_left !<
	816	REAL(wp) :: wall_north !<
	817	REAL(wp) :: wall_right !<
	818	REAL(wp) :: wall_south !<
	819	REAL(wp) :: wall_total !<
	820	REAL(wp) :: wall_top !<
[2586]	821
[4649]	822	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	823	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
	824	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	825	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[2586]	826
[4649]	827	l = grid_level
	828
[2586]	829	!
[4649]	830	!-- Choose flag array of this level
	831	SELECT CASE ( l )
	832	CASE ( 1 )
	833	flags => wall_flags_1
	834	CASE ( 2 )
	835	flags => wall_flags_2
	836	CASE ( 3 )
	837	flags => wall_flags_3
	838	CASE ( 4 )
	839	flags => wall_flags_4
	840	CASE ( 5 )
	841	flags => wall_flags_5
	842	CASE ( 6 )
	843	flags => wall_flags_6
	844	CASE ( 7 )
	845	flags => wall_flags_7
	846	CASE ( 8 )
	847	flags => wall_flags_8
	848	CASE ( 9 )
	849	flags => wall_flags_9
	850	CASE ( 10 )
	851	flags => wall_flags_10
	852	END SELECT
[2586]	853
[4649]	854	unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. &
	855	MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
[2586]	856
[4649]	857	DO n = 1, ngsrb
[2586]	858
[4649]	859	DO color = 1, 2
[2586]	860
[4649]	861	IF ( .NOT. unroll ) THEN
[2586]	862
[4649]	863	CALL cpu_log( log_point_s(36), 'redblack_no_unroll_noopt', 'start' )
	864
[2586]	865	!
[4649]	866	!-- Without unrolling of loops, no cache optimization
	867	DO i = nxl_mg(l), nxr_mg(l), 2
	868	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	869	DO k = nzb+1, nzt_mg(l), 2
	870	! p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
	871	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	872	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	873	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	874	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	875	! )
[2586]	876
[4649]	877	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	878	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	879	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	880	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	881	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	882	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	883	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	884	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	885	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	886	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	887	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	888	+ f2_mg(k,l) * &
	889	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	890	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	891	+ f3_mg(k,l) * &
	892	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	893	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	894	- f_mg(k,j,i) &
	895	)
[2586]	896	ENDDO
	897	ENDDO
[4649]	898	ENDDO
	899
	900	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	901	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	902	DO k = nzb+1, nzt_mg(l), 2
	903	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	904	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	905	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	906	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	907	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	908	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	909	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	910	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	911	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	912	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	913	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	914	+ f2_mg(k,l) * &
	915	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	916	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	917	+ f3_mg(k,l) * &
	918	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	919	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	920	- f_mg(k,j,i) &
	921	)
[2586]	922	ENDDO
	923	ENDDO
[4649]	924	ENDDO
	925
	926	DO i = nxl_mg(l), nxr_mg(l), 2
	927	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
	928	DO k = nzb+2, nzt_mg(l), 2
	929	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	930	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	931	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	932	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	933	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	934	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	935	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	936	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	937	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	938	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	939	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	940	+ f2_mg(k,l) * &
	941	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	942	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	943	+ f3_mg(k,l) * &
	944	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	945	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	946	- f_mg(k,j,i) &
	947	)
[2586]	948	ENDDO
	949	ENDDO
[4649]	950	ENDDO
[2586]	951
[4649]	952	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	953	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
	954	DO k = nzb+2, nzt_mg(l), 2
	955	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	956	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	957	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	958	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	959	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	960	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	961	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	962	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	963	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	964	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	965	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	966	+ f2_mg(k,l) * &
	967	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	968	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	969	+ f3_mg(k,l) * &
	970	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	971	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	972	- f_mg(k,j,i) &
	973	)
[2586]	974	ENDDO
	975	ENDDO
[4649]	976	ENDDO
	977	CALL cpu_log( log_point_s(36), 'redblack_no_unroll_noopt', 'stop' )
[2586]	978
[4649]	979	ELSE
[2586]	980
	981	!
[4649]	982	!-- Loop unrolling along y, only one i loop for better cache use
	983	CALL cpu_log( log_point_s(38), 'redblack_unroll_noopt', 'start' )
	984	DO ic = nxl_mg(l), nxr_mg(l), 2
	985	DO jc = nys_mg(l), nyn_mg(l), 4
	986	i = ic
	987	jj = jc+2-color
	988	DO k = nzb+1, nzt_mg(l), 2
	989	j = jj
	990	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	991	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	992	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	993	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	994	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	995	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	996	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	997	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	998	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	999	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1000	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1001	+ f2_mg(k,l) * &
	1002	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1003	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1004	+ f3_mg(k,l) * &
	1005	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1006	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1007	- f_mg(k,j,i) &
	1008	)
[2586]	1009
[4649]	1010	j = jj+2
	1011	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	1012	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	1013	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1014	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1015	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1016	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1017	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	1018	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1019	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1020	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1021	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1022	+ f2_mg(k,l) * &
	1023	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1024	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1025	+ f3_mg(k,l) * &
	1026	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1027	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1028	- f_mg(k,j,i) &
	1029	)
	1030	ENDDO
[2586]	1031
[4649]	1032	i = ic+1
	1033	jj = jc+color-1
	1034	DO k = nzb+1, nzt_mg(l), 2
	1035	j =jj
	1036	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	1037	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	1038	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1039	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1040	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1041	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1042	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	1043	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1044	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1045	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1046	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1047	+ f2_mg(k,l) * &
	1048	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1049	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1050	+ f3_mg(k,l) * &
	1051	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1052	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1053	- f_mg(k,j,i) &
	1054	)
[2586]	1055
[4649]	1056	j = jj+2
	1057	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	1058	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	1059	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1060	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1061	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1062	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1063	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	1064	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1065	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1066	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1067	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1068	+ f2_mg(k,l) * &
	1069	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1070	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1071	+ f3_mg(k,l) * &
	1072	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1073	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1074	- f_mg(k,j,i) &
	1075	)
[2586]	1076	ENDDO
[4649]	1077
	1078	i = ic
	1079	jj = jc+color-1
	1080	DO k = nzb+2, nzt_mg(l), 2
	1081	j =jj
	1082	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	1083	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	1084	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1085	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1086	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1087	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1088	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	1089	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1090	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1091	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1092	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1093	+ f2_mg(k,l) * &
	1094	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1095	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1096	+ f3_mg(k,l) * &
	1097	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1098	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1099	- f_mg(k,j,i) &
	1100	)
	1101
	1102	j = jj+2
	1103	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	1104	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	1105	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1106	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1107	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1108	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1109	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	1110	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1111	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1112	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1113	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1114	+ f2_mg(k,l) * &
	1115	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1116	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1117	+ f3_mg(k,l) * &
	1118	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1119	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1120	- f_mg(k,j,i) &
	1121	)
	1122	ENDDO
	1123
	1124	i = ic+1
	1125	jj = jc+2-color
	1126	DO k = nzb+2, nzt_mg(l), 2
	1127	j =jj
	1128	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	1129	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	1130	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1131	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1132	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1133	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1134	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	1135	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1136	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1137	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1138	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1139	+ f2_mg(k,l) * &
	1140	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1141	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1142	+ f3_mg(k,l) * &
	1143	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1144	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1145	- f_mg(k,j,i) &
	1146	)
	1147
	1148	j = jj+2
	1149	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) &
	1150	* ( rho_air_mg(k,l) * ddx2_mg(l) * &
	1151	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1152	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1153	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1154	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1155	+ rho_air_mg(k,l) * ddy2_mg(l) * &
	1156	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1157	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1158	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1159	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1160	+ f2_mg(k,l) * &
	1161	( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
	1162	( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
	1163	+ f3_mg(k,l) * &
	1164	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1165	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1166	- f_mg(k,j,i) &
	1167	)
	1168	ENDDO
	1169
[2586]	1170	ENDDO
[4649]	1171	ENDDO
	1172	CALL cpu_log( log_point_s(38), 'redblack_unroll_noopt', 'stop' )
[2586]	1173
[4649]	1174	ENDIF
[2586]	1175
	1176	!
[4649]	1177	!-- Horizontal boundary conditions
	1178	CALL exchange_horiz( p_mg, 1 )
[2586]	1179
[4649]	1180	IF ( .NOT. bc_lr_cyc ) THEN
	1181	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
	1182	p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
[2586]	1183	ENDIF
[4649]	1184	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
	1185	p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
	1186	ENDIF
	1187	ENDIF
[2586]	1188
[4649]	1189	IF ( .NOT. bc_ns_cyc ) THEN
	1190	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
	1191	p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
[2586]	1192	ENDIF
[4649]	1193	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
	1194	p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
	1195	ENDIF
	1196	ENDIF
[2586]	1197
	1198	!
[4649]	1199	!-- Bottom and top boundary conditions
	1200	IF ( ibc_p_b == 1 ) THEN
	1201	p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
	1202	ELSE
	1203	p_mg(nzb,:,: ) = 0.0_wp
	1204	ENDIF
[2586]	1205
[4649]	1206	IF ( ibc_p_t == 1 ) THEN
	1207	p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
	1208	ELSE
	1209	p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
	1210	ENDIF
[2586]	1211
	1212	ENDDO
	1213
[4649]	1214	ENDDO
	1215
[2586]	1216	!
	1217	!-- Set pressure within topography and at the topography surfaces
	1218	!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
	1219	!$OMP DO
[4649]	1220	DO i = nxl_mg(l), nxr_mg(l)
	1221	DO j = nys_mg(l), nyn_mg(l)
	1222	DO k = nzb, nzt_mg(l)
[2586]	1223	!
[4649]	1224	!-- First, set pressure inside topography to zero
	1225	p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
[2586]	1226	!
[4649]	1227	!-- Second, determine if the grid point inside topography is adjacent to a wall and set its
	1228	!-- value to a value given by the average of those values obtained from Neumann boundary
	1229	!-- condition.
	1230	wall_left = IBITS( flags(k,j,i-1), 5, 1 )
	1231	wall_right = IBITS( flags(k,j,i+1), 4, 1 )
	1232	wall_south = IBITS( flags(k,j-1,i), 3, 1 )
	1233	wall_north = IBITS( flags(k,j+1,i), 2, 1 )
	1234	wall_top = IBITS( flags(k+1,j,i), 0, 1 )
	1235	wall_total = wall_left + wall_right + wall_south + wall_north + wall_top
[2586]	1236
[4649]	1237	IF ( wall_total > 0.0_wp ) THEN
	1238	p_mg(k,j,i) = 1.0_wp / wall_total * ( wall_left * p_mg(k,j,i-1) + &
	1239	wall_right * p_mg(k,j,i+1) + &
	1240	wall_south * p_mg(k,j-1,i) + &
	1241	wall_north * p_mg(k,j+1,i) + &
	1242	wall_top * p_mg(k+1,j,i) )
	1243	ENDIF
[2586]	1244	ENDDO
	1245	ENDDO
[4649]	1246	ENDDO
[2586]	1247	!$OMP END PARALLEL
	1248
	1249	!
[4649]	1250	!-- One more time horizontal boundary conditions
	1251	CALL exchange_horiz( p_mg, 1)
[2586]	1252
	1253
[4649]	1254	END SUBROUTINE redblack_noopt
[2586]	1255
	1256
	1257
[4649]	1258	!--------------------------------------------------------------------------------------------------!
[2586]	1259	! Description:
	1260	! ------------
	1261	!> Gather subdomain data from all PEs.
[4649]	1262	!--------------------------------------------------------------------------------------------------!
[4429]	1263	#if defined( __parallel )
[4649]	1264	SUBROUTINE mg_gather_noopt( f2, f2_sub )
[2586]	1265
[4649]	1266	USE cpulog, &
	1267	ONLY: cpu_log, &
	1268	log_point_s
[2586]	1269
[4649]	1270	USE indices, &
	1271	ONLY: mg_loc_ind, &
	1272	nxl_mg, &
	1273	nxr_mg, &
	1274	nys_mg, &
	1275	nyn_mg, &
	1276	nzb, &
	1277	nzt_mg
[2586]	1278
[4649]	1279	USE kinds
[2586]	1280
[4649]	1281	USE pegrid
[2586]	1282
[4649]	1283	IMPLICIT NONE
[2586]	1284
[4649]	1285	INTEGER(iwp) :: i !<
	1286	INTEGER(iwp) :: il !<
	1287	INTEGER(iwp) :: ir !<
	1288	INTEGER(iwp) :: j !<
	1289	INTEGER(iwp) :: jn !<
	1290	INTEGER(iwp) :: js !<
	1291	INTEGER(iwp) :: k !<
	1292	INTEGER(iwp) :: nwords !<
[2586]	1293
[4649]	1294	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1295	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !<
	1296	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1297	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !<
[2586]	1298
[4649]	1299	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	1300	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !<
[2586]	1301
	1302
[4649]	1303	CALL cpu_log( log_point_s(34), 'mg_gather_noopt', 'start' )
[2586]	1304
[4649]	1305	f2_l = 0.0_wp
[2586]	1306
	1307	!
[4649]	1308	!-- Store the local subdomain array on the total array
	1309	js = mg_loc_ind(3,myid)
	1310	IF ( south_border_pe ) js = js - 1
	1311	jn = mg_loc_ind(4,myid)
	1312	IF ( north_border_pe ) jn = jn + 1
	1313	il = mg_loc_ind(1,myid)
	1314	IF ( left_border_pe ) il = il - 1
	1315	ir = mg_loc_ind(2,myid)
	1316	IF ( right_border_pe ) ir = ir + 1
	1317	DO i = il, ir
	1318	DO j = js, jn
	1319	DO k = nzb, nzt_mg(grid_level)+1
	1320	f2_l(k,j,i) = f2_sub(k,j,i)
[2586]	1321	ENDDO
	1322	ENDDO
[4649]	1323	ENDDO
[2586]	1324
	1325	!
[4649]	1326	!-- Find out the number of array elements of the total array
	1327	nwords = SIZE( f2 )
[2586]	1328
	1329	!
[4649]	1330	!-- Gather subdomain data from all PEs
	1331	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1332	CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
	1333	f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), nwords, MPI_REAL, &
	1334	MPI_SUM, comm2d, ierr )
[2586]	1335
[4649]	1336	CALL cpu_log( log_point_s(34), 'mg_gather_noopt', 'stop' )
	1337
	1338	END SUBROUTINE mg_gather_noopt
[4429]	1339	#endif
[2586]	1340
	1341
[4649]	1342	!--------------------------------------------------------------------------------------------------!
[2586]	1343	! Description:
	1344	! ------------
[4649]	1345	!> @todo It may be possible to improve the speed of this routine by using non-blocking communication
	1346	!--------------------------------------------------------------------------------------------------!
[4429]	1347	#if defined( __parallel )
[4649]	1348	SUBROUTINE mg_scatter_noopt( p2, p2_sub )
[2586]	1349
[4649]	1350	USE cpulog, &
	1351	ONLY: cpu_log, &
	1352	log_point_s
[2586]	1353
[4649]	1354	USE indices, &
	1355	ONLY: mg_loc_ind, &
	1356	nxl_mg, &
	1357	nxr_mg, &
	1358	nys_mg, &
	1359	nyn_mg, &
	1360	nzb, &
	1361	nzt_mg
[2586]	1362
[4649]	1363	USE kinds
[2586]	1364
[4649]	1365	USE pegrid
[2586]	1366
[4649]	1367	IMPLICIT NONE
[2586]	1368
[4649]	1369	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
	1370	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[2586]	1371
[4649]	1372	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	1373	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !<
[2586]	1374
	1375
[4649]	1376	CALL cpu_log( log_point_s(35), 'mg_scatter_noopt', 'start' )
[2586]	1377
[4649]	1378	p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	1379	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
[2586]	1380
[4649]	1381	CALL cpu_log( log_point_s(35), 'mg_scatter_noopt', 'stop' )
	1382
	1383	END SUBROUTINE mg_scatter_noopt
[4429]	1384	#endif
[2586]	1385
	1386
[4649]	1387	!--------------------------------------------------------------------------------------------------!
[2586]	1388	! Description:
	1389	! ------------
[4649]	1390	!> This is where the multigrid technique takes place. V- and W- Cycle are implemented and steered by
	1391	!> the parameter "gamma". Parameter "nue" determines the convergence of the multigrid iterative
	1392	!> solution. There are nue times RB-GS iterations. It should be set to "1" or "2", considering the
	1393	!> time effort one would like to invest. Last choice shows a very good converging factor, but leads
	1394	!> to an increase in computing time.
	1395	!--------------------------------------------------------------------------------------------------!
	1396	RECURSIVE SUBROUTINE next_mg_level_noopt( f_mg, p_mg, p3, r )
[2586]	1397
[4649]	1398	USE control_parameters, &
	1399	ONLY: bc_lr_dirrad, &
	1400	bc_lr_raddir, &
	1401	bc_ns_dirrad, &
	1402	bc_ns_raddir, &
	1403	gamma_mg, &
	1404	grid_level_count, &
	1405	maximum_grid_level, &
	1406	mg_switch_to_pe0_level, &
	1407	mg_switch_to_pe0, &
	1408	ngsrb
[2586]	1409
	1410
[4649]	1411	USE indices, &
	1412	ONLY: mg_loc_ind, &
	1413	nxl, &
	1414	nxl_mg, &
	1415	nxr, &
	1416	nxr_mg, &
	1417	nys, &
	1418	nys_mg, &
	1419	nyn, &
	1420	nyn_mg, &
	1421	nzb, &
	1422	nzt, &
	1423	nzt_mg
[2586]	1424
[4649]	1425	USE kinds
[2586]	1426
[4649]	1427	USE pegrid
[2586]	1428
[4649]	1429	IMPLICIT NONE
[2586]	1430
[4649]	1431	INTEGER(iwp) :: i !<
	1432	INTEGER(iwp) :: j !<
	1433	INTEGER(iwp) :: k !<
	1434	INTEGER(iwp) :: nxl_mg_save !<
	1435	INTEGER(iwp) :: nxr_mg_save !<
	1436	INTEGER(iwp) :: nyn_mg_save !<
	1437	INTEGER(iwp) :: nys_mg_save !<
	1438	INTEGER(iwp) :: nzt_mg_save !<
[2586]	1439
[4649]	1440	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1441	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
	1442	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1443	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
	1444	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1445	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !<
	1446	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1447	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
[2586]	1448
[4649]	1449	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
	1450	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !<
	1451	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
	1452	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[2586]	1453
[4649]	1454	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !<
[4429]	1455
	1456	#if defined( __parallel )
[4649]	1457	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p2_sub !<
[4429]	1458	#endif
[2586]	1459
	1460	!
[4649]	1461	!-- Restriction to the coarsest grid
	1462	10 IF ( grid_level == 1 ) THEN
[2586]	1463
	1464	!
[4649]	1465	!-- Solution on the coarsest grid. Double the number of Gauss-Seidel iterations in order to get a
	1466	!-- more accurate solution.
	1467	ngsrb = 2 * ngsrb
[2586]	1468
[4649]	1469	CALL redblack_noopt( f_mg, p_mg )
[2586]	1470
[4649]	1471	ngsrb = ngsrb / 2
[2586]	1472
	1473
[4649]	1474	ELSEIF ( grid_level /= 1 ) THEN
[2586]	1475
[4649]	1476	grid_level_count(grid_level) = grid_level_count(grid_level) + 1
[2586]	1477
	1478	!
[4649]	1479	!-- Solution on the actual grid level
	1480	CALL redblack_noopt( f_mg, p_mg )
[2586]	1481
	1482	!
[4649]	1483	!-- Determination of the actual residual
	1484	CALL resid_noopt( f_mg, p_mg, r )
[2586]	1485
	1486	!
[4649]	1487	!-- Restriction of the residual (finer grid values!) to the next coarser grid. Therefore, the
	1488	!-- grid level has to be decremented now. nxl..nzt have to be set to the coarse grid values,
	1489	!-- because these variables are needed for the exchange of ghost points in routine exchange_horiz
	1490	grid_level = grid_level - 1
	1491	nxl = nxl_mg(grid_level)
	1492	nys = nys_mg(grid_level)
	1493	nxr = nxr_mg(grid_level)
	1494	nyn = nyn_mg(grid_level)
	1495	nzt = nzt_mg(grid_level)
[2586]	1496
[4649]	1497	IF ( grid_level == mg_switch_to_pe0_level ) THEN
[2586]	1498
	1499	!
[4649]	1500	!-- From this level on, calculations are done on PE0 only. First, carry out restriction on the
	1501	!-- subdomain. Therefore, indices of the level have to be changed to subdomain values in
	1502	!-- between (otherwise, the restrict routine would expect the gathered array)
[2586]	1503
[4649]	1504	nxl_mg_save = nxl_mg(grid_level)
	1505	nxr_mg_save = nxr_mg(grid_level)
	1506	nys_mg_save = nys_mg(grid_level)
	1507	nyn_mg_save = nyn_mg(grid_level)
	1508	nzt_mg_save = nzt_mg(grid_level)
	1509	nxl_mg(grid_level) = mg_loc_ind(1,myid)
	1510	nxr_mg(grid_level) = mg_loc_ind(2,myid)
	1511	nys_mg(grid_level) = mg_loc_ind(3,myid)
	1512	nyn_mg(grid_level) = mg_loc_ind(4,myid)
	1513	nzt_mg(grid_level) = mg_loc_ind(5,myid)
	1514	nxl = mg_loc_ind(1,myid)
	1515	nxr = mg_loc_ind(2,myid)
	1516	nys = mg_loc_ind(3,myid)
	1517	nyn = mg_loc_ind(4,myid)
	1518	nzt = mg_loc_ind(5,myid)
[2586]	1519
[4649]	1520	ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1521	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
[2586]	1522
[4649]	1523	CALL restrict_noopt( f2_sub, r )
[2586]	1524
	1525	!
[4649]	1526	!-- Restore the correct indices of this level
	1527	nxl_mg(grid_level) = nxl_mg_save
	1528	nxr_mg(grid_level) = nxr_mg_save
	1529	nys_mg(grid_level) = nys_mg_save
	1530	nyn_mg(grid_level) = nyn_mg_save
	1531	nzt_mg(grid_level) = nzt_mg_save
	1532	nxl = nxl_mg(grid_level)
	1533	nxr = nxr_mg(grid_level)
	1534	nys = nys_mg(grid_level)
	1535	nyn = nyn_mg(grid_level)
	1536	nzt = nzt_mg(grid_level)
[2586]	1537	!
[4649]	1538	!-- Gather all arrays from the subdomains on PE0
[4429]	1539	#if defined( __parallel )
[4649]	1540	CALL mg_gather_noopt( f2, f2_sub )
[4429]	1541	#endif
[2586]	1542
	1543	!
[4649]	1544	!-- Set switch for routine exchange_horiz, that no ghostpoint exchange has to be carried out
	1545	!-- from now on
	1546	mg_switch_to_pe0 = .TRUE.
[2586]	1547
	1548	!
[4649]	1549	!-- In case of non-cyclic lateral boundary conditions, both in- and outflow conditions have to
	1550	!-- be used on all PEs after the switch, because then they have the total domain.
	1551	IF ( bc_lr_dirrad ) THEN
	1552	bc_dirichlet_l = .TRUE.
	1553	bc_dirichlet_r = .FALSE.
	1554	bc_radiation_l = .FALSE.
	1555	bc_radiation_r = .TRUE.
	1556	ELSEIF ( bc_lr_raddir ) THEN
	1557	bc_dirichlet_l = .FALSE.
	1558	bc_dirichlet_r = .TRUE.
	1559	bc_radiation_l = .TRUE.
	1560	bc_radiation_r = .FALSE.
	1561	ELSEIF ( child_domain .OR. nesting_offline ) THEN
	1562	bc_dirichlet_l = .TRUE.
	1563	bc_dirichlet_r = .TRUE.
	1564	ENDIF
[2586]	1565
[4649]	1566	IF ( bc_ns_dirrad ) THEN
	1567	bc_dirichlet_n = .TRUE.
	1568	bc_dirichlet_s = .FALSE.
	1569	bc_radiation_n = .FALSE.
	1570	bc_radiation_s = .TRUE.
	1571	ELSEIF ( bc_ns_raddir ) THEN
	1572	bc_dirichlet_n = .FALSE.
	1573	bc_dirichlet_s = .TRUE.
	1574	bc_radiation_n = .TRUE.
	1575	bc_radiation_s = .FALSE.
	1576	ELSEIF ( child_domain .OR. nesting_offline ) THEN
	1577	bc_dirichlet_s = .TRUE.
	1578	bc_dirichlet_n = .TRUE.
	1579	ENDIF
[2586]	1580
[4649]	1581	DEALLOCATE( f2_sub )
[2586]	1582
[4649]	1583	ELSE
[2586]	1584
[4649]	1585	CALL restrict_noopt( f2, r )
[2586]	1586
[4649]	1587	ENDIF
[2586]	1588
[4649]	1589	p2 = 0.0_wp
[2586]	1590
	1591	!
[4649]	1592	!-- Repeat the same procedure till the coarsest grid is reached
	1593	CALL next_mg_level_noopt( f2, p2, p3, r )
[2586]	1594
[4649]	1595	ENDIF
[2586]	1596
	1597	!
[4649]	1598	!-- Now follows the prolongation
	1599	IF ( grid_level >= 2 ) THEN
[2586]	1600
	1601	!
[4649]	1602	!-- Prolongation of the new residual. The values are transferred from the coarse to the next
	1603	!-- finer grid.
	1604	IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN
[2586]	1605
	1606	#if defined( __parallel )
	1607	!
[4649]	1608	!-- At this level, the new residual first has to be scattered from PE0 to the other PEs
	1609	ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,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) )
[2586]	1611
[4649]	1612	CALL mg_scatter_noopt( p2, p2_sub )
[2586]	1613
	1614	!
[4649]	1615	!-- Therefore, indices of the previous level have to be changed to subdomain values in between
	1616	!-- (otherwise, the prolong routine would expect the gathered array)
	1617	nxl_mg_save = nxl_mg(grid_level-1)
	1618	nxr_mg_save = nxr_mg(grid_level-1)
	1619	nys_mg_save = nys_mg(grid_level-1)
	1620	nyn_mg_save = nyn_mg(grid_level-1)
	1621	nzt_mg_save = nzt_mg(grid_level-1)
	1622	nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
	1623	nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
	1624	nys_mg(grid_level-1) = mg_loc_ind(3,myid)
	1625	nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
	1626	nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
[2586]	1627
	1628	!
[4649]	1629	!-- Set switch for routine exchange_horiz, that ghostpoint exchange has to be carried again
	1630	!-- out from now on.
	1631	mg_switch_to_pe0 = .FALSE.
[2586]	1632
	1633	!
[4649]	1634	!-- For non-cyclic lateral boundary conditions and in case of nesting, restore the in-/outflow
	1635	!-- conditions.
	1636	bc_dirichlet_l = .FALSE.; bc_dirichlet_r = .FALSE.
	1637	bc_dirichlet_n = .FALSE.; bc_dirichlet_s = .FALSE.
	1638	bc_radiation_l = .FALSE.; bc_radiation_r = .FALSE.
	1639	bc_radiation_n = .FALSE.; bc_radiation_s = .FALSE.
[2586]	1640
[4649]	1641	IF ( pleft == MPI_PROC_NULL ) THEN
	1642	IF ( bc_lr_dirrad .OR. child_domain .OR. nesting_offline ) THEN
	1643	bc_dirichlet_l = .TRUE.
	1644	ELSEIF ( bc_lr_raddir ) THEN
	1645	bc_radiation_l = .TRUE.
[2586]	1646	ENDIF
[4649]	1647	ENDIF
[2586]	1648
[4649]	1649	IF ( pright == MPI_PROC_NULL ) THEN
	1650	IF ( bc_lr_dirrad ) THEN
	1651	bc_radiation_r = .TRUE.
	1652	ELSEIF ( bc_lr_raddir .OR. child_domain .OR. nesting_offline ) THEN
	1653	bc_dirichlet_r = .TRUE.
[2586]	1654	ENDIF
[4649]	1655	ENDIF
[2586]	1656
[4649]	1657	IF ( psouth == MPI_PROC_NULL ) THEN
	1658	IF ( bc_ns_dirrad ) THEN
	1659	bc_radiation_s = .TRUE.
	1660	ELSEIF ( bc_ns_raddir .OR. child_domain .OR. nesting_offline ) THEN
	1661	bc_dirichlet_s = .TRUE.
[2586]	1662	ENDIF
[4649]	1663	ENDIF
[2586]	1664
[4649]	1665	IF ( pnorth == MPI_PROC_NULL ) THEN
	1666	IF ( bc_ns_dirrad .OR. child_domain .OR. nesting_offline ) THEN
	1667	bc_dirichlet_n = .TRUE.
	1668	ELSEIF ( bc_ns_raddir ) THEN
	1669	bc_radiation_n = .TRUE.
[2586]	1670	ENDIF
[4649]	1671	ENDIF
[2586]	1672
[4649]	1673	CALL prolong_noopt( p2_sub, p3 )
[2586]	1674
	1675	!
[4649]	1676	!-- Restore the correct indices of the previous level
	1677	nxl_mg(grid_level-1) = nxl_mg_save
	1678	nxr_mg(grid_level-1) = nxr_mg_save
	1679	nys_mg(grid_level-1) = nys_mg_save
	1680	nyn_mg(grid_level-1) = nyn_mg_save
	1681	nzt_mg(grid_level-1) = nzt_mg_save
[2586]	1682
[4649]	1683	DEALLOCATE( p2_sub )
[2586]	1684	#endif
	1685
[4649]	1686	ELSE
[2586]	1687
[4649]	1688	CALL prolong_noopt( p2, p3 )
[2586]	1689
[4649]	1690	ENDIF
[2586]	1691
	1692	!
[4649]	1693	!-- Computation of the new pressure correction. Therefore, values from prior grids are added up
	1694	!-- automatically stage by stage.
	1695	DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
	1696	DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
	1697	DO k = nzb, nzt_mg(grid_level)+1
	1698	p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
[2586]	1699	ENDDO
	1700	ENDDO
[4649]	1701	ENDDO
[2586]	1702
	1703	!
[4649]	1704	!-- Relaxation of the new solution
	1705	CALL redblack_noopt( f_mg, p_mg )
[2586]	1706
[4649]	1707	ENDIF
[2586]	1708
	1709
	1710	!
[4649]	1711	!-- The following few lines serve the steering of the multigrid scheme
	1712	IF ( grid_level == maximum_grid_level ) THEN
[2586]	1713
[4649]	1714	GOTO 20
[2586]	1715
[4649]	1716	ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. &
	1717	grid_level_count(grid_level) /= gamma_mg ) THEN
[2586]	1718
[4649]	1719	GOTO 10
[2586]	1720
[4649]	1721	ENDIF
[2586]	1722
	1723	!
[4649]	1724	!-- Reset counter for the next call of poismg_noopt
	1725	grid_level_count(grid_level) = 0
[2586]	1726
	1727	!
[4649]	1728	!-- Continue with the next finer level. nxl..nzt have to be set to the finer grid values, because
	1729	!-- these variables are needed for the exchange of ghost points in routine exchange_horiz.
	1730	grid_level = grid_level + 1
	1731	nxl = nxl_mg(grid_level)
	1732	nxr = nxr_mg(grid_level)
	1733	nys = nys_mg(grid_level)
	1734	nyn = nyn_mg(grid_level)
	1735	nzt = nzt_mg(grid_level)
[2586]	1736
[4649]	1737	20 CONTINUE
[2586]	1738
[4649]	1739	END SUBROUTINE next_mg_level_noopt
[2586]	1740
	1741
	1742
[4649]	1743	!--------------------------------------------------------------------------------------------------!
	1744	! Description:
	1745	! ------------
	1746	!> @TODO: Missing subroutine description
	1747	!--------------------------------------------------------------------------------------------------!
	1748	SUBROUTINE poismg_noopt_init
[2586]	1749
[4649]	1750	USE control_parameters, &
	1751	ONLY: bc_lr_cyc, &
	1752	bc_ns_cyc, &
	1753	masking_method, &
	1754	maximum_grid_level, &
	1755	psolver
[2586]	1756
[4649]	1757	USE indices, &
	1758	ONLY: flags, &
	1759	nxl_mg, &
	1760	nxr_mg, &
	1761	nyn_mg, &
	1762	nys_mg, &
	1763	nzb, &
	1764	nzt_mg, &
	1765	wall_flags_total_0, &
	1766	wall_flags_1, &
	1767	wall_flags_10, &
	1768	wall_flags_2, &
	1769	wall_flags_3, &
	1770	wall_flags_4, &
	1771	wall_flags_5, &
	1772	wall_flags_6, &
	1773	wall_flags_7, &
	1774	wall_flags_8, &
	1775	wall_flags_9
[2586]	1776
[4649]	1777	IMPLICIT NONE
[2586]	1778
[4649]	1779	INTEGER(iwp) :: i !< index variable along x
	1780	INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level
	1781	INTEGER(iwp) :: j !< index variable along y
	1782	INTEGER(iwp) :: k !< index variable along z
	1783	INTEGER(iwp) :: l !< loop variable indication current grid level
	1784	INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level
	1785	INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level
	1786	INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level
	1787	INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level
	1788	INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level
[2586]	1789
[4649]	1790	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo_tmp !<
[4414]	1791
[4649]	1792	IF ( psolver /= 'multigrid_noopt' ) RETURN
[2586]	1793	!
[4649]	1794	!-- Grid point increment of the current level.
	1795	inc = 1
	1796	DO l = maximum_grid_level, 1 , -1
[2586]	1797	!
[4649]	1798	!-- Set grid_level as it is required for exchange_horiz_2d_int
	1799	grid_level = l
[2586]	1800
[4649]	1801	nxl_l = nxl_mg(l)
	1802	nxr_l = nxr_mg(l)
	1803	nys_l = nys_mg(l)
	1804	nyn_l = nyn_mg(l)
	1805	nzt_l = nzt_mg(l)
[2586]	1806	!
[4649]	1807	!-- Assign the flag level to be calculated
	1808	SELECT CASE ( l )
	1809	CASE ( 1 )
	1810	flags => wall_flags_1
	1811	CASE ( 2 )
	1812	flags => wall_flags_2
	1813	CASE ( 3 )
	1814	flags => wall_flags_3
	1815	CASE ( 4 )
	1816	flags => wall_flags_4
	1817	CASE ( 5 )
	1818	flags => wall_flags_5
	1819	CASE ( 6 )
	1820	flags => wall_flags_6
	1821	CASE ( 7 )
	1822	flags => wall_flags_7
	1823	CASE ( 8 )
	1824	flags => wall_flags_8
	1825	CASE ( 9 )
	1826	flags => wall_flags_9
	1827	CASE ( 10 )
	1828	flags => wall_flags_10
	1829	END SELECT
[2586]	1830
	1831	!
[4649]	1832	!-- Depending on the grid level, set the respective bits in case of neighbouring walls
	1833	!-- Bit 0: wall to the bottom
	1834	!-- Bit 1: wall to the top (not realized in remaining PALM code so far)
	1835	!-- Bit 2: wall to the south
	1836	!-- Bit 3: wall to the north
	1837	!-- Bit 4: wall to the left
	1838	!-- Bit 5: wall to the right
	1839	!-- Bit 6: inside building
[2586]	1840
[4649]	1841	flags = 0
[2586]	1842	!
[4649]	1843	!-- In case of masking method, flags are not set and multigrid method works like FFT-solver
	1844	IF ( .NOT. masking_method ) THEN
[2586]	1845
	1846	!
[4649]	1847	!-- Allocate temporary array for topography heights on coarser grid level. Please note, 2
	1848	!-- ghost points are required, in order to calculate flags() on the interior ghost point.
	1849	ALLOCATE( topo_tmp(nzb:nzt_l+1,nys_l-1:nyn_l+1,nxl_l-1:nxr_l+1) )
	1850	topo_tmp = 0
	1851
	1852	DO i = nxl_l, nxr_l
	1853	DO j = nys_l, nyn_l
	1854	DO k = nzb, nzt_l
	1855	topo_tmp(k,j,i) = wall_flags_total_0(kinc,jinc,i*inc)
[2586]	1856	ENDDO
	1857	ENDDO
[4649]	1858	ENDDO
	1859	topo_tmp(nzt_l+1,:,:) = topo_tmp(nzt_l,:,:)
[2586]	1860	!
[4649]	1861	!-- Exchange ghost points on respective multigrid level. 2 ghost points are required, in order
	1862	!-- to calculate flags on.
	1863	!-- nys_l-1 / nyn_l+1 / nxl_l-1 / nxr_l+1.
	1864	CALL exchange_horiz_int( topo_tmp, nys_l, nyn_l, nxl_l, nxr_l, nzt_l, 1 )
[2586]	1865	!
[4649]	1866	!-- Set non-cyclic boundary conditions on respective multigrid level
	1867	IF ( .NOT. bc_ns_cyc ) THEN
	1868	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
	1869	! topo_tmp(:,-2,:) = topo_tmp(:,0,:)
	1870	topo_tmp(:,-1,:) = topo_tmp(:,0,:)
[2586]	1871	ENDIF
[4649]	1872	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
	1873	! topo_tmp(:,nyn_l+2,:) = topo_tmp(:,nyn_l,:)
	1874	topo_tmp(:,nyn_l+1,:) = topo_tmp(:,nyn_l,:)
[2586]	1875	ENDIF
[4649]	1876	ENDIF
	1877	IF ( .NOT. bc_lr_cyc ) THEN
	1878	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
	1879	! topo_tmp(:,:,-2) = topo_tmp(:,:,0)
	1880	topo_tmp(:,:,-1) = topo_tmp(:,:,0)
	1881	ENDIF
	1882	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
	1883	! topo_tmp(:,:,nxr_l+2) = topo_tmp(:,:,nxr_l)
	1884	topo_tmp(:,:,nxr_l+1) = topo_tmp(:,:,nxr_l)
	1885	ENDIF
	1886	ENDIF
	1887
	1888	DO i = nxl_l, nxr_l
	1889	DO j = nys_l, nyn_l
	1890	DO k = nzb, nzt_l
[2586]	1891	!
[4649]	1892	!-- Inside/outside building (inside building does not need further tests for walls)
	1893	IF ( .NOT. BTEST( topo_tmp(k,j,i), 0 ) ) THEN
[2586]	1894
[4649]	1895	flags(k,j,i) = IBSET( flags(k,j,i), 6 )
[2586]	1896
[4649]	1897	ELSE
[2586]	1898	!
[4649]	1899	!-- Bottom wall
	1900	IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) ) THEN
	1901	flags(k,j,i) = IBSET( flags(k,j,i), 0 )
	1902	ENDIF
[2586]	1903	!
[4649]	1904	!-- South wall
	1905	IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) ) THEN
	1906	flags(k,j,i) = IBSET( flags(k,j,i), 2 )
	1907	ENDIF
[2586]	1908	!
[4649]	1909	!-- North wall
	1910	IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) ) THEN
	1911	flags(k,j,i) = IBSET( flags(k,j,i), 3 )
	1912	ENDIF
[2586]	1913	!
[4649]	1914	!-- Left wall
	1915	IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) ) THEN
	1916	flags(k,j,i) = IBSET( flags(k,j,i), 4 )
	1917	ENDIF
[2586]	1918	!
[4649]	1919	!-- Right wall
	1920	IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) ) THEN
	1921	flags(k,j,i) = IBSET( flags(k,j,i), 5 )
	1922	ENDIF
[2586]	1923	!
[4649]	1924	!-- Top wall
	1925	IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) ) THEN
	1926	flags(k,j,i) = IBSET( flags(k,j,i), 7 )
	1927	ENDIF
[2586]	1928
[4649]	1929	ENDIF
	1930
[2586]	1931	ENDDO
	1932	ENDDO
[4649]	1933	ENDDO
	1934	flags(nzt_l+1,:,:) = flags(nzt_l,:,:)
[2586]	1935
[4649]	1936	CALL exchange_horiz_int( flags, nys_l, nyn_l, nxl_l, nxr_l, nzt_l, 1 )
[2586]	1937
[4649]	1938	DEALLOCATE( topo_tmp )
[2586]	1939
[4649]	1940	ENDIF
[2586]	1941
[4649]	1942	inc = inc * 2
[2586]	1943
[4649]	1944	ENDDO
[2586]	1945	!
[4649]	1946	!-- Reset grid_level to "normal" grid
	1947	grid_level = 0
[2586]	1948
[4649]	1949	END SUBROUTINE poismg_noopt_init
[2586]	1950
[4649]	1951	END MODULE poismg_noopt_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |