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

Last change on this file since 77 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

  • Property svn:keywords set to Id
File size: 24.4 KB
Line 
1 SUBROUTINE init_pegrid
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: init_pegrid.f90 77 2007-03-29 04:26:56Z raasch $
11!
12! 75 2007-03-22 09:54:05Z raasch
13! uxrp, vynp eliminated,
14! dirichlet/neumann changed to dirichlet/radiation, etc.,
15! poisfft_init is only called if fft-solver is switched on
16!
17! RCS Log replace by Id keyword, revision history cleaned up
18!
19! Revision 1.28  2006/04/26 13:23:32  raasch
20! lcmuk does not understand the !$ comment so a cpp-directive is required
21!
22! Revision 1.1  1997/07/24 11:15:09  raasch
23! Initial revision
24!
25!
26! Description:
27! ------------
28! Determination of the virtual processor topology (if not prescribed by the
29! user)and computation of the grid point number and array bounds of the local
30! domains.
31!------------------------------------------------------------------------------!
32
33    USE control_parameters
34    USE fft_xy
35    USE indices
36    USE pegrid
37    USE poisfft_mod
38    USE poisfft_hybrid_mod
39    USE statistics
40    USE transpose_indices
41
42
43    IMPLICIT NONE
44
45    INTEGER ::  gathered_size, i, ind(5), j, k, maximum_grid_level_l,     &
46                mg_switch_to_pe0_level_l, mg_levels_x, mg_levels_y,       &
47                mg_levels_z, nnx_y, nnx_z, nny_x, nny_z, nnz_x, nnz_y,    &
48                numproc_sqr, nx_total, nxl_l, nxr_l, nyn_l, nys_l, nzb_l, &
49                nzt_l, omp_get_num_threads, subdomain_size
50
51    INTEGER, DIMENSION(:), ALLOCATABLE ::  ind_all, nxlf, nxrf, nynf, nysf
52
53    LOGICAL ::  found
54
55!
56!-- Get the number of OpenMP threads
57    !$OMP PARALLEL
58#if defined( __lcmuk )
59    threads_per_task = omp_get_num_threads()
60#else
61!$  threads_per_task = omp_get_num_threads()
62#endif
63    !$OMP END PARALLEL
64
65
66#if defined( __parallel )
67!
68!-- Determine the processor topology or check it, if prescribed by the user
69    IF ( npex == -1  .AND.  npey == -1 )  THEN
70
71!
72!--    Automatic determination of the topology
73!--    The default on SMP- and cluster-hosts is a 1d-decomposition along x
74#if defined( __lcmuk )
75       host = 'lcmuk'
76#endif
77       IF ( host(1:3) == 'ibm'  .OR.  host(1:3) == 'nec'  .OR. &
78            host(1:2) == 'lc'   .OR.  host(1:3) == 'dec' )  THEN
79
80          pdims(1) = numprocs
81          pdims(2) = 1
82
83       ELSE
84
85          numproc_sqr = SQRT( REAL( numprocs ) )
86          pdims(1)    = MAX( numproc_sqr , 1 )
87          DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
88             pdims(1) = pdims(1) - 1
89          ENDDO
90          pdims(2) = numprocs / pdims(1)
91
92       ENDIF
93
94    ELSEIF ( npex /= -1  .AND.  npey /= -1 )  THEN
95
96!
97!--    Prescribed by user. Number of processors on the prescribed topology
98!--    must be equal to the number of PEs available to the job
99       IF ( ( npex * npey ) /= numprocs )  THEN
100          PRINT*, '+++ init_pegrid:'
101          PRINT*, '    number of PEs of the prescribed topology (', npex*npey, &
102                      ') does not match the number of PEs available to the ',  &
103                      'job (', numprocs, ')'
104          CALL local_stop
105       ENDIF
106       pdims(1) = npex
107       pdims(2) = npey
108
109    ELSE
110!
111!--    If the processor topology is prescribed by the user, the number of
112!--    PEs must be given in both directions
113       PRINT*, '+++ init_pegrid:'
114       PRINT*, '    if the processor topology is prescribed by the user, ',   &
115                    'both values of "npex" and "npey" must be given in the ', &
116                    'NAMELIST-parameter file'
117       CALL local_stop
118
119    ENDIF
120
121!
122!-- The hybrid solver can only be used in case of a 1d-decomposition along x
123    IF ( pdims(2) /= 1  .AND.  psolver == 'poisfft_hybrid' )  THEN
124       IF ( myid == 0 )  THEN
125          PRINT*, '*** init_pegrid: psolver = "poisfft_hybrid" can only be'
126          PRINT*, '                 used in case of a 1d-decomposition along x'
127       ENDIF
128    ENDIF
129
130!
131!-- If necessary, set horizontal boundary conditions to non-cyclic
132    IF ( bc_lr /= 'cyclic' )  cyclic(1) = .FALSE.
133    IF ( bc_ns /= 'cyclic' )  cyclic(2) = .FALSE.
134
135!
136!-- Create the virtual processor grid
137    CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
138                          comm2d, ierr )
139    CALL MPI_COMM_RANK( comm2d, myid, ierr )
140    WRITE (myid_char,'(''_'',I4.4)')  myid
141
142    CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
143    CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
144    CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
145
146!
147!-- Determine sub-topologies for transpositions
148!-- Transposition from z to x:
149    remain_dims(1) = .TRUE.
150    remain_dims(2) = .FALSE.
151    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
152    CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
153!
154!-- Transposition from x to y
155    remain_dims(1) = .FALSE.
156    remain_dims(2) = .TRUE.
157    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
158    CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
159
160
161!
162!-- Find a grid (used for array d) which will match the transposition demands
163    IF ( grid_matching == 'strict' )  THEN
164
165       nxa = nx;  nya = ny;  nza = nz
166
167    ELSE
168
169       found = .FALSE.
170   xn: DO  nxa = nx, 2*nx
171!
172!--       Meet conditions for nx
173          IF ( MOD( nxa+1, pdims(1) ) /= 0 .OR. &
174               MOD( nxa+1, pdims(2) ) /= 0 )  CYCLE xn
175
176      yn: DO  nya = ny, 2*ny
177!
178!--          Meet conditions for ny
179             IF ( MOD( nya+1, pdims(2) ) /= 0 .OR. &
180                  MOD( nya+1, pdims(1) ) /= 0 )  CYCLE yn
181
182
183         zn: DO  nza = nz, 2*nz
184!
185!--             Meet conditions for nz
186                IF ( ( MOD( nza, pdims(1) ) /= 0  .AND.  pdims(1) /= 1  .AND. &
187                       pdims(2) /= 1 )  .OR.                                  &
188                     ( MOD( nza, pdims(2) ) /= 0  .AND.  dt_dosp /= 9999999.9 &
189                     ) )  THEN
190                   CYCLE zn
191                ELSE
192                   found = .TRUE.
193                   EXIT xn
194                ENDIF
195
196             ENDDO zn
197
198          ENDDO yn
199
200       ENDDO xn
201
202       IF ( .NOT. found )  THEN
203          IF ( myid == 0 )  THEN
204             PRINT*,'+++ init_pegrid: no matching grid for transpositions found'
205          ENDIF
206          CALL local_stop
207       ENDIF
208
209    ENDIF
210
211!
212!-- Calculate array bounds in x-direction for every PE.
213!-- The last PE along x may get less grid points than the others
214    ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
215              nysf(0:pdims(2)-1), nnx_pe(0:pdims(1)-1), nny_pe(0:pdims(2)-1) )
216
217    IF ( MOD( nxa+1 , pdims(1) ) /= 0 )  THEN
218       IF ( myid == 0 )  THEN
219          PRINT*,'+++ x-direction:  gridpoint number (',nx+1,') is not an'
220          PRINT*,'                  integral divisor of the number of proces', &
221                                   &'sors (', pdims(1),')'
222       ENDIF
223       CALL local_stop
224    ELSE
225       nnx  = ( nxa + 1 ) / pdims(1)
226       IF ( nnx*pdims(1) - ( nx + 1) > nnx )  THEN
227          IF ( myid == 0 )  THEN
228             PRINT*,'+++ x-direction: nx does not match the requirements ', &
229                         'given by the number of PEs'
230             PRINT*,'                 used'
231             PRINT*,'    please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
232                         - ( nx + 1 ) ) ), ' instead of nx =', nx
233          ENDIF
234          CALL local_stop
235       ENDIF
236    ENDIF   
237
238!
239!-- Left and right array bounds, number of gridpoints
240    DO  i = 0, pdims(1)-1
241       nxlf(i)   = i * nnx
242       nxrf(i)   = ( i + 1 ) * nnx - 1
243       nnx_pe(i) = MIN( nx, nxrf(i) ) - nxlf(i) + 1
244    ENDDO
245
246!
247!-- Calculate array bounds in y-direction for every PE.
248    IF ( MOD( nya+1 , pdims(2) ) /= 0 )  THEN
249       IF ( myid == 0 )  THEN
250          PRINT*,'+++ y-direction:  gridpoint number (',ny+1,') is not an'
251          PRINT*,'                  integral divisor of the number of proces', &
252                                   &'sors (', pdims(2),')'
253       ENDIF
254       CALL local_stop
255    ELSE
256       nny  = ( nya + 1 ) / pdims(2)
257       IF ( nny*pdims(2) - ( ny + 1) > nny )  THEN
258          IF ( myid == 0 )  THEN
259             PRINT*,'+++ x-direction: nx does not match the requirements ', &
260                         'given by the number of PEs'
261             PRINT*,'                 used'
262             PRINT*,'    please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
263                         - ( nx + 1 ) ) ), ' instead of nx =', nx
264          ENDIF
265          CALL local_stop
266       ENDIF
267    ENDIF   
268
269!
270!-- South and north array bounds
271    DO  j = 0, pdims(2)-1
272       nysf(j)   = j * nny
273       nynf(j)   = ( j + 1 ) * nny - 1
274       nny_pe(j) = MIN( ny, nynf(j) ) - nysf(j) + 1
275    ENDDO
276
277!
278!-- Local array bounds of the respective PEs
279    nxl  = nxlf(pcoord(1))
280    nxra = nxrf(pcoord(1))
281    nxr  = MIN( nx, nxra )
282    nys  = nysf(pcoord(2))
283    nyna = nynf(pcoord(2))
284    nyn  = MIN( ny, nyna )
285    nzb  = 0
286    nzta = nza
287    nzt  = MIN( nz, nzta )
288    nnz  = nza
289
290!
291!-- Calculate array bounds and gridpoint numbers for the transposed arrays
292!-- (needed in the pressure solver)
293!-- For the transposed arrays, cyclic boundaries as well as top and bottom
294!-- boundaries are omitted, because they are obstructive to the transposition
295
296!
297!-- 1. transposition  z --> x
298!-- This transposition is not neccessary in case of a 1d-decomposition along x,
299!-- except that the uptream-spline method is switched on
300    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
301         scalar_advec == 'ups-scheme' )  THEN
302
303       IF ( pdims(2) == 1  .AND. ( momentum_advec == 'ups-scheme'  .OR. &
304            scalar_advec == 'ups-scheme' ) )  THEN
305          IF ( myid == 0 )  THEN
306             PRINT*,'+++ WARNING: init_pegrid: 1d-decomposition along x ', &
307                                &'chosen but nz restrictions may occur'
308             PRINT*,'             since ups-scheme is activated'
309          ENDIF
310       ENDIF
311       nys_x  = nys
312       nyn_xa = nyna
313       nyn_x  = nyn
314       nny_x  = nny
315       IF ( MOD( nza , pdims(1) ) /= 0 )  THEN
316          IF ( myid == 0 )  THEN
317             PRINT*,'+++ transposition z --> x:'
318             PRINT*,'    nz=',nz,' is not an integral divisior of pdims(1)=', &
319                    &pdims(1)
320          ENDIF
321          CALL local_stop
322       ENDIF
323       nnz_x  = nza / pdims(1)
324       nzb_x  = 1 + myidx * nnz_x
325       nzt_xa = ( myidx + 1 ) * nnz_x
326       nzt_x  = MIN( nzt, nzt_xa )
327
328       sendrecvcount_zx = nnx * nny * nnz_x
329
330    ENDIF
331
332!
333!-- 2. transposition  x --> y
334    nnz_y  = nnz_x
335    nzb_y  = nzb_x
336    nzt_ya = nzt_xa
337    nzt_y  = nzt_x
338    IF ( MOD( nxa+1 , pdims(2) ) /= 0 )  THEN
339       IF ( myid == 0 )  THEN
340          PRINT*,'+++ transposition x --> y:'
341          PRINT*,'    nx+1=',nx+1,' is not an integral divisor of ',&
342                 &'pdims(2)=',pdims(2)
343       ENDIF
344       CALL local_stop
345    ENDIF
346    nnx_y = (nxa+1) / pdims(2)
347    nxl_y = myidy * nnx_y
348    nxr_ya = ( myidy + 1 ) * nnx_y - 1
349    nxr_y  = MIN( nx, nxr_ya )
350
351    sendrecvcount_xy = nnx_y * nny_x * nnz_y
352
353!
354!-- 3. transposition  y --> z  (ELSE:  x --> y  in case of 1D-decomposition
355!-- along x)
356    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
357         scalar_advec == 'ups-scheme' )  THEN
358!
359!--    y --> z
360!--    This transposition is not neccessary in case of a 1d-decomposition
361!--    along x, except that the uptream-spline method is switched on
362       nnx_z  = nnx_y
363       nxl_z  = nxl_y
364       nxr_za = nxr_ya
365       nxr_z  = nxr_y
366       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
367          IF ( myid == 0 )  THEN
368             PRINT*,'+++ Transposition y --> z:'
369             PRINT*,'    ny+1=',ny+1,' is not an integral divisor of ',&
370                    &'pdims(1)=',pdims(1)
371          ENDIF
372          CALL local_stop
373       ENDIF
374       nny_z  = (nya+1) / pdims(1)
375       nys_z  = myidx * nny_z
376       nyn_za = ( myidx + 1 ) * nny_z - 1
377       nyn_z  = MIN( ny, nyn_za )
378
379       sendrecvcount_yz = nnx_y * nny_z * nnz_y
380
381    ELSE
382!
383!--    x --> y. This condition must be fulfilled for a 1D-decomposition along x
384       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
385          IF ( myid == 0 )  THEN
386             PRINT*,'+++ Transposition x --> y:'
387             PRINT*,'    ny+1=',ny+1,' is not an integral divisor of ',&
388                    &'pdims(1)=',pdims(1)
389          ENDIF
390          CALL local_stop
391       ENDIF
392
393    ENDIF
394
395!
396!-- Indices for direct transpositions z --> y (used for calculating spectra)
397    IF ( dt_dosp /= 9999999.9 )  THEN
398       IF ( MOD( nza, pdims(2) ) /= 0 )  THEN
399          IF ( myid == 0 )  THEN
400             PRINT*,'+++ Direct transposition z --> y (needed for spectra):'
401             PRINT*,'    nz=',nz,' is not an integral divisor of ',&
402                    &'pdims(2)=',pdims(2)
403          ENDIF
404          CALL local_stop
405       ELSE
406          nxl_yd  = nxl
407          nxr_yda = nxra
408          nxr_yd  = nxr
409          nzb_yd  = 1 + myidy * ( nza / pdims(2) )
410          nzt_yda = ( myidy + 1 ) * ( nza / pdims(2) )
411          nzt_yd  = MIN( nzt, nzt_yda )
412
413          sendrecvcount_zyd = nnx * nny * ( nza / pdims(2) )
414       ENDIF
415    ENDIF
416
417!
418!-- Indices for direct transpositions y --> x (they are only possible in case
419!-- of a 1d-decomposition along x)
420    IF ( pdims(2) == 1 )  THEN
421       nny_x  = nny / pdims(1)
422       nys_x  = myid * nny_x
423       nyn_xa = ( myid + 1 ) * nny_x - 1
424       nyn_x  = MIN( ny, nyn_xa )
425       nzb_x  = 1
426       nzt_xa = nza
427       nzt_x  = nz
428       sendrecvcount_xy = nnx * nny_x * nza
429    ENDIF
430
431!
432!-- Indices for direct transpositions x --> y (they are only possible in case
433!-- of a 1d-decomposition along y)
434    IF ( pdims(1) == 1 )  THEN
435       nnx_y  = nnx / pdims(2)
436       nxl_y  = myid * nnx_y
437       nxr_ya = ( myid + 1 ) * nnx_y - 1
438       nxr_y  = MIN( nx, nxr_ya )
439       nzb_y  = 1
440       nzt_ya = nza
441       nzt_y  = nz
442       sendrecvcount_xy = nnx_y * nny * nza
443    ENDIF
444
445!
446!-- Arrays for storing the array bounds are needed any more
447    DEALLOCATE( nxlf , nxrf , nynf , nysf )
448
449#if defined( __print )
450!
451!-- Control output
452    IF ( myid == 0 )  THEN
453       PRINT*, '*** processor topology ***'
454       PRINT*, ' '
455       PRINT*, 'myid   pcoord    left right  south north  idx idy   nxl: nxr',&
456               &'   nys: nyn'
457       PRINT*, '------------------------------------------------------------',&
458               &'-----------'
459       WRITE (*,1000)  0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
460                       myidx, myidy, nxl, nxr, nys, nyn
4611000   FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
462               2(2X,I4,':',I4))
463
464!
465!--    Recieve data from the other PEs
466       DO  i = 1,numprocs-1
467          CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
468                         ierr )
469          WRITE (*,1000)  i, ( ibuf(j) , j = 1,12 )
470       ENDDO
471    ELSE
472
473!
474!--    Send data to PE0
475       ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
476       ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
477       ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
478       ibuf(12) = nyn
479       CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
480    ENDIF
481#endif
482
483#else
484
485!
486!-- Array bounds when running on a single PE (respectively a non-parallel
487!-- machine)
488    nxl  = 0
489    nxr  = nx
490    nxra = nx
491    nnx  = nxr - nxl + 1
492    nys  = 0
493    nyn  = ny
494    nyna = ny
495    nny  = nyn - nys + 1
496    nzb  = 0
497    nzt  = nz
498    nzta = nz
499    nnz  = nz
500
501!
502!-- Array bounds for the pressure solver (in the parallel code, these bounds
503!-- are the ones for the transposed arrays)
504    nys_x  = nys
505    nyn_x  = nyn
506    nyn_xa = nyn
507    nzb_x  = nzb + 1
508    nzt_x  = nzt
509    nzt_xa = nzt
510
511    nxl_y  = nxl
512    nxr_y  = nxr
513    nxr_ya = nxr
514    nzb_y  = nzb + 1
515    nzt_y  = nzt
516    nzt_ya = nzt
517
518    nxl_z  = nxl
519    nxr_z  = nxr
520    nxr_za = nxr
521    nys_z  = nys
522    nyn_z  = nyn
523    nyn_za = nyn
524
525#endif
526
527!
528!-- Calculate number of grid levels necessary for the multigrid poisson solver
529!-- as well as the gridpoint indices on each level
530    IF ( psolver == 'multigrid' )  THEN
531
532!
533!--    First calculate number of possible grid levels for the subdomains
534       mg_levels_x = 1
535       mg_levels_y = 1
536       mg_levels_z = 1
537
538       i = nnx
539       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
540          i = i / 2
541          mg_levels_x = mg_levels_x + 1
542       ENDDO
543
544       j = nny
545       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
546          j = j / 2
547          mg_levels_y = mg_levels_y + 1
548       ENDDO
549
550       k = nnz
551       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
552          k = k / 2
553          mg_levels_z = mg_levels_z + 1
554       ENDDO
555
556       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
557
558!
559!--    Find out, if the total domain allows more levels. These additional
560!--    levels are processed on PE0 only.
561       IF ( numprocs > 1 )  THEN
562          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
563             mg_switch_to_pe0_level_l = maximum_grid_level
564
565             mg_levels_x = 1
566             mg_levels_y = 1
567
568             i = nx+1
569             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
570                i = i / 2
571                mg_levels_x = mg_levels_x + 1
572             ENDDO
573
574             j = ny+1
575             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
576                j = j / 2
577                mg_levels_y = mg_levels_y + 1
578             ENDDO
579
580             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
581
582             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
583                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
584                                           mg_switch_to_pe0_level_l + 1
585             ELSE
586                mg_switch_to_pe0_level_l = 0
587             ENDIF
588          ELSE
589             mg_switch_to_pe0_level_l = 0
590             maximum_grid_level_l = maximum_grid_level
591          ENDIF
592
593!
594!--       Use switch level calculated above only if it is not pre-defined
595!--       by user
596          IF ( mg_switch_to_pe0_level == 0 )  THEN
597
598             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
599                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
600                maximum_grid_level     = maximum_grid_level_l
601             ENDIF
602
603          ELSE
604!
605!--          Check pre-defined value and reset to default, if neccessary
606             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.  &
607                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
608                IF ( myid == 0 )  THEN
609                   PRINT*, '+++ WARNING init_pegrid: mg_switch_to_pe0_level ', &
610                               'out of range and reset to default (=0)'
611                ENDIF
612                mg_switch_to_pe0_level = 0
613             ELSE
614!
615!--             Use the largest number of possible levels anyway and recalculate
616!--             the switch level to this largest number of possible values
617                maximum_grid_level = maximum_grid_level_l
618
619             ENDIF
620          ENDIF
621
622       ENDIF
623
624       ALLOCATE( grid_level_count(maximum_grid_level),                   &
625                 nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
626                 nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
627                 nzt_mg(maximum_grid_level) )
628
629       grid_level_count = 0
630       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
631
632       DO  i = maximum_grid_level, 1 , -1
633
634          IF ( i == mg_switch_to_pe0_level )  THEN
635#if defined( __parallel )
636!
637!--          Save the grid size of the subdomain at the switch level, because
638!--          it is needed in poismg.
639!--          Array bounds of the local subdomain grids are gathered on PE0
640             ind(1) = nxl_l; ind(2) = nxr_l
641             ind(3) = nys_l; ind(4) = nyn_l
642             ind(5) = nzt_l
643             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
644             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
645                                 MPI_INTEGER, comm2d, ierr )
646             DO  j = 0, numprocs-1
647                DO  k = 1, 5
648                   mg_loc_ind(k,j) = ind_all(k+j*5)
649                ENDDO
650             ENDDO
651             DEALLOCATE( ind_all )
652!
653!--          Calculate the grid size of the total domain gathered on PE0
654             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
655             nxl_l = 0
656             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
657             nys_l = 0
658!
659!--          The size of this gathered array must not be larger than the
660!--          array tend, which is used in the multigrid scheme as a temporary
661!--          array
662             subdomain_size = ( nxr - nxl + 3 )     * ( nyn - nys + 3 )     * &
663                              ( nzt - nzb + 2 )
664             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
665                              ( nzt_l - nzb + 2 )
666
667             IF ( gathered_size > subdomain_size )  THEN
668                IF ( myid == 0 )  THEN
669                   PRINT*, '+++ init_pegrid: not enough memory for storing ', &
670                               'gathered multigrid data on PE0'
671                ENDIF
672                CALL local_stop
673             ENDIF
674#else
675             PRINT*, '+++ init_pegrid: multigrid gather/scatter impossible ', &
676                          'in non parallel mode'
677             CALL local_stop
678#endif
679          ENDIF
680
681          nxl_mg(i) = nxl_l
682          nxr_mg(i) = nxr_l
683          nys_mg(i) = nys_l
684          nyn_mg(i) = nyn_l
685          nzt_mg(i) = nzt_l
686
687          nxl_l = nxl_l / 2 
688          nxr_l = nxr_l / 2
689          nys_l = nys_l / 2 
690          nyn_l = nyn_l / 2 
691          nzt_l = nzt_l / 2 
692       ENDDO
693
694    ELSE
695
696       maximum_grid_level = 1
697
698    ENDIF
699
700    grid_level = maximum_grid_level
701
702#if defined( __parallel )
703!
704!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
705    ngp_y  = nyn - nys + 1
706
707!
708!-- Define a new MPI derived datatype for the exchange of ghost points in
709!-- y-direction for 2D-arrays (line)
710    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_REAL, type_x, ierr )
711    CALL MPI_TYPE_COMMIT( type_x, ierr )
712    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_INTEGER, type_x_int, ierr )
713    CALL MPI_TYPE_COMMIT( type_x_int, ierr )
714
715!
716!-- Calculate gridpoint numbers for the exchange of ghost points along x
717!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
718!-- exchange of ghost points in y-direction (xz-plane).
719!-- Do these calculations for the model grid and (if necessary) also
720!-- for the coarser grid levels used in the multigrid method
721    ALLOCATE ( ngp_yz(maximum_grid_level), type_xz(maximum_grid_level) )
722
723    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
724         
725    DO i = maximum_grid_level, 1 , -1
726       ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
727
728       CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
729                             MPI_REAL, type_xz(i), ierr )
730       CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
731
732       nxl_l = nxl_l / 2 
733       nxr_l = nxr_l / 2
734       nys_l = nys_l / 2 
735       nyn_l = nyn_l / 2 
736       nzt_l = nzt_l / 2 
737    ENDDO
738#endif
739
740#if defined( __parallel )
741!
742!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
743!-- horizontal boundary conditions. Set variables for extending array u (v)
744!-- by one gridpoint on the left/rightmost (northest/southest) processor
745    IF ( pleft == MPI_PROC_NULL )  THEN
746       IF ( bc_lr == 'dirichlet/radiation' )  THEN
747          inflow_l  = .TRUE.
748       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
749          outflow_l = .TRUE.
750       ENDIF
751    ENDIF
752
753    IF ( pright == MPI_PROC_NULL )  THEN
754       IF ( bc_lr == 'dirichlet/radiation' )  THEN
755          outflow_r = .TRUE.
756       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
757          inflow_r  = .TRUE.
758       ENDIF
759    ENDIF
760
761    IF ( psouth == MPI_PROC_NULL )  THEN
762       IF ( bc_ns == 'dirichlet/radiation' )  THEN
763          outflow_s = .TRUE.
764       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
765          inflow_s  = .TRUE.
766       ENDIF
767    ENDIF
768
769    IF ( pnorth == MPI_PROC_NULL )  THEN
770       IF ( bc_ns == 'dirichlet/radiation' )  THEN
771          inflow_n  = .TRUE.
772       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
773          outflow_n = .TRUE.
774       ENDIF
775    ENDIF
776
777#else
778    IF ( bc_lr == 'dirichlet/radiation' )  THEN
779       inflow_l  = .TRUE.
780       outflow_r = .TRUE.
781    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
782       outflow_l = .TRUE.
783       inflow_r  = .TRUE.
784    ENDIF
785
786    IF ( bc_ns == 'dirichlet/radiation' )  THEN
787       inflow_n  = .TRUE.
788       outflow_s = .TRUE.
789    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
790       outflow_n = .TRUE.
791       inflow_s  = .TRUE.
792    ENDIF
793#endif
794
795    IF ( psolver == 'poisfft_hybrid' )  THEN
796       CALL poisfft_hybrid_ini
797    ELSEIF ( psolver == 'poisfft' )  THEN
798       CALL poisfft_init
799    ENDIF
800
801 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.