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

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

preliminary updates for implementing buildings in poismg

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