Changeset 1871 for palm/trunk/SOURCE/lpm_init_mod.f90
- Timestamp:
- Apr 15, 2016 11:46:09 AM (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.