Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |