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

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

preliminary version for coupled runs

  • Property svn:keywords set to Id
File size: 26.5 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 102 2007-07-27 09:09:17Z 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       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
498
499          CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
500
501       ENDIF
502
503    ENDIF
504
505!
506!-- In case of coupled runs, establish the connection between the atmosphere
507!-- and the ocean model and define the intercommunicator (comm_inter)
508    CALL MPI_BARRIER( comm2d, ierr )
509    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
510
511       print*, '... before COMM_ACCEPT'
512       CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
513                             comm_inter, ierr )
514       print*, '--- ierr = ', ierr
515       print*, '--- comm_inter atmosphere = ', comm_inter
516
517    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
518
519       IF ( myid == 0 )  PRINT*, '*** read: ', port_name, '  ierr = ', ierr
520       print*, '... before COMM_CONNECT'
521       CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
522                              comm_inter, ierr )
523       print*, '--- ierr = ', ierr
524       print*, '--- comm_inter ocean      = ', comm_inter
525
526    ENDIF
527
528!
529!-- In case of coupled runs, create a new MPI derived datatype for the
530!-- exchange of surface (xy) data .
531!-- Gridpoint number for the exchange of ghost points (xy-plane)
532    ngp_xy  = ( nxr - nxl + 3 ) * ( nyn - nys + 3 )
533
534!
535!-- Define a new MPI derived datatype for the exchange of ghost points in
536!-- y-direction for 2D-arrays (line)
537    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
538    CALL MPI_TYPE_COMMIT( type_xy, ierr )
539#endif
540
541#else
542
543!
544!-- Array bounds when running on a single PE (respectively a non-parallel
545!-- machine)
546    nxl  = 0
547    nxr  = nx
548    nxra = nx
549    nnx  = nxr - nxl + 1
550    nys  = 0
551    nyn  = ny
552    nyna = ny
553    nny  = nyn - nys + 1
554    nzb  = 0
555    nzt  = nz
556    nzta = nz
557    nnz  = nz
558
559!
560!-- Array bounds for the pressure solver (in the parallel code, these bounds
561!-- are the ones for the transposed arrays)
562    nys_x  = nys
563    nyn_x  = nyn
564    nyn_xa = nyn
565    nzb_x  = nzb + 1
566    nzt_x  = nzt
567    nzt_xa = nzt
568
569    nxl_y  = nxl
570    nxr_y  = nxr
571    nxr_ya = nxr
572    nzb_y  = nzb + 1
573    nzt_y  = nzt
574    nzt_ya = nzt
575
576    nxl_z  = nxl
577    nxr_z  = nxr
578    nxr_za = nxr
579    nys_z  = nys
580    nyn_z  = nyn
581    nyn_za = nyn
582
583#endif
584
585!
586!-- Calculate number of grid levels necessary for the multigrid poisson solver
587!-- as well as the gridpoint indices on each level
588    IF ( psolver == 'multigrid' )  THEN
589
590!
591!--    First calculate number of possible grid levels for the subdomains
592       mg_levels_x = 1
593       mg_levels_y = 1
594       mg_levels_z = 1
595
596       i = nnx
597       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
598          i = i / 2
599          mg_levels_x = mg_levels_x + 1
600       ENDDO
601
602       j = nny
603       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
604          j = j / 2
605          mg_levels_y = mg_levels_y + 1
606       ENDDO
607
608       k = nnz
609       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
610          k = k / 2
611          mg_levels_z = mg_levels_z + 1
612       ENDDO
613
614       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
615
616!
617!--    Find out, if the total domain allows more levels. These additional
618!--    levels are processed on PE0 only.
619       IF ( numprocs > 1 )  THEN
620          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
621             mg_switch_to_pe0_level_l = maximum_grid_level
622
623             mg_levels_x = 1
624             mg_levels_y = 1
625
626             i = nx+1
627             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
628                i = i / 2
629                mg_levels_x = mg_levels_x + 1
630             ENDDO
631
632             j = ny+1
633             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
634                j = j / 2
635                mg_levels_y = mg_levels_y + 1
636             ENDDO
637
638             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
639
640             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
641                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
642                                           mg_switch_to_pe0_level_l + 1
643             ELSE
644                mg_switch_to_pe0_level_l = 0
645             ENDIF
646          ELSE
647             mg_switch_to_pe0_level_l = 0
648             maximum_grid_level_l = maximum_grid_level
649          ENDIF
650
651!
652!--       Use switch level calculated above only if it is not pre-defined
653!--       by user
654          IF ( mg_switch_to_pe0_level == 0 )  THEN
655
656             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
657                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
658                maximum_grid_level     = maximum_grid_level_l
659             ENDIF
660
661          ELSE
662!
663!--          Check pre-defined value and reset to default, if neccessary
664             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.  &
665                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
666                IF ( myid == 0 )  THEN
667                   PRINT*, '+++ WARNING init_pegrid: mg_switch_to_pe0_level ', &
668                               'out of range and reset to default (=0)'
669                ENDIF
670                mg_switch_to_pe0_level = 0
671             ELSE
672!
673!--             Use the largest number of possible levels anyway and recalculate
674!--             the switch level to this largest number of possible values
675                maximum_grid_level = maximum_grid_level_l
676
677             ENDIF
678          ENDIF
679
680       ENDIF
681
682       ALLOCATE( grid_level_count(maximum_grid_level),                   &
683                 nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
684                 nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
685                 nzt_mg(maximum_grid_level) )
686
687       grid_level_count = 0
688       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
689
690       DO  i = maximum_grid_level, 1 , -1
691
692          IF ( i == mg_switch_to_pe0_level )  THEN
693#if defined( __parallel )
694!
695!--          Save the grid size of the subdomain at the switch level, because
696!--          it is needed in poismg.
697!--          Array bounds of the local subdomain grids are gathered on PE0
698             ind(1) = nxl_l; ind(2) = nxr_l
699             ind(3) = nys_l; ind(4) = nyn_l
700             ind(5) = nzt_l
701             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
702             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
703                                 MPI_INTEGER, comm2d, ierr )
704             DO  j = 0, numprocs-1
705                DO  k = 1, 5
706                   mg_loc_ind(k,j) = ind_all(k+j*5)
707                ENDDO
708             ENDDO
709             DEALLOCATE( ind_all )
710!
711!--          Calculate the grid size of the total domain gathered on PE0
712             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
713             nxl_l = 0
714             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
715             nys_l = 0
716!
717!--          The size of this gathered array must not be larger than the
718!--          array tend, which is used in the multigrid scheme as a temporary
719!--          array
720             subdomain_size = ( nxr - nxl + 3 )     * ( nyn - nys + 3 )     * &
721                              ( nzt - nzb + 2 )
722             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
723                              ( nzt_l - nzb + 2 )
724
725             IF ( gathered_size > subdomain_size )  THEN
726                IF ( myid == 0 )  THEN
727                   PRINT*, '+++ init_pegrid: not enough memory for storing ', &
728                               'gathered multigrid data on PE0'
729                ENDIF
730                CALL local_stop
731             ENDIF
732#else
733             PRINT*, '+++ init_pegrid: multigrid gather/scatter impossible ', &
734                          'in non parallel mode'
735             CALL local_stop
736#endif
737          ENDIF
738
739          nxl_mg(i) = nxl_l
740          nxr_mg(i) = nxr_l
741          nys_mg(i) = nys_l
742          nyn_mg(i) = nyn_l
743          nzt_mg(i) = nzt_l
744
745          nxl_l = nxl_l / 2 
746          nxr_l = nxr_l / 2
747          nys_l = nys_l / 2 
748          nyn_l = nyn_l / 2 
749          nzt_l = nzt_l / 2 
750       ENDDO
751
752    ELSE
753
754       maximum_grid_level = 1
755
756    ENDIF
757
758    grid_level = maximum_grid_level
759
760#if defined( __parallel )
761!
762!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
763    ngp_y  = nyn - nys + 1
764
765!
766!-- Define a new MPI derived datatype for the exchange of ghost points in
767!-- y-direction for 2D-arrays (line)
768    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_REAL, type_x, ierr )
769    CALL MPI_TYPE_COMMIT( type_x, ierr )
770    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_INTEGER, type_x_int, ierr )
771    CALL MPI_TYPE_COMMIT( type_x_int, ierr )
772
773!
774!-- Calculate gridpoint numbers for the exchange of ghost points along x
775!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
776!-- exchange of ghost points in y-direction (xz-plane).
777!-- Do these calculations for the model grid and (if necessary) also
778!-- for the coarser grid levels used in the multigrid method
779    ALLOCATE ( ngp_yz(maximum_grid_level), type_xz(maximum_grid_level) )
780
781    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
782         
783    DO i = maximum_grid_level, 1 , -1
784       ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
785
786       CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
787                             MPI_REAL, type_xz(i), ierr )
788       CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
789
790       nxl_l = nxl_l / 2 
791       nxr_l = nxr_l / 2
792       nys_l = nys_l / 2 
793       nyn_l = nyn_l / 2 
794       nzt_l = nzt_l / 2 
795    ENDDO
796#endif
797
798#if defined( __parallel )
799!
800!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
801!-- horizontal boundary conditions. Set variables for extending array u (v)
802!-- by one gridpoint on the left/rightmost (northest/southest) processor
803    IF ( pleft == MPI_PROC_NULL )  THEN
804       IF ( bc_lr == 'dirichlet/radiation' )  THEN
805          inflow_l  = .TRUE.
806       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
807          outflow_l = .TRUE.
808       ENDIF
809    ENDIF
810
811    IF ( pright == MPI_PROC_NULL )  THEN
812       IF ( bc_lr == 'dirichlet/radiation' )  THEN
813          outflow_r = .TRUE.
814       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
815          inflow_r  = .TRUE.
816       ENDIF
817    ENDIF
818
819    IF ( psouth == MPI_PROC_NULL )  THEN
820       IF ( bc_ns == 'dirichlet/radiation' )  THEN
821          outflow_s = .TRUE.
822       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
823          inflow_s  = .TRUE.
824       ENDIF
825    ENDIF
826
827    IF ( pnorth == MPI_PROC_NULL )  THEN
828       IF ( bc_ns == 'dirichlet/radiation' )  THEN
829          inflow_n  = .TRUE.
830       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
831          outflow_n = .TRUE.
832       ENDIF
833    ENDIF
834
835#else
836    IF ( bc_lr == 'dirichlet/radiation' )  THEN
837       inflow_l  = .TRUE.
838       outflow_r = .TRUE.
839    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
840       outflow_l = .TRUE.
841       inflow_r  = .TRUE.
842    ENDIF
843
844    IF ( bc_ns == 'dirichlet/radiation' )  THEN
845       inflow_n  = .TRUE.
846       outflow_s = .TRUE.
847    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
848       outflow_n = .TRUE.
849       inflow_s  = .TRUE.
850    ENDIF
851#endif
852
853    IF ( psolver == 'poisfft_hybrid' )  THEN
854       CALL poisfft_hybrid_ini
855    ELSEIF ( psolver == 'poisfft' )  THEN
856       CALL poisfft_init
857    ENDIF
858
859 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.