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

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

trying to make creation of intercommunicator more stable

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