Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |