Changeset 828
- Timestamp:
- Feb 21, 2012 12:00:36 PM (13 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r826 r828 4 4 # Current revisions: 5 5 # ------------------ 6 # 6 # init_particles depends on lpm_collision_kernels 7 7 # 8 8 # Former revisions: … … 66 66 67 67 RCS = advec_particles.f90 advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \ 68 69 70 71 72 73 74 75 76 77 data_output_ptseries.f90 data_output_spectra.f90 data_output_tseries.f90 \78 data_output_2d.f90 data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 \79 diffusion_u.f90 diffusion_v.f90 diffusion_w.f90 diffusivities.f90 \80 disturb_field.f90 disturb_heatflux.f90 eqn_state_seawater.f90 \81 exchange_horiz.f90 exchange_horiz_2d.f90 \82 fft_xy.f90 flow_statistics.f90 global_min_max.f90 \83 header.f90 impact_of_latent_heat.f90 inflow_turbulence.f90 \84 in it_1d_model.f90 init_3d_model.f90 init_advec.f90 \85 init_ cloud_physics.f90 init_coupling.f90 init_dvrp.f90 init_grid.f90 \86 init_ masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \87 init_p t_anomaly.f90 init_rankine.f90 init_slope.f90 \68 advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \ 69 advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \ 70 advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \ 71 average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \ 72 calc_liquid_water_content.f90 calc_precipitation.f90 \ 73 calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \ 74 check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \ 75 coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \ 76 data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \ 77 data_output_ptseries.f90 data_output_spectra.f90 \ 78 data_output_tseries.f90 data_output_2d.f90 data_output_3d.f90 \ 79 diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 diffusion_v.f90 \ 80 diffusion_w.f90 diffusivities.f90 disturb_field.f90 \ 81 disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \ 82 exchange_horiz_2d.f90 fft_xy.f90 flow_statistics.f90 \ 83 global_min_max.f90 header.f90 impact_of_latent_heat.f90 \ 84 inflow_turbulence.f90 init_1d_model.f90 init_3d_model.f90 \ 85 init_advec.f90 init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 \ 86 init_grid.f90 init_masks.f90 init_ocean.f90 init_particles.f90 \ 87 init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \ 88 88 interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \ 89 local_stop.f90 local_system.f90 local_tremain.f90 local_tremain_ini.f90\90 lpm_collision_kernels.f90 message.f90 modules.f90 netcdf.f90 \91 package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \92 89 local_stop.f90 local_system.f90 local_tremain.f90 \ 90 local_tremain_ini.f90 lpm_collision_kernels.f90 message.f90 \ 91 modules.f90 netcdf.f90 package_parin.f90 palm.f90 parin.f90 \ 92 particle_boundary_conds.f90 plant_canopy_model.f90 poisfft.f90 \ 93 93 poisfft_hybrid.f90 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \ 94 94 production_e.f90 prognostic_equations.f90 random_function.f90 \ 95 95 random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \ 96 96 set_particle_attributes.f90 set_slicer_attributes_dvrp.f90 \ 97 97 singleton.f90 sor.f90 spline_x.f90 \ 98 98 spline_y.f90 spline_z.f90 subsidence.f90 \ 99 99 sum_up_3d_data.f90 surface_coupler.f90 \ 100 100 swap_timelevel.f90 temperton_fft.f90 time_integration.f90 \ 101 102 101 time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \ 102 transpose.f90 user_3d_data_averaging.f90 user_actions.f90 \ 103 103 user_additional_routines.f90 user_advec_particles.f90 \ 104 104 user_check_data_output.f90 user_check_data_output_pr.f90 \ 105 105 user_check_parameters.f90 user_data_output_2d.f90 \ 106 user_data_output_3d.f90 user_data_output_dvrp.f90 user_data_output_mask.f90\107 user_d efine_netcdf_grid.f90 user_dvrp_coltab.f90 user_header.f90 \108 user_ init.f90 user_init_3d_model.f90 user_init_grid.f90 \109 user_init_ particles.f90 user_init_plant_canopy.f90 \110 user_ last_actions.f90 user_module.f90 user_parin.f90 \111 user_par ticle_attributes.f90 user_read_restart_data.f90 \106 user_data_output_3d.f90 user_data_output_dvrp.f90 \ 107 user_data_output_mask.f90 user_define_netcdf_grid.f90 \ 108 user_dvrp_coltab.f90 user_header.f90 user_init.f90 \ 109 user_init_3d_model.f90 user_init_grid.f90 user_init_particles.f90 \ 110 user_init_plant_canopy.f90 user_last_actions.f90 user_module.f90 \ 111 user_parin.f90 user_particle_attributes.f90 user_read_restart_data.f90 \ 112 112 user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \ 113 write_3d_binary.f90 write_compressed.f90 \ 114 write_var_list.f90 113 write_3d_binary.f90 write_compressed.f90 write_var_list.f90 115 114 116 115 OBJS = advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \ … … 252 251 init_masks.o: modules.o 253 252 init_ocean.o: modules.o eqn_state_seawater.o 254 init_particles.o: modules.o random_function.o253 init_particles.o: modules.o lpm_collision_kernels.o random_function.o 255 254 init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o 256 255 init_pt_anomaly.o: modules.o -
palm/trunk/SOURCE/advec_particles.f90
r826 r828 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 6 ! fast hall/wang kernels with fixed radius/dissipation classes added, 7 ! particle feature color renamed class, routine colker renamed 8 ! recalculate_kernel, 9 ! lower limit for droplet radius changed from 1E-7 to 1E-8 10 ! 11 ! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4 7 12 ! 8 13 ! Former revisions: … … 144 149 IMPLICIT NONE 145 150 146 INTEGER :: agp, deleted_particles, deleted_tails, i, ie, ii, inc,&147 internal_timestep_count, is, j, &148 jj, js, jtry, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi,&151 INTEGER :: agp, deleted_particles, deleted_tails, eclass, i, ie, ii, inc, & 152 internal_timestep_count, is, j, jj, js, jtry, k, kk, kw, m, n, & 153 nc, nd, nn, num_gp, pse, psi, rclass_l, rclass_s, & 149 154 tlength, trlp_count, trlp_count_sum, trlp_count_recv, & 150 155 trlp_count_recv_sum, trlpt_count, trlpt_count_recv, & … … 173 178 lagr_timescale, lw_max, mean_r, new_r, p_int, pt_int, & 174 179 pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, & 175 ql_int_u, random_gauss, r_ros, r_ros_ini, sigma, sl_r3, sl_r4,&176 t_int, u_int, u_int_l, u_int_u,vv_int, v_int, v_int_l, &177 v_int _u, w_int, w_int_l, w_int_u, x, y180 ql_int_u, random_gauss, r_ros, r_ros_ini, & 181 sigma, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,vv_int, & 182 v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u, x, y 178 183 ! 179 184 !-- Parameters for Rosenbrock method … … 190 195 191 196 REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei 192 REAL, DIMENSION(:,:), ALLOCATABLE :: kern193 197 194 198 REAL :: location(1:30,1:3) … … 356 360 thermal_conductivity_l + & 357 361 rho_l * r_v * t_int / diff_coeff_l / e_s ) 358 IF ( arg < 1.0E-1 4) THEN359 new_r = 1.0E- 7362 IF ( arg < 1.0E-16 ) THEN 363 new_r = 1.0E-8 360 364 ELSE 361 365 new_r = SQRT( arg ) … … 550 554 CALL message( 'advec_particles', 'PA0144', 2, 2, -1, 6, 1 ) 551 555 ENDIF 552 particles(n)%radius = new_r553 556 554 557 ! 555 558 !-- Sum up the total volume of liquid water (needed below for 556 559 !-- re-calculating the weighting factors) 557 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * & 558 particles(n)%radius**3 560 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r**3 561 562 particles(n)%radius = new_r 563 564 ! 565 !-- Determine radius class of the particle needed for collision 566 IF ( ( hall_kernel .OR. wang_kernel ) .AND. use_kernel_tables ) & 567 THEN 568 particles(n)%class = ( LOG( new_r ) - rclass_lbound ) / & 569 ( rclass_ubound - rclass_lbound ) * & 570 radius_classes 571 particles(n)%class = MIN( particles(n)%class, radius_classes ) 572 particles(n)%class = MAX( particles(n)%class, 1 ) 573 ENDIF 574 559 575 ENDDO 560 576 CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' ) … … 601 617 pse = psi + prt_count(k,j,i)-1 602 618 603 IF ( hall_kernel .OR. wang_kernel ) THEN 604 605 ALLOCATE( kern(psi:pse,psi:pse) ) 606 607 ! 608 !-- Calculate collision kernel for all particles in grid box 609 CALL colker( i, j, k, kern ) 610 ! 611 !-- Collison kernel is calculated in cm**3/s but needed 612 !-- in m**3/s 613 kern = kern * 1.0E-6 ! to be moved to colker 619 ! 620 !-- Now apply the different kernels 621 IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 622 use_kernel_tables ) THEN 623 ! 624 !-- Fast method with pre-calculated efficiencies for 625 !-- discrete radius- and dissipation-classes. 626 ! 627 !-- Determine dissipation class index of this gridbox 628 eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * & 629 dissipation_classes ) + 1 630 IF ( hall_kernel .OR. diss(k,j,i) * 1.0E4 < 0.001 ) THEN 631 eclass = 0 ! Hall kernel is used 632 ELSE 633 eclass = MIN( dissipation_classes, eclass ) 634 ENDIF 614 635 615 636 DO n = pse, psi+1, -1 … … 617 638 integral = 0.0 618 639 lw_max = 0.0 619 640 rclass_l = particles(n)%class 620 641 ! 621 642 !-- Calculate growth of collector particle radius using all … … 623 644 DO is = psi, n-1 624 645 625 integral = integral + & 626 ( particles(is)%radius**3 * kern(n,is) *& 646 rclass_s = particles(is)%class 647 integral = integral + & 648 ( particles(is)%radius**3 * & 649 ckernel(rclass_l,rclass_s,eclass) * & 627 650 particles(is)%weight_factor ) 628 651 ! … … 653 676 654 677 ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN 655 678 !-- can this case really happen?? 656 679 DO is = psi, n-1 657 680 … … 667 690 DO is = psi, n-1 668 691 692 rclass_s = particles(is)%class 669 693 particles(is)%weight_factor = & 670 694 particles(is)%weight_factor - & 671 ( ( kern(n,is) * particles(is)%weight_factor &695 ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor & 672 696 / integral ) * delta_v ) 673 697 … … 699 723 ENDDO 700 724 701 DEALLOCATE( kern ) 725 ELSEIF ( ( hall_kernel .OR. wang_kernel ) .AND. & 726 .NOT. use_kernel_tables ) THEN 727 ! 728 !-- Collision efficiencies are calculated for every new 729 !-- grid box. First, allocate memory for kernel table. 730 !-- Third dimension is 1, because table is re-calculated for 731 !-- every new dissipation value. 732 ALLOCATE( ckernel(prt_start_index(k,j,i): & 733 prt_start_index(k,j,i)+prt_count(k,j,i)-1, & 734 prt_start_index(k,j,i): & 735 prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) ) 736 ! 737 !-- Now calculate collision efficiencies for this box 738 CALL recalculate_kernel( i, j, k ) 739 740 DO n = pse, psi+1, -1 741 742 integral = 0.0 743 lw_max = 0.0 744 ! 745 !-- Calculate growth of collector particle radius using all 746 !-- droplets smaller than current droplet 747 DO is = psi, n-1 748 749 integral = integral + particles(is)%radius**3 * & 750 ckernel(n,is,1) * & 751 particles(is)%weight_factor 752 ! 753 !-- Calculate volume of liquid water of the collected 754 !-- droplets which is the maximum liquid water available 755 !-- for droplet growth 756 lw_max = lw_max + ( particles(is)%radius**3 * & 757 particles(is)%weight_factor ) 758 759 ENDDO 760 761 ! 762 !-- Change in radius of collector droplet due to collision 763 delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & 764 integral * dt_3d * ddx * ddy / dz 765 766 ! 767 !-- Change in volume of collector droplet due to collision 768 delta_v = particles(n)%weight_factor & 769 * ( ( particles(n)%radius + delta_r )**3 & 770 - particles(n)%radius**3 ) 771 772 IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN 773 !-- replace by message call 774 PRINT*, 'Particle has grown to much because the', & 775 ' change of volume of particle is larger', & 776 ' than liquid water available!' 777 778 ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN 779 !-- can this case really happen?? 780 DO is = psi, n-1 781 782 particles(is)%weight_factor = 0.0 783 particle_mask(is) = .FALSE. 784 deleted_particles = deleted_particles + 1 785 786 ENDDO 787 788 ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN 789 ! 790 !-- Calculate new weighting factor of collected droplets 791 DO is = psi, n-1 792 793 particles(is)%weight_factor = & 794 particles(is)%weight_factor - & 795 ( ckernel(n,is,1) / integral * delta_v * & 796 particles(is)%weight_factor ) 797 798 IF ( particles(is)%weight_factor < 0.0 ) THEN 799 WRITE( message_string, * ) 'negative ', & 800 'weighting factor: ', & 801 particles(is)%weight_factor 802 CALL message( 'advec_particles', '', 2, 2,& 803 -1, 6, 1 ) 804 805 ELSEIF ( particles(is)%weight_factor == 0.0 ) & 806 THEN 807 808 particles(is)%weight_factor = 0.0 809 particle_mask(is) = .FALSE. 810 deleted_particles = deleted_particles + 1 811 812 ENDIF 813 814 ENDDO 815 816 ENDIF 817 818 particles(n)%radius = particles(n)%radius + delta_r 819 ql_vp(k,j,i) = ql_vp(k,j,i) + & 820 particles(n)%weight_factor & 821 * particles(n)%radius**3 822 823 ENDDO 824 825 DEALLOCATE( ckernel ) 702 826 703 827 ELSEIF ( palm_kernel ) THEN … … 2376 2500 particle_tail_coordinates(1,2,nn) = particles(n)%y 2377 2501 particle_tail_coordinates(1,3,nn) = particles(n)%z 2378 particle_tail_coordinates(1,4,nn) = particles(n)%c olor2502 particle_tail_coordinates(1,4,nn) = particles(n)%class 2379 2503 particles(n)%tailpoints = 1 2380 2504 IF ( minimum_tailpoint_distance /= 0.0 ) THEN … … 2382 2506 particle_tail_coordinates(2,2,nn) = particles(n)%y 2383 2507 particle_tail_coordinates(2,3,nn) = particles(n)%z 2384 particle_tail_coordinates(2,4,nn) = particles(n)%c olor2508 particle_tail_coordinates(2,4,nn) = particles(n)%class 2385 2509 particle_tail_coordinates(1:2,5,nn) = 0.0 2386 2510 particles(n)%tailpoints = 2 … … 3743 3867 3744 3868 ! 3745 !-- Set particle attributes 3746 CALL set_particle_attributes 3869 !-- Set particle attributes. 3870 !-- Feature is not available if collision is activated, because the respective 3871 !-- particle attribute (class) is then used for storing the particle radius 3872 !-- class. 3873 IF ( collision_kernel == 'none' ) CALL set_particle_attributes 3747 3874 3748 3875 ! 3749 3876 !-- Set particle attributes defined by the user 3750 3877 CALL user_particle_attributes 3751 ! WRITE ( 9, * ) '*** advec_particles: ##10'3752 ! CALL local_flush( 9 )3753 ! DO n = 1, number_of_particles3754 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &3755 ! THEN3756 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'3757 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'3758 ! CALL local_flush( 9 )3759 ! CALL MPI_ABORT( comm2d, 9999, ierr )3760 ! ENDIF3761 ! ENDDO3762 3878 3763 3879 ! … … 3810 3926 particle_tail_coordinates(1,2,nn) = particles(n)%y 3811 3927 particle_tail_coordinates(1,3,nn) = particles(n)%z 3812 particle_tail_coordinates(1,4,nn) = particles(n)%c olor3928 particle_tail_coordinates(1,4,nn) = particles(n)%class 3813 3929 ! WRITE ( 9, * ) '*** advec_particles: ##10.4' 3814 3930 ! CALL local_flush( 9 ) … … 4058 4174 CALL handle_netcdf_error( 'output_particles_netcdf', 15 ) 4059 4175 4060 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%c olor, &4176 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class, & 4061 4177 start = (/ 1, prt_time_count /), & 4062 4178 count = (/ maximum_number_of_particles /) ) -
palm/trunk/SOURCE/check_parameters.f90
r826 r828 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! check of collision_kernel extended 7 7 ! 8 8 ! Former revisions: … … 716 716 SELECT CASE ( TRIM( collision_kernel ) ) 717 717 718 CASE ( 'hall' )718 CASE ( 'hall', 'hall_fast' ) 719 719 hall_kernel = .TRUE. 720 720 … … 722 722 palm_kernel = .TRUE. 723 723 724 CASE ( 'wang' )724 CASE ( 'wang', 'wang_fast' ) 725 725 wang_kernel = .TRUE. 726 726 … … 734 734 735 735 END SELECT 736 IF ( collision_kernel(6:9) == 'fast' ) use_kernel_tables = .TRUE. 736 737 737 738 -
palm/trunk/SOURCE/data_output_dvrp.f90
r392 r828 32 32 ! Current revisions: 33 33 ! ----------------- 34 ! 34 ! particle feature color renamed class 35 35 ! 36 36 ! Former revisions: … … 226 226 p_y(k) = particles(n)%y * superelevation_y 227 227 p_z(k) = particles(n)%z * superelevation 228 p_c(k) = particles(n)%c olor228 p_c(k) = particles(n)%class 229 229 ENDIF 230 230 ENDDO -
palm/trunk/SOURCE/header.f90
r826 r828 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! output of dissipation_classes + radius_classes 7 7 ! 8 8 ! Former revisions: … … 1483 1483 IF ( collision_kernel /= 'none' ) THEN 1484 1484 WRITE ( io, 435 ) TRIM( collision_kernel ) 1485 ELSE 1486 WRITE ( io, 436 ) 1485 IF ( collision_kernel(6:9) == 'fast' ) THEN 1486 WRITE ( io, 436 ) radius_classes, dissipation_classes 1487 ENDIF 1488 ELSE 1489 WRITE ( io, 437 ) 1487 1490 ENDIF 1488 1491 ENDIF … … 1964 1967 ' droplets < 1.0E-6 m') 1965 1968 435 FORMAT (' Droplet collision is handled by ',A,'-kernel') 1966 436 FORMAT (' Droplet collision is switched off') 1969 436 FORMAT (' Fast kernel with fixed radius- and dissipation classes ', & 1970 'are used'/ & 1971 ' number of radius classes: ',I3,' interval ', & 1972 '[1.0E-6,2.0E-4] m'/ & 1973 ' number of dissipation classes: ',I2,' interval ', & 1974 '[0,1000] cm**2/s**3') 1975 437 FORMAT (' Droplet collision is switched off') 1967 1976 450 FORMAT (//' LES / Turbulence quantities:'/ & 1968 1977 ' ---------------------------'/) -
palm/trunk/SOURCE/init_particles.f90
r826 r828 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! call of init_kernels, particle feature color renamed class 7 7 ! 8 8 ! Former revisions: … … 72 72 USE grid_variables 73 73 USE indices 74 USE lpm_collision_kernels_mod 74 75 USE particle_attributes 75 76 USE pegrid … … 148 149 IF ( pdz(j) == 9999999.9 .OR. pdz(j) == 0.0 ) pdz(j) = pdz(j-1) 149 150 ENDDO 151 152 ! 153 !-- Initialize collision kernels 154 IF ( collision_kernel /= 'none' ) CALL init_kernels 150 155 151 156 ! … … 343 348 particles(n)%radius = particle_groups(i)%radius 344 349 particles(n)%weight_factor =initial_weighting_factor 345 particles(n)%c olor= 1350 particles(n)%class = 1 346 351 particles(n)%group = i 347 352 particles(n)%tailpoints = 0 … … 470 475 IF ( MINVAL( particles(1:n)%dvrp_psize ) == & 471 476 MAXVAL( particles(1:n)%dvrp_psize ) .AND. & 472 MINVAL( particles(1:n)%c olor) == &473 MAXVAL( particles(1:n)%c olor) ) THEN477 MINVAL( particles(1:n)%class ) == & 478 MAXVAL( particles(1:n)%class ) ) THEN 474 479 uniform_particles_l = .TRUE. 475 480 ELSE … … 532 537 particle_tail_coordinates(1,2,nn) = particles(n)%y 533 538 particle_tail_coordinates(1,3,nn) = particles(n)%z 534 particle_tail_coordinates(1,4,nn) = particles(n)%c olor539 particle_tail_coordinates(1,4,nn) = particles(n)%class 535 540 particles(n)%tailpoints = 1 536 541 IF ( minimum_tailpoint_distance /= 0.0 ) THEN … … 538 543 particle_tail_coordinates(2,2,nn) = particles(n)%y 539 544 particle_tail_coordinates(2,3,nn) = particles(n)%z 540 particle_tail_coordinates(2,4,nn) = particles(n)%c olor545 particle_tail_coordinates(2,4,nn) = particles(n)%class 541 546 particle_tail_coordinates(1:2,5,nn) = 0.0 542 547 particles(n)%tailpoints = 2 -
palm/trunk/SOURCE/lpm_collision_kernels.f90
r826 r828 1 MODULE lpm_collision_kernels_mod1 MODULE lpm_collision_kernels_mod 2 2 3 3 !------------------------------------------------------------------------------! 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! code has been completely reformatted, routine colker renamed 7 ! recalculate_kernel, 8 ! routine init_kernels added, radius is now communicated to the collision 9 ! routines by array radclass 10 ! 11 ! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4 7 12 ! 8 13 ! Former revisions: … … 24 29 ! Description: 25 30 ! ------------ 26 ! This routine calculates the effect of (SGS) turbulence on the collision 27 ! efficiency of droplets. 28 ! It is based on the original kernel developed by Wang (...) 31 ! This module calculates collision efficiencies either due to pure gravitational 32 ! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or 33 ! including the effects of (SGS) turbulence (Wang kernel, see Wang and 34 ! Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been 35 ! provided by L.-P. Wang but is substantially reformatted and speed optimized 36 ! here. 37 ! 38 ! ATTENTION: 39 ! Physical quantities (like g, densities, etc.) used in this module still 40 ! have to be adjusted to those values used in the main PALM code. 41 ! Also, quantities in CGS-units should be converted to SI-units eventually. 29 42 !------------------------------------------------------------------------------! 30 43 … … 33 46 USE constants 34 47 USE particle_attributes 48 USE pegrid 49 35 50 36 51 IMPLICIT NONE … … 38 53 PRIVATE 39 54 40 PUBLIC colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 41 42 INTEGER, SAVE :: ip, jp, kp, pend, pstart, psum 43 44 REAL, SAVE :: epsilon, eps2, urms, urms2 45 46 REAL, DIMENSION(:), ALLOCATABLE, SAVE :: winf 47 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ec, ecf, gck 55 PUBLIC ckernel, init_kernels, rclass_lbound, rclass_ubound, & 56 recalculate_kernel 57 58 REAL :: epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2 59 60 REAL, DIMENSION(:), ALLOCATABLE :: epsclass, radclass, winf 61 REAL, DIMENSION(:,:), ALLOCATABLE :: ec, ecf, gck, hkernel, hwratio 62 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ckernel 63 64 SAVE 48 65 49 66 ! 50 67 !-- Public interfaces 51 INTERFACE colker 52 MODULE PROCEDURE colker 53 END INTERFACE colker 54 55 INTERFACE effic 56 MODULE PROCEDURE effic 57 END INTERFACE effic 58 59 INTERFACE fallg 60 MODULE PROCEDURE fallg 61 END INTERFACE fallg 62 63 INTERFACE phi 64 MODULE PROCEDURE phi 65 END INTERFACE phi 66 67 INTERFACE turbsd 68 MODULE PROCEDURE turbsd 69 END INTERFACE turbsd 70 71 INTERFACE turb_enhance_eff 72 MODULE PROCEDURE turb_enhance_eff 73 END INTERFACE turb_enhance_eff 74 75 INTERFACE zhi 76 MODULE PROCEDURE zhi 77 END INTERFACE zhi 68 INTERFACE init_kernels 69 MODULE PROCEDURE init_kernels 70 END INTERFACE init_kernels 71 72 INTERFACE recalculate_kernel 73 MODULE PROCEDURE recalculate_kernel 74 END INTERFACE recalculate_kernel 75 78 76 79 77 CONTAINS 80 78 81 !------------------------------------------------------------------------------! 82 ! SUBROUTINE for calculation of collision kernel 83 !------------------------------------------------------------------------------! 84 SUBROUTINE colker( i1, j1, k1, kernel ) 79 80 SUBROUTINE init_kernels 81 !------------------------------------------------------------------------------! 82 ! Initialization of the collision efficiency matrix with fixed radius and 83 ! dissipation classes, calculated at simulation start only. 84 !------------------------------------------------------------------------------! 85 86 IMPLICIT NONE 87 88 INTEGER :: i, j, k 89 90 91 ! 92 !-- Calculate collision efficiencies for fixed radius- and dissipation 93 !-- classes 94 IF ( collision_kernel(6:9) == 'fast' ) THEN 95 96 ALLOCATE( ckernel(1:radius_classes,1:radius_classes, & 97 0:dissipation_classes), epsclass(1:dissipation_classes), & 98 radclass(1:radius_classes) ) 99 100 ! 101 !-- Calculate the radius class bounds with logarithmic distances 102 !-- in the interval [1.0E-6, 2.0E-4] m 103 rclass_lbound = LOG( 1.0E-6 ) 104 rclass_ubound = LOG( 2.0E-4 ) 105 radclass(1) = 1.0E-6 106 DO i = 2, radius_classes 107 radclass(i) = EXP( rclass_lbound + & 108 ( rclass_ubound - rclass_lbound ) * ( i-1.0 ) /& 109 ( radius_classes - 1.0 ) ) 110 ! IF ( myid == 0 ) THEN 111 ! PRINT*, 'i=', i, ' r = ', radclass(i)*1.0E6 112 ! ENDIF 113 ENDDO 114 ! 115 !-- Collision routines expect radius to be in cm 116 radclass = radclass * 100.0 117 118 ! 119 !-- Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3 120 DO i = 1, dissipation_classes 121 epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes 122 ! IF ( myid == 0 ) THEN 123 ! PRINT*, 'i=', i, ' eps = ', epsclass(i) 124 ! ENDIF 125 ENDDO 126 ! 127 !-- Calculate collision efficiencies of the Wang/ayala kernel 128 ALLOCATE( ec(1:radius_classes,1:radius_classes), & 129 ecf(1:radius_classes,1:radius_classes), & 130 gck(1:radius_classes,1:radius_classes), & 131 winf(1:radius_classes) ) 132 133 DO k = 1, dissipation_classes 134 135 epsilon = epsclass(k) 136 urms = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 ) 137 138 CALL turbsd 139 CALL turb_enhance_eff 140 CALL effic 141 142 DO j = 1, radius_classes 143 DO i = 1, radius_classes 144 ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j) 145 ENDDO 146 ENDDO 147 148 ENDDO 149 150 ! 151 !-- Calculate collision efficiencies of the Hall kernel 152 ALLOCATE( hkernel(1:radius_classes,1:radius_classes), & 153 hwratio(1:radius_classes,1:radius_classes) ) 154 155 CALL fallg 156 CALL effic 157 158 DO j = 1, radius_classes 159 DO i = 1, radius_classes 160 hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 & 161 * ec(i,j) * ABS( winf(j) - winf(i) ) 162 ckernel(i,j,0) = hkernel(i,j) ! hall kernel stored on index 0 163 ENDDO 164 ENDDO 165 166 ! 167 !-- Test output of efficiencies 168 IF ( j == -1 ) THEN 169 170 PRINT*, '*** Hall kernel' 171 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes ) 172 DO j = 1, radius_classes 173 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes ) 174 ENDDO 175 176 DO k = 1, dissipation_classes 177 DO i = 1, radius_classes 178 DO j = 1, radius_classes 179 IF ( hkernel(i,j) == 0.0 ) THEN 180 hwratio(i,j) = 9999999.9 181 ELSE 182 hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j) 183 ENDIF 184 ENDDO 185 ENDDO 186 187 PRINT*, '*** epsilon = ', epsclass(k) 188 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes ) 189 DO j = 1, radius_classes 190 ! WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes ) 191 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes ) 192 ENDDO 193 ENDDO 194 195 ENDIF 196 197 DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf ) 198 199 ckernel = ckernel * 1.0E-6 ! kernel is needed in m**3/s 200 201 ELSEIF( collision_kernel == 'hall' .OR. collision_kernel == 'wang' ) & 202 THEN 203 ! 204 !-- Initial settings for Hall- and Wang-Kernel 205 !-- To be done: move here parts from turbsd, fallg, ecoll, etc. 206 ENDIF 207 208 END SUBROUTINE init_kernels 209 210 211 !------------------------------------------------------------------------------! 212 ! Calculation of collision kernels during each timestep and for each grid box 213 !------------------------------------------------------------------------------! 214 SUBROUTINE recalculate_kernel( i1, j1, k1 ) 85 215 86 216 USE arrays_3d … … 94 224 IMPLICIT NONE 95 225 96 INTEGER :: i, i1, j, j1, k1 97 98 REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel 99 100 ! CALL cpu_log( log_point_s(46), 'colker', 'start' ) 101 102 ip = i1 103 jp = j1 104 kp = k1 105 106 pstart = prt_start_index(kp,jp,ip) 107 pend = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1 108 psum = prt_count(kp,jp,ip) 109 110 ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) ) 111 112 IF ( wang_kernel ) THEN 113 114 ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) ) 115 116 epsilon = diss(kp,jp,ip) * 1.0E5 !dissipation rate in cm**2/s**-3 117 urms = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 ) 118 119 IF ( epsilon <= 0.001 ) THEN 120 121 CALL fallg 122 CALL effic 123 124 DO j = pstart, pend 125 DO i = pstart, pend 126 kernel(i,j) = pi * ( particles(j)%radius * 100.0 + & 127 particles(i)%radius * 100.0 )**2 & 128 * ec(i,j) * ABS( winf(j) - winf(i) ) 129 ENDDO 226 INTEGER :: i, i1, j, j1, k1, pend, pstart 227 228 229 pstart = prt_start_index(k1,j1,i1) 230 pend = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1) - 1 231 radius_classes = prt_count(k1,j1,i1) 232 233 ALLOCATE( ec(1:radius_classes,1:radius_classes), & 234 radclass(1:radius_classes), winf(1:radius_classes) ) 235 236 ! 237 !-- Store particle radii on the radclass array. Collision routines 238 !-- expect radii to be in cm. 239 radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0 240 241 epsilon = diss(k1,j1,i1) * 1.0E4 ! dissipation rate in cm**2/s**-3 242 urms = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 ) 243 244 IF ( wang_kernel .AND. epsilon > 0.001 ) THEN 245 ! 246 !-- Call routines to calculate efficiencies for the Wang kernel 247 ALLOCATE( gck(1:radius_classes,1:radius_classes), & 248 ecf(1:radius_classes,1:radius_classes) ) 249 250 CALL turbsd 251 CALL turb_enhance_eff 252 CALL effic 253 254 DO j = 1, radius_classes 255 DO i = 1, radius_classes 256 ckernel(pstart+i-1,pstart+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j) 130 257 ENDDO 131 132 ELSE 133 134 CALL turbsd 135 CALL turb_enhance_eff 136 CALL effic 137 138 DO j = pstart, pend, 1 139 DO i = pstart, pend, 1 140 kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j) 141 ENDDO 258 ENDDO 259 260 DEALLOCATE( gck, ecf ) 261 262 ELSE 263 ! 264 !-- Call routines to calculate efficiencies for the Hall kernel 265 CALL fallg 266 CALL effic 267 268 DO j = 1, radius_classes 269 DO i = 1, radius_classes 270 ckernel(pstart+i-1,pstart+j-1,1) = pi * & 271 ( radclass(j) + radclass(i) )**2 & 272 * ec(i,j) * ABS( winf(j) - winf(i) ) 142 273 ENDDO 143 144 ENDIF145 146 DEALLOCATE(gck, ecf)147 148 ELSE149 150 ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )151 CALL fallg152 ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )153 ! CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )154 CALL effic155 ! CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )156 157 DO j = pstart, pend158 DO i = pstart, pend159 kernel(i,j) = pi * ( particles(j)%radius * 100.0 + &160 particles(i)%radius * 100.0 )**2 &161 * ec(i,j) * ABS( winf(j) - winf(i) )162 ENDDO163 274 ENDDO 164 275 165 276 ENDIF 166 277 167 DEALLOCATE( ec, winf ) 168 169 ! CALL cpu_log( log_point_s(46), 'colker', 'stop' ) 170 171 END SUBROUTINE colker 172 173 !------------------------------------------------------------------------------! 174 ! SUBROUTINE for calculation of w, g and gck 278 ckernel = ckernel * 1.0E-6 ! kernel is needed in m**3/s 279 280 DEALLOCATE( ec, radclass, winf ) 281 282 END SUBROUTINE recalculate_kernel 283 284 285 !------------------------------------------------------------------------------! 286 ! Calculation of gck 287 ! This is from Aayala 2008b, page 37ff. 288 ! Necessary input parameters: water density, radii of droplets, air density, 289 ! air viscosity, turbulent dissipation rate, taylor microscale reynolds number, 290 ! gravitational acceleration --> to be replaced by PALM parameters 175 291 !------------------------------------------------------------------------------! 176 292 SUBROUTINE turbsd 177 ! from Aayala 2008b, page 37ff, necessary input parameter water density, radii178 ! of droplets, air density, air viscosity, turbulent dissipation rate,179 ! taylor microscale reynolds number, gravitational acceleration180 293 181 294 USE constants … … 186 299 IMPLICIT NONE 187 300 188 REAL :: Relamda, & 189 Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be, & 190 bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq, & 191 vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, & 192 wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp 193 194 REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens 195 196 REAL, DIMENSION(pstart:pend) :: St, tau 301 INTEGER :: i, j 197 302 198 303 LOGICAL, SAVE :: first = .TRUE. 199 304 200 INTEGER :: i, j 201 202 ! 203 !-- Initial assignment of constants 305 REAL :: ao, ao_gr, bbb, be, b1, b2, ccc, c1, c1_gr, c2, d1, d2, eta, & 306 e1, e2, fao_gr, fr, grfin, lambda, lambda_re, lf, rc, rrp, & 307 sst, tauk, tl, t2, tt, t1, vk, vrms1xy, vrms2xy, v1, v1v2xy, & 308 v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z 309 310 REAL, SAVE :: airdens, airvisc, anu, gravity, waterdens 311 312 REAL, DIMENSION(1:radius_classes) :: st, tau 313 314 315 ! 316 !-- Initial assignment of constants 204 317 IF ( first ) THEN 205 318 206 319 first = .FALSE. 207 airvisc = 0.1818 !dynamic viscosity in mg/cm*s208 airdens = 1.2250 !air density in mg/cm**3209 waterdens = 1000.0 !water density in mg/cm**3210 gravity = 980.6650 !in cm/s**2320 airvisc = 0.1818 ! dynamic viscosity in mg/cm*s 321 airdens = 1.2250 ! air density in mg/cm**3 322 waterdens = 1000.0 ! water density in mg/cm**3 323 gravity = 980.6650 ! in cm/s**2 211 324 anu = airvisc/airdens ! kinetic viscosity in cm**2/s 212 325 213 326 ENDIF 214 327 215 Relamda = urms**2*sqrt(15.0/epsilon/anu) 216 217 Tl = urms**2/epsilon !in s 218 Lf = 0.5 * (urms**3)/epsilon !in cm 219 220 tauk = (anu/epsilon)**0.5 !in s 221 eta = (anu**3/epsilon)**0.25 !in cm 222 vk = eta/tauk !in cm/s 223 224 ao = (11.+7.*Relamda)/(205.+Relamda) 225 226 lambda = urms * sqrt(15.*anu/epsilon) !in cm 227 228 tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s 229 230 CALL fallg !gives winf in cm/s 231 232 DO i = pstart, pend 233 tau(i) = winf(i)/gravity !in s 234 St(i) = tau(i)/tauk 328 lambda = urms * SQRT( 15.0 * anu / epsilon ) ! in cm 329 lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu ) 330 tl = urms**2 / epsilon ! in s 331 lf = 0.5 * urms**3 / epsilon ! in cm 332 tauk = SQRT( anu / epsilon ) ! in s 333 eta = ( anu**3 / epsilon )**0.25 ! in cm 334 vk = eta / tauk ! in cm/s 335 336 ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re ) 337 tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk ! in s 338 339 CALL fallg ! gives winf in cm/s 340 341 DO i = 1, radius_classes 342 tau(i) = winf(i) / gravity ! in s 343 st(i) = tau(i) / tauk 235 344 ENDDO 236 345 237 !***** TO CALCULATE wr ******************************** 238 !from Aayala 2008b, page 38f 239 240 z = tt/Tl 241 be = sqrt(2.0)*lambda/Lf 242 243 bbb = sqrt(1.0-2.0*be**2) 244 d1 = (1.+bbb)/2.0/bbb 245 e1 = Lf*(1.0+bbb)/2.0 !in cm 246 d2 = (1.0-bbb)/2.0/bbb 247 e2 = Lf*(1.0-bbb)/2.0 !in cm 248 249 ccc = sqrt(1.0-2.0*z**2) 250 b1 = (1.+ccc)/2./ccc 251 c1 = Tl*(1.+ccc)/2. !in s 252 b2 = (1.-ccc)/2./ccc 253 c2 = Tl*(1.-ccc)/2. !in s 254 255 DO i = pstart, pend 256 v1 = winf(i) !in cm/s 257 t1 = tau(i) !in s 258 259 DO j = pstart,i 260 rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm 261 v2 = winf(j) !in cm/s 262 t2 = tau(j) !in s 263 264 v1xysq = b1*d1*PHI(c1,e1,v1,t1) & 265 - b1*d2*PHI(c1,e2,v1,t1) & 266 - b2*d1*PHI(c2,e1,v1,t1) & 267 + b2*d2*PHI(c2,e2,v1,t1) 268 v1xysq = v1xysq * urms**2/t1 !in cm**2/s**2 269 270 vrms1xy= sqrt(v1xysq) !in cm/s 271 272 v2xysq = b1*d1*PHI(c1,e1,v2,t2) & 273 - b1*d2*PHI(c1,e2,v2,t2) & 274 - b2*d1*PHI(c2,e1,v2,t2) & 275 + b2*d2*PHI(c2,e2,v2,t2) 276 v2xysq = v2xysq * urms**2/t2 !in cm**2/s**2 277 278 vrms2xy= sqrt(v2xysq) !in cm/s 279 280 IF(winf(i).ge.winf(j)) THEN 346 ! 347 !-- Calculate wr (from Aayala 2008b, page 38f) 348 z = tt / tl 349 be = SQRT( 2.0 ) * lambda / lf 350 bbb = SQRT( 1.0 - 2.0 * be**2 ) 351 d1 = ( 1.0 + bbb ) / ( 2.0 * bbb ) 352 e1 = lf * ( 1.0 + bbb ) * 0.5 ! in cm 353 d2 = ( 1.0 - bbb ) * 0.5 / bbb 354 e2 = lf * ( 1.0 - bbb ) * 0.5 ! in cm 355 ccc = SQRT( 1.0 - 2.0 * z**2 ) 356 b1 = ( 1.0 + ccc ) * 0.5 / ccc 357 c1 = tl * ( 1.0 + ccc ) * 0.5 ! in s 358 b2 = ( 1.0 - ccc ) * 0.5 / ccc 359 c2 = tl * ( 1.0 - ccc ) * 0.5 ! in s 360 361 DO i = 1, radius_classes 362 363 v1 = winf(i) ! in cm/s 364 t1 = tau(i) ! in s 365 366 DO j = 1, i 367 rrp = radclass(i) + radclass(j) ! radius in cm 368 v2 = winf(j) ! in cm/s 369 t2 = tau(j) ! in s 370 371 v1xysq = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) & 372 - b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1) 373 v1xysq = v1xysq * urms**2 / t1 ! in cm**2/s**2 374 vrms1xy = SQRT( v1xysq ) ! in cm/s 375 376 v2xysq = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) & 377 - b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2) 378 v2xysq = v2xysq * urms**2 / t2 ! in cm**2/s**2 379 vrms2xy = SQRT( v2xysq ) ! in cm/s 380 381 IF ( winf(i) >= winf(j) ) THEN 281 382 v1 = winf(i) 282 383 t1 = tau(i) … … 290 391 ENDIF 291 392 292 v1v2xy = b1*d1*ZHI(c1,e1,v1,t1,v2,t2) & 293 - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) & 294 - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) & 295 + b2*d2*ZHI(c2,e2,v1,t1,v2,t2) 296 fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2) 297 v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j) !in cm**2/s**2 298 299 wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy !in cm**2/s**2 300 301 IF (wrtur2xy.le.0.0) wrtur2xy=0.0 302 303 wrgrav2=pi/8.*(winf(j)-winf(i))**2 304 305 wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2)) !in cm/s 306 307 308 !***** TO CALCULATE gr ******************************** 309 310 IF(St(j).gt.St(i)) THEN 311 SSt = St(j) 393 v1v2xy = b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - & 394 b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - & 395 b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + & 396 b2 * d2* zhi(c2,e2,v1,t1,v2,t2) 397 fr = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 ) 398 v1v2xy = v1v2xy * fr * urms**2 / tau(i) / tau(j) ! in cm**2/s**2 399 wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy ! in cm**2/s**2 400 IF ( wrtur2xy < 0.0 ) wrtur2xy = 0.0 401 wrgrav2 = pi / 8.0 * ( winf(j) - winf(i) )**2 402 wrfin = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) ) ! in cm/s 403 404 ! 405 !-- Calculate gr 406 IF ( st(j) > st(i) ) THEN 407 sst = st(j) 312 408 ELSE 313 SSt = St(i)409 sst = st(i) 314 410 ENDIF 315 411 316 XX = -0.1988*SSt**4 + 1.5275*SSt**3 & 317 -4.2942*SSt**2 + 5.3406*SSt 318 319 IF(XX.le.0.0) XX = 0.0 320 321 YY = 0.1886*exp(20.306/Relamda) 322 323 c1_gr = XX/(gravity/(vk/tauk))**YY 324 325 ao_gr = ao + (pi/8.)*(gravity/(vk/tauk))**2 326 fao_gr = 20.115 * (ao_gr/Relamda)**0.5 327 rc = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta !in cm 328 329 grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.) 330 IF (grFIN.lt.1.0) grFIN = 1.0 331 332 gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN ! in cm**3/s 412 xx = -0.1988 * sst**4 + 1.5275 * sst**3 - 4.2942 * sst**2 + & 413 5.3406 * sst 414 IF ( xx < 0.0 ) xx = 0.0 415 yy = 0.1886 * EXP( 20.306 / lambda_re ) 416 417 c1_gr = xx / ( gravity / vk * tauk )**yy 418 419 ao_gr = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2 420 fao_gr = 20.115 * SQRT( ao_gr / lambda_re ) 421 rc = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta ! in cm 422 423 grfin = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5 ) 424 IF ( grfin < 1.0 ) grfin = 1.0 425 426 gck(i,j) = 2.0 * pi * rrp**2 * wrfin * grfin ! in cm**3/s 333 427 gck(j,i) = gck(i,j) 334 428 … … 336 430 ENDDO 337 431 338 END SUBROUTINE TurbSD 339 340 !------------------------------------------------------------------------------! 341 ! PHI as a function 342 !------------------------------------------------------------------------------! 343 REAL FUNCTION PHI(a,b,vsett,tau0) 432 END SUBROUTINE turbsd 433 434 435 !------------------------------------------------------------------------------! 436 ! phi as a function 437 !------------------------------------------------------------------------------! 438 REAL FUNCTION phi( a, b, vsett, tau0 ) 344 439 345 440 IMPLICIT NONE 346 441 347 REAL :: a, aa1, b, vsett, tau0348 349 aa1 = 1. /tau0 + 1./a + vsett/b350 351 PHI = 1./aa1 - vsett/2.0/b/aa1**2 !in s 352 353 END FUNCTION PHI 354 355 !------------------------------------------------------------------------------! 356 ! ZETAas a function357 !------------------------------------------------------------------------------! 358 REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)442 REAL :: a, aa1, b, tau0, vsett 443 444 aa1 = 1.0 / tau0 + 1.0 / a + vsett / b 445 phi = 1.0 / aa1 - 0.5 * vsett / b / aa1**2 ! in s 446 447 END FUNCTION phi 448 449 450 !------------------------------------------------------------------------------! 451 ! zeta as a function 452 !------------------------------------------------------------------------------! 453 REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 ) 359 454 360 455 IMPLICIT NONE 361 456 362 REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2 363 364 aa1 = vsett2/b - 1./tau2 - 1./a 365 aa2 = vsett1/b + 1./tau1 + 1./a 366 aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2 367 aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2 368 aa5 = vsett2/b + 1./tau2 + 1./a 369 aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2 370 ZHI = (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2 & 371 + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6 & 372 + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2) & 373 * 1./2./b/aa3 ! in s**2 374 375 END FUNCTION ZHI 376 377 !------------------------------------------------------------------------------! 378 ! SUBROUTINE for calculation of terminal velocity winf 379 !------------------------------------------------------------------------------! 380 SUBROUTINE fallg 457 REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, tau1, tau2, vsett1, vsett2 458 459 aa1 = vsett2 / b - 1.0 / tau2 - 1.0 / a 460 aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a 461 aa3 = ( vsett1 - vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2 462 aa4 = ( vsett2 / b )**2 - ( 1.0 / tau2 + 1.0 / a )**2 463 aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a 464 aa6 = 1.0 / tau1 - 1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2 465 zhi = (1.0 / aa1 - 1.0 / aa2 ) * ( vsett1 - vsett2 ) * 0.5 / b / aa3**2 & 466 + (4.0 / aa4 - 1.0 / aa5**2 - 1.0 / aa1**2 ) * vsett2 * 0.5 / b /aa6& 467 + (2.0 * ( b / aa2 - b / aa1 ) - vsett1 / aa2**2 + vsett2 / aa1**2 )& 468 * 0.5 / b / aa3 ! in s**2 469 470 END FUNCTION zhi 471 472 473 !------------------------------------------------------------------------------! 474 ! Calculation of terminal velocity winf 475 !------------------------------------------------------------------------------! 476 SUBROUTINE fallg 381 477 382 478 USE constants … … 385 481 USE arrays_3d 386 482 387 IMPLICIT NONE 388 389 INTEGER :: i, j 390 391 LOGICAL, SAVE :: first = .TRUE. 392 393 REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py 394 395 REAL :: bond, x, xrey, y 396 397 REAL, DIMENSION(1:7), SAVE :: b 398 REAL, DIMENSION(1:6), SAVE :: c 399 400 ! 401 !-- Initial assignment of constants 402 IF ( first ) THEN 483 IMPLICIT NONE 484 485 INTEGER :: i, j 486 487 LOGICAL, SAVE :: first = .TRUE. 488 489 REAL, SAVE :: cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, & 490 t0, xlamb 491 492 REAL :: bond, x, xrey, y 493 494 REAL, DIMENSION(1:7), SAVE :: b 495 REAL, DIMENSION(1:6), SAVE :: c 496 497 ! 498 !-- Initial assignment of constants 499 IF ( first ) THEN 500 501 first = .FALSE. 502 b = (/ -0.318657E1, 0.992696E0, -0.153193E-2, -0.987059E-3, & 503 -0.578878E-3, 0.855176E-4, -0.327815E-5 /) 504 c = (/ -0.500015E1, 0.523778E1, -0.204914E1, 0.475294E0, & 505 -0.542819E-1, 0.238449E-2 /) 506 507 eta = 1.818E-4 ! in poise = g/(cm s) 508 xlamb = 6.62E-6 ! in cm 509 rhow = 1.0 ! in g/cm**3 510 rhoa = 1.225E-3 ! in g/cm**3 511 grav = 980.665 ! in cm/s**2 512 cunh = 1.257 * xlamb ! in cm 513 t0 = 273.15 ! in K 514 sigma = 76.1 - 0.155 * ( 293.15 - t0 ) ! in N/m = g/s**2 515 stok = 2.0 * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s) 516 stb = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta) 517 ! in 1/cm**3 518 phy = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) ) 519 py = phy**( 1.0 / 6.0 ) 520 521 ENDIF 522 523 DO j = 1, radius_classes 524 525 IF ( radclass(j) <= 1.0E-3 ) THEN 526 527 winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s 528 529 ELSEIF ( radclass(j) > 1.0E-3 .AND. radclass(j) <= 5.35E-2 ) THEN 530 531 x = LOG( stb * radclass(j)**3 ) 532 y = 0.0 533 534 DO i = 1, 7 535 y = y + b(i) * x**(i-1) 536 ENDDO 537 538 xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y ) 539 winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) ) ! in cm/s 540 541 ELSEIF ( radclass(j) > 5.35E-2 ) THEN 542 543 IF ( radclass(j) > 0.35 ) THEN 544 bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma 545 ELSE 546 bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma 547 ENDIF 548 549 x = LOG( 16.0 * bond * py / 3.0 ) 550 y = 0.0 551 552 DO i = 1, 6 553 y = y + c(i) * x**(i-1) 554 ENDDO 555 556 xrey = py * EXP( y ) 557 558 IF ( radclass(j) > 0.35 ) THEN 559 winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 ) ! in cm/s 560 ELSE 561 winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) ) ! in cm/s 562 ENDIF 563 564 ENDIF 565 566 ENDDO 567 568 END SUBROUTINE fallg 569 570 571 !------------------------------------------------------------------------------! 572 ! Calculation of collision efficencies for the Hall kernel 573 !------------------------------------------------------------------------------! 574 SUBROUTINE effic 575 576 USE arrays_3d 577 USE cloud_parameters 578 USE constants 579 USE particle_attributes 580 581 IMPLICIT NONE 582 583 INTEGER :: i, iq, ir, j, k, kk 584 585 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 586 587 LOGICAL, SAVE :: first = .TRUE. 588 589 REAL :: ek, particle_radius, pp, qq, rq 590 591 REAL, DIMENSION(1:21), SAVE :: rat 592 REAL, DIMENSION(1:15), SAVE :: r0 593 REAL, DIMENSION(1:15,1:21), SAVE :: ecoll 594 595 ! 596 !-- Initial assignment of constants 597 IF ( first ) THEN 403 598 404 599 first = .FALSE. 405 b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, & 406 -0.578878e-3,0.855176e-4,-0.327815e-5/) 407 c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0, & 408 -0.542819e-1, 0.238449e-2/) 409 410 eta = 1.818e-4 !in poise = g/(cm s) 411 xlamb = 6.62e-6 !in cm 412 413 rhow = 1.0 !in g/cm**3 414 rhoa = 1.225e-3 !in g/cm**3 415 416 grav = 980.665 !in cm/s**2 417 cunh = 1.257 * xlamb !in cm 418 t0 = 273.15 !in K 419 sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2 420 stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s) 421 stb = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3 422 phy = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa)) 423 py = phy**(1.0/6.0) 424 425 ENDIF 426 427 !particle radius has to be in cm 428 DO j = pstart, pend 429 430 IF (particles(j)%radius*100.0 .le. 1.e-3) THEN 431 432 winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s 433 434 ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN 435 436 x = log(stb*(particles(j)%radius*100.0)**3) 437 y = 0.0 438 439 DO i = 1, 7 440 y = y + b(i) * (x**(i-1)) 441 ENDDO 442 443 xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y) 444 winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s 445 446 ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN 447 448 IF (particles(j)%radius*100.0.gt.0.35) THEN 449 bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma 450 ELSE 451 bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma 452 ENDIF 453 454 x = log(16.0*bond*py/3.0) 455 y = 0.0 456 457 DO i = 1, 6 458 y = y + c(i) * (x**(i-1)) 459 ENDDO 460 461 xrey = py*exp(y) 462 463 IF (particles(j)%radius*100.0 .gt.0.35) THEN 464 winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s 465 ELSE 466 winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s 467 ENDIF 468 469 ENDIF 470 ENDDO 471 RETURN 472 END SUBROUTINE fallg 473 474 !------------------------------------------------------------------------------! 475 ! SUBROUTINE for calculation of collision efficencies 476 !------------------------------------------------------------------------------! 477 478 SUBROUTINE effic 479 480 USE arrays_3d 481 USE constants 482 USE cloud_parameters 483 USE particle_attributes 484 485 !collision efficiencies of hall kernel 486 IMPLICIT NONE 487 488 INTEGER :: i, ir, iq, j, k, kk 489 490 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 491 492 LOGICAL, SAVE :: first = .TRUE. 493 494 REAL :: ek, particle_radius, pp, qq, rq 495 496 REAL, DIMENSION(1:21), SAVE :: rat 497 REAL, DIMENSION(1:15), SAVE :: r0 498 REAL, DIMENSION(1:15,1:21), SAVE :: ecoll 499 500 ! 501 !-- Initial assignment of constants 502 IF ( first ) THEN 503 504 first = .FALSE. 505 r0 = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., & 506 300./) 507 rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6, & 508 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/) 509 510 ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, & 511 0.001,0.001,0.001,0.001,0.001,0.001/) 512 ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, & 513 0.050,0.200,0.500,0.770,0.870,0.970/) 514 ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, & 515 0.430,0.580,0.790,0.930,0.960,1.000/) 516 ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, & 517 0.640,0.750,0.910,0.970,0.980,1.000/) 518 ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, & 519 0.770,0.840,0.950,0.970,1.000,1.000/) 520 ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, & 521 0.840,0.880,0.950,1.000,1.000,1.000/) 522 ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, & 523 0.870,0.900,0.950,1.000,1.000,1.000/) 524 ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, & 525 0.890,0.920,1.000,1.000,1.000,1.000/) 526 ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, & 527 0.900,0.940,1.000,1.000,1.000,1.000/) 528 ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, & 529 0.910,0.950,1.000,1.000,1.000,1.000/) 530 ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, & 531 0.910,0.950,1.000,1.000,1.000,1.000/) 532 ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, & 533 0.910,0.950,1.000,1.000,1.000,1.000/) 534 ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, & 535 0.910,0.950,1.000,1.000,1.000,1.000/) 536 ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, & 537 0.910,0.950,1.000,1.000,1.000,1.000/) 538 ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, & 539 0.920,0.950,1.000,1.000,1.000,1.000/) 540 ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, & 541 0.930,0.970,1.000,1.000,1.000,1.000/) 542 ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, & 543 0.950,1.000,1.000,1.000,1.000,1.000/) 544 ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, & 545 1.000,1.020,1.020,1.020,1.020,1.020/) 546 ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, & 547 1.030,1.040,1.040,1.040,1.040,1.040/) 548 ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, & 549 1.700,2.300,2.300,2.300,2.300,2.300/) 550 ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, & 551 3.000,4.000,4.000,4.000,4.000,4.000/) 552 ENDIF 553 554 ! 555 !-- Calculate the radius class index of particles with respect to array r 556 ALLOCATE( ira(pstart:pend) ) 557 DO j = pstart, pend 558 particle_radius = particles(j)%radius * 1.0E6 559 DO k = 1, 15 560 IF ( particle_radius < r0(k) ) THEN 561 ira(j) = k 562 EXIT 563 ENDIF 564 ENDDO 565 IF ( particle_radius >= r0(15) ) ira(j) = 16 566 ENDDO 567 568 ! 569 !-- Two-dimensional linear interpolation of the collision efficiency. 570 !-- Radius has to be in µm 571 DO j = pstart, pend 572 DO i = pstart, j 573 574 ir = ira(j) 575 576 rq = particles(i)%radius / particles(j)%radius 577 578 ! DO kk = 2, 21 579 ! IF ( rq <= rat(kk) ) THEN 580 ! iq = kk 581 ! EXIT 582 ! ENDIF 583 ! ENDDO 584 585 iq = INT( rq * 20 ) + 1 586 iq = MAX(iq , 2) 587 588 IF ( ir < 16 ) THEN 589 590 IF ( ir >= 2 ) THEN 591 pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / & 592 ( r0(ir) - r0(ir-1) ) 593 qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 594 ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1) & 595 + pp * ( 1.0-qq ) * ecoll(ir,iq-1) & 596 + qq * ( 1.0-pp ) * ecoll(ir-1,iq) & 597 + pp * qq * ecoll(ir,iq) 598 ELSE 599 qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) 600 ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq) 601 ENDIF 602 603 ELSE 604 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 605 ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq) 606 ec(j,i) = MIN( ek, 1.0 ) 607 ENDIF 608 609 ec(i,j) = ec(j,i) 610 IF ( ec(i,j) < 1.0E-20 ) ec(i,j) = 0.0 611 612 ENDDO 613 ENDDO 614 615 DEALLOCATE( ira ) 616 617 END SUBROUTINE effic 618 619 !------------------------------------------------------------------------------! 620 ! SUBROUTINE for calculation of enhancement factor collision efficencies 621 !------------------------------------------------------------------------------! 622 SUBROUTINE turb_enhance_eff 600 r0 = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., & 601 70.0, 100.0, 150.0, 200.0, 300.0 /) 602 rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, & 603 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, & 604 1.00 /) 605 606 ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, & 607 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/) 608 ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, & 609 0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/) 610 ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, & 611 0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/) 612 ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, & 613 0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/) 614 ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, & 615 0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/) 616 ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, & 617 0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/) 618 ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, & 619 0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/) 620 ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, & 621 0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/) 622 ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, & 623 0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/) 624 ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, & 625 0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 626 ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, & 627 0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 628 ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, & 629 0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 630 ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, & 631 0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 632 ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, & 633 0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 634 ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, & 635 0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/) 636 ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, & 637 0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/) 638 ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, & 639 0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/) 640 ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, & 641 0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/) 642 ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, & 643 0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/) 644 ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, & 645 0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/) 646 ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, & 647 1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/) 648 ENDIF 649 650 ! 651 !-- Calculate the radius class index of particles with respect to array r 652 ALLOCATE( ira(1:radius_classes) ) 653 DO j = 1, radius_classes 654 particle_radius = radclass(j) * 1.0E4 ! in microm 655 DO k = 1, 15 656 IF ( particle_radius < r0(k) ) THEN 657 ira(j) = k 658 EXIT 659 ENDIF 660 ENDDO 661 IF ( particle_radius >= r0(15) ) ira(j) = 16 662 ENDDO 663 664 ! 665 !-- Two-dimensional linear interpolation of the collision efficiency. 666 !-- Radius has to be in µm 667 DO j = 1, radius_classes 668 DO i = 1, j 669 670 ir = ira(j) 671 rq = radclass(i) / radclass(j) 672 iq = INT( rq * 20 ) + 1 673 iq = MAX( iq , 2) 674 675 IF ( ir < 16 ) THEN 676 IF ( ir >= 2 ) THEN 677 pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / & 678 ( r0(ir) - r0(ir-1) ) 679 qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 680 ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1) & 681 + pp * ( 1.0-qq ) * ecoll(ir,iq-1) & 682 + qq * ( 1.0-pp ) * ecoll(ir-1,iq) & 683 + pp * qq * ecoll(ir,iq) 684 ELSE 685 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 686 ec(j,i) = (1.0-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq) 687 ENDIF 688 ELSE 689 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 690 ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq) 691 ec(j,i) = MIN( ek, 1.0 ) 692 ENDIF 693 694 ec(i,j) = ec(j,i) 695 IF ( ec(i,j) < 1.0E-20 ) ec(i,j) = 0.0 696 697 ENDDO 698 ENDDO 699 700 DEALLOCATE( ira ) 701 702 END SUBROUTINE effic 703 704 705 !------------------------------------------------------------------------------! 706 ! Calculation of enhancement factor for collision efficencies due to turbulence 707 !------------------------------------------------------------------------------! 708 SUBROUTINE turb_enhance_eff 623 709 624 710 USE constants … … 627 713 USE arrays_3d 628 714 629 IMPLICIT NONE 630 631 INTEGER :: i, ik, ir, iq, j, k, kk 632 633 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 634 635 REAL :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3 636 637 LOGICAL, SAVE :: first = .TRUE. 638 639 REAL, DIMENSION(1:11), SAVE :: rat 640 REAL, DIMENSION(1:7), SAVE :: r0 641 REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400 642 643 ! 644 !-- Initial assignment of constants 645 IF ( first ) THEN 646 647 first = .FALSE. 648 649 r0 = (/10., 20., 30.,40., 50., 60.,100./) 650 rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/) 651 652 ! 100 cm^2/s^3 653 ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /) 654 ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /) 655 ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /) 656 ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /) 657 ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /) 658 ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /) 659 ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /) 660 ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /) 661 ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /) 662 ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /) 663 ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /) 664 665 ! 400 cm^2/s^3 666 ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) 667 ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) 668 ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /) 669 ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /) 670 ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /) 671 ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /) 672 ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /) 673 ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /) 674 ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /) 675 ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /) 676 ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /) 677 678 ENDIF 679 680 ! 681 !-- Calculate the radius class index of particles with respect to array r 682 ALLOCATE( ira(pstart:pend) ) 683 684 DO j = pstart, pend 685 particle_radius = particles(j)%radius * 1.0E6 686 DO k = 1, 7 687 IF ( particle_radius < r0(k) ) THEN 688 ira(j) = k 689 EXIT 690 ENDIF 691 ENDDO 692 IF ( particle_radius >= r0(7) ) ira(j) = 8 693 ENDDO 694 695 ! two-dimensional linear interpolation of the collision efficiency 696 DO j = pstart, pend 697 DO i = pstart, j 698 699 ir = ira(j) 700 701 rq = particles(i)%radius/particles(j)%radius 702 703 DO kk = 2, 11 704 IF ( rq <= rat(kk) ) THEN 705 iq = kk 706 EXIT 707 ENDIF 708 ENDDO 709 710 ! 0 cm2/s3 711 y1 = 1.0 712 ! 100 cm2/s3, 400 cm2/s3 713 IF (ir.lt.8) THEN 714 IF (ir.ge.2) THEN 715 pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1)) 716 qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) 717 y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+ & 718 pp*(1.-qq)*ecoll_100(ir,iq-1)+ & 719 qq*(1.-pp)*ecoll_100(ir-1,iq)+ & 720 pp*qq*ecoll_100(ir,iq) 721 y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+ & 722 pp*(1.-qq)*ecoll_400(ir,iq-1)+ & 723 qq*(1.-pp)*ecoll_400(ir-1,iq)+ & 724 pp*qq*ecoll_400(ir,iq) 725 ELSE 726 qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) 727 y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq) 728 y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq) 729 ENDIF 730 ELSE 731 qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) 732 y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq) 733 y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq) 734 ENDIF 735 ! linear interpolation 736 ! dissipation rate in cm2/s3 737 x1 = 0.0 738 x2 = 100.0 739 x3 = 400.0 740 741 IF (epsilon.le.100.) THEN 742 ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 & 743 + (epsilon-0.)/(100.-0.) * y2 744 ELSE IF(epsilon.le.600.)THEN 745 ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 & 746 + (epsilon-100.)/(400.-100.) * y3 747 748 ELSE 749 ecf(j,i) = (600.-400.)/(100.-400.) * y2 & 750 + (600.-100.)/(400.-100.) * y3 751 ENDIF 752 753 IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0 754 755 ecf(i,j)=ecf(j,i) 756 ENDDO 757 ENDDO 758 759 RETURN 760 END SUBROUTINE turb_enhance_eff 715 IMPLICIT NONE 716 717 INTEGER :: i, ik, iq, ir, j, k, kk 718 719 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 720 721 REAL :: particle_radius, pp, qq, rq, x1, x2, x3, y1, y2, y3 722 723 LOGICAL, SAVE :: first = .TRUE. 724 725 REAL, DIMENSION(1:11), SAVE :: rat 726 REAL, DIMENSION(1:7), SAVE :: r0 727 REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400 728 729 ! 730 !-- Initial assignment of constants 731 IF ( first ) THEN 732 733 first = .FALSE. 734 735 r0 = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /) 736 rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /) 737 ! 738 !-- In 100 cm^2/s^3 739 ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /) 740 ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /) 741 ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /) 742 ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /) 743 ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /) 744 ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /) 745 ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /) 746 ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /) 747 ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /) 748 ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /) 749 ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /) 750 ! 751 !-- In 400 cm^2/s^3 752 ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) 753 ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) 754 ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /) 755 ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /) 756 ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /) 757 ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /) 758 ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /) 759 ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /) 760 ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /) 761 ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /) 762 ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /) 763 764 ENDIF 765 766 ! 767 !-- Calculate the radius class index of particles with respect to array r0 768 ALLOCATE( ira(1:radius_classes) ) 769 770 DO j = 1, radius_classes 771 particle_radius = radclass(j) * 1.0E4 ! in microm 772 DO k = 1, 7 773 IF ( particle_radius < r0(k) ) THEN 774 ira(j) = k 775 EXIT 776 ENDIF 777 ENDDO 778 IF ( particle_radius >= r0(7) ) ira(j) = 8 779 ENDDO 780 781 ! 782 !-- Two-dimensional linear interpolation of the collision efficiencies 783 DO j = 1, radius_classes 784 DO i = 1, j 785 786 ir = ira(j) 787 rq = radclass(i) / radclass(j) 788 789 DO kk = 2, 11 790 IF ( rq <= rat(kk) ) THEN 791 iq = kk 792 EXIT 793 ENDIF 794 ENDDO 795 796 y1 = 1.0 ! 0 cm2/s3 797 ! 798 !-- 100 cm2/s3, 400 cm2/s3 799 IF ( ir < 8 ) THEN 800 IF ( ir >= 2 ) THEN 801 pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) ) 802 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 803 y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) + & 804 pp * ( 1.0-qq ) * ecoll_100(ir,iq-1) + & 805 qq * ( 10.-pp ) * ecoll_100(ir-1,iq) + & 806 pp * qq * ecoll_100(ir,iq) 807 y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) + & 808 pp * ( 1.0-qq ) * ecoll_400(ir,iq-1) + & 809 qq * ( 1.0-pp ) * ecoll_400(ir-1,iq) + & 810 pp * qq * ecoll_400(ir,iq) 811 ELSE 812 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 813 y2 = ( 1.0-qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq) 814 y3 = ( 1.0-qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq) 815 ENDIF 816 ELSE 817 qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) 818 y2 = ( 1.0-qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq) 819 y3 = ( 1.0-qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq) 820 ENDIF 821 ! 822 !-- Linear interpolation of dissipation rate in cm2/s3 823 IF ( epsilon <= 100.0 ) THEN 824 ecf(j,i) = ( epsilon - 100.0 ) / ( 0.0 - 100.0 ) * y1 & 825 + ( epsilon - 0.0 ) / ( 100.0 - 0.0 ) * y2 826 ELSEIF ( epsilon <= 600.0 ) THEN 827 ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 & 828 + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3 829 ELSE 830 ecf(j,i) = ( 600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 & 831 + ( 600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3 832 ENDIF 833 834 IF ( ecf(j,i) < 1.0 ) ecf(j,i) = 1.0 835 836 ecf(i,j) = ecf(j,i) 837 838 ENDDO 839 ENDDO 840 841 END SUBROUTINE turb_enhance_eff 761 842 762 843 END MODULE lpm_collision_kernels_mod -
palm/trunk/SOURCE/modules.f90
r826 r828 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +dissipation_classes, radius_classes, use_kernel_tables, 7 ! particle feature color renamed class 7 8 ! 8 9 ! Former revisions: … … 1195 1196 INTEGER :: mpi_particle_type 1196 1197 #endif 1197 INTEGER :: ibc_par_lr, ibc_par_ns, ibc_par_b, ibc_par_t, & 1198 iran_part = -1234567, maximum_number_of_particles = 1000, & 1198 INTEGER :: dissipation_classes = 10, ibc_par_lr, ibc_par_ns, ibc_par_b, & 1199 ibc_par_t, iran_part = -1234567, & 1200 maximum_number_of_particles = 1000, & 1199 1201 maximum_number_of_tailpoints = 100, & 1200 1202 maximum_number_of_tails = 0, & … … 1203 1205 number_of_initial_tails = 0, offset_ocean_nzt = 0, & 1204 1206 offset_ocean_nzt_m1 = 0, particles_per_point = 1, & 1205 particle_file_count = 0, skip_particles_for_tail = 100,&1206 s ort_count = 0, total_number_of_particles,&1207 total_number_of_ tails = 01207 particle_file_count = 0, radius_classes = 20, & 1208 skip_particles_for_tail = 100, sort_count = 0, & 1209 total_number_of_particles, total_number_of_tails = 0 1208 1210 1209 1211 INTEGER, PARAMETER :: max_number_of_particle_groups = 10 … … 1215 1217 particle_advection = .FALSE., random_start_position = .FALSE., & 1216 1218 read_particles_from_restartfile = .TRUE., & 1217 uniform_particles = .TRUE., use_ particle_tails = .FALSE.,&1218 use_ sgs_for_particles = .FALSE.,&1219 uniform_particles = .TRUE., use_kernel_tables = .FALSE., & 1220 use_particle_tails = .FALSE., use_sgs_for_particles = .FALSE., & 1219 1221 wang_kernel = .FALSE., write_particle_statistics = .FALSE. 1220 1222 … … 1247 1249 origin_z, radius, rvar1, rvar2, rvar3, speed_x, speed_y, & 1248 1250 speed_z, weight_factor, x, y, z 1249 INTEGER :: c olor, group, tailpoints, tail_id1251 INTEGER :: class, group, tailpoints, tail_id 1250 1252 END TYPE particle_type 1251 1253 -
palm/trunk/SOURCE/package_parin.f90
r826 r828 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +dissipation_classes, radius_classes in parpar 7 7 ! 8 8 ! Former revisions: … … 75 75 76 76 NAMELIST /particles_par/ bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, & 77 collision_kernel, 78 d ensity_ratio, radius, dt_dopts,&77 collision_kernel, density_ratio, & 78 dissipation_classes, dt_dopts, & 79 79 dt_min_part, dt_prel, dt_sort_particles, & 80 80 dt_write_particle_data, & … … 88 88 particle_advection_start, & 89 89 particle_maximum_age, pdx, pdy, pdz, psb, & 90 psl, psn, psr, pss, pst, 91 ra ndom_start_position,&90 psl, psn, psr, pss, pst, radius, & 91 radius_classes, random_start_position, & 92 92 read_particles_from_restartfile, & 93 93 skip_particles_for_tail, & -
palm/trunk/SOURCE/set_particle_attributes.f90
r623 r828 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! particle feature color renamed class 7 7 ! 8 8 ! Former revisions: … … 121 121 ! 122 122 !-- Number of available colors is defined in init_dvrp 123 particles(n)%c olor= 1 + absuv * ( dvrp_colortable_entries_prt - 1 )123 particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 ) 124 124 125 125 ENDDO … … 192 192 ! 193 193 !-- Number of available colors is defined in init_dvrp 194 particles(n)%c olor= 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )194 particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 ) 195 195 196 196 ENDDO … … 213 213 ! 214 214 !-- Number of available colors is defined in init_dvrp 215 particles(n)%c olor= 1 + height * ( dvrp_colortable_entries_prt - 1 )215 particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 ) 216 216 217 217 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.