source: palm/tags/release-3.2b/SOURCE/data_output_3d.f90 @ 237

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