source: palm/trunk/SOURCE/init_pegrid.f90 @ 75

Last change on this file since 75 was 75, checked in by raasch, 14 years ago

preliminary update for changes concerning non-cyclic boundary conditions

  • Property svn:keywords set to Id
File size: 24.4 KB
Line 
1 SUBROUTINE init_pegrid
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! uxrp, vynp eliminated,
7! dirichlet/neumann changed to dirichlet/radiation, etc.,
8! poisfft_init is only called if fft-solver is switched on
9!
10! Former revisions:
11! -----------------
12! $Id: init_pegrid.f90 75 2007-03-22 09:54:05Z raasch $
13! RCS Log replace by Id keyword, revision history cleaned up
14!
15! Revision 1.28  2006/04/26 13:23:32  raasch
16! lcmuk does not understand the !$ comment so a cpp-directive is required
17!
18! Revision 1.1  1997/07/24 11:15:09  raasch
19! Initial revision
20!
21!
22! Description:
23! ------------
24! Determination of the virtual processor topology (if not prescribed by the
25! user)and computation of the grid point number and array bounds of the local
26! domains.
27!------------------------------------------------------------------------------!
28
29    USE control_parameters
30    USE fft_xy
31    USE indices
32    USE pegrid
33    USE poisfft_mod
34    USE poisfft_hybrid_mod
35    USE statistics
36    USE transpose_indices
37
38
39    IMPLICIT NONE
40
41    INTEGER ::  gathered_size, i, ind(5), j, k, maximum_grid_level_l,     &
42                mg_switch_to_pe0_level_l, mg_levels_x, mg_levels_y,       &
43                mg_levels_z, nnx_y, nnx_z, nny_x, nny_z, nnz_x, nnz_y,    &
44                numproc_sqr, nx_total, nxl_l, nxr_l, nyn_l, nys_l, nzb_l, &
45                nzt_l, omp_get_num_threads, subdomain_size
46
47    INTEGER, DIMENSION(:), ALLOCATABLE ::  ind_all, nxlf, nxrf, nynf, nysf
48
49    LOGICAL ::  found
50
51!
52!-- Get the number of OpenMP threads
53    !$OMP PARALLEL
54#if defined( __lcmuk )
55    threads_per_task = omp_get_num_threads()
56#else
57!$  threads_per_task = omp_get_num_threads()
58#endif
59    !$OMP END PARALLEL
60
61
62#if defined( __parallel )
63!
64!-- Determine the processor topology or check it, if prescribed by the user
65    IF ( npex == -1  .AND.  npey == -1 )  THEN
66
67!
68!--    Automatic determination of the topology
69!--    The default on SMP- and cluster-hosts is a 1d-decomposition along x
70#if defined( __lcmuk )
71       host = 'lcmuk'
72#endif
73       IF ( host(1:3) == 'ibm'  .OR.  host(1:3) == 'nec'  .OR. &
74            host(1:2) == 'lc'   .OR.  host(1:3) == 'dec' )  THEN
75
76          pdims(1) = numprocs
77          pdims(2) = 1
78
79       ELSE
80
81          numproc_sqr = SQRT( REAL( numprocs ) )
82          pdims(1)    = MAX( numproc_sqr , 1 )
83          DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
84             pdims(1) = pdims(1) - 1
85          ENDDO
86          pdims(2) = numprocs / pdims(1)
87
88       ENDIF
89
90    ELSEIF ( npex /= -1  .AND.  npey /= -1 )  THEN
91
92!
93!--    Prescribed by user. Number of processors on the prescribed topology
94!--    must be equal to the number of PEs available to the job
95       IF ( ( npex * npey ) /= numprocs )  THEN
96          PRINT*, '+++ init_pegrid:'
97          PRINT*, '    number of PEs of the prescribed topology (', npex*npey, &
98                      ') does not match the number of PEs available to the ',  &
99                      'job (', numprocs, ')'
100          CALL local_stop
101       ENDIF
102       pdims(1) = npex
103       pdims(2) = npey
104
105    ELSE
106!
107!--    If the processor topology is prescribed by the user, the number of
108!--    PEs must be given in both directions
109       PRINT*, '+++ init_pegrid:'
110       PRINT*, '    if the processor topology is prescribed by the user, ',   &
111                    'both values of "npex" and "npey" must be given in the ', &
112                    'NAMELIST-parameter file'
113       CALL local_stop
114
115    ENDIF
116
117!
118!-- The hybrid solver can only be used in case of a 1d-decomposition along x
119    IF ( pdims(2) /= 1  .AND.  psolver == 'poisfft_hybrid' )  THEN
120       IF ( myid == 0 )  THEN
121          PRINT*, '*** init_pegrid: psolver = "poisfft_hybrid" can only be'
122          PRINT*, '                 used in case of a 1d-decomposition along x'
123       ENDIF
124    ENDIF
125
126!
127!-- If necessary, set horizontal boundary conditions to non-cyclic
128    IF ( bc_lr /= 'cyclic' )  cyclic(1) = .FALSE.
129    IF ( bc_ns /= 'cyclic' )  cyclic(2) = .FALSE.
130
131!
132!-- Create the virtual processor grid
133    CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
134                          comm2d, ierr )
135    CALL MPI_COMM_RANK( comm2d, myid, ierr )
136    WRITE (myid_char,'(''_'',I4.4)')  myid
137
138    CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
139    CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
140    CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
141
142!
143!-- Determine sub-topologies for transpositions
144!-- Transposition from z to x:
145    remain_dims(1) = .TRUE.
146    remain_dims(2) = .FALSE.
147    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
148    CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
149!
150!-- Transposition from x to y
151    remain_dims(1) = .FALSE.
152    remain_dims(2) = .TRUE.
153    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
154    CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
155
156
157!
158!-- Find a grid (used for array d) which will match the transposition demands
159    IF ( grid_matching == 'strict' )  THEN
160
161       nxa = nx;  nya = ny;  nza = nz
162
163    ELSE
164
165       found = .FALSE.
166   xn: DO  nxa = nx, 2*nx
167!
168!--       Meet conditions for nx
169          IF ( MOD( nxa+1, pdims(1) ) /= 0 .OR. &
170               MOD( nxa+1, pdims(2) ) /= 0 )  CYCLE xn
171
172      yn: DO  nya = ny, 2*ny
173!
174!--          Meet conditions for ny
175             IF ( MOD( nya+1, pdims(2) ) /= 0 .OR. &
176                  MOD( nya+1, pdims(1) ) /= 0 )  CYCLE yn
177
178
179         zn: DO  nza = nz, 2*nz
180!
181!--             Meet conditions for nz
182                IF ( ( MOD( nza, pdims(1) ) /= 0  .AND.  pdims(1) /= 1  .AND. &
183                       pdims(2) /= 1 )  .OR.                                  &
184                     ( MOD( nza, pdims(2) ) /= 0  .AND.  dt_dosp /= 9999999.9 &
185                     ) )  THEN
186                   CYCLE zn
187                ELSE
188                   found = .TRUE.
189                   EXIT xn
190                ENDIF
191
192             ENDDO zn
193
194          ENDDO yn
195
196       ENDDO xn
197
198       IF ( .NOT. found )  THEN
199          IF ( myid == 0 )  THEN
200             PRINT*,'+++ init_pegrid: no matching grid for transpositions found'
201          ENDIF
202          CALL local_stop
203       ENDIF
204
205    ENDIF
206
207!
208!-- Calculate array bounds in x-direction for every PE.
209!-- The last PE along x may get less grid points than the others
210    ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
211              nysf(0:pdims(2)-1), nnx_pe(0:pdims(1)-1), nny_pe(0:pdims(2)-1) )
212
213    IF ( MOD( nxa+1 , pdims(1) ) /= 0 )  THEN
214       IF ( myid == 0 )  THEN
215          PRINT*,'+++ x-direction:  gridpoint number (',nx+1,') is not an'
216          PRINT*,'                  integral divisor of the number of proces', &
217                                   &'sors (', pdims(1),')'
218       ENDIF
219       CALL local_stop
220    ELSE
221       nnx  = ( nxa + 1 ) / pdims(1)
222       IF ( nnx*pdims(1) - ( nx + 1) > nnx )  THEN
223          IF ( myid == 0 )  THEN
224             PRINT*,'+++ x-direction: nx does not match the requirements ', &
225                         'given by the number of PEs'
226             PRINT*,'                 used'
227             PRINT*,'    please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
228                         - ( nx + 1 ) ) ), ' instead of nx =', nx
229          ENDIF
230          CALL local_stop
231       ENDIF
232    ENDIF   
233
234!
235!-- Left and right array bounds, number of gridpoints
236    DO  i = 0, pdims(1)-1
237       nxlf(i)   = i * nnx
238       nxrf(i)   = ( i + 1 ) * nnx - 1
239       nnx_pe(i) = MIN( nx, nxrf(i) ) - nxlf(i) + 1
240    ENDDO
241
242!
243!-- Calculate array bounds in y-direction for every PE.
244    IF ( MOD( nya+1 , pdims(2) ) /= 0 )  THEN
245       IF ( myid == 0 )  THEN
246          PRINT*,'+++ y-direction:  gridpoint number (',ny+1,') is not an'
247          PRINT*,'                  integral divisor of the number of proces', &
248                                   &'sors (', pdims(2),')'
249       ENDIF
250       CALL local_stop
251    ELSE
252       nny  = ( nya + 1 ) / pdims(2)
253       IF ( nny*pdims(2) - ( ny + 1) > nny )  THEN
254          IF ( myid == 0 )  THEN
255             PRINT*,'+++ x-direction: nx does not match the requirements ', &
256                         'given by the number of PEs'
257             PRINT*,'                 used'
258             PRINT*,'    please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
259                         - ( nx + 1 ) ) ), ' instead of nx =', nx
260          ENDIF
261          CALL local_stop
262       ENDIF
263    ENDIF   
264
265!
266!-- South and north array bounds
267    DO  j = 0, pdims(2)-1
268       nysf(j)   = j * nny
269       nynf(j)   = ( j + 1 ) * nny - 1
270       nny_pe(j) = MIN( ny, nynf(j) ) - nysf(j) + 1
271    ENDDO
272
273!
274!-- Local array bounds of the respective PEs
275    nxl  = nxlf(pcoord(1))
276    nxra = nxrf(pcoord(1))
277    nxr  = MIN( nx, nxra )
278    nys  = nysf(pcoord(2))
279    nyna = nynf(pcoord(2))
280    nyn  = MIN( ny, nyna )
281    nzb  = 0
282    nzta = nza
283    nzt  = MIN( nz, nzta )
284    nnz  = nza
285
286!
287!-- Calculate array bounds and gridpoint numbers for the transposed arrays
288!-- (needed in the pressure solver)
289!-- For the transposed arrays, cyclic boundaries as well as top and bottom
290!-- boundaries are omitted, because they are obstructive to the transposition
291
292!
293!-- 1. transposition  z --> x
294!-- This transposition is not neccessary in case of a 1d-decomposition along x,
295!-- except that the uptream-spline method is switched on
296    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
297         scalar_advec == 'ups-scheme' )  THEN
298
299       IF ( pdims(2) == 1  .AND. ( momentum_advec == 'ups-scheme'  .OR. &
300            scalar_advec == 'ups-scheme' ) )  THEN
301          IF ( myid == 0 )  THEN
302             PRINT*,'+++ WARNING: init_pegrid: 1d-decomposition along x ', &
303                                &'chosen but nz restrictions may occur'
304             PRINT*,'             since ups-scheme is activated'
305          ENDIF
306       ENDIF
307       nys_x  = nys
308       nyn_xa = nyna
309       nyn_x  = nyn
310       nny_x  = nny
311       IF ( MOD( nza , pdims(1) ) /= 0 )  THEN
312          IF ( myid == 0 )  THEN
313             PRINT*,'+++ transposition z --> x:'
314             PRINT*,'    nz=',nz,' is not an integral divisior of pdims(1)=', &
315                    &pdims(1)
316          ENDIF
317          CALL local_stop
318       ENDIF
319       nnz_x  = nza / pdims(1)
320       nzb_x  = 1 + myidx * nnz_x
321       nzt_xa = ( myidx + 1 ) * nnz_x
322       nzt_x  = MIN( nzt, nzt_xa )
323
324       sendrecvcount_zx = nnx * nny * nnz_x
325
326    ENDIF
327
328!
329!-- 2. transposition  x --> y
330    nnz_y  = nnz_x
331    nzb_y  = nzb_x
332    nzt_ya = nzt_xa
333    nzt_y  = nzt_x
334    IF ( MOD( nxa+1 , pdims(2) ) /= 0 )  THEN
335       IF ( myid == 0 )  THEN
336          PRINT*,'+++ transposition x --> y:'
337          PRINT*,'    nx+1=',nx+1,' is not an integral divisor of ',&
338                 &'pdims(2)=',pdims(2)
339       ENDIF
340       CALL local_stop
341    ENDIF
342    nnx_y = (nxa+1) / pdims(2)
343    nxl_y = myidy * nnx_y
344    nxr_ya = ( myidy + 1 ) * nnx_y - 1
345    nxr_y  = MIN( nx, nxr_ya )
346
347    sendrecvcount_xy = nnx_y * nny_x * nnz_y
348
349!
350!-- 3. transposition  y --> z  (ELSE:  x --> y  in case of 1D-decomposition
351!-- along x)
352    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
353         scalar_advec == 'ups-scheme' )  THEN
354!
355!--    y --> z
356!--    This transposition is not neccessary in case of a 1d-decomposition
357!--    along x, except that the uptream-spline method is switched on
358       nnx_z  = nnx_y
359       nxl_z  = nxl_y
360       nxr_za = nxr_ya
361       nxr_z  = nxr_y
362       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
363          IF ( myid == 0 )  THEN
364             PRINT*,'+++ Transposition y --> z:'
365             PRINT*,'    ny+1=',ny+1,' is not an integral divisor of ',&
366                    &'pdims(1)=',pdims(1)
367          ENDIF
368          CALL local_stop
369       ENDIF
370       nny_z  = (nya+1) / pdims(1)
371       nys_z  = myidx * nny_z
372       nyn_za = ( myidx + 1 ) * nny_z - 1
373       nyn_z  = MIN( ny, nyn_za )
374
375       sendrecvcount_yz = nnx_y * nny_z * nnz_y
376
377    ELSE
378!
379!--    x --> y. This condition must be fulfilled for a 1D-decomposition along x
380       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
381          IF ( myid == 0 )  THEN
382             PRINT*,'+++ Transposition x --> y:'
383             PRINT*,'    ny+1=',ny+1,' is not an integral divisor of ',&
384                    &'pdims(1)=',pdims(1)
385          ENDIF
386          CALL local_stop
387       ENDIF
388
389    ENDIF
390
391!
392!-- Indices for direct transpositions z --> y (used for calculating spectra)
393    IF ( dt_dosp /= 9999999.9 )  THEN
394       IF ( MOD( nza, pdims(2) ) /= 0 )  THEN
395          IF ( myid == 0 )  THEN
396             PRINT*,'+++ Direct transposition z --> y (needed for spectra):'
397             PRINT*,'    nz=',nz,' is not an integral divisor of ',&
398                    &'pdims(2)=',pdims(2)
399          ENDIF
400          CALL local_stop
401       ELSE
402          nxl_yd  = nxl
403          nxr_yda = nxra
404          nxr_yd  = nxr
405          nzb_yd  = 1 + myidy * ( nza / pdims(2) )
406          nzt_yda = ( myidy + 1 ) * ( nza / pdims(2) )
407          nzt_yd  = MIN( nzt, nzt_yda )
408
409          sendrecvcount_zyd = nnx * nny * ( nza / pdims(2) )
410       ENDIF
411    ENDIF
412
413!
414!-- Indices for direct transpositions y --> x (they are only possible in case
415!-- of a 1d-decomposition along x)
416    IF ( pdims(2) == 1 )  THEN
417       nny_x  = nny / pdims(1)
418       nys_x  = myid * nny_x
419       nyn_xa = ( myid + 1 ) * nny_x - 1
420       nyn_x  = MIN( ny, nyn_xa )
421       nzb_x  = 1
422       nzt_xa = nza
423       nzt_x  = nz
424       sendrecvcount_xy = nnx * nny_x * nza
425    ENDIF
426
427!
428!-- Indices for direct transpositions x --> y (they are only possible in case
429!-- of a 1d-decomposition along y)
430    IF ( pdims(1) == 1 )  THEN
431       nnx_y  = nnx / pdims(2)
432       nxl_y  = myid * nnx_y
433       nxr_ya = ( myid + 1 ) * nnx_y - 1
434       nxr_y  = MIN( nx, nxr_ya )
435       nzb_y  = 1
436       nzt_ya = nza
437       nzt_y  = nz
438       sendrecvcount_xy = nnx_y * nny * nza
439    ENDIF
440
441!
442!-- Arrays for storing the array bounds are needed any more
443    DEALLOCATE( nxlf , nxrf , nynf , nysf )
444
445#if defined( __print )
446!
447!-- Control output
448    IF ( myid == 0 )  THEN
449       PRINT*, '*** processor topology ***'
450       PRINT*, ' '
451       PRINT*, 'myid   pcoord    left right  south north  idx idy   nxl: nxr',&
452               &'   nys: nyn'
453       PRINT*, '------------------------------------------------------------',&
454               &'-----------'
455       WRITE (*,1000)  0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
456                       myidx, myidy, nxl, nxr, nys, nyn
4571000   FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
458               2(2X,I4,':',I4))
459
460!
461!--    Recieve data from the other PEs
462       DO  i = 1,numprocs-1
463          CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
464                         ierr )
465          WRITE (*,1000)  i, ( ibuf(j) , j = 1,12 )
466       ENDDO
467    ELSE
468
469!
470!--    Send data to PE0
471       ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
472       ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
473       ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
474       ibuf(12) = nyn
475       CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
476    ENDIF
477#endif
478
479#else
480
481!
482!-- Array bounds when running on a single PE (respectively a non-parallel
483!-- machine)
484    nxl  = 0
485    nxr  = nx
486    nxra = nx
487    nnx  = nxr - nxl + 1
488    nys  = 0
489    nyn  = ny
490    nyna = ny
491    nny  = nyn - nys + 1
492    nzb  = 0
493    nzt  = nz
494    nzta = nz
495    nnz  = nz
496
497!
498!-- Array bounds for the pressure solver (in the parallel code, these bounds
499!-- are the ones for the transposed arrays)
500    nys_x  = nys
501    nyn_x  = nyn
502    nyn_xa = nyn
503    nzb_x  = nzb + 1
504    nzt_x  = nzt
505    nzt_xa = nzt
506
507    nxl_y  = nxl
508    nxr_y  = nxr
509    nxr_ya = nxr
510    nzb_y  = nzb + 1
511    nzt_y  = nzt
512    nzt_ya = nzt
513
514    nxl_z  = nxl
515    nxr_z  = nxr
516    nxr_za = nxr
517    nys_z  = nys
518    nyn_z  = nyn
519    nyn_za = nyn
520
521#endif
522
523!
524!-- Calculate number of grid levels necessary for the multigrid poisson solver
525!-- as well as the gridpoint indices on each level
526    IF ( psolver == 'multigrid' )  THEN
527
528!
529!--    First calculate number of possible grid levels for the subdomains
530       mg_levels_x = 1
531       mg_levels_y = 1
532       mg_levels_z = 1
533
534       i = nnx
535       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
536          i = i / 2
537          mg_levels_x = mg_levels_x + 1
538       ENDDO
539
540       j = nny
541       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
542          j = j / 2
543          mg_levels_y = mg_levels_y + 1
544       ENDDO
545
546       k = nnz
547       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
548          k = k / 2
549          mg_levels_z = mg_levels_z + 1
550       ENDDO
551
552       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
553
554!
555!--    Find out, if the total domain allows more levels. These additional
556!--    levels are processed on PE0 only.
557       IF ( numprocs > 1 )  THEN
558          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
559             mg_switch_to_pe0_level_l = maximum_grid_level
560
561             mg_levels_x = 1
562             mg_levels_y = 1
563
564             i = nx+1
565             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
566                i = i / 2
567                mg_levels_x = mg_levels_x + 1
568             ENDDO
569
570             j = ny+1
571             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
572                j = j / 2
573                mg_levels_y = mg_levels_y + 1
574             ENDDO
575
576             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
577
578             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
579                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
580                                           mg_switch_to_pe0_level_l + 1
581             ELSE
582                mg_switch_to_pe0_level_l = 0
583             ENDIF
584          ELSE
585             mg_switch_to_pe0_level_l = 0
586             maximum_grid_level_l = maximum_grid_level
587          ENDIF
588
589!
590!--       Use switch level calculated above only if it is not pre-defined
591!--       by user
592          IF ( mg_switch_to_pe0_level == 0 )  THEN
593
594             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
595                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
596                maximum_grid_level     = maximum_grid_level_l
597             ENDIF
598
599          ELSE
600!
601!--          Check pre-defined value and reset to default, if neccessary
602             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.  &
603                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
604                IF ( myid == 0 )  THEN
605                   PRINT*, '+++ WARNING init_pegrid: mg_switch_to_pe0_level ', &
606                               'out of range and reset to default (=0)'
607                ENDIF
608                mg_switch_to_pe0_level = 0
609             ELSE
610!
611!--             Use the largest number of possible levels anyway and recalculate
612!--             the switch level to this largest number of possible values
613                maximum_grid_level = maximum_grid_level_l
614
615             ENDIF
616          ENDIF
617
618       ENDIF
619
620       ALLOCATE( grid_level_count(maximum_grid_level),                   &
621                 nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
622                 nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
623                 nzt_mg(maximum_grid_level) )
624
625       grid_level_count = 0
626       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
627
628       DO  i = maximum_grid_level, 1 , -1
629
630          IF ( i == mg_switch_to_pe0_level )  THEN
631#if defined( __parallel )
632!
633!--          Save the grid size of the subdomain at the switch level, because
634!--          it is needed in poismg.
635!--          Array bounds of the local subdomain grids are gathered on PE0
636             ind(1) = nxl_l; ind(2) = nxr_l
637             ind(3) = nys_l; ind(4) = nyn_l
638             ind(5) = nzt_l
639             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
640             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
641                                 MPI_INTEGER, comm2d, ierr )
642             DO  j = 0, numprocs-1
643                DO  k = 1, 5
644                   mg_loc_ind(k,j) = ind_all(k+j*5)
645                ENDDO
646             ENDDO
647             DEALLOCATE( ind_all )
648!
649!--          Calculate the grid size of the total domain gathered on PE0
650             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
651             nxl_l = 0
652             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
653             nys_l = 0
654!
655!--          The size of this gathered array must not be larger than the
656!--          array tend, which is used in the multigrid scheme as a temporary
657!--          array
658             subdomain_size = ( nxr - nxl + 3 )     * ( nyn - nys + 3 )     * &
659                              ( nzt - nzb + 2 )
660             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
661                              ( nzt_l - nzb + 2 )
662
663             IF ( gathered_size > subdomain_size )  THEN
664                IF ( myid == 0 )  THEN
665                   PRINT*, '+++ init_pegrid: not enough memory for storing ', &
666                               'gathered multigrid data on PE0'
667                ENDIF
668                CALL local_stop
669             ENDIF
670#else
671             PRINT*, '+++ init_pegrid: multigrid gather/scatter impossible ', &
672                          'in non parallel mode'
673             CALL local_stop
674#endif
675          ENDIF
676
677          nxl_mg(i) = nxl_l
678          nxr_mg(i) = nxr_l
679          nys_mg(i) = nys_l
680          nyn_mg(i) = nyn_l
681          nzt_mg(i) = nzt_l
682
683          nxl_l = nxl_l / 2 
684          nxr_l = nxr_l / 2
685          nys_l = nys_l / 2 
686          nyn_l = nyn_l / 2 
687          nzt_l = nzt_l / 2 
688       ENDDO
689
690    ELSE
691
692       maximum_grid_level = 1
693
694    ENDIF
695
696    grid_level = maximum_grid_level
697
698#if defined( __parallel )
699!
700!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
701    ngp_y  = nyn - nys + 1
702
703!
704!-- Define a new MPI derived datatype for the exchange of ghost points in
705!-- y-direction for 2D-arrays (line)
706    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_REAL, type_x, ierr )
707    CALL MPI_TYPE_COMMIT( type_x, ierr )
708    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_INTEGER, type_x_int, ierr )
709    CALL MPI_TYPE_COMMIT( type_x_int, ierr )
710
711!
712!-- Calculate gridpoint numbers for the exchange of ghost points along x
713!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
714!-- exchange of ghost points in y-direction (xz-plane).
715!-- Do these calculations for the model grid and (if necessary) also
716!-- for the coarser grid levels used in the multigrid method
717    ALLOCATE ( ngp_yz(maximum_grid_level), type_xz(maximum_grid_level) )
718
719    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
720         
721    DO i = maximum_grid_level, 1 , -1
722       ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
723
724       CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
725                             MPI_REAL, type_xz(i), ierr )
726       CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
727
728       nxl_l = nxl_l / 2 
729       nxr_l = nxr_l / 2
730       nys_l = nys_l / 2 
731       nyn_l = nyn_l / 2 
732       nzt_l = nzt_l / 2 
733    ENDDO
734#endif
735
736#if defined( __parallel )
737!
738!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
739!-- horizontal boundary conditions. Set variables for extending array u (v)
740!-- by one gridpoint on the left/rightmost (northest/southest) processor
741    IF ( pleft == MPI_PROC_NULL )  THEN
742       IF ( bc_lr == 'dirichlet/radiation' )  THEN
743          inflow_l  = .TRUE.
744       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
745          outflow_l = .TRUE.
746       ENDIF
747    ENDIF
748
749    IF ( pright == MPI_PROC_NULL )  THEN
750       IF ( bc_lr == 'dirichlet/radiation' )  THEN
751          outflow_r = .TRUE.
752       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
753          inflow_r  = .TRUE.
754       ENDIF
755    ENDIF
756
757    IF ( psouth == MPI_PROC_NULL )  THEN
758       IF ( bc_ns == 'dirichlet/radiation' )  THEN
759          outflow_s = .TRUE.
760       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
761          inflow_s  = .TRUE.
762       ENDIF
763    ENDIF
764
765    IF ( pnorth == MPI_PROC_NULL )  THEN
766       IF ( bc_ns == 'dirichlet/radiation' )  THEN
767          inflow_n  = .TRUE.
768       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
769          outflow_n = .TRUE.
770       ENDIF
771    ENDIF
772
773#else
774    IF ( bc_lr == 'dirichlet/radiation' )  THEN
775       inflow_l  = .TRUE.
776       outflow_r = .TRUE.
777    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
778       outflow_l = .TRUE.
779       inflow_r  = .TRUE.
780    ENDIF
781
782    IF ( bc_ns == 'dirichlet/radiation' )  THEN
783       inflow_n  = .TRUE.
784       outflow_s = .TRUE.
785    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
786       outflow_n = .TRUE.
787       inflow_s  = .TRUE.
788    ENDIF
789#endif
790
791    IF ( psolver == 'poisfft_hybrid' )  THEN
792       CALL poisfft_hybrid_ini
793    ELSEIF ( psolver == 'poisfft' )  THEN
794       CALL poisfft_init
795    ENDIF
796
797 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.