Home

Context Navigation

poismg.f90 @ 108

Last change on this file since 108 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 39.5 KB

Line
1	SUBROUTINE poismg( r )
2
3	!------------------------------------------------------------------------------!
4	! Attention: Loop unrolling and cache optimization in SOR-Red/Black method
5	! still does not bring the expected speedup on ibm! Further work
6	! is required.
7	!
8	! Actual revisions:
9	! -----------------
10	!
11	!
12	! Former revisions:
13	! -----------------
14	! $Id: poismg.f90 77 2007-03-29 04:26:56Z letzel $
15	!
16	! 75 2007-03-22 09:54:05Z raasch
17	! 2nd+3rd argument removed from exchange horiz
18	!
19	! RCS Log replace by Id keyword, revision history cleaned up
20	!
21	! Revision 1.6 2005/03/26 20:55:54 raasch
22	! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
23	! routine prolong simplified (one call of exchange_horiz spared)
24	!
25	! Revision 1.1 2001/07/20 13:10:51 raasch
26	! Initial revision
27	!
28	!
29	! Description:
30	! ------------
31	! Solves the Poisson equation for the perturbation pressure with a multigrid
32	! V- or W-Cycle scheme.
33	!
34	! This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
35	! September 2000 - July 2001.
36	!------------------------------------------------------------------------------!
37
38	USE arrays_3d
39	USE control_parameters
40	USE cpulog
41	USE grid_variables
42	USE indices
43	USE interfaces
44	USE pegrid
45
46	IMPLICIT NONE
47
48	REAL :: maxerror, maximum_mgcycles, residual_norm
49
50	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r
51
52	REAL, DIMENSION(:,:,:), ALLOCATABLE :: p3
53
54
55	CALL cpu_log( log_point_s(29), 'poismg', 'start' )
56
57
58	!
59	!-- Initialize arrays and variables used in this subroutine
60	ALLOCATE ( p3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
61
62
63	!
64	!-- Some boundaries have to be added to divergence array
65	CALL exchange_horiz( d )
66	d(nzb,:,:) = d(nzb+1,:,:)
67
68	!
69	!-- Initiation of the multigrid scheme. Does n cycles until the
70	!-- residual is smaller than the given limit. The accuracy of the solution
71	!-- of the poisson equation will increase with the number of cycles.
72	!-- If the number of cycles is preset by the user, this number will be
73	!-- carried out regardless of the accuracy.
74	grid_level_count = 0
75	mgcycles = 0
76	IF ( mg_cycles == -1 ) THEN
77	maximum_mgcycles = 0
78	residual_norm = 1.0
79	ELSE
80	maximum_mgcycles = mg_cycles
81	residual_norm = 0.0
82	ENDIF
83
84	DO WHILE ( residual_norm > residual_limit .OR. &
85	mgcycles < maximum_mgcycles )
86
87	CALL next_mg_level( d, p, p3, r)
88
89	!
90	!-- Calculate the residual if the user has not preset the number of
91	!-- cycles to be performed
92	IF ( maximum_mgcycles == 0 ) THEN
93	CALL resid( d, p, r )
94	maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
95	#if defined( __parallel )
96	CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
97	comm2d, ierr)
98	#else
99	residual_norm = maxerror
100	#endif
101	residual_norm = SQRT( residual_norm )
102	ENDIF
103
104	mgcycles = mgcycles + 1
105
106	!
107	!-- If the user has not limited the number of cycles, stop the run in case
108	!-- of insufficient convergence
109	IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN
110	IF ( myid == 0 ) THEN
111	PRINT*, '+++ poismg: no sufficient convergence within 1000 cycles'
112	ENDIF
113	CALL local_stop
114	ENDIF
115
116	ENDDO
117
118	DEALLOCATE( p3 )
119
120	CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
121
122	END SUBROUTINE poismg
123
124
125
126	SUBROUTINE resid( f_mg, p_mg, r )
127
128	!------------------------------------------------------------------------------!
129	! Description:
130	! ------------
131	! Computes the residual of the perturbation pressure.
132	!------------------------------------------------------------------------------!
133
134	USE arrays_3d
135	USE control_parameters
136	USE grid_variables
137	USE indices
138	USE pegrid
139
140	IMPLICIT NONE
141
142	INTEGER :: i, j, k, l
143
144	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
145	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
146	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg, r
147
148	!
149	!-- Calculate the residual
150	l = grid_level
151
152	!$OMP PARALLEL PRIVATE (i,j,k)
153	!$OMP DO
154	DO i = nxl_mg(l), nxr_mg(l)
155	DO j = nys_mg(l), nyn_mg(l)
156	DO k = nzb+1, nzt_mg(l)
157	r(k,j,i) = f_mg(k,j,i) &
158	- ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
159	- ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
160	- f2_mg(k,l) * p_mg(k+1,j,i) &
161	- f3_mg(k,l) * p_mg(k-1,j,i) &
162	+ f1_mg(k,l) * p_mg(k,j,i)
163	ENDDO
164	ENDDO
165	ENDDO
166	!$OMP END PARALLEL
167
168	!
169	!-- Horizontal boundary conditions
170	CALL exchange_horiz( r )
171
172	IF ( bc_lr /= 'cyclic' ) THEN
173	IF ( inflow_l .OR. outflow_l ) r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
174	IF ( inflow_r .OR. outflow_r ) r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
175	ENDIF
176
177	IF ( bc_ns /= 'cyclic' ) THEN
178	IF ( inflow_n .OR. outflow_n ) r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
179	IF ( inflow_s .OR. outflow_s ) r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
180	ENDIF
181
182	!
183	!-- Bottom and top boundary conditions
184	IF ( ibc_p_b == 1 ) THEN
185	r(nzb,:,: ) = r(nzb+1,:,:)
186	ELSE
187	r(nzb,:,: ) = 0.0
188	ENDIF
189
190	IF ( ibc_p_t == 1 ) THEN
191	r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
192	ELSE
193	r(nzt_mg(l)+1,:,: ) = 0.0
194	ENDIF
195
196
197	END SUBROUTINE resid
198
199
200
201	SUBROUTINE restrict( f_mg, r )
202
203	!------------------------------------------------------------------------------!
204	! Description:
205	! ------------
206	! Interpolates the residual on the next coarser grid with "full weighting"
207	! scheme
208	!------------------------------------------------------------------------------!
209
210	USE control_parameters
211	USE grid_variables
212	USE indices
213	USE pegrid
214
215	IMPLICIT NONE
216
217	INTEGER :: i, ic, j, jc, k, kc, l
218
219	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
220	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
221	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg
222
223	REAL, DIMENSION(nzb:nzt_mg(grid_level+1)+1, &
224	nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, &
225	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r
226
227	!
228	!-- Interpolate the residual
229	l = grid_level
230
231	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
232	!$OMP DO
233	DO ic = nxl_mg(l), nxr_mg(l)
234	i = 2*ic
235	DO jc = nys_mg(l), nyn_mg(l)
236	j = 2*jc
237	DO kc = nzb+1, nzt_mg(l)
238	k = 2*kc-1
239	f_mg(kc,jc,ic) = 1.0 / 64.0 * ( &
240	8.0 * r(k,j,i) &
241	+ 4.0 * ( r(k,j,i-1) + r(k,j,i+1) + &
242	r(k,j+1,i) + r(k,j-1,i) ) &
243	+ 2.0 * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
244	r(k,j-1,i+1) + r(k,j+1,i+1) ) &
245	+ 4.0 * r(k-1,j,i) &
246	+ 2.0 * ( r(k-1,j,i-1) + r(k-1,j,i+1) + &
247	r(k-1,j+1,i) + r(k-1,j-1,i) ) &
248	+ ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
249	r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
250	+ 4.0 * r(k+1,j,i) &
251	+ 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
252	r(k+1,j+1,i) + r(k+1,j-1,i) ) &
253	+ ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
254	r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
255	)
256	ENDDO
257	ENDDO
258	ENDDO
259	!$OMP END PARALLEL
260
261	!
262	!-- Horizontal boundary conditions
263	CALL exchange_horiz( f_mg )
264
265	IF ( bc_lr /= 'cyclic' ) THEN
266	IF (inflow_l .OR. outflow_l) f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
267	IF (inflow_r .OR. outflow_r) f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
268	ENDIF
269
270	IF ( bc_ns /= 'cyclic' ) THEN
271	IF (inflow_n .OR. outflow_n) f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
272	IF (inflow_s .OR. outflow_s) f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
273	ENDIF
274
275	!
276	!-- Bottom and top boundary conditions
277	IF ( ibc_p_b == 1 ) THEN
278	f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
279	ELSE
280	f_mg(nzb,:,: ) = 0.0
281	ENDIF
282
283	IF ( ibc_p_t == 1 ) THEN
284	f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
285	ELSE
286	f_mg(nzt_mg(l)+1,:,: ) = 0.0
287	ENDIF
288
289
290	END SUBROUTINE restrict
291
292
293
294	SUBROUTINE prolong( p, temp )
295
296	!------------------------------------------------------------------------------!
297	! Description:
298	! ------------
299	! Interpolates the correction of the perturbation pressure
300	! to the next finer grid.
301	!------------------------------------------------------------------------------!
302
303	USE control_parameters
304	USE pegrid
305	USE indices
306
307	IMPLICIT NONE
308
309	INTEGER :: i, j, k, l
310
311	REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
312	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
313	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p
314
315	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
316	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
317	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp
318
319
320	!
321	!-- First, store elements of the coarser grid on the next finer grid
322	l = grid_level
323
324	!$OMP PARALLEL PRIVATE (i,j,k)
325	!$OMP DO
326	DO i = nxl_mg(l-1), nxr_mg(l-1)
327	DO j = nys_mg(l-1), nyn_mg(l-1)
328	!CDIR NODEP
329	DO k = nzb+1, nzt_mg(l-1)
330	!
331	!-- Points of the coarse grid are directly stored on the next finer
332	!-- grid
333	temp(2k-1,2j,2*i) = p(k,j,i)
334	!
335	!-- Points between two coarse-grid points
336	temp(2k-1,2j,2i+1) = 0.5 ( p(k,j,i) + p(k,j,i+1) )
337	temp(2k-1,2j+1,2i) = 0.5 ( p(k,j,i) + p(k,j+1,i) )
338	temp(2k,2j,2i) = 0.5 ( p(k,j,i) + p(k+1,j,i) )
339	!
340	!-- Points in the center of the planes stretched by four points
341	!-- of the coarse grid cube
342	temp(2k-1,2j+1,2i+1) = 0.25 ( p(k,j,i) + p(k,j,i+1) + &
343	p(k,j+1,i) + p(k,j+1,i+1) )
344	temp(2k,2j,2i+1) = 0.25 ( p(k,j,i) + p(k,j,i+1) + &
345	p(k+1,j,i) + p(k+1,j,i+1) )
346	temp(2k,2j+1,2i) = 0.25 ( p(k,j,i) + p(k,j+1,i) + &
347	p(k+1,j,i) + p(k+1,j+1,i) )
348	!
349	!-- Points in the middle of coarse grid cube
350	temp(2k,2j+1,2i+1) = 0.125 ( p(k,j,i) + p(k,j,i+1) + &
351	p(k,j+1,i) + p(k,j+1,i+1) + &
352	p(k+1,j,i) + p(k+1,j,i+1) + &
353	p(k+1,j+1,i) + p(k+1,j+1,i+1) )
354	ENDDO
355	ENDDO
356	ENDDO
357	!$OMP END PARALLEL
358
359	!
360	!-- Horizontal boundary conditions
361	CALL exchange_horiz( temp )
362
363	IF ( bc_lr /= 'cyclic' ) THEN
364	IF (inflow_l .OR. outflow_l) temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
365	IF (inflow_r .OR. outflow_r) temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
366	ENDIF
367
368	IF ( bc_ns /= 'cyclic' ) THEN
369	IF (inflow_n .OR. outflow_n) temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
370	IF (inflow_s .OR. outflow_s) temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
371	ENDIF
372
373	!
374	!-- Bottom and top boundary conditions
375	IF ( ibc_p_b == 1 ) THEN
376	temp(nzb,:,: ) = temp(nzb+1,:,:)
377	ELSE
378	temp(nzb,:,: ) = 0.0
379	ENDIF
380
381	IF ( ibc_p_t == 1 ) THEN
382	temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
383	ELSE
384	temp(nzt_mg(l)+1,:,: ) = 0.0
385	ENDIF
386
387
388	END SUBROUTINE prolong
389
390
391	SUBROUTINE redblack( f_mg, p_mg )
392
393	!------------------------------------------------------------------------------!
394	! Description:
395	! ------------
396	! Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
397	! 3D-Red-Black decomposition (GS-RB) is used.
398	!------------------------------------------------------------------------------!
399
400	USE arrays_3d
401	USE control_parameters
402	USE cpulog
403	USE grid_variables
404	USE indices
405	USE interfaces
406	USE pegrid
407
408	IMPLICIT NONE
409
410	INTEGER :: colour, i, ic, j, jc, jj, k, l, n
411
412	LOGICAL :: unroll
413
414	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
415	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
416	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg
417
418
419	l = grid_level
420
421	unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. &
422	MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
423
424	DO n = 1, ngsrb
425
426	DO colour = 1, 2
427
428	IF ( .NOT. unroll ) THEN
429	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
430
431	!
432	!-- Without unrolling of loops, no cache optimization
433	DO i = nxl_mg(l), nxr_mg(l), 2
434	DO j = nys_mg(l) + 2 - colour, nyn_mg(l), 2
435	DO k = nzb+1, nzt_mg(l), 2
436	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
437	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
438	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
439	+ f2_mg(k,l) * p_mg(k+1,j,i) &
440	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
441	)
442	ENDDO
443	ENDDO
444	ENDDO
445
446	DO i = nxl_mg(l)+1, nxr_mg(l), 2
447	DO j = nys_mg(l) + (colour-1), nyn_mg(l), 2
448	DO k = nzb+1, nzt_mg(l), 2
449	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
450	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
451	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
452	+ f2_mg(k,l) * p_mg(k+1,j,i) &
453	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
454	)
455	ENDDO
456	ENDDO
457	ENDDO
458
459	DO i = nxl_mg(l), nxr_mg(l), 2
460	DO j = nys_mg(l) + (colour-1), nyn_mg(l), 2
461	DO k = nzb+2, nzt_mg(l), 2
462	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
463	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
464	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
465	+ f2_mg(k,l) * p_mg(k+1,j,i) &
466	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
467	)
468	ENDDO
469	ENDDO
470	ENDDO
471
472	DO i = nxl_mg(l)+1, nxr_mg(l), 2
473	DO j = nys_mg(l) + 2 - colour, nyn_mg(l), 2
474	DO k = nzb+2, nzt_mg(l), 2
475	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
476	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
477	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
478	+ f2_mg(k,l) * p_mg(k+1,j,i) &
479	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
480	)
481	ENDDO
482	ENDDO
483	ENDDO
484	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
485
486	ELSE
487
488	!
489	!-- Loop unrolling along y, only one i loop for better cache use
490	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
491	DO ic = nxl_mg(l), nxr_mg(l), 2
492	DO jc = nys_mg(l), nyn_mg(l), 4
493	i = ic
494	jj = jc+2-colour
495	DO k = nzb+1, nzt_mg(l), 2
496	j = jj
497	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
498	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
499	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
500	+ f2_mg(k,l) * p_mg(k+1,j,i) &
501	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
502	)
503	j = jj+2
504	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
505	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
506	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
507	+ f2_mg(k,l) * p_mg(k+1,j,i) &
508	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
509	)
510	! j = jj+4
511	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
512	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
513	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
514	! + f2_mg(k,l) * p_mg(k+1,j,i) &
515	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
516	! )
517	! j = jj+6
518	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
519	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
520	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
521	! + f2_mg(k,l) * p_mg(k+1,j,i) &
522	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
523	! )
524	ENDDO
525
526	i = ic+1
527	jj = jc+colour-1
528	DO k = nzb+1, nzt_mg(l), 2
529	j =jj
530	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
531	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
532	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
533	+ f2_mg(k,l) * p_mg(k+1,j,i) &
534	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
535	)
536	j = jj+2
537	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
538	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
539	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
540	+ f2_mg(k,l) * p_mg(k+1,j,i) &
541	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
542	)
543	! j = jj+4
544	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
545	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
546	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
547	! + f2_mg(k,l) * p_mg(k+1,j,i) &
548	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
549	! )
550	! j = jj+6
551	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
552	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
553	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
554	! + f2_mg(k,l) * p_mg(k+1,j,i) &
555	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
556	! )
557	ENDDO
558
559	i = ic
560	jj = jc+colour-1
561	DO k = nzb+2, nzt_mg(l), 2
562	j =jj
563	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
564	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
565	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
566	+ f2_mg(k,l) * p_mg(k+1,j,i) &
567	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
568	)
569	j = jj+2
570	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
571	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
572	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
573	+ f2_mg(k,l) * p_mg(k+1,j,i) &
574	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
575	)
576	! j = jj+4
577	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
578	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
579	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
580	! + f2_mg(k,l) * p_mg(k+1,j,i) &
581	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
582	! )
583	! j = jj+6
584	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
585	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
586	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
587	! + f2_mg(k,l) * p_mg(k+1,j,i) &
588	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
589	! )
590	ENDDO
591
592	i = ic+1
593	jj = jc+2-colour
594	DO k = nzb+2, nzt_mg(l), 2
595	j =jj
596	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
597	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
598	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
599	+ f2_mg(k,l) * p_mg(k+1,j,i) &
600	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
601	)
602	j = jj+2
603	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
604	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
605	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
606	+ f2_mg(k,l) * p_mg(k+1,j,i) &
607	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
608	)
609	! j = jj+4
610	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
611	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
612	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
613	! + f2_mg(k,l) * p_mg(k+1,j,i) &
614	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
615	! )
616	! j = jj+6
617	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
618	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
619	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
620	! + f2_mg(k,l) * p_mg(k+1,j,i) &
621	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
622	! )
623	ENDDO
624
625	ENDDO
626	ENDDO
627	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
628
629	ENDIF
630
631	!
632	!-- Horizontal boundary conditions
633	CALL exchange_horiz( p_mg )
634
635	IF ( bc_lr /= 'cyclic' ) THEN
636	IF ( inflow_l .OR. outflow_l ) THEN
637	p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
638	ENDIF
639	IF ( inflow_r .OR. outflow_r ) THEN
640	p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
641	ENDIF
642	ENDIF
643
644	IF ( bc_ns /= 'cyclic' ) THEN
645	IF ( inflow_n .OR. outflow_n ) THEN
646	p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
647	ENDIF
648	IF ( inflow_s .OR. outflow_s ) THEN
649	p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
650	ENDIF
651	ENDIF
652
653	!
654	!-- Bottom and top boundary conditions
655	IF ( ibc_p_b == 1 ) THEN
656	p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
657	ELSE
658	p_mg(nzb,:,: ) = 0.0
659	ENDIF
660
661	IF ( ibc_p_t == 1 ) THEN
662	p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
663	ELSE
664	p_mg(nzt_mg(l)+1,:,: ) = 0.0
665	ENDIF
666
667	ENDDO
668
669	ENDDO
670
671
672	END SUBROUTINE redblack
673
674
675
676	SUBROUTINE mg_gather( f2, f2_sub )
677
678	USE control_parameters
679	USE cpulog
680	USE indices
681	USE interfaces
682	USE pegrid
683
684	IMPLICIT NONE
685
686	INTEGER :: n, nwords, sender
687
688	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
689	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
690	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2
691
692	REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
693	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
694	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub
695
696	!
697	!-- Find out the number of array elements of the subdomain array
698	nwords = SIZE( f2_sub )
699
700	#if defined( __parallel )
701	CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
702
703	IF ( myid == 0 ) THEN
704	!
705	!-- Store the local subdomain array on the total array
706	f2(:,mg_loc_ind(3,0)-1:mg_loc_ind(4,0)+1, &
707	mg_loc_ind(1,0)-1:mg_loc_ind(2,0)+1) = f2_sub
708
709	!
710	!-- Receive the subdomain arrays from all other PEs and store them on the
711	!-- total array
712	DO n = 1, numprocs-1
713	!
714	!-- Receive the arrays in arbitrary order from the PEs.
715	CALL MPI_RECV( f2_sub(nzb,mg_loc_ind(3,0)-1,mg_loc_ind(1,0)-1), &
716	nwords, MPI_REAL, MPI_ANY_SOURCE, 1, comm2d, status, &
717	ierr )
718	sender = status(MPI_SOURCE)
719	f2(:,mg_loc_ind(3,sender)-1:mg_loc_ind(4,sender)+1, &
720	mg_loc_ind(1,sender)-1:mg_loc_ind(2,sender)+1) = f2_sub
721	ENDDO
722
723	ELSE
724	!
725	!-- Send subdomain array to PE0
726	CALL MPI_SEND( f2_sub(nzb,mg_loc_ind(3,myid)-1,mg_loc_ind(1,myid)-1), &
727	nwords, MPI_REAL, 0, 1, comm2d, ierr )
728	ENDIF
729
730	CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
731	#endif
732
733	END SUBROUTINE mg_gather
734
735
736
737	SUBROUTINE mg_scatter( p2, p2_sub )
738	!
739	!-- TODO: It may be possible to improve the speed of this routine by using
740	!-- non-blocking communication
741
742	USE control_parameters
743	USE cpulog
744	USE indices
745	USE interfaces
746	USE pegrid
747
748	IMPLICIT NONE
749
750	INTEGER :: n, nwords, sender
751
752	REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
753	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
754	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2
755
756	REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
757	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
758	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub
759
760	!
761	!-- Find out the number of array elements of the subdomain array
762	nwords = SIZE( p2_sub )
763
764	#if defined( __parallel )
765	CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
766
767	IF ( myid == 0 ) THEN
768	!
769	!-- Scatter the subdomain arrays to the other PEs by blocking
770	!-- communication
771	DO n = 1, numprocs-1
772
773	p2_sub = p2(:,mg_loc_ind(3,n)-1:mg_loc_ind(4,n)+1, &
774	mg_loc_ind(1,n)-1:mg_loc_ind(2,n)+1)
775
776	CALL MPI_SEND( p2_sub(nzb,mg_loc_ind(3,0)-1,mg_loc_ind(1,0)-1), &
777	nwords, MPI_REAL, n, 1, comm2d, ierr )
778
779	ENDDO
780
781	!
782	!-- Store data from the total array to the local subdomain array
783	p2_sub = p2(:,mg_loc_ind(3,0)-1:mg_loc_ind(4,0)+1, &
784	mg_loc_ind(1,0)-1:mg_loc_ind(2,0)+1)
785
786	ELSE
787	!
788	!-- Receive subdomain array from PE0
789	CALL MPI_RECV( p2_sub(nzb,mg_loc_ind(3,myid)-1,mg_loc_ind(1,myid)-1), &
790	nwords, MPI_REAL, 0, 1, comm2d, status, ierr )
791
792	ENDIF
793
794	CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
795	#endif
796
797	END SUBROUTINE mg_scatter
798
799
800
801	RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
802
803	!------------------------------------------------------------------------------!
804	! Description:
805	! ------------
806	! This is where the multigrid technique takes place. V- and W- Cycle are
807	! implemented and steered by the parameter "gamma". Parameter "nue" determines
808	! the convergence of the multigrid iterative solution. There are nue times
809	! RB-GS iterations. It should be set to "1" or "2", considering the time effort
810	! one would like to invest. Last choice shows a very good converging factor,
811	! but leads to an increase in computing time.
812	!------------------------------------------------------------------------------!
813
814	USE arrays_3d
815	USE control_parameters
816	USE grid_variables
817	USE indices
818	USE pegrid
819
820	IMPLICIT NONE
821
822	INTEGER :: i, j, k, nxl_mg_save, nxr_mg_save, nyn_mg_save, nys_mg_save, &
823	nzt_mg_save
824
825	LOGICAL :: restore_boundary_lr_on_pe0, restore_boundary_ns_on_pe0
826
827	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
828	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
829	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg, p3, r
830
831	REAL, DIMENSION(:,:,:), ALLOCATABLE :: f2, f2_sub, p2, p2_sub
832
833	!
834	!-- Restriction to the coarsest grid
835	10 IF ( grid_level == 1 ) THEN
836
837	!
838	!-- Solution on the coarsest grid. Double the number of Gauss-Seidel
839	!-- iterations in order to get a more accurate solution.
840	ngsrb = 2 * ngsrb
841	CALL redblack( f_mg, p_mg )
842	ngsrb = ngsrb / 2
843
844	ELSEIF ( grid_level /= 1 ) THEN
845
846	grid_level_count(grid_level) = grid_level_count(grid_level) + 1
847
848	!
849	!-- Solution on the actual grid level
850	CALL redblack( f_mg, p_mg )
851
852	!
853	!-- Determination of the actual residual
854	CALL resid( f_mg, p_mg, r )
855
856	!
857	!-- Restriction of the residual (finer grid values!) to the next coarser
858	!-- grid. Therefore, the grid level has to be decremented now. nxl..nzt have
859	!-- to be set to the coarse grid values, because these variables are needed
860	!-- for the exchange of ghost points in routine exchange_horiz
861	grid_level = grid_level - 1
862	nxl = nxl_mg(grid_level)
863	nxr = nxr_mg(grid_level)
864	nys = nys_mg(grid_level)
865	nyn = nyn_mg(grid_level)
866	nzt = nzt_mg(grid_level)
867
868	ALLOCATE( f2(nzb:nzt_mg(grid_level)+1, &
869	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
870	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1), &
871	p2(nzb:nzt_mg(grid_level)+1, &
872	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
873	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
874
875	IF ( grid_level == mg_switch_to_pe0_level ) THEN
876	! print*, 'myid=',myid, ' restrict and switch to PE0. level=', grid_level
877	!
878	!-- From this level on, calculations are done on PE0 only.
879	!-- First, carry out restriction on the subdomain.
880	!-- Therefore, indices of the level have to be changed to subdomain values
881	!-- in between (otherwise, the restrict routine would expect
882	!-- the gathered array)
883	nxl_mg_save = nxl_mg(grid_level)
884	nxr_mg_save = nxr_mg(grid_level)
885	nys_mg_save = nys_mg(grid_level)
886	nyn_mg_save = nyn_mg(grid_level)
887	nzt_mg_save = nzt_mg(grid_level)
888	nxl_mg(grid_level) = mg_loc_ind(1,myid)
889	nxr_mg(grid_level) = mg_loc_ind(2,myid)
890	nys_mg(grid_level) = mg_loc_ind(3,myid)
891	nyn_mg(grid_level) = mg_loc_ind(4,myid)
892	nzt_mg(grid_level) = mg_loc_ind(5,myid)
893	nxl = mg_loc_ind(1,myid)
894	nxr = mg_loc_ind(2,myid)
895	nys = mg_loc_ind(3,myid)
896	nyn = mg_loc_ind(4,myid)
897	nzt = mg_loc_ind(5,myid)
898
899	ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1, &
900	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
901	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
902
903	CALL restrict( f2_sub, r )
904
905	!
906	!-- Restore the correct indices of this level
907	nxl_mg(grid_level) = nxl_mg_save
908	nxr_mg(grid_level) = nxr_mg_save
909	nys_mg(grid_level) = nys_mg_save
910	nyn_mg(grid_level) = nyn_mg_save
911	nzt_mg(grid_level) = nzt_mg_save
912	nxl = nxl_mg(grid_level)
913	nxr = nxr_mg(grid_level)
914	nys = nys_mg(grid_level)
915	nyn = nyn_mg(grid_level)
916	nzt = nzt_mg(grid_level)
917
918	!
919	!-- Gather all arrays from the subdomains on PE0
920	CALL mg_gather( f2, f2_sub )
921
922	!
923	!-- Set switch for routine exchange_horiz, that no ghostpoint exchange
924	!-- has to be carried out from now on
925	mg_switch_to_pe0 = .TRUE.
926
927	!
928	!-- In case of non-cyclic lateral boundary conditions, both in- and
929	!-- outflow conditions have to be used on PE0 after the switch, because
930	!-- it then contains the total domain. Due to the virtual processor
931	!-- grid, before the switch, PE0 can have in-/outflow at the left
932	!-- and south wall only (or on opposite walls in case of a 1d
933	!-- decomposition).
934	restore_boundary_lr_on_pe0 = .FALSE.
935	restore_boundary_ns_on_pe0 = .FALSE.
936	IF ( myid == 0 ) THEN
937	IF ( inflow_l .AND. .NOT. outflow_r ) THEN
938	outflow_r = .TRUE.
939	restore_boundary_lr_on_pe0 = .TRUE.
940	ENDIF
941	IF ( outflow_l .AND. .NOT. inflow_r ) THEN
942	inflow_r = .TRUE.
943	restore_boundary_lr_on_pe0 = .TRUE.
944	ENDIF
945	IF ( inflow_s .AND. .NOT. outflow_n ) THEN
946	outflow_n = .TRUE.
947	restore_boundary_ns_on_pe0 = .TRUE.
948	ENDIF
949	IF ( outflow_s .AND. .NOT. inflow_n ) THEN
950	inflow_n = .TRUE.
951	restore_boundary_ns_on_pe0 = .TRUE.
952	ENDIF
953	ENDIF
954
955	DEALLOCATE( f2_sub )
956
957	ELSE
958
959	CALL restrict( f2, r )
960
961	ENDIF
962	p2 = 0.0
963
964	!
965	!-- Repeat the same procedure till the coarsest grid is reached
966	IF ( myid == 0 .OR. grid_level > mg_switch_to_pe0_level ) THEN
967	CALL next_mg_level( f2, p2, p3, r )
968	ENDIF
969
970	ENDIF
971
972	!
973	!-- Now follows the prolongation
974	IF ( grid_level >= 2 ) THEN
975
976	!
977	!-- Grid level has to be incremented on the PEs where next_mg_level
978	!-- has not been called before (normally it is incremented at the end
979	!-- of next_mg_level)
980	IF ( myid /= 0 .AND. grid_level == mg_switch_to_pe0_level ) THEN
981	grid_level = grid_level + 1
982	nxl = nxl_mg(grid_level)
983	nxr = nxr_mg(grid_level)
984	nys = nys_mg(grid_level)
985	nyn = nyn_mg(grid_level)
986	nzt = nzt_mg(grid_level)
987	ENDIF
988
989	!
990	!-- Prolongation of the new residual. The values are transferred
991	!-- from the coarse to the next finer grid.
992	IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN
993	!
994	!-- At this level, the new residual first has to be scattered from
995	!-- PE0 to the other PEs
996	ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1, &
997	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
998	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
999
1000	CALL mg_scatter( p2, p2_sub )
1001
1002	!
1003	!-- Therefore, indices of the previous level have to be changed to
1004	!-- subdomain values in between (otherwise, the prolong routine would
1005	!-- expect the gathered array)
1006	nxl_mg_save = nxl_mg(grid_level-1)
1007	nxr_mg_save = nxr_mg(grid_level-1)
1008	nys_mg_save = nys_mg(grid_level-1)
1009	nyn_mg_save = nyn_mg(grid_level-1)
1010	nzt_mg_save = nzt_mg(grid_level-1)
1011	nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1012	nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1013	nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1014	nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1015	nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1016
1017	!
1018	!-- Set switch for routine exchange_horiz, that ghostpoint exchange
1019	!-- has to be carried again out from now on
1020	mg_switch_to_pe0 = .FALSE.
1021
1022	!
1023	!-- In case of non-cyclic lateral boundary conditions, restore the
1024	!-- in-/outflow conditions on PE0
1025	IF ( myid == 0 ) THEN
1026	IF ( restore_boundary_lr_on_pe0 ) THEN
1027	IF ( inflow_l ) outflow_r = .FALSE.
1028	IF ( outflow_l ) inflow_r = .FALSE.
1029	ENDIF
1030	IF ( restore_boundary_ns_on_pe0 ) THEN
1031	IF ( inflow_s ) outflow_n = .FALSE.
1032	IF ( outflow_s ) inflow_n = .FALSE.
1033	ENDIF
1034	ENDIF
1035
1036	CALL prolong( p2_sub, p3 )
1037
1038	!
1039	!-- Restore the correct indices of the previous level
1040	nxl_mg(grid_level-1) = nxl_mg_save
1041	nxr_mg(grid_level-1) = nxr_mg_save
1042	nys_mg(grid_level-1) = nys_mg_save
1043	nyn_mg(grid_level-1) = nyn_mg_save
1044	nzt_mg(grid_level-1) = nzt_mg_save
1045
1046	DEALLOCATE( p2_sub )
1047
1048	ELSE
1049
1050	CALL prolong( p2, p3 )
1051
1052	ENDIF
1053
1054	!
1055	!-- Temporary arrays for the actual grid are not needed any more
1056	DEALLOCATE( p2, f2 )
1057
1058	!
1059	!-- Computation of the new pressure correction. Therefore,
1060	!-- values from prior grids are added up automatically stage by stage.
1061	DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1062	DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1063	DO k = nzb, nzt_mg(grid_level)+1
1064	p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1065	ENDDO
1066	ENDDO
1067	ENDDO
1068
1069	!
1070	!-- Relaxation of the new solution
1071	CALL redblack( f_mg, p_mg )
1072
1073	ENDIF
1074
1075	!
1076	!-- The following few lines serve the steering of the multigrid scheme
1077	IF ( grid_level == maximum_grid_level ) THEN
1078
1079	GOTO 20
1080
1081	ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. &
1082	grid_level_count(grid_level) /= gamma_mg ) THEN
1083
1084	GOTO 10
1085
1086	ENDIF
1087
1088	!
1089	!-- Reset counter for the next call of poismg
1090	grid_level_count(grid_level) = 0
1091
1092	!
1093	!-- Continue with the next finer level. nxl..nzt have to be
1094	!-- set to the finer grid values, because these variables are needed for the
1095	!-- exchange of ghost points in routine exchange_horiz
1096	grid_level = grid_level + 1
1097	nxl = nxl_mg(grid_level)
1098	nxr = nxr_mg(grid_level)
1099	nys = nys_mg(grid_level)
1100	nyn = nyn_mg(grid_level)
1101	nzt = nzt_mg(grid_level)
1102
1103	20 CONTINUE
1104
1105	END SUBROUTINE next_mg_level

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

Context Navigation

source: palm/trunk/SOURCE/poismg.f90 @ 108

Download in other formats: