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

Last change on this file since 16 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
286END 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(2*k-1,2*j,2*i) = p(k,j,i) 
330!
331!--          Points between two coarse-grid points
332             temp(2*k-1,2*j,2*i+1) = 0.5 * ( p(k,j,i) + p(k,j,i+1) )
333             temp(2*k-1,2*j+1,2*i) = 0.5 * ( p(k,j,i) + p(k,j+1,i) )
334             temp(2*k,2*j,2*i)     = 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(2*k-1,2*j+1,2*i+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(2*k,2*j,2*i+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(2*k,2*j+1,2*i)     = 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(2*k,2*j+1,2*i+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.