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

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

vorlaeufige Standalone-Version fuer Linux-Cluster

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