source: palm/tags/release-3.4/SOURCE/data_output_dvrp.f90 @ 4404

Last change on this file since 4404 was 106, checked in by raasch, 17 years ago

preliminary update of bugfixes and extensions for non-cyclic BCs

  • Property svn:keywords set to Id
File size: 18.5 KB
Line 
1 MODULE dvrp_color
2
3    USE dvrp_variables
4
5    IMPLICIT NONE
6
7 CONTAINS
8
9    SUBROUTINE color_dvrp( value, color )
10
11       REAL, INTENT(IN)  ::  value
12       REAL, INTENT(OUT) ::  color(4)
13
14       REAL              ::  scale
15
16       scale = ( value - slicer_range_limits_dvrp(1,islice_dvrp) ) / &
17               ( slicer_range_limits_dvrp(2,islice_dvrp) -           &
18                 slicer_range_limits_dvrp(1,islice_dvrp) )
19
20       scale = MODULO( 180.0 + 180.0 * scale, 360.0 )
21
22       color = (/ scale, 0.5, 1.0, 0.0 /)
23
24    END SUBROUTINE color_dvrp
25
26 END MODULE dvrp_color
27
28
29 RECURSIVE SUBROUTINE data_output_dvrp
30
31!------------------------------------------------------------------------------!
32! Actual revisions:
33! -----------------
34! TEST: different colours for isosurfaces
35! TEST: write statements
36!
37! Former revisions:
38! -----------------
39! $Id: data_output_dvrp.f90 106 2007-08-16 14:30:26Z suehring $
40!
41! 82 2007-04-16 15:40:52Z raasch
42! Preprocessor strings for different linux clusters changed to "lc",
43! routine local_flush is used for buffer flushing
44!
45! 75 2007-03-22 09:54:05Z raasch
46! Particles-package is now part of the default code,
47! moisture renamed humidity
48!
49! RCS Log replace by Id keyword, revision history cleaned up
50!
51! Revision 1.13  2006/02/23 10:25:12  raasch
52! Former routine plot_dvrp renamed data_output_dvrp,
53! Only a fraction of the particles may have a tail,
54! pl.. replaced by do.., %size renamed %dvrp_psize
55!
56! Revision 1.1  2000/04/27 06:27:17  raasch
57! Initial revision
58!
59!
60! Description:
61! ------------
62! Plot of isosurface, particles and slicers with dvrp-software
63!------------------------------------------------------------------------------!
64#if defined( __dvrp_graphics )
65
66    USE arrays_3d
67    USE cloud_parameters
68    USE cpulog
69    USE DVRP
70    USE dvrp_color
71    USE dvrp_variables
72    USE grid_variables
73    USE indices
74    USE interfaces
75    USE particle_attributes
76    USE pegrid
77    USE control_parameters
78
79    IMPLICIT NONE
80
81    CHARACTER (LEN=2) ::  section_chr
82    CHARACTER (LEN=6) ::  output_variable
83    INTEGER ::  i, j, k, l, m, n, nn, section_mode, tv, vn
84    INTEGER, DIMENSION(:), ALLOCATABLE ::  p_c, p_t
85    REAL    ::  center(3), distance, slicer_position, surface_value
86    REAL, DIMENSION(:),     ALLOCATABLE ::  psize, p_x, p_y, p_z
87    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  local_pf
88
89
90    WRITE ( 9, * ) '*** myid=', myid, ' Anfang data_output_dvrp'
91    CALL local_flush( 9 )
92    CALL cpu_log( log_point(27), 'data_output_dvrp', 'start' )
93
94!
95!-- Loop over all output modes choosed
96    m           = 1
97    tv          = 0  ! threshold counter
98    islice_dvrp = 0  ! slice plane counter
99    DO WHILE ( mode_dvrp(m) /= ' ' )
100!
101!--    Update of the steering variables
102       IF ( .NOT. lock_steering_update )  THEN
103!
104!--       Set lock to avoid recursive calls of DVRP_STEERING_UPDATE
105          lock_steering_update = .TRUE.
106!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: vor steering_update'
107!   CALL local_flush( 9 )
108          CALL DVRP_STEERING_UPDATE( m-1, data_output_dvrp )
109!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: nach steering_update'
110!   CALL local_flush( 9 )
111          lock_steering_update = .FALSE.
112       ENDIF
113
114!
115!--    Determine the variable which shall be plotted (in case of slicers or
116!--    isosurfaces)
117       IF ( mode_dvrp(m)(1:10) == 'isosurface' )  THEN
118          READ ( mode_dvrp(m), '(10X,I1)' )  vn
119          output_variable = do3d(0,vn)
120          tv = tv + 1
121       ELSEIF ( mode_dvrp(m)(1:6) == 'slicer' )  THEN
122          READ ( mode_dvrp(m), '(6X,I1)' )  vn
123          output_variable = do2d(0,vn)
124          l = MAX( 2, LEN_TRIM( do2d(0,vn) ) )
125          section_chr = do2d(0,vn)(l-1:l)
126          SELECT CASE ( section_chr )
127             CASE ( 'xy' )
128                section_mode = 2
129                slicer_position = zu(MIN( slicer_position_dvrp(m), nz_do3d ))
130             CASE ( 'xz' )
131                section_mode = 1
132                slicer_position = slicer_position_dvrp(m) * dy
133             CASE ( 'yz' )
134                section_mode = 0
135                slicer_position = slicer_position_dvrp(m) * dx
136          END SELECT
137       ENDIF
138
139!
140!--    Select the plot mode (in case of isosurface or slicer only if user has
141!--    defined a variable which shall be plotted; otherwise do nothing)
142       IF ( mode_dvrp(m)(1:9) == 'particles'  .AND.  particle_advection  .AND. &
143            simulated_time >= particle_advection_start )  THEN
144
145!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: anfang particles'
146!   CALL local_flush( 9 )
147!
148!--       DVRP-Calls for plotting particles:
149          CALL cpu_log( log_point_s(28), 'dvrp_particles', 'start' )
150
151!
152!--       Definition of characteristics of particle material
153!          CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.7, 0.1, 0.0 )
154          CALL DVRP_MATERIAL_RGB( m-1, 1, 0.0, 0.0, 0.0, 0.0 )
155
156!
157!--       Move particle coordinates to one-dimensional arrays
158          IF ( .NOT. use_particle_tails )  THEN
159!
160!--          All particles are output
161             ALLOCATE( psize(number_of_particles), p_t(number_of_particles), &
162                       p_c(number_of_particles), p_x(number_of_particles),   &
163                       p_y(number_of_particles), p_z(number_of_particles) )
164             psize = 0.0;  p_t = 0;  p_c = 0.0;  p_x = 0.0;  p_y = 0.0
165             p_z   = 0.0;
166             psize = particles(1:number_of_particles)%dvrp_psize
167             p_x   = particles(1:number_of_particles)%x * superelevation_x
168             p_y   = particles(1:number_of_particles)%y * superelevation_y
169             p_z   = particles(1:number_of_particles)%z * superelevation
170             p_c   = particles(1:number_of_particles)%color
171          ELSE
172!
173!--          Particles have a tail
174!            WRITE (9,*) '--- before ALLOCATE  simtime=',simulated_time,' #of_tails=', number_of_tails, &
175!                          ' max#of_tp=', maximum_number_of_tailpoints
176!   CALL local_flush( 9 )
177             ALLOCATE( psize(number_of_tails), p_t(number_of_tails),      &
178                       p_c(number_of_tails*maximum_number_of_tailpoints), &
179                       p_x(number_of_tails*maximum_number_of_tailpoints), &
180                       p_y(number_of_tails*maximum_number_of_tailpoints), &
181                       p_z(number_of_tails*maximum_number_of_tailpoints) )
182!            WRITE (9,*) '--- after ALLOCATE'
183!   CALL local_flush( 9 )
184             psize = 0.0;  p_t = 0;  p_c = 0.0;  p_x = 0.0;  p_y = 0.0
185             p_z   = 0.0;
186             i = 0
187             k = 0
188             DO  n = 1, number_of_particles
189                nn = particles(n)%tail_id
190                IF ( nn /= 0 )  THEN
191                   k = k + 1
192!                  IF ( simulated_time > 1338.0 )  THEN
193!                     WRITE (9,*) '--- particle ',n,' tail_id=',nn,' #of_tp=',particles(n)%tailpoints
194!   CALL local_flush( 9 )
195!                  ENDIF
196                   DO  j = 1, particles(n)%tailpoints
197                      i = i + 1
198                      p_x(i) = particle_tail_coordinates(j,1,nn) * &
199                                                                superelevation_x
200                      p_y(i) = particle_tail_coordinates(j,2,nn) * &
201                                                                superelevation_y
202                      p_z(i) = particle_tail_coordinates(j,3,nn) * &
203                                                                superelevation
204                      p_c(i) = particle_tail_coordinates(j,4,nn)
205!                     IF ( simulated_time > 1338.0 )  THEN
206!                        WRITE (9,*) '--- tp= ',i,' x=',p_x(i),' y=',p_y(i), &
207!                                                ' z=',p_z(i),' c=',p_c(i)
208!   CALL local_flush( 9 )
209!                     ENDIF
210                   ENDDO
211                   psize(k) = particles(n)%dvrp_psize
212                   p_t(k)   = particles(n)%tailpoints - 1
213!                  IF ( simulated_time > 1338.0 )  THEN
214!                     WRITE (9,*) '--- t= ',k,' psize=',psize(k),' p_t=',p_t(k)
215!   CALL local_flush( 9 )
216!                  ENDIF
217                ENDIF               
218             ENDDO
219!            WRITE (9,*) '--- after locally storing the particle attributes'
220!   CALL local_flush( 9 )
221          ENDIF
222
223!
224!--       Compute and plot particles in dvr-format
225          IF ( uniform_particles  .AND.  .NOT. use_particle_tails )  THEN
226!
227!--          All particles have the same color. Use simple routine to set
228!--          the particle attributes (produces less output data)
229             CALL DVRP_PARTICLES( m-1, p_x, p_y, p_z, psize )
230          ELSE
231!
232!--          Set color definitions
233             CALL user_dvrp_coltab( 'particles', 'none' )
234
235             CALL DVRP_COLORTABLE_HLS( m-1, 0, interval_values_dvrp,     &
236                                       interval_h_dvrp, interval_l_dvrp, &
237                                       interval_s_dvrp, interval_a_dvrp )
238
239             IF ( .NOT. use_particle_tails )  THEN
240                CALL DVRP_PARTICLES( m-1, number_of_particles, p_x, p_y, p_z, &
241                                     3, psize, p_c, p_t )
242             ELSE
243!               WRITE (9,*) '--- before DVRP_PARTICLES'
244!   CALL local_flush( 9 )
245                CALL DVRP_PARTICLES( m-1, number_of_tails, p_x, p_y, p_z, 15, &
246                                     psize, p_c, p_t )
247!               WRITE (9,*) '--- after DVRP_PARTICLES'
248!               WRITE (9,*) 'm-1 = ',m-1
249!               WRITE (9,*) 'number_of_tails=', number_of_tails
250!               WRITE (9,*) 'p_x =', p_x
251!               WRITE (9,*) 'p_y =', p_y
252!               WRITE (9,*) 'p_z =', p_z
253!               WRITE (9,*) 'psize =', psize
254!               WRITE (9,*) 'p_c =', p_c
255!               WRITE (9,*) 'p_t =', p_t
256
257!   CALL local_flush( 9 )
258             ENDIF
259          ENDIF
260
261          CALL DVRP_VISUALIZE( m-1, 3, dvrp_filecount )
262!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: ende particles'
263!   CALL local_flush( 9 )
264
265          DEALLOCATE( psize, p_c, p_t, p_x, p_y, p_z )
266
267          CALL cpu_log( log_point_s(28), 'dvrp_particles', 'stop' )
268
269
270       ELSEIF ( ( mode_dvrp(m)(1:10) == 'isosurface'  .OR.   &
271                  mode_dvrp(m)(1:6)  == 'slicer'           ) &
272                  .AND.  output_variable /= ' ' )  THEN
273
274!
275!--       Create an intermediate array, properly dimensioned for plot-output
276          ALLOCATE( local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d) )
277
278!
279!--       Move original array to intermediate array
280          SELECT CASE ( output_variable )
281
282             CASE ( 'u', 'u_xy', 'u_xz', 'u_yz' )
283                DO  i = nxl, nxr+1
284                   DO  j = nys, nyn+1
285                      DO  k = nzb, nz_do3d
286                         local_pf(i,j,k) = u(k,j,i)
287                      ENDDO
288                   ENDDO
289                ENDDO
290!
291!--             Replace mirrored values at lower surface by real surface values
292                IF ( output_variable == 'u_xz'  .OR. &
293                     output_variable == 'u_yz' )  THEN
294                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
295                ENDIF
296
297
298             CASE ( 'v', 'v_xy', 'v_xz', 'v_yz' )
299                DO  i = nxl, nxr+1
300                   DO  j = nys, nyn+1
301                      DO  k = nzb, nz_do3d
302                         local_pf(i,j,k) = v(k,j,i)
303                      ENDDO
304                   ENDDO
305                ENDDO
306!
307!--             Replace mirrored values at lower surface by real surface values
308                IF ( output_variable == 'v_xz'  .OR. &
309                     output_variable == 'v_yz' )  THEN
310                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
311                ENDIF
312
313             CASE ( 'w', 'w_xy', 'w_xz', 'w_yz' )
314                DO  i = nxl, nxr+1
315                   DO  j = nys, nyn+1
316                      DO  k = nzb, nz_do3d
317                         local_pf(i,j,k) = w(k,j,i)
318                      ENDDO
319                   ENDDO
320                ENDDO
321! Averaging for Langmuir circulation
322!                DO  k = nzb, nz_do3d
323!                   DO  j = nys+1, nyn
324!                      DO  i = nxl, nxr+1
325!                         local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + &
326!                                           0.50 * local_pf(i,j,k)   + &
327!                                           0.25 * local_pf(i,j+1,k)
328!                      ENDDO
329!                   ENDDO
330!                ENDDO
331
332             CASE ( 'p', 'p_xy', 'p_xz', 'p_yz' )
333                DO  i = nxl, nxr+1
334                   DO  j = nys, nyn+1
335                      DO  k = nzb, nz_do3d
336                         local_pf(i,j,k) = p(k,j,i)
337                      ENDDO
338                   ENDDO
339                ENDDO
340
341             CASE ( 'pt', 'pt_xy', 'pt_xz', 'pt_yz' )
342                IF ( .NOT. cloud_physics ) THEN
343                   DO  i = nxl, nxr+1
344                      DO  j = nys, nyn+1
345                         DO  k = nzb, nz_do3d
346                            local_pf(i,j,k) = pt(k,j,i)
347                         ENDDO
348                      ENDDO
349                   ENDDO
350                ELSE
351                   DO  i = nxl, nxr+1
352                      DO  j = nys, nyn+1
353                         DO  k = nzb, nz_do3d
354                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp * pt_d_t(k) * &
355                                                          ql(k,j,i)
356                         ENDDO
357                      ENDDO
358                   ENDDO
359                ENDIF
360
361             CASE ( 'q', 'q_xy', 'q_xz', 'q_yz' )
362                IF ( humidity  .OR.  passive_scalar )  THEN
363                   DO  i = nxl, nxr+1
364                      DO  j = nys, nyn+1
365                         DO  k = nzb, nz_do3d
366                            local_pf(i,j,k) = q(k,j,i)
367                         ENDDO
368                      ENDDO
369                   ENDDO           
370                ELSE
371                   IF ( myid == 0 )  THEN
372                      PRINT*, '+++ data_output_dvrp: if humidity/passive_scalar = ', & 
373                              'FALSE output of ', output_variable,            &
374                              'is not provided' 
375                   ENDIF
376                ENDIF
377             
378             CASE ( 'ql', 'ql_xy', 'ql_xz', 'ql_yz' )
379                IF ( cloud_physics  .OR.  cloud_droplets )  THEN
380                   DO  i = nxl, nxr+1
381                      DO  j = nys, nyn+1
382                         DO  k = nzb, nz_do3d
383                            local_pf(i,j,k) = ql(k,j,i)
384                         ENDDO
385                      ENDDO
386                   ENDDO
387                ELSE
388                   IF ( myid == 0 ) THEN
389                      PRINT*, '+++ data_output_dvrp: if cloud_physics = FALSE ', & 
390                              'output of ', output_variable, 'is not provided' 
391                   ENDIF
392                ENDIF
393
394             CASE ( 'u*_xy' )
395                DO  i = nxl, nxr+1
396                   DO  j = nys, nyn+1
397                      local_pf(i,j,nzb+1) = us(j,i)
398                   ENDDO
399                ENDDO
400                slicer_position = zu(nzb+1)
401
402             CASE ( 't*_xy' )
403                DO  i = nxl, nxr+1
404                   DO  j = nys, nyn+1
405                      local_pf(i,j,nzb+1) = ts(j,i)
406                   ENDDO
407                ENDDO
408                slicer_position = zu(nzb+1)
409
410
411             CASE DEFAULT
412                IF ( myid == 0 )  THEN
413                   PRINT*,'+++ data_output_dvrp: no output possible for: ', &
414                          output_variable
415                ENDIF
416
417          END SELECT
418
419
420          IF ( mode_dvrp(m)(1:10) == 'isosurface' )  THEN
421
422!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: anfang isosurface'
423!   CALL local_flush( 9 )
424!
425!--          DVRP-Calls for plotting isosurfaces:
426             CALL cpu_log( log_point_s(26), 'dvrp_isosurface', 'start' )
427
428!
429!--          Definition of characteristics of isosurface material
430!--          Preliminary settings for w!
431             IF ( output_variable == 'w' )  THEN
432                IF ( tv == 1 )  THEN
433                   CALL DVRP_MATERIAL_RGB( m-1, 1, 0.8, 0.1, 0.1, 0.0 )
434                ELSE
435                   CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.1, 0.8, 0.0 )
436                ENDIF
437             ELSE
438                CALL DVRP_MATERIAL_RGB( m-1, 1, 0.9, 0.9, 0.9, 0.0 )
439             ENDIF
440
441!
442!--          Compute and plot isosurface in dvr-format
443             CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
444                             cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
445             CALL DVRP_THRESHOLD( m-1, threshold(tv) )
446             CALL DVRP_VISUALIZE( m-1, 1, dvrp_filecount )
447!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: ende isosurface'
448!   CALL local_flush( 9 )
449
450             CALL cpu_log( log_point_s(26), 'dvrp_isosurface', 'stop' )
451
452          ELSEIF ( mode_dvrp(m)(1:6) == 'slicer' )  THEN
453
454!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: anfang slicer'
455!   CALL local_flush( 9 )
456!
457!--          DVRP-Calls for plotting slicers:
458             CALL cpu_log( log_point_s(27), 'dvrp_slicer', 'start' )
459
460!
461!--          Material and color definitions
462             CALL DVRP_MATERIAL_RGB( m-1, 1, 0.0, 0.0, 0.0, 0.0 )
463
464             islice_dvrp = islice_dvrp + 1
465!             CALL DVRP_COLORFUNCTION( m-1, DVRP_CM_HLS, 25,                    &
466!                                      slicer_range_limits_dvrp(:,islice_dvrp), &
467!                                      color_dvrp )
468
469             CALL user_dvrp_coltab( 'slicer', output_variable )
470
471             CALL DVRP_COLORTABLE_HLS( m-1, 1, interval_values_dvrp,     &
472                                       interval_h_dvrp, interval_l_dvrp, &
473                                       interval_s_dvrp, interval_a_dvrp )
474
475!
476!--          Compute and plot slicer in dvr-format
477             CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
478                             cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
479!             CALL DVRP_SLICER( m-1, section_mode, slicer_position )
480             CALL DVRP_SLICER( m-1, 2, 1.0 )
481             WRITE (9,*) 'nx_dvrp=', nx_dvrp
482             WRITE (9,*) 'ny_dvrp=', ny_dvrp
483             WRITE (9,*) 'nz_dvrp=', nz_dvrp
484             WRITE (9,*) 'section_mode=', section_mode
485             WRITE (9,*) 'slicer_position=', slicer_position
486             CALL local_flush( 9 )
487
488             CALL DVRP_VISUALIZE( m-1, 2, dvrp_filecount )
489
490             CALL cpu_log( log_point_s(27), 'dvrp_slicer', 'stop' )
491
492!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: ende slicer'
493!   CALL local_flush( 9 )
494          ENDIF
495
496          DEALLOCATE( local_pf )
497
498       ENDIF
499
500       m = m + 1
501
502    ENDDO
503
504    dvrp_filecount = dvrp_filecount + 1
505
506    CALL cpu_log( log_point(27), 'data_output_dvrp', 'stop' )
507!   WRITE ( 9, * ) '*** myid=', myid, ' Ende data_output_dvrp'
508!   CALL local_flush( 9 )
509
510#endif
511 END SUBROUTINE data_output_dvrp
Note: See TracBrowser for help on using the repository browser.