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

Last change on this file since 66 was 60, checked in by raasch, 18 years ago

preliminary update of further changes, running

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