Ignore:
Timestamp:
Jan 18, 2017 12:22:54 PM (5 years ago)
Author:
hoffmann
Message:

introduction of a particle ID, improvement of equilibrium radius calculation, and reformatting

File:
1 edited

Legend:

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

    r2101 r2122  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Improved initialization of equilibrium aerosol radii
     23! Calculation of particle ID
    2324!
    2425! Former revisions:
     
    4445! 1873 2016-04-18 14:50:06Z maronga
    4546! Module renamed (removed _mod
    46 
    4747!
    4848! 1871 2016-04-15 11:46:09Z hoffmann
     
    236236!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
    237237!-- particle_attributes). Integer length is 4 byte, Real is 8 byte
     238!    blocklengths(1)  = 20;  blocklengths(2)  =   5;  blocklengths(3)  =   1
     239!    displacements(1) =  0;  displacements(2) = 160;  displacements(3) = 180
    238240    blocklengths(1)  = 19;  blocklengths(2)  =   6;  blocklengths(3)  =   1
    239241    displacements(1) =  0;  displacements(2) = 152;  displacements(3) = 176
     242
     243!
     244!-- WARNING: the double precision integer 'id' is treated as a double precision
     245!-- integer here, i.e., there are 20 real variables and only 5 integer variables
     246!-- instead of 19 and 6, respectively.
    240247
    241248    types(1) = MPI_REAL
     
    452459
    453460!
     461!--    initialize counter for particle IDs
     462       grid_particles%id_counter = 0
     463
     464!
    454465!--    Initialize all particles with dummy values (otherwise errors may
    455466!--    occur within restart runs). The reason for this is still not clear
     
    458469                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    459470                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    460                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
    461                                       0, .FALSE., -1 )
     471                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,         &
     472                                      0, 0, 0, 0, .FALSE., -1 )
    462473
    463474       particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
     
    602613                            tmp_particle%age_m         = 0.0_wp
    603614                            tmp_particle%dt_sum        = 0.0_wp
    604                             tmp_particle%dvrp_psize    = 0.0_wp !unused
     615                            tmp_particle%user          = 0.0_wp !unused, free for the user
    605616                            tmp_particle%e_m           = 0.0_wp
    606617                            IF ( curvature_solution_effects )  THEN
     
    610621                               tmp_particle%rvar1      = 1.0E-6_wp     !last Rosenbrock timestep
    611622                               tmp_particle%rvar2      = 0.1E-6_wp     !dry aerosol radius
    612                                tmp_particle%rvar3      = -9999999.9_wp !unused
     623                               tmp_particle%rvar3      = -9999999.9_wp !unused in this configuration
    613624                            ELSE
    614625!
     
    628639                            tmp_particle%class         = 1
    629640                            tmp_particle%group         = i
    630                             tmp_particle%tailpoints    = 0     !unused
     641                            tmp_particle%id1           = 0
     642                            tmp_particle%id2           = 0
    631643                            tmp_particle%particle_mask = .TRUE.
    632                             tmp_particle%tail_id       = 0     !unused
    633 
    634 
     644                            tmp_particle%block_nr      = -1
    635645!
    636646!--                         Determine the grid indices of the particle position
     
    722732    local_start = prt_count+1
    723733    prt_count   = local_count
     734
     735!
     736!-- Calculate particle IDs
     737    DO  ip = nxl, nxr
     738       DO  jp = nys, nyn
     739          DO  kp = nzb+1, nzt
     740             number_of_particles = prt_count(kp,jp,ip)
     741             IF ( number_of_particles <= 0 )  CYCLE
     742             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     743
     744             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
     745
     746                particles(n)%id1 = 10000 * grid_particles(kp,jp,ip)%id_counter + kp
     747                particles(n)%id2 = 10000 * jp + ip
     748
     749                grid_particles(kp,jp,ip)%id_counter =                          &
     750                                         grid_particles(kp,jp,ip)%id_counter + 1
     751
     752             ENDDO
     753
     754          ENDDO
     755       ENDDO
     756    ENDDO
    724757
    725758!
     
    837870
    838871    USE cloud_parameters,                                                      &
    839         ONLY: l_d_rv, rho_l
     872        ONLY: l_d_rv, rho_l, r_v
    840873
    841874    USE constants,                                                             &
     
    854887    REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_temp  !< dry aerosol radius spectrum
    855888
     889    REAL(wp)  :: afactor            !< curvature effects
    856890    REAL(wp)  :: bfactor            !< solute effects
    857891    REAL(wp)  :: dr                 !< width of radius bin
     
    864898    REAL(wp)  :: rs_rand            !< random number
    865899    REAL(wp)  :: r_mid              !< mean radius
     900    REAL(wp)  :: sigma              !< surface tension
    866901    REAL(wp)  :: t_int              !< temperature
    867902    REAL(wp)  :: weight_sum         !< sum of all weighting factors
     
    10351070             e_a = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp )
    10361071
     1072             sigma   = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
     1073             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
     1074
     1075             bfactor = vanthoff * molecular_weight_of_water *    &
     1076                       rho_s / ( molecular_weight_of_solute * rho_l )
    10371077!
    10381078!--          The formula is only valid for subsaturated environments. For
    1039 !--          supersaturations higher than -1 %, the supersaturation is set to -1%.
    1040              IF ( e_a / e_s < 0.99_wp )  THEN
    1041 
    1042                 DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
    1043 
    1044                    bfactor             = vanthoff * molecular_weight_of_water *    &
    1045                                          rho_s * particles(n)%rvar2**3 /           &
    1046                                          ( molecular_weight_of_solute * rho_l )
    1047                    particles(n)%radius = particles(n)%rvar2 * ( bfactor /          &
    1048                                          particles(n)%rvar2**3 )**(1.0_wp/3.0_wp) *&
    1049                                          ( 1.0_wp - e_a / e_s )**(-1.0_wp/3.0_wp)
    1050 
    1051                 ENDDO
    1052 
    1053              ELSE
    1054 
    1055                 DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
    1056 
    1057                    bfactor             = vanthoff * molecular_weight_of_water *    &
    1058                                          rho_s * particles(n)%rvar2**3 /           &
    1059                                          ( molecular_weight_of_solute * rho_l )
    1060                    particles(n)%radius = particles(n)%rvar2 * ( bfactor /          &
    1061                                          particles(n)%rvar2**3 )**(1.0_wp/3.0_wp) *&
    1062                                          0.01_wp**(-1.0_wp/3.0_wp)
    1063 
    1064                 ENDDO
    1065 
    1066              ENDIF
     1079!--          supersaturations higher than -5 %, the supersaturation is set to -5%.
     1080             IF ( e_a / e_s >= 0.95_wp )  e_a = 0.95_wp * e_s
     1081
     1082             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
     1083!
     1084!--             For details on this equation, see Eq. (14) of Khvorostyanov and
     1085!--             Curry (2007, JGR)
     1086                particles(n)%radius = bfactor**0.3333333_wp *                  &
     1087                   particles(n)%rvar2 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
     1088                   ( 1.0_wp + ( afactor / ( 3.0_wp * bfactor**0.3333333_wp *   &
     1089                     particles(n)%rvar2 ) ) /                                  &
     1090                     ( 1.0_wp - e_a / e_s )**0.6666666_wp                      &
     1091                   )
     1092
     1093             ENDDO
    10671094
    10681095          ENDDO
Note: See TracChangeset for help on using the changeset viewer.