source: palm/trunk/SOURCE/data_output_dvrp.f90 @ 86

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

small bugfixes

  • Property svn:keywords set to Id
File size: 18.0 KB
RevLine 
[1]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! -----------------
[86]34! TEST: different colours for isosurfaces
[1]35! TEST: write statements
36!
37! Former revisions:
38! -----------------
[3]39! $Id: data_output_dvrp.f90 86 2007-05-16 16:00:18Z raasch $
[77]40!
[83]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!
[77]45! 75 2007-03-22 09:54:05Z raasch
46! Particles-package is now part of the default code,
47! moisture renamed humidity
48!
[3]49! RCS Log replace by Id keyword, revision history cleaned up
50!
[1]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'
[82]91    CALL local_flush( 9 )
[1]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.
[86]106!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: vor steering_update'
107!   CALL local_flush( 9 )
[1]108          CALL DVRP_STEERING_UPDATE( m-1, data_output_dvrp )
[86]109!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: nach steering_update'
110!   CALL local_flush( 9 )
[1]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)
[86]142       IF ( mode_dvrp(m)(1:9) == 'particles'  .AND.  particle_advection  .AND. &
143            simulated_time >= particle_advection_start )  THEN
[1]144
[86]145!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: anfang particles'
146!   CALL local_flush( 9 )
[1]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
155!
156!--       Move particle coordinates to one-dimensional arrays
157          IF ( .NOT. use_particle_tails )  THEN
158!
159!--          All particles are output
160             ALLOCATE( psize(number_of_particles), p_t(number_of_particles), &
161                       p_c(number_of_particles), p_x(number_of_particles),   &
162                       p_y(number_of_particles), p_z(number_of_particles) )
163             psize = 0.0;  p_t = 0;  p_c = 0.0;  p_x = 0.0;  p_y = 0.0
164             p_z   = 0.0;
165             psize = particles(1:number_of_particles)%dvrp_psize
166             p_x   = particles(1:number_of_particles)%x * superelevation_x
167             p_y   = particles(1:number_of_particles)%y * superelevation_y
168             p_z   = particles(1:number_of_particles)%z * superelevation
169             p_c   = particles(1:number_of_particles)%color
170          ELSE
171!
172!--          Particles have a tail
[86]173!            WRITE (9,*) '--- before ALLOCATE  simtime=',simulated_time,' #of_tails=', number_of_tails, &
174!                          ' max#of_tp=', maximum_number_of_tailpoints
175!   CALL local_flush( 9 )
[1]176             ALLOCATE( psize(number_of_tails), p_t(number_of_tails),      &
177                       p_c(number_of_tails*maximum_number_of_tailpoints), &
178                       p_x(number_of_tails*maximum_number_of_tailpoints), &
179                       p_y(number_of_tails*maximum_number_of_tailpoints), &
180                       p_z(number_of_tails*maximum_number_of_tailpoints) )
[86]181!            WRITE (9,*) '--- after ALLOCATE'
182!   CALL local_flush( 9 )
[1]183             psize = 0.0;  p_t = 0;  p_c = 0.0;  p_x = 0.0;  p_y = 0.0
184             p_z   = 0.0;
185             i = 0
186             k = 0
187             DO  n = 1, number_of_particles
188                nn = particles(n)%tail_id
189                IF ( nn /= 0 )  THEN
190                   k = k + 1
[86]191!                  IF ( simulated_time > 1338.0 )  THEN
192!                     WRITE (9,*) '--- particle ',n,' tail_id=',nn,' #of_tp=',particles(n)%tailpoints
193!   CALL local_flush( 9 )
194!                  ENDIF
[1]195                   DO  j = 1, particles(n)%tailpoints
196                      i = i + 1
197                      p_x(i) = particle_tail_coordinates(j,1,nn) * &
198                                                                superelevation_x
199                      p_y(i) = particle_tail_coordinates(j,2,nn) * &
200                                                                superelevation_y
201                      p_z(i) = particle_tail_coordinates(j,3,nn) * &
202                                                                superelevation
203                      p_c(i) = particle_tail_coordinates(j,4,nn)
[86]204!                     IF ( simulated_time > 1338.0 )  THEN
205!                        WRITE (9,*) '--- tp= ',i,' x=',p_x(i),' y=',p_y(i), &
206!                                                ' z=',p_z(i),' c=',p_c(i)
207!   CALL local_flush( 9 )
208!                     ENDIF
[1]209                   ENDDO
210                   psize(k) = particles(n)%dvrp_psize
211                   p_t(k)   = particles(n)%tailpoints - 1
[86]212!                  IF ( simulated_time > 1338.0 )  THEN
213!                     WRITE (9,*) '--- t= ',k,' psize=',psize(k),' p_t=',p_t(k)
214!   CALL local_flush( 9 )
215!                  ENDIF
[1]216                ENDIF               
217             ENDDO
[86]218!            WRITE (9,*) '--- after locally storing the particle attributes'
219!   CALL local_flush( 9 )
[1]220          ENDIF
221
222!
223!--       Compute and plot particles in dvr-format
224          IF ( uniform_particles  .AND.  .NOT. use_particle_tails )  THEN
225!
226!--          All particles have the same color. Use simple routine to set
227!--          the particle attributes (produces less output data)
228             CALL DVRP_PARTICLES( m-1, p_x, p_y, p_z, psize )
229          ELSE
230!
231!--          Set color definitions
232             CALL user_dvrp_coltab( 'particles', 'none' )
233
234             CALL DVRP_COLORTABLE_HLS( m-1, 0, interval_values_dvrp,     &
235                                       interval_h_dvrp, interval_l_dvrp, &
236                                       interval_s_dvrp, interval_a_dvrp )
237
238             IF ( .NOT. use_particle_tails )  THEN
239                CALL DVRP_PARTICLES( m-1, number_of_particles, p_x, p_y, p_z, &
240                                     3, psize, p_c, p_t )
241             ELSE
[86]242!               WRITE (9,*) '--- before DVRP_PARTICLES'
243!   CALL local_flush( 9 )
[1]244                CALL DVRP_PARTICLES( m-1, number_of_tails, p_x, p_y, p_z, 15, &
245                                     psize, p_c, p_t )
[86]246!               WRITE (9,*) '--- after DVRP_PARTICLES'
247!               WRITE (9,*) 'm-1 = ',m-1
248!               WRITE (9,*) 'number_of_tails=', number_of_tails
249!               WRITE (9,*) 'p_x =', p_x
250!               WRITE (9,*) 'p_y =', p_y
251!               WRITE (9,*) 'p_z =', p_z
252!               WRITE (9,*) 'psize =', psize
253!               WRITE (9,*) 'p_c =', p_c
254!               WRITE (9,*) 'p_t =', p_t
[60]255
[86]256!   CALL local_flush( 9 )
[1]257             ENDIF
258          ENDIF
259
260          CALL DVRP_VISUALIZE( m-1, 3, dvrp_filecount )
[86]261!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: ende particles'
262!   CALL local_flush( 9 )
[1]263
264          DEALLOCATE( psize, p_c, p_t, p_x, p_y, p_z )
265
266          CALL cpu_log( log_point_s(28), 'dvrp_particles', 'stop' )
267
268
269       ELSEIF ( ( mode_dvrp(m)(1:10) == 'isosurface'  .OR.   &
270                  mode_dvrp(m)(1:6)  == 'slicer'           ) &
271                  .AND.  output_variable /= ' ' )  THEN
272
273!
274!--       Create an intermediate array, properly dimensioned for plot-output
275          ALLOCATE( local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d) )
276
277!
278!--       Move original array to intermediate array
279          SELECT CASE ( output_variable )
280
281             CASE ( 'u', 'u_xy', 'u_xz', 'u_yz' )
282                DO  i = nxl, nxr+1
283                   DO  j = nys, nyn+1
284                      DO  k = nzb, nz_do3d
285                         local_pf(i,j,k) = u(k,j,i)
286                      ENDDO
287                   ENDDO
288                ENDDO
289!
290!--             Replace mirrored values at lower surface by real surface values
291                IF ( output_variable == 'u_xz'  .OR. &
292                     output_variable == 'u_yz' )  THEN
293                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
294                ENDIF
295
296
297             CASE ( 'v', 'v_xy', 'v_xz', 'v_yz' )
298                DO  i = nxl, nxr+1
299                   DO  j = nys, nyn+1
300                      DO  k = nzb, nz_do3d
301                         local_pf(i,j,k) = v(k,j,i)
302                      ENDDO
303                   ENDDO
304                ENDDO
305!
306!--             Replace mirrored values at lower surface by real surface values
307                IF ( output_variable == 'v_xz'  .OR. &
308                     output_variable == 'v_yz' )  THEN
309                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
310                ENDIF
311
312             CASE ( 'w', 'w_xy', 'w_xz', 'w_yz' )
313                DO  i = nxl, nxr+1
314                   DO  j = nys, nyn+1
315                      DO  k = nzb, nz_do3d
316                         local_pf(i,j,k) = w(k,j,i)
317                      ENDDO
318                   ENDDO
319                ENDDO
[86]320                DO  k = nzb, nz_do3d
321                   DO  j = nys+1, nyn
322                      DO  i = nxl, nxr+1
323                         local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + &
324                                           0.50 * local_pf(i,j,k)   + &
325                                           0.25 * local_pf(i,j+1,k)
326                      ENDDO
327                   ENDDO
328                ENDDO
[1]329
330             CASE ( 'p', 'p_xy', 'p_xz', 'p_yz' )
331                DO  i = nxl, nxr+1
332                   DO  j = nys, nyn+1
333                      DO  k = nzb, nz_do3d
334                         local_pf(i,j,k) = p(k,j,i)
335                      ENDDO
336                   ENDDO
337                ENDDO
338
339             CASE ( 'pt', 'pt_xy', 'pt_xz', 'pt_yz' )
340                IF ( .NOT. cloud_physics ) THEN
341                   DO  i = nxl, nxr+1
342                      DO  j = nys, nyn+1
343                         DO  k = nzb, nz_do3d
344                            local_pf(i,j,k) = pt(k,j,i)
345                         ENDDO
346                      ENDDO
347                   ENDDO
348                ELSE
349                   DO  i = nxl, nxr+1
350                      DO  j = nys, nyn+1
351                         DO  k = nzb, nz_do3d
352                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp * pt_d_t(k) * &
353                                                          ql(k,j,i)
354                         ENDDO
355                      ENDDO
356                   ENDDO
357                ENDIF
358
359             CASE ( 'q', 'q_xy', 'q_xz', 'q_yz' )
[75]360                IF ( humidity  .OR.  passive_scalar )  THEN
[1]361                   DO  i = nxl, nxr+1
362                      DO  j = nys, nyn+1
363                         DO  k = nzb, nz_do3d
364                            local_pf(i,j,k) = q(k,j,i)
365                         ENDDO
366                      ENDDO
367                   ENDDO           
368                ELSE
369                   IF ( myid == 0 )  THEN
[75]370                      PRINT*, '+++ data_output_dvrp: if humidity/passive_scalar = ', & 
[1]371                              'FALSE output of ', output_variable,            &
372                              'is not provided' 
373                   ENDIF
374                ENDIF
375             
376             CASE ( 'ql', 'ql_xy', 'ql_xz', 'ql_yz' )
377                IF ( cloud_physics  .OR.  cloud_droplets )  THEN
378                   DO  i = nxl, nxr+1
379                      DO  j = nys, nyn+1
380                         DO  k = nzb, nz_do3d
381                            local_pf(i,j,k) = ql(k,j,i)
382                         ENDDO
383                      ENDDO
384                   ENDDO
385                ELSE
386                   IF ( myid == 0 ) THEN
387                      PRINT*, '+++ data_output_dvrp: if cloud_physics = FALSE ', & 
388                              'output of ', output_variable, 'is not provided' 
389                   ENDIF
390                ENDIF
391
392             CASE ( 'u*_xy' )
393                DO  i = nxl, nxr+1
394                   DO  j = nys, nyn+1
395                      local_pf(i,j,nzb+1) = us(j,i)
396                   ENDDO
397                ENDDO
398                slicer_position = zu(nzb+1)
399
400             CASE ( 't*_xy' )
401                DO  i = nxl, nxr+1
402                   DO  j = nys, nyn+1
403                      local_pf(i,j,nzb+1) = ts(j,i)
404                   ENDDO
405                ENDDO
406                slicer_position = zu(nzb+1)
407
408
409             CASE DEFAULT
410                IF ( myid == 0 )  THEN
411                   PRINT*,'+++ data_output_dvrp: no output possible for: ', &
412                          output_variable
413                ENDIF
414
415          END SELECT
416
417
418          IF ( mode_dvrp(m)(1:10) == 'isosurface' )  THEN
419
[86]420!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: anfang isosurface'
421!   CALL local_flush( 9 )
[1]422!
423!--          DVRP-Calls for plotting isosurfaces:
424             CALL cpu_log( log_point_s(26), 'dvrp_isosurface', 'start' )
425
426!
427!--          Definition of characteristics of isosurface material
428!--          Preliminary settings for w!
429             IF ( output_variable == 'w' )  THEN
[86]430                IF ( tv == 1 )  THEN
431                   CALL DVRP_MATERIAL_RGB( m-1, 1, 0.8, 0.1, 0.1, 0.0 )
432                ELSE
433                   CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.1, 0.8, 0.0 )
434                ENDIF
[1]435             ELSE
436                CALL DVRP_MATERIAL_RGB( m-1, 1, 0.9, 0.9, 0.9, 0.0 )
437             ENDIF
438
439!
440!--          Compute and plot isosurface in dvr-format
441             CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
442                             cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
443             CALL DVRP_THRESHOLD( m-1, threshold(tv) )
444             CALL DVRP_VISUALIZE( m-1, 1, dvrp_filecount )
[86]445!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: ende isosurface'
446!   CALL local_flush( 9 )
[1]447
448             CALL cpu_log( log_point_s(26), 'dvrp_isosurface', 'stop' )
449
450          ELSEIF ( mode_dvrp(m)(1:6) == 'slicer' )  THEN
451
[86]452!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: anfang slicer'
453!   CALL local_flush( 9 )
[1]454!
455!--          DVRP-Calls for plotting slicers:
456             CALL cpu_log( log_point_s(27), 'dvrp_slicer', 'start' )
457
458!
459!--          Material and color definitions
460             CALL DVRP_MATERIAL_RGB( m-1, 1, 0.0, 0.0, 0.0, 0.0 )
461
462             islice_dvrp = islice_dvrp + 1
463!             CALL DVRP_COLORFUNCTION( m-1, DVRP_CM_HLS, 25,                    &
464!                                      slicer_range_limits_dvrp(:,islice_dvrp), &
465!                                      color_dvrp )
466
467             CALL user_dvrp_coltab( 'slicer', output_variable )
468
469             CALL DVRP_COLORTABLE_HLS( m-1, 1, interval_values_dvrp,     &
470                                       interval_h_dvrp, interval_l_dvrp, &
471                                       interval_s_dvrp, interval_a_dvrp )
472
473!
474!--          Compute and plot slicer in dvr-format
475             CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
476                             cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
477             CALL DVRP_SLICER( m-1, section_mode, slicer_position )
478             CALL DVRP_VISUALIZE( m-1, 2, dvrp_filecount )
479
480             CALL cpu_log( log_point_s(27), 'dvrp_slicer', 'stop' )
481
[86]482!   WRITE ( 9, * ) '*** myid=', myid, ' data_output_dvrp: ende slicer'
483!   CALL local_flush( 9 )
[1]484          ENDIF
485
486          DEALLOCATE( local_pf )
487
488       ENDIF
489
490       m = m + 1
491
492    ENDDO
493
494    dvrp_filecount = dvrp_filecount + 1
495
496    CALL cpu_log( log_point(27), 'data_output_dvrp', 'stop' )
[86]497!   WRITE ( 9, * ) '*** myid=', myid, ' Ende data_output_dvrp'
498!   CALL local_flush( 9 )
[1]499
500#endif
501 END SUBROUTINE data_output_dvrp
Note: See TracBrowser for help on using the repository browser.