Home

Context Navigation

init_pegrid.f90 @ 4894

Last change on this file since 4894 was 4893, checked in by raasch, 4 years ago
revised output of surface data via MPI-IO for better performance
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 4893 2021-03-02 16:39:14Z raasch $
[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	!
[1003]	361	!-- Calculate array bounds along x-direction for every PE.
[4648]	362	ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), nysf(0:pdims(2)-1) )
[1]	363
[1003]	364	IF ( MOD( nx+1 , pdims(1) ) /= 0 ) THEN
[4648]	365	WRITE( message_string, * ) 'x-direction: gridpoint number (' ,nx+1, ') ', &
	366	'is not an& integral multiple of the number', &
	367	' of processors (', pdims(1), ')'
[254]	368	CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
[1]	369	ELSE
[1003]	370	nnx = ( nx + 1 ) / pdims(1)
[4264]	371	ENDIF
[1]	372
	373	!
	374	!-- Left and right array bounds, number of gridpoints
	375	DO i = 0, pdims(1)-1
	376	nxlf(i) = i * nnx
	377	nxrf(i) = ( i + 1 ) * nnx - 1
	378	ENDDO
	379
	380	!
	381	!-- Calculate array bounds in y-direction for every PE.
[1003]	382	IF ( MOD( ny+1 , pdims(2) ) /= 0 ) THEN
[4648]	383	WRITE( message_string, * ) 'y-direction: gridpoint number (', ny+1, ') ', &
	384	'is not an& integral multiple of the number', &
	385	' of processors (', pdims(2), ')'
[254]	386	CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
[1]	387	ELSE
[1003]	388	nny = ( ny + 1 ) / pdims(2)
[4264]	389	ENDIF
[1]	390
	391	!
	392	!-- South and north array bounds
	393	DO j = 0, pdims(2)-1
	394	nysf(j) = j * nny
	395	nynf(j) = ( j + 1 ) * nny - 1
	396	ENDDO
	397
	398	!
	399	!-- Local array bounds of the respective PEs
[1003]	400	nxl = nxlf(pcoord(1))
	401	nxr = nxrf(pcoord(1))
	402	nys = nysf(pcoord(2))
	403	nyn = nynf(pcoord(2))
	404	nzb = 0
	405	nzt = nz
	406	nnz = nz
[1]	407
	408	!
[4648]	409	!-- Set switches to define if the PE is situated at the border of the virtual processor grid
[707]	410	IF ( nxl == 0 ) left_border_pe = .TRUE.
	411	IF ( nxr == nx ) right_border_pe = .TRUE.
	412	IF ( nys == 0 ) south_border_pe = .TRUE.
	413	IF ( nyn == ny ) north_border_pe = .TRUE.
	414
	415	!
[4648]	416	!-- Calculate array bounds and gridpoint numbers for the transposed arrays (needed in the pressure
	417	!-- solver)
	418	!-- For the transposed arrays, cyclic boundaries as well as top and bottom boundaries are omitted,
	419	!-- because they are obstructive to the transposition
[1]	420
	421	!
	422	!-- 1. transposition z --> x
[1001]	423	!-- This transposition is not neccessary in case of a 1d-decomposition along x
[2938]	424	IF ( psolver == 'poisfft' .OR. calculate_spectra ) THEN
[1304]	425
[1922]	426	IF ( pdims(2) /= 1 ) THEN
	427	IF ( MOD( nz , pdims(1) ) /= 0 ) THEN
[4648]	428	WRITE( message_string, * ) 'transposition z --> x:& ', &
	429	'nz=', nz, ' is not an integral multiple ', &
	430	'of pdims(1)=', pdims(1)
[1922]	431	CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
	432	ENDIF
[1]	433	ENDIF
[1922]	434
	435	nys_x = nys
	436	nyn_x = nyn
	437	nny_x = nny
	438	nnz_x = nz / pdims(1)
	439	nzb_x = 1 + myidx * nnz_x
	440	nzt_x = ( myidx + 1 ) * nnz_x
	441	sendrecvcount_zx = nnx * nny * nnz_x
	442
[1]	443	ENDIF
	444
[1922]	445
[4264]	446	IF ( psolver == 'poisfft' ) THEN
[1]	447	!
[1922]	448	!-- 2. transposition x --> y
	449	IF ( MOD( nx+1 , pdims(2) ) /= 0 ) THEN
[4648]	450	WRITE( message_string, * ) 'transposition x --> y:& ', &
	451	'nx+1=', nx+1, ' is not an integral ', &
	452	'multiple of pdims(2)=', pdims(2)
[1922]	453	CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
	454	ENDIF
[1]	455
[1922]	456	nnz_y = nnz_x
	457	nzb_y = nzb_x
	458	nzt_y = nzt_x
	459	nnx_y = (nx+1) / pdims(2)
	460	nxl_y = myidy * nnx_y
	461	nxr_y = ( myidy + 1 ) * nnx_y - 1
	462	sendrecvcount_xy = nnx_y * nny_x * nnz_y
[1]	463	!
[4264]	464	!-- 3. transposition y --> z
[1922]	465	!-- (ELSE: x --> y in case of 1D-decomposition along x)
	466	nxl_z = nxl_y
	467	nxr_z = nxr_y
	468	nny_z = (ny+1) / pdims(1)
	469	nys_z = myidx * nny_z
	470	nyn_z = ( myidx + 1 ) * nny_z - 1
	471	sendrecvcount_yz = nnx_y * nny_z * nnz_y
[1304]	472
[1922]	473	IF ( pdims(2) /= 1 ) THEN
[1]	474	!
[1922]	475	!-- y --> z
	476	!-- This transposition is not neccessary in case of a 1d-decomposition
	477	!-- along x, except that the uptream-spline method is switched on
	478	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[4648]	479	WRITE( message_string, * ) 'transposition y --> z:& ', &
	480	'ny+1=', ny+1, ' is not an integral ', &
	481	'multiple of pdims(1)=', pdims(1)
[1922]	482	CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
	483	ENDIF
[1]	484
[1922]	485	ELSE
[1]	486	!
[1922]	487	!-- x --> y
	488	!-- This condition must be fulfilled for a 1D-decomposition along x
	489	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[4648]	490	WRITE( message_string, * ) 'transposition x --> y:& ', &
	491	'ny+1=', ny+1, ' is not an integral ', &
	492	'multiple of pdims(1)=', pdims(1)
[1922]	493	CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
	494	ENDIF
	495
[1]	496	ENDIF
	497
	498	ENDIF
	499
	500	!
	501	!-- Indices for direct transpositions z --> y (used for calculating spectra)
[1922]	502	IF ( calculate_spectra ) THEN
[1003]	503	IF ( MOD( nz, pdims(2) ) /= 0 ) THEN
[4648]	504	WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
	505	'for spectra):& nz=', nz, ' is not an ', &
	506	'integral multiple of pdims(2)=', pdims(2)
[254]	507	CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
[1]	508	ELSE
[1003]	509	nxl_yd = nxl
	510	nxr_yd = nxr
	511	nzb_yd = 1 + myidy * ( nz / pdims(2) )
	512	nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
	513	sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
[1]	514	ENDIF
	515	ENDIF
	516
[1922]	517	IF ( psolver == 'poisfft' .OR. calculate_spectra ) THEN
[1]	518	!
[4648]	519	!-- Indices for direct transpositions y --> x
[1922]	520	!-- (they are only possible in case of a 1d-decomposition along x)
	521	IF ( pdims(2) == 1 ) THEN
	522	nny_x = nny / pdims(1)
	523	nys_x = myid * nny_x
	524	nyn_x = ( myid + 1 ) * nny_x - 1
	525	nzb_x = 1
	526	nzt_x = nz
	527	sendrecvcount_xy = nnx * nny_x * nz
	528	ENDIF
	529
[1]	530	ENDIF
	531
[1922]	532	IF ( psolver == 'poisfft' ) THEN
[1]	533	!
[4648]	534	!-- Indices for direct transpositions x --> y
[1922]	535	!-- (they are only possible in case of a 1d-decomposition along y)
	536	IF ( pdims(1) == 1 ) THEN
	537	nnx_y = nnx / pdims(2)
	538	nxl_y = myid * nnx_y
	539	nxr_y = ( myid + 1 ) * nnx_y - 1
	540	nzb_y = 1
	541	nzt_y = nz
	542	sendrecvcount_xy = nnx_y * nny * nz
	543	ENDIF
	544
[1]	545	ENDIF
	546
	547	!
	548	!-- Arrays for storing the array bounds are needed any more
	549	DEALLOCATE( nxlf , nxrf , nynf , nysf )
	550
[807]	551
[145]	552	!
	553	!-- Collect index bounds from other PEs (to be written to restart file later)
	554	ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
	555
	556	IF ( myid == 0 ) THEN
	557
	558	hor_index_bounds(1,0) = nxl
	559	hor_index_bounds(2,0) = nxr
	560	hor_index_bounds(3,0) = nys
	561	hor_index_bounds(4,0) = nyn
	562
	563	!
	564	!-- Receive data from all other PEs
	565	DO i = 1, numprocs-1
[4648]	566	CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, ierr )
[145]	567	hor_index_bounds(:,i) = ibuf(1:4)
	568	ENDDO
	569
	570	ELSE
	571	!
	572	!-- Send index bounds to PE0
	573	ibuf(1) = nxl
	574	ibuf(2) = nxr
	575	ibuf(3) = nys
	576	ibuf(4) = nyn
	577	CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
	578
	579	ENDIF
	580
[807]	581
[1]	582	#if defined( __print )
	583	!
	584	!-- Control output
	585	IF ( myid == 0 ) THEN
	586	PRINT, ' processor topology *'
	587	PRINT*, ' '
[4648]	588	PRINT*, 'myid pcoord left right south north idx idy nxl: nxr',' nys: nyn'
	589	PRINT*, '------------------------------------------------------------','-----------'
	590	WRITE (*,1000) 0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, myidx, myidy, nxl, &
	591	nxr, nys, nyn
	592	1000 FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3,2(2X,I4,':',I4))
[1]	593
	594	!
[108]	595	!-- Receive data from the other PEs
[1]	596	DO i = 1,numprocs-1
[4648]	597	CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, ierr )
[1]	598	WRITE (*,1000) i, ( ibuf(j) , j = 1,12 )
	599	ENDDO
	600	ELSE
	601
	602	!
	603	!-- Send data to PE0
	604	ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
	605	ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
	606	ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
	607	ibuf(12) = nyn
[4648]	608	CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )
[1]	609	ENDIF
	610	#endif
	611
[4648]	612	!
[709]	613	!-- Determine the number of ghost point layers
[4648]	614	IF ( scalar_advec == 'ws-scheme' .OR. &
[3999]	615	momentum_advec == 'ws-scheme' .OR. nested_run ) THEN
[667]	616	nbgp = 3
	617	ELSE
	618	nbgp = 1
[4648]	619	ENDIF
[667]	620
[102]	621	!
[4648]	622	!-- Check that the number of computational grid points is not smaller than the number of ghost
	623	!-- points.
[4241]	624	IF ( nnx < nbgp ) THEN
	625	WRITE( message_string, * ) 'number of subdomain grid points along x (', nnx, ') is smaller',&
	626	'than the number of ghost points (', nbgp, ')'
	627	CALL message( 'init_pegrid', 'PA0682', 1, 2, 0, 6, 0 )
	628	ENDIF
	629	IF ( nny < nbgp ) THEN
	630	WRITE( message_string, * ) 'number of subdomain grid points along y (', nny, ') is smaller',&
	631	'than the number of ghost points (', nbgp, ')'
	632	CALL message( 'init_pegrid', 'PA0683', 1, 2, 0, 6, 0 )
	633	ENDIF
	634
	635	!
[4648]	636	!-- Create a new MPI derived datatype for the exchange of surface (xy) data, which is needed for
	637	!-- coupled atmosphere-ocean runs.
[709]	638	!-- First, calculate number of grid points of an xy-plane.
[667]	639	ngp_xy = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
[102]	640	CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
	641	CALL MPI_TYPE_COMMIT( type_xy, ierr )
[667]	642
[4564]	643	IF ( TRIM( coupling_mode ) /= 'uncoupled' ) THEN
[4648]	644
[667]	645	!
	646	!-- Pass the number of grid points of the atmosphere model to
	647	!-- the ocean model and vice versa
	648	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	649
	650	nx_a = nx
	651	ny_a = ny
	652
[709]	653	IF ( myid == 0 ) THEN
	654
[4648]	655	CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter, ierr )
	656	CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter, ierr )
	657	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, ierr )
	658	CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter, status, ierr )
	659	CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter, status, ierr )
	660	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, comm_inter, status, ierr )
[667]	661	ENDIF
	662
[709]	663	CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
[4648]	664	CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr )
[709]	665	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
[4648]	666
[667]	667	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	668
	669	nx_o = nx
[4648]	670	ny_o = ny
[667]	671
[4648]	672	IF ( myid == 0 ) THEN
[709]	673
[4648]	674	CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, ierr )
	675	CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, ierr )
	676	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, status, ierr )
[709]	677	CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
	678	CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
	679	CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
[667]	680	ENDIF
	681
	682	CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
[4648]	683	CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	684	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
[667]	685
	686	ENDIF
[4648]	687
[709]	688	ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
	689	ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
[667]	690
	691	!
[4648]	692	!-- Determine if the horizontal grid and the number of PEs in ocean and atmosphere is same or not.
	693	IF ( nx_o == nx_a .AND. ny_o == ny_a .AND. &
	694	pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) THEN
[667]	695	coupling_topology = 0
	696	ELSE
	697	coupling_topology = 1
[4648]	698	ENDIF
[667]	699
	700	!
[4648]	701	!-- Determine the target PEs for the exchange between ocean and atmosphere (comm2d)
[709]	702	IF ( coupling_topology == 0 ) THEN
	703	!
[4648]	704	!-- In case of identical topologies, every atmosphere PE has exactly one ocean PE counterpart
	705	!-- and vice versa
	706	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	707	target_id = myid + numprocs
	708	ELSE
[4648]	709	target_id = myid
[667]	710	ENDIF
	711
	712	ELSE
	713	!
[4648]	714	!-- In case of nonequivalent topology in ocean and atmosphere only for PE0 in ocean and PE0 in
	715	!-- atmosphere a target_id is needed, since data echxchange between ocean and atmosphere will
	716	!-- be done only between these PEs.
[709]	717	IF ( myid == 0 ) THEN
	718
	719	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[4648]	720	target_id = numprocs
[667]	721	ELSE
	722	target_id = 0
	723	ENDIF
[709]	724
[667]	725	ENDIF
[709]	726
[667]	727	ENDIF
	728
	729	ENDIF
	730
[1]	731	#else
	732
	733	!
[4648]	734	!-- Array bounds when running on a single PE (respectively a non-parallel machine)
[1003]	735	nxl = 0
	736	nxr = nx
	737	nnx = nxr - nxl + 1
	738	nys = 0
	739	nyn = ny
	740	nny = nyn - nys + 1
	741	nzb = 0
	742	nzt = nz
	743	nnz = nz
[1]	744
[145]	745	ALLOCATE( hor_index_bounds(4,0:0) )
	746	hor_index_bounds(1,0) = nxl
	747	hor_index_bounds(2,0) = nxr
	748	hor_index_bounds(3,0) = nys
	749	hor_index_bounds(4,0) = nyn
	750
[1]	751	!
[4648]	752	!-- Array bounds for the pressure solver (in the parallel code, these bounds are the ones for the
	753	!-- transposed arrays)
[1003]	754	nys_x = nys
	755	nyn_x = nyn
	756	nzb_x = nzb + 1
	757	nzt_x = nzt
[1]	758
[1003]	759	nxl_y = nxl
	760	nxr_y = nxr
	761	nzb_y = nzb + 1
	762	nzt_y = nzt
[1]	763
[1003]	764	nxl_z = nxl
	765	nxr_z = nxr
	766	nys_z = nys
	767	nyn_z = nyn
[1]	768
	769	#endif
	770
	771	!
[4648]	772	!-- Calculate number of grid levels necessary for the multigrid poisson solver as well as the
	773	!-- gridpoint indices on each level
[1575]	774	IF ( psolver(1:9) == 'multigrid' ) THEN
[1]	775
	776	!
	777	!-- First calculate number of possible grid levels for the subdomains
	778	mg_levels_x = 1
	779	mg_levels_y = 1
	780	mg_levels_z = 1
	781
	782	i = nnx
	783	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	784	i = i / 2
	785	mg_levels_x = mg_levels_x + 1
	786	ENDDO
	787
	788	j = nny
	789	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	790	j = j / 2
	791	mg_levels_y = mg_levels_y + 1
	792	ENDDO
	793
[181]	794	k = nz ! do not use nnz because it might be > nz due to transposition
	795	! requirements
[1]	796	DO WHILE ( MOD( k, 2 ) == 0 .AND. k /= 2 )
	797	k = k / 2
	798	mg_levels_z = mg_levels_z + 1
	799	ENDDO
[2197]	800	!
[4648]	801	!-- The optimized MG-solver does not allow odd values for nz at the coarsest grid level
[2197]	802	IF ( TRIM( psolver ) /= 'multigrid_noopt' ) THEN
	803	IF ( MOD( k, 2 ) /= 0 ) mg_levels_z = mg_levels_z - 1
[3057]	804	!
[4648]	805	!-- An odd value of nz does not work. The finest level must have an even value.
[3057]	806	IF ( mg_levels_z == 0 ) THEN
	807	message_string = 'optimized multigrid method requires nz to be even'
[3058]	808	CALL message( 'init_pegrid', 'PA0495', 1, 2, 0, 6, 0 )
[3057]	809	ENDIF
[2197]	810	ENDIF
[1]	811
	812	maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
[4648]	813	!
	814	!-- Check if subdomain sizes prevents any coarsening.
	815	!-- This case, the maximum number of grid levels is 1, i.e. effectively a Gauss-Seidel scheme is
	816	!-- applied rather than a multigrid approach.
[3897]	817	!-- Give a warning in this case.
	818	IF ( maximum_grid_level == 1 .AND. mg_switch_to_pe0_level == -1 ) THEN
[4648]	819	message_string = 'No grid coarsening possible, multigrid ' // &
	820	'approach effectively reduces to a Gauss-Seidel ' // &
[3897]	821	'scheme.'
[4648]	822
[3897]	823	CALL message( 'poismg', 'PA0648', 0, 1, 0, 6, 0 )
[3890]	824	ENDIF
[1]	825
	826	!
[4648]	827	!-- Find out, if the total domain allows more levels. These additional levels are identically
	828	!-- processed on all PEs.
[197]	829	IF ( numprocs > 1 .AND. mg_switch_to_pe0_level /= -1 ) THEN
[709]	830
[1]	831	IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) ) THEN
[709]	832
[1]	833	mg_switch_to_pe0_level_l = maximum_grid_level
	834
	835	mg_levels_x = 1
	836	mg_levels_y = 1
	837
	838	i = nx+1
	839	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	840	i = i / 2
	841	mg_levels_x = mg_levels_x + 1
	842	ENDDO
	843
	844	j = ny+1
	845	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	846	j = j / 2
	847	mg_levels_y = mg_levels_y + 1
	848	ENDDO
	849
	850	maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
	851
	852	IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l ) THEN
[4648]	853	mg_switch_to_pe0_level_l = maximum_grid_level_l - mg_switch_to_pe0_level_l + 1
[1]	854	ELSE
	855	mg_switch_to_pe0_level_l = 0
	856	ENDIF
[709]	857
[1]	858	ELSE
[3057]	859
[1]	860	mg_switch_to_pe0_level_l = 0
	861	maximum_grid_level_l = maximum_grid_level
[709]	862
[1]	863	ENDIF
	864
	865	!
[4648]	866	!-- Use switch level calculated above only if it is not pre-defined by user
[1]	867	IF ( mg_switch_to_pe0_level == 0 ) THEN
	868	IF ( mg_switch_to_pe0_level_l /= 0 ) THEN
	869	mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
	870	maximum_grid_level = maximum_grid_level_l
	871	ENDIF
	872
	873	ELSE
	874	!
	875	!-- Check pre-defined value and reset to default, if neccessary
[4648]	876	IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l .OR. &
[1]	877	mg_switch_to_pe0_level >= maximum_grid_level_l ) THEN
[4648]	878	message_string = 'mg_switch_to_pe0_level ' // &
[2271]	879	'out of range and reset to 0'
[254]	880	CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
[1]	881	mg_switch_to_pe0_level = 0
	882	ELSE
	883	!
[4648]	884	!-- Use the largest number of possible levels anyway and recalculate the switch level to
	885	!-- this largest number of possible values
[1]	886	maximum_grid_level = maximum_grid_level_l
	887
	888	ENDIF
[709]	889
[1]	890	ENDIF
	891
	892	ENDIF
	893
[4648]	894	ALLOCATE( grid_level_count(maximum_grid_level), &
	895	nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
	896	nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
[1056]	897	nzt_mg(0:maximum_grid_level) )
[1]	898
	899	grid_level_count = 0
[1056]	900	!
[4648]	901	!-- Index zero required as dummy due to definition of arrays f2 and p2 in recursive subroutine
	902	!-- next_mg_level
[1056]	903	nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
[778]	904
[1]	905	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
	906
	907	DO i = maximum_grid_level, 1 , -1
	908
	909	IF ( i == mg_switch_to_pe0_level ) THEN
[1804]	910	#if defined( __parallel )
[1]	911	!
[4648]	912	!-- Save the grid size of the subdomain at the switch level, because it is needed in poismg.
[1]	913	ind(1) = nxl_l; ind(2) = nxr_l
	914	ind(3) = nys_l; ind(4) = nyn_l
	915	ind(5) = nzt_l
	916	ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
	917	CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
	918	MPI_INTEGER, comm2d, ierr )
	919	DO j = 0, numprocs-1
	920	DO k = 1, 5
	921	mg_loc_ind(k,j) = ind_all(k+j*5)
	922	ENDDO
	923	ENDDO
	924	DEALLOCATE( ind_all )
	925	!
[709]	926	!-- Calculate the grid size of the total domain
[1]	927	nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
	928	nxl_l = 0
	929	nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
	930	nys_l = 0
	931	!
[4648]	932	!-- The size of this gathered array must not be larger than the array tend, which is used
	933	!-- in the multigrid scheme as a temporary array. Therefore the subdomain size of an PE is
	934	!-- calculated and the size of the gathered grid. These values are used in routines pres
	935	!-- and poismg.
	936	subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
[778]	937	( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
[4648]	938	gathered_size = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * ( nzt_l - nzb + 2 )
[1]	939
[1804]	940	#else
[4648]	941	message_string = 'multigrid gather/scatter impossible ' // &
[1]	942	'in non parallel mode'
[254]	943	CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
[1]	944	#endif
	945	ENDIF
	946
	947	nxl_mg(i) = nxl_l
	948	nxr_mg(i) = nxr_l
	949	nys_mg(i) = nys_l
	950	nyn_mg(i) = nyn_l
	951	nzt_mg(i) = nzt_l
	952
[4648]	953	nxl_l = nxl_l / 2
[1]	954	nxr_l = nxr_l / 2
[4648]	955	nys_l = nys_l / 2
	956	nyn_l = nyn_l / 2
	957	nzt_l = nzt_l / 2
	958
[1]	959	ENDDO
	960
[780]	961	!
[4648]	962	!-- Temporary problem: Currently calculation of maxerror in routine poismg crashes if grid data
	963	!-- are collected on PE0 already on the finest grid level.
[780]	964	!-- To be solved later.
	965	IF ( maximum_grid_level == mg_switch_to_pe0_level ) THEN
	966	message_string = 'grid coarsening on subdomain level cannot be performed'
	967	CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
	968	ENDIF
	969
[1]	970	ELSE
	971
[667]	972	maximum_grid_level = 0
[1]	973
	974	ENDIF
	975
[722]	976	!
[4648]	977	!-- Default level 0 tells exchange_horiz that all ghost planes have to be exchanged. grid_level is
	978	!-- adjusted in poismg, where only one ghost plane is required.
[722]	979	grid_level = 0
[1]	980
[1804]	981	#if defined( __parallel )
[1]	982	!
	983	!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
[667]	984	ngp_y = nyn - nys + 1 + 2 * nbgp
[1]	985
	986	!
[4648]	987	!-- Define new MPI derived datatypes for the exchange of ghost points in x- and y-direction for
	988	!-- 2D-arrays (line)
	989	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, ierr )
[1]	990	CALL MPI_TYPE_COMMIT( type_x, ierr )
	991
[667]	992	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
	993	CALL MPI_TYPE_COMMIT( type_y, ierr )
[1968]	994	!
[4648]	995	!-- Define new MPI derived datatypes for the exchange of ghost points in x- and y-direction for
	996	!-- 2D-INTEGER arrays (line) - on normal grid.
	997	!-- Define types for 32-bit and 8-bit Integer. The 8-bit Integer are only required on normal grid,
	998	!-- while 32-bit Integer may be also required on coarser grid level in case of multigrid solver.
[3542]	999	!
	1000	!-- 8-bit Integer
[4648]	1001	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_BYTE, type_x_byte, ierr )
[3542]	1002	CALL MPI_TYPE_COMMIT( type_x_byte, ierr )
	1003
[4648]	1004	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_BYTE, type_y_byte, ierr )
[3542]	1005	CALL MPI_TYPE_COMMIT( type_y_byte, ierr )
	1006	!
	1007	!-- 32-bit Integer
[4648]	1008	ALLOCATE( type_x_int(0:maximum_grid_level), type_y_int(0:maximum_grid_level) )
	1009
	1010	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, type_x_int(0), ierr )
[1968]	1011	CALL MPI_TYPE_COMMIT( type_x_int(0), ierr )
[667]	1012
[1968]	1013	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int(0), ierr )
	1014	CALL MPI_TYPE_COMMIT( type_y_int(0), ierr )
[1]	1015	!
[4648]	1016	!-- Calculate gridpoint numbers for the exchange of ghost points along x (yz-plane for 3D-arrays)
	1017	!-- and define MPI derived data type(s) for the exchange of ghost points in y-direction (xz-plane).
	1018	!-- Do these calculations for the model grid and (if necessary) also for the coarser grid levels
	1019	!-- used in the multigrid method
	1020	ALLOCATE ( ngp_xz(0:maximum_grid_level), &
	1021	ngp_xz_int(0:maximum_grid_level), &
	1022	ngp_yz(0:maximum_grid_level), &
	1023	ngp_yz_int(0:maximum_grid_level), &
	1024	type_xz(0:maximum_grid_level), &
	1025	type_xz_int(0:maximum_grid_level), &
	1026	type_yz(0:maximum_grid_level), &
[2696]	1027	type_yz_int(0:maximum_grid_level) )
[1]	1028
	1029	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
[709]	1030
[667]	1031	!
[4648]	1032	!-- Discern between the model grid, which needs nbgp ghost points and grid levels for the multigrid
	1033	!-- scheme. In the latter case only one ghost point is necessary.
	1034	!-- First definition of MPI-datatypes for exchange of ghost layers on normal grid. The following
	1035	!-- loop is needed for data exchange in poismg.f90.
[667]	1036	!
	1037	!-- Determine number of grid points of yz-layer for exchange
	1038	ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
[709]	1039
[667]	1040	!
[709]	1041	!-- Define an MPI-datatype for the exchange of left/right boundaries.
[4648]	1042	!-- Although data are contiguous in physical memory (which does not necessarily require an
	1043	!-- MPI-derived datatype), the data exchange between left and right PE's using the MPI-derived type
	1044	!-- is 10% faster than without.
	1045	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz(0), MPI_REAL, type_xz(0), &
	1046	ierr )
[667]	1047	CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
[1]	1048
[4648]	1049	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), ierr )
[667]	1050	CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
[709]	1051
[667]	1052	!
[2696]	1053	!-- Define data types for exchange of 3D Integer arrays.
	1054	ngp_yz_int(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
	1055
[4648]	1056	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz_int(0), MPI_INTEGER, &
	1057	type_xz_int(0), ierr )
[2696]	1058	CALL MPI_TYPE_COMMIT( type_xz_int(0), ierr )
	1059
[4648]	1060	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz_int(0), ngp_yz_int(0), MPI_INTEGER, type_yz_int(0), ierr )
[2696]	1061	CALL MPI_TYPE_COMMIT( type_yz_int(0), ierr )
	1062
	1063	!
[709]	1064	!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
[1575]	1065	IF ( psolver(1:9) == 'multigrid' ) THEN
[4648]	1066	!
[709]	1067	!-- Definition of MPI-datatyoe as above, but only 1 ghost level is used
	1068	DO i = maximum_grid_level, 1 , -1
[1968]	1069	!
[4648]	1070	!-- For 3D-exchange on different multigrid level, one ghost point for REAL arrays, two ghost
	1071	!-- points for INTEGER arrays
[1575]	1072	ngp_xz(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
[667]	1073	ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	1074
[2696]	1075	ngp_xz_int(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
	1076	ngp_yz_int(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	1077	!
	1078	!-- MPI data type for REAL arrays, for xz-layers
[4648]	1079	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), MPI_REAL, type_xz(i), &
	1080	ierr )
[667]	1081	CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
[1]	1082
[2696]	1083	!
	1084	!-- MPI data type for INTEGER arrays, for xz-layers
[4648]	1085	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz_int(i), MPI_INTEGER, &
	1086	type_xz_int(i), ierr )
[2696]	1087	CALL MPI_TYPE_COMMIT( type_xz_int(i), ierr )
	1088
	1089	!
	1090	!-- MPI data type for REAL arrays, for yz-layers
[4648]	1091	CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), ierr )
[667]	1092	CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
[2696]	1093	!
	1094	!-- MPI data type for INTEGER arrays, for yz-layers
[4648]	1095	CALL MPI_TYPE_VECTOR( 1, ngp_yz_int(i), ngp_yz_int(i), MPI_INTEGER, type_yz_int(i), ierr )
[2696]	1096	CALL MPI_TYPE_COMMIT( type_yz_int(i), ierr )
[667]	1097
[1968]	1098
[4648]	1099	!-- For 2D-exchange of INTEGER arrays on coarser grid level, where 2 ghost points need to be
	1100	!-- exchanged. Only required for 32-bit Integer arrays.
	1101	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+5, 2, nyn_l-nys_l+5, MPI_INTEGER, type_x_int(i), ierr )
[1968]	1102	CALL MPI_TYPE_COMMIT( type_x_int(i), ierr )
	1103
	1104
[4648]	1105	CALL MPI_TYPE_VECTOR( 2, nyn_l-nys_l+5, nyn_l-nys_l+5, MPI_INTEGER, type_y_int(i), ierr )
[1968]	1106	CALL MPI_TYPE_COMMIT( type_y_int(i), ierr )
	1107
[667]	1108	nxl_l = nxl_l / 2
	1109	nxr_l = nxr_l / 2
	1110	nys_l = nys_l / 2
	1111	nyn_l = nyn_l / 2
	1112	nzt_l = nzt_l / 2
[709]	1113
[667]	1114	ENDDO
[709]	1115
	1116	ENDIF
[1677]	1117
[1]	1118	#endif
	1119
[1804]	1120	#if defined( __parallel )
[1]	1121	!
[1933]	1122	!-- Setting of flags for inflow/outflow/nesting conditions.
[1964]	1123	IF ( pleft == MPI_PROC_NULL ) THEN
[4648]	1124	IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'nested' .OR. &
[3182]	1125	bc_lr == 'nesting_offline' ) THEN
	1126	bc_dirichlet_l = .TRUE.
[1964]	1127	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[3182]	1128	bc_radiation_l = .TRUE.
[1]	1129	ENDIF
	1130	ENDIF
[4648]	1131
[1964]	1132	IF ( pright == MPI_PROC_NULL ) THEN
	1133	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[3182]	1134	bc_radiation_r = .TRUE.
[4648]	1135	ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'nested' .OR. &
[3182]	1136	bc_lr == 'nesting_offline' ) THEN
	1137	bc_dirichlet_r = .TRUE.
[1]	1138	ENDIF
	1139	ENDIF
	1140
[1964]	1141	IF ( psouth == MPI_PROC_NULL ) THEN
	1142	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[3182]	1143	bc_radiation_s = .TRUE.
[4648]	1144	ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'nested' .OR. &
[3182]	1145	bc_ns == 'nesting_offline' ) THEN
	1146	bc_dirichlet_s = .TRUE.
[1]	1147	ENDIF
	1148	ENDIF
	1149
[1964]	1150	IF ( pnorth == MPI_PROC_NULL ) THEN
[4648]	1151	IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'nested' .OR. &
[3182]	1152	bc_ns == 'nesting_offline' ) THEN
	1153	bc_dirichlet_n = .TRUE.
[1964]	1154	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[3182]	1155	bc_radiation_n = .TRUE.
[1]	1156	ENDIF
	1157	ENDIF
[2938]	1158	!
[4648]	1159	!-- In case of synthetic turbulence geneartor determine ids.
	1160	!-- Please note, if no forcing or nesting is applied, the generator is applied only at the left
	1161	!-- lateral boundary.
[4848]	1162	IF ( syn_turb_gen ) THEN
[3182]	1163	IF ( bc_dirichlet_l ) THEN
[2938]	1164	id_stg_left_l = myidx
	1165	ELSE
	1166	id_stg_left_l = 0
	1167	ENDIF
[3182]	1168	IF ( bc_dirichlet_r ) THEN
[2938]	1169	id_stg_right_l = myidx
	1170	ELSE
	1171	id_stg_right_l = 0
	1172	ENDIF
[3182]	1173	IF ( bc_dirichlet_s ) THEN
[2938]	1174	id_stg_south_l = myidy
	1175	ELSE
	1176	id_stg_south_l = 0
	1177	ENDIF
[3182]	1178	IF ( bc_dirichlet_n ) THEN
[2938]	1179	id_stg_north_l = myidy
	1180	ELSE
	1181	id_stg_north_l = 0
	1182	ENDIF
[1968]	1183
[2938]	1184	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1185	CALL MPI_ALLREDUCE( id_stg_left_l, id_stg_left, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
[2938]	1186
	1187	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1188	CALL MPI_ALLREDUCE( id_stg_right_l, id_stg_right, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
[2938]	1189
	1190	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1191	CALL MPI_ALLREDUCE( id_stg_south_l, id_stg_south, 1, MPI_INTEGER, MPI_SUM, comm1dy, ierr )
[2938]	1192
	1193	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1194	CALL MPI_ALLREDUCE( id_stg_north_l, id_stg_north, 1, MPI_INTEGER, MPI_SUM, comm1dy, ierr )
[2938]	1195
[4648]	1196	ENDIF
	1197
[151]	1198	!
	1199	!-- Broadcast the id of the inflow PE
[3182]	1200	IF ( bc_dirichlet_l ) THEN
[163]	1201	id_inflow_l = myidx
[151]	1202	ELSE
	1203	id_inflow_l = 0
	1204	ENDIF
[622]	1205	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1206	CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
[151]	1207
[163]	1208	!
	1209	!-- Broadcast the id of the recycling plane
	1210	!-- WARNING: needs to be adjusted in case of inflows other than from left side!
[3884]	1211	IF ( turbulent_inflow ) THEN
[4648]	1212
[4045]	1213	IF ( NINT( recycling_width / dx, KIND=idp ) >= nxl .AND. &
	1214	NINT( recycling_width / dx, KIND=idp ) <= nxr ) THEN
[3884]	1215	id_recycling_l = myidx
	1216	ELSE
	1217	id_recycling_l = 0
	1218	ENDIF
	1219	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1220	CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
	1221
[163]	1222	ENDIF
	1223
[2050]	1224	!
	1225	!-- Broadcast the id of the outflow PE and outflow-source plane
	1226	IF ( turbulent_outflow ) THEN
	1227
[3182]	1228	IF ( bc_radiation_r ) THEN
[2050]	1229	id_outflow_l = myidx
	1230	ELSE
	1231	id_outflow_l = 0
	1232	ENDIF
	1233	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1234	CALL MPI_ALLREDUCE( id_outflow_l, id_outflow, 1, MPI_INTEGER, MPI_SUM, &
	1235	comm1dx, ierr )
	1236
[4648]	1237	IF ( NINT( outflow_source_plane / dx ) >= nxl .AND. &
[2050]	1238	NINT( outflow_source_plane / dx ) <= nxr ) THEN
	1239	id_outflow_source_l = myidx
	1240	ELSE
	1241	id_outflow_source_l = 0
	1242	ENDIF
	1243	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4648]	1244	CALL MPI_ALLREDUCE( id_outflow_source_l, id_outflow_source, 1, MPI_INTEGER, MPI_SUM, &
	1245	comm1dx, ierr )
[2050]	1246
	1247	ENDIF
	1248
[3885]	1249	CALL location_message( 'creating virtual PE grids + MPI derived data types', 'finished' )
[1384]	1250
[1804]	1251	#else
[1159]	1252	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[3182]	1253	bc_dirichlet_l = .TRUE.
	1254	bc_radiation_r = .TRUE.
[1159]	1255	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[3182]	1256	bc_radiation_l = .TRUE.
	1257	bc_dirichlet_r = .TRUE.
[1]	1258	ENDIF
	1259
[1159]	1260	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[3182]	1261	bc_dirichlet_n = .TRUE.
	1262	bc_radiation_s = .TRUE.
[1159]	1263	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[3182]	1264	bc_radiation_n = .TRUE.
	1265	bc_dirichlet_s = .TRUE.
[1]	1266	ENDIF
	1267	#endif
[807]	1268
[106]	1269	!
[4648]	1270	!-- At the inflow or outflow, u or v, respectively, have to be calculated for one more grid point.
[3182]	1271	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[106]	1272	nxlu = nxl + 1
	1273	ELSE
	1274	nxlu = nxl
	1275	ENDIF
[3182]	1276	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[106]	1277	nysv = nys + 1
	1278	ELSE
	1279	nysv = nys
	1280	ENDIF
[1]	1281
[114]	1282	!
	1283	!-- Allocate wall flag arrays used in the multigrid solver
[1575]	1284	IF ( psolver(1:9) == 'multigrid' ) THEN
[114]	1285
	1286	DO i = maximum_grid_level, 1, -1
	1287
	1288	SELECT CASE ( i )
	1289
	1290	CASE ( 1 )
[4648]	1291	ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1, &
	1292	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1293	nxl_mg(i)-1:nxr_mg(i)+1) )
	1294
	1295	CASE ( 2 )
[4648]	1296	ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1, &
	1297	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1298	nxl_mg(i)-1:nxr_mg(i)+1) )
	1299
	1300	CASE ( 3 )
[4648]	1301	ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1, &
	1302	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1303	nxl_mg(i)-1:nxr_mg(i)+1) )
	1304
	1305	CASE ( 4 )
[4648]	1306	ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1, &
	1307	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1308	nxl_mg(i)-1:nxr_mg(i)+1) )
	1309
	1310	CASE ( 5 )
[4648]	1311	ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1, &
	1312	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1313	nxl_mg(i)-1:nxr_mg(i)+1) )
	1314
	1315	CASE ( 6 )
[4648]	1316	ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1, &
	1317	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1318	nxl_mg(i)-1:nxr_mg(i)+1) )
	1319
	1320	CASE ( 7 )
[4648]	1321	ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1, &
	1322	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1323	nxl_mg(i)-1:nxr_mg(i)+1) )
	1324
	1325	CASE ( 8 )
[4648]	1326	ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1, &
	1327	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1328	nxl_mg(i)-1:nxr_mg(i)+1) )
	1329
	1330	CASE ( 9 )
[4648]	1331	ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1, &
	1332	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1333	nxl_mg(i)-1:nxr_mg(i)+1) )
	1334
	1335	CASE ( 10 )
[4648]	1336	ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1, &
	1337	nys_mg(i)-1:nyn_mg(i)+1, &
[114]	1338	nxl_mg(i)-1:nxr_mg(i)+1) )
	1339
	1340	CASE DEFAULT
[254]	1341	message_string = 'more than 10 multigrid levels'
	1342	CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
[114]	1343
	1344	END SELECT
	1345
	1346	ENDDO
	1347
	1348	ENDIF
	1349
[1]	1350	END SUBROUTINE init_pegrid

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

Context Navigation

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

Download in other formats: