Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |