source: palm/trunk/SOURCE/data_output_3d.f90 @ 438

Last change on this file since 438 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: 12.5 KB
RevLine 
[1]1 SUBROUTINE data_output_3d( av )
2
3!------------------------------------------------------------------------------!
[254]4! Current revisions:
[1]5! -----------------
[392]6!
7!
8! Former revisions:
9! -----------------
10! $Id: data_output_3d.f90 392 2009-09-24 10:39:14Z raasch $
11!
12! 355 2009-07-17 01:03:01Z letzel
[291]13! simulated_time in NetCDF output replaced by time_since_reference_point.
[263]14! Output of NetCDF messages with aid of message handling routine.
[254]15! Output of messages replaced by message handling routine.
[355]16! Bugfix: to_be_resorted => s_av for time-averaged scalars
[1]17!
[98]18! 96 2007-06-04 08:07:41Z raasch
19! Output of density and salinity
20!
[77]21! 75 2007-03-22 09:54:05Z raasch
22! 2nd+3rd argument removed from exchange horiz
23!
[3]24! RCS Log replace by Id keyword, revision history cleaned up
25!
[1]26! Revision 1.3  2006/06/02 15:18:59  raasch
27! +argument "found", -argument grid in call of routine user_data_output_3d
28!
29! Revision 1.2  2006/02/23 10:23:07  raasch
30! Former subroutine plot_3d renamed data_output_3d, pl.. renamed do..,
31! .._anz renamed .._n,
32! output extended to (almost) all quantities, output of user-defined quantities
33!
34! Revision 1.1  1997/09/03 06:29:36  raasch
35! Initial revision
36!
37!
38! Description:
39! ------------
40! Output of the 3D-arrays in NetCDF and/or AVS format.
41!------------------------------------------------------------------------------!
42
43    USE array_kind
44    USE arrays_3d
45    USE averaging
46    USE cloud_parameters
47    USE control_parameters
48    USE cpulog
49    USE indices
50    USE interfaces
51    USE netcdf_control
52    USE particle_attributes
53    USE pegrid
54
55    IMPLICIT NONE
56
57    CHARACTER (LEN=9) ::  simulated_time_mod
58
59    INTEGER           ::  av, i, if, j, k, n, pos, prec, psi
60
61    LOGICAL           ::  found, resorted
62
63    REAL              ::  mean_r, s_r3, s_r4
64
65    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE  ::  local_pf
66
67    REAL, DIMENSION(:,:,:), POINTER ::  to_be_resorted
68
69!
70!-- Return, if nothing to output
71    IF ( do3d_no(av) == 0 )  RETURN
72
73    CALL cpu_log (log_point(14),'data_output_3d','start')
74
75!
76!-- Open output file.
77!-- Also creates coordinate and fld-file for AVS.
78!-- In case of a run on more than one PE, each PE opens its own file and
79!-- writes the data of its subdomain in binary format (regardless of the format
80!-- the user has requested). After the run, these files are combined to one
81!-- file by combine_plot_fields in the format requested by the user (netcdf
82!-- and/or avs).
83    IF ( avs_output  .OR.  ( numprocs > 1 ) )  CALL check_open( 30 )
84
85#if defined( __netcdf )
86    IF ( myid == 0  .AND.  netcdf_output )  CALL check_open( 106+av*10 )
87#endif
88
89!
90!-- Allocate a temporary array with the desired output dimensions.
91    ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nz_do3d) )
92
93!
94!-- Update the NetCDF time axis
95#if defined( __netcdf )
96    do3d_time_count(av) = do3d_time_count(av) + 1
97    IF ( myid == 0  .AND.  netcdf_output )  THEN
98       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av), &
[291]99                               (/ time_since_reference_point /),  &
[1]100                               start = (/ do3d_time_count(av) /), &
101                               count = (/ 1 /) )
[314]102       CALL handle_netcdf_error( 'data_output_3d', 376 )
[1]103    ENDIF
104#endif
105
106!
107!-- Loop over all variables to be written.
108    if = 1
109
110    DO  WHILE ( do3d(av,if)(1:1) /= ' ' )
111!
112!--    Set the precision for data compression.
113       IF ( do3d_compress )  THEN
114          DO  i = 1, 100
115             IF ( plot_3d_precision(i)%variable == do3d(av,if) )  THEN
116                prec = plot_3d_precision(i)%precision
117                EXIT
118             ENDIF
119          ENDDO
120       ENDIF
121
122!
123!--    Store the array chosen on the temporary array.
124       resorted = .FALSE.
125       SELECT CASE ( TRIM( do3d(av,if) ) )
126
127          CASE ( 'e' )
128             IF ( av == 0 )  THEN
129                to_be_resorted => e
130             ELSE
131                to_be_resorted => e_av
132             ENDIF
133
134          CASE ( 'p' )
135             IF ( av == 0 )  THEN
136                to_be_resorted => p
137             ELSE
138                to_be_resorted => p_av
139             ENDIF
140
141          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
142             IF ( av == 0 )  THEN
143                tend = prt_count
[75]144                CALL exchange_horiz( tend )
[1]145                DO  i = nxl-1, nxr+1
146                   DO  j = nys-1, nyn+1
147                      DO  k = nzb, nz_do3d
148                         local_pf(i,j,k) = tend(k,j,i)
149                      ENDDO
150                   ENDDO
151                ENDDO
152                resorted = .TRUE.
153             ELSE
[75]154                CALL exchange_horiz( pc_av )
[1]155                to_be_resorted => pc_av
156             ENDIF
157
158          CASE ( 'pr' )  ! mean particle radius
159             IF ( av == 0 )  THEN
160                DO  i = nxl, nxr
161                   DO  j = nys, nyn
162                      DO  k = nzb, nzt+1
163                         psi = prt_start_index(k,j,i)
164                         s_r3 = 0.0
165                         s_r4 = 0.0
166                         DO  n = psi, psi+prt_count(k,j,i)-1
167                            s_r3 = s_r3 + particles(n)%radius**3
168                            s_r4 = s_r4 + particles(n)%radius**4
169                         ENDDO
170                         IF ( s_r3 /= 0.0 )  THEN
171                            mean_r = s_r4 / s_r3
172                         ELSE
173                            mean_r = 0.0
174                         ENDIF
175                         tend(k,j,i) = mean_r
176                      ENDDO
177                   ENDDO
178                ENDDO
[75]179                CALL exchange_horiz( tend )
[1]180                DO  i = nxl-1, nxr+1
181                   DO  j = nys-1, nyn+1
182                      DO  k = nzb, nzt+1
183                         local_pf(i,j,k) = tend(k,j,i)
184                      ENDDO
185                   ENDDO
186                ENDDO
187                resorted = .TRUE.
188             ELSE
[75]189                CALL exchange_horiz( pr_av )
[1]190                to_be_resorted => pr_av
191             ENDIF
192
193          CASE ( 'pt' )
194             IF ( av == 0 )  THEN
195                IF ( .NOT. cloud_physics ) THEN
196                   to_be_resorted => pt
197                ELSE
198                   DO  i = nxl-1, nxr+1
199                      DO  j = nys-1, nyn+1
200                         DO  k = nzb, nz_do3d
201                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
202                                                          pt_d_t(k) * &
203                                                          ql(k,j,i)
204                         ENDDO
205                      ENDDO
206                   ENDDO
207                   resorted = .TRUE.
208                ENDIF
209             ELSE
210                to_be_resorted => pt_av
211             ENDIF
212
213          CASE ( 'q' )
214             IF ( av == 0 )  THEN
215                to_be_resorted => q
216             ELSE
217                to_be_resorted => q_av
218             ENDIF
219             
220          CASE ( 'ql' )
221             IF ( av == 0 )  THEN
222                to_be_resorted => ql
223             ELSE
224                to_be_resorted => ql_av
225             ENDIF
226
227          CASE ( 'ql_c' )
228             IF ( av == 0 )  THEN
229                to_be_resorted => ql_c
230             ELSE
231                to_be_resorted => ql_c_av
232             ENDIF
233
234          CASE ( 'ql_v' )
235             IF ( av == 0 )  THEN
236                to_be_resorted => ql_v
237             ELSE
238                to_be_resorted => ql_v_av
239             ENDIF
240
241          CASE ( 'ql_vp' )
242             IF ( av == 0 )  THEN
243                to_be_resorted => ql_vp
244             ELSE
245                to_be_resorted => ql_vp_av
246             ENDIF
247
248          CASE ( 'qv' )
249             IF ( av == 0 )  THEN
250                DO  i = nxl-1, nxr+1
251                   DO  j = nys-1, nyn+1
252                      DO  k = nzb, nz_do3d
253                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
254                      ENDDO
255                   ENDDO
256                ENDDO
257                resorted = .TRUE.
258             ELSE
259                to_be_resorted => qv_av
260             ENDIF
261
[96]262          CASE ( 'rho' )
263             IF ( av == 0 )  THEN
264                to_be_resorted => rho
265             ELSE
266                to_be_resorted => rho_av
267             ENDIF
268             
[1]269          CASE ( 's' )
270             IF ( av == 0 )  THEN
271                to_be_resorted => q
272             ELSE
[355]273                to_be_resorted => s_av
[1]274             ENDIF
275             
[96]276          CASE ( 'sa' )
277             IF ( av == 0 )  THEN
278                to_be_resorted => sa
279             ELSE
280                to_be_resorted => sa_av
281             ENDIF
282             
[1]283          CASE ( 'u' )
284             IF ( av == 0 )  THEN
285                to_be_resorted => u
286             ELSE
287                to_be_resorted => u_av
288             ENDIF
289
290          CASE ( 'v' )
291             IF ( av == 0 )  THEN
292                to_be_resorted => v
293             ELSE
294                to_be_resorted => v_av
295             ENDIF
296
297          CASE ( 'vpt' )
298             IF ( av == 0 )  THEN
299                to_be_resorted => vpt
300             ELSE
301                to_be_resorted => vpt_av
302             ENDIF
303
304          CASE ( 'w' )
305             IF ( av == 0 )  THEN
306                to_be_resorted => w
307             ELSE
308                to_be_resorted => w_av
309             ENDIF
310
311          CASE DEFAULT
312!
313!--          User defined quantity
314             CALL user_data_output_3d( av, do3d(av,if), found, local_pf, &
315                                       nz_do3d )
316             resorted = .TRUE.
317
[254]318             IF ( .NOT. found )  THEN
[274]319                message_string =  'no output available for: ' //   &
320                                  TRIM( do3d(av,if) )
[254]321                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
[1]322             ENDIF
323
324       END SELECT
325
326!
327!--    Resort the array to be output, if not done above
328       IF ( .NOT. resorted )  THEN
329          DO  i = nxl-1, nxr+1
330             DO  j = nys-1, nyn+1
331                DO  k = nzb, nz_do3d
332                   local_pf(i,j,k) = to_be_resorted(k,j,i)
333                ENDDO
334             ENDDO
335          ENDDO
336       ENDIF
337
338!
339!--    Output of the volume data information for the AVS-FLD-file.
340       do3d_avs_n = do3d_avs_n + 1
341       IF ( myid == 0  .AND.  avs_output )  THEN
342!
343!--       AVS-labels must not contain any colons. Hence they must be removed
344!--       from the time character string.
345          simulated_time_mod = simulated_time_chr
346          DO  WHILE ( SCAN( simulated_time_mod, ':' ) /= 0 )
347             pos = SCAN( simulated_time_mod, ':' )
348             simulated_time_mod(pos:pos) = '/'
349          ENDDO
350
351          IF ( av == 0 )  THEN
352             WRITE ( 33, 3300 )  do3d_avs_n, TRIM( avs_data_file ), &
353                                 skip_do_avs, TRIM( do3d(av,if) ), &
354                                 TRIM( simulated_time_mod )
355          ELSE
356             WRITE ( 33, 3300 )  do3d_avs_n, TRIM( avs_data_file ), &
357                                 skip_do_avs, TRIM( do3d(av,if) ) // &
358                                 ' averaged', TRIM( simulated_time_mod )
359          ENDIF
360!
361!--       Determine the Skip-value for the next array. Record end and start
362!--       require 4 byte each.
363          skip_do_avs = skip_do_avs + ( ((nx+2)*(ny+2)*(nz_do3d+1)) * 4 + 8 )
364       ENDIF
365
366!
367!--    Output of the 3D-array. (compressed/uncompressed)
368       IF ( do3d_compress )  THEN
369!
370!--       Compression, output of compression information on FLD-file and output
371!--       of compressed data.
372          CALL write_compressed( local_pf, 30, 33, myid, nxl, nxr, nyn, nys, &
373                                 nzb, nz_do3d, prec )
374       ELSE
375!
376!--       Uncompressed output.
377#if defined( __parallel )
378          IF ( netcdf_output  .AND.  myid == 0 )  THEN
379             WRITE ( 30 )  simulated_time, do3d_time_count(av), av
380          ENDIF
381          WRITE ( 30 )  nxl-1, nxr+1, nys-1, nyn+1, nzb, nz_do3d
382          WRITE ( 30 )  local_pf
383#else
384          IF ( avs_output )  THEN
385             WRITE ( 30 )  local_pf(nxl:nxr+1,nys:nyn+1,:)
386          ENDIF
387#if defined( __netcdf )
388          IF ( netcdf_output )  THEN
389             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),       &
390                                  local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
391                                  start = (/ 1, 1, 1, do3d_time_count(av) /), &
392                                  count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
[314]393             CALL handle_netcdf_error( 'data_output_3d', 386 )
[1]394          ENDIF
395#endif
396#endif
397       ENDIF
398
399       if = if + 1
400
401    ENDDO
402
403!
404!-- Deallocate temporary array.
405    DEALLOCATE( local_pf )
406
407
408    CALL cpu_log (log_point(14),'data_output_3d','stop','nobarrier')
409
410!
411!-- Formats.
4123300 FORMAT ('variable ',I4,'  file=',A,'  filetype=unformatted  skip=',I12/ &
413             'label = ',A,A)
414
415 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.