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

Last change on this file since 154 was 151, checked in by raasch, 16 years ago

preliminary update for the turbulence recycling method

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