Home

Context Navigation

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

Last change on this file since 4656 was 4648, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 53.4 KB

Rev	Line
[1682]	1	!> @file init_pegrid.f90
[4648]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[4648]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1036]	8	!
[4648]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1036]	12	!
[4648]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1036]	15	!
[254]	16	! Current revisions:
[1322]	17	! ------------------
[4648]	18	!
	19	!
[1321]	20	! Former revisions:
	21	! -----------------
	22	! $Id: init_pegrid.f90 4648 2020-08-25 07:52:08Z pavelkrc $
[4648]	23	! file re-formatted to follow the PALM coding standard
	24	!
	25	! 4564 2020-06-12 14:03:36Z raasch
[4564]	26	! Vertical nesting method of Huq et al. (2019) removed
[4648]	27	!
[4564]	28	! 4461 2020-03-12 16:51:59Z raasch
[4461]	29	! communicator configurations for four virtual pe grids defined
[4648]	30	!
[4461]	31	! 4444 2020-03-05 15:59:50Z raasch
[4444]	32	! bugfix: cpp-directives for serial mode added
[4648]	33	!
[4444]	34	! 4360 2020-01-07 11:25:50Z suehring
[4301]	35	! changed message PA0467
[4648]	36	!
[4301]	37	! 4264 2019-10-15 16:00:23Z scharf
[4264]	38	! corrected error message string
[4648]	39	!
[4264]	40	! 4241 2019-09-27 06:32:47Z raasch
[4648]	41	! Check added to ensure that subdomain grid has at least the size as given by the number of ghost
	42	! points
	43	!
[4241]	44	! 4182 2019-08-22 15:20:23Z scharf
[4182]	45	! Corrected "Former revisions" section
[4648]	46	!
[4182]	47	! 4045 2019-06-21 10:58:47Z raasch
[4648]	48	! bugfix: kind attribute added to nint function to allow for large integers which may appear in case
	49	! of default recycling width and small grid spacings
	50	!
[4045]	51	! 3999 2019-05-23 16:09:37Z suehring
[3999]	52	! Spend 3 ghost points also in case of pw-scheme when nesting is applied
[4648]	53	!
[3999]	54	! 3897 2019-04-15 11:51:14Z suehring
[3897]	55	! Minor revision of multigrid check; give warning instead of an abort.
[4648]	56	!
[3897]	57	! 3890 2019-04-12 15:59:20Z suehring
[4648]	58	! Check if grid coarsening is possible on subdomain, in order to avoid that multigrid approach
	59	! effectively reduces to a Gauss-Seidel scheme.
	60	!
[3890]	61	! 3885 2019-04-11 11:29:34Z kanani
[4648]	62	! Changes related to global restructuring of location messages and introduction of additional debug
	63	! messages
	64	!
[3885]	65	! 3884 2019-04-10 13:31:55Z Giersch
[3884]	66	! id_recycling is only calculated in case of tubulent inflow
[4648]	67	!
[3884]	68	! 3761 2019-02-25 15:31:42Z raasch
[3761]	69	! unused variable removed
[4648]	70	!
[3761]	71	! 3655 2019-01-07 16:51:22Z knoop
[3553]	72	! variables documented
[2259]	73	!
[4182]	74	! Revision 1.1 1997/07/24 11:15:09 raasch
	75	! Initial revision
	76	!
	77	!
[1]	78	! Description:
	79	! ------------
[4648]	80	!> Determination of the virtual processor topology (if not prescribed by the user) and computation
	81	!> of the grid point number and array bounds of the local domains.
	82	!> @todo: remove MPI-data types for 2D exchange on coarse multigrid level (not used any more)
	83	!--------------------------------------------------------------------------------------------------!
[1682]	84	SUBROUTINE init_pegrid
[1]	85
	86
[4648]	87	USE control_parameters, &
	88	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, bc_lr, bc_ns, &
	89	bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s, grid_level, &
	90	grid_level_count, maximum_grid_level, message_string, mg_switch_to_pe0_level, psolver
[4444]	91
[4648]	92
[4444]	93	#if defined( __parallel )
[4648]	94	USE control_parameters, &
	95	ONLY: coupling_mode, coupling_topology, gathered_size, momentum_advec, &
	96	outflow_source_plane, recycling_width, scalar_advec, subdomain_size, &
[4444]	97	turbulent_inflow, turbulent_outflow, y_shift
	98
[4648]	99	USE grid_variables, &
[1320]	100	ONLY: dx
[4444]	101	#endif
[1]	102
[4648]	103	USE indices, &
	104	ONLY: nnx, nny, nnz, nx, nxl, nxl_mg, nxlu, nxr, nxr_mg, ny, nyn, nyn_mg, nys, nys_mg, &
	105	nysv, nz, nzb, nzt, nzt_mg, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4, &
	106	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, wall_flags_9, wall_flags_10
	107
[4444]	108	#if defined( __parallel )
[4648]	109	USE indices, &
[4444]	110	ONLY: mg_loc_ind, nbgp, nx_a, nx_o, ny_a, ny_o
	111	#endif
	112
[1320]	113	USE kinds
[4648]	114
[1320]	115	USE pegrid
[4648]	116
[4444]	117	#if defined( __parallel )
[4648]	118	USE pmc_interface, &
[3999]	119	ONLY: nested_run
[4444]	120
[4648]	121	USE spectra_mod, &
[4444]	122	ONLY: calculate_spectra
[1833]	123
[4648]	124	USE synthetic_turbulence_generator_mod, &
	125	ONLY: id_stg_left, id_stg_north, id_stg_right, id_stg_south, use_syn_turb_gen
[4444]	126	#endif
[2259]	127
[4648]	128	USE transpose_indices, &
	129	ONLY: nxl_y, nxl_z, nxr_y, nxr_z, nyn_x, nyn_z, nys_x, nys_z, nzb_x, nzb_y, nzt_x, nzt_y
[667]	130
[4444]	131	#if defined( __parallel )
[4648]	132	USE transpose_indices, &
[4444]	133	ONLY: nxl_yd, nxr_yd, nzb_yd, nzt_yd
	134	#endif
[2365]	135
[1]	136	IMPLICIT NONE
	137
[3552]	138	INTEGER(iwp) :: i !< running index over number of processors or number of multigrid level
[4444]	139	#if defined( __parallel )
[3552]	140	INTEGER(iwp) :: id_inflow_l !< ID indicating processors located at the left inflow boundary
[2050]	141	INTEGER(iwp) :: id_outflow_l !< local value of id_outflow
	142	INTEGER(iwp) :: id_outflow_source_l !< local value of id_outflow_source
[3552]	143	INTEGER(iwp) :: id_recycling_l !< ID indicating processors located at the recycling plane
[4648]	144	INTEGER(iwp) :: id_stg_left_l !< left lateral boundary local core id in case of turbulence generator
	145	INTEGER(iwp) :: id_stg_north_l !< north lateral boundary local core id in case of turbulence generator
	146	INTEGER(iwp) :: id_stg_right_l !< right lateral boundary local core id in case of turbulence generator
	147	INTEGER(iwp) :: id_stg_south_l !< south lateral boundary local core id in case of turbulence generator
[3552]	148	INTEGER(iwp) :: ind(5) !< array containing the subdomain bounds
[4444]	149	#endif
[3552]	150	INTEGER(iwp) :: j !< running index, used for various loops
	151	INTEGER(iwp) :: k !< number of vertical grid points in different multigrid level
	152	INTEGER(iwp) :: maximum_grid_level_l !< maximum number of grid level without switching to PE 0
	153	INTEGER(iwp) :: mg_levels_x !< maximum number of grid level allowed along x-direction
	154	INTEGER(iwp) :: mg_levels_y !< maximum number of grid level allowed along y-direction
	155	INTEGER(iwp) :: mg_levels_z !< maximum number of grid level allowed along z-direction
	156	INTEGER(iwp) :: mg_switch_to_pe0_level_l !< maximum number of grid level with switching to PE 0
[4444]	157	#if defined( __parallel )
[3552]	158	INTEGER(iwp) :: nnx_y !< quotient of number of grid points along x-direction and number of PEs used along y-direction
	159	INTEGER(iwp) :: nny_x !< quotient of number of grid points along y-direction and number of PEs used along x-direction
	160	INTEGER(iwp) :: nny_z !< quotient of number of grid points along y-direction and number of PEs used along x-direction
	161	INTEGER(iwp) :: nnz_x !< quotient of number of grid points along z-direction and number of PEs used along x-direction
	162	INTEGER(iwp) :: nnz_y !< quotient of number of grid points along z-direction and number of PEs used along x-direction
	163	INTEGER(iwp) :: numproc_sqr !< square root of the number of processors
[4444]	164	#endif
[3552]	165	INTEGER(iwp) :: nxl_l !< lower index bound along x-direction on subdomain and different multigrid level
	166	INTEGER(iwp) :: nxr_l !< upper index bound along x-direction on subdomain and different multigrid level
	167	INTEGER(iwp) :: nyn_l !< lower index bound along y-direction on subdomain and different multigrid level
	168	INTEGER(iwp) :: nys_l !< upper index bound along y-direction on subdomain and different multigrid level
[4444]	169	#if defined( __parallel )
[3552]	170	INTEGER(iwp) :: nzb_l !< lower index bound along z-direction on subdomain and different multigrid level
[4444]	171	#endif
[3552]	172	INTEGER(iwp) :: nzt_l !< upper index bound along z-direction on subdomain and different multigrid level
	173	!$ INTEGER(iwp) :: omp_get_num_threads !< number of OpenMP threads
[1]	174
[4444]	175	#if defined( __parallel )
[3552]	176	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ind_all !< dummy array containing index bounds on subdomain, used for gathering
	177	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxlf !< lower index bound allong x-direction for every PE
	178	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxrf !< upper index bound allong x-direction for every PE
	179	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nynf !< lower index bound allong y-direction for every PE
	180	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nysf !< lower index bound allong y-direction for every PE
[1]	181
[3552]	182	INTEGER(iwp), DIMENSION(2) :: pdims_remote !< number of PEs used for coupled model (only in atmospher-ocean coupling)
[2372]	183	INTEGER(iwp) :: lcoord(2) !< PE coordinates of left neighbor along x and y
	184	INTEGER(iwp) :: rcoord(2) !< PE coordinates of right neighbor along x and y
[4444]	185	#endif
[667]	186
[1]	187	!
	188	!-- Get the number of OpenMP threads
	189	!$OMP PARALLEL
	190	!$ threads_per_task = omp_get_num_threads()
	191	!$OMP END PARALLEL
	192
	193
	194	#if defined( __parallel )
[667]	195
[3885]	196	CALL location_message( 'creating virtual PE grids + MPI derived data types', 'start' )
[1764]	197
[1]	198	!
[2177]	199	!-- Determine the processor topology or check it, if prescribed by the user
[1779]	200	IF ( npex == -1 .AND. npey == -1 ) THEN
[1]	201
	202	!
[2177]	203	!-- Automatic determination of the topology
[1779]	204	numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
	205	pdims(1) = MAX( numproc_sqr , 1 )
[2180]	206	DO WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
[1779]	207	pdims(1) = pdims(1) - 1
	208	ENDDO
[2180]	209	pdims(2) = numprocs / pdims(1)
[1]	210
[1779]	211	ELSEIF ( npex /= -1 .AND. npey /= -1 ) THEN
[1]	212
	213	!
[4648]	214	!-- Prescribed by user. Number of processors on the prescribed topology must be equal to the
	215	!-- number of PEs available to the job
[1779]	216	IF ( ( npex * npey ) /= numprocs ) THEN
[4648]	217	WRITE( message_string, * ) 'number of PEs of the prescribed ', 'topology (', npex*npey, &
	218	') does not match & the number of ', &
	219	'PEs available to the job (', numprocs, ')'
[1779]	220	CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
	221	ENDIF
	222	pdims(1) = npex
	223	pdims(2) = npey
[1]	224
[1779]	225	ELSE
[1]	226	!
[1779]	227	!-- If the processor topology is prescribed by the user, the number of
	228	!-- PEs must be given in both directions
[4648]	229	message_string = 'if the processor topology is prescribed by th' // &
	230	'e user & both values of "npex" and "npey" must be given' // &
	231	' in the &NAMELIST-parameter file'
[1779]	232	CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
[1]	233
	234	ENDIF
	235
	236	!
[4648]	237	!-- Create four default MPI communicators for the 2d virtual PE grid. One of them will be used as
	238	!-- the main communicator for this run, while others might be used for specific quantities like
[4461]	239	!-- aerosol, chemical species, or passive scalars), if their horizontal boundary conditions shall
	240	!-- be different from those of the other quantities (e.g. non-cyclic conditions for aerosols, and
	241	!-- cyclic conditions for all others).
	242	DO i = 1, 4
	243
	244	IF ( i == 1 ) cyclic = (/ .TRUE., .TRUE. /) ! cyclic along x and y
	245	IF ( i == 2 ) cyclic = (/ .TRUE., .FALSE. /) ! cyclic along x
	246	IF ( i == 3 ) cyclic = (/ .FALSE., .TRUE. /) ! cyllic along y
	247	IF ( i == 4 ) cyclic = (/ .FALSE., .FALSE. /) ! non-cyclic
	248
	249	CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
	250	communicator_configurations(i)%mpi_communicator, ierr )
	251
	252	CALL MPI_CART_SHIFT( communicator_configurations(i)%mpi_communicator, 0, 1, &
	253	communicator_configurations(i)%pleft, &
	254	communicator_configurations(i)%pright, ierr )
	255
	256	CALL MPI_CART_SHIFT( communicator_configurations(i)%mpi_communicator, 1, 1, &
	257	communicator_configurations(i)%psouth, &
	258	communicator_configurations(i)%pnorth, ierr )
	259
	260	ENDDO
	261
	262	!
[1]	263	!-- If necessary, set horizontal boundary conditions to non-cyclic
[722]	264	IF ( bc_lr /= 'cyclic' ) cyclic(1) = .FALSE.
	265	IF ( bc_ns /= 'cyclic' ) cyclic(2) = .FALSE.
[1]	266
[807]	267
[1]	268	!
[4461]	269	!-- Set the main communicator (virtual pe grid) for this run
	270	IF ( bc_lr == 'cyclic' .AND. bc_ns == 'cyclic' ) i = 1
	271	IF ( bc_lr == 'cyclic' .AND. bc_ns /= 'cyclic' ) i = 2
	272	IF ( bc_lr /= 'cyclic' .AND. bc_ns == 'cyclic' ) i = 3
	273	IF ( bc_lr /= 'cyclic' .AND. bc_ns /= 'cyclic' ) i = 4
	274
	275	comm2d = communicator_configurations(i)%mpi_communicator
	276	pleft = communicator_configurations(i)%pleft
	277	pright = communicator_configurations(i)%pright
	278	psouth = communicator_configurations(i)%psouth
	279	pnorth = communicator_configurations(i)%pnorth
	280
	281	!
	282	!-- Set rank and coordinates of the main communicator
[1]	283	CALL MPI_COMM_RANK( comm2d, myid, ierr )
[1468]	284	WRITE (myid_char,'(''_'',I6.6)') myid
[1]	285
	286	CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
[4461]	287
[1]	288	!
[4648]	289	!-- In case of cyclic boundary conditions, a y-shift at the boundaries in x-direction can be
	290	!-- introduced via parameter y_shift. The shift is done by modifying the processor grid in such a
	291	!-- way that processors located at the x-boundary communicate across it to processors with
	292	!-- y-coordinate shifted by y_shift relative to their own. This feature can not be used in
	293	!-- combination with an fft pressure solver. It has been implemented to counter the effect of streak
	294	!-- structures in case of cyclic boundary conditions. For a description of these see Munters
[2372]	295	!-- (2016; dx.doi.org/10.1063/1.4941912)
	296	!--
	297	!-- Get coordinates of left and right neighbor on PE grid
	298	IF ( y_shift /= 0 ) THEN
[4301]	299	IF ( bc_lr == 'cyclic' ) THEN
[4648]	300	IF ( TRIM( psolver ) /= 'multigrid' .AND. TRIM( psolver ) /= 'multigrid_noopt') THEN
[4301]	301	message_string = 'y_shift /= 0 requires a multigrid pressure solver '
	302	CALL message( 'check_parameters', 'PA0468', 1, 2, 0, 6, 0 )
	303	ENDIF
[2372]	304
[4301]	305	CALL MPI_CART_COORDS( comm2d, pright, ndim, rcoord, ierr )
	306	CALL MPI_CART_COORDS( comm2d, pleft, ndim, lcoord, ierr )
[2372]	307
[4648]	308	!
	309	!-- If the x(y)-coordinate of the right (left) neighbor is smaller (greater) than that of the
	310	!-- calling process, then the calling process is located on the right (left) boundary of the
	311	!-- processor grid. In that case, the y-coordinate of that neighbor is increased (decreased)
	312	!-- by y_shift.
	313	!-- The rank of the process with that coordinate is then inquired and the neighbor rank for
	314	!-- MPI_SENDRECV, pright (pleft) is set to it.
	315	!-- In this way, the calling process receives a new right (left) neighbor for all future
	316	!-- MPI_SENDRECV calls. That neighbor has a y-coordinate of y+(-)y_shift, where y is the
	317	!-- original right (left) neighbor's y-coordinate. The modulo-operation ensures that if the
	318	!-- neighbor's y-coordinate exceeds the grid-boundary, it will be relocated to the opposite
	319	!-- part of the grid cyclicly.
	320	IF ( rcoord(1) < pcoord(1) ) THEN
[4301]	321	rcoord(2) = MODULO( rcoord(2) + y_shift, pdims(2) )
	322	CALL MPI_CART_RANK( comm2d, rcoord, pright, ierr )
	323	ENDIF
[2372]	324
[4648]	325	IF ( lcoord(1) > pcoord(1) ) THEN
[4301]	326	lcoord(2) = MODULO( lcoord(2) - y_shift, pdims(2) )
	327	CALL MPI_CART_RANK( comm2d, lcoord, pleft, ierr )
	328	ENDIF
[4648]	329
[4301]	330	ELSE
[2372]	331	!
[4648]	332	!-- y-shift for non-cyclic boundary conditions is only implemented
[4301]	333	!-- for the turbulence recycling method in inflow_turbulence.f90
	334	IF ( .NOT. turbulent_inflow ) THEN
[4648]	335	message_string = 'y_shift /= 0 is only allowed for cyclic ' // &
	336	'boundary conditions in both directions ' // &
[4301]	337	'or with turbulent_inflow == .TRUE.'
	338	CALL message( 'check_parameters', 'PA0467', 1, 2, 0, 6, 0 )
	339	ENDIF
[2372]	340	ENDIF
	341	ENDIF
[4461]	342
[2372]	343	!
[1]	344	!-- Determine sub-topologies for transpositions
	345	!-- Transposition from z to x:
	346	remain_dims(1) = .TRUE.
	347	remain_dims(2) = .FALSE.
	348	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
	349	CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
	350	!
	351	!-- Transposition from x to y
	352	remain_dims(1) = .FALSE.
	353	remain_dims(2) = .TRUE.
	354	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
	355	CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
	356
	357
	358	!
[1003]	359	!-- Calculate array bounds along x-direction for every PE.
[4648]	360	ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), nysf(0:pdims(2)-1) )
[1]	361
[1003]	362	IF ( MOD( nx+1 , pdims(1) ) /= 0 ) THEN
[4648]	363	WRITE( message_string, * ) 'x-direction: gridpoint number (' ,nx+1, ') ', &
	364	'is not an& integral multiple of the number', &
	365	' of processors (', pdims(1), ')'
[254]	366	CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
[1]	367	ELSE
[1003]	368	nnx = ( nx + 1 ) / pdims(1)
[4264]	369	ENDIF
[1]	370
	371	!
	372	!-- Left and right array bounds, number of gridpoints
	373	DO i = 0, pdims(1)-1
	374	nxlf(i) = i * nnx
	375	nxrf(i) = ( i + 1 ) * nnx - 1
	376	ENDDO
	377
	378	!
	379	!-- Calculate array bounds in y-direction for every PE.
[1003]	380	IF ( MOD( ny+1 , pdims(2) ) /= 0 ) THEN
[4648]	381	WRITE( message_string, * ) 'y-direction: gridpoint number (', ny+1, ') ', &
	382	'is not an& integral multiple of the number', &
	383	' of processors (', pdims(2), ')'
[254]	384	CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
[1]	385	ELSE
[1003]	386	nny = ( ny + 1 ) / pdims(2)
[4264]	387	ENDIF
[1]	388
	389	!
	390	!-- South and north array bounds
	391	DO j = 0, pdims(2)-1
	392	nysf(j) = j * nny
	393	nynf(j) = ( j + 1 ) * nny - 1
	394	ENDDO
	395
	396	!
	397	!-- Local array bounds of the respective PEs
[1003]	398	nxl = nxlf(pcoord(1))
	399	nxr = nxrf(pcoord(1))
	400	nys = nysf(pcoord(2))
	401	nyn = nynf(pcoord(2))
	402	nzb = 0
	403	nzt = nz
	404	nnz = nz
[1]	405
	406	!
[4648]	407	!-- Set switches to define if the PE is situated at the border of the virtual processor grid
[707]	408	IF ( nxl == 0 ) left_border_pe = .TRUE.
	409	IF ( nxr == nx ) right_border_pe = .TRUE.
	410	IF ( nys == 0 ) south_border_pe = .TRUE.
	411	IF ( nyn == ny ) north_border_pe = .TRUE.
	412
	413	!
[4648]	414	!-- Calculate array bounds and gridpoint numbers for the transposed arrays (needed in the pressure
	415	!-- solver)
	416	!-- For the transposed arrays, cyclic boundaries as well as top and bottom boundaries are omitted,
	417	!-- because they are obstructive to the transposition
[1]	418
	419	!
	420	!-- 1. transposition z --> x
[1001]	421	!-- This transposition is not neccessary in case of a 1d-decomposition along x
[2938]	422	IF ( psolver == 'poisfft' .OR. calculate_spectra ) THEN
[1304]	423
[1922]	424	IF ( pdims(2) /= 1 ) THEN
	425	IF ( MOD( nz , pdims(1) ) /= 0 ) THEN
[4648]	426	WRITE( message_string, * ) 'transposition z --> x:& ', &
	427	'nz=', nz, ' is not an integral multiple ', &
	428	'of pdims(1)=', pdims(1)
[1922]	429	CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
	430	ENDIF
[1]	431	ENDIF
[1922]	432
	433	nys_x = nys
	434	nyn_x = nyn
	435	nny_x = nny
	436	nnz_x = nz / pdims(1)
	437	nzb_x = 1 + myidx * nnz_x
	438	nzt_x = ( myidx + 1 ) * nnz_x
	439	sendrecvcount_zx = nnx * nny * nnz_x
	440
[1]	441	ENDIF
	442
[1922]	443
[4264]	444	IF ( psolver == 'poisfft' ) THEN
[1]	445	!
[1922]	446	!-- 2. transposition x --> y
	447	IF ( MOD( nx+1 , pdims(2) ) /= 0 ) THEN
[4648]	448	WRITE( message_string, * ) 'transposition x --> y:& ', &
	449	'nx+1=', nx+1, ' is not an integral ', &
	450	'multiple of pdims(2)=', pdims(2)
[1922]	451	CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
	452	ENDIF
[1]	453
[1922]	454	nnz_y = nnz_x
	455	nzb_y = nzb_x
	456	nzt_y = nzt_x
	457	nnx_y = (nx+1) / pdims(2)
	458	nxl_y = myidy * nnx_y
	459	nxr_y = ( myidy + 1 ) * nnx_y - 1
	460	sendrecvcount_xy = nnx_y * nny_x * nnz_y
[1]	461	!
[4264]	462	!-- 3. transposition y --> z
[1922]	463	!-- (ELSE: x --> y in case of 1D-decomposition along x)
	464	nxl_z = nxl_y
	465	nxr_z = nxr_y
	466	nny_z = (ny+1) / pdims(1)
	467	nys_z = myidx * nny_z
	468	nyn_z = ( myidx + 1 ) * nny_z - 1
	469	sendrecvcount_yz = nnx_y * nny_z * nnz_y
[1304]	470
[1922]	471	IF ( pdims(2) /= 1 ) THEN
[1]	472	!
[1922]	473	!-- y --> z
	474	!-- This transposition is not neccessary in case of a 1d-decomposition
	475	!-- along x, except that the uptream-spline method is switched on
	476	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[4648]	477	WRITE( message_string, * ) 'transposition y --> z:& ', &
	478	'ny+1=', ny+1, ' is not an integral ', &
	479	'multiple of pdims(1)=', pdims(1)
[1922]	480	CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
	481	ENDIF
[1]	482
[1922]	483	ELSE
[1]	484	!
[1922]	485	!-- x --> y
	486	!-- This condition must be fulfilled for a 1D-decomposition along x
	487	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[4648]	488	WRITE( message_string, * ) 'transposition x --> y:& ', &
	489	'ny+1=', ny+1, ' is not an integral ', &
	490	'multiple of pdims(1)=', pdims(1)
[1922]	491	CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
	492	ENDIF
	493
[1]	494	ENDIF
	495
	496	ENDIF
	497
	498	!
	499	!-- Indices for direct transpositions z --> y (used for calculating spectra)
[1922]	500	IF ( calculate_spectra ) THEN
[1003]	501	IF ( MOD( nz, pdims(2) ) /= 0 ) THEN
[4648]	502	WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
	503	'for spectra):& nz=', nz, ' is not an ', &
	504	'integral multiple of pdims(2)=', pdims(2)
[254]	505	CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
[1]	506	ELSE
[1003]	507	nxl_yd = nxl
	508	nxr_yd = nxr
	509	nzb_yd = 1 + myidy * ( nz / pdims(2) )
	510	nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
	511	sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
[1]	512	ENDIF
	513	ENDIF
	514
[1922]	515	IF ( psolver == 'poisfft' .OR. calculate_spectra ) THEN
[1]	516	!
[4648]	517	!-- Indices for direct transpositions y --> x
[1922]	518	!-- (they are only possible in case of a 1d-decomposition along x)
	519	IF ( pdims(2) == 1 ) THEN
	520	nny_x = nny / pdims(1)
	521	nys_x = myid * nny_x
	522	nyn_x = ( myid + 1 ) * nny_x - 1
	523	nzb_x = 1
	524	nzt_x = nz
	525	sendrecvcount_xy = nnx * nny_x * nz
	526	ENDIF
	527
[1]	528	ENDIF
	529
[1922]	530	IF ( psolver == 'poisfft' ) THEN
[1]	531	!
[4648]	532	!-- Indices for direct transpositions x --> y
[1922]	533	!-- (they are only possible in case of a 1d-decomposition along y)
	534	IF ( pdims(1) == 1 ) THEN
	535	nnx_y = nnx / pdims(2)
	536	nxl_y = myid * nnx_y
	537	nxr_y = ( myid + 1 ) * nnx_y - 1
	538	nzb_y = 1
	539	nzt_y = nz
	540	sendrecvcount_xy = nnx_y * nny * nz
	541	ENDIF
	542
[1]	543	ENDIF
	544
	545	!
	546	!-- Arrays for storing the array bounds are needed any more
	547	DEALLOCATE( nxlf , nxrf , nynf , nysf )
	548
[807]	549
[145]	550	!
	551	!-- Collect index bounds from other PEs (to be written to restart file later)
	552	ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
	553
	554	IF ( myid == 0 ) THEN
	555
	556	hor_index_bounds(1,0) = nxl
	557	hor_index_bounds(2,0) = nxr
	558	hor_index_bounds(3,0) = nys
	559	hor_index_bounds(4,0) = nyn
	560
	561	!
	562	!-- Receive data from all other PEs
	563	DO i = 1, numprocs-1
[4648]	564	CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, ierr )
[145]	565	hor_index_bounds(:,i) = ibuf(1:4)
	566	ENDDO
	567
	568	ELSE
	569	!
	570	!-- Send index bounds to PE0
	571	ibuf(1) = nxl
	572	ibuf(2) = nxr
	573	ibuf(3) = nys
	574	ibuf(4) = nyn
	575	CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
	576
	577	ENDIF
	578
[807]	579
[1]	580	#if defined( __print )
	581	!
	582	!-- Control output
	583	IF ( myid == 0 ) THEN
	584	PRINT, ' processor topology *'
	585	PRINT*, ' '
[4648]	586	PRINT*, 'myid pcoord left right south north idx idy nxl: nxr',' nys: nyn'
	587	PRINT*, '------------------------------------------------------------','-----------'
	588	WRITE (*,1000) 0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, myidx, myidy, nxl, &
	589	nxr, nys, nyn
	590	1000 FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3,2(2X,I4,':',I4))
[1]	591
	592	!
[108]	593	!-- Receive data from the other PEs
[1]	594	DO i = 1,numprocs-1
[4648]	595	CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, ierr )
[1]	596	WRITE (*,1000) i, ( ibuf(j) , j = 1,12 )
	597	ENDDO
	598	ELSE
	599
	600	!
	601	!-- Send data to PE0
	602	ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
	603	ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
	604	ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
	605	ibuf(12) = nyn
[4648]	606	CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )
[1]	607	ENDIF
	608	#endif
	609
[4648]	610	!
[709]	611	!-- Determine the number of ghost point layers
[4648]	612	IF ( scalar_advec == 'ws-scheme' .OR. &
[3999]	613	momentum_advec == 'ws-scheme' .OR. nested_run ) THEN
[667]	614	nbgp = 3
	615	ELSE
	616	nbgp = 1
[4648]	617	ENDIF
[667]	618
[102]	619	!
[4648]	620	!-- Check that the number of computational grid points is not smaller than the number of ghost
	621	!-- points.
[4241]	622	IF ( nnx < nbgp ) THEN
	623	WRITE( message_string, * ) 'number of subdomain grid points along x (', nnx, ') is smaller',&
	624	'than the number of ghost points (', nbgp, ')'
	625	CALL message( 'init_pegrid', 'PA0682', 1, 2, 0, 6, 0 )
	626	ENDIF
	627	IF ( nny < nbgp ) THEN
	628	WRITE( message_string, * ) 'number of subdomain grid points along y (', nny, ') is smaller',&
	629	'than the number of ghost points (', nbgp, ')'
	630	CALL message( 'init_pegrid', 'PA0683', 1, 2, 0, 6, 0 )
	631	ENDIF
	632
	633	!
[4648]	634	!-- Create a new MPI derived datatype for the exchange of surface (xy) data, which is needed for
	635	!-- coupled atmosphere-ocean runs.
[709]	636	!-- First, calculate number of grid points of an xy-plane.
[667]	637	ngp_xy = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
[102]	638	CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
	639	CALL MPI_TYPE_COMMIT( type_xy, ierr )
[667]	640
[4564]	641	IF ( TRIM( coupling_mode ) /= 'uncoupled' ) THEN
[4648]	642
[667]	643	!
	644	!-- Pass the number of grid points of the atmosphere model to
	645	!-- the ocean model and vice versa
	646	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	647
	648	nx_a = nx
	649	ny_a = ny
	650
[709]	651	IF ( myid == 0 ) THEN
	652
[4648]	653	CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter, ierr )
	654	CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter, ierr )
	655	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, ierr )
	656	CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter, status, ierr )
	657	CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter, status, ierr )
	658	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, comm_inter, status, ierr )
[667]	659	ENDIF
	660
[709]	661	CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
[4648]	662	CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr )
[709]	663	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
[4648]	664
[667]	665	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	666
	667	nx_o = nx
[4648]	668	ny_o = ny
[667]	669
[4648]	670	IF ( myid == 0 ) THEN
[709]	671
[4648]	672	CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, ierr )
	673	CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, ierr )
	674	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, status, ierr )
[709]	675	CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
	676	CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
	677	CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
[667]	678	ENDIF
	679
	680	CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
[4648]	681	CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	682	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
[667]	683
	684	ENDIF
[4648]	685
[709]	686	ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
	687	ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
[667]	688
	689	!
[4648]	690	!-- Determine if the horizontal grid and the number of PEs in ocean and atmosphere is same or not.
	691	IF ( nx_o == nx_a .AND. ny_o == ny_a .AND. &
	692	pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) THEN
[667]	693	coupling_topology = 0
	694	ELSE
	695	coupling_topology = 1
[4648]	696	ENDIF
[667]	697
	698	!
[4648]	699	!-- Determine the target PEs for the exchange between ocean and atmosphere (comm2d)
[709]	700	IF ( coupling_topology == 0 ) THEN
	701	!
[4648]	702	!-- In case of identical topologies, every atmosphere PE has exactly one ocean PE counterpart
	703	!-- and vice versa
	704	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	705	target_id = myid + numprocs
	706	ELSE
[4648]	707	target_id = myid
[667]	708	ENDIF
	709
	710	ELSE
	711	!
[4648]	712	!-- In case of nonequivalent topology in ocean and atmosphere only for PE0 in ocean and PE0 in
	713	!-- atmosphere a target_id is needed, since data echxchange between ocean and atmosphere will
	714	!-- be done only between these PEs.
[709]	715	IF ( myid == 0 ) THEN
	716
	717	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[4648]	718	target_id = numprocs
[667]	719	ELSE
	720	target_id = 0
	721	ENDIF
[709]	722
[667]	723	ENDIF
[709]	724
[667]	725	ENDIF
	726
	727	ENDIF
	728
[1]	729	#else
	730
	731	!
[4648]	732	!-- Array bounds when running on a single PE (respectively a non-parallel machine)
[1003]	733	nxl = 0
	734	nxr = nx
	735	nnx = nxr - nxl + 1
	736	nys = 0
	737	nyn = ny
	738	nny = nyn - nys + 1
	739	nzb = 0
	740	nzt = nz
	741	nnz = nz
[1]	742
[145]	743	ALLOCATE( hor_index_bounds(4,0:0) )
	744	hor_index_bounds(1,0) = nxl
	745	hor_index_bounds(2,0) = nxr
	746	hor_index_bounds(3,0) = nys
	747	hor_index_bounds(4,0) = nyn
	748
[1]	749	!
[4648]	750	!-- Array bounds for the pressure solver (in the parallel code, these bounds are the ones for the
	751	!-- transposed arrays)
[1003]	752	nys_x = nys
	753	nyn_x = nyn
	754	nzb_x = nzb + 1
	755	nzt_x = nzt
[1]	756
[1003]	757	nxl_y = nxl
	758	nxr_y = nxr
	759	nzb_y = nzb + 1
	760	nzt_y = nzt
[1]	761
[1003]	762	nxl_z = nxl
	763	nxr_z = nxr
	764	nys_z = nys
	765	nyn_z = nyn
[1]	766
	767	#endif
	768
	769	!
[4648]	770	!-- Calculate number of grid levels necessary for the multigrid poisson solver as well as the
	771	!-- gridpoint indices on each level
[1575]	772	IF ( psolver(1:9) == 'multigrid' ) THEN
[1]	773
	774	!
	775	!-- First calculate number of possible grid levels for the subdomains
	776	mg_levels_x = 1
	777	mg_levels_y = 1
	778	mg_levels_z = 1
	779
	780	i = nnx
	781	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	782	i = i / 2
	783	mg_levels_x = mg_levels_x + 1
	784	ENDDO
	785
	786	j = nny
	787	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	788	j = j / 2
	789	mg_levels_y = mg_levels_y + 1
	790	ENDDO
	791
[181]	792	k = nz ! do not use nnz because it might be > nz due to transposition
	793	! requirements
[1]	794	DO WHILE ( MOD( k, 2 ) == 0 .AND. k /= 2 )
	795	k = k / 2
	796	mg_levels_z = mg_levels_z + 1
	797	ENDDO
[2197]	798	!
[4648]	799	!-- The optimized MG-solver does not allow odd values for nz at the coarsest grid level
[2197]	800	IF ( TRIM( psolver ) /= 'multigrid_noopt' ) THEN
	801	IF ( MOD( k, 2 ) /= 0 ) mg_levels_z = mg_levels_z - 1
[3057]	802	!
[4648]	803	!-- An odd value of nz does not work. The finest level must have an even value.
[3057]	804	IF ( mg_levels_z == 0 ) THEN
	805	message_string = 'optimized multigrid method requires nz to be even'
[3058]	806	CALL message( 'init_pegrid', 'PA0495', 1, 2, 0, 6, 0 )
[3057]	807	ENDIF
[2197]	808	ENDIF
[1]	809
	810	maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
[4648]	811	!
	812	!-- Check if subdomain sizes prevents any coarsening.
	813	!-- This case, the maximum number of grid levels is 1, i.e. effectively a Gauss-Seidel scheme is
	814	!-- applied rather than a multigrid approach.
[3897]	815	!-- Give a warning in this case.
	816	IF ( maximum_grid_level == 1 .AND. mg_switch_to_pe0_level == -1 ) THEN
[4648]	817	message_string = 'No grid coarsening possible, multigrid ' // &
	818	'approach effectively reduces to a Gauss-Seidel ' // &
[3897]	819	'scheme.'
[4648]	820
[3897]	821	CALL message( 'poismg', 'PA0648', 0, 1, 0, 6, 0 )
[3890]	822	ENDIF
[1]	823
	824	!
[4648]	825	!-- Find out, if the total domain allows more levels. These additional levels are identically
	826	!-- processed on all PEs.
[197]	827	IF ( numprocs > 1 .AND. mg_switch_to_pe0_level /= -1 ) THEN
[709]	828
[1]	829	IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) ) THEN
[709]	830
[1]	831	mg_switch_to_pe0_level_l = maximum_grid_level
	832
	833	mg_levels_x = 1
	834	mg_levels_y = 1
	835
	836	i = nx+1
	837	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	838	i = i / 2
	839	mg_levels_x = mg_levels_x + 1
	840	ENDDO
	841
	842	j = ny+1
	843	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	844	j = j / 2
	845	mg_levels_y = mg_levels_y + 1
	846	ENDDO
	847
	848	maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
	849
	850	IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l ) THEN
[4648]	851	mg_switch_to_pe0_level_l = maximum_grid_level_l - mg_switch_to_pe0_level_l + 1
[1]	852	ELSE
	853	mg_switch_to_pe0_level_l = 0
	854	ENDIF
[709]	855
[1]	856	ELSE
[3057]	857
[1]	858	mg_switch_to_pe0_level_l = 0
	859	maximum_grid_level_l = maximum_grid_level
[709]	860
[1]	861	ENDIF
	862
	863	!
[4648]	864	!-- Use switch level calculated above only if it is not pre-defined by user
[1]	865	IF ( mg_switch_to_pe0_level == 0 ) THEN
	866	IF ( mg_switch_to_pe0_level_l /= 0 ) THEN
	867	mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
	868	maximum_grid_level = maximum_grid_level_l
	869	ENDIF
	870
	871	ELSE
	872	!
	873	!-- Check pre-defined value and reset to default, if neccessary
[4648]	874	IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l .OR. &
[1]	875	mg_switch_to_pe0_level >= maximum_grid_level_l ) THEN
[4648]	876	message_string = 'mg_switch_to_pe0_level ' // &
[2271]	877	'out of range and reset to 0'
[254]	878	CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
[1]	879	mg_switch_to_pe0_level = 0
	880	ELSE
	881	!
[4648]	882	!-- Use the largest number of possible levels anyway and recalculate the switch level to
	883	!-- this largest number of possible values
[1]	884	maximum_grid_level = maximum_grid_level_l
	885
	886	ENDIF
[709]	887
[1]	888	ENDIF
	889
	890	ENDIF
	891
[4648]	892	ALLOCATE( grid_level_count(maximum_grid_level), &
	893	nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
	894	nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
[1056]	895	nzt_mg(0:maximum_grid_level) )
[1]	896
	897	grid_level_count = 0
[1056]	898	!
[4648]	899	!-- Index zero required as dummy due to definition of arrays f2 and p2 in recursive subroutine
	900	!-- next_mg_level
[1056]	901	nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
[778]	902
[1]	903	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
	904
	905	DO i = maximum_grid_level, 1 , -1
	906
	907	IF ( i == mg_switch_to_pe0_level ) THEN
[1804]	908	#if defined( __parallel )
[1]	909	!
[4648]	910	!-- Save the grid size of the subdomain at the switch level, because it is needed in poismg.
[1]	911	ind(1) = nxl_l; ind(2) = nxr_l
	912	ind(3) = nys_l; ind(4) = nyn_l
	913	ind(5) = nzt_l
	914	ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
	915	CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
	916	MPI_INTEGER, comm2d, ierr )
	917	DO j = 0, numprocs-1
	918	DO k = 1, 5
	919	mg_loc_ind(k,j) = ind_all(k+j*5)
	920	ENDDO
	921	ENDDO
	922	DEALLOCATE( ind_all )
	923	!
[709]	924	!-- Calculate the grid size of the total domain
[1]	925	nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
	926	nxl_l = 0
	927	nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
	928	nys_l = 0
	929	!
[4648]	930	!-- The size of this gathered array must not be larger than the array tend, which is used
	931	!-- in the multigrid scheme as a temporary array. Therefore the subdomain size of an PE is
	932	!-- calculated and the size of the gathered grid. These values are used in routines pres
	933	!-- and poismg.
	934	subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
[778]	935	( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
[4648]	936	gathered_size = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * ( nzt_l - nzb + 2 )
[1]	937
[1804]	938	#else
[4648]	939	message_string = 'multigrid gather/scatter impossible ' // &
[1]	940	'in non parallel mode'
[254]	941	CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
[1]	942	#endif
	943	ENDIF
	944
	945	nxl_mg(i) = nxl_l
	946	nxr_mg(i) = nxr_l
	947	nys_mg(i) = nys_l
	948	nyn_mg(i) = nyn_l
	949	nzt_mg(i) = nzt_l
	950
[4648]	951	nxl_l = nxl_l / 2
[1]	952	nxr_l = nxr_l / 2
[4648]	953	nys_l = nys_l / 2
	954	nyn_l = nyn_l / 2
	955	nzt_l = nzt_l / 2
	956
[1]	957	ENDDO
	958
[780]	959	!
[4648]	960	!-- Temporary problem: Currently calculation of maxerror in routine poismg crashes if grid data
	961	!-- are collected on PE0 already on the finest grid level.
[780]	962	!-- To be solved later.
	963	IF ( maximum_grid_level == mg_switch_to_pe0_level ) THEN
	964	message_string = 'grid coarsening on subdomain level cannot be performed'
	965	CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
	966	ENDIF
	967
[1]	968	ELSE
	969
[667]	970	maximum_grid_level = 0
[1]	971
	972	ENDIF
	973
[722]	974	!
[4648]	975	!-- Default level 0 tells exchange_horiz that all ghost planes have to be exchanged. grid_level is
	976	!-- adjusted in poismg, where only one ghost plane is required.
[722]	977	grid_level = 0
[1]	978
[1804]	979	#if defined( __parallel )
[1]	980	!
	981	!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
[667]	982	ngp_y = nyn - nys + 1 + 2 * nbgp
[1]	983
	984	!
[4648]	985	!-- Define new MPI derived datatypes for the exchange of ghost points in x- and y-direction for
	986	!-- 2D-arrays (line)
	987	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, ierr )
[1]	988	CALL MPI_TYPE_COMMIT( type_x, ierr )
	989
[667]	990	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
	991	CALL MPI_TYPE_COMMIT( type_y, ierr )
[1968]	992	!
[4648]	993	!-- Define new MPI derived datatypes for the exchange of ghost points in x- and y-direction for
	994	!-- 2D-INTEGER arrays (line) - on normal grid.
	995	!-- Define types for 32-bit and 8-bit Integer. The 8-bit Integer are only required on normal grid,
	996	!-- while 32-bit Integer may be also required on coarser grid level in case of multigrid solver.
[3542]	997	!
	998	!-- 8-bit Integer
[4648]	999	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_BYTE, type_x_byte, ierr )
[3542]	1000	CALL MPI_TYPE_COMMIT( type_x_byte, ierr )
	1001
[4648]	1002	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_BYTE, type_y_byte, ierr )
[3542]	1003	CALL MPI_TYPE_COMMIT( type_y_byte, ierr )
	1004	!
	1005	!-- 32-bit Integer
[4648]	1006	ALLOCATE( type_x_int(0:maximum_grid_level), type_y_int(0:maximum_grid_level) )
	1007
	1008	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, type_x_int(0), ierr )
[1968]	1009	CALL MPI_TYPE_COMMIT( type_x_int(0), ierr )
[667]	1010
[1968]	1011	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int(0), ierr )
	1012	CALL MPI_TYPE_COMMIT( type_y_int(0), ierr )
[1]	1013	!
[4648]	1014	!-- Calculate gridpoint numbers for the exchange of ghost points along x (yz-plane for 3D-arrays)
	1015	!-- and define MPI derived data type(s) for the exchange of ghost points in y-direction (xz-plane).
	1016	!-- Do these calculations for the model grid and (if necessary) also for the coarser grid levels
	1017	!-- used in the multigrid method
	1018	ALLOCATE ( ngp_xz(0:maximum_grid_level), &
	1019	ngp_xz_int(0:maximum_grid_level), &
	1020	ngp_yz(0:maximum_grid_level), &
	1021	ngp_yz_int(0:maximum_grid_level), &
	1022	type_xz(0:maximum_grid_level), &
	1023	type_xz_int(0:maximum_grid_level), &
	1024	type_yz(0:maximum_grid_level), &
[2696]	1025	type_yz_int(0:maximum_grid_level) )
[1]	1026
	1027	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
[709]	1028
[667]	1029	!
[4648]	1030	!-- Discern between the model grid, which needs nbgp ghost points and grid levels for the multigrid
	1031	!-- scheme. In the latter case only one ghost point is necessary.
	1032	!-- First definition of MPI-datatypes for exchange of ghost layers on normal grid. The following
	1033	!-- loop is needed for data exchange in poismg.f90.
[667]	1034	!
	1035	!-- Determine number of grid points of yz-layer for exchange
	1036	ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
[709]	1037
[667]	1038	!
[709]	1039	!-- Define an MPI-datatype for the exchange of left/right boundaries.
[4648]	1040	!-- Although data are contiguous in physical memory (which does not necessarily require an
	1041	!-- MPI-derived datatype), the data exchange between left and right PE's using the MPI-derived type
	1042	!-- is 10% faster than without.
	1043	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz(0), MPI_REAL, type_xz(0), &
	1044	ierr )
[667]	1045	CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
[1]	1046
[4648]	1047	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), ierr )
[667]	1048	CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
[709]	1049
[667]	1050	!
[2696]	1051	!-- Define data types for exchange of 3D Integer arrays.
	1052	ngp_yz_int(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
	1053
[4648]	1054	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz_int(0), MPI_INTEGER, &
	1055	type_xz_int(0), ierr )
[2696]	1056	CALL MPI_TYPE_COMMIT( type_xz_int(0), ierr )
	1057
[4648]	1058	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz_int(0), ngp_yz_int(0), MPI_INTEGER, type_yz_int(0), ierr )
[2696]	1059	CALL MPI_TYPE_COMMIT( type_yz_int(0), ierr )
	1060
	1061	!
[709]	1062	!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
[1575]	1063	IF ( psolver(1:9) == 'multigrid' ) THEN
[4648]	1064	!
[709]	1065	!-- Definition of MPI-datatyoe as above, but only 1 ghost level is used
	1066	DO i = maximum_grid_level, 1 , -1
[1968]	1067	!
[4648]	1068	!-- For 3D-exchange on different multigrid level, one ghost point for REAL arrays, two ghost
	1069	!-- points for INTEGER arrays
[1575]	1070	ngp_xz(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
[667]	1071	ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	1072
[2696]	1073	ngp_xz_int(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
	1074	ngp_yz_int(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	1075	!
	1076	!-- MPI data type for REAL arrays, for xz-layers
[4648]	1077	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), MPI_REAL, type_xz(i), &
	1078	ierr )
[667]	1079	CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
[1]	1080
[2696]	1081	!
	1082	!-- MPI data type for INTEGER arrays, for xz-layers
[4648]	1083	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz_int(i), MPI_INTEGER, &
	1084	type_xz_int(i), ierr )
[2696]	1085	CALL MPI_TYPE_COMMIT( type_xz_int(i), ierr )
	1086
	1087	!
	1088	!-- MPI data type for REAL arrays, for yz-layers
[4648]	1089	CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), ierr )
[667]	1090	CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
[2696]	1091	!
	1092	!-- MPI data type for INTEGER arrays, for yz-layers
[4648]	1093	CALL MPI_TYPE_VECTOR( 1, ngp_yz_int(i), ngp_yz_int(i), MPI_INTEGER, type_yz_int(i), ierr )
[2696]	1094	CALL MPI_TYPE_COMMIT( type_yz_int(i), ierr )
[667]	1095
[1968]	1096
[4648]	1097	!-- For 2D-exchange of INTEGER arrays on coarser grid level, where 2 ghost points need to be
	1098	!-- exchanged. Only required for 32-bit Integer arrays.
	1099	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+5, 2, nyn_l-nys_l+5, MPI_INTEGER, type_x_int(i), ierr )
[1968]	1100	CALL MPI_TYPE_COMMIT( type_x_int(i), ierr )
	1101
	1102
[4648]	1103	CALL MPI_TYPE_VECTOR( 2, nyn_l-nys_l+5, nyn_l-nys_l+5, MPI_INTEGER, type_y_int(i), ierr )
[1968]	1104	CALL MPI_TYPE_COMMIT( type_y_int(i), ierr )
	1105
[667]	1106	nxl_l = nxl_l / 2
	1107	nxr_l = nxr_l / 2
	1108	nys_l = nys_l / 2
	1109	nyn_l = nyn_l / 2
	1110	nzt_l = nzt_l / 2
[709]	1111
[667]	1112	ENDDO
[709]	1113
	1114	ENDIF
[1677]	1115
[1]	1116	#endif
	1117
[1804]	1118	#if defined( __parallel )
[1]	1119	!
[1933]	1120	!-- Setting of flags for inflow/outflow/nesting conditions.
[1964]	1121	IF ( pleft == MPI_PROC_NULL ) THEN
[4648]	1122	IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'nested' .OR. &
[3182]	1123	bc_lr == 'nesting_offline' ) THEN
	1124	bc_dirichlet_l = .TRUE.
[1964]	1125	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[3182]	1126	bc_radiation_l = .TRUE.
[1]	1127	ENDIF
	1128	ENDIF
[4648]	1129
[1964]	1130	IF ( pright == MPI_PROC_NULL ) THEN
	1131	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[3182]	1132	bc_radiation_r = .TRUE.
[4648]	1133	ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'nested' .OR. &
[3182]	1134	bc_lr == 'nesting_offline' ) THEN
	1135	bc_dirichlet_r = .TRUE.
[1]	1136	ENDIF
	1137	ENDIF
	1138
[1964]	1139	IF ( psouth == MPI_PROC_NULL ) THEN
	1140	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[3182]	1141	bc_radiation_s = .TRUE.
[4648]	1142	ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'nested' .OR. &
[3182]	1143	bc_ns == 'nesting_offline' ) THEN
	1144	bc_dirichlet_s = .TRUE.
[1]	1145	ENDIF
	1146	ENDIF
	1147
[1964]	1148	IF ( pnorth == MPI_PROC_NULL ) THEN
[4648]	1149	IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'nested' .OR. &
[3182]	1150	bc_ns == 'nesting_offline' ) THEN
	1151	bc_dirichlet_n = .TRUE.
[1964]	1152	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[3182]	1153	bc_radiation_n = .TRUE.
[1]	1154	ENDIF
	1155	ENDIF
[2938]	1156	!
[4648]	1157	!-- In case of synthetic turbulence geneartor determine ids.
	1158	!-- Please note, if no forcing or nesting is applied, the generator is applied only at the left
	1159	!-- lateral boundary.
[2938]	1160	IF ( use_syn_turb_gen ) THEN
[3182]	1161	IF ( bc_dirichlet_l ) THEN
[2938]	1162	id_stg_left_l = myidx
	1163	ELSE
	1164	id_stg_left_l = 0
	1165	ENDIF
[3182]	1166	IF ( bc_dirichlet_r ) THEN
[2938]	1167	id_stg_right_l = myidx
	1168	ELSE
	1169	id_stg_right_l = 0
	1170	ENDIF
[3182]	1171	IF ( bc_dirichlet_s ) THEN
[2938]	1172	id_stg_south_l = myidy
	1173	ELSE
	1174	id_stg_south_l = 0
	1175	ENDIF
[3182]	1176	IF ( bc_dirichlet_n ) THEN
[2938]	1177	id_stg_north_l = myidy
	1178	ELSE
	1179	id_stg_north_l = 0
	1180	ENDIF
[1968]	1181
[2938]	1182	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1183	CALL MPI_ALLREDUCE( id_stg_left_l, id_stg_left, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
[2938]	1184
	1185	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1186	CALL MPI_ALLREDUCE( id_stg_right_l, id_stg_right, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
[2938]	1187
	1188	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1189	CALL MPI_ALLREDUCE( id_stg_south_l, id_stg_south, 1, MPI_INTEGER, MPI_SUM, comm1dy, ierr )
[2938]	1190
	1191	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1192	CALL MPI_ALLREDUCE( id_stg_north_l, id_stg_north, 1, MPI_INTEGER, MPI_SUM, comm1dy, ierr )
[2938]	1193
[4648]	1194	ENDIF
	1195
[151]	1196	!
	1197	!-- Broadcast the id of the inflow PE
[3182]	1198	IF ( bc_dirichlet_l ) THEN
[163]	1199	id_inflow_l = myidx
[151]	1200	ELSE
	1201	id_inflow_l = 0
	1202	ENDIF
[622]	1203	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1204	CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
[151]	1205
[163]	1206	!
	1207	!-- Broadcast the id of the recycling plane
	1208	!-- WARNING: needs to be adjusted in case of inflows other than from left side!
[3884]	1209	IF ( turbulent_inflow ) THEN
[4648]	1210
[4045]	1211	IF ( NINT( recycling_width / dx, KIND=idp ) >= nxl .AND. &
	1212	NINT( recycling_width / dx, KIND=idp ) <= nxr ) THEN
[3884]	1213	id_recycling_l = myidx
	1214	ELSE
	1215	id_recycling_l = 0
	1216	ENDIF
	1217	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1218	CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
	1219
[163]	1220	ENDIF
	1221
[2050]	1222	!
	1223	!-- Broadcast the id of the outflow PE and outflow-source plane
	1224	IF ( turbulent_outflow ) THEN
	1225
[3182]	1226	IF ( bc_radiation_r ) THEN
[2050]	1227	id_outflow_l = myidx
	1228	ELSE
	1229	id_outflow_l = 0
	1230	ENDIF
	1231	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1232	CALL MPI_ALLREDUCE( id_outflow_l, id_outflow, 1, MPI_INTEGER, MPI_SUM, &
	1233	comm1dx, ierr )
	1234
[4648]	1235	IF ( NINT( outflow_source_plane / dx ) >= nxl .AND. &
[2050]	1236	NINT( outflow_source_plane / dx ) <= nxr ) THEN
	1237	id_outflow_source_l = myidx
	1238	ELSE
	1239	id_outflow_source_l = 0
	1240	ENDIF
	1241	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1242	CALL MPI_ALLREDUCE( id_outflow_source_l, id_outflow_source, 1, MPI_INTEGER, MPI_SUM, &
	1243	comm1dx, ierr )
[2050]	1244
	1245	ENDIF
	1246
[3885]	1247	CALL location_message( 'creating virtual PE grids + MPI derived data types', 'finished' )
[1384]	1248
[1804]	1249	#else
[1159]	1250	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[3182]	1251	bc_dirichlet_l = .TRUE.
	1252	bc_radiation_r = .TRUE.
[1159]	1253	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[3182]	1254	bc_radiation_l = .TRUE.
	1255	bc_dirichlet_r = .TRUE.
[1]	1256	ENDIF
	1257
[1159]	1258	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[3182]	1259	bc_dirichlet_n = .TRUE.
	1260	bc_radiation_s = .TRUE.
[1159]	1261	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[3182]	1262	bc_radiation_n = .TRUE.
	1263	bc_dirichlet_s = .TRUE.
[1]	1264	ENDIF
	1265	#endif
[807]	1266
[106]	1267	!
[4648]	1268	!-- At the inflow or outflow, u or v, respectively, have to be calculated for one more grid point.
[3182]	1269	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[106]	1270	nxlu = nxl + 1
	1271	ELSE
	1272	nxlu = nxl
	1273	ENDIF
[3182]	1274	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[106]	1275	nysv = nys + 1
	1276	ELSE
	1277	nysv = nys
	1278	ENDIF
[1]	1279
[114]	1280	!
	1281	!-- Allocate wall flag arrays used in the multigrid solver
[1575]	1282	IF ( psolver(1:9) == 'multigrid' ) THEN
[114]	1283
	1284	DO i = maximum_grid_level, 1, -1
	1285
	1286	SELECT CASE ( i )
	1287
	1288	CASE ( 1 )
[4648]	1289	ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1, &
	1290	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1291	nxl_mg(i)-1:nxr_mg(i)+1) )
	1292
	1293	CASE ( 2 )
[4648]	1294	ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1, &
	1295	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1296	nxl_mg(i)-1:nxr_mg(i)+1) )
	1297
	1298	CASE ( 3 )
[4648]	1299	ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1, &
	1300	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1301	nxl_mg(i)-1:nxr_mg(i)+1) )
	1302
	1303	CASE ( 4 )
[4648]	1304	ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1, &
	1305	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1306	nxl_mg(i)-1:nxr_mg(i)+1) )
	1307
	1308	CASE ( 5 )
[4648]	1309	ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1, &
	1310	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1311	nxl_mg(i)-1:nxr_mg(i)+1) )
	1312
	1313	CASE ( 6 )
[4648]	1314	ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1, &
	1315	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1316	nxl_mg(i)-1:nxr_mg(i)+1) )
	1317
	1318	CASE ( 7 )
[4648]	1319	ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1, &
	1320	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1321	nxl_mg(i)-1:nxr_mg(i)+1) )
	1322
	1323	CASE ( 8 )
[4648]	1324	ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1, &
	1325	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1326	nxl_mg(i)-1:nxr_mg(i)+1) )
	1327
	1328	CASE ( 9 )
[4648]	1329	ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1, &
	1330	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1331	nxl_mg(i)-1:nxr_mg(i)+1) )
	1332
	1333	CASE ( 10 )
[4648]	1334	ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1, &
	1335	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1336	nxl_mg(i)-1:nxr_mg(i)+1) )
	1337
	1338	CASE DEFAULT
[254]	1339	message_string = 'more than 10 multigrid levels'
	1340	CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
[114]	1341
	1342	END SELECT
	1343
	1344	ENDDO
	1345
	1346	ENDIF
	1347
[1]	1348	END SUBROUTINE init_pegrid

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |