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

Last change on this file since 667 was 667, checked in by suehring, 11 years ago

summary:


Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90


Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90


Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90


Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90


When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90


Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90


Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90


Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90


Steering parameter dissipation_control added in inipar.

Makefile


Module advec_ws added.

Modules


Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90


New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90


Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90


Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90


Output of total array size was adapted to nbgp.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90


Calls of exchange_horiz are modified.

diffusion_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90


Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90


Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90


determination of target_id's moved to init_pegrid

init_pt_anomaly.f90


Call of exchange_horiz are modified.

init_rankine.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90


Output of advection scheme.

poismg.f90


Calls of exchange_horiz are modified.

prandtl_fluxes.f90


Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90


+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90


Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90


Allocation with nbgp.

wall_fluxes.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90


Array bounds and nx, ny adapted with nbgp.

sor.f90


bug removed in declaration of ddzw(), nz replaced by nzt+1

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