Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |