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

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

further adjustments for SGI and other small changes

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