Home

Context Navigation

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

Last change on this file since 1468 was 1468, checked in by maronga, 10 years ago
New flag files allow to force unscheduled termination/restarts of batch jobs, progress output is made for batch runs, small adjustments for lxce6 and lccrayh/lccrayb
Property svn:keywords set to `Id`
File size: 41.1 KB

Rev	Line
[1]	1	SUBROUTINE init_pegrid
[1036]	2
	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[254]	20	! Current revisions:
[1322]	21	! ------------------
[1468]	22	! Adapted for use on up to 6-digit processor cores
[1354]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: init_pegrid.f90 1468 2014-09-24 14:06:57Z maronga $
	27	!
[1436]	28	! 1435 2014-07-21 10:37:02Z keck
	29	! bugfix: added missing parameter coupling_mode_remote to ONLY-attribute
	30	!
[1403]	31	! 1402 2014-05-09 14:25:13Z raasch
	32	! location messages modified
	33	!
[1385]	34	! 1384 2014-05-02 14:31:06Z raasch
	35	! location messages added
	36	!
[1354]	37	! 1353 2014-04-08 15:21:23Z heinze
	38	! REAL constants provided with KIND-attribute
	39	!
[1323]	40	! 1322 2014-03-20 16:38:49Z raasch
	41	! REAL functions provided with KIND-attribute
	42	!
[1321]	43	! 1320 2014-03-20 08:40:49Z raasch
[1320]	44	! ONLY-attribute added to USE-statements,
	45	! kind-parameters added to all INTEGER and REAL declaration statements,
	46	! kinds are defined in new module kinds,
	47	! revision history before 2012 removed,
	48	! comment fields (!:) to be used for variable explanations added to
	49	! all variable declaration statements
[760]	50	!
[1305]	51	! 1304 2014-03-12 10:29:42Z raasch
	52	! bugfix: single core MPI runs missed some settings of transpose indices
	53	!
[1213]	54	! 1212 2013-08-15 08:46:27Z raasch
	55	! error message for poisfft_hybrid removed
	56	!
[1160]	57	! 1159 2013-05-21 11:58:22Z fricke
	58	! dirichlet/neumann and neumann/dirichlet removed
	59	!
[1140]	60	! 1139 2013-04-18 07:25:03Z raasch
	61	! bugfix for calculating the id of the PE carrying the recycling plane
	62	!
[1112]	63	! 1111 2013-03-08 23:54:10Z raasch
	64	! initialization of poisfft moved to module poisfft
	65	!
[1093]	66	! 1092 2013-02-02 11:24:22Z raasch
	67	! unused variables removed
	68	!
[1057]	69	! 1056 2012-11-16 15:28:04Z raasch
	70	! Indices for arrays n.._mg start from zero due to definition of arrays f2 and
	71	! p2 as automatic arrays in recursive subroutine next_mg_level
	72	!
[1042]	73	! 1041 2012-11-06 02:36:29Z raasch
	74	! a 2d virtual processor topology is used by default for all machines
	75	!
[1037]	76	! 1036 2012-10-22 13:43:42Z raasch
	77	! code put under GPL (PALM 3.9)
	78	!
[1004]	79	! 1003 2012-09-14 14:35:53Z raasch
	80	! subdomains must have identical size (grid matching = "match" removed)
	81	!
[1002]	82	! 1001 2012-09-13 14:08:46Z raasch
	83	! all actions concerning upstream-spline-method removed
	84	!
[979]	85	! 978 2012-08-09 08:28:32Z fricke
	86	! dirichlet/neumann and neumann/dirichlet added
	87	! nxlu and nysv are also calculated for inflow boundary
	88	!
[810]	89	! 809 2012-01-30 13:32:58Z maronga
	90	! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
	91	!
[808]	92	! 807 2012-01-25 11:53:51Z maronga
	93	! New cpp directive "__check" implemented which is used by check_namelist_files
	94	!
[1]	95	! Revision 1.1 1997/07/24 11:15:09 raasch
	96	! Initial revision
	97	!
	98	!
	99	! Description:
	100	! ------------
	101	! Determination of the virtual processor topology (if not prescribed by the
	102	! user)and computation of the grid point number and array bounds of the local
	103	! domains.
	104	!------------------------------------------------------------------------------!
	105
[1320]	106	USE control_parameters, &
[1435]	107	ONLY: bc_lr, bc_ns, coupling_mode, coupling_mode_remote, &
	108	coupling_topology, dt_dosp, gathered_size, grid_level, &
	109	grid_level_count, host, inflow_l, inflow_n, inflow_r, inflow_s, &
	110	io_blocks, io_group, maximum_grid_level, &
	111	maximum_parallel_io_streams, message_string, &
[1320]	112	mg_switch_to_pe0_level, momentum_advec, psolver, outflow_l, &
	113	outflow_n, outflow_r, outflow_s, recycling_width, scalar_advec, &
	114	subdomain_size
[1]	115
[1320]	116	USE grid_variables, &
	117	ONLY: dx
	118
	119	USE indices, &
	120	ONLY: mg_loc_ind, nbgp, nnx, nny, nnz, nx, nx_a, nx_o, nxl, nxl_mg, &
	121	nxlu, nxr, nxr_mg, ny, ny_a, ny_o, nyn, nyn_mg, nys, nys_mg, &
	122	nysv, nz, nzb, nzt, nzt_mg, wall_flags_1, wall_flags_2, &
	123	wall_flags_3, wall_flags_4, wall_flags_5, wall_flags_6, &
	124	wall_flags_7, wall_flags_8, wall_flags_9, wall_flags_10
[1]	125
[1320]	126	USE kinds
	127
	128	USE pegrid
	129
	130	USE transpose_indices, &
	131	ONLY: nxl_y, nxl_yd, nxl_z, nxr_y, nxr_yd, nxr_z, nyn_x, nyn_z, nys_x,&
	132	nys_z, nzb_x, nzb_y, nzb_yd, nzt_x, nzt_yd, nzt_y
[667]	133
[1]	134	IMPLICIT NONE
	135
[1320]	136	INTEGER(iwp) :: i !:
	137	INTEGER(iwp) :: id_inflow_l !:
	138	INTEGER(iwp) :: id_recycling_l !:
	139	INTEGER(iwp) :: ind(5) !:
	140	INTEGER(iwp) :: j !:
	141	INTEGER(iwp) :: k !:
	142	INTEGER(iwp) :: maximum_grid_level_l !:
	143	INTEGER(iwp) :: mg_levels_x !:
	144	INTEGER(iwp) :: mg_levels_y !:
	145	INTEGER(iwp) :: mg_levels_z !:
	146	INTEGER(iwp) :: mg_switch_to_pe0_level_l !:
	147	INTEGER(iwp) :: nnx_y !:
	148	INTEGER(iwp) :: nnx_z !:
	149	INTEGER(iwp) :: nny_x !:
	150	INTEGER(iwp) :: nny_z !:
	151	INTEGER(iwp) :: nnz_x !:
	152	INTEGER(iwp) :: nnz_y !:
	153	INTEGER(iwp) :: numproc_sqr !:
	154	INTEGER(iwp) :: nxl_l !:
	155	INTEGER(iwp) :: nxr_l !:
	156	INTEGER(iwp) :: nyn_l !:
	157	INTEGER(iwp) :: nys_l !:
	158	INTEGER(iwp) :: nzb_l !:
	159	INTEGER(iwp) :: nzt_l !:
	160	INTEGER(iwp) :: omp_get_num_threads !:
[1]	161
[1320]	162	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ind_all !:
	163	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxlf !:
	164	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxrf !:
	165	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nynf !:
	166	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nysf !:
[1]	167
[1320]	168	INTEGER(iwp), DIMENSION(2) :: pdims_remote !:
[667]	169
[1092]	170	#if defined( __mpi2 )
[1320]	171	LOGICAL :: found !:
[1092]	172	#endif
[1]	173
	174	!
	175	!-- Get the number of OpenMP threads
	176	!$OMP PARALLEL
[82]	177	#if defined( __intel_openmp_bug )
[1]	178	threads_per_task = omp_get_num_threads()
	179	#else
	180	!$ threads_per_task = omp_get_num_threads()
	181	#endif
	182	!$OMP END PARALLEL
	183
	184
	185	#if defined( __parallel )
[667]	186
[1402]	187	CALL location_message( 'creating virtual PE grids + MPI derived data types', &
	188	.FALSE. )
[1]	189	!
	190	!-- Determine the processor topology or check it, if prescribed by the user
	191	IF ( npex == -1 .AND. npey == -1 ) THEN
	192
	193	!
	194	!-- Automatic determination of the topology
[1322]	195	numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
[1041]	196	pdims(1) = MAX( numproc_sqr , 1 )
	197	DO WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
	198	pdims(1) = pdims(1) - 1
	199	ENDDO
	200	pdims(2) = numprocs / pdims(1)
[1]	201
	202	ELSEIF ( npex /= -1 .AND. npey /= -1 ) THEN
	203
	204	!
	205	!-- Prescribed by user. Number of processors on the prescribed topology
	206	!-- must be equal to the number of PEs available to the job
	207	IF ( ( npex * npey ) /= numprocs ) THEN
[274]	208	WRITE( message_string, * ) 'number of PEs of the prescribed ', &
	209	'topology (', npex*npey,') does not match & the number of ', &
	210	'PEs available to the job (', numprocs, ')'
[254]	211	CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
[1]	212	ENDIF
	213	pdims(1) = npex
	214	pdims(2) = npey
	215
	216	ELSE
	217	!
	218	!-- If the processor topology is prescribed by the user, the number of
	219	!-- PEs must be given in both directions
[274]	220	message_string = 'if the processor topology is prescribed by the, ' // &
	221	' user& both values of "npex" and "npey" must be given ' // &
	222	'in the &NAMELIST-parameter file'
[254]	223	CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
[1]	224
	225	ENDIF
	226
	227	!
[622]	228	!-- For communication speedup, set barriers in front of collective
	229	!-- communications by default on SGI-type systems
	230	IF ( host(3:5) == 'sgi' ) collective_wait = .TRUE.
	231
	232	!
[1]	233	!-- If necessary, set horizontal boundary conditions to non-cyclic
[722]	234	IF ( bc_lr /= 'cyclic' ) cyclic(1) = .FALSE.
	235	IF ( bc_ns /= 'cyclic' ) cyclic(2) = .FALSE.
[1]	236
[807]	237
[809]	238	#if ! defined( __check)
[1]	239	!
	240	!-- Create the virtual processor grid
	241	CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
	242	comm2d, ierr )
	243	CALL MPI_COMM_RANK( comm2d, myid, ierr )
[1468]	244	WRITE (myid_char,'(''_'',I6.6)') myid
[1]	245
	246	CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
	247	CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
	248	CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
	249
	250	!
	251	!-- Determine sub-topologies for transpositions
	252	!-- Transposition from z to x:
	253	remain_dims(1) = .TRUE.
	254	remain_dims(2) = .FALSE.
	255	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
	256	CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
	257	!
	258	!-- Transposition from x to y
	259	remain_dims(1) = .FALSE.
	260	remain_dims(2) = .TRUE.
	261	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
	262	CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
	263
[807]	264	#endif
[1]	265
	266	!
[1003]	267	!-- Calculate array bounds along x-direction for every PE.
[1]	268	ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
[1003]	269	nysf(0:pdims(2)-1) )
[1]	270
[1003]	271	IF ( MOD( nx+1 , pdims(1) ) /= 0 ) THEN
[274]	272	WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ',&
	273	'is not an& integral divisor of the number ', &
	274	'processors (', pdims(1),')'
[254]	275	CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
[1]	276	ELSE
[1003]	277	nnx = ( nx + 1 ) / pdims(1)
[1]	278	IF ( nnx*pdims(1) - ( nx + 1) > nnx ) THEN
[274]	279	WRITE( message_string, * ) 'x-direction: nx does not match the', &
	280	'requirements given by the number of PEs &used', &
	281	'& please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
	282	- ( nx + 1 ) ) ), ' instead of nx =', nx
[254]	283	CALL message( 'init_pegrid', 'PA0226', 1, 2, 0, 6, 0 )
[1]	284	ENDIF
	285	ENDIF
	286
	287	!
	288	!-- Left and right array bounds, number of gridpoints
	289	DO i = 0, pdims(1)-1
	290	nxlf(i) = i * nnx
	291	nxrf(i) = ( i + 1 ) * nnx - 1
	292	ENDDO
	293
	294	!
	295	!-- Calculate array bounds in y-direction for every PE.
[1003]	296	IF ( MOD( ny+1 , pdims(2) ) /= 0 ) THEN
[274]	297	WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
	298	'is not an& integral divisor of the number of', &
	299	'processors (', pdims(2),')'
[254]	300	CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
[1]	301	ELSE
[1003]	302	nny = ( ny + 1 ) / pdims(2)
[1]	303	IF ( nny*pdims(2) - ( ny + 1) > nny ) THEN
[274]	304	WRITE( message_string, * ) 'y-direction: ny does not match the', &
	305	'requirements given by the number of PEs &used ', &
	306	'& please use ny = ', ny - ( pdims(2) - ( nnx*pdims(2) &
[254]	307	- ( ny + 1 ) ) ), ' instead of ny =', ny
	308	CALL message( 'init_pegrid', 'PA0228', 1, 2, 0, 6, 0 )
[1]	309	ENDIF
	310	ENDIF
	311
	312	!
	313	!-- South and north array bounds
	314	DO j = 0, pdims(2)-1
	315	nysf(j) = j * nny
	316	nynf(j) = ( j + 1 ) * nny - 1
	317	ENDDO
	318
	319	!
	320	!-- Local array bounds of the respective PEs
[1003]	321	nxl = nxlf(pcoord(1))
	322	nxr = nxrf(pcoord(1))
	323	nys = nysf(pcoord(2))
	324	nyn = nynf(pcoord(2))
	325	nzb = 0
	326	nzt = nz
	327	nnz = nz
[1]	328
	329	!
[707]	330	!-- Set switches to define if the PE is situated at the border of the virtual
	331	!-- processor grid
	332	IF ( nxl == 0 ) left_border_pe = .TRUE.
	333	IF ( nxr == nx ) right_border_pe = .TRUE.
	334	IF ( nys == 0 ) south_border_pe = .TRUE.
	335	IF ( nyn == ny ) north_border_pe = .TRUE.
	336
	337	!
[1]	338	!-- Calculate array bounds and gridpoint numbers for the transposed arrays
	339	!-- (needed in the pressure solver)
	340	!-- For the transposed arrays, cyclic boundaries as well as top and bottom
	341	!-- boundaries are omitted, because they are obstructive to the transposition
	342
	343	!
	344	!-- 1. transposition z --> x
[1001]	345	!-- This transposition is not neccessary in case of a 1d-decomposition along x
[1304]	346	nys_x = nys
	347	nyn_x = nyn
	348	nny_x = nny
	349	nnz_x = nz / pdims(1)
	350	nzb_x = 1 + myidx * nnz_x
	351	nzt_x = ( myidx + 1 ) * nnz_x
	352	sendrecvcount_zx = nnx * nny * nnz_x
	353
[1001]	354	IF ( pdims(2) /= 1 ) THEN
[1003]	355	IF ( MOD( nz , pdims(1) ) /= 0 ) THEN
[274]	356	WRITE( message_string, * ) 'transposition z --> x:', &
	357	'&nz=',nz,' is not an integral divisior of pdims(1)=', &
	358	pdims(1)
[254]	359	CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
[1]	360	ENDIF
	361	ENDIF
	362
	363	!
	364	!-- 2. transposition x --> y
[1003]	365	nnz_y = nnz_x
	366	nzb_y = nzb_x
	367	nzt_y = nzt_x
	368	IF ( MOD( nx+1 , pdims(2) ) /= 0 ) THEN
[274]	369	WRITE( message_string, * ) 'transposition x --> y:', &
	370	'&nx+1=',nx+1,' is not an integral divisor of ',&
	371	'pdims(2)=',pdims(2)
[254]	372	CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
[1]	373	ENDIF
[1003]	374	nnx_y = (nx+1) / pdims(2)
[1]	375	nxl_y = myidy * nnx_y
[1003]	376	nxr_y = ( myidy + 1 ) * nnx_y - 1
[1]	377	sendrecvcount_xy = nnx_y * nny_x * nnz_y
	378
	379	!
	380	!-- 3. transposition y --> z (ELSE: x --> y in case of 1D-decomposition
	381	!-- along x)
[1304]	382	nnx_z = nnx_y
	383	nxl_z = nxl_y
	384	nxr_z = nxr_y
	385	nny_z = (ny+1) / pdims(1)
	386	nys_z = myidx * nny_z
	387	nyn_z = ( myidx + 1 ) * nny_z - 1
	388	sendrecvcount_yz = nnx_y * nny_z * nnz_y
	389
[1001]	390	IF ( pdims(2) /= 1 ) THEN
[1]	391	!
	392	!-- y --> z
	393	!-- This transposition is not neccessary in case of a 1d-decomposition
	394	!-- along x, except that the uptream-spline method is switched on
[1003]	395	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[274]	396	WRITE( message_string, * ) 'transposition y --> z:', &
	397	'& ny+1=',ny+1,' is not an integral divisor of ',&
	398	'pdims(1)=',pdims(1)
[254]	399	CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
[1]	400	ENDIF
	401
	402	ELSE
	403	!
	404	!-- x --> y. This condition must be fulfilled for a 1D-decomposition along x
[1003]	405	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[274]	406	WRITE( message_string, * ) 'transposition x --> y:', &
	407	'& ny+1=',ny+1,' is not an integral divisor of ',&
	408	'pdims(1)=',pdims(1)
[254]	409	CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
[1]	410	ENDIF
	411
	412	ENDIF
	413
	414	!
	415	!-- Indices for direct transpositions z --> y (used for calculating spectra)
[1353]	416	IF ( dt_dosp /= 9999999.9_wp ) THEN
[1003]	417	IF ( MOD( nz, pdims(2) ) /= 0 ) THEN
[274]	418	WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
	419	'for spectra):& nz=',nz,' is not an integral divisor of ',&
	420	'pdims(2)=',pdims(2)
[254]	421	CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
[1]	422	ELSE
[1003]	423	nxl_yd = nxl
	424	nxr_yd = nxr
	425	nzb_yd = 1 + myidy * ( nz / pdims(2) )
	426	nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
	427	sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
[1]	428	ENDIF
	429	ENDIF
	430
	431	!
	432	!-- Indices for direct transpositions y --> x (they are only possible in case
	433	!-- of a 1d-decomposition along x)
	434	IF ( pdims(2) == 1 ) THEN
[1003]	435	nny_x = nny / pdims(1)
	436	nys_x = myid * nny_x
	437	nyn_x = ( myid + 1 ) * nny_x - 1
	438	nzb_x = 1
	439	nzt_x = nz
	440	sendrecvcount_xy = nnx * nny_x * nz
[1]	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
[1003]	447	nnx_y = nnx / pdims(2)
	448	nxl_y = myid * nnx_y
	449	nxr_y = ( myid + 1 ) * nnx_y - 1
	450	nzb_y = 1
	451	nzt_y = nz
	452	sendrecvcount_xy = nnx_y * nny * nz
[1]	453	ENDIF
	454
	455	!
	456	!-- Arrays for storing the array bounds are needed any more
	457	DEALLOCATE( nxlf , nxrf , nynf , nysf )
	458
[807]	459
[809]	460	#if ! defined( __check)
[145]	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
[807]	491	#endif
	492
[1]	493	#if defined( __print )
	494	!
	495	!-- Control output
	496	IF ( myid == 0 ) THEN
	497	PRINT, ' processor topology *'
	498	PRINT*, ' '
	499	PRINT*, 'myid pcoord left right south north idx idy nxl: nxr',&
	500	&' nys: nyn'
	501	PRINT*, '------------------------------------------------------------',&
	502	&'-----------'
	503	WRITE (*,1000) 0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
	504	myidx, myidy, nxl, nxr, nys, nyn
	505	1000 FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
	506	2(2X,I4,':',I4))
	507
	508	!
[108]	509	!-- Receive data from the other PEs
[1]	510	DO i = 1,numprocs-1
	511	CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
	512	ierr )
	513	WRITE (*,1000) i, ( ibuf(j) , j = 1,12 )
	514	ENDDO
	515	ELSE
	516
	517	!
	518	!-- Send data to PE0
	519	ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
	520	ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
	521	ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
	522	ibuf(12) = nyn
	523	CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )
	524	ENDIF
	525	#endif
	526
[809]	527	#if defined( __parallel ) && ! defined( __check)
[102]	528	#if defined( __mpi2 )
	529	!
	530	!-- In case of coupled runs, get the port name on PE0 of the atmosphere model
	531	!-- and pass it to PE0 of the ocean model
	532	IF ( myid == 0 ) THEN
	533
	534	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	535
	536	CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
[108]	537
[102]	538	CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
	539	ierr )
[108]	540
	541	!
[104]	542	!-- Write a flag file for the ocean model and the other atmosphere
	543	!-- processes.
	544	!-- There seems to be a bug in MPICH2 which causes hanging processes
	545	!-- in case that execution of LOOKUP_NAME is continued too early
	546	!-- (i.e. before the port has been created)
	547	OPEN( 90, FILE='COUPLING_PORT_OPENED', FORM='FORMATTED' )
	548	WRITE ( 90, '(''TRUE'')' )
	549	CLOSE ( 90 )
[102]	550
	551	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	552
[104]	553	!
	554	!-- Continue only if the atmosphere model has created the port.
	555	!-- There seems to be a bug in MPICH2 which causes hanging processes
	556	!-- in case that execution of LOOKUP_NAME is continued too early
	557	!-- (i.e. before the port has been created)
	558	INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
	559	DO WHILE ( .NOT. found )
	560	INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
	561	ENDDO
	562
[102]	563	CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
	564
	565	ENDIF
	566
	567	ENDIF
	568
	569	!
	570	!-- In case of coupled runs, establish the connection between the atmosphere
	571	!-- and the ocean model and define the intercommunicator (comm_inter)
	572	CALL MPI_BARRIER( comm2d, ierr )
	573	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	574
	575	CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
	576	comm_inter, ierr )
[108]	577	coupling_mode_remote = 'ocean_to_atmosphere'
	578
[102]	579	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	580
	581	CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
	582	comm_inter, ierr )
[108]	583	coupling_mode_remote = 'atmosphere_to_ocean'
	584
[102]	585	ENDIF
[206]	586	#endif
[102]	587
[667]	588	!
[709]	589	!-- Determine the number of ghost point layers
	590	IF ( scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) THEN
[667]	591	nbgp = 3
	592	ELSE
	593	nbgp = 1
[709]	594	ENDIF
[667]	595
[102]	596	!
[709]	597	!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
	598	!-- which is needed for coupled atmosphere-ocean runs.
	599	!-- First, calculate number of grid points of an xy-plane.
[667]	600	ngp_xy = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
[102]	601	CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
	602	CALL MPI_TYPE_COMMIT( type_xy, ierr )
[667]	603
[709]	604	IF ( TRIM( coupling_mode ) /= 'uncoupled' ) THEN
[667]	605
	606	!
	607	!-- Pass the number of grid points of the atmosphere model to
	608	!-- the ocean model and vice versa
	609	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	610
	611	nx_a = nx
	612	ny_a = ny
	613
[709]	614	IF ( myid == 0 ) THEN
	615
	616	CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter, &
	617	ierr )
	618	CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter, &
	619	ierr )
	620	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
	621	ierr )
	622	CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter, &
	623	status, ierr )
	624	CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter, &
	625	status, ierr )
	626	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
[667]	627	comm_inter, status, ierr )
	628	ENDIF
	629
[709]	630	CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
	631	CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr )
	632	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
[667]	633
	634	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	635
	636	nx_o = nx
	637	ny_o = ny
	638
	639	IF ( myid == 0 ) THEN
[709]	640
	641	CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
	642	ierr )
	643	CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
	644	ierr )
	645	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
	646	status, ierr )
	647	CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
	648	CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
	649	CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
[667]	650	ENDIF
	651
	652	CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	653	CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	654	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
	655
	656	ENDIF
	657
[709]	658	ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
	659	ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
[667]	660
	661	!
[709]	662	!-- Determine if the horizontal grid and the number of PEs in ocean and
	663	!-- atmosphere is same or not
	664	IF ( nx_o == nx_a .AND. ny_o == ny_a .AND. &
[667]	665	pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
	666	THEN
	667	coupling_topology = 0
	668	ELSE
	669	coupling_topology = 1
	670	ENDIF
	671
	672	!
	673	!-- Determine the target PEs for the exchange between ocean and
	674	!-- atmosphere (comm2d)
[709]	675	IF ( coupling_topology == 0 ) THEN
	676	!
	677	!-- In case of identical topologies, every atmosphere PE has exactly one
	678	!-- ocean PE counterpart and vice versa
	679	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	680	target_id = myid + numprocs
	681	ELSE
	682	target_id = myid
	683	ENDIF
	684
	685	ELSE
	686	!
	687	!-- In case of nonequivalent topology in ocean and atmosphere only for
	688	!-- PE0 in ocean and PE0 in atmosphere a target_id is needed, since
[709]	689	!-- data echxchange between ocean and atmosphere will be done only
	690	!-- between these PEs.
	691	IF ( myid == 0 ) THEN
	692
	693	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	694	target_id = numprocs
	695	ELSE
	696	target_id = 0
	697	ENDIF
[709]	698
[667]	699	ENDIF
[709]	700
[667]	701	ENDIF
	702
	703	ENDIF
	704
	705
[102]	706	#endif
	707
[1]	708	#else
	709
	710	!
	711	!-- Array bounds when running on a single PE (respectively a non-parallel
	712	!-- machine)
[1003]	713	nxl = 0
	714	nxr = nx
	715	nnx = nxr - nxl + 1
	716	nys = 0
	717	nyn = ny
	718	nny = nyn - nys + 1
	719	nzb = 0
	720	nzt = nz
	721	nnz = nz
[1]	722
[145]	723	ALLOCATE( hor_index_bounds(4,0:0) )
	724	hor_index_bounds(1,0) = nxl
	725	hor_index_bounds(2,0) = nxr
	726	hor_index_bounds(3,0) = nys
	727	hor_index_bounds(4,0) = nyn
	728
[1]	729	!
	730	!-- Array bounds for the pressure solver (in the parallel code, these bounds
	731	!-- are the ones for the transposed arrays)
[1003]	732	nys_x = nys
	733	nyn_x = nyn
	734	nzb_x = nzb + 1
	735	nzt_x = nzt
[1]	736
[1003]	737	nxl_y = nxl
	738	nxr_y = nxr
	739	nzb_y = nzb + 1
	740	nzt_y = nzt
[1]	741
[1003]	742	nxl_z = nxl
	743	nxr_z = nxr
	744	nys_z = nys
	745	nyn_z = nyn
[1]	746
	747	#endif
	748
	749	!
	750	!-- Calculate number of grid levels necessary for the multigrid poisson solver
	751	!-- as well as the gridpoint indices on each level
	752	IF ( psolver == 'multigrid' ) THEN
	753
	754	!
	755	!-- First calculate number of possible grid levels for the subdomains
	756	mg_levels_x = 1
	757	mg_levels_y = 1
	758	mg_levels_z = 1
	759
	760	i = nnx
	761	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	762	i = i / 2
	763	mg_levels_x = mg_levels_x + 1
	764	ENDDO
	765
	766	j = nny
	767	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	768	j = j / 2
	769	mg_levels_y = mg_levels_y + 1
	770	ENDDO
	771
[181]	772	k = nz ! do not use nnz because it might be > nz due to transposition
	773	! requirements
[1]	774	DO WHILE ( MOD( k, 2 ) == 0 .AND. k /= 2 )
	775	k = k / 2
	776	mg_levels_z = mg_levels_z + 1
	777	ENDDO
	778
	779	maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
	780
	781	!
	782	!-- Find out, if the total domain allows more levels. These additional
[709]	783	!-- levels are identically processed on all PEs.
[197]	784	IF ( numprocs > 1 .AND. mg_switch_to_pe0_level /= -1 ) THEN
[709]	785
[1]	786	IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) ) THEN
[709]	787
[1]	788	mg_switch_to_pe0_level_l = maximum_grid_level
	789
	790	mg_levels_x = 1
	791	mg_levels_y = 1
	792
	793	i = nx+1
	794	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	795	i = i / 2
	796	mg_levels_x = mg_levels_x + 1
	797	ENDDO
	798
	799	j = ny+1
	800	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	801	j = j / 2
	802	mg_levels_y = mg_levels_y + 1
	803	ENDDO
	804
	805	maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
	806
	807	IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l ) THEN
	808	mg_switch_to_pe0_level_l = maximum_grid_level_l - &
	809	mg_switch_to_pe0_level_l + 1
	810	ELSE
	811	mg_switch_to_pe0_level_l = 0
	812	ENDIF
[709]	813
[1]	814	ELSE
	815	mg_switch_to_pe0_level_l = 0
	816	maximum_grid_level_l = maximum_grid_level
[709]	817
[1]	818	ENDIF
	819
	820	!
	821	!-- Use switch level calculated above only if it is not pre-defined
	822	!-- by user
	823	IF ( mg_switch_to_pe0_level == 0 ) THEN
	824	IF ( mg_switch_to_pe0_level_l /= 0 ) THEN
	825	mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
	826	maximum_grid_level = maximum_grid_level_l
	827	ENDIF
	828
	829	ELSE
	830	!
	831	!-- Check pre-defined value and reset to default, if neccessary
	832	IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l .OR. &
	833	mg_switch_to_pe0_level >= maximum_grid_level_l ) THEN
[254]	834	message_string = 'mg_switch_to_pe0_level ' // &
	835	'out of range and reset to default (=0)'
	836	CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
[1]	837	mg_switch_to_pe0_level = 0
	838	ELSE
	839	!
	840	!-- Use the largest number of possible levels anyway and recalculate
	841	!-- the switch level to this largest number of possible values
	842	maximum_grid_level = maximum_grid_level_l
	843
	844	ENDIF
[709]	845
[1]	846	ENDIF
	847
	848	ENDIF
	849
[1056]	850	ALLOCATE( grid_level_count(maximum_grid_level), &
	851	nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
	852	nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
	853	nzt_mg(0:maximum_grid_level) )
[1]	854
	855	grid_level_count = 0
[1056]	856	!
	857	!-- Index zero required as dummy due to definition of arrays f2 and p2 in
	858	!-- recursive subroutine next_mg_level
	859	nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
[778]	860
[1]	861	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
	862
	863	DO i = maximum_grid_level, 1 , -1
	864
	865	IF ( i == mg_switch_to_pe0_level ) THEN
[809]	866	#if defined( __parallel ) && ! defined( __check )
[1]	867	!
	868	!-- Save the grid size of the subdomain at the switch level, because
	869	!-- it is needed in poismg.
	870	ind(1) = nxl_l; ind(2) = nxr_l
	871	ind(3) = nys_l; ind(4) = nyn_l
	872	ind(5) = nzt_l
	873	ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
	874	CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
	875	MPI_INTEGER, comm2d, ierr )
	876	DO j = 0, numprocs-1
	877	DO k = 1, 5
	878	mg_loc_ind(k,j) = ind_all(k+j*5)
	879	ENDDO
	880	ENDDO
	881	DEALLOCATE( ind_all )
	882	!
[709]	883	!-- Calculate the grid size of the total domain
[1]	884	nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
	885	nxl_l = 0
	886	nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
	887	nys_l = 0
	888	!
	889	!-- The size of this gathered array must not be larger than the
	890	!-- array tend, which is used in the multigrid scheme as a temporary
[778]	891	!-- array. Therefore the subdomain size of an PE is calculated and
	892	!-- the size of the gathered grid. These values are used in
	893	!-- routines pres and poismg
	894	subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
	895	( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
[1]	896	gathered_size = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
	897	( nzt_l - nzb + 2 )
	898
[809]	899	#elif ! defined ( __parallel )
[254]	900	message_string = 'multigrid gather/scatter impossible ' // &
[1]	901	'in non parallel mode'
[254]	902	CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
[1]	903	#endif
	904	ENDIF
	905
	906	nxl_mg(i) = nxl_l
	907	nxr_mg(i) = nxr_l
	908	nys_mg(i) = nys_l
	909	nyn_mg(i) = nyn_l
	910	nzt_mg(i) = nzt_l
	911
	912	nxl_l = nxl_l / 2
	913	nxr_l = nxr_l / 2
	914	nys_l = nys_l / 2
	915	nyn_l = nyn_l / 2
	916	nzt_l = nzt_l / 2
[778]	917
[1]	918	ENDDO
	919
[780]	920	!
	921	!-- Temporary problem: Currently calculation of maxerror iin routine poismg crashes
	922	!-- if grid data are collected on PE0 already on the finest grid level.
	923	!-- To be solved later.
	924	IF ( maximum_grid_level == mg_switch_to_pe0_level ) THEN
	925	message_string = 'grid coarsening on subdomain level cannot be performed'
	926	CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
	927	ENDIF
	928
[1]	929	ELSE
	930
[667]	931	maximum_grid_level = 0
[1]	932
	933	ENDIF
	934
[722]	935	!
	936	!-- Default level 0 tells exchange_horiz that all ghost planes have to be
	937	!-- exchanged. grid_level is adjusted in poismg, where only one ghost plane
	938	!-- is required.
	939	grid_level = 0
[1]	940
[809]	941	#if defined( __parallel ) && ! defined ( __check )
[1]	942	!
	943	!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
[667]	944	ngp_y = nyn - nys + 1 + 2 * nbgp
[1]	945
	946	!
[709]	947	!-- Define new MPI derived datatypes for the exchange of ghost points in
	948	!-- x- and y-direction for 2D-arrays (line)
	949	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, &
	950	ierr )
[1]	951	CALL MPI_TYPE_COMMIT( type_x, ierr )
[709]	952	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, &
	953	type_x_int, ierr )
[1]	954	CALL MPI_TYPE_COMMIT( type_x_int, ierr )
	955
[667]	956	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
	957	CALL MPI_TYPE_COMMIT( type_y, ierr )
	958	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int, ierr )
	959	CALL MPI_TYPE_COMMIT( type_y_int, ierr )
	960
	961
[1]	962	!
	963	!-- Calculate gridpoint numbers for the exchange of ghost points along x
	964	!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
	965	!-- exchange of ghost points in y-direction (xz-plane).
	966	!-- Do these calculations for the model grid and (if necessary) also
	967	!-- for the coarser grid levels used in the multigrid method
[667]	968	ALLOCATE ( ngp_yz(0:maximum_grid_level), type_xz(0:maximum_grid_level),&
	969	type_yz(0:maximum_grid_level) )
[1]	970
	971	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
[709]	972
[667]	973	!
	974	!-- Discern between the model grid, which needs nbgp ghost points and
	975	!-- grid levels for the multigrid scheme. In the latter case only one
	976	!-- ghost point is necessary.
[709]	977	!-- First definition of MPI-datatypes for exchange of ghost layers on normal
[667]	978	!-- grid. The following loop is needed for data exchange in poismg.f90.
	979	!
	980	!-- Determine number of grid points of yz-layer for exchange
	981	ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
[709]	982
[667]	983	!
[709]	984	!-- Define an MPI-datatype for the exchange of left/right boundaries.
	985	!-- Although data are contiguous in physical memory (which does not
	986	!-- necessarily require an MPI-derived datatype), the data exchange between
	987	!-- left and right PE's using the MPI-derived type is 10% faster than without.
[667]	988	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz(0), &
[709]	989	MPI_REAL, type_xz(0), ierr )
[667]	990	CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
[1]	991
[709]	992	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
	993	ierr )
[667]	994	CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
[709]	995
[667]	996	!
[709]	997	!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
[667]	998	IF ( psolver == 'multigrid' ) THEN
	999	!
[709]	1000	!-- Definition of MPI-datatyoe as above, but only 1 ghost level is used
	1001	DO i = maximum_grid_level, 1 , -1
	1002
[667]	1003	ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	1004
	1005	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
[709]	1006	MPI_REAL, type_xz(i), ierr )
[667]	1007	CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
[1]	1008
[709]	1009	CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
	1010	ierr )
[667]	1011	CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
	1012
	1013	nxl_l = nxl_l / 2
	1014	nxr_l = nxr_l / 2
	1015	nys_l = nys_l / 2
	1016	nyn_l = nyn_l / 2
	1017	nzt_l = nzt_l / 2
[709]	1018
[667]	1019	ENDDO
[709]	1020
	1021	ENDIF
[1]	1022	#endif
	1023
[809]	1024	#if defined( __parallel ) && ! defined ( __check )
[1]	1025	!
	1026	!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
[106]	1027	!-- horizontal boundary conditions.
[1]	1028	IF ( pleft == MPI_PROC_NULL ) THEN
[1159]	1029	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	1030	inflow_l = .TRUE.
[1159]	1031	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	1032	outflow_l = .TRUE.
	1033	ENDIF
	1034	ENDIF
	1035
	1036	IF ( pright == MPI_PROC_NULL ) THEN
[1159]	1037	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	1038	outflow_r = .TRUE.
[1159]	1039	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	1040	inflow_r = .TRUE.
	1041	ENDIF
	1042	ENDIF
	1043
	1044	IF ( psouth == MPI_PROC_NULL ) THEN
[1159]	1045	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	1046	outflow_s = .TRUE.
[1159]	1047	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	1048	inflow_s = .TRUE.
	1049	ENDIF
	1050	ENDIF
	1051
	1052	IF ( pnorth == MPI_PROC_NULL ) THEN
[1159]	1053	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	1054	inflow_n = .TRUE.
[1159]	1055	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	1056	outflow_n = .TRUE.
	1057	ENDIF
	1058	ENDIF
	1059
[151]	1060	!
	1061	!-- Broadcast the id of the inflow PE
	1062	IF ( inflow_l ) THEN
[163]	1063	id_inflow_l = myidx
[151]	1064	ELSE
	1065	id_inflow_l = 0
	1066	ENDIF
[622]	1067	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[151]	1068	CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
	1069	comm1dx, ierr )
	1070
[163]	1071	!
	1072	!-- Broadcast the id of the recycling plane
	1073	!-- WARNING: needs to be adjusted in case of inflows other than from left side!
[1139]	1074	IF ( NINT( recycling_width / dx ) >= nxl .AND. &
	1075	NINT( recycling_width / dx ) <= nxr ) THEN
[163]	1076	id_recycling_l = myidx
	1077	ELSE
	1078	id_recycling_l = 0
	1079	ENDIF
[622]	1080	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[163]	1081	CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
	1082	comm1dx, ierr )
	1083
[1402]	1084	CALL location_message( 'finished', .TRUE. )
[1384]	1085
[809]	1086	#elif ! defined ( __parallel )
[1159]	1087	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	1088	inflow_l = .TRUE.
	1089	outflow_r = .TRUE.
[1159]	1090	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	1091	outflow_l = .TRUE.
	1092	inflow_r = .TRUE.
	1093	ENDIF
	1094
[1159]	1095	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	1096	inflow_n = .TRUE.
	1097	outflow_s = .TRUE.
[1159]	1098	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	1099	outflow_n = .TRUE.
	1100	inflow_s = .TRUE.
	1101	ENDIF
	1102	#endif
[807]	1103
[106]	1104	!
[978]	1105	!-- At the inflow or outflow, u or v, respectively, have to be calculated for
	1106	!-- one more grid point.
	1107	IF ( inflow_l .OR. outflow_l ) THEN
[106]	1108	nxlu = nxl + 1
	1109	ELSE
	1110	nxlu = nxl
	1111	ENDIF
[978]	1112	IF ( inflow_s .OR. outflow_s ) THEN
[106]	1113	nysv = nys + 1
	1114	ELSE
	1115	nysv = nys
	1116	ENDIF
[1]	1117
[114]	1118	!
	1119	!-- Allocate wall flag arrays used in the multigrid solver
	1120	IF ( psolver == 'multigrid' ) THEN
	1121
	1122	DO i = maximum_grid_level, 1, -1
	1123
	1124	SELECT CASE ( i )
	1125
	1126	CASE ( 1 )
	1127	ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1, &
	1128	nys_mg(i)-1:nyn_mg(i)+1, &
	1129	nxl_mg(i)-1:nxr_mg(i)+1) )
	1130
	1131	CASE ( 2 )
	1132	ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1, &
	1133	nys_mg(i)-1:nyn_mg(i)+1, &
	1134	nxl_mg(i)-1:nxr_mg(i)+1) )
	1135
	1136	CASE ( 3 )
	1137	ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1, &
	1138	nys_mg(i)-1:nyn_mg(i)+1, &
	1139	nxl_mg(i)-1:nxr_mg(i)+1) )
	1140
	1141	CASE ( 4 )
	1142	ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1, &
	1143	nys_mg(i)-1:nyn_mg(i)+1, &
	1144	nxl_mg(i)-1:nxr_mg(i)+1) )
	1145
	1146	CASE ( 5 )
	1147	ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1, &
	1148	nys_mg(i)-1:nyn_mg(i)+1, &
	1149	nxl_mg(i)-1:nxr_mg(i)+1) )
	1150
	1151	CASE ( 6 )
	1152	ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1, &
	1153	nys_mg(i)-1:nyn_mg(i)+1, &
	1154	nxl_mg(i)-1:nxr_mg(i)+1) )
	1155
	1156	CASE ( 7 )
	1157	ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1, &
	1158	nys_mg(i)-1:nyn_mg(i)+1, &
	1159	nxl_mg(i)-1:nxr_mg(i)+1) )
	1160
	1161	CASE ( 8 )
	1162	ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1, &
	1163	nys_mg(i)-1:nyn_mg(i)+1, &
	1164	nxl_mg(i)-1:nxr_mg(i)+1) )
	1165
	1166	CASE ( 9 )
	1167	ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1, &
	1168	nys_mg(i)-1:nyn_mg(i)+1, &
	1169	nxl_mg(i)-1:nxr_mg(i)+1) )
	1170
	1171	CASE ( 10 )
	1172	ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1, &
	1173	nys_mg(i)-1:nyn_mg(i)+1, &
	1174	nxl_mg(i)-1:nxr_mg(i)+1) )
	1175
	1176	CASE DEFAULT
[254]	1177	message_string = 'more than 10 multigrid levels'
	1178	CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
[114]	1179
	1180	END SELECT
	1181
	1182	ENDDO
	1183
	1184	ENDIF
	1185
[759]	1186	!
	1187	!-- Calculate the number of groups into which parallel I/O is split.
	1188	!-- The default for files which are opened by all PEs (or where each
	1189	!-- PE opens his own independent file) is, that all PEs are doing input/output
	1190	!-- in parallel at the same time. This might cause performance or even more
	1191	!-- severe problems depending on the configuration of the underlying file
	1192	!-- system.
	1193	!-- First, set the default:
	1194	IF ( maximum_parallel_io_streams == -1 .OR. &
	1195	maximum_parallel_io_streams > numprocs ) THEN
	1196	maximum_parallel_io_streams = numprocs
	1197	ENDIF
	1198
	1199	!
	1200	!-- Now calculate the number of io_blocks and the io_group to which the
	1201	!-- respective PE belongs. I/O of the groups is done in serial, but in parallel
	1202	!-- for all PEs belonging to the same group. A preliminary setting with myid
	1203	!-- based on MPI_COMM_WORLD has been done in parin.
	1204	io_blocks = numprocs / maximum_parallel_io_streams
	1205	io_group = MOD( myid+1, io_blocks )
	1206
	1207
[1]	1208	END SUBROUTINE init_pegrid

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |