Changeset 4145
- Timestamp:
- Aug 6, 2019 9:55:22 AM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
r4144 r4145 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Some reformatting 28 ! 29 ! 4144 2019-08-06 09:11:47Z raasch 27 30 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc. 28 31 ! … … 387 390 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ckernel !< 388 391 389 INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !<390 INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2 !<392 INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !< 393 INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2 !< 391 394 392 395 SAVE … … 871 874 !-- it can be combined as the sgs-velocites for active particles are 872 875 !-- calculated differently, i.e. no subboxes are needed. 873 IF ( .NOT. TRIM( particle_advection_interpolation) == 'trilinear' .AND.&874 use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN876 IF ( .NOT. TRIM( particle_advection_interpolation ) == 'trilinear' .AND. & 877 use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 875 878 message_string = 'subrgrid scale velocities in combination with ' // & 876 879 'simple interpolation method is not ' // & … … 1076 1079 ! 1077 1080 !-- Check which particle interpolation method should be used 1078 IF ( TRIM( particle_advection_interpolation) == 'trilinear' ) THEN1081 IF ( TRIM( particle_advection_interpolation ) == 'trilinear' ) THEN 1079 1082 interpolation_simple_corrector = .FALSE. 1080 1083 interpolation_simple_predictor = .FALSE. 1081 1084 interpolation_trilinear = .TRUE. 1082 ELSEIF ( TRIM( particle_advection_interpolation) == 'simple_corrector' ) THEN1085 ELSEIF ( TRIM( particle_advection_interpolation ) == 'simple_corrector' ) THEN 1083 1086 interpolation_simple_corrector = .TRUE. 1084 1087 interpolation_simple_predictor = .FALSE. 1085 1088 interpolation_trilinear = .FALSE. 1086 ELSEIF ( TRIM( particle_advection_interpolation) == 'simple_predictor' ) THEN1089 ELSEIF ( TRIM( particle_advection_interpolation ) == 'simple_predictor' ) THEN 1087 1090 interpolation_simple_corrector = .FALSE. 1088 1091 interpolation_simple_predictor = .TRUE. … … 1334 1337 ! 1335 1338 !-- Calculate initial_weighting_factor diagnostically 1336 IF ( number_concentration /= -1.0_wp .AND.number_concentration > 0.0_wp ) THEN1337 initial_weighting_factor = number_concentration * &1339 IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN 1340 initial_weighting_factor = number_concentration * & 1338 1341 pdx(1) * pdy(1) * pdz(1) 1339 1342 END IF … … 1471 1474 alloc_factor / 100.0_wp ) ), 1 ) 1472 1475 IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) ) THEN 1473 CALL realloc_particles_array( ip,jp,kp,alloc_size)1476 CALL realloc_particles_array( ip, jp, kp, alloc_size ) 1474 1477 ENDIF 1475 1478 ENDIF … … 1510 1513 !-- Initialize aerosol background spectrum 1511 1514 IF ( curvature_solution_effects ) THEN 1512 CALL lpm_init_aerosols( local_start)1515 CALL lpm_init_aerosols( local_start ) 1513 1516 ENDIF 1514 1517 ! … … 1629 1632 SUBROUTINE lpm_init_aerosols(local_start) 1630 1633 1631 REAL(wp) ::afactor !< curvature effects1632 REAL(wp) ::bfactor !< solute effects1633 REAL(wp) ::dlogr !< logarithmic width of radius bin1634 REAL(wp) ::e_a !< vapor pressure1635 REAL(wp) ::e_s !< saturation vapor pressure1636 REAL(wp) ::rmin = 0.005e-6_wp !< minimum aerosol radius1637 REAL(wp) ::rmax = 10.0e-6_wp !< maximum aerosol radius1638 REAL(wp) ::r_mid !< mean radius of bin1639 REAL(wp) ::r_l !< left radius of bin1640 REAL(wp) ::r_r !< right radius of bin1641 REAL(wp) ::sigma !< surface tension1642 REAL(wp) ::t_int !< temperature1634 REAL(wp) :: afactor !< curvature effects 1635 REAL(wp) :: bfactor !< solute effects 1636 REAL(wp) :: dlogr !< logarithmic width of radius bin 1637 REAL(wp) :: e_a !< vapor pressure 1638 REAL(wp) :: e_s !< saturation vapor pressure 1639 REAL(wp) :: rmin = 0.005e-6_wp !< minimum aerosol radius 1640 REAL(wp) :: rmax = 10.0e-6_wp !< maximum aerosol radius 1641 REAL(wp) :: r_mid !< mean radius of bin 1642 REAL(wp) :: r_l !< left radius of bin 1643 REAL(wp) :: r_r !< right radius of bin 1644 REAL(wp) :: sigma !< surface tension 1645 REAL(wp) :: t_int !< temperature 1643 1646 1644 1647 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: local_start !< 1645 1648 1646 INTEGER(iwp) ::n !<1647 INTEGER(iwp) ::ip !<1648 INTEGER(iwp) ::jp !<1649 INTEGER(iwp) ::kp !<1649 INTEGER(iwp) :: n !< 1650 INTEGER(iwp) :: ip !< 1651 INTEGER(iwp) :: jp !< 1652 INTEGER(iwp) :: kp !< 1650 1653 1651 1654 ! 1652 1655 !-- Set constants for different aerosol species 1653 IF ( TRIM( aero_species) == 'nacl' ) THEN1656 IF ( TRIM( aero_species ) == 'nacl' ) THEN 1654 1657 molecular_weight_of_solute = 0.05844_wp 1655 1658 rho_s = 2165.0_wp 1656 1659 vanthoff = 2.0_wp 1657 ELSEIF ( TRIM( aero_species) == 'c3h4o4' ) THEN1660 ELSEIF ( TRIM( aero_species ) == 'c3h4o4' ) THEN 1658 1661 molecular_weight_of_solute = 0.10406_wp 1659 1662 rho_s = 1600.0_wp 1660 1663 vanthoff = 1.37_wp 1661 ELSEIF ( TRIM( aero_species) == 'nh4o3' ) THEN1664 ELSEIF ( TRIM( aero_species ) == 'nh4o3' ) THEN 1662 1665 molecular_weight_of_solute = 0.08004_wp 1663 1666 rho_s = 1720.0_wp … … 1671 1674 !-- The following typical aerosol spectra are taken from Jaenicke (1993): 1672 1675 !-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions. 1673 IF ( TRIM( aero_type) == 'polar' ) THEN1676 IF ( TRIM( aero_type ) == 'polar' ) THEN 1674 1677 na = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6_wp 1675 1678 rm = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6_wp 1676 1679 log_sigma = (/ 0.245, 0.300, 0.291 /) 1677 ELSEIF ( TRIM( aero_type) == 'background' ) THEN1680 ELSEIF ( TRIM( aero_type ) == 'background' ) THEN 1678 1681 na = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6_wp 1679 1682 rm = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6_wp 1680 1683 log_sigma = (/ 0.645, 0.253, 0.425 /) 1681 ELSEIF ( TRIM( aero_type) == 'maritime' ) THEN1684 ELSEIF ( TRIM( aero_type ) == 'maritime' ) THEN 1682 1685 na = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6_wp 1683 1686 rm = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6_wp 1684 1687 log_sigma = (/ 0.657, 0.210, 0.396 /) 1685 ELSEIF ( TRIM( aero_type) == 'continental' ) THEN1688 ELSEIF ( TRIM( aero_type ) == 'continental' ) THEN 1686 1689 na = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6_wp 1687 1690 rm = (/ 0.01, 0.058, 0.9 /) * 1.0E-6_wp 1688 1691 log_sigma = (/ 0.161, 0.217, 0.380 /) 1689 ELSEIF ( TRIM( aero_type) == 'desert' ) THEN1692 ELSEIF ( TRIM( aero_type ) == 'desert' ) THEN 1690 1693 na = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6_wp 1691 1694 rm = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6_wp 1692 1695 log_sigma = (/ 0.247, 0.770, 0.438 /) 1693 ELSEIF ( TRIM( aero_type) == 'rural' ) THEN1696 ELSEIF ( TRIM( aero_type ) == 'rural' ) THEN 1694 1697 na = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6_wp 1695 1698 rm = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6_wp 1696 1699 log_sigma = (/ 0.225, 0.557, 0.266 /) 1697 ELSEIF ( TRIM( aero_type) == 'urban' ) THEN1700 ELSEIF ( TRIM( aero_type ) == 'urban' ) THEN 1698 1701 na = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6_wp 1699 1702 rm = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6_wp 1700 1703 log_sigma = (/ 0.245, 0.666, 0.337 /) 1701 ELSEIF ( TRIM( aero_type) == 'user' ) THEN1704 ELSEIF ( TRIM( aero_type ) == 'user' ) THEN 1702 1705 CONTINUE 1703 1706 ELSE … … 2117 2120 ! 2118 2121 !-- If necessary, release new set of particles 2119 IF ( ( simulated_time - last_particle_release_time ) >= dt_prel .AND. end_time_prel > simulated_time )&2120 THEN2122 IF ( ( simulated_time - last_particle_release_time ) >= dt_prel .AND. & 2123 end_time_prel > simulated_time ) THEN 2121 2124 DO WHILE ( ( simulated_time - last_particle_release_time ) >= dt_prel ) 2122 2125 CALL lpm_create_particle( PHASE_RELEASE ) … … 2212 2215 ! 2213 2216 !-- Particle advection 2214 CALL lpm_advec( i,j,k)2217 CALL lpm_advec( i, j, k ) 2215 2218 ! 2216 2219 !-- Particle reflection from walls. Only applied if the particles 2217 2220 !-- are in the vertical range of the topography. (Here, some 2218 2221 !-- optimization is still possible.) 2219 IF ( topography /= 'flat' .AND.k < nzb_max + 2 ) THEN2222 IF ( topography /= 'flat' .AND. k < nzb_max + 2 ) THEN 2220 2223 CALL lpm_boundary_conds( 'walls', i, j, k ) 2221 2224 ENDIF … … 2223 2226 !-- User-defined actions after the calculation of the new particle 2224 2227 !-- position 2225 CALL user_lpm_advec( i,j,k)2228 CALL user_lpm_advec( i, j, k ) 2226 2229 ! 2227 2230 !-- Apply boundary conditions to those particles that have crossed … … 2284 2287 ! 2285 2288 !-- IF .FALSE., lpm_sort_and_delete is done inside pcmp 2286 IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run ) THEN2289 IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run ) THEN 2287 2290 ! 2288 2291 !-- Pack particles (eliminate those marked for deletion), … … 2386 2389 SUBROUTINE lpm_write_exchange_statistics 2387 2390 2388 INTEGER(iwp) :: ip !<2389 INTEGER(iwp) :: jp !<2390 INTEGER(iwp) :: kp !<2391 INTEGER(iwp) :: tot_number_of_particles2391 INTEGER(iwp) :: ip !< 2392 INTEGER(iwp) :: jp !< 2393 INTEGER(iwp) :: kp !< 2394 INTEGER(iwp) :: tot_number_of_particles !< 2392 2395 2393 2396 ! … … 2462 2465 WRITE ( 85 ) simulated_time 2463 2466 WRITE ( 85 ) prt_count 2464 2467 2465 2468 DO ip = nxl, nxr 2466 2469 DO jp = nys, nyn … … 2898 2901 CHARACTER (LEN=10) :: version_on_file !< 2899 2902 2900 INTEGER(iwp) :: alloc_size !<2901 INTEGER(iwp) :: ip !<2902 INTEGER(iwp) :: jp !<2903 INTEGER(iwp) :: kp !<2904 2905 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles !<2903 INTEGER(iwp) :: alloc_size !< 2904 INTEGER(iwp) :: ip !< 2905 INTEGER(iwp) :: jp !< 2906 INTEGER(iwp) :: kp !< 2907 2908 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles !< 2906 2909 2907 2910 ! … … 3016 3019 LOGICAL, INTENT(OUT) :: found 3017 3020 3018 REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<3021 REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< 3019 3022 3020 3023 … … 3359 3362 !-- Predictor step 3360 3363 kkw = kp - 1 3361 DO n = 1, number_of_particles3364 DO n = 1, number_of_particles 3362 3365 3363 3366 alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) … … 3374 3377 ! 3375 3378 !-- Corrector step 3376 DO n = 1, number_of_particles3379 DO n = 1, number_of_particles 3377 3380 3378 3381 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3379 3382 3380 DO nn = 1, iteration_steps3383 DO nn = 1, iteration_steps 3381 3384 3382 3385 ! … … 3428 3431 !-- The particle position for the w velociy is based on the value of kp and kp-1 3429 3432 kkw = kp - 1 3430 DO n = 1, number_of_particles3433 DO n = 1, number_of_particles 3431 3434 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3432 3435 … … 3830 3833 ENDDO 3831 3834 3832 DO nb = 0,73835 DO nb = 0,7 3833 3836 i = ip + block_offset(nb)%i_off 3834 3837 j = jp + block_offset(nb)%j_off … … 3963 3966 ENDIF 3964 3967 3965 CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n),&3968 CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n), & 3966 3969 de_dx_int(n), de_dt, diss_int(n), & 3967 3970 dt_particle(n), rg(n,1), term_1_2(n) ) 3968 3971 3969 CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n),&3972 CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n), & 3970 3973 de_dy_int(n), de_dt, diss_int(n), & 3971 3974 dt_particle(n), rg(n,2), term_1_2(n) ) 3972 3975 3973 CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n),&3976 CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n), & 3974 3977 de_dz_int(n), de_dt, diss_int(n), & 3975 3978 dt_particle(n), rg(n,3), term_1_2(n) ) … … 4011 4014 ENDIF 4012 4015 4013 CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n),&4016 CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n), & 4014 4017 de_dx_int(n), de_dt, diss_int(n), & 4015 4018 dt_particle(n), rg(n,1), term_1_2(n) ) 4016 4019 4017 CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n),&4020 CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n), & 4018 4021 de_dy_int(n), de_dt, diss_int(n), & 4019 4022 dt_particle(n), rg(n,2), term_1_2(n) ) 4020 4023 4021 CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n),&4024 CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n), & 4022 4025 de_dz_int(n), de_dt, diss_int(n), & 4023 4026 dt_particle(n), rg(n,3), term_1_2(n) ) … … 4420 4423 INTEGER(iwp) :: tmp_z !< dummy for sorting algorithm 4421 4424 4422 INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions4423 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions4424 INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions4425 INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 4426 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions 4427 INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions 4425 4428 4426 4429 LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary … … 4556 4559 i2 = particles(n)%x * ddx 4557 4560 j2 = particles(n)%y * ddy 4558 IF ( zw(k) < particles(n)%z ) k2 = k + 14559 IF ( zw(k) > particles(n)%z .AND.zw(k-1) < particles(n)%z ) k2 = k4560 IF ( zw(k-1) > particles(n)%z ) k2 = k - 14561 IF ( zw(k) < particles(n)%z ) k2 = k + 1 4562 IF ( zw(k) > particles(n)%z .AND. zw(k-1) < particles(n)%z ) k2 = k 4563 IF ( zw(k-1) > particles(n)%z ) k2 = k - 1 4561 4564 ! 4562 4565 !-- Save current particle positions … … 4590 4593 !-- Wall to the north 4591 4594 IF ( prt_y > pos_y_old ) THEN 4592 ywall = ( j1 + 1 ) * dy4595 ywall = ( j1 + 1 ) * dy 4593 4596 !-- Wall to the south 4594 4597 ELSE … … 4645 4648 !-- Store these values only if particle really reaches any wall. t must 4646 4649 !-- be in a interval between [0:1]. 4647 IF ( t(t_index) <= 1.0_wp .AND.t(t_index) >= 0.0_wp ) THEN4650 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN 4648 4651 t_index = t_index + 1 4649 4652 cross_wall_x = .TRUE. … … 4661 4664 reach_y(t_index) = .TRUE. 4662 4665 reach_z(t_index) = .FALSE. 4663 IF ( t(t_index) <= 1.0_wp .AND.t(t_index) >= 0.0_wp ) THEN4666 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN 4664 4667 t_index = t_index + 1 4665 4668 cross_wall_y = .TRUE. … … 4678 4681 reach_y(t_index) = .FALSE. 4679 4682 reach_z(t_index) = .TRUE. 4680 IF( t(t_index) <= 1.0_wp .AND.t(t_index) >= 0.0_wp) THEN4683 IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN 4681 4684 t_index = t_index + 1 4682 4685 cross_wall_z = .TRUE. … … 4686 4689 ! 4687 4690 !-- Carry out reflection only if particle reaches any wall 4688 IF ( cross_wall_x .OR. cross_wall_y .OR.cross_wall_z ) THEN4691 IF ( cross_wall_x .OR. cross_wall_y .OR. cross_wall_z ) THEN 4689 4692 ! 4690 4693 !-- Sort fractional timesteps in ascending order. Also sort the … … 4885 4888 !-- If a particle was reflected, calculate final position from last 4886 4889 !-- intermediate position. 4887 IF ( reflect_x .OR. reflect_y .OR.reflect_z ) THEN4890 IF ( reflect_x .OR. reflect_y .OR. reflect_z ) THEN 4888 4891 4889 4892 particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle & … … 4922 4925 SUBROUTINE lpm_droplet_condensation (i,j,k) 4923 4926 4924 INTEGER(iwp), INTENT(IN) :: i !<4925 INTEGER(iwp), INTENT(IN) :: j !<4926 INTEGER(iwp), INTENT(IN) :: k !<4927 INTEGER(iwp) :: n !<4927 INTEGER(iwp), INTENT(IN) :: i !< 4928 INTEGER(iwp), INTENT(IN) :: j !< 4929 INTEGER(iwp), INTENT(IN) :: k !< 4930 INTEGER(iwp) :: n !< 4928 4931 4929 4932 REAL(wp) :: afactor !< curvature effects … … 4954 4957 ! 4955 4958 !-- Parameters for Rosenbrock method (see Verwer et al., 1999) 4956 REAL(wp), PARAMETER :: prec = 1.0E-3_wp !< precision of Rosenbrock solution4957 REAL(wp), PARAMETER :: q_increase = 1.5_wp !< increase factor in timestep4958 REAL(wp), PARAMETER :: q_decrease = 0.9_wp !< decrease factor in timestep4959 REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)4959 REAL(wp), PARAMETER :: prec = 1.0E-3_wp !< precision of Rosenbrock solution 4960 REAL(wp), PARAMETER :: q_increase = 1.5_wp !< increase factor in timestep 4961 REAL(wp), PARAMETER :: q_decrease = 0.9_wp !< decrease factor in timestep 4962 REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0) 4960 4963 ! 4961 4964 !-- Parameters for terminal velocity … … 5152 5155 !-- case, the model timestep might be too long. 5153 5156 delta_r = new_r(n) - particles(n)%radius 5154 IF ( delta_r < 0.0_wp .AND. new_r(n) < 0.0_wp ) THEN5155 WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, &5157 IF ( delta_r < 0.0_wp .AND. new_r(n) < 0.0_wp ) THEN 5158 WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, & 5156 5159 ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, & 5157 5160 '&delta_r=',delta_r, & … … 6158 6161 SUBROUTINE turb_enhance_eff 6159 6162 6160 INTEGER(iwp) :: i !<6161 INTEGER(iwp) :: iq !<6162 INTEGER(iwp) :: ir !<6163 INTEGER(iwp) :: j !<6164 INTEGER(iwp) :: k !<6165 INTEGER(iwp) :: kk !<6163 INTEGER(iwp) :: i !< 6164 INTEGER(iwp) :: iq !< 6165 INTEGER(iwp) :: ir !< 6166 INTEGER(iwp) :: j !< 6167 INTEGER(iwp) :: k !< 6168 INTEGER(iwp) :: kk !< 6166 6169 6167 6170 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira !< … … 6348 6351 INTEGER(iwp) :: np !< 6349 6352 INTEGER(iwp) :: old_size !< old particle array size 6350 6353 6351 6354 INTEGER(iwp), PARAMETER :: n_max = 100 !< number of radii bin for splitting functions 6352 6355 6353 6356 LOGICAL :: first_loop_stride_sp = .TRUE. !< flag to calculate constants only once 6354 6357 … … 6384 6387 REAL(wp), DIMENSION(0:n_max-1) :: r_bin_mid !< mass weighted mean radius of a bin 6385 6388 REAL(wp), DIMENSION(0:n_max) :: r_bin !< boundaries of a radius bin 6386 6389 6387 6390 TYPE(particle_type) :: tmp_particle !< temporary particle TYPE 6388 6391 … … 6442 6445 !-- Reallocate particle array if necessary. 6443 6446 IF ( new_size > SIZE(particles) ) THEN 6444 CALL realloc_particles_array( i,j,k,new_size)6447 CALL realloc_particles_array( i, j, k, new_size ) 6445 6448 ENDIF 6446 6449 old_size = prt_count(k,j,i) … … 6652 6655 !-- Reallocate particle array if necessary. 6653 6656 IF ( new_size > SIZE(particles) ) THEN 6654 CALL realloc_particles_array( i,j,k,new_size)6657 CALL realloc_particles_array( i, j, k, new_size ) 6655 6658 ENDIF 6656 6659 old_size = prt_count(k,j,i) … … 6819 6822 !-- Reallocate particle array if necessary. 6820 6823 IF ( new_size > SIZE(particles) ) THEN 6821 CALL realloc_particles_array( i,j,k,new_size)6824 CALL realloc_particles_array( i, j, k, new_size ) 6822 6825 ENDIF 6823 6826 ! … … 6829 6832 ! 6830 6833 !-- Create splitting_factor-1 new particles. 6831 DO jpp = 1, splitting_factor-16832 grid_particles(k,j,i)%particles( jpp+old_size) =&6833 tmp_particle 6834 DO jpp = 1, splitting_factor-1 6835 grid_particles(k,j,i)%particles( jpp + old_size ) = & 6836 tmp_particle 6834 6837 ENDDO 6835 6838 ! … … 6840 6843 splitting_factor - 1 6841 6844 ENDIF 6842 ENDDO 6845 ENDDO 6843 6846 ENDDO 6844 6847 ENDDO … … 6846 6849 ENDDO 6847 6850 ENDIF 6848 6851 6849 6852 CALL cpu_log( log_point_s(80), 'lpm_splitting', 'stop' ) 6850 6853 … … 7610 7613 pack_done = .FALSE. 7611 7614 7612 DO n = 1, SIZE(particle_array)7615 DO n = 1, SIZE(particle_array) 7613 7616 7614 7617 IF ( .NOT. particle_array(n)%particle_mask ) CYCLE … … 7633 7636 IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) ) THEN 7634 7637 IF ( pack_done ) THEN 7635 CALL realloc_particles_array ( ip,jp,kp)7638 CALL realloc_particles_array ( ip, jp, kp ) 7636 7639 ELSE 7637 7640 CALL lpm_pack … … 7639 7642 pindex = prt_count(kp,jp,ip)+1 7640 7643 IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) ) THEN 7641 CALL realloc_particles_array ( ip,jp,kp)7644 CALL realloc_particles_array ( ip, jp, kp ) 7642 7645 ENDIF 7643 7646 pack_done = .TRUE. … … 7813 7816 IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) & 7814 7817 THEN 7815 CALL realloc_particles_array( i,j,k)7818 CALL realloc_particles_array( i, j, k ) 7816 7819 ENDIF 7817 7820 … … 7855 7858 IF ( number_of_particles <= 0 ) CYCLE 7856 7859 particles => grid_particles(k,j,i)%particles(1:number_of_particles) 7857 DO n = 1, number_of_particles7860 DO n = 1, number_of_particles 7858 7861 ! 7859 7862 !-- Note, check for CFL does not work at first particle timestep … … 7950 7953 7951 7954 7952 INTEGER(iwp) :: i 7953 INTEGER(iwp) :: j 7954 INTEGER(iwp) :: k 7955 INTEGER(iwp) :: old_size !<7956 INTEGER(iwp) :: new_size !<7957 7958 LOGICAL :: dealloc7959 7960 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<7961 TYPE(particle_type), DIMENSION(500) :: tmp_particles_s !<7955 INTEGER(iwp) :: i !< 7956 INTEGER(iwp) :: j !< 7957 INTEGER(iwp) :: k !< 7958 INTEGER(iwp) :: old_size !< 7959 INTEGER(iwp) :: new_size !< 7960 7961 LOGICAL :: dealloc 7962 7963 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !< 7964 TYPE(particle_type), DIMENSION(500) :: tmp_particles_s !< 7962 7965 7963 7966 DO i = nxl, nxr … … 8082 8085 !-- the number_of_particles to the actual value. 8083 8086 nn = 0 8084 DO is = 0,78087 DO is = 0,7 8085 8088 grid_particles(kp,jp,ip)%start_index(is) = nn + 1 8086 DO n = 1,sort_count(is)8089 DO n = 1,sort_count(is) 8087 8090 nn = nn + 1 8088 8091 particles(nn) = sort_particles(n,is) … … 8175 8178 SUBROUTINE lpm_sort_timeloop_done 8176 8179 8177 INTEGER(iwp) :: end_index !< particle end index for each sub-box8178 INTEGER(iwp) :: i !< index of particle grid box in x-direction8179 INTEGER(iwp) :: j !< index of particle grid box in y-direction8180 INTEGER(iwp) :: k !< index of particle grid box in z-direction8181 INTEGER(iwp) :: n !< running index for number of particles8182 INTEGER(iwp) :: nb !< index of subgrid boux8183 INTEGER(iwp) :: nf !< indices for particles in each sub-box that already finalized their substeps8184 INTEGER(iwp) :: nnf !< indices for particles in each sub-box that need further treatment8185 INTEGER(iwp) :: num_finalized !< number of particles in each sub-box that already finalized their substeps8186 INTEGER(iwp) :: start_index !< particle start index for each sub-box8187 8188 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: sort_particles !< temporary particle array8180 INTEGER(iwp) :: end_index !< particle end index for each sub-box 8181 INTEGER(iwp) :: i !< index of particle grid box in x-direction 8182 INTEGER(iwp) :: j !< index of particle grid box in y-direction 8183 INTEGER(iwp) :: k !< index of particle grid box in z-direction 8184 INTEGER(iwp) :: n !< running index for number of particles 8185 INTEGER(iwp) :: nb !< index of subgrid boux 8186 INTEGER(iwp) :: nf !< indices for particles in each sub-box that already finalized their substeps 8187 INTEGER(iwp) :: nnf !< indices for particles in each sub-box that need further treatment 8188 INTEGER(iwp) :: num_finalized !< number of particles in each sub-box that already finalized their substeps 8189 INTEGER(iwp) :: start_index !< particle start index for each sub-box 8190 8191 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: sort_particles !< temporary particle array 8189 8192 8190 8193 DO i = nxl, nxr
Note: See TracChangeset
for help on using the changeset viewer.
