source: palm/tags/release-3.1c/SOURCE/data_output_3d.f90 @ 710

Last change on this file since 710 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

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