- Timestamp:
- Apr 15, 2016 11:46:09 AM (9 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 5 edited
-
Makefile (modified) (4 diffs)
-
lpm_droplet_condensation.f90 (modified) (15 diffs)
-
lpm_init_mod.f90 (modified) (5 diffs)
-
mod_particle_attributes.f90 (modified) (3 diffs)
-
package_parin.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r1852 r1871 20 20 # Current revisions: 21 21 # ------------------ 22 # 22 # dependencies for package_parin corrected 23 23 # 24 24 # Former revisions: … … 399 399 data_output_3d.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o\ 400 400 netcdf_interface_mod.o land_surface_model_mod.o 401 diffusion_e_mod.o: modules.o mod_kinds.o microphysics_mod.o 401 diffusion_e_mod.o: modules.o mod_kinds.o microphysics_mod.o \ 402 mod_particle_attributes.o 402 403 diffusion_s_mod.o: modules.o mod_kinds.o 403 404 diffusion_u_mod.o: modules.o mod_kinds.o wall_fluxes_mod.o … … 423 424 lpm_init_mod.o ls_forcing_mod.o netcdf_interface_mod.o plant_canopy_model_mod.o \ 424 425 radiation_model_mod.o random_function_mod.o random_generator_parallel_mod.o \ 425 surface_layer_fluxes_mod.o microphysics_mod.o 426 surface_layer_fluxes_mod.o microphysics_mod.o mod_particle_attributes.o 426 427 init_advec.o: modules.o mod_kinds.o 427 428 init_cloud_physics.o: modules.o mod_kinds.o … … 475 476 land_surface_model_mod.o spectra_mod.o 476 477 nudging_mod.o: modules.o cpulog_mod.o mod_kinds.o 477 package_parin.o: modules.o mod_kinds.o 478 package_parin.o: modules.o mod_kinds.o mod_particle_attributes.o 478 479 palm.o: modules.o cpulog_mod.o ls_forcing_mod.o mod_kinds.o nudging_mod.o\ 479 480 pmc_interface_mod.o surface_layer_fluxes_mod.o -
palm/trunk/SOURCE/lpm_droplet_condensation.f90
r1852 r1871 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Initialization of aerosols added. 22 22 ! 23 23 ! Former revisions: 24 24 ! ----------------- 25 25 ! $Id$ 26 !27 26 ! 28 27 ! 1849 2016-04-08 11:33:18Z hoffmann … … 121 120 122 121 USE particle_attributes, & 123 ONLY: curvature_solution_effects, hall_kernel, mass_of_solute,&122 ONLY: curvature_solution_effects, hall_kernel, & 124 123 molecular_weight_of_solute, molecular_weight_of_water, & 125 number_of_particles, particles, radius_classes, 124 number_of_particles, particles, radius_classes, rho_s, & 126 125 use_kernel_tables, vanthoff, wang_kernel 127 126 … … 140 139 INTEGER(iwp) :: ros_count !< 141 140 142 INTEGER(iwp), PARAMETER :: maxtry = 40!<141 INTEGER(iwp), PARAMETER :: maxtry = 100 !< 143 142 144 143 LOGICAL :: repeat !< … … 195 194 REAL(wp), PARAMETER :: e3 = 0.0_wp !< 196 195 REAL(wp), PARAMETER :: e4 = 125.0_wp / 108.0_wp !< 196 REAL(wp), PARAMETER :: eps_ros = 1.0E-3_wp !< accuracy of Rosenbrock method 197 197 REAL(wp), PARAMETER :: gam = 0.5_wp !< 198 198 REAL(wp), PARAMETER :: grow = 1.5_wp !< … … 200 200 REAL(wp), PARAMETER :: pshrnk = -1.0_wp /3.0_wp !< 201 201 REAL(wp), PARAMETER :: shrnk = 0.5_wp !< 202 REAL(wp), PARAMETER :: eps_ros = 1.0E-4_wp !< accuracy of Rosenbrock method203 202 204 203 ! … … 311 310 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int ) 312 311 ! 313 !-- Solute effect (bfactor), mass of solute to be replaced by variable 314 !-- aerosol radius 315 bfactor = 3.0_wp * vanthoff * mass_of_solute * & 316 molecular_weight_of_water / ( 4.0_wp * pi * rho_l * & 317 molecular_weight_of_solute & 318 ) 312 !-- Solute effect (bfactor) 313 bfactor = vanthoff * rho_s * particles(n)%rvar2**3 * & 314 molecular_weight_of_water / & 315 ( rho_l * molecular_weight_of_solute ) 319 316 320 317 r_ros = particles(n)%radius … … 326 323 327 324 ! 328 !-- Internal time step should not be > 1.0E-2 in case of evaporation 329 !-- because larger values may lead to secondary solutions which are 330 !-- physically unlikely 331 IF ( dt_ros_next > 1.0E-2_wp .AND. e_a / e_s < 1.0_wp ) THEN 325 !-- Internal time step should not be > 1.0E-2 and < 1.0E-6 326 IF ( dt_ros_next > 1.0E-2_wp ) THEN 332 327 dt_ros_next = 1.0E-3_wp 328 ELSEIF ( dt_ros_next < 1.0E-6_wp ) THEN 329 dt_ros_next = 1.0E-6_wp 333 330 ENDIF 331 334 332 ! 335 333 !-- If calculation of Rosenbrock method is repeated due to unreasonalble … … 337 335 !-- reduced 338 336 IF ( ros_count > 1 ) THEN 339 dt_ros_next = dt_ros_next - ( 0.2_wp * dt_ros_next )337 dt_ros_next = dt_ros_next * 0.1_wp 340 338 ELSEIF ( ros_count > 5 ) THEN 341 339 ! … … 349 347 !-- Internal time step must not be larger than PALM time step 350 348 dt_ros = MIN( dt_ros_next, dt_3d ) 349 351 350 ! 352 351 !-- Integrate growth equation in time unless PALM time step is reached … … 364 363 drdt_ini = drdt 365 364 dt_ros_sum_ini = dt_ros_sum 366 r_ros_ini = r_ros365 r_ros_ini = MAX( r_ros, particles(n)%rvar2 ) 367 366 368 367 ! … … 378 377 379 378 IF ( jtry == maxtry+1 ) THEN 380 message_string = 'maxtry > 40 in Rosenbrock method'379 message_string = 'maxtry > 100 in Rosenbrock method' 381 380 CALL message( 'lpm_droplet_condensation', 'PA0347', 2, & 382 381 2, 0, 6, 0 ) … … 385 384 aa = 1.0_wp / ( gam * dt_ros ) - d2rdtdr 386 385 g1 = drdt_ini / aa 387 r_ros = r_ros_ini + a21 * g1386 r_ros = MAX( r_ros_ini + a21 * g1, particles(n)%rvar2 ) 388 387 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 389 388 afactor / r_ros + & … … 393 392 g2 = ( drdt + c21 * g1 / dt_ros )& 394 393 / aa 395 r_ros = r_ros_ini + a31 * g1 + a32 * g2394 r_ros = MAX( r_ros_ini + a31 * g1 + a32 * g2, particles(n)%rvar2 ) 396 395 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - & 397 396 afactor / r_ros + & … … 403 402 g4 = ( drdt + & 404 403 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 405 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 404 r_ros = MAX( r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + & 405 b4 * g4, particles(n)%rvar2 ) 406 406 407 407 dt_ros_sum = dt_ros_sum_ini + dt_ros … … 451 451 particles(n)%rvar1 = dt_ros_next 452 452 453 new_r(n) = r_ros454 453 ! 455 454 !-- Radius should not fall below 1E-8 because Rosenbrock method may 456 455 !-- lead to errors otherwise 457 new_r(n) = MAX( new_r(n), 1.0E-8_wp)456 new_r(n) = MAX( r_ros, particles(n)%rvar2 ) 458 457 ! 459 458 !-- Check if calculated droplet radius change is reasonable since in -
palm/trunk/SOURCE/lpm_init_mod.f90
r1851 r1871 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Initialization of aerosols added. 22 22 ! 23 23 ! Former revisions: … … 509 509 ONLY: lpm_pack_all_arrays 510 510 511 USE particle_attributes, & 512 ONLY: monodisperse_aerosols 513 511 514 IMPLICIT NONE 512 515 … … 581 584 !-- Initial values (internal timesteps, derivative) 582 585 !-- for Rosenbrock method 583 tmp_particle%rvar1 = 1.0E- 12_wp584 tmp_particle%rvar2 = 1.0E-3_wp585 tmp_particle%rvar3 = -9999999.9_wp 586 tmp_particle%rvar1 = 1.0E-6_wp !last Rosenbrock timestep 587 tmp_particle%rvar2 = 0.1E-6_wp !dry aerosol radius 588 tmp_particle%rvar3 = -9999999.9_wp !unused 586 589 ELSE 587 590 ! … … 686 689 local_start = prt_count+1 687 690 prt_count = local_count 691 692 ! 693 !-- Initialize aerosol background spectrum 694 IF ( curvature_solution_effects .AND. .NOT. monodisperse_aerosols ) THEN 695 CALL lpm_init_aerosols(local_start) 696 ENDIF 697 688 698 ! 689 699 !-- Add random fluctuation to particle positions … … 759 769 END SUBROUTINE lpm_create_particle 760 770 771 SUBROUTINE lpm_init_aerosols(local_start) 772 773 USE arrays_3d, & 774 ONLY: hyp, pt, q 775 776 USE cloud_parameters, & 777 ONLY: l_d_rv, rho_l 778 779 USE constants, & 780 ONLY: pi 781 782 USE kinds 783 784 USE particle_attributes, & 785 ONLY: init_aerosol_probabilistic, molecular_weight_of_solute, & 786 molecular_weight_of_water, n1, n2, n3, rho_s, rm1, rm2, rm3, & 787 s1, s2, s3, vanthoff 788 789 IMPLICIT NONE 790 791 REAL(wp), DIMENSION(:), ALLOCATABLE :: cdf !< CDF of aerosol spectrum 792 REAL(wp), DIMENSION(:), ALLOCATABLE :: r_temp !< dry aerosol radius spectrum 793 794 REAL(wp) :: bfactor !< solute effects 795 REAL(wp) :: dr !< width of radius bin 796 REAL(wp) :: e_a !< vapor pressure 797 REAL(wp) :: e_s !< saturation vapor pressure 798 REAL(wp) :: n_init !< sum of all aerosol concentrations 799 REAL(wp) :: pdf !< PDF of aerosol spectrum 800 REAL(wp) :: rmin = 1.0e-8_wp !< minimum aerosol radius 801 REAL(wp) :: rmax = 1.0e-6_wp !< maximum aerosol radius 802 REAL(wp) :: rs_rand !< random number 803 REAL(wp) :: r_mid !< mean radius 804 REAL(wp) :: t_int !< temperature 805 REAL(wp) :: weight_sum !< sum of all weighting factors 806 807 INTEGER(iwp), DIMENSION(:,:,:), INTENT(IN) :: local_start !< 808 809 INTEGER(iwp) :: n !< 810 INTEGER(iwp) :: nn !< 811 INTEGER(iwp) :: no_bins = 999 !< number of bins 812 INTEGER(iwp) :: ip !< 813 INTEGER(iwp) :: jp !< 814 INTEGER(iwp) :: kp !< 815 816 LOGICAL :: new_pdf = .FALSE. !< check if aerosol PDF has to be recalculated 817 818 ! 819 !-- Compute aerosol background distribution 820 IF ( init_aerosol_probabilistic ) THEN 821 ALLOCATE( cdf(0:no_bins), r_temp(0:no_bins) ) 822 DO n = 0, no_bins 823 r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / & 824 REAL(no_bins, KIND=wp) * REAL(n, KIND=wp) ) 825 826 cdf(n) = 0.0_wp 827 n_init = n1 + n2 + n3 828 IF ( n1 > 0.0_wp ) THEN 829 cdf(n) = cdf(n) + n1 / n_init * ( 0.5_wp + 0.5_wp * & 830 ERF( LOG( r_temp(n) / rm1 ) / & 831 ( SQRT(2.0_wp) * LOG(s1) ) & 832 ) ) 833 ENDIF 834 IF ( n2 > 0.0_wp ) THEN 835 cdf(n) = cdf(n) + n2 / n_init * ( 0.5_wp + 0.5_wp * & 836 ERF( LOG( r_temp(n) / rm2 ) / & 837 ( SQRT(2.0_wp) * LOG(s2) ) & 838 ) ) 839 ENDIF 840 IF ( n3 > 0.0_wp ) THEN 841 cdf(n) = cdf(n) + n3 / n_init * ( 0.5_wp + 0.5_wp * & 842 ERF( LOG( r_temp(n) / rm3 ) / & 843 ( SQRT(2.0_wp) * LOG(s3) ) & 844 ) ) 845 ENDIF 846 847 ENDDO 848 ENDIF 849 850 DO ip = nxl, nxr 851 DO jp = nys, nyn 852 DO kp = nzb+1, nzt 853 854 number_of_particles = prt_count(kp,jp,ip) 855 IF ( number_of_particles <= 0 ) CYCLE 856 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 857 ! 858 !-- Initialize the aerosols with a predefined spectral distribution 859 !-- of the dry radius (logarithmically increasing bins) and a varying 860 !-- weighting factor 861 IF ( .NOT. init_aerosol_probabilistic ) THEN 862 863 new_pdf = .FALSE. 864 IF ( .NOT. ALLOCATED( r_temp ) ) THEN 865 new_pdf = .TRUE. 866 ELSE 867 IF ( SIZE( r_temp ) .NE. & 868 number_of_particles - local_start(kp,jp,ip) + 2 ) THEN 869 new_pdf = .TRUE. 870 DEALLOCATE( r_temp ) 871 ENDIF 872 ENDIF 873 874 IF ( new_pdf ) THEN 875 876 no_bins = number_of_particles + 1 - local_start(kp,jp,ip) 877 ALLOCATE( r_temp(0:no_bins) ) 878 879 DO n = 0, no_bins 880 r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / & 881 REAL(no_bins, KIND=wp) * & 882 REAL(n, KIND=wp) ) 883 ENDDO 884 885 ENDIF 886 887 ! 888 !-- Calculate radius and concentration of each aerosol 889 DO n = local_start(kp,jp,ip), number_of_particles 890 891 nn = n - local_start(kp,jp,ip) 892 893 r_mid = SQRT( r_temp(nn) * r_temp(nn+1) ) 894 dr = r_temp(nn+1) - r_temp(nn) 895 896 pdf = 0.0_wp 897 n_init = n1 + n2 + n3 898 IF ( n1 > 0.0_wp ) THEN 899 pdf = pdf + n1 / n_init * ( 1.0_wp / ( r_mid * LOG(s1) * & 900 SQRT( 2.0_wp * pi ) & 901 ) * & 902 EXP( -( LOG( r_mid / rm1 ) )**2 / & 903 ( 2.0_wp * LOG(s1)**2 ) & 904 ) & 905 ) 906 ENDIF 907 IF ( n2 > 0.0_wp ) THEN 908 pdf = pdf + n2 / n_init * ( 1.0_wp / ( r_mid * LOG(s2) * & 909 SQRT( 2.0_wp * pi ) & 910 ) * & 911 EXP( -( LOG( r_mid / rm2 ) )**2 / & 912 ( 2.0_wp * LOG(s2)**2 ) & 913 ) & 914 ) 915 ENDIF 916 IF ( n3 > 0.0_wp ) THEN 917 pdf = pdf + n3 / n_init * ( 1.0_wp / ( r_mid * LOG(s3) * & 918 SQRT( 2.0_wp * pi ) & 919 ) * & 920 EXP( -( LOG( r_mid / rm3 ) )**2 / & 921 ( 2.0_wp * LOG(s3)**2 ) & 922 ) & 923 ) 924 ENDIF 925 926 particles(n)%rvar2 = r_mid 927 particles(n)%weight_factor = pdf * dr 928 929 END DO 930 ! 931 !-- Adjust weighting factors to initialize the same number of aerosols 932 !-- in every grid box 933 weight_sum = SUM(particles(local_start(kp,jp,ip):number_of_particles)%weight_factor) 934 935 particles(local_start(kp,jp,ip):number_of_particles)%weight_factor = & 936 particles(local_start(kp,jp,ip):number_of_particles)%weight_factor / & 937 weight_sum * initial_weighting_factor * ( no_bins + 1 ) 938 939 ENDIF 940 ! 941 !-- Initialize the aerosols with a predefined weighting factor but 942 !-- a randomly choosen dry radius 943 IF ( init_aerosol_probabilistic ) THEN 944 945 DO n = local_start(kp,jp,ip), number_of_particles !only new particles 946 947 rs_rand = -1.0_wp 948 DO WHILE ( rs_rand .LT. cdf(0) .OR. rs_rand .GE. cdf(no_bins) ) 949 rs_rand = random_function( iran_part ) 950 ENDDO 951 ! 952 !-- Determine aerosol dry radius by a random number generator 953 DO nn = 0, no_bins-1 954 IF ( cdf(nn) .LE. rs_rand .AND. cdf(nn+1) .GT. rs_rand ) THEN 955 particles(n)%rvar2 = r_temp(nn) + ( r_temp(nn+1) - r_temp(nn) ) / & 956 ( cdf(nn+1) - cdf(nn) ) * ( rs_rand - cdf(nn) ) 957 EXIT 958 ENDIF 959 ENDDO 960 961 ENDDO 962 963 ENDIF 964 965 ! 966 !-- Set particle radius to equilibrium radius based on the environmental 967 !-- supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids 968 !-- the sometimes lengthy growth toward their equilibrium radius within 969 !-- the simulation. 970 t_int = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp 971 972 e_s = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) ) 973 e_a = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp ) 974 975 ! 976 !-- The formula is only valid for subsaturated environments. In (super-) 977 !-- saturated air, the inital radius is used. 978 IF ( e_a / e_s < 1.0_wp ) THEN 979 980 DO n = local_start(kp,jp,ip), number_of_particles !only new particles 981 982 bfactor = vanthoff * molecular_weight_of_water * & 983 rho_s * particles(n)%rvar2**3 / & 984 ( molecular_weight_of_solute * rho_l ) 985 particles(n)%radius = particles(n)%rvar2 * ( bfactor / & 986 particles(n)%rvar2**3 )**(1.0_wp/3.0_wp) *& 987 ( 1.0_wp - e_a / e_s )**(-1.0_wp/3.0_wp) 988 989 ENDDO 990 991 ENDIF 992 993 ENDDO 994 ENDDO 995 ENDDO 996 ! 997 !-- Deallocate used arrays 998 IF ( ALLOCATED(r_temp) ) DEALLOCATE( r_temp ) 999 IF ( ALLOCATED(cdf) ) DEALLOCATE( cdf ) 1000 1001 END SUBROUTINE lpm_init_aerosols 1002 761 1003 END MODULE lpm_init_mod -
palm/trunk/SOURCE/mod_particle_attributes.f90
r1852 r1871 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Initialization of aerosols added. 22 22 ! 23 23 ! Former revisions: … … 119 119 REAL(wp), DIMENSION(:), ALLOCATABLE :: log_z_z0 120 120 121 !122 !-- Lagrangian cloud model constants123 REAL(wp) :: mass_of_solute = 1.0E-17_wp !< soluted NaCl (kg)124 REAL(wp) :: molecular_weight_of_solute = 0.05844_wp !< mol. m. NaCl (kg mol-1)125 REAL(wp) :: molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1)126 REAL(wp) :: vanthoff = 2.0_wp !< van't Hoff factor for NaCl127 128 121 TYPE particle_type 129 122 SEQUENCE … … 164 157 TYPE(block_offset_def), DIMENSION(0:7) :: block_offset 165 158 159 ! 160 !-- Lagrangian cloud model constants (especially for steering aerosols) 161 REAL(wp) :: molecular_weight_of_solute = 0.05844_wp !< mol. m. NaCl (kg mol-1) 162 REAL(wp) :: molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1) 163 REAL(wp) :: rho_s = 2165.0_wp !< density of NaCl (kg m-3) 164 REAL(wp) :: vanthoff = 2.0_wp !< van't Hoff factor for NaCl 165 166 REAL(wp) :: n1 = 100.0_wp, s1 = 2.0_wp, rm1 = 0.05E-6_wp, & 167 n2 = 0.0_wp, s2 = 2.0_wp, rm2 = 0.05E-6_wp, & 168 n3 = 0.0_wp, s3 = 2.0_wp, rm3 = 0.05E-6_wp 169 170 LOGICAL :: monodisperse_aerosols = .FALSE. !< initialize monodisperse aerosols 171 LOGICAL :: init_aerosol_probabilistic = .FALSE. !< how to initialize aerosol spectra 172 166 173 SAVE 167 174 -
palm/trunk/SOURCE/package_parin.f90
r1834 r1871 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Initialization of aerosols added. 22 22 ! 23 23 ! Former revisions: … … 151 151 ONLY: alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, & 152 152 collision_algorithm, collision_kernel, & 153 curvature_solution_effects, density_ratio, & 154 dissipation_classes, dt_min_part, dt_prel, & 155 dt_write_particle_data, & 153 curvature_solution_effects, density_ratio, dissipation_classes, & 154 dt_min_part, dt_prel, dt_write_particle_data, & 156 155 end_time_prel, initial_weighting_factor, & 157 min_nr_particle, number_of_particle_groups, particles_per_point,& 158 particle_advection, particle_advection_start, pdx, pdy, pdz, & 159 psb, psl, psn, psr, pss, pst, radius, radius_classes, & 160 random_start_position, read_particles_from_restartfile, & 161 seed_follows_topography, use_sgs_for_particles, & 156 init_aerosol_probabilistic, min_nr_particle, & 157 monodisperse_aerosols, n1, n2, n3, number_of_particle_groups, & 158 particles_per_point, particle_advection, & 159 particle_advection_start, & 160 pdx, pdy, pdz, psb, psl, psn, psr, pss, pst, radius, & 161 radius_classes, random_start_position, & 162 read_particles_from_restartfile, rm1, rm2, rm3, & 163 seed_follows_topography, s1, s2, s3, use_sgs_for_particles, & 162 164 vertical_particle_advection, write_particle_statistics 163 165 … … 191 193 dt_write_particle_data, & 192 194 end_time_prel, initial_weighting_factor, & 193 min_nr_particle, & 194 number_of_particle_groups, & 195 init_aerosol_probabilistic, & 196 min_nr_particle, monodisperse_aerosols, & 197 n1, n2, n3, number_of_particle_groups, & 195 198 particles_per_point, & 196 199 particle_advection_start, & … … 199 202 radius_classes, random_start_position, & 200 203 read_particles_from_restartfile, & 201 seed_follows_topography, & 204 rm1, rm2, rm3, & 205 seed_follows_topography, s1, s2, s3, & 202 206 use_sgs_for_particles, & 203 207 vertical_particle_advection, &
Note: See TracChangeset
for help on using the changeset viewer.
