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

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

New:
---

Changed:


PALM can be generally installed on any kind of Linux-, IBM-AIX-, or NEC-SX-system by adding appropriate settings to the configuration file.

Scripts are also running under the public domain ksh.

All system relevant compile and link options as well as the host identifier (local_host) are specified in the configuration file.

Filetransfer by ftp removed (options -f removed from mrun and mbuild).

Call of (system-)FLUSH routine moved to new routine local_flush.

return_addres and return_username are read from ENVPAR-NAMELIST-file instead of using local_getenv.

Preprocessor strings for different linux clusters changed to "lc", some preprocessor directives renamed (new: intel_openmp_bug), preprocessor directives for old systems removed

advec_particles, check_open, cpu_log, cpu_statistics, data_output_dvrp, flow_statistics, header, init_dvrp, init_particles, init_1d_model, init_dvrp, init_pegrid, local_getenv, local_system, local_tremain, local_tremain_ini, modules, palm, parin, run_control

new:
local_flush

mbuild, mrun

Errors:


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