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

Last change on this file since 198 was 198, checked in by raasch, 13 years ago

file headers updated for the next release 3.5

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