Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |