Changeset 2122


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

Location:
palm/trunk/SOURCE
Files:
5 edited

Legend:

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

    r2101 r2122  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Calculation of particle ID
    2323!
    2424! Former revisions:
     
    136136!     CALL netcdf_handle_error( 'lpm_data_output_particles', 3 )
    137137!
    138 !     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize,&
    139 !                             start = (/ 1, prt_time_count /),                &
     138!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%user,     &
     139!                             start = (/ 1, prt_time_count /),               &
    140140!                             count = (/ maximum_number_of_particles /) )
    141141!     CALL netcdf_handle_error( 'lpm_data_output_particles', 4 )
     
    208208!
    209209!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16),                    &
    210 !                             particles%tailpoints,                          &
     210!                             particles%id2,                                 &
    211211!                             start = (/ 1, prt_time_count /),               &
    212212!                             count = (/ maximum_number_of_particles /) )
    213213!     CALL netcdf_handle_error( 'lpm_data_output_particles', 18 )
    214214!
    215 !     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id, &
     215!     nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%id1,    &
    216216!                             start = (/ 1, prt_time_count /),               &
    217217!                             count = (/ maximum_number_of_particles /) )
  • palm/trunk/SOURCE/lpm_droplet_collision.f90

    r2101 r2122  
    2020! Current revisions:
    2121! ------------------
    22 !
    23 ! 
     22! Some reformatting of the code.
     23!
    2424! Former revisions:
    2525! -----------------
     
    445445
    446446                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    447                          mass(n) = mass(n) + weight(n) * xm 
    448                          weight(m)    = weight(m)    - weight(n)
    449                          mass(m) = mass(m) - weight(n) * xm
     447                         mass(n)   = mass(n)  + weight(n) * xm 
     448                         weight(m) = weight(m) - weight(n)
     449                         mass(m)   = mass(m)  - weight(n) * xm
    450450                      ENDIF
    451451
     
    457457
    458458                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
    459                          mass(m) = mass(m) + weight(m) * xn
    460                          weight(n)    = weight(n)    - weight(m)
    461                          mass(n) = mass(n) - weight(m) * xn
     459                         mass(m)   = mass(m)  + weight(m) * xn
     460                         weight(n) = weight(n) - weight(m)
     461                         mass(n)   = mass(n)  - weight(m) * xn
    462462                      ENDIF
    463463                   ELSE
  • 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
  • palm/trunk/SOURCE/lpm_set_attributes.f90

    r2101 r2122  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! DVRP routine removed
    2323!
    2424! Former revisions:
     
    397397    ENDIF
    398398
    399 !
    400 !-- Set particle size for dvrp graphics
    401     IF ( particle_dvrpsize == 'absw' )  THEN
    402 
    403        DO  ip = nxl, nxr
    404           DO  jp = nys, nyn
    405              DO  kp = nzb+1, nzt
    406 
    407                 number_of_particles = prt_count(kp,jp,ip)
    408                 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    409                 IF ( number_of_particles <= 0 )  CYCLE
    410                 start_index = grid_particles(kp,jp,ip)%start_index
    411                 end_index   = grid_particles(kp,jp,ip)%end_index
    412 
    413                 ALLOCATE( xv(1:number_of_particles), &
    414                           yv(1:number_of_particles) )
    415 
    416                 xv = particles(1:number_of_particles)%x
    417                 yv = particles(1:number_of_particles)%y
    418                 zv = particles(1:number_of_particles)%z
    419 
    420                 DO  nb = 0,7
    421 
    422                    i = ip + block_offset(nb)%i_off
    423                    j = jp + block_offset(nb)%j_off
    424                    k = kp-1
    425 
    426                    DO  n = start_index(nb), end_index(nb)
    427 !
    428 !--                   Interpolate w-component to the current particle position
    429                       x  = xv(n) - i * dx
    430                       y  = yv(n) - j * dy
    431                       aa = x**2          + y**2
    432                       bb = ( dx - x )**2 + y**2
    433                       cc = x**2          + ( dy - y )**2
    434                       dd = ( dx - x )**2 + ( dy - y )**2
    435                       gg = aa + bb + cc + dd
    436 
    437                       w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) *     &
    438                                   w(k,j,i+1) + ( gg - cc ) * w(k,j+1,i) +      &
    439                                   ( gg - dd ) * w(k,j+1,i+1)                   &
    440                                 ) / ( 3.0_wp * gg )
    441 
    442                       IF ( k+1 == nzt+1 )  THEN
    443                          w_int = w_int_l
    444                       ELSE
    445                          w_int_u = ( ( gg - aa ) * w(k+1,j,i)   + ( gg - bb ) *  &
    446                                      w(k+1,j,i+1) + ( gg - cc ) * w(k+1,j+1,i) + &
    447                                      ( gg - dd ) * w(k+1,j+1,i+1)                &
    448                                    ) / ( 3.0_wp * gg )
    449                          w_int = w_int_l + ( zv(n) - zw(k) ) / dz *     &
    450                                            ( w_int_u - w_int_l )
    451                       ENDIF
    452 
    453 !
    454 !--                   Limit values by the given interval and normalize to
    455 !--                   interval [0,1]
    456                       w_int = ABS( w_int )
    457                       w_int = MIN( w_int, dvrpsize_interval(2) )
    458                       w_int = MAX( w_int, dvrpsize_interval(1) )
    459 
    460                       w_int = ( w_int - dvrpsize_interval(1) ) / &
    461                               ( dvrpsize_interval(2) - dvrpsize_interval(1) )
    462 
    463                       particles(n)%dvrp_psize = ( 0.25_wp + w_int * 0.6_wp ) * &
    464                                                 dx
    465 
    466                    ENDDO
    467                 ENDDO
    468 
    469                 DEALLOCATE( xv, yv, zv )
    470 
    471              ENDDO
    472           ENDDO
    473        ENDDO
    474 
    475     ENDIF
    476 
    477399    CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'stop' )
    478400
  • palm/trunk/SOURCE/mod_particle_attributes.f90

    r2101 r2122  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Calculation of particle ID
     23! Particle attribute dvrp_psize renamed to user: this attribute can be used by
     24! by the user to store any variable
    2325!
    2426! Former revisions:
     
    134136    TYPE particle_type
    135137        SEQUENCE
    136         REAL(wp)     ::  radius, age, age_m, dt_sum, dvrp_psize, e_m,          &
     138        REAL(wp)     ::  radius, age, age_m, dt_sum, user, e_m,                &
    137139                         origin_x, origin_y, origin_z, rvar1, rvar2, rvar3,    &
    138140                         speed_x, speed_y, speed_z, weight_factor, x, y, z
    139         INTEGER(iwp) ::  class, group, tailpoints, tail_id
     141        INTEGER(iwp) ::  class, group, id1, id2
    140142        LOGICAL      ::  particle_mask
    141143        INTEGER(iwp) ::  block_nr
     
    156158        INTEGER(iwp), DIMENSION(0:7)               ::  start_index
    157159        INTEGER(iwp), DIMENSION(0:7)               ::  end_index
     160        INTEGER(iwp)                               ::  id_counter
    158161        LOGICAL                                    ::  time_loop_done
    159162        TYPE(particle_type), POINTER, DIMENSION(:) ::  particles                !Particle array for this grid cell
Note: See TracChangeset for help on using the changeset viewer.