Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Jul 14, 2017 8:26:51 PM (7 years ago)
Author:
hoffmann
Message:

various improvements of the LCM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/lpm_init.f90

    r2305 r2312  
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Extended particle data type. Aerosol initialization improved.
     28!
     29! 2305 2017-07-06 11:18:47Z hoffmann
    2730! Improved calculation of particle IDs.
    28 ! 
     31!
    2932! 2274 2017-06-09 13:27:48Z Giersch
    3033!  Changed error messages
    31 ! 
     34!
    3235! 2265 2017-06-08 16:58:28Z schwenkel
    3336! Unused variables removed.
    34 ! 
     37!
    3538! 2263 2017-06-08 14:59:01Z schwenkel
    3639! Implemented splitting and merging algorithm
    37 ! 
     40!
    3841! 2233 2017-05-30 18:08:54Z suehring
    3942!
    4043! 2232 2017-05-30 17:47:52Z suehring
    4144! Adjustments according to new topography realization
    42 ! 
    43 ! 
     45!
     46!
    4447! 2223 2017-05-15 16:38:09Z suehring
    4548! Add check for particle release at model top
    46 ! 
     49!
    4750! 2182 2017-03-17 14:27:40Z schwenkel
    4851! Added parameters for simplified particle initialization.
    4952!
    5053! 2122 2017-01-18 12:22:54Z hoffmann
    51 ! Improved initialization of equilibrium aerosol radii 
     54! Improved initialization of equilibrium aerosol radii
    5255! Calculation of particle ID
    5356!
    5457! 2000 2016-08-20 18:09:15Z knoop
    5558! Forced header and separation lines into 80 columns
    56 ! 
     59!
    5760! 2016-06-09 16:25:25Z suehring
    58 ! Bugfix in determining initial particle height and grid index in case of 
     61! Bugfix in determining initial particle height and grid index in case of
    5962! seed_follows_topography.
    60 ! Bugfix concerning random positions, ensure that particles do not move more 
     63! Bugfix concerning random positions, ensure that particles do not move more
    6164! than one grid length.
    6265! Bugfix logarithmic interpolation.
    … …  
    6467!
    6568! 1890 2016-04-22 08:52:11Z hoffmann
    66 ! Initialization of aerosol equilibrium radius not possible in supersaturated 
     69! Initialization of aerosol equilibrium radius not possible in supersaturated
    6770! environments. Therefore, a maximum supersaturation of -1 % is assumed during
    6871! initialization.
    … …  
    7073! 1873 2016-04-18 14:50:06Z maronga
    7174! Module renamed (removed _mod
    72 ! 
     75!
    7376! 1871 2016-04-15 11:46:09Z hoffmann
    7477! Initialization of aerosols added.
    … …  
    8689! netcdf module added
    8790!
    88 ! 1725 2015-11-17 13:01:51Z hoffmann 
    89 ! Bugfix: Processor-dependent seed for random function is generated before it is 
     91! 1725 2015-11-17 13:01:51Z hoffmann
     92! Bugfix: Processor-dependent seed for random function is generated before it is
    9093! used.
    9194!
    … …  
    97100!
    98101! 1682 2015-10-07 23:56:08Z knoop
    99 ! Code annotations made doxygen readable 
     102! Code annotations made doxygen readable
    100103!
    101104! 1575 2015-03-27 09:56:27Z raasch
    … …  
    103106!
    104107! 1359 2014-04-11 17:15:14Z hoffmann
    105 ! New particle structure integrated. 
     108! New particle structure integrated.
    106109! Kind definition added to all floating point numbers.
    107110! lpm_init changed form a subroutine to a module.
    108 ! 
     111!
    109112! 1327 2014-03-21 11:00:16Z raasch
    110113! -netcdf_output
    … …  
    123126!
    124127! 1314 2014-03-14 18:25:17Z suehring
    125 ! Vertical logarithmic interpolation of horizontal particle speed for particles 
     128! Vertical logarithmic interpolation of horizontal particle speed for particles
    126129! between roughness height and first vertical grid level.
    127130!
    … …  
    146149!
    147150! 667 2010-12-23 12:06:00Z suehring/gryschka
    148 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation 
     151! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
    149152! of arrays.
    150153!
    … …  
    159162!------------------------------------------------------------------------------!
    160163 MODULE lpm_init_mod
    161  
     164
    162165    USE, INTRINSIC ::  ISO_C_BINDING
    163166
    … …  
    167170    USE control_parameters,                                                    &
    168171        ONLY:  cloud_droplets, constant_flux_layer, current_timestep_number,   &
    169                dz, initializing_actions, message_string, ocean, simulated_time
     172               dt_3d, dz, initializing_actions, message_string, ocean,         &
     173               simulated_time
    170174
    171175    USE grid_variables,                                                        &
    … …  
    197201                particles, particle_advection_start, particle_groups,          &
    198202                particle_groups_type, particles_per_point,                     &
    199                 particle_type, pdx, pdy, pdz,  prt_count, psb, psl, psn, psr,  & 
    200                 pss, pst, radius, random_start_position,                       & 
     203                particle_type, pdx, pdy, pdz,  prt_count, psb, psl, psn, psr,  &
     204                pss, pst, radius, random_start_position,                       &
    201205                read_particles_from_restartfile, seed_follows_topography,      &
    202206                sgs_wf_part, sort_count, splitting_function, splitting_mode,   &
    … …  
    247251    INTEGER(iwp) ::  k                           !<
    248252
    249     REAL(wp) ::  div                             !< 
     253    REAL(wp) ::  div                             !<
    250254    REAL(wp) ::  height_int                      !<
    251255    REAL(wp) ::  height_p                        !<
    … …  
    286290!-- Check if downward-facing walls exist. This case, reflection boundary
    287291!-- conditions (as well as subgrid-scale velocities) may do not work
    288 !-- propably (not realized so far). 
     292!-- propably (not realized so far).
    289293    IF ( surf_def_h(1)%ns >= 1 )  THEN
    290        WRITE( message_string, * ) 'Overhanging topography do not work '// &   
    291                                   'with particles' 
     294       WRITE( message_string, * ) 'Overhanging topography do not work '// &
     295                                  'with particles'
    292296       CALL message( 'lpm_init', 'PA0212', 0, 1, 0, 6, 0 )
    293297
    294     ENDIF 
     298    ENDIF
    295299
    296300!
    … …  
    310314!-- If number_particles_per_gridbox is set, the parametres pdx, pdy and pdz are
    311315!-- calculated diagnostically. Therfore an isotropic distribution is prescribed.
    312     IF ( number_particles_per_gridbox /= -1 .AND.   & 
     316    IF ( number_particles_per_gridbox /= -1 .AND.   &
    313317         number_particles_per_gridbox >= 1 )    THEN
    314        pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) /  & 
     318       pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) /  &
    315319             REAL(number_particles_per_gridbox))**0.3333333_wp
    316320!
    317 !--    Ensure a smooth value (two significant digits) of distance between 
     321!--    Ensure a smooth value (two significant digits) of distance between
    318322!--    particles (pdx, pdy, pdz).
    319323       div = 1000.0_wp
    … …  
    340344
    341345!
    342 !-- Allocate arrays required for calculating particle SGS velocities. 
     346!-- Allocate arrays required for calculating particle SGS velocities.
    343347!-- Initialize prefactor required for stoachastic Weil equation.
    344348    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    … …  
    347351                 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    348352
    349        sgs_wf_part = 1.0_wp / 3.0_wp   
     353       sgs_wf_part = 1.0_wp / 3.0_wp
    350354    ENDIF
    351355
    352356!
    353 !-- Allocate array required for logarithmic vertical interpolation of 
     357!-- Allocate array required for logarithmic vertical interpolation of
    354358!-- horizontal particle velocities between the surface and the first vertical
    355 !-- grid level. In order to avoid repeated CPU cost-intensive CALLS of 
    356 !-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for 
    357 !-- several heights. Splitting into 20 sublayers turned out to be sufficient. 
     359!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
     360!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
     361!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
    358362!-- To obtain exact height levels of particles, linear interpolation is applied
    359363!-- (see lpm_advec.f90).
    360364    IF ( constant_flux_layer )  THEN
    361        
    362        ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 
     365
     366       ALLOCATE ( log_z_z0(0:number_of_sublayers) )
    363367       z_p = zu(nzb+1) - zw(nzb)
    364368
    365369!
    366370!--    Calculate horizontal mean value of z0 used for logartihmic
    367 !--    interpolation. Note: this is not exact for heterogeneous z0. 
    368 !--    However, sensitivity studies showed that the effect is 
    369 !--    negligible. 
     371!--    interpolation. Note: this is not exact for heterogeneous z0.
     372!--    However, sensitivity studies showed that the effect is
     373!--    negligible.
    370374       z0_av_local  = SUM( surf_def_h(0)%z0 ) + SUM( surf_lsm_h%z0 ) +         &
    371375                      SUM( surf_usm_h%z0 )
    … …  
    401405!-- Check boundary condition and set internal variables
    402406    SELECT CASE ( bc_par_b )
    403    
     407
    404408       CASE ( 'absorb' )
    405409          ibc_par_b = 1
    … …  
    407411       CASE ( 'reflect' )
    408412          ibc_par_b = 2
    409          
     413
    410414       CASE DEFAULT
    411415          WRITE( message_string, * )  'unknown boundary condition ',   &
    412416                                       'bc_par_b = "', TRIM( bc_par_b ), '"'
    413417          CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
    414          
     418
    415419    END SELECT
    416420    SELECT CASE ( bc_par_t )
    417    
     421
    418422       CASE ( 'absorb' )
    419423          ibc_par_t = 1
    … …  
    421425       CASE ( 'reflect' )
    422426          ibc_par_t = 2
    423          
     427
    424428       CASE DEFAULT
    425429          WRITE( message_string, * ) 'unknown boundary condition ',   &
    426430                                     'bc_par_t = "', TRIM( bc_par_t ), '"'
    427431          CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
    428          
     432
    429433    END SELECT
    430434    SELECT CASE ( bc_par_lr )
    … …  
    438442       CASE ( 'reflect' )
    439443          ibc_par_lr = 2
    440          
     444
    441445       CASE DEFAULT
    442446          WRITE( message_string, * ) 'unknown boundary condition ',   &
    443447                                     'bc_par_lr = "', TRIM( bc_par_lr ), '"'
    444448          CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
    445          
     449
    446450    END SELECT
    447451    SELECT CASE ( bc_par_ns )
    … …  
    455459       CASE ( 'reflect' )
    456460          ibc_par_ns = 2
    457          
     461
    458462       CASE DEFAULT
    459463          WRITE( message_string, * ) 'unknown boundary condition ',   &
    460464                                     'bc_par_ns = "', TRIM( bc_par_ns ), '"'
    461465          CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
    462          
     466
    463467    END SELECT
    464468    SELECT CASE ( splitting_mode )
    465    
     469
    466470       CASE ( 'const' )
    467471          i_splitting_mode = 1
    … …  
    472476       CASE ( 'gb_av' )
    473477          i_splitting_mode = 3
    474          
     478
    475479       CASE DEFAULT
    476480          WRITE( message_string, * )  'unknown splitting condition ',   &
    477481                                       'splitting_mode = "', TRIM( splitting_mode ), '"'
    478482          CALL message( 'lpm_init', 'PA0146', 1, 2, 0, 6, 0 )
    479          
     483
    480484    END SELECT
    481485    SELECT CASE ( splitting_function )
    482    
     486
    483487       CASE ( 'gamma' )
    484488          isf = 1
    … …  
    489493       CASE ( 'exp' )
    490494          isf = 3
    491          
     495
    492496       CASE DEFAULT
    493497          WRITE( message_string, * )  'unknown splitting function ',   &
    494498                                       'splitting_function = "', TRIM( splitting_function ), '"'
    495499          CALL message( 'lpm_init', 'PA0147', 1, 2, 0, 6, 0 )
    496          
     500
    497501    END SELECT
    498    
     502
    499503
    500504!
    … …  
    524528
    525529!
    526 !--    initialize counter for particle IDs 
     530!--    initialize counter for particle IDs
    527531       grid_particles%id_counter = 1
    528532
    … …  
    534538                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    535539                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    536                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
     540                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    537541                                      0, 0, 0_idp, .FALSE., -1 )
    538542
    … …  
    577581          CALL check_open( 80 )
    578582          WRITE ( 80, 8000 )  current_timestep_number, simulated_time,         &
    579                               number_of_particles                             
     583                              number_of_particles
    580584          CALL close_file( 80 )
    581585       ENDIF
    … …  
    584588
    585589!
    586 !-- To avoid programm abort, assign particles array to the local version of 
     590!-- To avoid programm abort, assign particles array to the local version of
    587591!-- first grid cell
    588592    number_of_particles = prt_count(nzb+1,nys,nxl)
    … …  
    608612
    609613    USE particle_attributes,                                                   &
    610         ONLY: deleted_particles, monodisperse_aerosols
     614        ONLY: deleted_particles
    611615
    612616    IMPLICIT  NONE
    … …  
    631635    LOGICAL                    ::  first_stride !< flag for initialization
    632636
    633     REAL(wp)                   ::  pos_x      !< increment for particle position in x     
    634     REAL(wp)                   ::  pos_y      !< increment for particle position in y 
    635     REAL(wp)                   ::  pos_z      !< increment for particle position in z     
     637    REAL(wp)                   ::  pos_x      !< increment for particle position in x
     638    REAL(wp)                   ::  pos_y      !< increment for particle position in y
     639    REAL(wp)                   ::  pos_z      !< increment for particle position in z
    636640    REAL(wp)                   ::  rand_contr !< dummy argument for random position
    637641
    … …  
    652656!--    Calculate initial_weighting_factor diagnostically
    653657       IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN
    654           initial_weighting_factor =  number_concentration * 1.0E6_wp *             & 
    655                                       pdx(1) * pdy(1) * pdz(1) 
     658          initial_weighting_factor =  number_concentration * 1.0E6_wp *             &
     659                                      pdx(1) * pdy(1) * pdz(1)
    656660       END IF
    657661
    … …  
    677681               xloop: DO WHILE ( pos_x <= psr(i) )
    678682
    679                          IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            & 
     683                         IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            &
    680684                              pos_x <  ( nxr + 0.5_wp ) * dx )  THEN
    681685
    … …  
    690694                               tmp_particle%age_m         = 0.0_wp
    691695                               tmp_particle%dt_sum        = 0.0_wp
    692                                tmp_particle%user          = 0.0_wp !unused, free for the user
    693696                               tmp_particle%e_m           = 0.0_wp
    694                                IF ( curvature_solution_effects )  THEN
    695 !
    696 !--                               Initial values (internal timesteps, derivative)
    697 !--                               for Rosenbrock method
    698                                   tmp_particle%rvar1      = 1.0E-6_wp     !last Rosenbrock timestep
    699                                   tmp_particle%rvar2      = 0.1E-6_wp     !dry aerosol radius
    700                                   tmp_particle%rvar3      = -9999999.9_wp !unused in this configuration
    701                                ELSE
    702 !
    703 !--                               Initial values for SGS velocities
    704                                   tmp_particle%rvar1      = 0.0_wp
    705                                   tmp_particle%rvar2      = 0.0_wp
    706                                   tmp_particle%rvar3      = 0.0_wp
    707                                ENDIF
     697                               tmp_particle%rvar1         = 0.0_wp
     698                               tmp_particle%rvar2         = 0.0_wp
     699                               tmp_particle%rvar3         = 0.0_wp
    708700                               tmp_particle%speed_x       = 0.0_wp
    709701                               tmp_particle%speed_y       = 0.0_wp
    … …  
    712704                               tmp_particle%origin_y      = pos_y
    713705                               tmp_particle%origin_z      = pos_z
     706                               IF ( curvature_solution_effects )  THEN
     707                                  tmp_particle%aux1      = 0.0_wp    ! dry aerosol radius
     708                                  tmp_particle%aux2      = dt_3d     ! last Rosenbrock timestep
     709                               ELSE
     710                                  tmp_particle%aux1      = 0.0_wp    ! free to use
     711                                  tmp_particle%aux2      = 0.0_wp    ! free to use
     712                               ENDIF
    714713                               tmp_particle%radius        = particle_groups(i)%radius
    715714                               tmp_particle%weight_factor = initial_weighting_factor
    … …  
    726725!
    727726!--                            Determine surface level. Therefore, check for
    728 !--                            upward-facing wall on w-grid. MAXLOC will return 
     727!--                            upward-facing wall on w-grid. MAXLOC will return
    729728!--                            the index of the lowest upward-facing wall.
    730729                               k_surf = MAXLOC(                                &
    … …  
    748747                                  ENDIF
    749748!
    750 !--                            Skip particle release if particle position is 
     749!--                            Skip particle release if particle position is
    751750!--                            below surface, or within topography in case
    752751!--                            of overhanging structures.
    … …  
    756755                               THEN
    757756                                  pos_x = pos_x + pdx(i)
    758                                   CYCLE xloop                               
     757                                  CYCLE xloop
    759758                               ENDIF
    760759
    … …  
    840839             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    841840
    842              DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
     841             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
    843842
    844843                particles(n)%id = 10000_idp**3 * grid_particles(kp,jp,ip)%id_counter + &
    … …  
    857856!
    858857!-- Initialize aerosol background spectrum
    859     IF ( curvature_solution_effects  .AND.  .NOT. monodisperse_aerosols )  THEN
     858    IF ( curvature_solution_effects )  THEN
    860859       CALL lpm_init_aerosols(local_start)
    861860    ENDIF
    … …  
    871870                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    872871!
    873 !--             Move only new particles. Moreover, limit random fluctuation 
    874 !--             in order to prevent that particles move more than one grid box, 
    875 !--             which would lead to problems concerning particle exchange 
    876 !--             between processors in case pdx/pdy are larger than dx/dy, 
    877 !--             respectively. 
     872!--             Move only new particles. Moreover, limit random fluctuation
     873!--             in order to prevent that particles move more than one grid box,
     874!--             which would lead to problems concerning particle exchange
     875!--             between processors in case pdx/pdy are larger than dx/dy,
     876!--             respectively.
    878877                DO  n = local_start(kp,jp,ip), number_of_particles
    879878                   IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
    … …  
    883882                              MERGE( rand_contr, SIGN( dx, rand_contr ),       &
    884883                                     ABS( rand_contr ) < dx                    &
    885                                    ) 
     884                                   )
    886885                   ENDIF
    887886                   IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
    … …  
    891890                              MERGE( rand_contr, SIGN( dy, rand_contr ),       &
    892891                                     ABS( rand_contr ) < dy                    &
    893                                    ) 
     892                                   )
    894893                   ENDIF
    895894                   IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
    … …  
    899898                              MERGE( rand_contr, SIGN( dz, rand_contr ),       &
    900899                                     ABS( rand_contr ) < dz                    &
    901                                    ) 
     900                                   )
    902901                   ENDIF
    903902                ENDDO
    … …  
    908907!
    909908!--             Furthermore, remove particles located in topography. Note, as
    910 !--             the particle speed is still zero at this point, wall 
     909!--             the particle speed is still zero at this point, wall
    911910!--             reflection boundary conditions will not work in this case.
    912911                particles =>                                                   &
    … …  
    934933    ENDIF
    935934!
    936 !-- In case of random_start_position, delete particles identified by 
    937 !-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks, 
    938 !-- which is needed for a fast interpolation of the LES fields on the particle 
     935!-- In case of random_start_position, delete particles identified by
     936!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks,
     937!-- which is needed for a fast interpolation of the LES fields on the particle
    939938!-- position.
    940939    CALL lpm_pack_all_arrays
    … …  
    967966
    968967    USE arrays_3d,                                                             &
    969         ONLY: hyp, pt, q 
     968        ONLY: hyp, pt, q
    970969
    971970    USE cloud_parameters,                                                      &
    … …  
    978977
    979978    USE particle_attributes,                                                   &
    980         ONLY: init_aerosol_probabilistic, molecular_weight_of_solute,          &
    981               molecular_weight_of_water, n1, n2, n3, rho_s, rm1, rm2, rm3,     &
    982               s1, s2, s3, vanthoff
     979        ONLY: aero_type, aero_weight, log_sigma, molecular_weight_of_solute,   &
     980              molecular_weight_of_water, na, rho_s, rm, vanthoff
    983981
    984982    IMPLICIT NONE
    985 
    986     REAL(wp), DIMENSION(:), ALLOCATABLE ::  cdf     !< CDF of aerosol spectrum
    987     REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_temp  !< dry aerosol radius spectrum
    988983
    989984    REAL(wp)  :: afactor            !< curvature effects
    990985    REAL(wp)  :: bfactor            !< solute effects
    991     REAL(wp)  :: dr                 !< width of radius bin
     986    REAL(wp)  :: dlogr              !< logarithmic width of radius bin
    992987    REAL(wp)  :: e_a                !< vapor pressure
    993988    REAL(wp)  :: e_s                !< saturation vapor pressure
    994     REAL(wp)  :: n_init             !< sum of all aerosol concentrations
    995     REAL(wp)  :: pdf                !< PDF of aerosol spectrum
    996     REAL(wp)  :: rmin = 1.0e-8_wp   !< minimum aerosol radius
    997     REAL(wp)  :: rmax = 1.0e-6_wp   !< maximum aerosol radius
    998     REAL(wp)  :: rs_rand            !< random number
    999     REAL(wp)  :: r_mid              !< mean radius
     989    REAL(wp)  :: rmin = 0.005e-6_wp !< minimum aerosol radius
     990    REAL(wp)  :: rmax = 10.0e-6_wp  !< maximum aerosol radius
     991    REAL(wp)  :: r_mid              !< mean radius of bin
     992    REAL(wp)  :: r_l                !< left radius of bin
     993    REAL(wp)  :: r_r                !< right radius of bin
    1000994    REAL(wp)  :: sigma              !< surface tension
    1001995    REAL(wp)  :: t_int              !< temperature
    1002     REAL(wp)  :: weight_sum         !< sum of all weighting factors
    1003996
    1004997    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  local_start !<
    1005998
    1006999    INTEGER(iwp)  :: n              !<
    1007     INTEGER(iwp)  :: nn             !<
    1008     INTEGER(iwp)  :: no_bins = 999  !< number of bins
    10091000    INTEGER(iwp)  :: ip             !<
    10101001    INTEGER(iwp)  :: jp             !<
    10111002    INTEGER(iwp)  :: kp             !<
    10121003
    1013     LOGICAL ::  new_pdf = .FALSE.   !< check if aerosol PDF has to be recalculated
    1014 
    1015 !
    1016 !-- Compute aerosol background distribution
    1017     IF ( init_aerosol_probabilistic )  THEN
    1018        ALLOCATE( cdf(0:no_bins), r_temp(0:no_bins) )
    1019        DO n = 0, no_bins
    1020           r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) /            &
    1021                            REAL(no_bins, KIND=wp) * REAL(n, KIND=wp) )
    1022 
    1023           cdf(n) = 0.0_wp
    1024           n_init = n1 + n2 + n3
    1025           IF ( n1 > 0.0_wp )  THEN
    1026              cdf(n) = cdf(n) + n1 / n_init * ( 0.5_wp + 0.5_wp *        &
    1027                                   ERF( LOG( r_temp(n) / rm1 ) /         &
    1028                                        ( SQRT(2.0_wp) * LOG(s1) )       &
    1029                                      ) )
    1030           ENDIF
    1031           IF ( n2 > 0.0_wp )  THEN
    1032              cdf(n) = cdf(n) + n2 / n_init * ( 0.5_wp + 0.5_wp *        &
    1033                                   ERF( LOG( r_temp(n) / rm2 ) /         &
    1034                                        ( SQRT(2.0_wp) * LOG(s2) )       &
    1035                                      ) )
    1036           ENDIF
    1037           IF ( n3 > 0.0_wp )  THEN
    1038              cdf(n) = cdf(n) + n3 / n_init * ( 0.5_wp + 0.5_wp *        &
    1039                                   ERF( LOG( r_temp(n) / rm3 ) /         &
    1040                                        ( SQRT(2.0_wp) * LOG(s3) )       &
    1041                                      ) )
    1042           ENDIF
    1043 
    1044        ENDDO
     1004!
     1005!-- The following typical aerosol spectra are taken from Jaenicke (1993):
     1006!-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
     1007    IF ( TRIM(aero_type) .EQ. 'polar' )  THEN
     1008       na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6
     1009       rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6
     1010       log_sigma = (/ 0.245, 0.300, 0.291 /)
     1011    ELSEIF ( TRIM(aero_type) .EQ. 'background' )  THEN
     1012       na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6
     1013       rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6
     1014       log_sigma = (/ 0.645, 0.253, 0.425 /)
     1015    ELSEIF ( TRIM(aero_type) .EQ. 'maritime' )  THEN
     1016       na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6
     1017       rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6
     1018       log_sigma = (/ 0.657, 0.210, 0.396 /)
     1019    ELSEIF ( TRIM(aero_type) .EQ. 'continental' )  THEN
     1020       na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6
     1021       rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6
     1022       log_sigma = (/ 0.161, 0.217, 0.380 /)
     1023    ELSEIF ( TRIM(aero_type) .EQ. 'desert' )  THEN
     1024       na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6
     1025       rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6
     1026       log_sigma = (/ 0.247, 0.770, 0.438 /)
     1027    ELSEIF ( TRIM(aero_type) .EQ. 'rural' )  THEN
     1028       na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6
     1029       rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6
     1030       log_sigma = (/ 0.225, 0.557, 0.266 /)
     1031    ELSEIF ( TRIM(aero_type) .EQ. 'urban' )  THEN
     1032       na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6
     1033       rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6
     1034       log_sigma = (/ 0.245, 0.666, 0.337 /)
     1035    ELSEIF ( TRIM(aero_type) .EQ. 'user' )  THEN
     1036       CONTINUE
     1037    ELSE
     1038       WRITE( message_string, * ) 'unknown aerosol type ',   &
     1039                                'aero_type = "', TRIM( aero_type ), '"'
     1040       CALL message( 'lpm_init', 'PA0459', 1, 2, 0, 6, 0 )
    10451041    ENDIF
    10461042
    … …  
    10521048             IF ( number_of_particles <= 0 )  CYCLE
    10531049             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     1050
     1051             dlogr   = ( LOG10(rmax) - LOG10(rmin) ) / ( number_of_particles - local_start(kp,jp,ip) + 1 )
    10541052!
    10551053!--          Initialize the aerosols with a predefined spectral distribution
    1056 !--          of the dry radius (logarithmically increasing bins) and a varying 
     1054!--          of the dry radius (logarithmically increasing bins) and a varying
    10571055!--          weighting factor
    1058              IF ( .NOT. init_aerosol_probabilistic )  THEN
    1059 
    1060                 new_pdf = .FALSE.
    1061                 IF ( .NOT. ALLOCATED( r_temp ) )  THEN
    1062                    new_pdf = .TRUE.
     1056             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
     1057
     1058                r_l   = 10.0**( LOG10( rmin ) + (n-1) * dlogr )
     1059                r_r   = 10.0**( LOG10( rmin ) + n * dlogr )
     1060                r_mid = SQRT( r_l * r_r )
     1061
     1062                particles(n)%aux1          = r_mid
     1063                particles(n)%weight_factor =                                           &
     1064                   ( na(1) / ( SQRT( 2.0 * pi ) * log_sigma(1) ) *                     &
     1065                     EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0 * log_sigma(1)**2 ) ) +  &
     1066                     na(2) / ( SQRT( 2.0 * pi ) * log_sigma(2) ) *                     &
     1067                     EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0 * log_sigma(2)**2 ) ) +  &
     1068                     na(3) / ( SQRT( 2.0 * pi ) * log_sigma(3) ) *                     &
     1069                     EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0 * log_sigma(3)**2 ) )    &
     1070                   ) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dz )
     1071
     1072!
     1073!--             Multiply weight_factor with the namelist parameter aero_weight
     1074!--             to increase or decrease the number of simulated aerosols
     1075                particles(n)%weight_factor = particles(n)%weight_factor * aero_weight
     1076
     1077                IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) &
     1078                     .GT. random_function( iran_part ) )  THEN
     1079                   particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0
    10631080                ELSE
    1064                    IF ( SIZE( r_temp ) .NE. &
    1065                         number_of_particles - local_start(kp,jp,ip) + 2 )  THEN
    1066                       new_pdf = .TRUE.
    1067                       DEALLOCATE( r_temp )
    1068                    ENDIF
     1081                   particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp)
    10691082                ENDIF
    1070 
    1071                 IF ( new_pdf )  THEN
    1072 
    1073                    no_bins = number_of_particles + 1 - local_start(kp,jp,ip)
    1074                    ALLOCATE( r_temp(0:no_bins) )
    1075 
    1076                    DO n = 0, no_bins
    1077                       r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / &
    1078                                        REAL(no_bins, KIND=wp) *                 &
    1079                                        REAL(n, KIND=wp) )
    1080                    ENDDO
    1081 
    1082                 ENDIF
    1083 
    1084 !
    1085 !--             Calculate radius and concentration of each aerosol
    1086                 DO n = local_start(kp,jp,ip), number_of_particles
    1087 
    1088                    nn = n - local_start(kp,jp,ip)
    1089 
    1090                    r_mid = SQRT( r_temp(nn) * r_temp(nn+1) )
    1091                    dr    = r_temp(nn+1) - r_temp(nn)
    1092 
    1093                    pdf    = 0.0_wp
    1094                    n_init = n1 + n2 + n3
    1095                    IF ( n1 > 0.0_wp )  THEN
    1096                       pdf = pdf + n1 / n_init * ( 1.0_wp / ( r_mid * LOG(s1) *      &
    1097                                                              SQRT( 2.0_wp * pi )    &
    1098                                                            ) *                      &
    1099                                                   EXP( -( LOG( r_mid / rm1 ) )**2 / &
    1100                                                        ( 2.0_wp * LOG(s1)**2 )      &
    1101                                                      )                              &
    1102                                                 )
    1103                    ENDIF
    1104                    IF ( n2 > 0.0_wp )  THEN
    1105                       pdf = pdf + n2 / n_init * ( 1.0_wp / ( r_mid * LOG(s2) *      &
    1106                                                              SQRT( 2.0_wp * pi )    &
    1107                                                            ) *                      &
    1108                                                   EXP( -( LOG( r_mid / rm2 ) )**2 / &
    1109                                                        ( 2.0_wp * LOG(s2)**2 )      &
    1110                                                      )                              &
    1111                                                 )
    1112                    ENDIF
    1113                    IF ( n3 > 0.0_wp )  THEN
    1114                       pdf = pdf + n3 / n_init * ( 1.0_wp / ( r_mid * LOG(s3) *      &
    1115                                                              SQRT( 2.0_wp * pi )    &
    1116                                                            ) *                      &
    1117                                                   EXP( -( LOG( r_mid / rm3 ) )**2 / &
    1118                                                        ( 2.0_wp * LOG(s3)**2 )      &
    1119                                                      )                              &
    1120                                                 )
    1121                    ENDIF
    1122 
    1123                    particles(n)%rvar2         = r_mid
    1124                    particles(n)%weight_factor = pdf * dr
    1125 
    1126                 END DO
    1127 !
    1128 !--             Adjust weighting factors to initialize the same number of aerosols
    1129 !--             in every grid box
    1130                 weight_sum = SUM(particles(local_start(kp,jp,ip):number_of_particles)%weight_factor)
    1131 
    1132                 particles(local_start(kp,jp,ip):number_of_particles)%weight_factor =     &
    1133                    particles(local_start(kp,jp,ip):number_of_particles)%weight_factor /  &
    1134                    weight_sum * initial_weighting_factor * ( no_bins + 1 )
    1135 
    1136              ENDIF
    1137 !
    1138 !--          Initialize the aerosols with a predefined weighting factor but
    1139 !--          a randomly choosen dry radius
    1140              IF ( init_aerosol_probabilistic )  THEN
    1141 
    1142                 DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
    1143 
    1144                    rs_rand = -1.0_wp
    1145                    DO WHILE ( rs_rand .LT. cdf(0)  .OR.  rs_rand .GE. cdf(no_bins)  )
    1146                       rs_rand = random_function( iran_part )
    1147                    ENDDO
    1148 !
    1149 !--                Determine aerosol dry radius by a random number generator
    1150                    DO nn = 0, no_bins-1
    1151                       IF ( cdf(nn) .LE. rs_rand  .AND.  cdf(nn+1) .GT. rs_rand )  THEN
    1152                          particles(n)%rvar2 = r_temp(nn) + ( r_temp(nn+1) - r_temp(nn) ) / &
    1153                                               ( cdf(nn+1) - cdf(nn) ) * ( rs_rand - cdf(nn) )
    1154                          EXIT
    1155                       ENDIF
    1156                    ENDDO
    1157 
    1158                 ENDDO
    1159 
    1160              ENDIF
    1161 
     1083!
     1084!--             Unnecessary particles will be deleted
     1085                IF ( particles(n)%weight_factor .LE. 0.0 )  particles(n)%particle_mask = .FALSE.
     1086
     1087             ENDDO
    11621088!
    11631089!--          Set particle radius to equilibrium radius based on the environmental
    11641090!--          supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids
    1165 !--          the sometimes lengthy growth toward their equilibrium radius within 
     1091!--          the sometimes lengthy growth toward their equilibrium radius within
    11661092!--          the simulation.
    11671093             t_int  = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
    … …  
    11761102                       rho_s / ( molecular_weight_of_solute * rho_l )
    11771103!
    1178 !--          The formula is only valid for subsaturated environments. For 
     1104!--          The formula is only valid for subsaturated environments. For
    11791105!--          supersaturations higher than -5 %, the supersaturation is set to -5%.
    11801106             IF ( e_a / e_s >= 0.95_wp )  e_a = 0.95_wp * e_s
    11811107
    1182              DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
    1183 !
    1184 !--             For details on this equation, see Eq. (14) of Khvorostyanov and 
    1185 !--             Curry (2007, JGR) 
     1108             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
     1109!
     1110!--             For details on this equation, see Eq. (14) of Khvorostyanov and
     1111!--             Curry (2007, JGR)
    11861112                particles(n)%radius = bfactor**0.3333333_wp *                  &
    1187                    particles(n)%rvar2 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
     1113                   particles(n)%aux1 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
    11881114                   ( 1.0_wp + ( afactor / ( 3.0_wp * bfactor**0.3333333_wp *   &
    1189                      particles(n)%rvar2 ) ) /                                  &
     1115                     particles(n)%aux1 ) ) /                                  &
    11901116                     ( 1.0_wp - e_a / e_s )**0.6666666_wp                      &
    11911117                   )
    … …  
    11961122       ENDDO
    11971123    ENDDO
    1198 !
    1199 !-- Deallocate used arrays
    1200     IF ( ALLOCATED(r_temp) )  DEALLOCATE( r_temp )
    1201     IF ( ALLOCATED(cdf) )     DEALLOCATE( cdf )
    12021124
    12031125 END SUBROUTINE lpm_init_aerosols
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |