Home

Context Navigation

source: palm/trunk/SOURCE/poismg_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: 87.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |