Home

Context Navigation

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

Last change on this file since 197 was 197, checked in by raasch, 16 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 (', npexnpey, &
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) - ( nnxpdims(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) - ( nnxpdims(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
518	1000 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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |