source: palm/trunk/SOURCE/init_dvrp.f90 @ 590

Last change on this file since 590 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

  • Property svn:keywords set to Id
File size: 26.9 KB
RevLine 
[1]1  SUBROUTINE init_dvrp
2
3!------------------------------------------------------------------------------!
[254]4! Current revisions:
[82]5! -----------------
[392]6!
7!
8! Former revisions:
9! -----------------
10! $Id: init_dvrp.f90 392 2009-09-24 10:39:14Z heinze $
11!
12! 284 2009-04-06 06:36:10Z raasch
[264]13! Definition of a colortable to be used for particles.
[260]14! Output names are changed: surface=groundplate, buildings=topography
[254]15! Output of messages replaced by message handling routine.
[242]16! Clipping implemented.
[237]17! Polygon reduction for building and ground plate isosurface. Reduction level
18! for buildings can be chosen with parameter cluster_size.
[273]19! Steering, splitting, and rtsp routines not used on nec.
[1]20! ToDo: checking of mode_dvrp for legal values is not correct
[206]21! Implementation of a MPI-1 coupling: __mpi2 adjustments for MPI_COMM_WORLD
[226]22!
23! 210 2008-11-06 08:54:02Z raasch
24! DVRP arguments changed to single precision, mode pathlines added
25!
[198]26! 155 2008-03-28 10:56:30Z letzel
27! introduce prefix_chr to ensure unique dvrp_file path
28!
[139]29! 130 2007-11-13 14:08:40Z letzel
30! allow two instead of one digit to specify isosurface and slicer variables
31! Test output of isosurface on camera file
32!
[83]33! 82 2007-04-16 15:40:52Z raasch
34! Preprocessor strings for different linux clusters changed to "lc",
35! routine local_flush is used for buffer flushing
36!
[39]37! 17 2007-02-19 01:57:39Z raasch
[17]38! dvrp_output_local activated for all streams
39!
40! 13 2007-02-14 12:15:07Z raasch
[3]41! RCS Log replace by Id keyword, revision history cleaned up
42!
[1]43! Revision 1.12  2006/02/23 12:30:22  raasch
44! ebene renamed section, pl.. replaced by do..,
45!
46! Revision 1.1  2000/04/27 06:24:39  raasch
47! Initial revision
48!
49!
50! Description:
51! ------------
52! Initializing actions needed when using dvrp-software
53!------------------------------------------------------------------------------!
54#if defined( __dvrp_graphics )
55
56    USE arrays_3d
57    USE DVRP
58    USE dvrp_variables
59    USE grid_variables
60    USE indices
61    USE pegrid
62    USE control_parameters
63
64    IMPLICIT NONE
65
66    CHARACTER (LEN=2)  ::  section_chr
[155]67    CHARACTER (LEN=3)  ::  prefix_chr
[1]68    CHARACTER (LEN=80) ::  dvrp_file_local
[237]69    INTEGER ::  cluster_mode, cluster_size_x, cluster_size_y, cluster_size_z, &
[242]70                gradient_normals, i, j, k, l, m, nx_dvrp_l, nx_dvrp_r,        &
71                ny_dvrp_n, ny_dvrp_s, pn, tv, vn
[1]72    LOGICAL ::  allocated
[237]73    REAL(4) ::  center(3), cluster_alpha, distance, tmp_b, tmp_g, tmp_r, &
74                tmp_t, tmp_th, tmp_thr, tmp_x1, tmp_x2, tmp_y1, tmp_y2,  &
75                tmp_z1, tmp_z2, tmp_1, tmp_2, tmp_3, tmp_4, tmp_5, tmp_6, tmp_7
[1]76
[210]77    REAL(4), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf
[1]78
79    TYPE(CSTRING), SAVE ::  dvrp_directory_c, dvrp_file_c, &
80                            dvrp_file_local_c,dvrp_host_c, &
81                            dvrp_password_c, dvrp_username_c, name_c
82
83!
[242]84!-- Set clipping to default (total domain), if not set by user
85    IF ( clip_dvrp_l == 9999999.9 )  clip_dvrp_l = 0.0
86    IF ( clip_dvrp_r == 9999999.9 )  clip_dvrp_r = ( nx + 1 ) * dx
87    IF ( clip_dvrp_s == 9999999.9 )  clip_dvrp_s = 0.0
88    IF ( clip_dvrp_n == 9999999.9 )  clip_dvrp_n = ( ny + 1 ) * dy
89
90!
91!-- Calculate the clipping index limits
92    nx_dvrp_l = clip_dvrp_l / dx
93    nx_dvrp_r = clip_dvrp_r / dx
94    ny_dvrp_s = clip_dvrp_s / dy
95    ny_dvrp_n = clip_dvrp_n / dy
96
97    IF ( nx_dvrp_l < nxr  .AND.  nx_dvrp_r > nxl  .AND. &
98         ny_dvrp_s < nyn  .AND.  ny_dvrp_n > nys )  THEN
99
100       dvrp_overlap = .TRUE.
101       nxl_dvrp = MAX( nxl, nx_dvrp_l )
102       nxr_dvrp = MIN( nxr, nx_dvrp_r )
103       nys_dvrp = MAX( nys, ny_dvrp_s )
104       nyn_dvrp = MIN( nyn, ny_dvrp_n )
105
106       IF ( nxl_dvrp == nxl  .AND.  nxr_dvrp == nxr  .AND.  &
107            nys_dvrp == nys  .AND.  nyn_dvrp == nyn )  THEN
108          dvrp_total_overlap = .TRUE.
109       ELSE
110          dvrp_total_overlap = .FALSE.
111       ENDIF
112
113    ELSE
114!
115!--    This subdomain does not overlap with the clipping area. Define an
116!--    arbitrary (small) domain within in the clipping area.
117       dvrp_overlap       = .FALSE.
118       dvrp_total_overlap = .FALSE.
[264]119!       nxl_dvrp = nx_dvrp_l
120!       nxr_dvrp = nxl_dvrp + 4
121!       nys_dvrp = ny_dvrp_s
122!       nyn_dvrp = nys_dvrp + 4
123       nxl_dvrp = nxl
124       nxr_dvrp = MIN( nxl+4, nxr )
125       nys_dvrp = nys
126       nyn_dvrp = MIN( nys+4, nyn )
[242]127
128    ENDIF
129
130!
[1]131!-- Set the maximum time the program can be suspended on user request (by
132!-- dvrp steering). This variable is defined in module DVRP.
133    DVRP_MAX_SUSPEND_TIME = 7200
134
135!
136!-- Allocate array holding the names and limits of the steering variables
137!-- (must have the same number of elements as array mode_dvrp!)
138    ALLOCATE( steering_dvrp(10) )
139
140!
141!-- Check, if output parameters are given and/or allowed
142!-- and set default-values, where necessary
143    IF ( dvrp_username == ' ' )  THEN
[254]144        message_string = 'dvrp_username is undefined'
145        CALL message( 'init_dvrp', 'PA0195', 1, 2, 0, 6, 0 )         
[1]146    ENDIF
147
148    IF ( dvrp_output /= 'ftp'  .AND.  dvrp_output /= 'rtsp'  .AND. &
149         dvrp_output /= 'local' )  THEN
[274]150       message_string = 'dvrp_output="' // TRIM( dvrp_output ) // &
151                        '" not allowed'
[254]152       CALL message( 'init_dvrp', 'PA0196', 1, 2, 0, 6, 0 )
[1]153    ENDIF
154
155    IF ( dvrp_directory == 'default' )  THEN
156       dvrp_directory = TRIM( dvrp_username ) // '/' // TRIM( run_identifier )
157    ENDIF
158
[260]159!
160!-- A local dvrserver running always outputs on temporary directory DATA_DVR
161    IF ( local_dvrserver_running )  THEN
162       dvrp_directory = 'DATA_DVR'
163    ENDIF
164
[13]165    IF ( dvrp_output /= 'local' )  THEN
166       IF ( dvrp_file /= 'default'  .AND.  dvrp_file /= '/dev/null' )  THEN
[254]167          message_string = 'dvrp_file="' // TRIM( dvrp_file ) // '" not allowed'
168          CALL message( 'init_dvrp', 'PA0197', 1, 2, 0, 6, 0 )
[1]169       ENDIF
170    ENDIF
171
172!
173!-- Strings are assigned to strings of special type which have a CHAR( 0 )
174!-- (C end-of-character symbol) at their end. This is needed when strings are
175!-- passed to C routines.
176    dvrp_directory_c = dvrp_directory
177    dvrp_file_c      = dvrp_file
178    dvrp_host_c      = dvrp_host
179    dvrp_password_c  = dvrp_password
180    dvrp_username_c  = dvrp_username
181
182!
183!-- Loop over all output modes choosed
184    m = 1
185    allocated = .FALSE.
186    DO WHILE ( mode_dvrp(m) /= ' ' )
187   
188!
189!--    Check, if mode is allowed
190       IF ( mode_dvrp(m)(1:10) /= 'isosurface'  .AND. &
191            mode_dvrp(m)(1:6)  /= 'slicer'      .AND. &
[210]192            mode_dvrp(m)(1:9)  /= 'particles'   .AND. &
193            mode_dvrp(m)(1:9)  /= 'pathlines' )  THEN
[1]194
[274]195          message_string = 'mode_dvrp="' // TRIM( mode_dvrp(m) ) // &
196                           '" not allowed'
[254]197          CALL message( 'init_dvrp', 'PA0198', 1, 2, 0, 6, 0 )
[1]198          CALL local_stop
199
200       ENDIF
201!
[155]202!--    Determine prefix for dvrp_file
203       WRITE ( prefix_chr, '(I2.2,''_'')' )  m
204!
[1]205!--    Camera position must be computed and written on file when no dvrp-output
206!--    has been generated so far (in former runs)
207!       IF ( dvrp_filecount == 0 )  THEN
208!
209!--       Compute center of domain and distance of camera from center
[242]210          center(1) = ( clip_dvrp_l + clip_dvrp_r ) * 0.5 * superelevation_x
211          center(2) = ( clip_dvrp_s + clip_dvrp_n ) * 0.5 * superelevation_y
[1]212          center(3) = ( zu(nz_do3d) - zu(nzb) ) * 0.5 * superelevation
[242]213          distance  = 1.5 * MAX( (clip_dvrp_r-clip_dvrp_l) * superelevation_x, &
214                                 (clip_dvrp_n-clip_dvrp_s) * superelevation_y, &
[1]215                                 ( zu(nz_do3d) - zu(nzb) ) * superelevation )
216
217!
218!--       Write camera position on file
219          CALL DVRP_INIT( m-1, 0 )
220
221!
222!--       Create filename for camera
223          IF ( dvrp_output == 'rtsp' )  THEN
224
[155]225             dvrp_file   = prefix_chr // TRIM( mode_dvrp(m) ) // '/camera.dvr'
[1]226             dvrp_file_c = dvrp_file
227             CALL DVRP_OUTPUT_RTSP( m-1, dvrp_host_c, dvrp_username_c, &
228                                    dvrp_password_c, dvrp_directory_c, &
229                                    dvrp_file_c )
230
231          ELSEIF ( dvrp_output == 'ftp' )  THEN
232
[155]233             dvrp_file   = prefix_chr // TRIM( mode_dvrp(m) ) // '.camera.dvr'
[1]234             dvrp_file_c = dvrp_file
235!             CALL DVRP_OUTPUT_FTP( m-1, 0, dvrp_host_c, dvrp_username_c, &
236!                                   dvrp_password_c, dvrp_directory_c,    &
237!                                   dvrp_file_c )
238
239          ELSE
240
241             IF ( dvrp_file(1:9) /= '/dev/null' )  THEN
[155]242                dvrp_file_local   = prefix_chr // TRIM( mode_dvrp(m) )  &
243                     // '.camera.dvr'
[1]244                dvrp_file_local_c = dvrp_file_local
245             ELSE
246                dvrp_file_local_c = dvrp_file_c
247             ENDIF
[13]248             CALL DVRP_OUTPUT_LOCAL( m-1, 0, dvrp_file_local_c )
[1]249
250          ENDIF
251
252          CALL DVRP_CAMERA( m-1, center, distance )
253
254!
255!--       Define bounding box material and create a bounding box
[210]256          tmp_r = 0.5;  tmp_g = 0.5;  tmp_b = 0.5;  tmp_t = 0.0
257          CALL DVRP_MATERIAL_RGB( m-1, 1, tmp_r, tmp_g, tmp_b, tmp_t )
[1]258
[242]259          tmp_1 = 0.01;
260          tmp_2 = clip_dvrp_l * superelevation_x
261          tmp_3 = clip_dvrp_s * superelevation_y
262          tmp_4 = 0.0
[264]263          tmp_5 = (clip_dvrp_r+dx) * superelevation_x
264          tmp_6 = (clip_dvrp_n+dy) * superelevation_y
[210]265          tmp_7 = zu(nz_do3d) * superelevation
266          CALL DVRP_BOUNDINGBOX( m-1, 1, tmp_1, tmp_2, tmp_3, tmp_4, tmp_5, &
267                                 tmp_6, tmp_7 )
268
[1]269          CALL DVRP_VISUALIZE( m-1, 0, 0 )
270          CALL DVRP_EXIT( m-1 )
271
272!
273!--       Write topography isosurface on file
[237]274          IF ( TRIM( topography ) /= 'flat' )  THEN
[1]275
[237]276             CALL DVRP_INIT( m-1, 0 )
277
[1]278!
[260]279!--          Create filename for topography
[237]280             IF ( dvrp_output == 'rtsp' )  THEN
[1]281
[237]282                dvrp_file   = prefix_chr // TRIM( mode_dvrp(m) )  &
[260]283                              // '/topography.dvr'
[237]284                dvrp_file_c = dvrp_file
285                CALL DVRP_OUTPUT_RTSP( m-1, dvrp_host_c, dvrp_username_c, &
286                                       dvrp_password_c, dvrp_directory_c, &
287                                       dvrp_file_c )
[1]288
[237]289             ELSEIF ( dvrp_output == 'ftp' )  THEN
[17]290
[237]291                dvrp_file   = prefix_chr // TRIM( mode_dvrp(m) )  &
[260]292                              // '.topography.dvr'
[237]293                dvrp_file_c = dvrp_file
294!                CALL DVRP_OUTPUT_FTP( m-1, 0, dvrp_host_c, dvrp_username_c, &
295!                                      dvrp_password_c, dvrp_directory_c,    &
296!                                      dvrp_file_c )
[17]297
[237]298             ELSE
[17]299
[237]300                IF ( dvrp_file(1:9) /= '/dev/null' )  THEN
301                   dvrp_file_local   = prefix_chr // TRIM( mode_dvrp(m) )  &
[260]302                                       // '.topography.dvr'
[237]303                   dvrp_file_local_c = dvrp_file_local
304                ELSE
305                   dvrp_file_local_c = dvrp_file_c
306                ENDIF
307                CALL DVRP_OUTPUT_LOCAL( m-1, 0, dvrp_file_local_c )
308
[17]309             ENDIF
310
[237]311!
312!--          Determine local gridpoint coordinates
313             IF ( .NOT. allocated )  THEN
[242]314                ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
315                          ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
[237]316                          zcoor_dvrp(nzb:nz_do3d) )
317                allocated = .TRUE.
[1]318
[242]319                DO  i = nxl_dvrp, nxr_dvrp+1
[237]320                   xcoor_dvrp(i) = i * dx * superelevation_x
321                ENDDO
[242]322                DO  j = nys_dvrp, nyn_dvrp+1
[237]323                   ycoor_dvrp(j) = j * dy * superelevation_y
324                ENDDO
325                zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
[242]326                nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
327                ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
[237]328                nz_dvrp    = nz_do3d - nzb + 1
329             ENDIF
330
[1]331!
[237]332!--          Define the grid used by dvrp
333             CALL DVRP_NO_GLOBAL_GRID( m-1, 1 )
334             CALL DVRP_GRID( m-1, nx_dvrp, ny_dvrp, nz_dvrp, xcoor_dvrp, &
335                             ycoor_dvrp, zcoor_dvrp )
[1]336
[284]337             tmp_r = topography_color(1)
338             tmp_g = topography_color(2)
339             tmp_b = topography_color(3)
340             tmp_t = 0.0
[237]341             CALL DVRP_MATERIAL_RGB( m-1, 1, tmp_r, tmp_g, tmp_b, tmp_t )
342
343!
344!--          Compute and plot isosurface in dvr-format
[242]345             ALLOCATE( local_pf(nxl_dvrp:nxr_dvrp+1,nys_dvrp:nyn_dvrp+1, &
346                                nzb:nz_do3d) )
[237]347             local_pf = 0.0
[242]348             IF ( dvrp_overlap )  THEN
349                DO  i = nxl_dvrp, nxr_dvrp+1
350                   DO  j = nys_dvrp, nyn_dvrp+1
351                      IF ( nzb_s_inner(j,i) > 0 )  THEN
[237]352                         local_pf(i,j,nzb:nzb_s_inner(j,i)) = 1.0
353                      ENDIF
[242]354                   ENDDO
[237]355                ENDDO
[242]356             ENDIF
[1]357
[237]358             CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
359                             cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
360
361             tmp_th = 1.0
362             CALL DVRP_THRESHOLD( m-1, tmp_th )
363
[1]364!
[237]365!--          Reduce the number of polygones, if required
366             IF ( cluster_size > 1 )  THEN
[210]367
[237]368                cluster_size_x = cluster_size
369                cluster_size_y = cluster_size
370                cluster_size_z = cluster_size
371                cluster_mode     = 4    ! vertex clustering mode
372                gradient_normals = 0    ! use flat-shading
[1]373
[237]374                CALL DVRP_CLUSTER_SIZE( m-1, cluster_size_x, cluster_size_y, &
375                                        cluster_size_z )
376                CALL DVRP_CLUSTERING_MODE( m-1, cluster_mode )
377                CALL DVRP_GRADIENTNORMALS( m-1, gradient_normals )
[1]378!
[237]379!--             Set parameter for vertex clustering mode 4.
380!--             ATTENTION: A seperate procedure for setting cluster_alpha will
381!--                        be in the next version of libDVRP (Feb 09)
382                cluster_alpha = 38.0
383                CALL DVRP_THRESHOLD( -(m-1)-1, cluster_alpha )
[1]384
[237]385                CALL DVRP_VISUALIZE( m-1, 21, 0 )
[1]386
[237]387             ELSE
388!
389!--             No polygon reduction
390                CALL DVRP_VISUALIZE( m-1, 1, 0 )
[1]391
[237]392             ENDIF
393
394             DEALLOCATE( local_pf )
395
396             CALL DVRP_EXIT( m-1 )
397
398          ENDIF
399
[1]400!
[260]401!--       Write the ground plate (z=0) isosurface on file
[1]402          CALL DVRP_INIT( m-1, 0 )
403
404!
405!--       Create filename for surface
406          IF ( dvrp_output == 'rtsp' )  THEN
407
[260]408             dvrp_file   = prefix_chr // TRIM( mode_dvrp(m) ) // &
409                           '/groundplate.dvr'
[1]410             dvrp_file_c = dvrp_file
411             CALL DVRP_OUTPUT_RTSP( m-1, dvrp_host_c, dvrp_username_c, &
412                                    dvrp_password_c, dvrp_directory_c, &
413                                    dvrp_file_c )
414
[17]415          ELSEIF ( dvrp_output == 'ftp' )  THEN
416
[260]417             dvrp_file   = prefix_chr // TRIM( mode_dvrp(m) ) // &
418                           '.groundplate.dvr'
[17]419             dvrp_file_c = dvrp_file
420!             CALL DVRP_OUTPUT_FTP( m-1, 0, dvrp_host_c, dvrp_username_c, &
421!                                   dvrp_password_c, dvrp_directory_c,    &
422!                                   dvrp_file_c )
423
424          ELSE
425
426             IF ( dvrp_file(1:9) /= '/dev/null' )  THEN
[155]427                dvrp_file_local   = prefix_chr // TRIM( mode_dvrp(m) )  &
[260]428                     // '.groundplate.dvr'
[17]429                dvrp_file_local_c = dvrp_file_local
430             ELSE
431                dvrp_file_local_c = dvrp_file_c
432             ENDIF
433             CALL DVRP_OUTPUT_LOCAL( m-1, 0, dvrp_file_local_c )
434
[1]435          ENDIF
436
437!
438!--       Determine local gridpoint coordinates
439          IF ( .NOT. allocated )  THEN
[242]440             ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
441                       ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
[1]442                       zcoor_dvrp(nzb:nz_do3d) )
443             allocated = .TRUE.
444
[242]445             DO  i = nxl_dvrp, nxr_dvrp+1
[1]446                xcoor_dvrp(i) = i * dx * superelevation_x
447             ENDDO
[242]448             DO  j = nys_dvrp, nyn_dvrp+1
[1]449                ycoor_dvrp(j) = j * dy * superelevation_y
450             ENDDO
451             zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
[242]452             nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
453             ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
[1]454             nz_dvrp    = nz_do3d - nzb + 1
455          ENDIF
456
457!
458!--       Define the grid used by dvrp
[210]459          CALL DVRP_NO_GLOBAL_GRID( m-1, 1 )
[1]460          CALL DVRP_GRID( m-1, nx_dvrp, ny_dvrp, nz_dvrp, xcoor_dvrp, &
461                          ycoor_dvrp, zcoor_dvrp )
[210]462
[284]463          tmp_r = groundplate_color(1)
464          tmp_g = groundplate_color(2)
465          tmp_b = groundplate_color(3)
466          tmp_t = 0.0
[210]467          CALL DVRP_MATERIAL_RGB( m-1, 1, tmp_r, tmp_g, tmp_b, tmp_t )
[1]468
469!
470!--       Compute and plot isosurface in dvr-format
[242]471          ALLOCATE( local_pf(nxl_dvrp:nxr_dvrp+1,nys_dvrp:nyn_dvrp+1, &
472                             nzb:nz_do3d) )
[1]473          local_pf = 0.0
[242]474          IF (dvrp_overlap )  local_pf(:,:,0) = 1.0
[1]475
476          CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
477                          cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
[210]478          tmp_th = 1.0
479          CALL DVRP_THRESHOLD( m-1, tmp_th )
[1]480
[237]481!
[242]482!--       Always reduce the number of polygones as much as possible
[237]483          cluster_size_x = 5
484          cluster_size_y = 5
485          cluster_size_z = 5
486          cluster_mode     = 4    ! vertex clustering mode
487          gradient_normals = 0    ! use flat-shading
488
489          CALL DVRP_CLUSTER_SIZE( m-1, cluster_size_x, cluster_size_y, &
490                                  cluster_size_z )
491          CALL DVRP_CLUSTERING_MODE( m-1, cluster_mode )
492          CALL DVRP_GRADIENTNORMALS( m-1, gradient_normals )
493!
494!--       Set parameter for vertex clustering mode 4.
495!--       ATTENTION: A seperate procedure for setting cluster_alpha will be in
496!--                  the next version of libDVRP (Feb 09)
497          cluster_alpha = 38.0
498          CALL DVRP_THRESHOLD( -(m-1)-1, cluster_alpha )
499
500          CALL DVRP_VISUALIZE( m-1, 21, 0 )
501
[1]502          DEALLOCATE( local_pf )
503
504          CALL DVRP_EXIT( m-1 )
505   
506!       ENDIF
507
508
509!
510!--    Initialize dvrp for all dvrp-calls during the run
511       CALL DVRP_INIT( m-1, 0 )
512
513!
514!--    Preliminary definition of filename for dvrp-output
515       IF ( dvrp_output == 'rtsp' )  THEN
516
517!
518!--       First initialize parameters for possible interactive steering.
519!--       Every parameter has to be passed to the respective stream.
520          pn = 1
521!
522!--       Initialize threshold counter needed for initialization of the
523!--       isosurface steering variables
524          tv = 0
525
526          DO WHILE ( mode_dvrp(pn) /= ' ' )
527
528             IF ( mode_dvrp(pn)(1:10) == 'isosurface' )  THEN
529
[130]530                READ ( mode_dvrp(pn), '(10X,I2)' )  vn
[1]531                steering_dvrp(pn)%name = do3d(0,vn)
532                tv = tv + 1
533
534                IF ( do3d(0,vn)(1:1) == 'w' )  THEN
535                   steering_dvrp(pn)%min  = -4.0
536                   steering_dvrp(pn)%max  =  5.0
537                ELSE
538                   steering_dvrp(pn)%min  = 288.0
539                   steering_dvrp(pn)%max  = 292.0
540                ENDIF
541
[237]542                name_c  = TRIM( do3d(0,vn) )
[210]543                tmp_thr = threshold(tv)
[1]544                CALL DVRP_STEERING_INIT( m-1, name_c, steering_dvrp(pn)%min, &
[210]545                                         steering_dvrp(pn)%max, tmp_thr )
[1]546
547             ELSEIF ( mode_dvrp(pn)(1:6) == 'slicer' )  THEN
548
[130]549                READ ( mode_dvrp(pn), '(6X,I2)' )  vn
[1]550                steering_dvrp(pn)%name = do2d(0,vn)
551                name_c = TRIM( do2d(0,vn) )
552
553                l = MAX( 2, LEN_TRIM( do2d(0,vn) ) )
554                section_chr = do2d(0,vn)(l-1:l)
555                SELECT CASE ( section_chr )
556                   CASE ( 'xy' )
557                      steering_dvrp(pn)%imin   = 0
558                      steering_dvrp(pn)%imax   = nz_do3d
559                      slicer_position_dvrp(pn) = section(1,1)
560                      CALL DVRP_STEERING_INIT( m-1, name_c,            &
561                                               steering_dvrp(pn)%imin, &
562                                               steering_dvrp(pn)%imax, &
563                                               slicer_position_dvrp(pn) )
564                   CASE ( 'xz' )
565                      steering_dvrp(pn)%imin   = 0
566                      steering_dvrp(pn)%imax   = ny
567                      slicer_position_dvrp(pn) = section(1,2)
568                      CALL DVRP_STEERING_INIT( m-1, name_c,            &
569                                               steering_dvrp(pn)%imin, &
570                                               steering_dvrp(pn)%imax, &
571                                               slicer_position_dvrp(pn) )
572                   CASE ( 'yz' )
573                      steering_dvrp(pn)%imin = 0
574                      steering_dvrp(pn)%imax = nx
575                      slicer_position_dvrp(pn) = section(1,3)
576                      CALL DVRP_STEERING_INIT( m-1, name_c,            &
577                                               steering_dvrp(pn)%imin, &
578                                               steering_dvrp(pn)%imax, &
579                                               slicer_position_dvrp(pn) )
580                END SELECT
581
582             ENDIF
583
584             pn = pn + 1
585
586          ENDDO
587
[155]588          dvrp_file = prefix_chr // TRIM( mode_dvrp(m) ) // '/*****.dvr'
[1]589          dvrp_file_c = dvrp_file
590          CALL DVRP_OUTPUT_RTSP( m-1, dvrp_host_c, dvrp_username_c, &
591                                 dvrp_password_c, dvrp_directory_c, &
592                                 dvrp_file_c )
593
594       ELSEIF ( dvrp_output == 'ftp' )  THEN
595
[155]596          dvrp_file   = prefix_chr // TRIM( mode_dvrp(m) ) // '.%05d.dvr'
[1]597          dvrp_file_c = dvrp_file
598!          CALL DVRP_OUTPUT_FTP( m-1, 0, dvrp_host_c, dvrp_username_c, &
599!                                dvrp_password_c, dvrp_directory_c, dvrp_file_c )
600
601       ELSE
602
603          IF ( dvrp_file(1:9) /= '/dev/null' )  THEN
[155]604             dvrp_file_local   = prefix_chr // TRIM( mode_dvrp(m) )  &
605                  // '_%05d.dvr'
[1]606             dvrp_file_local_c = dvrp_file_local
607          ELSE
608             dvrp_file_local_c = dvrp_file_c
609          ENDIF
[13]610          CALL DVRP_OUTPUT_LOCAL( m-1, 0, dvrp_file_local_c )
[1]611
612       ENDIF
613
614!
615!--    Determine local gridpoint coordinates
616       IF ( .NOT. allocated )  THEN
[242]617          ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
618                    ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
[1]619                    zcoor_dvrp(nzb:nz_do3d) )
620          allocated = .TRUE.
621
[242]622          DO  i = nxl_dvrp, nxr_dvrp+1
[1]623             xcoor_dvrp(i) = i * dx * superelevation_x
624          ENDDO
[242]625          DO  j = nys_dvrp, nyn_dvrp+1
[1]626             ycoor_dvrp(j) = j * dy * superelevation_y
627          ENDDO
628          zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
[242]629          nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
630          ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
[1]631          nz_dvrp    = nz_do3d - nzb + 1
632       ENDIF
633
634!
635!--    Define the grid used by dvrp
[210]636       IF ( mode_dvrp(m) /= 'pathlines' )  THEN
637          CALL DVRP_NO_GLOBAL_GRID( m-1, 1 )
638       ENDIF
[1]639       CALL DVRP_GRID( m-1, nx_dvrp, ny_dvrp, nz_dvrp, xcoor_dvrp, ycoor_dvrp, &
640                       zcoor_dvrp )
[210]641
642       IF ( mode_dvrp(m) == 'pathlines' )  THEN
643
644          tmp_x1 = 0.0;  tmp_y1 = 0.0;  tmp_z1 = 0.0
645          tmp_x2 = 1.0;  tmp_y2 = 1.0;  tmp_z2 = 0.3
646          CALL DVRP_CUBIC_SEEDING( m-1, tmp_x1, tmp_y1, tmp_z1, tmp_x2, tmp_y2,&
647                                   tmp_z2, pathlines_linecount, 2, 0 )
648!
649!--       Set wavecount and wavetime
650          CALL DVRP_PATHLINES_BEHAVIOUR_WAVE( m-1, pathlines_wavecount, &
651                                              pathlines_wavetime,       &
652                                              pathlines_fadeintime,     &
653                                              pathlines_fadeouttime )
654!
655!--       Set pathline length
656          CALL DVRP_PATHLINES_SETMAXHISTORY( m-1, pathlines_maxhistory )
657          CALL DVRP_PATHLINES_SETFADING( m-1, 1, 0.0 )
658
659          CALL DVRP_INIT_PATHLINES( m-1, 0 )
660
661       ENDIF
662
[271]663       IF ( mode_dvrp(m)(1:9) == 'particles' )  THEN
[264]664!
665!--       Define a default colourtable for particles
666          DO  i = 1, 11
667             interval_values_dvrp_prt(1,i) = i - 1.0
668             interval_values_dvrp_prt(2,i) = REAL( i )
669             interval_h_dvrp_prt(:,i) = 270.0 - ( i - 1.0 ) * 9.0
670          ENDDO
671
[271]672          DO  i = 12, 22
[264]673             interval_values_dvrp_prt(1,i) = i - 1.0
674             interval_values_dvrp_prt(2,i) = REAL( i )
675             interval_h_dvrp_prt(:,i) = 70.0 - ( i - 12.0 ) * 9.5
676          ENDDO
677
678          dvrp_colortable_entries_prt = 22
679
680       ENDIF
681
[1]682       m = m + 1
683
684    ENDDO
685
686#endif
687 END SUBROUTINE init_dvrp
688
689 
690 SUBROUTINE init_dvrp_logging
691
692!------------------------------------------------------------------------------!
693! Description:
694! ------------
695! Initializes logging events for time measurement with dvrp software
696! and splits one PE from the global communicator in case that dvrp output
697! shall be done by one single PE.
698!------------------------------------------------------------------------------!
699#if defined( __dvrp_graphics )
700
[210]701    USE control_parameters
[1]702    USE dvrp_variables
703    USE pegrid
704
705    IMPLICIT NONE
706
707    CHARACTER (LEN=4) ::  chr
708    INTEGER           ::  idummy
709
710!
711!-- Initialize logging of calls by DVRP graphic software
712    CALL DVRP_LOG_INIT( 'DVRP_LOG' // CHAR( 0 ), 0 )
713
714!
715!-- User-defined logging events: #1 (total time needed by PALM)
716    CALL DVRP_LOG_SYMBOL( 1, 'PALM_total' // CHAR( 0 ) )
717    CALL DVRP_LOG_SYMBOL( 2, 'PALM_timestep' // CHAR( 0 ) )
718    CALL DVRP_LOG_EVENT( 1, 1 )
719
720#if defined( __parallel )
721!
722!-- Find out, if dvrp output shall be done by a dedicated PE
723    CALL local_getenv( 'use_seperate_pe_for_dvrp_output', 31, chr, idummy )
724    IF ( chr == 'true' )  THEN
[237]725
[1]726       use_seperate_pe_for_dvrp_output = .TRUE.
[206]727
728!
[210]729!--    Adjustment for new MPI-1 coupling. This might be unnecessary.
[206]730#if defined( __mpi2 )
[1]731       CALL DVRP_SPLIT( MPI_COMM_WORLD, comm_palm )
[206]732#else
[210]733       IF ( coupling_mode /= 'uncoupled' ) THEN
[254]734          message_string = 'split of communicator not realized with' // &
[210]735                          ' MPI1 coupling atmosphere-ocean'
[254]736          CALL message( 'init_dvrp_logging', 'PA0199', 1, 2, 0, 6, 0 )
737 
[272]738          CALL DVRP_SPLIT( comm_inter, comm_palm )
[210]739       ELSE
740          CALL DVRP_SPLIT( MPI_COMM_WORLD, comm_palm )
741       ENDIF
[206]742#endif
743
[1]744       CALL MPI_COMM_SIZE( comm_palm, numprocs, ierr )
[237]745
[1]746    ENDIF
747#endif
748
749#endif
750 END SUBROUTINE init_dvrp_logging
751
752
753 SUBROUTINE close_dvrp
754
755!------------------------------------------------------------------------------!
756! Description:
757! ------------
758! Exit of dvrp software and finish dvrp logging
759!------------------------------------------------------------------------------!
760#if defined( __dvrp_graphics )
761
762    USE control_parameters
763    USE dvrp
764    USE dvrp_variables
765
766    INTEGER ::  m
767
768!
769!-- If required, close dvrp-software and logging of dvrp-calls
770    IF ( dt_dvrp /= 9999999.9 )  THEN
771       m = 1
772       DO WHILE ( mode_dvrp(m) /= ' ' )
773          CALL DVRP_EXIT( m-1 )
774          m = m + 1
775       ENDDO
776       CALL DVRP_LOG_EVENT( -1, 1 )   ! Logging of total cpu-time used by PALM
777       IF ( use_seperate_pe_for_dvrp_output )  THEN
[273]778#ifndef __nec
[1]779          CALL DVRP_SPLIT_EXIT( 1 )      ! Argument 0: reduced output
[272]780#endif
[1]781       ELSE
782          CALL DVRP_LOG_EXIT( 1 )        ! Argument 0: reduced output
783       ENDIF
784    ENDIF
785
786#endif
787 END SUBROUTINE close_dvrp
Note: See TracBrowser for help on using the repository browser.