Ignore:
Timestamp:
Apr 15, 2016 11:46:09 AM (8 years ago)
Author:
hoffmann
Message:

initialization of aerosol spectra added

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lpm_init_mod.f90

    r1851 r1871  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Initialization of aerosols added.
    2222!
    2323! Former revisions:
     
    509509        ONLY: lpm_pack_all_arrays
    510510
     511    USE particle_attributes,                                                   &
     512        ONLY: monodisperse_aerosols
     513
    511514    IMPLICIT  NONE
    512515
     
    581584!--                            Initial values (internal timesteps, derivative)
    582585!--                            for Rosenbrock method
    583                                tmp_particle%rvar1      = 1.0E-12_wp
    584                                tmp_particle%rvar2      = 1.0E-3_wp
    585                                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
    586589                            ELSE
    587590!
     
    686689    local_start = prt_count+1
    687690    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
    688698!
    689699!-- Add random fluctuation to particle positions
     
    759769 END SUBROUTINE lpm_create_particle
    760770
     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
    7611003END MODULE lpm_init_mod
Note: See TracChangeset for help on using the changeset viewer.