Changeset 4131 for palm/trunk


Ignore:
Timestamp:
Aug 2, 2019 11:06:18 AM (2 years ago)
Author:
monakurppa
Message:

Several changes in the salsa aerosol module:

  • Add "salsa_" before each salsa output variable
  • Add a possibility to output the number (salsa_N_UFP) and mass concentration (salsa_PM0.1) of ultrafine particles, i.e. particles with a diameter smaller than 100 nm
  • Implement aerosol emission mode "parameterized" which is based on the street type (similar to the chemistry module).
  • Remove unnecessary nucleation subroutines.
  • Add the z-dimension for gaseous emissions to correspond the implementation in the chemistry module
Location:
palm/trunk
Files:
14 edited

Legend:

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

    r4039 r4131  
    2525! -----------------
    2626! $Id$
     27! Allow profile output for salsa variables.
     28!
     29! 4039 2019-06-18 10:32:41Z suehring
    2730! Correct conversion to kinematic scalar fluxes in case of pw-scheme and
    2831! statistic regions
     
    335338                simulated_time, simulated_time_at_begin,                       &
    336339                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
    337                 ws_scheme_mom, ws_scheme_sca
     340                ws_scheme_mom, ws_scheme_sca, salsa, max_pr_salsa
    338341
    339342    USE cpulog,                                                                &
     
    19541957                ENDIF
    19551958             ENDIF
     1959             IF ( salsa )  THEN
     1960                IF ( max_pr_cs > 0 )  THEN
     1961                   sums_l(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0) =    &
     1962                      sums_l(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0) + &
     1963                      sums_l(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,i)
     1964
     1965                ENDIF
     1966             ENDIF
    19561967          ENDDO
    19571968       ENDIF
     
    19741985                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),       &
    19751986                                    sums(nzb,pr_palm+max_pr_user+i),           &
     1987                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
     1988             ENDDO
     1989       ENDIF
     1990
     1991       IF ( salsa  .AND.  max_pr_salsa > 0 )  THEN
     1992          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1993             DO  i = 1, max_pr_salsa
     1994                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+max_pr_cs+i,0),                 &
     1995                                    sums(nzb,pr_palm+max_pr_user+max_pr_user+i),                   &
    19761996                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
    19771997             ENDDO
     
    20512071                                 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) / &
    20522072                                 ngp_2dh_s_inner(k,sr)
     2073             ENDDO
     2074          ENDIF 
     2075       ENDIF
     2076
     2077       IF ( salsa ) THEN
     2078          IF ( max_pr_salsa > 0 )  THEN
     2079             DO k = nzb, nzt+1
     2080                sums(k,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa) = &
     2081                  sums(k,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa) &
     2082                  / ngp_2dh_s_inner(k,sr)
    20532083             ENDDO
    20542084          ENDIF 
     
    22022232          ENDIF
    22032233       ENDIF
     2234
     2235       IF ( salsa )  THEN
     2236          IF ( max_pr_salsa > 0 )  THEN    ! salsa profiles
     2237             hom(:,1,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa, sr) = &
     2238                  sums(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa)
     2239          ENDIF
     2240       ENDIF
    22042241!
    22052242!--    Determine the boundary layer height using two different schemes.
  • palm/trunk/SOURCE/init_3d_model.f90

    r4130 r4131  
    2525! -----------------
    2626! $Id$
     27! Allocate sums and sums_l to allow profile output for salsa variables.
     28!
     29! 4130 2019-08-01 13:04:13Z suehring
    2730! Effectively reduce 3D initialization to 1D initial profiles. This is because
    2831! 3D initialization produces structures in the w-component that are correlated
     
    765768              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
    766769              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
    767               sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs),                   &
    768               sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs,0:threads_per_task-1),      &
     770              sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa),      &
     771              sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:threads_per_task-1),   &
    769772              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
    770773              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions) )
  • palm/trunk/SOURCE/module_interface.f90

    r4127 r4131  
    397397        ONLY:  salsa_parin,                                                    &
    398398               salsa_check_parameters,                                         &
     399               salsa_check_data_output_pr,                                     &
    399400               salsa_check_data_output,                                        &
    400401               salsa_init_arrays,                                              &
     
    409410               salsa_data_output_2d,                                           &
    410411               salsa_data_output_3d,                                           &
     412               salsa_statistics,                                               &
    411413               salsa_rrd_local,                                                &
    412414               salsa_wrd_local
     
    783785    ENDIF
    784786
     787    IF ( unit == 'illegal'  .AND.  salsa )  THEN
     788       CALL salsa_check_data_output_pr( variable, var_count, unit, dopr_unit )
     789    ENDIF
     790
    785791    IF ( unit == 'illegal'  .AND.  user_module_enabled )  THEN
    786792       unit = '' ! ToDo: Seems like a hack. Find a general soultion!
     
    15261532    IF ( gust_module_enabled )  CALL gust_statistics( mode, sr, tn, dots_max )
    15271533    IF ( air_chemistry       )  CALL chem_statistics( mode, sr, tn )
     1534    IF ( salsa               )  CALL salsa_statistics( mode, sr, tn )
    15281535
    15291536    IF ( user_module_enabled )  CALL user_statistics( mode, sr, tn )
  • palm/trunk/SOURCE/modules.f90

    r4110 r4131  
    2525! -----------------
    2626! $Id$
     27! Add max_pr_salsa to control_parameters. Used in creating profile output for
     28! salsa.
     29!
     30! 4110 2019-07-22 17:05:21Z suehring
    2731! -advc_flags_1, advc_flags_2
    2832! +advc_flags_m, advc_flags_s
     
    12351239    INTEGER(iwp) ::  maximum_grid_level                !< number of grid levels that the multigrid solver is using
    12361240    INTEGER(iwp) ::  maximum_parallel_io_streams = -1  !< maximum number of parallel io streams that the underlying parallel file system allows, set with palmrun option -w, ENVPAR namelist parameter, provided by palmrun
     1241    INTEGER(iwp) ::  max_pr_salsa = 0                  !< number of salsa profiles (must not change within a job chain)
    12371242    INTEGER(iwp) ::  max_pr_user = 0                   !< number of user-defined profiles (must not change within a job chain)
    12381243    INTEGER(iwp) ::  max_pr_user_tmp = 0               !< number of user-defined profiles that is temporary stored to check it against max_pr_user in case of restarts
  • palm/trunk/SOURCE/parin.f90

    r4079 r4131  
    2525! -----------------
    2626! $Id$
     27! Allocate hom and hom_sum to allow profile output for salsa variables.
     28!
     29! 4079 2019-07-09 18:04:41Z suehring
    2730! +monotonic_limiter_z
    2831!
     
    10931096                       ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),         &
    10941097                       u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),             &
    1095                        hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions),  &
    1096                        hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions) )
     1098                    hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:statistic_regions),  &
     1099                    hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:statistic_regions) )
    10971100
    10981101             hom = 0.0_wp
  • palm/trunk/SOURCE/read_restart_data_mod.f90

    r4101 r4131  
    2525! -----------------
    2626! $Id$
     27! Allocate hom and hom_sum to allow profile output for salsa variables.
     28!
     29! 4101 2019-07-17 15:14:26Z gronemeier
    2730! remove old_dt
    2831!
     
    316319                    v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),           &
    317320                    ref_state(0:nz+1), s_init(0:nz+1), sa_init(0:nz+1),        &
    318                     hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions),     &
    319                     hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions) )
     321                    hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:statistic_regions),  &
     322                    hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:statistic_regions) )
    320323       ENDIF
    321324
  • palm/trunk/SOURCE/salsa_mod.f90

    r4118 r4131  
    2626! -----------------
    2727! $Id$
     28! - Add "salsa_" before each salsa output variable
     29! - Add a possibility to output the number (salsa_N_UFP) and mass concentration
     30!   (salsa_PM0.1) of ultrafine particles, i.e. particles with a diameter smaller
     31!   than 100 nm
     32! - Implement aerosol emission mode "parameterized" which is based on the street
     33!   type (similar to the chemistry module).
     34! - Remove unnecessary nucleation subroutines.
     35! - Add the z-dimension for gaseous emissions to correspond the implementation
     36!   in the chemistry module
     37!
     38! 4118 2019-07-25 16:11:45Z suehring
    2839! - When Dirichlet condition is applied in decycling, the boundary conditions are
    2940!   only set at the ghost points and not at the prognostic grid points as done
     
    211222               bc_lr, bc_lr_cyc, bc_ns, bc_ns_cyc, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
    212223               bc_radiation_s, coupling_char, debug_output, dt_3d, intermediate_timestep_count,    &
    213                intermediate_timestep_count_max, land_surface, message_string, monotonic_limiter_z, &
    214                plant_canopy, pt_surface, salsa, scalar_advec, surface_pressure,                    &
    215                time_since_reference_point, timestep_scheme, tsc, urban_surface, ws_scheme_sca
     224               intermediate_timestep_count_max, land_surface, max_pr_salsa, message_string,        &
     225               monotonic_limiter_z, plant_canopy, pt_surface, salsa, scalar_advec,                 &
     226               surface_pressure, time_since_reference_point, timestep_scheme, tsc, urban_surface,  &
     227               ws_scheme_sca
    216228
    217229    USE indices,                                                                                   &
     
    392404                                            !< 1 = read vertical profile from an input file
    393405    INTEGER(iwp) ::  lod_gas_emissions = 0  !< level of detail of the gaseous emission data
     406    INTEGER(iwp) ::  main_street_id = 0     !< lower bound of main street IDs (OpenStreetMaps) for parameterized mode
     407    INTEGER(iwp) ::  max_street_id = 0      !< upper bound of main street IDs (OpenStreetMaps) for parameterized mode
    394408    INTEGER(iwp) ::  nbins_aerosol = 1      !< total number of size bins
    395409    INTEGER(iwp) ::  ncc   = 1              !< number of chemical components used
     
    413427                                            !< 9 = homomolecular nucleation of H2SO4 and ORG
    414428                                            !<     + heteromolecular nucleation with H2SO4*ORG
     429    INTEGER(iwp) ::  salsa_pr_count = 0     !< counter for salsa variable profiles
     430    INTEGER(iwp) ::  side_street_id = 0     !< lower bound of side street IDs (OpenStreetMaps) for parameterized mode
    415431    INTEGER(iwp) ::  start_subrange_1a = 1  !< start index for bin subranges: subrange 1a
    416432    INTEGER(iwp) ::  start_subrange_2a = 1  !<                                subrange 2a
     
    419435    INTEGER(iwp), DIMENSION(nreg) ::  nbin = (/ 3, 7/)  !< Number of size bins per subrange: 1 & 2
    420436
    421     INTEGER(iwp), DIMENSION(ngases_salsa) ::  gas_index_chem = &
    422                                                  (/ 1, 1, 1, 1, 1/)  !< gas indices in chemistry_model_mod
    423                                                  !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV, 5 = OCSV
     437    INTEGER(iwp), DIMENSION(ngases_salsa) ::  gas_index_chem = (/ 1, 1, 1, 1, 1/)  !< gas indices in chemistry_model_mod
     438                                                                                   !< 1 = H2SO4, 2 = HNO3,
     439                                                                                   !< 3 = NH3,  4 = OCNV, 5 = OCSV
    424440    INTEGER(iwp), DIMENSION(ngases_salsa) ::  emission_index_chem  !< gas indices in the gas emission file
     441    INTEGER(iwp), DIMENSION(99) ::  salsa_pr_index  = 0            !< index for salsa profiles
    425442
    426443    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  k_topo_top  !< vertical index of the topography top
     
    466483    LOGICAL ::  lspartition  = .FALSE.  !< Partition of HNO3 and NH3
    467484
    468     REAL(wp) ::  act_coeff = 1.0E-7_wp               !< Activation coefficient
     485    REAL(wp) ::  act_coeff = 1.0E-7_wp               !< Activation coefficient (1/s)
    469486    REAL(wp) ::  dt_salsa  = 0.00001_wp              !< Time step of SALSA
     487    REAL(wp) ::  emiss_factor_main = 0.0_wp          !< relative emission factor for main streets
     488    REAL(wp) ::  emiss_factor_side = 0.0_wp          !< relative emission factor for side streets
    470489    REAL(wp) ::  h2so4_init = nclim                  !< Init value for sulphuric acid gas
    471490    REAL(wp) ::  hno3_init  = nclim                  !< Init value for nitric acid gas
     
    557576       INTEGER(iwp) ::  tind               !< time index for reference time in salsa emission data
    558577
    559        INTEGER(iwp), DIMENSION(maxspec) ::  cc_input_to_model   !<
     578       INTEGER(iwp), DIMENSION(maxspec) ::  cc_in2mod = 0   !<
    560579
    561580       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  cat_index  !< Index of emission categories
     
    630649    TYPE salsa_variable
    631650
    632        REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  init  !<
    633 
    634        REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  diss_s     !<
    635        REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  flux_s     !<
    636        REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  source     !<
    637        REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  sums_ws_l  !<
    638 
    639        REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  diss_l  !<
    640        REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  flux_l  !<
    641 
    642        REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  conc     !<
    643        REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  conc_p   !<
    644        REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS ::  tconc_m  !<
     651       REAL(wp), DIMENSION(:), ALLOCATABLE     ::  init  !<
     652
     653       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s     !<
     654       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s     !<
     655       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  source     !<
     656       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_ws_l  !<
     657
     658       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l  !<
     659       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l  !<
     660
     661       REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  conc     !<
     662       REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  conc_p   !<
     663       REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tconc_m  !<
    645664
    646665    END TYPE salsa_variable
     
    673692
    674693    TYPE(chem_emis_att_type) ::  chem_emission_att  !< chemistry emission attributes
    675     TYPE(chem_emis_val_type) ::  chem_emission      !< chemistry emission values
     694
     695    TYPE(chem_emis_val_type), DIMENSION(:), ALLOCATABLE ::  chem_emission  !< chemistry emissions
    676696
    677697    TYPE(t_section), DIMENSION(:), ALLOCATABLE ::  aero  !< local aerosol properties
     
    683703    TYPE(match_surface), DIMENSION(0:3) ::  usm_to_depo_v  !< to match the deposition mod. and vertical USM surfaces
    684704!
    685 !-- SALSA variables: as x = x(k,j,i,bin). 
     705!-- SALSA variables: as x = x(k,j,i,bin).
    686706!-- The 4th dimension contains all the size bins sequentially for each aerosol species  + water.
    687707!
     
    689709!
    690710!-- Number concentration (#/m3)
    691     TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_number  !<
    692     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_1  !<
    693     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_2  !<
    694     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nconc_3  !<
     711    TYPE(salsa_variable), DIMENSION(:), ALLOCATABLE, TARGET ::  aerosol_number  !<
     712    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  nconc_1  !<
     713    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  nconc_2  !<
     714    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  nconc_3  !<
    695715!
    696716!-- Mass concentration (kg/m3)
    697     TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  aerosol_mass  !<
    698     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_1  !<
    699     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_2  !<
    700     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mconc_3  !<
     717    TYPE(salsa_variable), DIMENSION(:), ALLOCATABLE, TARGET ::  aerosol_mass  !<
     718    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  mconc_1  !<
     719    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  mconc_2  !<
     720    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  mconc_3  !<
    701721!
    702722!-- Gaseous concentrations (#/m3)
    703     TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET ::  salsa_gas  !<
    704     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_1  !<
    705     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_2  !<
    706     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  gconc_3  !<
     723    TYPE(salsa_variable), DIMENSION(:), ALLOCATABLE, TARGET ::  salsa_gas  !<
     724    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  gconc_1  !<
     725    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  gconc_2  !<
     726    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  gconc_3  !<
    707727!
    708728!-- Diagnostic tracers
    709     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  sedim_vd  !< sedimentation velocity per bin (m/s)
    710     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  ra_dry    !< aerosol dry radius (m)
     729    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  sedim_vd  !< sedimentation velocity per bin (m/s)
     730    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  ra_dry    !< aerosol dry radius (m)
    711731
    712732!-- Particle component index tables
     
    717737!
    718738!-- Gases:
    719     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_h2so4_av  !< H2SO4
    720     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_hno3_av   !< HNO3
    721     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_nh3_av    !< NH3
    722     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_ocnv_av   !< non-volatile OC
    723     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  g_ocsv_av   !< semi-volatile OC
     739    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_h2so4_av  !< H2SO4
     740    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_hno3_av   !< HNO3
     741    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_nh3_av    !< NH3
     742    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_ocnv_av   !< non-volatile OC
     743    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  g_ocsv_av   !< semi-volatile OC
    724744!
    725745!-- Integrated:
    726     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ldsa_av  !< lung-deposited surface area
    727     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ntot_av  !< total number concentration
    728     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  pm25_av  !< PM2.5
    729     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  pm10_av  !< PM10
     746    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ldsa_av  !< lung-deposited surface area
     747    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ntot_av  !< total number concentration
     748    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nufp_av  !< ultrafine particles (UFP)
     749    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  pm01_av  !< PM0.1
     750    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  pm25_av  !< PM2.5
     751    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  pm10_av  !< PM10
    730752!
    731753!-- In the particle phase:
    732     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_bc_av   !< black carbon
    733     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_du_av   !< dust
    734     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_h2o_av  !< liquid water
    735     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_nh_av   !< ammonia
    736     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_no_av   !< nitrates
    737     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_oc_av   !< org. carbon
    738     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_so4_av  !< sulphates
    739     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  s_ss_av   !< sea salt
     754    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_bc_av   !< black carbon
     755    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_du_av   !< dust
     756    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_h2o_av  !< liquid water
     757    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_nh_av   !< ammonia
     758    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_no_av   !< nitrates
     759    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_oc_av   !< org. carbon
     760    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_so4_av  !< sulphates
     761    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  s_ss_av   !< sea salt
    740762!
    741763!-- Bin specific mass and number concentrations:
    742     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  mbins_av  !< bin mas
    743     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  nbins_av  !< bin number
     764    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  mbins_av  !< bin mas
     765    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  nbins_av  !< bin number
    744766
    745767!
    746768!-- PALM interfaces:
    747 !
    748 !-- Boundary conditions:
     769
     770    INTERFACE salsa_actions
     771       MODULE PROCEDURE salsa_actions
     772       MODULE PROCEDURE salsa_actions_ij
     773    END INTERFACE salsa_actions
     774
     775    INTERFACE salsa_3d_data_averaging
     776       MODULE PROCEDURE salsa_3d_data_averaging
     777    END INTERFACE salsa_3d_data_averaging
     778
    749779    INTERFACE salsa_boundary_conds
    750780       MODULE PROCEDURE salsa_boundary_conds
    751781       MODULE PROCEDURE salsa_boundary_conds_decycle
    752782    END INTERFACE salsa_boundary_conds
    753 !
    754 !-- Data output checks for 2D/3D data to be done in check_parameters
     783
    755784    INTERFACE salsa_check_data_output
    756785       MODULE PROCEDURE salsa_check_data_output
    757786    END INTERFACE salsa_check_data_output
    758 !
    759 !-- Input parameter checks to be done in check_parameters
     787
     788    INTERFACE salsa_check_data_output_pr
     789       MODULE PROCEDURE salsa_check_data_output_pr
     790    END INTERFACE salsa_check_data_output_pr
     791
    760792    INTERFACE salsa_check_parameters
    761793       MODULE PROCEDURE salsa_check_parameters
    762794    END INTERFACE salsa_check_parameters
    763 !
    764 !-- Averaging of 3D data for output
    765     INTERFACE salsa_3d_data_averaging
    766        MODULE PROCEDURE salsa_3d_data_averaging
    767     END INTERFACE salsa_3d_data_averaging
    768 !
    769 !-- Data output of 2D quantities
     795
    770796    INTERFACE salsa_data_output_2d
    771797       MODULE PROCEDURE salsa_data_output_2d
    772798    END INTERFACE salsa_data_output_2d
    773 !
    774 !-- Data output of 3D data
     799
    775800    INTERFACE salsa_data_output_3d
    776801       MODULE PROCEDURE salsa_data_output_3d
    777802    END INTERFACE salsa_data_output_3d
    778 !
    779 !-- Data output of 3D data
     803
    780804    INTERFACE salsa_data_output_mask
    781805       MODULE PROCEDURE salsa_data_output_mask
    782806    END INTERFACE salsa_data_output_mask
    783 !
    784 !-- Definition of data output quantities
     807
    785808    INTERFACE salsa_define_netcdf_grid
    786809       MODULE PROCEDURE salsa_define_netcdf_grid
    787810    END INTERFACE salsa_define_netcdf_grid
    788 !
    789 !-- Output of information to the header file
     811
     812    INTERFACE salsa_driver
     813       MODULE PROCEDURE salsa_driver
     814    END INTERFACE salsa_driver
     815
     816    INTERFACE salsa_emission_update
     817       MODULE PROCEDURE salsa_emission_update
     818    END INTERFACE salsa_emission_update
     819
     820    INTERFACE salsa_exchange_horiz_bounds
     821       MODULE PROCEDURE salsa_exchange_horiz_bounds
     822    END INTERFACE salsa_exchange_horiz_bounds
     823
    790824    INTERFACE salsa_header
    791825       MODULE PROCEDURE salsa_header
    792826    END INTERFACE salsa_header
    793 !
    794 !-- Initialization actions
     827
    795828    INTERFACE salsa_init
    796829       MODULE PROCEDURE salsa_init
    797830    END INTERFACE salsa_init
    798 !
    799 !-- Initialization of arrays
     831
    800832    INTERFACE salsa_init_arrays
    801833       MODULE PROCEDURE salsa_init_arrays
    802834    END INTERFACE salsa_init_arrays
    803 !
    804 !-- Writing of binary output for restart runs  !!! renaming?!
    805     INTERFACE salsa_wrd_local
    806        MODULE PROCEDURE salsa_wrd_local
    807     END INTERFACE salsa_wrd_local
    808 !
    809 !-- Reading of NAMELIST parameters
    810     INTERFACE salsa_parin
    811        MODULE PROCEDURE salsa_parin
    812     END INTERFACE salsa_parin
    813 !
    814 !-- Reading of parameters for restart runs
    815     INTERFACE salsa_rrd_local
    816        MODULE PROCEDURE salsa_rrd_local
    817     END INTERFACE salsa_rrd_local
    818 !
    819 !-- Swapping of time levels (required for prognostic variables)
    820     INTERFACE salsa_swap_timelevel
    821        MODULE PROCEDURE salsa_swap_timelevel
    822     END INTERFACE salsa_swap_timelevel
    823 !
    824 !-- Interface between PALM and salsa
    825     INTERFACE salsa_driver
    826        MODULE PROCEDURE salsa_driver
    827     END INTERFACE salsa_driver
    828 
    829 !-- Actions salsa variables
    830     INTERFACE salsa_actions
    831        MODULE PROCEDURE salsa_actions
    832        MODULE PROCEDURE salsa_actions_ij
    833     END INTERFACE salsa_actions
    834 !
    835 !-- Non-advective processes (i.e. aerosol dynamic reactions) for salsa variables
     835
    836836    INTERFACE salsa_non_advective_processes
    837837       MODULE PROCEDURE salsa_non_advective_processes
    838838       MODULE PROCEDURE salsa_non_advective_processes_ij
    839839    END INTERFACE salsa_non_advective_processes
    840 !
    841 !-- Exchange horiz for salsa variables
    842     INTERFACE salsa_exchange_horiz_bounds
    843        MODULE PROCEDURE salsa_exchange_horiz_bounds
    844     END INTERFACE salsa_exchange_horiz_bounds
    845 !
    846 !-- Prognostics equations for salsa variables
     840
     841    INTERFACE salsa_parin
     842       MODULE PROCEDURE salsa_parin
     843    END INTERFACE salsa_parin
     844
    847845    INTERFACE salsa_prognostic_equations
    848846       MODULE PROCEDURE salsa_prognostic_equations
    849847       MODULE PROCEDURE salsa_prognostic_equations_ij
    850848    END INTERFACE salsa_prognostic_equations
    851 !
    852 !-- Tendency salsa variables
     849
     850    INTERFACE salsa_rrd_local
     851       MODULE PROCEDURE salsa_rrd_local
     852    END INTERFACE salsa_rrd_local
     853
     854    INTERFACE salsa_statistics
     855       MODULE PROCEDURE salsa_statistics
     856    END INTERFACE salsa_statistics
     857
     858    INTERFACE salsa_swap_timelevel
     859       MODULE PROCEDURE salsa_swap_timelevel
     860    END INTERFACE salsa_swap_timelevel
     861
    853862    INTERFACE salsa_tendency
    854863       MODULE PROCEDURE salsa_tendency
    855864       MODULE PROCEDURE salsa_tendency_ij
    856865    END INTERFACE salsa_tendency
     866
     867    INTERFACE salsa_wrd_local
     868       MODULE PROCEDURE salsa_wrd_local
     869    END INTERFACE salsa_wrd_local
    857870
    858871
     
    867880           salsa_emission_update, salsa_header, salsa_init, salsa_init_arrays, salsa_parin,        &
    868881           salsa_rrd_local, salsa_swap_timelevel, salsa_prognostic_equations, salsa_wrd_local,     &
    869            salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds
     882           salsa_actions, salsa_non_advective_processes, salsa_exchange_horiz_bounds,              &
     883           salsa_check_data_output_pr, salsa_statistics
    870884!
    871885!-- Public parameters, constants and initial values
     
    888902 SUBROUTINE salsa_parin
    889903
     904    USE control_parameters,                                                                        &
     905        ONLY:  data_output_pr
     906
    890907    IMPLICIT NONE
    891908
    892     CHARACTER(LEN=80) ::  line   !< dummy string that contains the current line
    893                                   !< of the parameter file
    894 
    895     NAMELIST /salsa_parameters/      aerosol_flux_dpg, aerosol_flux_mass_fracs_a,                  &
    896                                      aerosol_flux_mass_fracs_b, aerosol_flux_sigmag,               &
    897                                      advect_particle_water, bc_salsa_b, bc_salsa_t,                &
    898                                      decycle_salsa_lr, decycle_method_salsa, decycle_salsa_ns,     &
    899                                      depo_pcm_par, depo_pcm_type,                                  &
    900                                      depo_surf_par, dpg, dt_salsa, feedback_to_palm, h2so4_init,   &
    901                                      hno3_init, init_gases_type, init_aerosol_type, listspec,      &
    902                                      mass_fracs_a, mass_fracs_b, n_lognorm, nbin, nest_salsa, nf2a,&
    903                                      nh3_init, nj3, nlcnd, nlcndgas, nlcndh2oae, nlcoag, nldepo,   &
    904                                      nldepo_pcm,  nldepo_surf, nldistupdate, nsnucl, ocnv_init,    &
    905                                      ocsv_init, read_restart_data_salsa, reglim, salsa,            &
    906                                      salsa_emission_mode, sigmag, skip_time_do_salsa,              &
    907                                      surface_aerosol_flux, van_der_waals_coagc, write_binary_salsa
     909    CHARACTER(LEN=80) ::  line   !< dummy string that contains the current line of parameter file
     910
     911    INTEGER(iwp) ::  i                 !< loop index
     912    INTEGER(iwp) ::  max_pr_salsa_tmp  !< dummy variable
     913
     914    NAMELIST /salsa_parameters/      aerosol_flux_dpg,                         &
     915                                     aerosol_flux_mass_fracs_a,                &
     916                                     aerosol_flux_mass_fracs_b,                &
     917                                     aerosol_flux_sigmag,                      &
     918                                     advect_particle_water,                    &
     919                                     bc_salsa_b,                               &
     920                                     bc_salsa_t,                               &
     921                                     decycle_salsa_lr,                         &
     922                                     decycle_method_salsa,                     &
     923                                     decycle_salsa_ns,                         &
     924                                     depo_pcm_par,                             &
     925                                     depo_pcm_type,                            &
     926                                     depo_surf_par,                            &
     927                                     dpg,                                      &
     928                                     dt_salsa,                                 &
     929                                     emiss_factor_main,                        &
     930                                     emiss_factor_side,                        &
     931                                     feedback_to_palm,                         &
     932                                     h2so4_init,                               &
     933                                     hno3_init,                                &
     934                                     init_gases_type,                          &
     935                                     init_aerosol_type,                        &
     936                                     listspec,                                 &
     937                                     main_street_id,                           &
     938                                     mass_fracs_a,                             &
     939                                     mass_fracs_b,                             &
     940                                     max_street_id,                            &
     941                                     n_lognorm,                                &
     942                                     nbin,                                     &
     943                                     nest_salsa,                               &
     944                                     nf2a,                                     &
     945                                     nh3_init,                                 &
     946                                     nj3,                                      &
     947                                     nlcnd,                                    &
     948                                     nlcndgas,                                 &
     949                                     nlcndh2oae,                               &
     950                                     nlcoag,                                   &
     951                                     nldepo,                                   &
     952                                     nldepo_pcm,                               &
     953                                     nldepo_surf,                              &
     954                                     nldistupdate,                             &
     955                                     nsnucl,                                   &
     956                                     ocnv_init,                                &
     957                                     ocsv_init,                                &
     958                                     read_restart_data_salsa,                  &
     959                                     reglim,                                   &
     960                                     salsa,                                    &
     961                                     salsa_emission_mode,                      &
     962                                     sigmag,                                   &
     963                                     side_street_id,                           &
     964                                     skip_time_do_salsa,                       &
     965                                     surface_aerosol_flux,                     &
     966                                     van_der_waals_coagc,                      &
     967                                     write_binary_salsa
    908968
    909969    line = ' '
     
    924984
    925985 10 CONTINUE
     986!
     987!-- Update the number of output profiles
     988    max_pr_salsa_tmp = 0
     989    i = 1
     990    DO WHILE ( data_output_pr(i) /= ' '  .AND.  i <= 100 )
     991       IF ( TRIM( data_output_pr(i)(1:6) ) == 'salsa_' )  max_pr_salsa_tmp = max_pr_salsa_tmp + 1
     992       i = i + 1
     993    ENDDO
     994    IF ( max_pr_salsa_tmp > 0 )  max_pr_salsa = max_pr_salsa_tmp
    926995
    927996 END SUBROUTINE salsa_parin
     
    10031072!-- Check for the grid
    10041073
    1005     IF ( var(1:2) == 'g_' )  THEN
    1006        grid_x = 'x'
    1007        grid_y = 'y'
    1008        grid_z = 'zu'
    1009     ELSEIF ( var(1:4) == 'LDSA' )  THEN
    1010        grid_x = 'x'
    1011        grid_y = 'y'
    1012        grid_z = 'zu'
    1013     ELSEIF ( var(1:5) == 'm_bin' )  THEN
    1014        grid_x = 'x'
    1015        grid_y = 'y'
    1016        grid_z = 'zu'
    1017     ELSEIF ( var(1:5) == 'N_bin' )  THEN
    1018        grid_x = 'x'
    1019        grid_y = 'y'
    1020        grid_z = 'zu'
    1021     ELSEIF ( var(1:4) == 'Ntot' ) THEN
    1022        grid_x = 'x'
    1023        grid_y = 'y'
    1024        grid_z = 'zu'
    1025     ELSEIF ( var(1:2) == 'PM' )  THEN
    1026        grid_x = 'x'
    1027        grid_y = 'y'
    1028        grid_z = 'zu'
    1029     ELSEIF ( var(1:2) == 's_' )  THEN
     1074    IF ( var(1:6) == 'salsa_' )  THEN  ! same grid for all salsa output variables
    10301075       grid_x = 'x'
    10311076       grid_y = 'y'
     
    11631208    INTEGER(iwp) ::  gases_available !< Number of available gas components in the chemistry model
    11641209    INTEGER(iwp) ::  i               !< loop index for allocating
     1210    INTEGER(iwp) ::  ii              !< index for indexing chemical components
    11651211    INTEGER(iwp) ::  l               !< loop index for allocating: surfaces
    11661212    INTEGER(iwp) ::  lsp             !< loop index for chem species in the chemistry model
     
    11941240
    11951241    ncomponents_mass = ncc
    1196     IF ( advect_particle_water )  ncomponents_mass = ncc + 1  ! Add water
    1197 
    1198 !
    1199 !-- Allocate:
    1200     ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass),                    &
    1201               bc_an_t_val(nbins_aerosol), bc_gt_t_val(ngases_salsa), bin_low_limits(nbins_aerosol),&
    1202               nsect(nbins_aerosol), massacc(nbins_aerosol) )
    1203     ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) )
    1204     IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
    1205     ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
    1206 
    1207 !
    1208 !-- Aerosol number concentration
    1209     ALLOCATE( aerosol_number(nbins_aerosol) )
    1210     ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
    1211               nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
    1212               nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
    1213     nconc_1 = 0.0_wp
    1214     nconc_2 = 0.0_wp
    1215     nconc_3 = 0.0_wp
    1216 
    1217     DO i = 1, nbins_aerosol
    1218        aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => nconc_1(:,:,:,i)
    1219        aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => nconc_2(:,:,:,i)
    1220        aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i)
    1221        ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                         &
    1222                  aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                         &
    1223                  aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
    1224                  aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
    1225                  aerosol_number(i)%init(nzb:nzt+1),                                                &
    1226                  aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    1227        aerosol_number(i)%init = nclim
    1228        IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
    1229           ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) )
    1230        ENDIF
    1231     ENDDO
    1232 
    1233 !
    1234 !-- Aerosol mass concentration
    1235     ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) )
    1236     ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
    1237               mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
    1238               mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) )
    1239     mconc_1 = 0.0_wp
    1240     mconc_2 = 0.0_wp
    1241     mconc_3 = 0.0_wp
    1242 
    1243     DO i = 1, ncomponents_mass*nbins_aerosol
    1244        aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => mconc_1(:,:,:,i)
    1245        aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => mconc_2(:,:,:,i)
    1246        aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i)
    1247        ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                           &
    1248                  aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                           &
    1249                  aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
    1250                  aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
    1251                  aerosol_mass(i)%init(nzb:nzt+1),                                                  &
    1252                  aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
    1253        aerosol_mass(i)%init = mclim
    1254        IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
    1255           ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) )
    1256        ENDIF
    1257     ENDDO
    1258 
    1259 !
    1260 !-- Surface fluxes: answs = aerosol number, amsws = aerosol mass
    1261 !
    1262 !-- Horizontal surfaces: default type
    1263     DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
    1264        ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) )
    1265        ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) )
    1266        surf_def_h(l)%answs = 0.0_wp
    1267        surf_def_h(l)%amsws = 0.0_wp
    1268     ENDDO
    1269 !
    1270 !-- Horizontal surfaces: natural type
    1271     ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) )
    1272     ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) )
    1273     surf_lsm_h%answs = 0.0_wp
    1274     surf_lsm_h%amsws = 0.0_wp
    1275 !
    1276 !-- Horizontal surfaces: urban type
    1277     ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) )
    1278     ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) )
    1279     surf_usm_h%answs = 0.0_wp
    1280     surf_usm_h%amsws = 0.0_wp
    1281 
    1282 !
    1283 !-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing
    1284     DO  l = 0, 3
    1285        ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) )
    1286        surf_def_v(l)%answs = 0.0_wp
    1287        ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
    1288        surf_def_v(l)%amsws = 0.0_wp
    1289 
    1290        ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) )
    1291        surf_lsm_v(l)%answs = 0.0_wp
    1292        ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
    1293        surf_lsm_v(l)%amsws = 0.0_wp
    1294 
    1295        ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) )
    1296        surf_usm_v(l)%answs = 0.0_wp
    1297        ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
    1298        surf_usm_v(l)%amsws = 0.0_wp
    1299 
    1300     ENDDO
    1301 
    1302 !
    1303 !-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV)
    1304 !-- (number concentration (#/m3) )
    1305 !
    1306 !-- If chemistry is on, read gas phase concentrations from there. Otherwise,
    1307 !-- allocate salsa_gas array.
    1308 
    1309     IF ( air_chemistry )  THEN
    1310        DO  lsp = 1, nvar
    1311           SELECT CASE ( TRIM( chem_species(lsp)%name ) )
    1312              CASE ( 'H2SO4', 'h2so4' )
    1313                 gases_available = gases_available + 1
    1314                 gas_index_chem(1) = lsp
    1315              CASE ( 'HNO3', 'hno3' )
    1316                 gases_available = gases_available + 1
    1317                 gas_index_chem(2) = lsp
    1318              CASE ( 'NH3', 'nh3' )
    1319                 gases_available = gases_available + 1
    1320                 gas_index_chem(3) = lsp
    1321              CASE ( 'OCNV', 'ocnv' )
    1322                 gases_available = gases_available + 1
    1323                 gas_index_chem(4) = lsp
    1324              CASE ( 'OCSV', 'ocsv' )
    1325                 gases_available = gases_available + 1
    1326                 gas_index_chem(5) = lsp
    1327           END SELECT
    1328        ENDDO
    1329 
    1330        IF ( gases_available == ngases_salsa )  THEN
    1331           salsa_gases_from_chem = .TRUE.
    1332        ELSE
    1333           WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// &
    1334                                      'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)'
    1335        CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 )
    1336        ENDIF
    1337 
    1338     ELSE
    1339 
    1340        ALLOCATE( salsa_gas(ngases_salsa) )
    1341        ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
    1342                  gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
    1343                  gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
    1344        gconc_1 = 0.0_wp
    1345        gconc_2 = 0.0_wp
    1346        gconc_3 = 0.0_wp
    1347 
    1348        DO i = 1, ngases_salsa
    1349           salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => gconc_1(:,:,:,i)
    1350           salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => gconc_2(:,:,:,i)
    1351           salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i)
    1352           ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),       &
    1353                     salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),       &
    1354                     salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
    1355                     salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
    1356                     salsa_gas(i)%init(nzb:nzt+1),                              &
    1357                     salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    1358           salsa_gas(i)%init = nclim
    1359           IF ( include_emission )  ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) )
    1360        ENDDO
    1361 !
    1362 !--    Surface fluxes: gtsws = gaseous tracer flux
    1363 !
    1364 !--    Horizontal surfaces: default type
    1365        DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
    1366           ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) )
    1367           surf_def_h(l)%gtsws = 0.0_wp
    1368        ENDDO
    1369 !--    Horizontal surfaces: natural type
    1370        ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) )
    1371        surf_lsm_h%gtsws = 0.0_wp
    1372 !--    Horizontal surfaces: urban type
    1373        ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) )
    1374        surf_usm_h%gtsws = 0.0_wp
    1375 !
    1376 !--    Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and
    1377 !--    westward (l=3) facing
    1378        DO  l = 0, 3
    1379           ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) )
    1380           surf_def_v(l)%gtsws = 0.0_wp
    1381           ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) )
    1382           surf_lsm_v(l)%gtsws = 0.0_wp
    1383           ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) )
    1384           surf_usm_v(l)%gtsws = 0.0_wp
    1385        ENDDO
    1386     ENDIF
    1387 
    1388     IF ( ws_scheme_sca )  THEN
    1389 
    1390        IF ( salsa )  THEN
    1391           ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    1392           sums_salsa_ws_l = 0.0_wp
    1393        ENDIF
    1394 
    1395     ENDIF
    1396 !
    1397 !-- Set control flags for decycling only at lateral boundary cores. Within the inner cores the
    1398 !-- decycle flag is set to .FALSE.. Even though it does not affect the setting of chemistry boundary
    1399 !-- conditions, this flag is used to set advection control flags appropriately.
    1400     decycle_salsa_lr = MERGE( decycle_salsa_lr, .FALSE., nxl == 0  .OR.  nxr == nx )
    1401     decycle_salsa_ns = MERGE( decycle_salsa_ns, .FALSE., nys == 0  .OR.  nyn == ny )
    1402 !
    1403 !-- Decycling can be applied separately for aerosol variables, while wind and other scalars may have
    1404 !-- cyclic or nested boundary conditions. However, large gradients near the boundaries may produce
    1405 !-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
    1406 !-- applied near these boundaries. To get rid-off this, set-up additional flags that control the
    1407 !-- order of the scalar advection scheme near the lateral boundaries for passive scalars with
    1408 !-- decycling.
    1409     IF ( scalar_advec == 'ws-scheme' )  THEN
    1410        ALLOCATE( salsa_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    1411 !
    1412 !--    In case of decycling, set Neuman boundary conditions for wall_flags_0 bit 31 instead of
    1413 !--    cyclic boundary conditions. Bit 31 is used to identify extended degradation zones (please see
    1414 !--    the following comment). Note, since several also other modules may access this bit but may
    1415 !--    have other boundary conditions, the original value of wall_flags_0 bit 31 must not be
    1416 !--    modified. Hence, store the boundary conditions directly on salsa_advc_flags_s.
    1417 !--    salsa_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used
    1418 !--    to control the numerical order.
    1419 !--    Initialize with flag 31 only.
    1420        salsa_advc_flags_s = 0
    1421        salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_0, 31 ) )
    1422 
    1423        IF ( decycle_salsa_ns )  THEN
    1424           IF ( nys == 0 )  THEN
    1425              DO  i = 1, nbgp
    1426                 salsa_advc_flags_s(:,nys-i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nys,:), 31 ),   &
    1427                                                        IBCLR( salsa_advc_flags_s(:,nys,:), 31 ),   &
    1428                                                        BTEST( salsa_advc_flags_s(:,nys,:), 31 ) )
    1429              ENDDO
    1430           ENDIF
    1431           IF ( nyn == ny )  THEN
    1432              DO  i = 1, nbgp
    1433                 salsa_advc_flags_s(:,nyn+i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nyn,:), 31 ),   &
    1434                                                        IBCLR( salsa_advc_flags_s(:,nyn,:), 31 ),   &
    1435                                                        BTEST( salsa_advc_flags_s(:,nyn,:), 31 ) )
    1436              ENDDO
    1437           ENDIF
    1438        ENDIF
    1439        IF ( decycle_salsa_lr )  THEN
    1440           IF ( nxl == 0 )  THEN
    1441              DO  i = 1, nbgp
    1442                 salsa_advc_flags_s(:,:,nxl-i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxl), 31 ),   &
    1443                                                        IBCLR( salsa_advc_flags_s(:,:,nxl), 31 ),   &
    1444                                                        BTEST( salsa_advc_flags_s(:,:,nxl), 31 ) )
    1445              ENDDO
    1446           ENDIF
    1447           IF ( nxr == nx )  THEN
    1448              DO  i = 1, nbgp
    1449                 salsa_advc_flags_s(:,:,nxr+i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxr), 31 ),   &
    1450                                                        IBCLR( salsa_advc_flags_s(:,:,nxr), 31 ),   &
    1451                                                        BTEST( salsa_advc_flags_s(:,:,nxr), 31 ) )
    1452              ENDDO
    1453           ENDIF
    1454        ENDIF
    1455 !
    1456 !--    To initialise the advection flags appropriately, pass the boundary flags to
    1457 !--    ws_init_flags_scalar. The last argument in ws_init_flags_scalar indicates that a passive
    1458 !--    scalar is being treated and the horizontal advection terms are degraded already 2 grid points
    1459 !--    before the lateral boundary. Also, extended degradation zones are applied, where
    1460 !--    horizontal advection of scalars is discretised by the first-order scheme at all grid points
    1461 !--    in the vicinity of buildings (<= 3 grid points). Even though no building is within the
    1462 !--    numerical stencil, the first-order scheme is used. At fourth and fifth grid points, the order
    1463 !--    of the horizontal advection scheme is successively upgraded.
    1464 !--    These degradations of the advection scheme are done to avoid stationary numerical
    1465 !--    oscillations, which are responsible for high concentration maxima that may appear e.g. under
    1466 !--    shear-free stable conditions.
    1467        CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_salsa_lr,    &
    1468                                   bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_salsa_ns,    &
    1469                                   bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_salsa_lr,    &
    1470                                   bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_salsa_ns,    &
    1471                                   salsa_advc_flags_s, .TRUE. )
    1472     ENDIF
    1473 
    1474 
    1475  END SUBROUTINE salsa_init_arrays
    1476 
    1477 !------------------------------------------------------------------------------!
    1478 ! Description:
    1479 ! ------------
    1480 !> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA.
    1481 !> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are
    1482 !> also merged here.
    1483 !------------------------------------------------------------------------------!
    1484  SUBROUTINE salsa_init
    1485 
    1486     IMPLICIT NONE
    1487 
    1488     INTEGER(iwp) :: i   !<
    1489     INTEGER(iwp) :: ib  !< loop index for aerosol number bins
    1490     INTEGER(iwp) :: ic  !< loop index for aerosol mass bins
    1491     INTEGER(iwp) :: ig  !< loop index for gases
    1492     INTEGER(iwp) :: ii  !< index for indexing
    1493     INTEGER(iwp) :: j   !<
    1494 
    1495     IF ( debug_output )  CALL debug_message( 'salsa_init', 'start' )
    1496 
    1497     bin_low_limits = 0.0_wp
    1498     k_topo_top     = 0
    1499     nsect          = 0.0_wp
    1500     massacc        = 1.0_wp
    1501 
     1242    IF ( advect_particle_water )  ncomponents_mass = ncc + 1  ! Add water
    15021243!
    15031244!-- Indices for chemical components used (-1 = not used)
     
    15381279    ENDIF
    15391280!
    1540 !-- Partition and dissolutional growth by gaseous HNO3 and NH3
    1541     IF ( index_no > 0  .AND.  index_nh > 0  .AND.  index_so4 > 0 )  lspartition = .TRUE.
     1281!-- Allocate:
     1282    ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass),                    &
     1283              bc_an_t_val(nbins_aerosol), bc_gt_t_val(ngases_salsa), bin_low_limits(nbins_aerosol),&
     1284              nsect(nbins_aerosol), massacc(nbins_aerosol) )
     1285    ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) )
     1286    IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     1287    ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     1288!
     1289!-- Initialise the sectional particle size distribution
     1290    CALL set_sizebins
     1291!
     1292!-- Aerosol number concentration
     1293    ALLOCATE( aerosol_number(nbins_aerosol) )
     1294    ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
     1295              nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol),                                &
     1296              nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     1297    nconc_1 = 0.0_wp
     1298    nconc_2 = 0.0_wp
     1299    nconc_3 = 0.0_wp
     1300
     1301    DO i = 1, nbins_aerosol
     1302       aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => nconc_1(:,:,:,i)
     1303       aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => nconc_2(:,:,:,i)
     1304       aerosol_number(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_3(:,:,:,i)
     1305       ALLOCATE( aerosol_number(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                         &
     1306                 aerosol_number(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                         &
     1307                 aerosol_number(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
     1308                 aerosol_number(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                 &
     1309                 aerosol_number(i)%init(nzb:nzt+1),                                                &
     1310                 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     1311       aerosol_number(i)%init = nclim
     1312       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
     1313          ALLOCATE( aerosol_number(i)%source(nys:nyn,nxl:nxr) )
     1314          aerosol_number(i)%source = 0.0_wp
     1315       ENDIF
     1316    ENDDO
     1317
     1318!
     1319!-- Aerosol mass concentration
     1320    ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) )
     1321    ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
     1322              mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol),               &
     1323              mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) )
     1324    mconc_1 = 0.0_wp
     1325    mconc_2 = 0.0_wp
     1326    mconc_3 = 0.0_wp
     1327
     1328    DO i = 1, ncomponents_mass*nbins_aerosol
     1329       aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => mconc_1(:,:,:,i)
     1330       aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => mconc_2(:,:,:,i)
     1331       aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i)
     1332       ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),                           &
     1333                 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),                           &
     1334                 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
     1335                 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                   &
     1336                 aerosol_mass(i)%init(nzb:nzt+1),                                                  &
     1337                 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1)  )
     1338       aerosol_mass(i)%init = mclim
     1339       IF ( include_emission  .OR.  ( nldepo  .AND.  nldepo_surf ) )  THEN
     1340          ALLOCATE( aerosol_mass(i)%source(nys:nyn,nxl:nxr) )
     1341          aerosol_mass(i)%source = 0.0_wp
     1342       ENDIF
     1343    ENDDO
     1344
     1345!
     1346!-- Surface fluxes: answs = aerosol number, amsws = aerosol mass
     1347!
     1348!-- Horizontal surfaces: default type
     1349    DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
     1350       ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) )
     1351       ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) )
     1352       surf_def_h(l)%answs = 0.0_wp
     1353       surf_def_h(l)%amsws = 0.0_wp
     1354    ENDDO
     1355!
     1356!-- Horizontal surfaces: natural type
     1357    ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) )
     1358    ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) )
     1359    surf_lsm_h%answs = 0.0_wp
     1360    surf_lsm_h%amsws = 0.0_wp
     1361!
     1362!-- Horizontal surfaces: urban type
     1363    ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) )
     1364    ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) )
     1365    surf_usm_h%answs = 0.0_wp
     1366    surf_usm_h%amsws = 0.0_wp
     1367
     1368!
     1369!-- Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing
     1370    DO  l = 0, 3
     1371       ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) )
     1372       surf_def_v(l)%answs = 0.0_wp
     1373       ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
     1374       surf_def_v(l)%amsws = 0.0_wp
     1375
     1376       ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) )
     1377       surf_lsm_v(l)%answs = 0.0_wp
     1378       ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
     1379       surf_lsm_v(l)%amsws = 0.0_wp
     1380
     1381       ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) )
     1382       surf_usm_v(l)%answs = 0.0_wp
     1383       ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) )
     1384       surf_usm_v(l)%amsws = 0.0_wp
     1385
     1386    ENDDO
     1387
     1388!
     1389!-- Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV)
     1390!-- (number concentration (#/m3) )
     1391!
     1392!-- If chemistry is on, read gas phase concentrations from there. Otherwise,
     1393!-- allocate salsa_gas array.
     1394
     1395    IF ( air_chemistry )  THEN
     1396       DO  lsp = 1, nvar
     1397          SELECT CASE ( TRIM( chem_species(lsp)%name ) )
     1398             CASE ( 'H2SO4', 'h2so4' )
     1399                gases_available = gases_available + 1
     1400                gas_index_chem(1) = lsp
     1401             CASE ( 'HNO3', 'hno3' )
     1402                gases_available = gases_available + 1
     1403                gas_index_chem(2) = lsp
     1404             CASE ( 'NH3', 'nh3' )
     1405                gases_available = gases_available + 1
     1406                gas_index_chem(3) = lsp
     1407             CASE ( 'OCNV', 'ocnv' )
     1408                gases_available = gases_available + 1
     1409                gas_index_chem(4) = lsp
     1410             CASE ( 'OCSV', 'ocsv' )
     1411                gases_available = gases_available + 1
     1412                gas_index_chem(5) = lsp
     1413          END SELECT
     1414       ENDDO
     1415
     1416       IF ( gases_available == ngases_salsa )  THEN
     1417          salsa_gases_from_chem = .TRUE.
     1418       ELSE
     1419          WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// &
     1420                                     'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)'
     1421       CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 )
     1422       ENDIF
     1423
     1424    ELSE
     1425
     1426       ALLOCATE( salsa_gas(ngases_salsa) )
     1427       ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
     1428                 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa),                 &
     1429                 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) )
     1430       gconc_1 = 0.0_wp
     1431       gconc_2 = 0.0_wp
     1432       gconc_3 = 0.0_wp
     1433
     1434       DO i = 1, ngases_salsa
     1435          salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => gconc_1(:,:,:,i)
     1436          salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => gconc_2(:,:,:,i)
     1437          salsa_gas(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_3(:,:,:,i)
     1438          ALLOCATE( salsa_gas(i)%flux_s(nzb+1:nzt,0:threads_per_task-1),       &
     1439                    salsa_gas(i)%diss_s(nzb+1:nzt,0:threads_per_task-1),       &
     1440                    salsa_gas(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
     1441                    salsa_gas(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task-1),&
     1442                    salsa_gas(i)%init(nzb:nzt+1),                              &
     1443                    salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     1444          salsa_gas(i)%init = nclim
     1445          IF ( include_emission )  THEN
     1446             ALLOCATE( salsa_gas(i)%source(nys:nys,nxl:nxr) )
     1447             salsa_gas(i)%source = 0.0_wp
     1448          ENDIF
     1449       ENDDO
     1450!
     1451!--    Surface fluxes: gtsws = gaseous tracer flux
     1452!
     1453!--    Horizontal surfaces: default type
     1454       DO  l = 0, 2   ! upward (l=0), downward (l=1) and model top (l=2)
     1455          ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) )
     1456          surf_def_h(l)%gtsws = 0.0_wp
     1457       ENDDO
     1458!--    Horizontal surfaces: natural type
     1459       ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) )
     1460       surf_lsm_h%gtsws = 0.0_wp
     1461!--    Horizontal surfaces: urban type
     1462       ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) )
     1463       surf_usm_h%gtsws = 0.0_wp
     1464!
     1465!--    Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and
     1466!--    westward (l=3) facing
     1467       DO  l = 0, 3
     1468          ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) )
     1469          surf_def_v(l)%gtsws = 0.0_wp
     1470          ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) )
     1471          surf_lsm_v(l)%gtsws = 0.0_wp
     1472          ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) )
     1473          surf_usm_v(l)%gtsws = 0.0_wp
     1474       ENDDO
     1475    ENDIF
     1476
     1477    IF ( ws_scheme_sca )  THEN
     1478
     1479       IF ( salsa )  THEN
     1480          ALLOCATE( sums_salsa_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     1481          sums_salsa_ws_l = 0.0_wp
     1482       ENDIF
     1483
     1484    ENDIF
     1485!
     1486!-- Set control flags for decycling only at lateral boundary cores. Within the inner cores the
     1487!-- decycle flag is set to .FALSE.. Even though it does not affect the setting of chemistry boundary
     1488!-- conditions, this flag is used to set advection control flags appropriately.
     1489    decycle_salsa_lr = MERGE( decycle_salsa_lr, .FALSE., nxl == 0  .OR.  nxr == nx )
     1490    decycle_salsa_ns = MERGE( decycle_salsa_ns, .FALSE., nys == 0  .OR.  nyn == ny )
     1491!
     1492!-- Decycling can be applied separately for aerosol variables, while wind and other scalars may have
     1493!-- cyclic or nested boundary conditions. However, large gradients near the boundaries may produce
     1494!-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
     1495!-- applied near these boundaries. To get rid-off this, set-up additional flags that control the
     1496!-- order of the scalar advection scheme near the lateral boundaries for passive scalars with
     1497!-- decycling.
     1498    IF ( scalar_advec == 'ws-scheme' )  THEN
     1499       ALLOCATE( salsa_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1500!
     1501!--    In case of decycling, set Neuman boundary conditions for wall_flags_0 bit 31 instead of
     1502!--    cyclic boundary conditions. Bit 31 is used to identify extended degradation zones (please see
     1503!--    the following comment). Note, since several also other modules may access this bit but may
     1504!--    have other boundary conditions, the original value of wall_flags_0 bit 31 must not be
     1505!--    modified. Hence, store the boundary conditions directly on salsa_advc_flags_s.
     1506!--    salsa_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used
     1507!--    to control the numerical order.
     1508!--    Initialize with flag 31 only.
     1509       salsa_advc_flags_s = 0
     1510       salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_0, 31 ) )
     1511
     1512       IF ( decycle_salsa_ns )  THEN
     1513          IF ( nys == 0 )  THEN
     1514             DO  i = 1, nbgp
     1515                salsa_advc_flags_s(:,nys-i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nys,:), 31 ),   &
     1516                                                       IBCLR( salsa_advc_flags_s(:,nys,:), 31 ),   &
     1517                                                       BTEST( salsa_advc_flags_s(:,nys,:), 31 ) )
     1518             ENDDO
     1519          ENDIF
     1520          IF ( nyn == ny )  THEN
     1521             DO  i = 1, nbgp
     1522                salsa_advc_flags_s(:,nyn+i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nyn,:), 31 ),   &
     1523                                                       IBCLR( salsa_advc_flags_s(:,nyn,:), 31 ),   &
     1524                                                       BTEST( salsa_advc_flags_s(:,nyn,:), 31 ) )
     1525             ENDDO
     1526          ENDIF
     1527       ENDIF
     1528       IF ( decycle_salsa_lr )  THEN
     1529          IF ( nxl == 0 )  THEN
     1530             DO  i = 1, nbgp
     1531                salsa_advc_flags_s(:,:,nxl-i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxl), 31 ),   &
     1532                                                       IBCLR( salsa_advc_flags_s(:,:,nxl), 31 ),   &
     1533                                                       BTEST( salsa_advc_flags_s(:,:,nxl), 31 ) )
     1534             ENDDO
     1535          ENDIF
     1536          IF ( nxr == nx )  THEN
     1537             DO  i = 1, nbgp
     1538                salsa_advc_flags_s(:,:,nxr+i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxr), 31 ),   &
     1539                                                       IBCLR( salsa_advc_flags_s(:,:,nxr), 31 ),   &
     1540                                                       BTEST( salsa_advc_flags_s(:,:,nxr), 31 ) )
     1541             ENDDO
     1542          ENDIF
     1543       ENDIF
     1544!
     1545!--    To initialise the advection flags appropriately, pass the boundary flags to
     1546!--    ws_init_flags_scalar. The last argument in ws_init_flags_scalar indicates that a passive
     1547!--    scalar is being treated and the horizontal advection terms are degraded already 2 grid points
     1548!--    before the lateral boundary. Also, extended degradation zones are applied, where
     1549!--    horizontal advection of scalars is discretised by the first-order scheme at all grid points
     1550!--    in the vicinity of buildings (<= 3 grid points). Even though no building is within the
     1551!--    numerical stencil, the first-order scheme is used. At fourth and fifth grid points, the order
     1552!--    of the horizontal advection scheme is successively upgraded.
     1553!--    These degradations of the advection scheme are done to avoid stationary numerical
     1554!--    oscillations, which are responsible for high concentration maxima that may appear e.g. under
     1555!--    shear-free stable conditions.
     1556       CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_salsa_lr,    &
     1557                                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_salsa_ns,    &
     1558                                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_salsa_lr,    &
     1559                                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_salsa_ns,    &
     1560                                  salsa_advc_flags_s, .TRUE. )
     1561    ENDIF
     1562
     1563
     1564 END SUBROUTINE salsa_init_arrays
     1565
     1566!------------------------------------------------------------------------------!
     1567! Description:
     1568! ------------
     1569!> Initialization of SALSA. Based on salsa_initialize in UCLALES-SALSA.
     1570!> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALES-SALSA are
     1571!> also merged here.
     1572!------------------------------------------------------------------------------!
     1573 SUBROUTINE salsa_init
     1574
     1575    IMPLICIT NONE
     1576
     1577    INTEGER(iwp) :: i   !<
     1578    INTEGER(iwp) :: ib  !< loop index for aerosol number bins
     1579    INTEGER(iwp) :: ic  !< loop index for aerosol mass bins
     1580    INTEGER(iwp) :: ig  !< loop index for gases
     1581    INTEGER(iwp) :: j   !<
     1582
     1583    IF ( debug_output )  CALL debug_message( 'salsa_init', 'start' )
     1584
     1585    bin_low_limits = 0.0_wp
     1586    k_topo_top     = 0
     1587    nsect          = 0.0_wp
     1588    massacc        = 1.0_wp
    15421589!
    15431590!-- Initialise
    1544 !
    1545 !-- Aerosol size distribution (TYPE t_section)
    1546     aero(:)%dwet     = 1.0E-10_wp
    1547     aero(:)%veqh2o   = 1.0E-10_wp
    1548     aero(:)%numc     = nclim
    1549     aero(:)%core     = 1.0E-10_wp
    1550     DO ic = 1, maxspec+1    ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
    1551        aero(:)%volc(ic) = 0.0_wp
    1552     ENDDO
    1553 
    15541591    IF ( nldepo )  sedim_vd = 0.0_wp
    15551592
     
    15831620!-- Aerosol radius in each bin: dry and wet (m)
    15841621    ra_dry = 1.0E-10_wp
    1585 !
    1586 !-- Initialise aerosol tracers
    1587     aero(:)%vhilim   = 0.0_wp
    1588     aero(:)%vlolim   = 0.0_wp
    1589     aero(:)%vratiohi = 0.0_wp
    1590     aero(:)%vratiolo = 0.0_wp
    1591     aero(:)%dmid     = 0.0_wp
    1592 !
    1593 !-- Initialise the sectional particle size distribution
    1594     CALL set_sizebins
    15951622!
    15961623!-- Initialise location-dependent aerosol size distributions and chemical compositions:
     
    16391666       ENDIF
    16401667    ENDIF
     1668!
     1669!-- Partition and dissolutional growth by gaseous HNO3 and NH3
     1670    IF ( index_no > 0  .AND.  index_nh > 0  .AND.  index_so4 > 0 )  lspartition = .TRUE.
    16411671
    16421672    IF ( debug_output )  CALL debug_message( 'salsa_init', 'end' )
     
    16761706
    16771707    REAL(wp) ::  ratio_d  !< ratio of the upper and lower diameter of subranges
     1708
     1709    aero(:)%dwet     = 1.0E-10_wp
     1710    aero(:)%veqh2o   = 1.0E-10_wp
     1711    aero(:)%numc     = nclim
     1712    aero(:)%core     = 1.0E-10_wp
     1713    DO  cc = 1, maxspec+1    ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O
     1714       aero(:)%volc(cc) = 0.0_wp
     1715    ENDDO
    16781716!
    16791717!-- vlolim&vhilim: min & max *dry* volumes [fxm]
     
    17371775
    17381776    USE netcdf_data_input_mod,                                                                     &
    1739         ONLY:  close_input_file, get_attribute, get_variable,                                      &
     1777        ONLY:  check_existence, close_input_file, get_attribute, get_variable,                     &
     1778               inquire_num_variables, inquire_variable_names,                                      &
    17401779               netcdf_data_input_get_dimension_length, open_read_file
    17411780
    17421781    IMPLICIT NONE
    17431782
    1744     CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cc_name  !< chemical component name
     1783    CHARACTER(LEN=25),  DIMENSION(:), ALLOCATABLE ::  cc_name    !< chemical component name
     1784    CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names
    17451785
    17461786    INTEGER(iwp) ::  ee        !< index: end
     
    17531793    INTEGER(iwp) ::  k         !< loop index: z-direction
    17541794    INTEGER(iwp) ::  lod_aero  !< level of detail of inital aerosol concentrations
    1755     INTEGER(iwp) ::  pr_nbins  !< Number of aerosol size bins in file
    1756     INTEGER(iwp) ::  pr_ncc    !< Number of aerosol chemical components in file
    1757     INTEGER(iwp) ::  pr_nz     !< Number of vertical grid-points in file
     1795    INTEGER(iwp) ::  num_vars  !< number of variables
     1796    INTEGER(iwp) ::  pr_nbins  !< number of aerosol size bins in file
     1797    INTEGER(iwp) ::  pr_ncc    !< number of aerosol chemical components in file
     1798    INTEGER(iwp) ::  pr_nz     !< number of vertical grid-points in file
    17581799    INTEGER(iwp) ::  prunmode  !< running mode of SALSA
    17591800    INTEGER(iwp) ::  ss        !< index: start
    17601801
    1761     INTEGER(iwp), DIMENSION(maxspec) ::  cc_input_to_model
     1802    INTEGER(iwp), DIMENSION(maxspec) ::  cc_in2mod
    17621803
    17631804    LOGICAL  ::  netcdf_extend = .FALSE. !< Flag: netcdf file exists
     
    17701811    REAL(wp), DIMENSION(0:nz+1) ::  pmfoc1a !< mass fraction of OC in 1a
    17711812
    1772     REAL(wp), DIMENSION(0:nz+1,nbins_aerosol)   ::  pndist  !< size dist as a function of height (#/m3)
     1813    REAL(wp), DIMENSION(0:nz+1,nbins_aerosol)   ::  pndist  !< vertical profile of size dist. (#/m3)
    17731814    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2a   !< mass distributions in subrange 2a
    17741815    REAL(wp), DIMENSION(0:nz+1,maxspec)         ::  pmf2b   !< mass distributions in subrange 2b
     
    17801821    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pr_mass_fracs_b  !< and b
    17811822
    1782     cc_input_to_model = 0
     1823    cc_in2mod = 0
    17831824    prunmode = 1
    17841825!
    17851826!-- Bin mean aerosol particle volume (m3)
    1786     core(:) = 0.0_wp
    17871827    core(1:nbins_aerosol) = api6 * aero(1:nbins_aerosol)%dmid**3
    17881828!
     
    18061846          CALL open_read_file( TRIM( input_file_dynamic ) // TRIM( coupling_char ), id_dyn )
    18071847!
    1808 !--       Inquire dimensions:
     1848!--       At first, inquire all variable names
     1849          CALL inquire_num_variables( id_dyn, num_vars )
     1850!
     1851!--       Allocate memory to store variable names
     1852          ALLOCATE( var_names(1:num_vars) )
     1853          CALL inquire_variable_names( id_dyn, var_names )
     1854!
     1855!--       Inquire vertical dimension and number of aerosol chemical components
    18091856          CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' )
    18101857          IF ( pr_nz /= nz )  THEN
     
    18241871          CALL get_variable( id_dyn, 'z', pr_z )
    18251872!
    1826 !--       Read name and index of chemical components
    1827           CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc )
     1873!--       Read the names of chemical components
     1874          IF ( check_existence( var_names, 'composition_name' ) )  THEN
     1875             CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc )
     1876          ELSE
     1877             WRITE( message_string, * ) 'Missing composition_name in ' // TRIM( input_file_dynamic )
     1878             CALL message( 'aerosol_init', 'PA0655', 1, 2, 0, 6, 0 )
     1879          ENDIF
     1880!
     1881!--       Define the index of each chemical component in the model
    18281882          DO  ic = 1, pr_ncc
    18291883             SELECT CASE ( TRIM( cc_name(ic) ) )
    18301884                CASE ( 'H2SO4', 'SO4', 'h2so4', 'so4' )
    1831                    cc_input_to_model(1) = ic
     1885                   cc_in2mod(1) = ic
    18321886                CASE ( 'OC', 'oc' )
    1833                    cc_input_to_model(2) = ic
     1887                   cc_in2mod(2) = ic
    18341888                CASE ( 'BC', 'bc' )
    1835                    cc_input_to_model(3) = ic
     1889                   cc_in2mod(3) = ic
    18361890                CASE ( 'DU', 'du' )
    1837                    cc_input_to_model(4) = ic
     1891                   cc_in2mod(4) = ic
    18381892                CASE ( 'SS', 'ss' )
    1839                    cc_input_to_model(5) = ic
     1893                   cc_in2mod(5) = ic
    18401894                CASE ( 'HNO3', 'hno3', 'NO', 'no' )
    1841                    cc_input_to_model(6) = ic
     1895                   cc_in2mod(6) = ic
    18421896                CASE ( 'NH3', 'nh3', 'NH', 'nh' )
    1843                    cc_input_to_model(7) = ic
     1897                   cc_in2mod(7) = ic
    18441898             END SELECT
    18451899          ENDDO
    18461900
    1847           IF ( SUM( cc_input_to_model ) == 0 )  THEN
     1901          IF ( SUM( cc_in2mod ) == 0 )  THEN
    18481902             message_string = 'None of the aerosol chemical components in ' // TRIM(               &
    18491903                              input_file_dynamic ) // ' correspond to ones applied in SALSA.'
     
    18521906!
    18531907!--       Vertical profiles of mass fractions of different chemical components:
    1854           CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a,              &
    1855                              0, pr_ncc-1, 0, pr_nz-1 )
     1908          IF ( check_existence( var_names, 'init_atmosphere_mass_fracs_a' ) )  THEN
     1909             CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a,           &
     1910                                0, pr_ncc-1, 0, pr_nz-1 )
     1911          ELSE
     1912             WRITE( message_string, * ) 'Missing init_atmosphere_mass_fracs_a in ' //              &
     1913                                        TRIM( input_file_dynamic )
     1914             CALL message( 'aerosol_init', 'PA0656', 1, 2, 0, 6, 0 )
     1915          ENDIF
    18561916          CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_b', pr_mass_fracs_b,              &
    18571917                             0, pr_ncc-1, 0, pr_nz-1  )
     
    18591919!--       Match the input data with the chemical composition applied in the model
    18601920          DO  ic = 1, maxspec
    1861              ss = cc_input_to_model(ic)
     1921             ss = cc_in2mod(ic)
    18621922             IF ( ss == 0 )  CYCLE
    18631923             pmf2a(nzb+1:nzt+1,ic) = pr_mass_fracs_a(nzb:nzt,ss)
     
    18651925          ENDDO
    18661926!
    1867 !--       Aerosol concentrations: lod=1 (total PM) or lod=2 (sectional number size distribution)
     1927!--       Aerosol concentrations: lod=1 (vertical profile of sectional number size distribution)
    18681928          CALL get_attribute( id_dyn, 'lod', lod_aero, .FALSE., 'init_atmosphere_aerosol' )
    18691929          IF ( lod_aero /= 1 )  THEN
     
    30203080!--    Condensation
    30213081       IF ( lscnd )   THEN
    3022           CALL condensation( lo_aero, zgso4, zgocnv, zgocsv,  zghno3, zgnh3, in_cw(k), in_cs(k),      &
     3082          CALL condensation( lo_aero, zgso4, zgocnv, zgocsv,  zghno3, zgnh3, in_cw(k), in_cs(k),   &
    30233083                             in_t(k), in_p(k), dt_salsa, prtcl )
    30243084       ENDIF
     
    30263086!--    Deposition
    30273087       IF ( lsdepo )  THEN
    3028           CALL deposition( lo_aero, in_t(k), in_adn(k), in_u(k), in_lad, kvis(k), schmidt_num(k,:),   &
     3088          CALL deposition( lo_aero, in_t(k), in_adn(k), in_u(k), in_lad, kvis(k), schmidt_num(k,:),&
    30293089                           vd(k,:) )
    30303090       ENDIF
     
    32583318! Description:
    32593319! ------------
    3260 !> Set logical switches according to the host model state and user-specified
    3261 !> NAMELIST options.
     3320!> Set logical switches according to the salsa_parameters options.
    32623321!> Juha Tonttila, FMI, 2014
    32633322!> Only aerosol processes included, Mona Kurppa, UHel, 2017
     
    38323891!
    38333892!--    Brownian diffusion
    3834        v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                             &
     3893       v_bd = mag_u * c_brownian_diff * schmidt_num**( -0.66666666_wp ) *                          &
    38353894              ( mag_u * par_l / kvis_a )**( -0.5_wp )
    38363895!
    38373896!--    Interception
    3838        v_in = mag_u * c_interception * diameter / par_l * ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) )
     3897       v_in = mag_u * c_interception * diameter / par_l *                                          &
     3898              ( 2.0_wp + LOG( 2.0_wp * par_l / diameter ) )
    38393899!
    38403900!--    Impaction: Petroff (2009) Eq. 18
     
    41794239    REAL(wp), DIMENSION(maxspec+1)     ::  zplusterm  !< coagulation gain in a bin (for each
    41804240                                                      !< chemical compound)
    4181     REAL(wp), DIMENSION(nbins_aerosol,nbins_aerosol) ::  zcc  !< updated coagulation coefficients (m3/s)
     4241    REAL(wp), DIMENSION(nbins_aerosol,nbins_aerosol) ::  zcc  !< updated coagulation coeff. (m3/s)
    41824242
    41834243    TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) ::  paero  !< Aerosol properties
     
    41864246    zdpart_nn = 0.0_wp
    41874247!
    4188 !-- 1) Coagulation to coarse mode calculated in a simplified way:
    4189 !--    CoagSink ~ Dp in continuum subrange, thus we calculate 'effective'
    4190 !--    number concentration of coarse particles
     4248!-- 1) Coagulation to coarse mode calculated in a simplified way:
     4249!--    CoagSink ~ Dp in continuum subrange --> 'effective' number conc. of coarse particles
    41914250
    41924251!-- 2) Updating coagulation coefficients
     
    44184477    fmdist = SQRT( tva(1)**2 + tva(2)**2 )
    44194478!
    4420 !-- 5) Coagulation coefficient = coalescence efficiency * collision kernel (m3/s) (eq. 15.33). 
     4479!-- 5) Coagulation coefficient = coalescence efficiency * collision kernel (m3/s) (eq. 15.33).
    44214480!--    Here assumed coalescence efficiency 1!!
    44224481    coagc = flux(1) / ( mdiam / ( mdiam + fmdist) + flux(1) / flux(2) )
     
    46324691!
    46334692!--    5.1.2) Semivolatile organic compound: all bins except subrange 1
    4634        zcs_ocsv = SUM( zcolrate(start_subrange_2a:end_subrange_2b) ) !< sink for semi-volatile organics
     4693       zcs_ocsv = SUM( zcolrate(start_subrange_2a:end_subrange_2b) ) !< sink for semi-volatile org.
    46354694       IF ( pcocsv > 1.0E+10_wp  .AND.  zcs_ocsv > 1.0E-30  .AND. is_used( prtcl,'OC') )  THEN
    46364695!
     
    48494908    zkocnv = 0.0_wp
    48504909
     4910    zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
     4911    zc_org   = pc_ocnv * 1.0E-6_wp   ! conc. of non-volatile OC to #/cm3
     4912    zmixnh3  = pc_nh3 * ptemp * argas / ( ppres * avo )
     4913
    48514914    SELECT CASE ( nsnucl )
    48524915!
     
    48544917       CASE(1)
    48554918
    4856           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    48574919          CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit,  zksa, zkocnv )
    48584920!
    4859 !--    Activation type nucleation
     4921!--    Activation type nucleation (See Riipinen et al. (2007), Atmos. Chem. Phys., 7(8), 1899-1914)
    48604922       CASE(2)
    4861 
    4862           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    4863           CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa,  znoc, zdcrit, zksa, zkocnv )
    4864           CALL actnucl( pc_sa, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv, act_coeff )
     4923!
     4924!--       Nucleation rate (#/(m3 s))
     4925          zc_h2so4  = MAX( zc_h2so4, 1.0E4_wp  )
     4926          zc_h2so4  = MIN( zc_h2so4, 1.0E11_wp )
     4927          zjnuc = act_coeff * pc_sa  ! (#/(m3 s))
     4928!
     4929!--       Organic compounds not involved when kinetic nucleation is assumed.
     4930          zdcrit  = 7.9375E-10_wp   ! (m)
     4931          zkocnv  = 0.0_wp
     4932          zksa    = 1.0_wp
     4933          znoc    = 0.0_wp
     4934          znsa    = 2.0_wp
    48654935!
    48664936!--    Kinetically limited nucleation of (NH4)HSO4 clusters
     4937!--    (See Sihto et al. (2006), Atmos. Chem. Phys., 6(12), 4079-4091.)
    48674938       CASE(3)
    4868 
    4869           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    4870           CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv )
    4871           CALL kinnucl( zc_h2so4, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv )
     4939!
     4940!--       Nucleation rate = coagcoeff*zpcsa**2 (#/(m3 s))
     4941          zc_h2so4  = MAX( zc_h2so4, 1.0E4_wp  )
     4942          zc_h2so4  = MIN( zc_h2so4, 1.0E11_wp )
     4943          zjnuc = 5.0E-13_wp * zc_h2so4**2.0_wp * 1.0E+6_wp
     4944!
     4945!--       Organic compounds not involved when kinetic nucleation is assumed.
     4946          zdcrit  = 7.9375E-10_wp   ! (m)
     4947          zkocnv  = 0.0_wp
     4948          zksa    = 1.0_wp
     4949          znoc    = 0.0_wp
     4950          znsa    = 2.0_wp
    48724951!
    48734952!--    Ternary H2SO4-H2O-NH3 nucleation
    48744953       CASE(4)
    48754954
    4876           zmixnh3 = pc_nh3 * ptemp * argas / ( ppres * avo )
    4877           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    48784955          CALL ternucl( zc_h2so4, zmixnh3, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv )
    48794956!
    48804957!--    Organic nucleation, J~[ORG] or J~[ORG]**2
     4958!--    (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.)
    48814959       CASE(5)
    4882 
    4883           zc_org = pc_ocnv * 1.0E-6_wp   ! conc. of non-volatile OC to #/cm3
    4884           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    4885           CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv )
    4886           CALL orgnucl( pc_ocnv, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv )
     4960!
     4961!--       Homomolecular nuleation rate
     4962          zjnuc = 1.3E-7_wp * pc_ocnv   ! (1/s) (Paasonen et al. Table 4: median a_org)
     4963!
     4964!--       H2SO4 not involved when pure organic nucleation is assumed.
     4965          zdcrit  = 1.5E-9  ! (m)
     4966          zkocnv  = 1.0_wp
     4967          zksa    = 0.0_wp
     4968          znoc    = 1.0_wp
     4969          znsa    = 0.0_wp
    48874970!
    48884971!--    Sum of H2SO4 and organic activation type nucleation, J~[H2SO4]+[ORG]
     4972!--    (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242)
    48894973       CASE(6)
    4890 
    4891           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    4892           CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv )
    4893           CALL sumnucl( pc_sa, pc_ocnv, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv )
    4894 !
    4895 !--    Heteromolecular nucleation, J~[H2SO4]*[ORG]
     4974!
     4975!--       Nucleation rate  (#/m3/s)
     4976          zjnuc = 6.1E-7_wp * pc_sa + 0.39E-7_wp * pc_ocnv   ! (Paasonen et al. Table 3.)
     4977!
     4978!--       Both organic compounds and H2SO4 are involved when sumnucleation is assumed.
     4979          zdcrit  = 1.5E-9_wp   ! (m)
     4980          zkocnv  = 1.0_wp
     4981          zksa    = 1.0_wp
     4982          znoc    = 1.0_wp
     4983          znsa    = 1.0_wp
     4984!
     4985!--    Heteromolecular nucleation, J~[H2SO4]*[ORG]
     4986!--    (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.)
    48964987       CASE(7)
    4897 
    4898           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    4899           zc_org = pc_ocnv * 1.0E-6_wp   ! conc. of non-volatile OC to #/cm3
    4900           CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv )
    4901           CALL hetnucl( zc_h2so4, zc_org, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv )
     4988!
     4989!--       Nucleation rate (#/m3/s)
     4990          zjnuc = 4.1E-14_wp * pc_sa * pc_ocnv * 1.0E6_wp   ! (Paasonen et al. Table 4: median)
     4991!
     4992!--       Both organic compounds and H2SO4 are involved when heteromolecular nucleation is assumed
     4993          zdcrit  = 1.5E-9_wp   ! (m)
     4994          zkocnv  = 1.0_wp
     4995          zksa    = 1.0_wp
     4996          znoc    = 1.0_wp
     4997          znsa    = 1.0_wp
    49024998!
    49034999!--    Homomolecular nucleation of H2SO4 and heteromolecular nucleation of H2SO4 and organic vapour,
    49045000!--    J~[H2SO4]**2 + [H2SO4]*[ORG] (EUCAARI project)
     5001!--    (See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242)
    49055002       CASE(8)
    4906           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    4907           zc_org = pc_ocnv * 1.0E-6_wp   ! conc. of non-volatile OC to #/cm3
    4908           CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv )
    4909           CALL SAnucl( zc_h2so4, zc_org, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv )
     5003!
     5004!--       Nucleation rate (#/m3/s)
     5005          zjnuc = ( 1.1E-14_wp * zc_h2so4**2 + 3.2E-14_wp * zc_h2so4 * zc_org ) * 1.0E+6_wp
     5006!
     5007!--       Both organic compounds and H2SO4 are involved when SAnucleation is assumed
     5008          zdcrit  = 1.5E-9_wp   ! (m)
     5009          zkocnv  = 1.0_wp
     5010          zksa    = 1.0_wp
     5011          znoc    = 1.0_wp
     5012          znsa    = 3.0_wp
    49105013!
    49115014!--    Homomolecular nucleation of H2SO4 and organic vapour and heteromolecular nucleation of H2SO4
    49125015!--    and organic vapour, J~[H2SO4]**2 + [H2SO4]*[ORG]+[ORG]**2 (EUCAARI project)
    49135016       CASE(9)
    4914 
    4915           zc_h2so4 = pc_sa * 1.0E-6_wp   ! sulphuric acid conc. to #/cm3
    4916           zc_org = pc_ocnv * 1.0E-6_wp   ! conc. of non-volatile OC to #/cm3
    4917           CALL binnucl( zc_h2so4, ptemp, prh, zjnuc, znsa, znoc, zdcrit, zksa, zkocnv )
    4918           CALL SAORGnucl( zc_h2so4, zc_org, zjnuc, zdcrit, znsa, znoc, zksa, zkocnv )
     5017!
     5018!--       Nucleation rate (#/m3/s)
     5019          zjnuc = ( 1.4E-14_wp * zc_h2so4**2 + 2.6E-14_wp * zc_h2so4 * zc_org + 0.037E-14_wp *     &
     5020                    zc_org**2 ) * 1.0E+6_wp
     5021!
     5022!--       Both organic compounds and H2SO4 are involved when SAORGnucleation is assumed
     5023          zdcrit  = 1.5E-9_wp   ! (m)
     5024          zkocnv  = 1.0_wp
     5025          zksa    = 1.0_wp
     5026          znoc    = 3.0_wp
     5027          znsa    = 3.0_wp
    49195028
    49205029    END SELECT
     
    50105119!
    50115120!--    Defining the parameter m (zm_para) for calculating the coagulation sink onto background
    5012 !--    particles (Eq. 5&6 in Lehtinen et al. 2007). The growth is investigated between 
     5121!--    particles (Eq. 5&6 in Lehtinen et al. 2007). The growth is investigated between
    50135122!--    [d1,reglim(1)] = [zdcrit,3nm] and m = LOG( CoagS_dx / CoagX_zdcrit ) / LOG( reglim / zdcrit )
    50145123!--    (Lehtinen et al. 2007, Eq. 6).
    5015 !--    The steps for the coagulation sink for reglim = 3nm and zdcrit ~= 1nm are explained in 
    5016 !--    Kulmala et al. (2001). The particles of diameter zdcrit ~1.14 nm  and reglim = 3nm are both 
     5124!--    The steps for the coagulation sink for reglim = 3nm and zdcrit ~= 1nm are explained in
     5125!--    Kulmala et al. (2001). The particles of diameter zdcrit ~1.14 nm  and reglim = 3nm are both
    50175126!--    in turn the "number 1" variables (Kulmala et al. 2001).
    50185127!--    c = critical (1nm), x = 3nm, 2 = wet or mean droplet
     
    51295238             z_gr_tot = z_gr_clust * 2.77777777E-7_wp +  1.5708E-6_wp * zlambda * zdcrit**3 *      &
    51305239                      ( z_n_nuc * 1.0E-6_wp ) * zcv_c * avo * 2.77777777E-7_wp ! (Eq. 3)
    5131              zeta = - zcoagstot / ( ( zm_para + 1.0_wp ) * z_gr_tot * ( zdcrit**zm_para ) ) ! (Eq. 7b)
     5240             zeta = - zcoagstot / ( ( zm_para + 1.0_wp ) * z_gr_tot * ( zdcrit**zm_para ) ) ! (Eq.7b)
    51325241!
    51335242!--          Calculate Eq. 7a (Taylor series for the number of particles between [d1,dx])
     
    51395248          z_gr_tot = z_gr_clust * 1.0E-9_wp / 3600.0_wp + 1.5708E-6_wp *  zlambda * zdcrit**3 *    &
    51405249                   ( z_n_nuc * 1.0E-6_wp ) * zcv_c * avo * 1.0E-9_wp / 3600.0_wp !< (m/s)
    5141           zj3 = zjnuc * EXP( MIN( 0.0_wp, -zgamma * zdcrit * zcoagstot / z_gr_tot ) ) ! (#/m3s, Eq. 5a)
     5250          zj3 = zjnuc * EXP( MIN( 0.0_wp, -zgamma * zdcrit * zcoagstot / z_gr_tot ) ) ! (#/m3s, Eq.5a)
    51425251
    51435252       ENDIF
     
    51455254!
    51465255!-- If J3 very small (< 1 #/cm3), neglect particle formation. In real atmosphere this would mean
    5147 !-- that clusters form but coagulate to pre-existing particles who gain sulphate. Since 
    5148 !-- CoagS ~ CS (4piD*CS'), we do *not* update H2SO4 concentration here but let condensation take 
     5256!-- that clusters form but coagulate to pre-existing particles who gain sulphate. Since
     5257!-- CoagS ~ CS (4piD*CS'), we do *not* update H2SO4 concentration here but let condensation take
    51495258!-- care of it. Formation mass rate of molecules (molec/m3s) for 1: H2SO4 and 2: organic vapour
    51505259    pj3n3(1) = zj3 * n3 * pxsa
     
    51565265! Description:
    51575266! ------------
    5158 !> Calculate the nucleation rate and the size of critical clusters assuming 
     5267!> Calculate the nucleation rate and the size of critical clusters assuming
    51595268!> binary nucleation.
    5160 !> Parametrisation according to Vehkamaki et al. (2002), J. Geophys. Res., 
    5161 !> 107(D22), 4622. Called from subroutine nucleation. 
     5269!> Parametrisation according to Vehkamaki et al. (2002), J. Geophys. Res.,
     5270!> 107(D22), 4622. Called from subroutine nucleation.
    51625271!------------------------------------------------------------------------------!
    51635272 SUBROUTINE binnucl( pc_sa, ptemp, prh, pnuc_rate, pn_crit_sa, pn_crit_ocnv, pd_crit, pk_sa,       &
     
    55195628
    55205629 END SUBROUTINE ternucl
    5521 
    5522 !------------------------------------------------------------------------------!
    5523 ! Description:
    5524 ! ------------
    5525 !> Calculate the nucleation rate and the size of critical clusters assuming
    5526 !> kinetic nucleation. Each sulphuric acid molecule forms an (NH4)HSO4 molecule
    5527 !> in the atmosphere and two colliding (NH4)HSO4 molecules form a stable
    5528 !> cluster. See Sihto et al. (2006), Atmos. Chem. Phys., 6(12), 4079-4091.
    5529 !>
    5530 !> Below the following assumption have been made:
    5531 !>  nucrate = coagcoeff*zpcsa**2
    5532 !>  coagcoeff = 8*sqrt(3*boltz*ptemp*r_abs/dens_abs)
    5533 !>  r_abs = 0.315d-9 radius of bisulphate molecule [m]
    5534 !>  dens_abs = 1465  density of - " - [kg/m3]
    5535 !------------------------------------------------------------------------------!
    5536  SUBROUTINE kinnucl( pc_sa, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )
    5537 
    5538     IMPLICIT NONE
    5539 
    5540     REAL(wp), INTENT(in) ::  pc_sa  !< H2SO4 conc. (#/m3)
    5541 
    5542     REAL(wp), INTENT(out) ::  pd_crit  !< critical diameter of clusters (m)
    5543     REAL(wp), INTENT(out) ::  pk_ocnv  !< if pk_ocnv = 1, organic compounds participate in nucleation
    5544     REAL(wp), INTENT(out) ::  pk_sa    !< if pk_sa = 1, H2SO4 is participate in nucleation
    5545     REAL(wp), INTENT(out) ::  pn_crit_ocnv  !< number of organic molecules in cluster (#)
    5546     REAL(wp), INTENT(out) ::  pn_crit_sa    !< number of H2SO4 molecules in cluster (#)
    5547     REAL(wp), INTENT(out) ::  pnuc_rate     !< nucl. rate (#/(m3 s))
    5548 !
    5549 !-- Nucleation rate (#/(m3 s))
    5550     pnuc_rate = 5.0E-13_wp * pc_sa**2.0_wp * 1.0E+6_wp
    5551 !
    5552 !-- Organic compounds not involved when kinetic nucleation is assumed.
    5553     pn_crit_sa   = 2.0_wp
    5554     pn_crit_ocnv = 0.0_wp
    5555     pk_sa        = 1.0_wp
    5556     pk_ocnv      = 0.0_wp
    5557     pd_crit      = 7.9375E-10_wp   ! (m)
    5558 
    5559  END SUBROUTINE kinnucl
    5560 
    5561 !------------------------------------------------------------------------------!
    5562 ! Description:
    5563 ! ------------
    5564 !> Calculate the nucleation rate and the size of critical clusters assuming
    5565 !> activation type nucleation.
    5566 !> See Riipinen et al. (2007), Atmos. Chem. Phys., 7(8), 1899-1914.
    5567 !------------------------------------------------------------------------------!
    5568  SUBROUTINE actnucl( psa_conc, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv, activ )
    5569 
    5570     IMPLICIT NONE
    5571 
    5572     REAL(wp), INTENT(in) ::  activ     !< activation coefficient (1e-7 by default)
    5573     REAL(wp), INTENT(in) ::  psa_conc  !< H2SO4 conc. (#/m3)
    5574 
    5575     REAL(wp), INTENT(out) ::  pd_crit  !< critical diameter of clusters (m)
    5576     REAL(wp), INTENT(out) ::  pk_ocnv  !< if pk_ocnv = 1, organic compounds participate in nucleation
    5577     REAL(wp), INTENT(out) ::  pk_sa    !< if pk_sa = 1, H2SO4 participate in nucleation
    5578     REAL(wp), INTENT(out) ::  pn_crit_ocnv  !< number of organic molecules in cluster (#)
    5579     REAL(wp), INTENT(out) ::  pn_crit_sa    !< number of H2SO4 molecules in cluster (#)
    5580     REAL(wp), INTENT(out) ::  pnuc_rate     !< nucl. rate (#/(m3 s))
    5581 !
    5582 !-- Nucleation rate (#/(m3 s))
    5583     pnuc_rate = activ * psa_conc   ! (#/(m3 s))
    5584 !
    5585 !-- Organic compounds not involved when kinetic nucleation is assumed.
    5586     pn_crit_sa   = 2.0_wp
    5587     pn_crit_ocnv = 0.0_wp
    5588     pk_sa        = 1.0_wp
    5589     pk_ocnv      = 0.0_wp
    5590     pd_crit      = 7.9375E-10_wp   ! (m)
    5591 
    5592  END SUBROUTINE actnucl
    5593 
    5594 !------------------------------------------------------------------------------!
    5595 ! Description:
    5596 ! ------------
    5597 !> Conciders only the organic matter in nucleation. Paasonen et al. (2010)
    5598 !> determined particle formation rates for 2 nm particles, J2, from different
    5599 !> kind of combinations of sulphuric acid and organic matter concentration.
    5600 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.
    5601 !------------------------------------------------------------------------------!
    5602  SUBROUTINE orgnucl( pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv,     &
    5603                      pk_sa, pk_ocnv )
    5604 
    5605     IMPLICIT NONE
    5606 
    5607     REAL(wp) ::  a_org = 1.3E-7_wp  !< (1/s) (Paasonen et al. Table 4: median)
    5608 
    5609     REAL(wp), INTENT(in) ::  pc_org   !< organic vapour concentration (#/m3)
    5610 
    5611     REAL(wp), INTENT(out) ::  pd_crit  !< critical diameter of clusters (m)
    5612     REAL(wp), INTENT(out) ::  pk_ocnv  !< if pk_ocnv = 1, organic compounds participate in nucleation
    5613     REAL(wp), INTENT(out) ::  pk_sa    !< if pk_sa = 1, H2SO4 participate in nucleation
    5614     REAL(wp), INTENT(out) ::  pn_crit_ocnv  !< number of organic molecules in cluster (#)
    5615     REAL(wp), INTENT(out) ::  pn_crit_sa    !< number of H2SO4 molecules in cluster (#)
    5616     REAL(wp), INTENT(out) ::  pnuc_rate     !< nucl. rate (#/(m3 s))
    5617 !
    5618 !-- Homomolecular nuleation rate
    5619     pnuc_rate = a_org * pc_org
    5620 !
    5621 !-- H2SO4 not involved when pure organic nucleation is assumed.
    5622     pn_crit_sa   = 0.0_wp
    5623     pn_crit_ocnv = 1.0_wp
    5624     pk_sa        = 0.0_wp
    5625     pk_ocnv      = 1.0_wp
    5626     pd_crit      = 1.5E-9_wp   ! (m)
    5627 
    5628  END SUBROUTINE orgnucl
    5629 
    5630 !------------------------------------------------------------------------------!
    5631 ! Description:
    5632 ! ------------
    5633 !> Conciders both the organic vapor and H2SO4 in nucleation - activation type
    5634 !> of nucleation.
    5635 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.
    5636 !------------------------------------------------------------------------------!
    5637  SUBROUTINE sumnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )
    5638 
    5639     IMPLICIT NONE
    5640 
    5641     REAL(wp) ::  a_s1 = 6.1E-7_wp   !< (1/s)
    5642     REAL(wp) ::  a_s2 = 0.39E-7_wp  !< (1/s) (Paasonen et al. Table 3.)
    5643 
    5644     REAL(wp), INTENT(in) ::  pc_org   !< organic vapour concentration (#/m3)
    5645     REAL(wp), INTENT(in) ::  pc_sa    !< H2SO4 conc. (#/m3)
    5646 
    5647     REAL(wp), INTENT(out) ::  pd_crit  !< critical diameter of clusters (m)
    5648     REAL(wp), INTENT(out) ::  pk_ocnv  !< if pk_ocnv = 1, organic compounds participate in nucleation
    5649     REAL(wp), INTENT(out) ::  pk_sa    !< if pk_sa = 1, H2SO4 participate in nucleation
    5650     REAL(wp), INTENT(out) ::  pn_crit_ocnv  !< number of organic molecules in cluster (#)
    5651     REAL(wp), INTENT(out) ::  pn_crit_sa    !< number of H2SO4 molecules in cluster (#)
    5652     REAL(wp), INTENT(out) ::  pnuc_rate     !< nucl. rate (#/(m3 s))
    5653 !
    5654 !-- Nucleation rate  (#/m3/s)
    5655     pnuc_rate = a_s1 * pc_sa + a_s2 * pc_org
    5656 !
    5657 !-- Both organic compounds and H2SO4 are involved when sumnucleation is assumed.
    5658     pn_crit_sa   = 1.0_wp
    5659     pn_crit_ocnv = 1.0_wp
    5660     pk_sa        = 1.0_wp
    5661     pk_ocnv      = 1.0_wp
    5662     pd_crit      = 1.5E-9_wp   ! (m)
    5663 
    5664  END SUBROUTINE sumnucl
    5665 
    5666 !------------------------------------------------------------------------------!
    5667 ! Description:
    5668 ! ------------
    5669 !> Conciders both the organic vapor and H2SO4 in nucleation - heteromolecular
    5670 !> nucleation.
    5671 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.
    5672 !------------------------------------------------------------------------------!
    5673  SUBROUTINE hetnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )
    5674 
    5675     IMPLICIT NONE
    5676 
    5677     REAL(wp) ::  z_k_het = 4.1E-14_wp  !< (cm3/s) (Paasonen et al. Table 4: median)
    5678 
    5679     REAL(wp), INTENT(in) ::  pc_org   !< organic vapour concentration (#/m3)
    5680     REAL(wp), INTENT(in) ::  pc_sa    !< H2SO4 conc. (#/m3)
    5681 
    5682     REAL(wp), INTENT(out) ::  pd_crit  !< critical diameter of clusters (m)
    5683     REAL(wp), INTENT(out) ::  pk_ocnv  !< if pk_ocnv = 1, organic compounds participate in nucleation
    5684     REAL(wp), INTENT(out) ::  pk_sa    !< if pk_sa = 1, H2SO4 participate in nucleation
    5685     REAL(wp), INTENT(out) ::  pn_crit_ocnv  !< number of organic molecules in cluster (#)
    5686     REAL(wp), INTENT(out) ::  pn_crit_sa    !< number of H2SO4 molecules in cluster (#)
    5687     REAL(wp), INTENT(out) ::  pnuc_rate     !< nucl. rate (#/(m3 s))
    5688 !
    5689 !-- Nucleation rate (#/m3/s)
    5690     pnuc_rate = z_k_het * pc_sa * pc_org * 1.0E6_wp
    5691 !
    5692 !-- Both organic compounds and H2SO4 are involved when heteromolecular
    5693 !-- nucleation is assumed.
    5694     pn_crit_sa   = 1.0_wp
    5695     pn_crit_ocnv = 1.0_wp
    5696     pk_sa        = 1.0_wp
    5697     pk_ocnv      = 1.0_wp
    5698     pd_crit      = 1.5E-9_wp   ! (m)
    5699 
    5700  END SUBROUTINE hetnucl
    5701 
    5702 !------------------------------------------------------------------------------!
    5703 ! Description:
    5704 ! ------------
    5705 !> Takes into account the homomolecular nucleation of sulphuric acid H2SO4 with
    5706 !> both of the available vapours.
    5707 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.
    5708 !------------------------------------------------------------------------------!
    5709  SUBROUTINE SAnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )
    5710 
    5711     IMPLICIT NONE
    5712 
    5713     REAL(wp) ::  z_k_sa1 = 1.1E-14_wp  !< (cm3/s)
    5714     REAL(wp) ::  z_k_sa2 = 3.2E-14_wp  !< (cm3/s) (Paasonen et al. Table 3.)
    5715 
    5716     REAL(wp), INTENT(in) ::  pc_org   !< organic vapour concentration (#/m3)
    5717     REAL(wp), INTENT(in) ::  pc_sa    !< H2SO4 conc. (#/m3)
    5718 
    5719     REAL(wp), INTENT(out) ::  pd_crit  !< critical diameter of clusters (m)
    5720     REAL(wp), INTENT(out) ::  pk_ocnv  !< if pk_ocnv = 1, organic compounds participate nucleation
    5721     REAL(wp), INTENT(out) ::  pk_sa    !< if pk_sa = 1, H2SO4 participate in nucleation
    5722     REAL(wp), INTENT(out) ::  pn_crit_ocnv  !< number of organic molecules in cluster (#)
    5723     REAL(wp), INTENT(out) ::  pn_crit_sa    !< number of H2SO4 molecules in cluster (#)
    5724     REAL(wp), INTENT(out) ::  pnuc_rate     !< nucleation rate (#/(m3 s))
    5725 !
    5726 !-- Nucleation rate (#/m3/s)
    5727     pnuc_rate = ( z_k_sa1 * pc_sa**2 + z_k_sa2 * pc_sa * pc_org ) * 1.0E+6_wp
    5728 !
    5729 !-- Both organic compounds and H2SO4 are involved when SAnucleation is assumed.
    5730     pn_crit_sa   = 3.0_wp
    5731     pn_crit_ocnv = 1.0_wp
    5732     pk_sa        = 1.0_wp
    5733     pk_ocnv      = 1.0_wp
    5734     pd_crit      = 1.5E-9_wp   ! (m)
    5735 
    5736  END SUBROUTINE SAnucl
    5737 
    5738 !------------------------------------------------------------------------------!
    5739 ! Description:
    5740 ! ------------
    5741 !> Takes into account the homomolecular nucleation of both sulphuric acid and
    5742 !> Lorganic with heteromolecular nucleation.
    5743 !> See Paasonen et al. (2010), Atmos. Chem. Phys., 10, 11223-11242.
    5744 !------------------------------------------------------------------------------!
    5745  SUBROUTINE SAORGnucl( pc_sa, pc_org, pnuc_rate, pd_crit, pn_crit_sa, pn_crit_ocnv, pk_sa, pk_ocnv )
    5746 
    5747     IMPLICIT NONE
    5748 
    5749     REAL(wp) ::  z_k_s1 = 1.4E-14_wp    !< (cm3/s])
    5750     REAL(wp) ::  z_k_s2 = 2.6E-14_wp    !< (cm3/s])
    5751     REAL(wp) ::  z_k_s3 = 0.037E-14_wp  !< (cm3/s]) (Paasonen et al. Table 3.)
    5752 
    5753     REAL(wp), INTENT(in) ::  pc_org   !< organic vapour concentration (#/m3)
    5754     REAL(wp), INTENT(in) ::  pc_sa    !< H2SO4 conc. (#/m3)
    5755 
    5756     REAL(wp), INTENT(out) ::  pd_crit  !< critical diameter of clusters (m)
    5757     REAL(wp), INTENT(out) ::  pk_ocnv  !< if pk_ocnv = 1, organic compounds participate in nucleation
    5758     REAL(wp), INTENT(out) ::  pk_sa    !< if pk_sa = 1, H2SO4 participate in nucleation
    5759     REAL(wp), INTENT(out) ::  pn_crit_ocnv  !< number of organic molecules in cluster (#)
    5760     REAL(wp), INTENT(out) ::  pn_crit_sa    !< number of H2SO4 molecules in cluster (#)
    5761     REAL(wp), INTENT(out) ::  pnuc_rate     !< nucl. rate (#/(m3 s))
    5762 !
    5763 !-- Nucleation rate (#/m3/s)
    5764     pnuc_rate = ( z_k_s1 * pc_sa**2 + z_k_s2 * pc_sa * pc_org + z_k_s3 * pc_org**2 ) * 1.0E+6_wp
    5765 !
    5766 !-- Organic compounds not involved when kinetic nucleation is assumed.
    5767     pn_crit_sa   = 3.0_wp
    5768     pn_crit_ocnv = 3.0_wp
    5769     pk_sa        = 1.0_wp
    5770     pk_ocnv      = 1.0_wp
    5771     pd_crit      = 1.5E-9_wp   ! (m)
    5772 
    5773  END SUBROUTINE SAORGnucl
    57745630
    57755631!------------------------------------------------------------------------------!
     
    82378093
    82388094          DO  ib = 1, nbins_aerosol
    8239              aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:) +              &
     8095             aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:) +           &
    82408096                                                    bc_an_t_val(ib) * dzu(nzt+1)
    82418097             DO  ic = 1, ncomponents_mass
    82428098                icc = ( ic - 1 ) * nbins_aerosol + ib
    8243                 aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:) +             &
     8099                aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:) +          &
    82448100                                                      bc_am_t_val(icc) * dzu(nzt+1)
    82458101             ENDDO
     
    82478103          IF ( .NOT. salsa_gases_from_chem )  THEN
    82488104             DO  ig = 1, ngases_salsa
    8249                 salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:) +                     &
     8105                salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:) +                  &
    82508106                                                  bc_gt_t_val(ig) * dzu(nzt+1)
    82518107             ENDDO
     
    85508406        ONLY:  check_existence, close_input_file, get_attribute, get_variable,                     &
    85518407               inquire_num_variables, inquire_variable_names,                                      &
    8552                netcdf_data_input_get_dimension_length, open_read_file
     8408               netcdf_data_input_get_dimension_length, open_read_file, street_type_f
    85538409
    85548410    USE surface_mod,                                                                               &
     
    85618417    CHARACTER(LEN=25) ::  mod_name             !< name in the input file
    85628418
     8419    INTEGER(iwp) ::  i         !< loop index
    85638420    INTEGER(iwp) ::  ib        !< loop index: aerosol number bins
    85648421    INTEGER(iwp) ::  ic        !< loop index: aerosol chemical components
     
    85668423    INTEGER(iwp) ::  in        !< loop index: emission category
    85678424    INTEGER(iwp) ::  inn       !< loop index
     8425    INTEGER(iwp) ::  j         !< loop index
    85688426    INTEGER(iwp) ::  ss        !< loop index
    85698427
     
    85798437
    85808438!
    8581 !-- Set source arrays to zero:
    8582     DO  ib = 1, nbins_aerosol
    8583        aerosol_number(ib)%source = 0.0_wp
    8584     ENDDO
    8585 
    8586     DO  ic = 1, ncomponents_mass * nbins_aerosol
    8587        aerosol_mass(ic)%source = 0.0_wp
    8588     ENDDO
    8589 
    8590 !
    85918439!-- Define emissions:
    8592 
    85938440    SELECT CASE ( salsa_emission_mode )
    85948441
    8595        CASE ( 'uniform' )
     8442       CASE ( 'uniform', 'parameterized' )
    85968443
    85978444          IF ( init )  THEN  ! Do only once
     
    86058452             CALL size_distribution( surface_aerosol_flux, aerosol_flux_dpg, aerosol_flux_sigmag,  &
    86068453                                     nsect_emission )
    8607              DO  ib = 1, nbins_aerosol
    8608                 source_array(:,:,ib) = nsect_emission(ib)
    8609              ENDDO
     8454             IF ( salsa_emission_mode == 'uniform' )  THEN
     8455                DO  ib = 1, nbins_aerosol
     8456                   source_array(:,:,ib) = nsect_emission(ib)
     8457                ENDDO
     8458             ELSE
     8459                IF ( street_type_f%from_file )  THEN
     8460                   DO  i = nxl, nxr
     8461                      DO  j = nys, nyn
     8462                         IF ( street_type_f%var(j,i) >= main_street_id  .AND.                      &
     8463                              street_type_f%var(j,i) < max_street_id )  THEN
     8464                            source_array(j,i,:) = nsect_emission(:) * emiss_factor_main
     8465                         ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.                  &
     8466                                  street_type_f%var(j,i) < main_street_id )  THEN
     8467                            source_array(j,i,:) = nsect_emission(:) * emiss_factor_side
     8468                         ENDIF
     8469                      ENDDO
     8470                   ENDDO
     8471                ELSE
     8472                   WRITE( message_string, * ) 'salsa_emission_mode = "parameterized" but the '//  &
     8473                                              'street_type data is missing.'
     8474                   CALL message( 'salsa_emission_setup', 'PA0661', 1, 2, 0, 6, 0 )
     8475                ENDIF
     8476             ENDIF
    86108477!
    86118478!--          Check which chemical components are used
     
    86228489             aerosol_flux_mass_fracs_a = aerosol_flux_mass_fracs_a /                               &
    86238490                                         SUM( aerosol_flux_mass_fracs_a(1:ncc ) )
    8624 !
    8625 !--          Set uniform fluxes of default horizontal surfaces
    8626              CALL set_flux( surf_def_h(0), cc_i2m, aerosol_flux_mass_fracs_a, source_array )
     8491             IF ( salsa_emission_mode ==  'uniform' )  THEN
     8492!
     8493!--             Set uniform fluxes of default horizontal surfaces
     8494                CALL set_flux( surf_def_h(0), cc_i2m, aerosol_flux_mass_fracs_a, source_array )
     8495             ELSE
     8496!
     8497!--             Set fluxes normalised based on the street type on land surfaces
     8498                CALL set_flux( surf_lsm_h, cc_i2m, aerosol_flux_mass_fracs_a, source_array )
     8499             ENDIF
    86278500
    86288501             DEALLOCATE( nsect_emission, source_array )
    86298502          ENDIF
    8630 
    8631        CASE ( 'parameterized' )
    8632 !
    8633 !--       TO DO
    86348503
    86358504       CASE ( 'read_from_file' )
     
    86438512          surf_usm_h%amsws = 0.0_wp
    86448513
     8514!
     8515!--       Reset source arrays:
     8516          DO  ib = 1, nbins_aerosol
     8517             aerosol_number(ib)%source = 0.0_wp
     8518          ENDDO
     8519
     8520          DO  ic = 1, ncomponents_mass * nbins_aerosol
     8521             aerosol_mass(ic)%source = 0.0_wp
     8522          ENDDO
     8523
    86458524#if defined( __netcdf )
    86468525!
     
    86588537          IF ( init )  THEN
    86598538!
     8539!--          Variable names
     8540             CALL inquire_num_variables( id_salsa, aero_emission_att%num_vars )
     8541             ALLOCATE( aero_emission_att%var_names(1:aero_emission_att%num_vars) )
     8542             CALL inquire_variable_names( id_salsa, aero_emission_att%var_names )
     8543!
    86608544!--          Read the index and name of chemical components
    86618545             CALL netcdf_data_input_get_dimension_length( id_salsa, aero_emission_att%ncc,         &
     
    86638547             ALLOCATE( aero_emission_att%cc_index(1:aero_emission_att%ncc) )
    86648548             CALL get_variable( id_salsa, 'composition_index', aero_emission_att%cc_index )
    8665              CALL get_variable( id_salsa, 'composition_name', aero_emission_att%cc_name,           &
    8666                                 aero_emission_att%ncc )
     8549
     8550             IF ( check_existence( aero_emission_att%var_names, 'composition_name' ) )  THEN
     8551                CALL get_variable( id_salsa, 'composition_name', aero_emission_att%cc_name,        &
     8552                                   aero_emission_att%ncc )
     8553             ELSE
     8554                message_string = 'Missing composition_name in ' // TRIM( input_file_salsa )
     8555                CALL message( 'salsa_emission_setup', 'PA0657', 1, 2, 0, 6, 0 )
     8556             ENDIF
    86678557!
    86688558!--          Find the corresponding chemical components in the model
    8669              aero_emission_att%cc_input_to_model = 0
     8559             aero_emission_att%cc_in2mod = 0
    86708560             DO  ic = 1, aero_emission_att%ncc
    86718561                in_name = aero_emission_att%cc_name(ic)
    86728562                SELECT CASE ( TRIM( in_name ) )
    86738563                   CASE ( 'H2SO4', 'h2so4', 'SO4', 'so4' )
    8674                       aero_emission_att%cc_input_to_model(1) = ic
     8564                      aero_emission_att%cc_in2mod(1) = ic
    86758565                   CASE ( 'OC', 'oc', 'organics' )
    8676                       aero_emission_att%cc_input_to_model(2) = ic
     8566                      aero_emission_att%cc_in2mod(2) = ic
    86778567                   CASE ( 'BC', 'bc' )
    8678                       aero_emission_att%cc_input_to_model(3) = ic
     8568                      aero_emission_att%cc_in2mod(3) = ic
    86798569                   CASE ( 'DU', 'du' )
    8680                       aero_emission_att%cc_input_to_model(4) = ic
     8570                      aero_emission_att%cc_in2mod(4) = ic
    86818571                   CASE ( 'SS', 'ss' )
    8682                       aero_emission_att%cc_input_to_model(5) = ic
     8572                      aero_emission_att%cc_in2mod(5) = ic
    86838573                   CASE ( 'HNO3', 'hno3', 'NO', 'no' )
    8684                       aero_emission_att%cc_input_to_model(6) = ic
     8574                      aero_emission_att%cc_in2mod(6) = ic
    86858575                   CASE ( 'NH3', 'nh3', 'NH', 'nh' )
    8686                       aero_emission_att%cc_input_to_model(7) = ic
     8576                      aero_emission_att%cc_in2mod(7) = ic
    86878577                END SELECT
    86888578
    86898579             ENDDO
    86908580
    8691              IF ( SUM( aero_emission_att%cc_input_to_model ) == 0 )  THEN
     8581             IF ( SUM( aero_emission_att%cc_in2mod ) == 0 )  THEN
    86928582                message_string = 'None of the aerosol chemical components in ' // TRIM(            &
    86938583                                 input_file_salsa ) // ' correspond to the ones applied in SALSA.'
     
    87068596             CALL get_attribute( id_salsa, 'lod', aero_emission_att%lod, .FALSE.,                  &
    87078597                                 'aerosol_emission_values' )
    8708 !
    8709 !--          Variable names
    8710              CALL inquire_num_variables( id_salsa, aero_emission_att%num_vars )
    8711              ALLOCATE( aero_emission_att%var_names(1:aero_emission_att%num_vars) )
    8712              CALL inquire_variable_names( id_salsa, aero_emission_att%var_names )
    87138598
    87148599!
     
    87388623!
    87398624!--             Get emission category names and indices
    8740                 CALL get_variable( id_salsa, 'emission_category_name', aero_emission_att%cat_name, &
    8741                                    aero_emission_att%ncat)
     8625                IF ( check_existence( aero_emission_att%var_names, 'emission_category_name' ) )  THEN
     8626                   CALL get_variable( id_salsa, 'emission_category_name',                          &
     8627                                      aero_emission_att%cat_name,  aero_emission_att%ncat )
     8628                ELSE
     8629                   message_string = 'Missing emission_category_name in ' // TRIM( input_file_salsa )
     8630                   CALL message( 'salsa_emission_setup', 'PA0658', 1, 2, 0, 6, 0 )
     8631                ENDIF
    87428632                CALL get_variable( id_salsa, 'emission_category_index', aero_emission_att%cat_index )
    87438633!
     
    87648654!--             For each hour of year:
    87658655                IF ( check_existence( aero_emission_att%var_names, 'nhoursyear' ) )  THEN
    8766                    CALL netcdf_data_input_get_dimension_length( id_salsa, aero_emission_att%nhoursyear,&
    8767                                                                 'nhoursyear' )
     8656                   CALL netcdf_data_input_get_dimension_length( id_salsa,                          &
     8657                                                        aero_emission_att%nhoursyear, 'nhoursyear' )
    87688658                   ALLOCATE( aero_emission_att%etf(1:aero_emission_att%ncat,                       &
    87698659                                                   1:aero_emission_att%nhoursyear) )
     
    87908680!
    87918681!--             Get chemical composition (i.e. mass fraction of different species) in aerosols
    8792                 ALLOCATE( aero_emission%def_mass_fracs(1:aero_emission_att%ncat,                   &
    8793                                                        1:aero_emission_att%ncc) )
    8794                 aero_emission%def_mass_fracs = 0.0_wp
    8795                 CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%def_mass_fracs,  &
    8796                                    0, aero_emission_att%ncc-1, 0, aero_emission_att%ncat-1 )
     8682                IF ( check_existence( aero_emission_att%var_names, 'emission_mass_fracs' ) )  THEN
     8683                   ALLOCATE( aero_emission%def_mass_fracs(1:aero_emission_att%ncat,                &
     8684                                                          1:aero_emission_att%ncc) )
     8685                   aero_emission%def_mass_fracs = 0.0_wp
     8686                   CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%def_mass_fracs,&
     8687                                      0, aero_emission_att%ncc-1, 0, aero_emission_att%ncat-1 )
     8688                ELSE
     8689                   message_string = 'Missing emission_mass_fracs in ' //  TRIM( input_file_salsa )
     8690                   CALL message( 'salsa_emission_setup', 'PA0659', 1, 2, 0, 6, 0 )
     8691                ENDIF
    87978692!
    87988693!--             If the chemical component is not activated, set its mass fraction to 0 to avoid
    87998694!--             inbalance between number and mass flux
    8800                 cc_i2m = aero_emission_att%cc_input_to_model
    8801                 IF ( index_so4 < 0  .AND.  cc_i2m(1) /= 0 )                                        &
     8695                cc_i2m = aero_emission_att%cc_in2mod
     8696                IF ( index_so4 < 0  .AND.  cc_i2m(1) > 0 )                                         &
    88028697                                                  aero_emission%def_mass_fracs(:,cc_i2m(1)) = 0.0_wp
    8803                 IF ( index_oc  < 0  .AND.  cc_i2m(2) /= 0 )                                        &
     8698                IF ( index_oc  < 0  .AND.  cc_i2m(2) > 0 )                                         &
    88048699                                                  aero_emission%def_mass_fracs(:,cc_i2m(2)) = 0.0_wp
    8805                 IF ( index_bc  < 0  .AND.  cc_i2m(3) /= 0 )                                        &
     8700                IF ( index_bc  < 0  .AND.  cc_i2m(3) > 0 )                                         &
    88068701                                                  aero_emission%def_mass_fracs(:,cc_i2m(3)) = 0.0_wp
    8807                 IF ( index_du  < 0  .AND.  cc_i2m(4) /= 0 )                                        &
     8702                IF ( index_du  < 0  .AND.  cc_i2m(4) > 0 )                                         &
    88088703                                                  aero_emission%def_mass_fracs(:,cc_i2m(4)) = 0.0_wp
    8809                 IF ( index_ss  < 0  .AND.  cc_i2m(5) /= 0 )                                        &
     8704                IF ( index_ss  < 0  .AND.  cc_i2m(5) > 0 )                                         &
    88108705                                                  aero_emission%def_mass_fracs(:,cc_i2m(5)) = 0.0_wp
    8811                 IF ( index_no  < 0  .AND.  cc_i2m(6) /= 0 )                                        &
     8706                IF ( index_no  < 0  .AND.  cc_i2m(6) > 0 )                                         &
    88128707                                                  aero_emission%def_mass_fracs(:,cc_i2m(6)) = 0.0_wp
    8813                 IF ( index_nh  < 0  .AND.  cc_i2m(7) /= 0 )                                        &
     8708                IF ( index_nh  < 0  .AND.  cc_i2m(7) > 0 )                                         &
    88148709                                                  aero_emission%def_mass_fracs(:,cc_i2m(7)) = 0.0_wp
    88158710!
     
    88868781!
    88878782!--             Read time stamps:
    8888                 CALL get_variable( id_salsa, 'time', aero_emission_att%time )
     8783                IF ( check_existence( aero_emission_att%var_names, 'time' ) )  THEN
     8784                   CALL get_variable( id_salsa, 'time', aero_emission_att%time )
     8785                ELSE
     8786                   message_string = 'Missing time in ' //  TRIM( input_file_salsa )
     8787                   CALL message( 'salsa_emission_setup', 'PA0660', 1, 2, 0, 6, 0 )
     8788                ENDIF
    88898789!
    88908790!--             Read emission mass fractions
    8891                 CALL get_variable( id_salsa, 'emission_mass_fracs', aero_emission%preproc_mass_fracs )
     8791                IF ( check_existence( aero_emission_att%var_names, 'emission_mass_fracs' ) )  THEN
     8792                   CALL get_variable( id_salsa, 'emission_mass_fracs',                             &
     8793                                      aero_emission%preproc_mass_fracs )
     8794                ELSE
     8795                   message_string = 'Missing emission_mass_fracs in ' //  TRIM( input_file_salsa )
     8796                   CALL message( 'salsa_emission_setup', 'PA0659', 1, 2, 0, 6, 0 )
     8797                ENDIF
    88928798!
    88938799!--             If the chemical component is not activated, set its mass fraction to 0
    8894                 cc_i2m = aero_emission_att%cc_input_to_model
     8800                cc_i2m = aero_emission_att%cc_in2mod
    88958801                IF ( index_so4 < 0  .AND.  cc_i2m(1) /= 0 )                                        &
    88968802                   aero_emission%preproc_mass_fracs(cc_i2m(1)) = 0.0_wp
     
    89708876!--             only for either default, land or urban surface.
    89718877                IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
    8972                    CALL set_flux( surf_def_h(0), aero_emission_att%cc_input_to_model,              &
     8878                   CALL set_flux( surf_def_h(0), aero_emission_att%cc_in2mod,                      &
    89738879                                  aero_emission%def_mass_fracs(in,:), source_array )
    89748880                ELSE
    8975                    CALL set_flux( surf_usm_h, aero_emission_att%cc_input_to_model,                 &
     8881                   CALL set_flux( surf_usm_h, aero_emission_att%cc_in2mod,                         &
    89768882                                  aero_emission%def_mass_fracs(in,:), source_array )
    8977                    CALL set_flux( surf_lsm_h, aero_emission_att%cc_input_to_model,                 &
     8883                   CALL set_flux( surf_lsm_h, aero_emission_att%cc_in2mod,                         &
    89788884                                  aero_emission%def_mass_fracs(in,:), source_array )
    89798885                ENDIF
     
    90048910!--          for either default, land and urban surface.
    90058911             IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
    9006                 CALL set_flux( surf_def_h(0), aero_emission_att%cc_input_to_model,                 &
     8912                CALL set_flux( surf_def_h(0), aero_emission_att%cc_in2mod,                         &
    90078913                               aero_emission%preproc_mass_fracs, aero_emission%preproc_data )
    90088914             ELSE
    9009                 CALL set_flux( surf_usm_h, aero_emission_att%cc_input_to_model,                    &
     8915                CALL set_flux( surf_usm_h, aero_emission_att%cc_in2mod,                            &
    90108916                               aero_emission%preproc_mass_fracs, aero_emission%preproc_data )
    9011                 CALL set_flux( surf_lsm_h, aero_emission_att%cc_input_to_model,                    &
     8917                CALL set_flux( surf_lsm_h, aero_emission_att%cc_in2mod,                            &
    90128918                               aero_emission%preproc_mass_fracs, aero_emission%preproc_data )
    90138919             ENDIF
     
    92499155
    92509156    INTEGER(iwp) ::  id_chem        !< NetCDF id of chemistry emission file
     9157    INTEGER(iwp) ::  i              !< loop index
    92519158    INTEGER(iwp) ::  ig             !< loop index
    92529159    INTEGER(iwp) ::  in             !< running index for emission categories
     9160    INTEGER(iwp) ::  j              !< loop index
    92539161    INTEGER(iwp) ::  num_vars       !< number of variables
    92549162
     
    92589166
    92599167    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_factor  !< emission time factor
     9168
     9169    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dum_var_3d  !<
     9170
     9171    REAL(wp), DIMENSION(:,:,:,:,:), ALLOCATABLE ::  dum_var_5d  !<
    92609172
    92619173!
     
    92809192!
    92819193!--    Read the index and name of chemical components
    9282        CALL netcdf_data_input_get_dimension_length( id_chem, chem_emission_att%nspec,              &
     9194       CALL netcdf_data_input_get_dimension_length( id_chem, chem_emission_att%n_emiss_species,    &
    92839195                                                    'nspecies' )
    9284        ALLOCATE( chem_emission_att%species_index(1:chem_emission_att%nspec) )
     9196       ALLOCATE( chem_emission_att%species_index(1:chem_emission_att%n_emiss_species) )
    92859197       CALL get_variable( id_chem, 'emission_index', chem_emission_att%species_index )
    92869198       CALL get_variable( id_chem, 'emission_name', chem_emission_att%species_name,                &
    9287                           chem_emission_att%nspec )
     9199                          chem_emission_att%n_emiss_species )
     9200!
     9201!--    Allocate emission data
     9202       ALLOCATE( chem_emission(1:chem_emission_att%n_emiss_species) )
    92889203!
    92899204!--    Find the corresponding indices in the model
    92909205       emission_index_chem = 0
    9291        DO  ig = 1, chem_emission_att%nspec
     9206       DO  ig = 1, chem_emission_att%n_emiss_species
    92929207          in_name = chem_emission_att%species_name(ig)
    92939208          SELECT CASE ( TRIM( in_name ) )
     
    93669281!--       Allocate and read surface emission data (in total PM) (NOTE that "preprocessed" input data
    93679282!--       array is applied now here)
    9368           ALLOCATE( chem_emission%preproc_emission_data(nys:nyn,nxl:nxr, 1:chem_emission_att%nspec,&
    9369                                                         1:chem_emission_att%ncat) )
    9370           CALL get_variable( id_chem, 'emission_values', chem_emission%preproc_emission_data,      &
    9371                              0, chem_emission_att%ncat-1, 0, chem_emission_att%nspec-1,            &
    9372                              nxl, nxr, nys, nyn )
     9283          ALLOCATE( dum_var_5d(1,nys:nyn,nxl:nxr,1:chem_emission_att%n_emiss_species,              &
     9284                               1:chem_emission_att%ncat) )
     9285          CALL get_variable( id_chem, 'emission_values', dum_var_5d, 0, chem_emission_att%ncat-1,  &
     9286                             0, chem_emission_att%n_emiss_species-1, nxl, nxr, nys, nyn, 0, 0 )
     9287          DO  ig = 1, chem_emission_att%n_emiss_species
     9288             ALLOCATE( chem_emission(ig)%default_emission_data(nys:nyn,nxl:nxr,                    &
     9289                                                               1:chem_emission_att%ncat) )
     9290             DO  in = 1, chem_emission_att%ncat
     9291                DO  i = nxl, nxr
     9292                   DO  j = nys, nyn
     9293                      chem_emission(ig)%default_emission_data(j,i,in) = dum_var_5d(1,j,i,ig,in)
     9294                   ENDDO
     9295                ENDDO
     9296             ENDDO
     9297          ENDDO
     9298          DEALLOCATE( dum_var_5d )
    93739299!
    93749300!--    Pre-processed mode:
     
    94129338!
    94139339!--    Set gas emissions for each emission category
     9340       ALLOCATE( dum_var_3d(nys:nyn,nxl:nxr,1:chem_emission_att%n_emiss_species) )
     9341
    94149342       DO  in = 1, chem_emission_att%ncat
     9343          DO  ig = 1, chem_emission_att%n_emiss_species
     9344             dum_var_3d(:,:,ig) = chem_emission(ig)%default_emission_data(:,:,in)
     9345          ENDDO
    94159346!
    94169347!--       Set surface fluxes only for either default, land or urban surface
    94179348          IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
    94189349             CALL set_gas_flux( surf_def_h(0), emission_index_chem, chem_emission_att%units,       &
    9419                                 chem_emission%preproc_emission_data(:,:,:,in), time_factor(in) )
     9350                                dum_var_3d(:,:,in), time_factor(in) )
    94209351          ELSE
    94219352             CALL set_gas_flux( surf_usm_h, emission_index_chem, chem_emission_att%units,          &
    9422                                 chem_emission%preproc_emission_data(:,:,:,in), time_factor(in) )
     9353                                dum_var_3d(:,:,in), time_factor(in) )
    94239354             CALL set_gas_flux( surf_lsm_h, emission_index_chem, chem_emission_att%units,          &
    9424                                 chem_emission%preproc_emission_data(:,:,:,in), time_factor(in) )
     9355                                dum_var_3d(:,:,in), time_factor(in) )
    94259356          ENDIF
    94269357       ENDDO
     9358       DEALLOCATE( dum_var_3d )
    94279359!
    94289360!--    The next emission update is again after one hour
     
    94389370!--    Allocate the data input array always before reading in the data and deallocate after (NOTE
    94399371!--    that "preprocessed" input data array is applied now here)
    9440        ALLOCATE( chem_emission%default_emission_data(nys:nyn,nxl:nxr,1:nbins_aerosol) )
     9372       ALLOCATE( dum_var_5d(1,1,nys:nyn,nxl:nxr,1:chem_emission_att%n_emiss_species) )
    94419373!
    94429374!--    Read in the next time step
    9443        CALL get_variable( id_chem, 'emission_values', chem_emission%default_emission_data,         &
    9444                           chem_emission_att%i_hour, 0, chem_emission_att%nspec-1, nxl, nxr, nys, nyn )
     9375       CALL get_variable( id_chem, 'emission_values', dum_var_5d,                                  &
     9376                          0, chem_emission_att%n_emiss_species-1, nxl, nxr, nys, nyn, 0, 0,        &
     9377                          chem_emission_att%i_hour, chem_emission_att%i_hour )
    94459378!
    94469379!--    Set surface fluxes only for either default, land or urban surface
    94479380       IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
    94489381          CALL set_gas_flux( surf_def_h(0), emission_index_chem, chem_emission_att%units,          &
    9449                              chem_emission%default_emission_data )
     9382                             dum_var_5d(1,1,:,:,:) )
    94509383       ELSE
    94519384          CALL set_gas_flux( surf_usm_h, emission_index_chem, chem_emission_att%units,             &
    9452                              chem_emission%default_emission_data )
     9385                             dum_var_5d(1,1,:,:,:) )
    94539386          CALL set_gas_flux( surf_lsm_h, emission_index_chem, chem_emission_att%units,             &
    9454                              chem_emission%default_emission_data )
    9455        ENDIF
     9387                             dum_var_5d(1,1,:,:,:) )
     9388       ENDIF
     9389       DEALLOCATE ( dum_var_5d )
    94569390!
    94579391!--    Determine the next emission update
    94589392       next_gas_emission_update = gas_emission_time(chem_emission_att%i_hour+2)
    94599393
    9460        DEALLOCATE( chem_emission%default_emission_data )
    94619394    ENDIF
    94629395!
     
    95069439       REAL(wp), DIMENSION(ngases_salsa) ::  conv     !< unit conversion factor
    95079440
    9508        REAL(wp), DIMENSION(nys:nyn,nxl:nxr,chem_emission_att%nspec), INTENT(in) ::  source_array  !<
     9441       REAL(wp), DIMENSION(nys:nyn,nxl:nxr,chem_emission_att%n_emiss_species), INTENT(in) ::  source_array  !<
    95099442
    95109443       TYPE(surf_type), INTENT(inout) :: surface  !< respective surface type
     
    95819514    CHARACTER(LEN=*) ::  var      !<
    95829515
    9583     SELECT CASE ( TRIM( var ) )
    9584 
    9585        CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
    9586           IF (  .NOT.  salsa )  THEN
    9587              message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
    9588              CALL message( 'check_parameters', 'PA0652', 1, 2, 0, 6, 0 )
    9589           ENDIF
    9590           IF (  salsa_gases_from_chem )  THEN
    9591              message_string = 'gases are imported from the chemistry module and thus output of "'  &
    9592                                // TRIM( var ) // '" is not allowed'
    9593              CALL message( 'check_parameters', 'PA0653', 1, 2, 0, 6, 0 )
    9594           ENDIF
     9516    INTEGER(iwp) ::  char_to_int   !< for converting character to integer
     9517
     9518    IF ( var(1:6) /= 'salsa_' )  THEN
     9519       unit = 'illegal'
     9520       RETURN
     9521    ENDIF
     9522!
     9523!-- Treat bin-specific outputs separately
     9524    IF ( var(7:11) ==  'N_bin' )  THEN
     9525       READ( var(12:),* ) char_to_int
     9526       IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
    95959527          unit = '#/m3'
    9596 
    9597        CASE ( 'LDSA' )
    9598           IF (  .NOT.  salsa )  THEN
    9599              message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
    9600              CALL message( 'check_parameters', 'PA0646', 1, 2, 0, 6, 0 )
    9601           ENDIF
    9602           unit = 'mum2/cm3'
    9603 
    9604        CASE ( 'm_bin1', 'm_bin2',  'm_bin3',  'm_bin4',  'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8',   &
    9605               'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12', 'PM2.5',  'PM10',   's_BC',   's_DU',     &
    9606               's_H2O',  's_NH',    's_NO',    's_OC',    's_SO4',  's_SS' )
    9607           IF (  .NOT.  salsa )  THEN
    9608              message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
    9609              CALL message( 'check_parameters', 'PA0647', 1, 2, 0, 6, 0 )
    9610           ENDIF
     9528       ELSE
     9529          unit = 'illegal'
     9530          RETURN
     9531       ENDIF
     9532
     9533    ELSEIF ( var(7:11) ==  'm_bin' )  THEN
     9534       READ( var(12:),* ) char_to_int
     9535       IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
    96119536          unit = 'kg/m3'
    9612 
    9613        CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4',  'N_bin5',  'N_bin6', 'N_bin7', 'N_bin8',    &
    9614               'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12', 'Ntot' )
    9615           IF (  .NOT.  salsa )  THEN
    9616              message_string = 'output of "' // TRIM( var ) // '" requires salsa = .TRUE.'
    9617              CALL message( 'check_parameters', 'PA0645', 1, 2, 0, 6, 0 )
    9618           ENDIF
    9619           unit = '#/m3'
     9537       ELSE
     9538          unit = 'illegal'
     9539          RETURN
     9540       ENDIF
     9541
     9542    ELSE
     9543       SELECT CASE ( TRIM( var(7:) ) )
     9544
     9545          CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
     9546             IF (  salsa_gases_from_chem )  THEN
     9547                message_string = 'gases are imported from the chemistry module and thus output '// &
     9548                                 'of "' // TRIM( var ) // '" is not allowed'
     9549                CALL message( 'check_parameters', 'PA0653', 1, 2, 0, 6, 0 )
     9550             ENDIF
     9551             unit = '#/m3'
     9552
     9553          CASE ( 'LDSA' )
     9554             unit = 'mum2/cm3'
     9555
     9556          CASE ( 'PM0.1', 'PM2.5', 'PM10', 's_BC', 's_DU', 's_H2O', 's_NH', 's_NO', 's_OC',        &
     9557                 's_SO4', 's_SS' )
     9558             unit = 'kg/m3'
     9559
     9560          CASE ( 'N_UFP', 'Ntot' )
     9561             unit = '#/m3'
     9562
     9563          CASE DEFAULT
     9564             unit = 'illegal'
     9565
     9566       END SELECT
     9567    ENDIF
     9568
     9569 END SUBROUTINE salsa_check_data_output
     9570
     9571!------------------------------------------------------------------------------!
     9572! Description:
     9573! ------------
     9574!> Check profile data output for salsa. Currently only for diagnostic variables
     9575!> Ntot, N_UFP, PM0.1, PM2.5, PM10 and LDSA
     9576!------------------------------------------------------------------------------!
     9577 SUBROUTINE salsa_check_data_output_pr( var, var_count, unit, dopr_unit )
     9578
     9579    USE arrays_3d,                                                                                 &
     9580        ONLY: zu
     9581
     9582    USE profil_parameter,                                                                          &
     9583        ONLY:  dopr_index
     9584
     9585    USE statistics,                                                                                &
     9586        ONLY:  hom, pr_palm, statistic_regions
     9587
     9588    IMPLICIT NONE
     9589
     9590    CHARACTER(LEN=*) ::  dopr_unit  !<
     9591    CHARACTER(LEN=*) ::  unit       !<
     9592    CHARACTER(LEN=*) ::  var        !<
     9593
     9594    INTEGER(iwp) ::  var_count     !<
     9595
     9596    IF ( var(1:6) /= 'salsa_' )  THEN
     9597       unit = 'illegal'
     9598       RETURN
     9599    ENDIF
     9600
     9601    SELECT CASE ( TRIM( var(7:) ) )
     9602
     9603       CASE( 'LDSA' )
     9604          salsa_pr_count = salsa_pr_count + 1
     9605          salsa_pr_index(salsa_pr_count) = 1
     9606          dopr_index(var_count) = pr_palm + salsa_pr_count
     9607          dopr_unit = 'mum2/cm3'
     9608          unit = dopr_unit
     9609          hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
     9610
     9611       CASE( 'N_UFP' )
     9612          salsa_pr_count = salsa_pr_count + 1
     9613          salsa_pr_index(salsa_pr_count) = 2
     9614          dopr_index(var_count) = pr_palm + salsa_pr_count
     9615          dopr_unit = '#/m3'
     9616          unit = dopr_unit
     9617          hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
     9618
     9619       CASE( 'Ntot' )
     9620          salsa_pr_count = salsa_pr_count + 1
     9621          salsa_pr_index(salsa_pr_count) = 3
     9622          dopr_index(var_count) = pr_palm + salsa_pr_count
     9623          dopr_unit = '#/m3'
     9624          unit = dopr_unit
     9625          hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
     9626
     9627       CASE( 'PM0.1' )
     9628          salsa_pr_count = salsa_pr_count + 1
     9629          salsa_pr_index(salsa_pr_count) = 4
     9630          dopr_index(var_count) = pr_palm + salsa_pr_count
     9631          dopr_unit = 'kg/m3'
     9632          unit = dopr_unit
     9633          hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
     9634
     9635       CASE( 'PM2.5' )
     9636          salsa_pr_count = salsa_pr_count + 1
     9637          salsa_pr_index(salsa_pr_count) = 5
     9638          dopr_index(var_count) = pr_palm + salsa_pr_count
     9639          dopr_unit = 'kg/m3'
     9640          unit = dopr_unit
     9641          hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
     9642
     9643       CASE( 'PM10' )
     9644          salsa_pr_count = salsa_pr_count + 1
     9645          salsa_pr_index(salsa_pr_count) = 6
     9646          dopr_index(var_count) = pr_palm + salsa_pr_count
     9647          dopr_unit = 'kg/m3'
     9648          unit = dopr_unit
     9649          hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
    96209650
    96219651       CASE DEFAULT
     
    96249654    END SELECT
    96259655
    9626  END SUBROUTINE salsa_check_data_output
     9656
     9657 END SUBROUTINE salsa_check_data_output_pr
     9658
     9659!-------------------------------------------------------------------------------!
     9660!> Description:
     9661!> Calculation of horizontally averaged profiles for salsa.
     9662!-------------------------------------------------------------------------------!
     9663 SUBROUTINE salsa_statistics( mode, sr, tn )
     9664
     9665    USE control_parameters,                                                                        &
     9666        ONLY:  max_pr_user
     9667
     9668    USE chem_modules,                                                                              &
     9669        ONLY:  max_pr_cs
     9670
     9671    USE statistics,                                                                                &
     9672        ONLY:  pr_palm, rmask, sums_l
     9673
     9674    IMPLICIT NONE
     9675
     9676    CHARACTER(LEN=*) ::  mode  !<
     9677
     9678    INTEGER(iwp) ::  i    !< loop index
     9679    INTEGER(iwp) ::  ib   !< loop index
     9680    INTEGER(iwp) ::  ic   !< loop index
     9681    INTEGER(iwp) ::  ii   !< loop index
     9682    INTEGER(iwp) ::  ind  !< index in the statistical output
     9683    INTEGER(iwp) ::  j    !< loop index
     9684    INTEGER(iwp) ::  k    !< loop index
     9685    INTEGER(iwp) ::  sr   !< statistical region
     9686    INTEGER(iwp) ::  tn   !< thread number
     9687
     9688    REAL(wp) ::  df        !< For calculating LDSA: fraction of particles depositing in the alveolar
     9689                           !< (or tracheobronchial) region of the lung. Depends on the particle size
     9690    REAL(wp) ::  mean_d    !< Particle diameter in micrometres
     9691    REAL(wp) ::  temp_bin  !< temporary variable
     9692
     9693    IF ( mode == 'profiles' )  THEN
     9694       !$OMP DO
     9695       DO  ii = 1, salsa_pr_count
     9696
     9697          ind = pr_palm + max_pr_user + max_pr_cs + ii
     9698
     9699          SELECT CASE( salsa_pr_index(ii) )
     9700
     9701             CASE( 1 )  ! LDSA
     9702                DO  i = nxl, nxr
     9703                   DO  j = nys, nyn
     9704                      DO  k = nzb, nzt+1
     9705                         temp_bin = 0.0_wp
     9706                         DO  ib = 1, nbins_aerosol
     9707   !
     9708   !--                      Diameter in micrometres
     9709                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
     9710   !
     9711   !--                      Deposition factor: alveolar
     9712                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
     9713                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
     9714                                   1.362_wp )**2 ) )
     9715   !
     9716   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
     9717                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
     9718                                       aerosol_number(ib)%conc(k,j,i)
     9719                         ENDDO
     9720                         sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr)  *         &
     9721                                           MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     9722                      ENDDO
     9723                   ENDDO
     9724                ENDDO
     9725
     9726             CASE( 2 )  ! N_UFP
     9727                DO  i = nxl, nxr
     9728                   DO  j = nys, nyn
     9729                      DO  k = nzb, nzt+1
     9730                         temp_bin = 0.0_wp
     9731                         DO  ib = 1, nbins_aerosol
     9732                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )                          &
     9733                               temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     9734                         ENDDO
     9735                         sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr)  *         &
     9736                                           MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     9737                      ENDDO
     9738                   ENDDO
     9739                ENDDO
     9740
     9741             CASE( 3 )  ! Ntot
     9742                DO  i = nxl, nxr
     9743                   DO  j = nys, nyn
     9744                      DO  k = nzb, nzt+1
     9745                         temp_bin = 0.0_wp
     9746                         DO  ib = 1, nbins_aerosol
     9747                            temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     9748                         ENDDO
     9749                         sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr)  *         &
     9750                                           MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     9751                      ENDDO
     9752                   ENDDO
     9753                ENDDO
     9754
     9755             CASE( 4 )  ! PM0.1
     9756                DO  i = nxl, nxr
     9757                   DO  j = nys, nyn
     9758                      DO  k = nzb, nzt+1
     9759                         temp_bin = 0.0_wp
     9760                         DO  ib = 1, nbins_aerosol
     9761                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     9762                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     9763                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     9764                               ENDDO
     9765                            ENDIF
     9766                         ENDDO
     9767                         sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr)  *         &
     9768                                           MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     9769                      ENDDO
     9770                   ENDDO
     9771                ENDDO
     9772
     9773             CASE( 5 )  ! PM2.5
     9774                DO  i = nxl, nxr
     9775                   DO  j = nys, nyn
     9776                      DO  k = nzb, nzt+1
     9777                         temp_bin = 0.0_wp
     9778                         DO  ib = 1, nbins_aerosol
     9779                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
     9780                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     9781                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     9782                               ENDDO
     9783                            ENDIF
     9784                         ENDDO
     9785                         sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr)  *         &
     9786                                           MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     9787                      ENDDO
     9788                   ENDDO
     9789                ENDDO
     9790
     9791             CASE( 6 )  ! PM10
     9792                DO  i = nxl, nxr
     9793                   DO  j = nys, nyn
     9794                      DO  k = nzb, nzt+1
     9795                         temp_bin = 0.0_wp
     9796                         DO  ib = 1, nbins_aerosol
     9797                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
     9798                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     9799                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     9800                               ENDDO
     9801                            ENDIF
     9802                         ENDDO
     9803                         sums_l(k,ind,tn) = sums_l(k,ind,tn) + temp_bin * rmask(j,i,sr)  *         &
     9804                                           MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     9805                      ENDDO
     9806                   ENDDO
     9807                ENDDO
     9808
     9809          END SELECT
     9810       ENDDO
     9811
     9812    ELSEIF ( mode == 'time_series' )  THEN
     9813!
     9814!--    TODO
     9815    ENDIF
     9816
     9817 END SUBROUTINE salsa_statistics
     9818
    96279819
    96289820!------------------------------------------------------------------------------!
     
    96439835    CHARACTER(LEN=*)  ::  variable   !<
    96449836
     9837    INTEGER(iwp) ::  char_to_int  !< for converting character to integer
    96459838    INTEGER(iwp) ::  found_index  !<
    96469839    INTEGER(iwp) ::  i            !<
     
    96539846                          !< (or tracheobronchial) region of the lung. Depends on the particle size
    96549847    REAL(wp) ::  mean_d   !< Particle diameter in micrometres
    9655     REAL(wp) ::  nc       !< Particle number concentration in units 1/cm**3
    96569848    REAL(wp) ::  temp_bin !< temporary variable
    96579849
     
    96629854    IF ( mode == 'allocate' )  THEN
    96639855
    9664        SELECT CASE ( TRIM( variable ) )
    9665 
    9666           CASE ( 'g_H2SO4' )
    9667              IF ( .NOT. ALLOCATED( g_h2so4_av ) )  THEN
    9668                 ALLOCATE( g_h2so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9669              ENDIF
    9670              g_h2so4_av = 0.0_wp
    9671 
    9672           CASE ( 'g_HNO3' )
    9673              IF ( .NOT. ALLOCATED( g_hno3_av ) )  THEN
    9674                 ALLOCATE( g_hno3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9675              ENDIF
    9676              g_hno3_av = 0.0_wp
    9677 
    9678           CASE ( 'g_NH3' )
    9679              IF ( .NOT. ALLOCATED( g_nh3_av ) )  THEN
    9680                 ALLOCATE( g_nh3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9681              ENDIF
    9682              g_nh3_av = 0.0_wp
    9683 
    9684           CASE ( 'g_OCNV' )
    9685              IF ( .NOT. ALLOCATED( g_ocnv_av ) )  THEN
    9686                 ALLOCATE( g_ocnv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9687              ENDIF
    9688              g_ocnv_av = 0.0_wp
    9689 
    9690           CASE ( 'g_OCSV' )
    9691              IF ( .NOT. ALLOCATED( g_ocsv_av ) )  THEN
    9692                 ALLOCATE( g_ocsv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9693              ENDIF
    9694              g_ocsv_av = 0.0_wp
    9695 
    9696           CASE ( 'LDSA' )
    9697              IF ( .NOT. ALLOCATED( ldsa_av ) )  THEN
    9698                 ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9699              ENDIF
    9700              ldsa_av = 0.0_wp
    9701 
    9702           CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8',   &
    9703                  'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' )
    9704              IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
    9705                 ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
    9706              ENDIF
    9707              nbins_av = 0.0_wp
    9708 
    9709           CASE ( 'Ntot' )
    9710              IF ( .NOT. ALLOCATED( ntot_av ) )  THEN
    9711                 ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9712              ENDIF
    9713              ntot_av = 0.0_wp
    9714 
    9715           CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8',   &
    9716                  'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12' )
    9717              IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
    9718                 ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
    9719              ENDIF
    9720              mbins_av = 0.0_wp
    9721 
    9722           CASE ( 'PM2.5' )
    9723              IF ( .NOT. ALLOCATED( pm25_av ) )  THEN
    9724                 ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9725              ENDIF
    9726              pm25_av = 0.0_wp
    9727 
    9728           CASE ( 'PM10' )
    9729              IF ( .NOT. ALLOCATED( pm10_av ) )  THEN
    9730                 ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9731              ENDIF
    9732              pm10_av = 0.0_wp
    9733 
    9734           CASE ( 's_BC' )
    9735              IF ( .NOT. ALLOCATED( s_bc_av ) )  THEN
    9736                 ALLOCATE( s_bc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9737              ENDIF
    9738              s_bc_av = 0.0_wp
    9739 
    9740           CASE ( 's_DU' )
    9741              IF ( .NOT. ALLOCATED( s_du_av ) )  THEN
    9742                 ALLOCATE( s_du_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9743              ENDIF
    9744              s_du_av = 0.0_wp
    9745 
    9746           CASE ( 's_H2O' )
    9747              IF ( .NOT. ALLOCATED( s_h2o_av ) )  THEN
    9748                 ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9749              ENDIF
    9750              s_h2o_av = 0.0_wp
    9751 
    9752           CASE ( 's_NH' )
    9753              IF ( .NOT. ALLOCATED( s_nh_av ) )  THEN
    9754                 ALLOCATE( s_nh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9755              ENDIF
    9756              s_nh_av = 0.0_wp
    9757 
    9758           CASE ( 's_NO' )
    9759              IF ( .NOT. ALLOCATED( s_no_av ) )  THEN
    9760                 ALLOCATE( s_no_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9761              ENDIF
    9762              s_no_av = 0.0_wp
    9763 
    9764           CASE ( 's_OC' )
    9765              IF ( .NOT. ALLOCATED( s_oc_av ) )  THEN
    9766                 ALLOCATE( s_oc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9767              ENDIF
    9768              s_oc_av = 0.0_wp
    9769 
    9770           CASE ( 's_SO4' )
    9771              IF ( .NOT. ALLOCATED( s_so4_av ) )  THEN
    9772                 ALLOCATE( s_so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9773              ENDIF
    9774              s_so4_av = 0.0_wp   
    9775 
    9776           CASE ( 's_SS' )
    9777              IF ( .NOT. ALLOCATED( s_ss_av ) )  THEN
    9778                 ALLOCATE( s_ss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    9779              ENDIF
    9780              s_ss_av = 0.0_wp
    9781 
    9782           CASE DEFAULT
    9783              CONTINUE
    9784 
    9785        END SELECT
     9856       IF ( variable(7:11) ==  'N_bin' )  THEN
     9857          IF ( .NOT. ALLOCATED( nbins_av ) )  THEN
     9858             ALLOCATE( nbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     9859          ENDIF
     9860          nbins_av = 0.0_wp
     9861
     9862       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
     9863          IF ( .NOT. ALLOCATED( mbins_av ) )  THEN
     9864             ALLOCATE( mbins_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) )
     9865          ENDIF
     9866          mbins_av = 0.0_wp
     9867
     9868       ELSE
     9869
     9870          SELECT CASE ( TRIM( variable(7:) ) )
     9871
     9872             CASE ( 'g_H2SO4' )
     9873                IF ( .NOT. ALLOCATED( g_h2so4_av ) )  THEN
     9874                   ALLOCATE( g_h2so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9875                ENDIF
     9876                g_h2so4_av = 0.0_wp
     9877
     9878             CASE ( 'g_HNO3' )
     9879                IF ( .NOT. ALLOCATED( g_hno3_av ) )  THEN
     9880                   ALLOCATE( g_hno3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9881                ENDIF
     9882                g_hno3_av = 0.0_wp
     9883
     9884             CASE ( 'g_NH3' )
     9885                IF ( .NOT. ALLOCATED( g_nh3_av ) )  THEN
     9886                   ALLOCATE( g_nh3_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9887                ENDIF
     9888                g_nh3_av = 0.0_wp
     9889
     9890             CASE ( 'g_OCNV' )
     9891                IF ( .NOT. ALLOCATED( g_ocnv_av ) )  THEN
     9892                   ALLOCATE( g_ocnv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9893                ENDIF
     9894                g_ocnv_av = 0.0_wp
     9895
     9896             CASE ( 'g_OCSV' )
     9897                IF ( .NOT. ALLOCATED( g_ocsv_av ) )  THEN
     9898                   ALLOCATE( g_ocsv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9899                ENDIF
     9900                g_ocsv_av = 0.0_wp
     9901
     9902             CASE ( 'LDSA' )
     9903                IF ( .NOT. ALLOCATED( ldsa_av ) )  THEN
     9904                   ALLOCATE( ldsa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9905                ENDIF
     9906                ldsa_av = 0.0_wp
     9907
     9908             CASE ( 'N_UFP' )
     9909                IF ( .NOT. ALLOCATED( nufp_av ) )  THEN
     9910                   ALLOCATE( nufp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9911                ENDIF
     9912                nufp_av = 0.0_wp
     9913
     9914             CASE ( 'Ntot' )
     9915                IF ( .NOT. ALLOCATED( ntot_av ) )  THEN
     9916                   ALLOCATE( ntot_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9917                ENDIF
     9918                ntot_av = 0.0_wp
     9919
     9920             CASE ( 'PM0.1' )
     9921                IF ( .NOT. ALLOCATED( pm01_av ) )  THEN
     9922                   ALLOCATE( pm01_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9923                ENDIF
     9924                pm01_av = 0.0_wp
     9925
     9926             CASE ( 'PM2.5' )
     9927                IF ( .NOT. ALLOCATED( pm25_av ) )  THEN
     9928                   ALLOCATE( pm25_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9929                ENDIF
     9930                pm25_av = 0.0_wp
     9931
     9932             CASE ( 'PM10' )
     9933                IF ( .NOT. ALLOCATED( pm10_av ) )  THEN
     9934                   ALLOCATE( pm10_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9935                ENDIF
     9936                pm10_av = 0.0_wp
     9937
     9938             CASE ( 's_BC' )
     9939                IF ( .NOT. ALLOCATED( s_bc_av ) )  THEN
     9940                   ALLOCATE( s_bc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9941                ENDIF
     9942                s_bc_av = 0.0_wp
     9943
     9944             CASE ( 's_DU' )
     9945                IF ( .NOT. ALLOCATED( s_du_av ) )  THEN
     9946                   ALLOCATE( s_du_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9947                ENDIF
     9948                s_du_av = 0.0_wp
     9949
     9950             CASE ( 's_H2O' )
     9951                IF ( .NOT. ALLOCATED( s_h2o_av ) )  THEN
     9952                   ALLOCATE( s_h2o_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9953                ENDIF
     9954                s_h2o_av = 0.0_wp
     9955
     9956             CASE ( 's_NH' )
     9957                IF ( .NOT. ALLOCATED( s_nh_av ) )  THEN
     9958                   ALLOCATE( s_nh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9959                ENDIF
     9960                s_nh_av = 0.0_wp
     9961
     9962             CASE ( 's_NO' )
     9963                IF ( .NOT. ALLOCATED( s_no_av ) )  THEN
     9964                   ALLOCATE( s_no_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9965                ENDIF
     9966                s_no_av = 0.0_wp
     9967
     9968             CASE ( 's_OC' )
     9969                IF ( .NOT. ALLOCATED( s_oc_av ) )  THEN
     9970                   ALLOCATE( s_oc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9971                ENDIF
     9972                s_oc_av = 0.0_wp
     9973
     9974             CASE ( 's_SO4' )
     9975                IF ( .NOT. ALLOCATED( s_so4_av ) )  THEN
     9976                   ALLOCATE( s_so4_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9977                ENDIF
     9978                s_so4_av = 0.0_wp
     9979
     9980             CASE ( 's_SS' )
     9981                IF ( .NOT. ALLOCATED( s_ss_av ) )  THEN
     9982                   ALLOCATE( s_ss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     9983                ENDIF
     9984                s_ss_av = 0.0_wp
     9985
     9986             CASE DEFAULT
     9987                CONTINUE
     9988
     9989          END SELECT
     9990
     9991       ENDIF
    97869992
    97879993    ELSEIF ( mode == 'sum' )  THEN
    97889994
    9789        SELECT CASE ( TRIM( variable ) )
    9790 
    9791           CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
    9792 
    9793              vari = TRIM( variable(3:) )
    9794 
    9795              SELECT CASE( vari )
    9796 
    9797                 CASE( 'H2SO4' )
    9798                    found_index = 1
    9799                    to_be_resorted => g_h2so4_av
    9800 
    9801                 CASE( 'HNO3' )
    9802                    found_index = 2
    9803                    to_be_resorted => g_hno3_av
    9804 
    9805                 CASE( 'NH3' )
    9806                    found_index = 3
    9807                    to_be_resorted => g_nh3_av
    9808 
    9809                 CASE( 'OCNV' )
    9810                    found_index = 4
    9811                    to_be_resorted => g_ocnv_av
    9812 
    9813                 CASE( 'OCSN' )
    9814                    found_index = 5
    9815                    to_be_resorted => g_ocsv_av
    9816 
    9817              END SELECT
    9818 
     9995       IF ( variable(7:11) ==  'N_bin' )  THEN
     9996          READ( variable(12:),* ) char_to_int
     9997          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     9998             ib = char_to_int
    98199999             DO  i = nxlg, nxrg
    982010000                DO  j = nysg, nyng
    982110001                   DO  k = nzb, nzt+1
    9822                       to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                              &
    9823                                               salsa_gas(found_index)%conc(k,j,i)
     10002                      nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i)
    982410003                   ENDDO
    982510004                ENDDO
    982610005             ENDDO
    9827 
    9828           CASE ( 'LDSA' )
     10006          ENDIF
     10007
     10008       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
     10009          READ( variable(12:),* ) char_to_int
     10010          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10011             ib = char_to_int
    982910012             DO  i = nxlg, nxrg
    983010013                DO  j = nysg, nyng
    983110014                   DO  k = nzb, nzt+1
    983210015                      temp_bin = 0.0_wp
    9833                       DO  ib = 1, nbins_aerosol
    9834 !
    9835 !--                      Diameter in micrometres
    9836                          mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
    9837 !
    9838 !--                      Deposition factor: alveolar (use ra_dry)
    9839                          df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +       &
    9840                                 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -     &
    9841                                 1.362_wp )**2 ) )
    9842 !
    9843 !--                      Number concentration in 1/cm3
    9844                          nc = 1.0E-6_wp * aerosol_number(ib)%conc(k,j,i)
    9845 !
    9846 !--                      Lung-deposited surface area LDSA (units mum2/cm3)
    9847                          temp_bin = temp_bin + pi * mean_d**2 * df * nc
     10016                      DO  ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol
     10017                         temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    984810018                      ENDDO
    9849                       ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin
     10019                      mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + temp_bin
    985010020                   ENDDO
    985110021                ENDDO
    985210022             ENDDO
    9853 
    9854           CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8',   &
    9855                  'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' )
    9856              DO  i = nxlg, nxrg
    9857                 DO  j = nysg, nyng
    9858                    DO  k = nzb, nzt+1
    9859                       DO  ib = 1, nbins_aerosol
    9860                          nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) + aerosol_number(ib)%conc(k,j,i)
     10023          ENDIF
     10024       ELSE
     10025
     10026          SELECT CASE ( TRIM( variable(7:) ) )
     10027
     10028             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
     10029
     10030                vari = TRIM( variable(9:) )  ! remove salsa_g_ from beginning
     10031
     10032                SELECT CASE( vari )
     10033
     10034                   CASE( 'H2SO4' )
     10035                      found_index = 1
     10036                      to_be_resorted => g_h2so4_av
     10037
     10038                   CASE( 'HNO3' )
     10039                      found_index = 2
     10040                      to_be_resorted => g_hno3_av
     10041
     10042                   CASE( 'NH3' )
     10043                      found_index = 3
     10044                      to_be_resorted => g_nh3_av
     10045
     10046                   CASE( 'OCNV' )
     10047                      found_index = 4
     10048                      to_be_resorted => g_ocnv_av
     10049
     10050                   CASE( 'OCSV' )
     10051                      found_index = 5
     10052                      to_be_resorted => g_ocsv_av
     10053
     10054                END SELECT
     10055
     10056                DO  i = nxlg, nxrg
     10057                   DO  j = nysg, nyng
     10058                      DO  k = nzb, nzt+1
     10059                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                           &
     10060                                                 salsa_gas(found_index)%conc(k,j,i)
    986110061                      ENDDO
    986210062                   ENDDO
    986310063                ENDDO
    9864              ENDDO
    9865 
    9866           CASE ( 'Ntot' )
    9867              DO  i = nxlg, nxrg
    9868                 DO  j = nysg, nyng
    9869                    DO  k = nzb, nzt+1
    9870                       DO  ib = 1, nbins_aerosol
    9871                          ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i)
     10064
     10065             CASE ( 'LDSA' )
     10066                DO  i = nxlg, nxrg
     10067                   DO  j = nysg, nyng
     10068                      DO  k = nzb, nzt+1
     10069                         temp_bin = 0.0_wp
     10070                         DO  ib = 1, nbins_aerosol
     10071   !
     10072   !--                      Diameter in micrometres
     10073                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
     10074   !
     10075   !--                      Deposition factor: alveolar (use ra_dry)
     10076                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
     10077                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
     10078                                   1.362_wp )**2 ) )
     10079   !
     10080   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
     10081                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
     10082                                       aerosol_number(ib)%conc(k,j,i)
     10083                         ENDDO
     10084                         ldsa_av(k,j,i) = ldsa_av(k,j,i) + temp_bin
    987210085                      ENDDO
    987310086                   ENDDO
    987410087                ENDDO
    9875              ENDDO
    9876 
    9877           CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8',   &
    9878                  'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12' )
    9879              DO  i = nxlg, nxrg
    9880                 DO  j = nysg, nyng
    9881                    DO  k = nzb, nzt+1
    9882                       DO  ib = 1, nbins_aerosol
    9883                          DO  ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol
    9884                             mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) + aerosol_mass(ic)%conc(k,j,i)
     10088
     10089             CASE ( 'N_UFP' )
     10090                DO  i = nxlg, nxrg
     10091                   DO  j = nysg, nyng
     10092                      DO  k = nzb, nzt+1
     10093                         temp_bin = 0.0_wp
     10094                         DO  ib = 1, nbins_aerosol
     10095                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     10096                               temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     10097                            ENDIF
     10098                         ENDDO
     10099                         nufp_av(k,j,i) = nufp_av(k,j,i) + temp_bin
     10100                      ENDDO
     10101                   ENDDO
     10102                ENDDO
     10103
     10104             CASE ( 'Ntot' )
     10105                DO  i = nxlg, nxrg
     10106                   DO  j = nysg, nyng
     10107                      DO  k = nzb, nzt+1
     10108                         DO  ib = 1, nbins_aerosol
     10109                            ntot_av(k,j,i) = ntot_av(k,j,i) + aerosol_number(ib)%conc(k,j,i)
    988510110                         ENDDO
    988610111                      ENDDO
    988710112                   ENDDO
    988810113                ENDDO
    9889              ENDDO
    9890 
    9891           CASE ( 'PM2.5' )
    9892              DO  i = nxlg, nxrg
    9893                 DO  j = nysg, nyng
    9894                    DO  k = nzb, nzt+1
    9895                       temp_bin = 0.0_wp
    9896                       DO  ib = 1, nbins_aerosol
    9897                          IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
    9898                             DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
    9899                                temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    9900                             ENDDO
    9901                          ENDIF
     10114
     10115             CASE ( 'PM0.1' )
     10116                DO  i = nxlg, nxrg
     10117                   DO  j = nysg, nyng
     10118                      DO  k = nzb, nzt+1
     10119                         temp_bin = 0.0_wp
     10120                         DO  ib = 1, nbins_aerosol
     10121                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     10122                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10123                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10124                               ENDDO
     10125                            ENDIF
     10126                         ENDDO
     10127                         pm01_av(k,j,i) = pm01_av(k,j,i) + temp_bin
    990210128                      ENDDO
    9903                       pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin
    990410129                   ENDDO
    990510130                ENDDO
    9906              ENDDO
    9907 
    9908           CASE ( 'PM10' )
    9909              DO  i = nxlg, nxrg
    9910                 DO  j = nysg, nyng
    9911                    DO  k = nzb, nzt+1
    9912                       temp_bin = 0.0_wp
    9913                       DO  ib = 1, nbins_aerosol
    9914                          IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
    9915                             DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
    9916                                temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    9917                             ENDDO
    9918                          ENDIF
     10131
     10132             CASE ( 'PM2.5' )
     10133                DO  i = nxlg, nxrg
     10134                   DO  j = nysg, nyng
     10135                      DO  k = nzb, nzt+1
     10136                         temp_bin = 0.0_wp
     10137                         DO  ib = 1, nbins_aerosol
     10138                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
     10139                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10140                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10141                               ENDDO
     10142                            ENDIF
     10143                         ENDDO
     10144                         pm25_av(k,j,i) = pm25_av(k,j,i) + temp_bin
    991910145                      ENDDO
    9920                       pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin
    992110146                   ENDDO
    992210147                ENDDO
    9923              ENDDO
    9924 
    9925           CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    9926              IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
    9927                 found_index = get_index( prtcl, TRIM( variable(3:) ) )
    9928                 IF ( TRIM( variable(3:) ) == 'BC' )   to_be_resorted => s_bc_av
    9929                 IF ( TRIM( variable(3:) ) == 'DU' )   to_be_resorted => s_du_av
    9930                 IF ( TRIM( variable(3:) ) == 'NH' )   to_be_resorted => s_nh_av
    9931                 IF ( TRIM( variable(3:) ) == 'NO' )   to_be_resorted => s_no_av
    9932                 IF ( TRIM( variable(3:) ) == 'OC' )   to_be_resorted => s_oc_av
    9933                 IF ( TRIM( variable(3:) ) == 'SO4' )  to_be_resorted => s_so4_av
    9934                 IF ( TRIM( variable(3:) ) == 'SS' )   to_be_resorted => s_ss_av
     10148
     10149             CASE ( 'PM10' )
     10150                DO  i = nxlg, nxrg
     10151                   DO  j = nysg, nyng
     10152                      DO  k = nzb, nzt+1
     10153                         temp_bin = 0.0_wp
     10154                         DO  ib = 1, nbins_aerosol
     10155                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
     10156                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10157                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10158                               ENDDO
     10159                            ENDIF
     10160                         ENDDO
     10161                         pm10_av(k,j,i) = pm10_av(k,j,i) + temp_bin
     10162                      ENDDO
     10163                   ENDDO
     10164                ENDDO
     10165
     10166             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10167                IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
     10168                   found_index = get_index( prtcl, TRIM( variable(9:) ) )
     10169                   IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
     10170                   IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
     10171                   IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
     10172                   IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
     10173                   IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
     10174                   IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
     10175                   IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
     10176                   DO  i = nxlg, nxrg
     10177                      DO  j = nysg, nyng
     10178                         DO  k = nzb, nzt+1
     10179                            DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
     10180                               to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                     &
     10181                                                       aerosol_mass(ic)%conc(k,j,i)
     10182                            ENDDO
     10183                         ENDDO
     10184                      ENDDO
     10185                   ENDDO
     10186                ENDIF
     10187
     10188             CASE ( 's_H2O' )
     10189                found_index = get_index( prtcl,'H2O' )
     10190                to_be_resorted => s_h2o_av
    993510191                DO  i = nxlg, nxrg
    993610192                   DO  j = nysg, nyng
    993710193                      DO  k = nzb, nzt+1
    993810194                         DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
    9939                             to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                        &
    9940                                                     aerosol_mass(ic)%conc(k,j,i)
     10195                            s_h2o_av(k,j,i) = s_h2o_av(k,j,i) + aerosol_mass(ic)%conc(k,j,i)
    994110196                         ENDDO
    994210197                      ENDDO
    994310198                   ENDDO
    994410199                ENDDO
    9945              ENDIF
    9946 
    9947           CASE ( 's_H2O' )
    9948              found_index = get_index( prtcl,'H2O' )
    9949              to_be_resorted => s_h2o_av
     10200
     10201             CASE DEFAULT
     10202                CONTINUE
     10203
     10204          END SELECT
     10205
     10206       ENDIF
     10207
     10208    ELSEIF ( mode == 'average' )  THEN
     10209
     10210       IF ( variable(7:11) ==  'N_bin' )  THEN
     10211          READ( variable(12:),* ) char_to_int
     10212          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10213             ib = char_to_int
    995010214             DO  i = nxlg, nxrg
    995110215                DO  j = nysg, nyng
    995210216                   DO  k = nzb, nzt+1
    9953                       DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
    9954                          to_be_resorted(k,j,i) = to_be_resorted(k,j,i) +                           &
    9955                                                  aerosol_mass(ic)%conc(k,j,i)
    9956                       ENDDO
     10217                      nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp )
    995710218                   ENDDO
    995810219                ENDDO
    995910220             ENDDO
    9960 
    9961           CASE DEFAULT
    9962              CONTINUE
    9963 
    9964        END SELECT
    9965 
    9966     ELSEIF ( mode == 'average' )  THEN
    9967 
    9968        SELECT CASE ( TRIM( variable ) )
    9969 
    9970           CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
    9971              IF ( TRIM( variable(3:) ) == 'H2SO4' )  THEN
    9972                 found_index = 1
    9973                 to_be_resorted => g_h2so4_av
    9974              ELSEIF ( TRIM( variable(3:) ) == 'HNO3' )  THEN
    9975                 found_index = 2
    9976                 to_be_resorted => g_hno3_av
    9977              ELSEIF ( TRIM( variable(3:) ) == 'NH3' )  THEN
    9978                 found_index = 3
    9979                 to_be_resorted => g_nh3_av
    9980              ELSEIF ( TRIM( variable(3:) ) == 'OCNV' )  THEN
    9981                 found_index = 4
    9982                 to_be_resorted => g_ocnv_av
    9983              ELSEIF ( TRIM( variable(3:) ) == 'OCSV' )  THEN
    9984                 found_index = 5
    9985                 to_be_resorted => g_ocsv_av
    9986              ENDIF
     10221          ENDIF
     10222
     10223       ELSEIF ( variable(7:11) ==  'm_bin' )  THEN
     10224          READ( variable(12:),* ) char_to_int
     10225          IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10226             ib = char_to_int
    998710227             DO  i = nxlg, nxrg
    998810228                DO  j = nysg, nyng
    998910229                   DO  k = nzb, nzt+1
    9990                       to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                              &
    9991                                               REAL( average_count_3d, KIND=wp )
     10230                      mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp)
    999210231                   ENDDO
    999310232                ENDDO
    999410233             ENDDO
    9995 
    9996           CASE ( 'LDSA' )
    9997              DO  i = nxlg, nxrg
    9998                 DO  j = nysg, nyng
    9999                    DO  k = nzb, nzt+1
    10000                       ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10001                    ENDDO
    10002                 ENDDO
    10003              ENDDO
    10004 
    10005           CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8',   &
    10006                  'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' )
    10007              DO  i = nxlg, nxrg
    10008                 DO  j = nysg, nyng
    10009                    DO  k = nzb, nzt+1
    10010                       DO  ib = 1, nbins_aerosol
    10011                          nbins_av(k,j,i,ib) = nbins_av(k,j,i,ib) / REAL( average_count_3d, KIND=wp )
    10012                       ENDDO
    10013                    ENDDO
    10014                 ENDDO
    10015              ENDDO
    10016 
    10017           CASE ( 'Ntot' )
    10018              DO  i = nxlg, nxrg
    10019                 DO  j = nysg, nyng
    10020                    DO  k = nzb, nzt+1
    10021                       ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10022                    ENDDO
    10023                 ENDDO
    10024              ENDDO
    10025 
    10026           CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8',   &
    10027                  'm_bin9', 'm_bin10', 'm_bin11', 'm_bin12' )
    10028              DO  i = nxlg, nxrg
    10029                 DO  j = nysg, nyng
    10030                    DO  k = nzb, nzt+1
    10031                       DO  ib = 1, nbins_aerosol
    10032                          DO  ic = ib, nbins_aerosol * ncomponents_mass, nbins_aerosol
    10033                             mbins_av(k,j,i,ib) = mbins_av(k,j,i,ib) /                              &
    10034                                                  REAL( average_count_3d, KIND=wp)
    10035                          ENDDO
    10036                       ENDDO
    10037                    ENDDO
    10038                 ENDDO
    10039              ENDDO
    10040 
    10041           CASE ( 'PM2.5' )
    10042              DO  i = nxlg, nxrg
    10043                 DO  j = nysg, nyng
    10044                    DO  k = nzb, nzt+1
    10045                       pm25_av(k,j,i) = pm25_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10046                    ENDDO
    10047                 ENDDO
    10048              ENDDO
    10049 
    10050           CASE ( 'PM10' )
    10051              DO  i = nxlg, nxrg
    10052                 DO  j = nysg, nyng
    10053                    DO  k = nzb, nzt+1
    10054                       pm10_av(k,j,i) = pm10_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    10055                    ENDDO
    10056                 ENDDO
    10057              ENDDO
    10058 
    10059           CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    10060              IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
    10061                 IF ( TRIM( variable(3:) ) == 'BC' )   to_be_resorted => s_bc_av
    10062                 IF ( TRIM( variable(3:) ) == 'DU' )   to_be_resorted => s_du_av
    10063                 IF ( TRIM( variable(3:) ) == 'NH' )   to_be_resorted => s_nh_av
    10064                 IF ( TRIM( variable(3:) ) == 'NO' )   to_be_resorted => s_no_av
    10065                 IF ( TRIM( variable(3:) ) == 'OC' )   to_be_resorted => s_oc_av
    10066                 IF ( TRIM( variable(3:) ) == 'SO4' )  to_be_resorted => s_so4_av
    10067                 IF ( TRIM( variable(3:) ) == 'SS' )   to_be_resorted => s_ss_av
     10234          ENDIF
     10235       ELSE
     10236
     10237          SELECT CASE ( TRIM( variable(7:) ) )
     10238
     10239             CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
     10240                IF ( TRIM( variable(9:) ) == 'H2SO4' )  THEN  ! 9: remove salsa_g_ from beginning
     10241                   found_index = 1
     10242                   to_be_resorted => g_h2so4_av
     10243                ELSEIF ( TRIM( variable(9:) ) == 'HNO3' )  THEN
     10244                   found_index = 2
     10245                   to_be_resorted => g_hno3_av
     10246                ELSEIF ( TRIM( variable(9:) ) == 'NH3' )  THEN
     10247                   found_index = 3
     10248                   to_be_resorted => g_nh3_av
     10249                ELSEIF ( TRIM( variable(9:) ) == 'OCNV' )  THEN
     10250                   found_index = 4
     10251                   to_be_resorted => g_ocnv_av
     10252                ELSEIF ( TRIM( variable(9:) ) == 'OCSV' )  THEN
     10253                   found_index = 5
     10254                   to_be_resorted => g_ocsv_av
     10255                ENDIF
    1006810256                DO  i = nxlg, nxrg
    1006910257                   DO  j = nysg, nyng
     
    1007410262                   ENDDO
    1007510263                ENDDO
    10076              ENDIF
    10077 
    10078           CASE ( 's_H2O' )
    10079              to_be_resorted => s_h2o_av
    10080              DO  i = nxlg, nxrg
    10081                 DO  j = nysg, nyng
    10082                    DO  k = nzb, nzt+1
    10083                       to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                              &
    10084                                               REAL( average_count_3d, KIND=wp )
     10264
     10265             CASE ( 'LDSA' )
     10266                DO  i = nxlg, nxrg
     10267                   DO  j = nysg, nyng
     10268                      DO  k = nzb, nzt+1
     10269                         ldsa_av(k,j,i) = ldsa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10270                      ENDDO
    1008510271                   ENDDO
    1008610272                ENDDO
    10087              ENDDO
    10088 
    10089        END SELECT
    10090 
     10273
     10274             CASE ( 'N_UFP' )
     10275                DO  i = nxlg, nxrg
     10276                   DO  j = nysg, nyng
     10277                      DO  k = nzb, nzt+1
     10278                         nufp_av(k,j,i) = nufp_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10279                      ENDDO
     10280                   ENDDO
     10281                ENDDO
     10282
     10283             CASE ( 'Ntot' )
     10284                DO  i = nxlg, nxrg
     10285                   DO  j = nysg, nyng
     10286                      DO  k = nzb, nzt+1
     10287                         ntot_av(k,j,i) = ntot_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10288                      ENDDO
     10289                   ENDDO
     10290                ENDDO
     10291
     10292
     10293             CASE ( 'PM0.1' )
     10294                DO  i = nxlg, nxrg
     10295                   DO  j = nysg, nyng
     10296                      DO  k = nzb, nzt+1
     10297                         pm01_av(k,j,i) = pm01_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10298                      ENDDO
     10299                   ENDDO
     10300                ENDDO
     10301
     10302             CASE ( 'PM2.5' )
     10303                DO  i = nxlg, nxrg
     10304                   DO  j = nysg, nyng
     10305                      DO  k = nzb, nzt+1
     10306                         pm25_av(k,j,i) = pm25_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10307                      ENDDO
     10308                   ENDDO
     10309                ENDDO
     10310
     10311             CASE ( 'PM10' )
     10312                DO  i = nxlg, nxrg
     10313                   DO  j = nysg, nyng
     10314                      DO  k = nzb, nzt+1
     10315                         pm10_av(k,j,i) = pm10_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     10316                      ENDDO
     10317                   ENDDO
     10318                ENDDO
     10319
     10320             CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10321                IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
     10322                   IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
     10323                   IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
     10324                   IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
     10325                   IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
     10326                   IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
     10327                   IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
     10328                   IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
     10329                   DO  i = nxlg, nxrg
     10330                      DO  j = nysg, nyng
     10331                         DO  k = nzb, nzt+1
     10332                            to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                        &
     10333                                                    REAL( average_count_3d, KIND=wp )
     10334                         ENDDO
     10335                      ENDDO
     10336                   ENDDO
     10337                ENDIF
     10338
     10339             CASE ( 's_H2O' )
     10340                to_be_resorted => s_h2o_av
     10341                DO  i = nxlg, nxrg
     10342                   DO  j = nysg, nyng
     10343                      DO  k = nzb, nzt+1
     10344                         to_be_resorted(k,j,i) = to_be_resorted(k,j,i) /                           &
     10345                                                 REAL( average_count_3d, KIND=wp )
     10346                      ENDDO
     10347                   ENDDO
     10348                ENDDO
     10349
     10350          END SELECT
     10351
     10352       ENDIF
    1009110353    ENDIF
    1009210354
     
    1011510377
    1011610378    INTEGER(iwp) ::  av           !<
     10379    INTEGER(iwp) ::  char_to_int  !< for converting character to integer
    1011710380    INTEGER(iwp) ::  found_index  !< index of a chemical compound
    1011810381    INTEGER(iwp) ::  i            !<
     
    1013210395    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
    1013310396    REAL(wp) ::  mean_d                   !< Particle diameter in micrometres
    10134     REAL(wp) ::  nc                       !< Particle number concentration in units 1/cm**3
    1013510397    REAL(wp) ::  temp_bin                 !< temporary array for calculating output variables
    1013610398
     
    1014510407    temp_bin  = 0.0_wp
    1014610408
    10147     SELECT CASE ( TRIM( variable( 1:LEN( TRIM( variable ) ) - 3 ) ) )  ! cut out _xy, _xz or _yz
    10148 
    10149        CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
    10150           vari = TRIM( variable( 3:LEN( TRIM( variable ) ) - 3 ) )
    10151           IF ( av == 0 )  THEN
    10152              IF ( vari == 'H2SO4')  found_index = 1
    10153              IF ( vari == 'HNO3')   found_index = 2
    10154              IF ( vari == 'NH3')    found_index = 3
    10155              IF ( vari == 'OCNV')   found_index = 4
    10156              IF ( vari == 'OCSV')   found_index = 5
    10157              DO  i = nxl, nxr
    10158                 DO  j = nys, nyn
    10159                    DO  k = nzb_do, nzt_do
    10160                       local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i), REAL( fill_value, &
    10161                                                KIND = wp ),  BTEST( wall_flags_0(k,j,i), 0 ) )
    10162                    ENDDO
    10163                 ENDDO
    10164              ENDDO
    10165           ELSE
    10166              IF ( vari == 'H2SO4' )  to_be_resorted => g_h2so4_av
    10167              IF ( vari == 'HNO3' )   to_be_resorted => g_hno3_av
    10168              IF ( vari == 'NH3' )    to_be_resorted => g_nh3_av
    10169              IF ( vari == 'OCNV' )   to_be_resorted => g_ocnv_av
    10170              IF ( vari == 'OCSV' )   to_be_resorted => g_ocsv_av
    10171              DO  i = nxl, nxr
    10172                 DO  j = nys, nyn
    10173                    DO  k = nzb_do, nzt_do
    10174                       local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,            &
    10175                                                KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
    10176                    ENDDO
    10177                 ENDDO
    10178              ENDDO
    10179           ENDIF
    10180 
    10181           IF ( mode == 'xy' )  grid = 'zu'
    10182 
    10183        CASE ( 'LDSA' )
    10184           IF ( av == 0 )  THEN
    10185              DO  i = nxl, nxr
    10186                 DO  j = nys, nyn
    10187                    DO  k = nzb_do, nzt_do
    10188                       temp_bin = 0.0_wp
    10189                       DO  ib = 1, nbins_aerosol
    10190 !
    10191 !--                      Diameter in micrometres
    10192                          mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
    10193 !
    10194 !--                      Deposition factor: alveolar
    10195                          df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +       &
    10196                                 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -     &
    10197                                 1.362_wp )**2 ) )
    10198 !
    10199 !--                      Number concentration in 1/cm3
    10200                          nc = 1.0E-6_wp * aerosol_number(ib)%conc(k,j,i)
    10201 !
    10202 !--                      Lung-deposited surface area LDSA (units mum2/cm3)
    10203                          temp_bin = temp_bin + pi * mean_d**2 * df * nc
    10204                       ENDDO
    10205 
    10206                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10207                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10208                    ENDDO
    10209                 ENDDO
    10210              ENDDO
    10211           ELSE
    10212              DO  i = nxl, nxr
    10213                 DO  j = nys, nyn
    10214                    DO  k = nzb_do, nzt_do
    10215                       local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10216                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10217                    ENDDO
    10218                 ENDDO
    10219              ENDDO
    10220           ENDIF
    10221 
    10222           IF ( mode == 'xy' )  grid = 'zu'
    10223 
    10224        CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4',   'N_bin5',  'N_bin6', 'N_bin7', 'N_bin8',   &
    10225               'N_bin9', 'N_bin10' , 'N_bin11', 'N_bin12' )
    10226           vari = TRIM( variable( 6:LEN( TRIM( variable ) ) - 3 ) )
    10227 
    10228           IF ( vari == '1' ) ib = 1
    10229           IF ( vari == '2' ) ib = 2
    10230           IF ( vari == '3' ) ib = 3
    10231           IF ( vari == '4' ) ib = 4
    10232           IF ( vari == '5' ) ib = 5
    10233           IF ( vari == '6' ) ib = 6
    10234           IF ( vari == '7' ) ib = 7
    10235           IF ( vari == '8' ) ib = 8
    10236           IF ( vari == '9' ) ib = 9
    10237           IF ( vari == '10' ) ib = 10
    10238           IF ( vari == '11' ) ib = 11
    10239           IF ( vari == '12' ) ib = 12
    10240 
     10409    IF ( variable(7:11)  == 'N_bin' )  THEN
     10410
     10411       READ( variable( 12:LEN( TRIM( variable ) ) - 3 ), * ) char_to_int
     10412       IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10413
     10414          ib = char_to_int
    1024110415          IF ( av == 0 )  THEN
    1024210416             DO  i = nxl, nxr
     
    1025810432             ENDDO
    1025910433          ENDIF
    10260 
    1026110434          IF ( mode == 'xy' )  grid = 'zu'
    10262 
    10263        CASE ( 'Ntot' )
    10264 
    10265           IF ( av == 0 )  THEN
    10266              DO  i = nxl, nxr
    10267                 DO  j = nys, nyn
    10268                    DO  k = nzb_do, nzt_do
    10269                       temp_bin = 0.0_wp
    10270                       DO  ib = 1, nbins_aerosol
    10271                          temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
    10272                       ENDDO
    10273                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10274                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10275                    ENDDO
    10276                 ENDDO
    10277              ENDDO
    10278           ELSE
    10279              DO  i = nxl, nxr
    10280                 DO  j = nys, nyn
    10281                    DO  k = nzb_do, nzt_do
    10282                       local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10283                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10284                    ENDDO
    10285                 ENDDO
    10286              ENDDO
    10287           ENDIF
    10288 
    10289           IF ( mode == 'xy' )  grid = 'zu'
    10290 
    10291        CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4',   'm_bin5',  'm_bin6', 'm_bin7', 'm_bin8',   &
    10292               'm_bin9', 'm_bin10' , 'm_bin11', 'm_bin12' )
    10293           vari = TRIM( variable( 6:LEN( TRIM( variable ) ) - 3 ) )
    10294 
    10295           IF ( vari == '1' ) ib = 1
    10296           IF ( vari == '2' ) ib = 2
    10297           IF ( vari == '3' ) ib = 3
    10298           IF ( vari == '4' ) ib = 4
    10299           IF ( vari == '5' ) ib = 5
    10300           IF ( vari == '6' ) ib = 6
    10301           IF ( vari == '7' ) ib = 7
    10302           IF ( vari == '8' ) ib = 8
    10303           IF ( vari == '9' ) ib = 9
    10304           IF ( vari == '10' ) ib = 10
    10305           IF ( vari == '11' ) ib = 11
    10306           IF ( vari == '12' ) ib = 12
    10307 
     10435       ENDIF
     10436
     10437    ELSEIF ( variable(7:11)  == 'm_bin' )  THEN
     10438
     10439       READ( variable( 12:LEN( TRIM( variable ) ) - 3 ), * ) char_to_int
     10440       IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10441
     10442          ib = char_to_int
    1030810443          IF ( av == 0 )  THEN
    1030910444             DO  i = nxl, nxr
     
    1032910464             ENDDO
    1033010465          ENDIF
    10331 
    1033210466          IF ( mode == 'xy' )  grid = 'zu'
    10333 
    10334        CASE ( 'PM2.5' )
    10335           IF ( av == 0 )  THEN
    10336              DO  i = nxl, nxr
    10337                 DO  j = nys, nyn
    10338                    DO  k = nzb_do, nzt_do
    10339                       temp_bin = 0.0_wp
    10340                       DO  ib = 1, nbins_aerosol
    10341                          IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
    10342                             DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10467       ENDIF
     10468
     10469    ELSE
     10470
     10471       SELECT CASE ( TRIM( variable( 7:LEN( TRIM( variable ) ) - 3 ) ) )  ! cut out _xy, _xz or _yz
     10472
     10473          CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV', 'g_OCSV' )
     10474             vari = TRIM( variable( 9:LEN( TRIM( variable ) ) - 3 ) )  ! 9: remove salsa_g_
     10475             IF ( av == 0 )  THEN
     10476                IF ( vari == 'H2SO4')  found_index = 1
     10477                IF ( vari == 'HNO3')   found_index = 2
     10478                IF ( vari == 'NH3')    found_index = 3
     10479                IF ( vari == 'OCNV')   found_index = 4
     10480                IF ( vari == 'OCSV')   found_index = 5
     10481                DO  i = nxl, nxr
     10482                   DO  j = nys, nyn
     10483                      DO  k = nzb_do, nzt_do
     10484                         local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i),              &
     10485                                                  REAL( fill_value,  KIND = wp ),                  &
     10486                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10487                      ENDDO
     10488                   ENDDO
     10489                ENDDO
     10490             ELSE
     10491                IF ( vari == 'H2SO4' )  to_be_resorted => g_h2so4_av
     10492                IF ( vari == 'HNO3' )   to_be_resorted => g_hno3_av
     10493                IF ( vari == 'NH3' )    to_be_resorted => g_nh3_av
     10494                IF ( vari == 'OCNV' )   to_be_resorted => g_ocnv_av
     10495                IF ( vari == 'OCSV' )   to_be_resorted => g_ocsv_av
     10496                DO  i = nxl, nxr
     10497                   DO  j = nys, nyn
     10498                      DO  k = nzb_do, nzt_do
     10499                         local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
     10500                                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
     10501                      ENDDO
     10502                   ENDDO
     10503                ENDDO
     10504             ENDIF
     10505
     10506             IF ( mode == 'xy' )  grid = 'zu'
     10507
     10508          CASE ( 'LDSA' )
     10509             IF ( av == 0 )  THEN
     10510                DO  i = nxl, nxr
     10511                   DO  j = nys, nyn
     10512                      DO  k = nzb_do, nzt_do
     10513                         temp_bin = 0.0_wp
     10514                         DO  ib = 1, nbins_aerosol
     10515   !
     10516   !--                      Diameter in micrometres
     10517                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
     10518   !
     10519   !--                      Deposition factor: alveolar
     10520                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
     10521                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
     10522                                   1.362_wp )**2 ) )
     10523   !
     10524   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
     10525                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
     10526                                       aerosol_number(ib)%conc(k,j,i)
     10527                         ENDDO
     10528
     10529                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10530                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10531                      ENDDO
     10532                   ENDDO
     10533                ENDDO
     10534             ELSE
     10535                DO  i = nxl, nxr
     10536                   DO  j = nys, nyn
     10537                      DO  k = nzb_do, nzt_do
     10538                         local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10539                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10540                      ENDDO
     10541                   ENDDO
     10542                ENDDO
     10543             ENDIF
     10544
     10545             IF ( mode == 'xy' )  grid = 'zu'
     10546
     10547          CASE ( 'N_UFP' )
     10548
     10549             IF ( av == 0 )  THEN
     10550                DO  i = nxl, nxr
     10551                   DO  j = nys, nyn
     10552                      DO  k = nzb_do, nzt_do
     10553                         temp_bin = 0.0_wp
     10554                         DO  ib = 1, nbins_aerosol
     10555                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     10556                               temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     10557                            ENDIF
     10558                         ENDDO
     10559                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10560                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10561                      ENDDO
     10562                   ENDDO
     10563                ENDDO
     10564             ELSE
     10565                DO  i = nxl, nxr
     10566                   DO  j = nys, nyn
     10567                      DO  k = nzb_do, nzt_do
     10568                         local_pf(i,j,k) = MERGE( nufp_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10569                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10570                      ENDDO
     10571                   ENDDO
     10572                ENDDO
     10573             ENDIF
     10574
     10575             IF ( mode == 'xy' )  grid = 'zu'
     10576
     10577          CASE ( 'Ntot' )
     10578
     10579             IF ( av == 0 )  THEN
     10580                DO  i = nxl, nxr
     10581                   DO  j = nys, nyn
     10582                      DO  k = nzb_do, nzt_do
     10583                         temp_bin = 0.0_wp
     10584                         DO  ib = 1, nbins_aerosol
     10585                            temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     10586                         ENDDO
     10587                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10588                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10589                      ENDDO
     10590                   ENDDO
     10591                ENDDO
     10592             ELSE
     10593                DO  i = nxl, nxr
     10594                   DO  j = nys, nyn
     10595                      DO  k = nzb_do, nzt_do
     10596                         local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10597                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10598                      ENDDO
     10599                   ENDDO
     10600                ENDDO
     10601             ENDIF
     10602
     10603             IF ( mode == 'xy' )  grid = 'zu'
     10604
     10605          CASE ( 'PM0.1' )
     10606             IF ( av == 0 )  THEN
     10607                DO  i = nxl, nxr
     10608                   DO  j = nys, nyn
     10609                      DO  k = nzb_do, nzt_do
     10610                         temp_bin = 0.0_wp
     10611                         DO  ib = 1, nbins_aerosol
     10612                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     10613                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10614                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10615                               ENDDO
     10616                            ENDIF
     10617                         ENDDO
     10618                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10619                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10620                      ENDDO
     10621                   ENDDO
     10622                ENDDO
     10623             ELSE
     10624                DO  i = nxl, nxr
     10625                   DO  j = nys, nyn
     10626                      DO  k = nzb_do, nzt_do
     10627                         local_pf(i,j,k) = MERGE( pm01_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10628                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10629                      ENDDO
     10630                   ENDDO
     10631                ENDDO
     10632             ENDIF
     10633
     10634             IF ( mode == 'xy' )  grid = 'zu'
     10635
     10636          CASE ( 'PM2.5' )
     10637             IF ( av == 0 )  THEN
     10638                DO  i = nxl, nxr
     10639                   DO  j = nys, nyn
     10640                      DO  k = nzb_do, nzt_do
     10641                         temp_bin = 0.0_wp
     10642                         DO  ib = 1, nbins_aerosol
     10643                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
     10644                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10645                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10646                               ENDDO
     10647                            ENDIF
     10648                         ENDDO
     10649                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10650                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10651                      ENDDO
     10652                   ENDDO
     10653                ENDDO
     10654             ELSE
     10655                DO  i = nxl, nxr
     10656                   DO  j = nys, nyn
     10657                      DO  k = nzb_do, nzt_do
     10658                         local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10659                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10660                      ENDDO
     10661                   ENDDO
     10662                ENDDO
     10663             ENDIF
     10664
     10665             IF ( mode == 'xy' )  grid = 'zu'
     10666
     10667          CASE ( 'PM10' )
     10668             IF ( av == 0 )  THEN
     10669                DO  i = nxl, nxr
     10670                   DO  j = nys, nyn
     10671                      DO  k = nzb_do, nzt_do
     10672                         temp_bin = 0.0_wp
     10673                         DO  ib = 1, nbins_aerosol
     10674                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
     10675                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10676                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     10677                               ENDDO
     10678                            ENDIF
     10679                         ENDDO
     10680                         local_pf(i,j,k) = MERGE( temp_bin,  REAL( fill_value, KIND = wp ),        &
     10681                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10682                      ENDDO
     10683                   ENDDO
     10684                ENDDO
     10685             ELSE
     10686                DO  i = nxl, nxr
     10687                   DO  j = nys, nyn
     10688                      DO  k = nzb_do, nzt_do
     10689                         local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10690                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10691                      ENDDO
     10692                   ENDDO
     10693                ENDDO
     10694             ENDIF
     10695
     10696             IF ( mode == 'xy' )  grid = 'zu'
     10697
     10698          CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     10699             vari = TRIM( variable( 9:LEN( TRIM( variable ) ) - 3 ) )  ! 9: remove salsa_s_
     10700             IF ( is_used( prtcl, vari ) )  THEN
     10701                found_index = get_index( prtcl, vari )
     10702                IF ( av == 0 )  THEN
     10703                   DO  i = nxl, nxr
     10704                      DO  j = nys, nyn
     10705                         DO  k = nzb_do, nzt_do
     10706                            temp_bin = 0.0_wp
     10707                            DO  ic = ( found_index-1 ) * nbins_aerosol+1, found_index * nbins_aerosol
    1034310708                               temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    1034410709                            ENDDO
    10345                          ENDIF
     10710                            local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),      &
     10711                                                     BTEST( wall_flags_0(k,j,i), 0 ) )
     10712                         ENDDO
    1034610713                      ENDDO
    10347                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10348                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    1034910714                   ENDDO
    10350                 ENDDO
    10351              ENDDO
    10352           ELSE
    10353              DO  i = nxl, nxr
    10354                 DO  j = nys, nyn
    10355                    DO  k = nzb_do, nzt_do
    10356                       local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10357                                                BTEST( wall_flags_0(k,j,i), 0 ) )
     10715                ELSE
     10716                   IF ( vari == 'BC' )   to_be_resorted => s_bc_av
     10717                   IF ( vari == 'DU' )   to_be_resorted => s_du_av
     10718                   IF ( vari == 'NH' )   to_be_resorted => s_nh_av
     10719                   IF ( vari == 'NO' )   to_be_resorted => s_no_av
     10720                   IF ( vari == 'OC' )   to_be_resorted => s_oc_av
     10721                   IF ( vari == 'SO4' )  to_be_resorted => s_so4_av
     10722                   IF ( vari == 'SS' )   to_be_resorted => s_ss_av
     10723                   DO  i = nxl, nxr
     10724                      DO  j = nys, nyn
     10725                         DO  k = nzb_do, nzt_do
     10726                            local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,      &
     10727                                                     KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
     10728                         ENDDO
     10729                      ENDDO
    1035810730                   ENDDO
    10359                 ENDDO
    10360              ENDDO
    10361           ENDIF
    10362 
    10363           IF ( mode == 'xy' )  grid = 'zu'
    10364 
    10365        CASE ( 'PM10' )
    10366           IF ( av == 0 )  THEN
    10367              DO  i = nxl, nxr
    10368                 DO  j = nys, nyn
    10369                    DO  k = nzb_do, nzt_do
    10370                       temp_bin = 0.0_wp
    10371                       DO  ib = 1, nbins_aerosol
    10372                          IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
    10373                             DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
    10374                                temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10375                             ENDDO
    10376                          ENDIF
    10377                       ENDDO
    10378                       local_pf(i,j,k) = MERGE( temp_bin,  REAL( fill_value, KIND = wp ),           &
    10379                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10380                    ENDDO
    10381                 ENDDO
    10382              ENDDO
    10383           ELSE
    10384              DO  i = nxl, nxr
    10385                 DO  j = nys, nyn
    10386                    DO  k = nzb_do, nzt_do
    10387                       local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10388                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10389                    ENDDO
    10390                 ENDDO
    10391              ENDDO
    10392           ENDIF
    10393 
    10394           IF ( mode == 'xy' )  grid = 'zu'
    10395 
    10396        CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    10397           vari = TRIM( variable( 3:LEN( TRIM( variable ) ) - 3 ) )
    10398           IF ( is_used( prtcl, vari ) )  THEN
    10399              found_index = get_index( prtcl, vari )
     10731                ENDIF
     10732             ELSE
     10733                local_pf = fill_value
     10734             ENDIF
     10735
     10736             IF ( mode == 'xy' )  grid = 'zu'
     10737
     10738          CASE ( 's_H2O' )
     10739             found_index = get_index( prtcl, 'H2O' )
    1040010740             IF ( av == 0 )  THEN
    1040110741                DO  i = nxl, nxr
     
    1041210752                ENDDO
    1041310753             ELSE
    10414                 IF ( vari == 'BC' )   to_be_resorted => s_bc_av
    10415                 IF ( vari == 'DU' )   to_be_resorted => s_du_av
    10416                 IF ( vari == 'NH' )   to_be_resorted => s_nh_av
    10417                 IF ( vari == 'NO' )   to_be_resorted => s_no_av
    10418                 IF ( vari == 'OC' )   to_be_resorted => s_oc_av
    10419                 IF ( vari == 'SO4' )  to_be_resorted => s_so4_av
    10420                 IF ( vari == 'SS' )   to_be_resorted => s_ss_av
     10754                to_be_resorted => s_h2o_av
    1042110755                DO  i = nxl, nxr
    1042210756                   DO  j = nys, nyn
     
    1042810762                ENDDO
    1042910763             ENDIF
    10430           ELSE
    10431              local_pf = fill_value
    10432           ENDIF
    10433 
    10434           IF ( mode == 'xy' )  grid = 'zu'
    10435 
    10436        CASE ( 's_H2O' )
    10437           found_index = get_index( prtcl, 'H2O' )
    10438           IF ( av == 0 )  THEN
    10439              DO  i = nxl, nxr
    10440                 DO  j = nys, nyn
    10441                    DO  k = nzb_do, nzt_do
    10442                       temp_bin = 0.0_wp
    10443                       DO  ic = ( found_index-1 ) * nbins_aerosol+1, found_index * nbins_aerosol
    10444                          temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10445                       ENDDO
    10446                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
    10447                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10448                    ENDDO
    10449                 ENDDO
    10450              ENDDO
    10451           ELSE
    10452              to_be_resorted => s_h2o_av
    10453              DO  i = nxl, nxr
    10454                 DO  j = nys, nyn
    10455                    DO  k = nzb_do, nzt_do
    10456                       local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
    10457                                                KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
    10458                    ENDDO
    10459                 ENDDO
    10460              ENDDO
    10461           ENDIF
    10462 
    10463           IF ( mode == 'xy' )  grid = 'zu'
    10464 
    10465        CASE DEFAULT
    10466           found = .FALSE.
    10467           grid  = 'none'
    10468 
    10469     END SELECT
     10764
     10765             IF ( mode == 'xy' )  grid = 'zu'
     10766
     10767          CASE DEFAULT
     10768             found = .FALSE.
     10769             grid  = 'none'
     10770
     10771       END SELECT
     10772
     10773    ENDIF
    1047010774
    1047110775 END SUBROUTINE salsa_data_output_2d
     
    1048910793
    1049010794    INTEGER(iwp) ::  av           !<
     10795    INTEGER(iwp) ::  char_to_int  !< for converting character to integer
    1049110796    INTEGER(iwp) ::  found_index  !< index of a chemical compound
    1049210797    INTEGER(iwp) ::  ib           !< running index: size bins
     
    1050510810    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
    1050610811    REAL(wp) ::  mean_d                   !< Particle diameter in micrometres
    10507     REAL(wp) ::  nc                       !< Particle number concentration in units 1/cm**3
    1050810812    REAL(wp) ::  temp_bin                 !< temporary array for calculating output variables
    1050910813
     
    1051510819    temp_bin  = 0.0_wp
    1051610820
    10517     SELECT CASE ( TRIM( variable ) )
    10518 
    10519        CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
    10520           IF ( av == 0 )  THEN
    10521              IF ( TRIM( variable ) == 'g_H2SO4')  found_index = 1
    10522              IF ( TRIM( variable ) == 'g_HNO3')   found_index = 2
    10523              IF ( TRIM( variable ) == 'g_NH3')    found_index = 3
    10524              IF ( TRIM( variable ) == 'g_OCNV')   found_index = 4
    10525              IF ( TRIM( variable ) == 'g_OCSV')   found_index = 5
    10526 
    10527              DO  i = nxl, nxr
    10528                 DO  j = nys, nyn
    10529                    DO  k = nzb_do, nzt_do
    10530                       local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i),                 &
    10531                                                REAL( fill_value, KIND = wp ),                      &
    10532                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10533                    ENDDO
    10534                 ENDDO
    10535              ENDDO
    10536           ELSE
    10537              IF ( TRIM( variable(3:) ) == 'H2SO4' ) to_be_resorted => g_h2so4_av
    10538              IF ( TRIM( variable(3:) ) == 'HNO3' )  to_be_resorted => g_hno3_av
    10539              IF ( TRIM( variable(3:) ) == 'NH3' )   to_be_resorted => g_nh3_av
    10540              IF ( TRIM( variable(3:) ) == 'OCNV' )  to_be_resorted => g_ocnv_av
    10541              IF ( TRIM( variable(3:) ) == 'OCSV' )  to_be_resorted => g_ocsv_av
    10542              DO  i = nxl, nxr
    10543                 DO  j = nys, nyn
    10544                    DO  k = nzb_do, nzt_do
    10545                       local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,            &
    10546                                                KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
    10547                    ENDDO
    10548                 ENDDO
    10549              ENDDO
    10550           ENDIF
    10551 
    10552        CASE ( 'LDSA' )
    10553           IF ( av == 0 )  THEN
    10554              DO  i = nxl, nxr
    10555                 DO  j = nys, nyn
    10556                    DO  k = nzb_do, nzt_do
    10557                       temp_bin = 0.0_wp
    10558                       DO  ib = 1, nbins_aerosol
    10559 !
    10560 !--                      Diameter in micrometres
    10561                          mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
    10562 !
    10563 !--                      Deposition factor: alveolar
    10564                          df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +       &
    10565                                 2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -     &
    10566                                 1.362_wp )**2 ) )
    10567 !
    10568 !--                      Number concentration in 1/cm3
    10569                          nc = 1.0E-6_wp * aerosol_number(ib)%conc(k,j,i)
    10570 !
    10571 !--                      Lung-deposited surface area LDSA (units mum2/cm3)
    10572                          temp_bin = temp_bin + pi * mean_d**2 * df * nc
    10573                       ENDDO
    10574                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10575                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10576                    ENDDO
    10577                 ENDDO
    10578              ENDDO
    10579           ELSE
    10580              DO  i = nxl, nxr
    10581                 DO  j = nys, nyn
    10582                    DO  k = nzb_do, nzt_do
    10583                       local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10584                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10585                    ENDDO
    10586                 ENDDO
    10587              ENDDO
    10588           ENDIF
    10589 
    10590        CASE ( 'N_bin1', 'N_bin2', 'N_bin3', 'N_bin4', 'N_bin5', 'N_bin6', 'N_bin7', 'N_bin8',      &
    10591               'N_bin9', 'N_bin10', 'N_bin11', 'N_bin12' )
    10592           IF ( TRIM( variable(6:) ) == '1' ) ib = 1
    10593           IF ( TRIM( variable(6:) ) == '2' ) ib = 2
    10594           IF ( TRIM( variable(6:) ) == '3' ) ib = 3
    10595           IF ( TRIM( variable(6:) ) == '4' ) ib = 4
    10596           IF ( TRIM( variable(6:) ) == '5' ) ib = 5
    10597           IF ( TRIM( variable(6:) ) == '6' ) ib = 6
    10598           IF ( TRIM( variable(6:) ) == '7' ) ib = 7
    10599           IF ( TRIM( variable(6:) ) == '8' ) ib = 8
    10600           IF ( TRIM( variable(6:) ) == '9' ) ib = 9
    10601           IF ( TRIM( variable(6:) ) == '10' ) ib = 10
    10602           IF ( TRIM( variable(6:) ) == '11' ) ib = 11
    10603           IF ( TRIM( variable(6:) ) == '12' ) ib = 12
    10604 
     10821    IF ( variable(7:11) == 'N_bin' )  THEN
     10822       READ( variable(12:),* ) char_to_int
     10823       IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10824
     10825          ib = char_to_int
    1060510826          IF ( av == 0 )  THEN
    1060610827             DO  i = nxl, nxr
     
    1062210843             ENDDO
    1062310844          ENDIF
    10624 
    10625        CASE ( 'Ntot' )
    10626           IF ( av == 0 )  THEN
    10627              DO  i = nxl, nxr
    10628                 DO  j = nys, nyn
    10629                    DO  k = nzb_do, nzt_do
    10630                       temp_bin = 0.0_wp
    10631                       DO  ib = 1, nbins_aerosol
    10632                          temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
    10633                       ENDDO
    10634                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10635                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10636                    ENDDO
    10637                 ENDDO
    10638              ENDDO
    10639           ELSE
    10640              DO  i = nxl, nxr
    10641                 DO  j = nys, nyn
    10642                    DO  k = nzb_do, nzt_do
    10643                       local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10644                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10645                    ENDDO
    10646                 ENDDO
    10647              ENDDO
    10648           ENDIF
    10649 
    10650        CASE ( 'm_bin1', 'm_bin2', 'm_bin3', 'm_bin4', 'm_bin5', 'm_bin6', 'm_bin7', 'm_bin8',      &
    10651               'm_bin9', 'm_bin10' , 'm_bin11', 'm_bin12' )
    10652           IF ( TRIM( variable(6:) ) == '1' ) ib = 1
    10653           IF ( TRIM( variable(6:) ) == '2' ) ib = 2
    10654           IF ( TRIM( variable(6:) ) == '3' ) ib = 3
    10655           IF ( TRIM( variable(6:) ) == '4' ) ib = 4
    10656           IF ( TRIM( variable(6:) ) == '5' ) ib = 5
    10657           IF ( TRIM( variable(6:) ) == '6' ) ib = 6
    10658           IF ( TRIM( variable(6:) ) == '7' ) ib = 7
    10659           IF ( TRIM( variable(6:) ) == '8' ) ib = 8
    10660           IF ( TRIM( variable(6:) ) == '9' ) ib = 9
    10661           IF ( TRIM( variable(6:) ) == '10' ) ib = 10
    10662           IF ( TRIM( variable(6:) ) == '11' ) ib = 11
    10663           IF ( TRIM( variable(6:) ) == '12' ) ib = 12
    10664 
     10845       ENDIF
     10846
     10847    ELSEIF ( variable(7:11) == 'm_bin' )  THEN
     10848       READ( variable(12:),* ) char_to_int
     10849       IF ( char_to_int >= 1  .AND. char_to_int <= SUM( nbin ) )  THEN
     10850
     10851          ib = char_to_int
    1066510852          IF ( av == 0 )  THEN
    1066610853             DO  i = nxl, nxr
     
    1068610873             ENDDO
    1068710874          ENDIF
    10688 
    10689        CASE ( 'PM2.5' )
    10690           IF ( av == 0 )  THEN
    10691              DO  i = nxl, nxr
    10692                 DO  j = nys, nyn
    10693                    DO  k = nzb_do, nzt_do
    10694                       temp_bin = 0.0_wp
    10695                       DO  ib = 1, nbins_aerosol
    10696                          IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
    10697                             DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     10875       ENDIF
     10876
     10877    ELSE
     10878       SELECT CASE ( TRIM( variable(7:) ) )
     10879
     10880          CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
     10881             IF ( av == 0 )  THEN
     10882                IF ( TRIM( variable(7:) ) == 'g_H2SO4')  found_index = 1
     10883                IF ( TRIM( variable(7:) ) == 'g_HNO3')   found_index = 2
     10884                IF ( TRIM( variable(7:) ) == 'g_NH3')    found_index = 3
     10885                IF ( TRIM( variable(7:) ) == 'g_OCNV')   found_index = 4
     10886                IF ( TRIM( variable(7:) ) == 'g_OCSV')   found_index = 5
     10887
     10888                DO  i = nxl, nxr
     10889                   DO  j = nys, nyn
     10890                      DO  k = nzb_do, nzt_do
     10891                         local_pf(i,j,k) = MERGE( salsa_gas(found_index)%conc(k,j,i),              &
     10892                                                  REAL( fill_value, KIND = wp ),                   &
     10893                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10894                      ENDDO
     10895                   ENDDO
     10896                ENDDO
     10897             ELSE
     10898!
     10899!--             9: remove salsa_g_ from the beginning
     10900                IF ( TRIM( variable(9:) ) == 'H2SO4' ) to_be_resorted => g_h2so4_av
     10901                IF ( TRIM( variable(9:) ) == 'HNO3' )  to_be_resorted => g_hno3_av
     10902                IF ( TRIM( variable(9:) ) == 'NH3' )   to_be_resorted => g_nh3_av
     10903                IF ( TRIM( variable(9:) ) == 'OCNV' )  to_be_resorted => g_ocnv_av
     10904                IF ( TRIM( variable(9:) ) == 'OCSV' )  to_be_resorted => g_ocsv_av
     10905                DO  i = nxl, nxr
     10906                   DO  j = nys, nyn
     10907                      DO  k = nzb_do, nzt_do
     10908                         local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,         &
     10909                                                  KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
     10910                      ENDDO
     10911                   ENDDO
     10912                ENDDO
     10913             ENDIF
     10914
     10915          CASE ( 'LDSA' )
     10916             IF ( av == 0 )  THEN
     10917                DO  i = nxl, nxr
     10918                   DO  j = nys, nyn
     10919                      DO  k = nzb_do, nzt_do
     10920                         temp_bin = 0.0_wp
     10921                         DO  ib = 1, nbins_aerosol
     10922   !
     10923   !--                      Diameter in micrometres
     10924                            mean_d = 1.0E+6_wp * ra_dry(k,j,i,ib) * 2.0_wp
     10925   !
     10926   !--                      Deposition factor: alveolar
     10927                            df = ( 0.01555_wp / mean_d ) * ( EXP( -0.416_wp * ( LOG( mean_d ) +    &
     10928                                   2.84_wp )**2 ) + 19.11_wp * EXP( -0.482_wp * ( LOG( mean_d ) -  &
     10929                                   1.362_wp )**2 ) )
     10930   !
     10931   !--                      Lung-deposited surface area LDSA (units mum2/cm3)
     10932                            temp_bin = temp_bin + pi * mean_d**2 * df * 1.0E-6_wp *                &
     10933                                       aerosol_number(ib)%conc(k,j,i)
     10934                         ENDDO
     10935                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10936                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10937                      ENDDO
     10938                   ENDDO
     10939                ENDDO
     10940             ELSE
     10941                DO  i = nxl, nxr
     10942                   DO  j = nys, nyn
     10943                      DO  k = nzb_do, nzt_do
     10944                         local_pf(i,j,k) = MERGE( ldsa_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10945                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10946                      ENDDO
     10947                   ENDDO
     10948                ENDDO
     10949             ENDIF
     10950
     10951          CASE ( 'N_UFP' )
     10952             IF ( av == 0 )  THEN
     10953                DO  i = nxl, nxr
     10954                   DO  j = nys, nyn
     10955                      DO  k = nzb_do, nzt_do
     10956                         temp_bin = 0.0_wp
     10957                         DO  ib = 1, nbins_aerosol
     10958                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     10959                               temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     10960                            ENDIF
     10961                         ENDDO
     10962                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10963                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10964                      ENDDO
     10965                   ENDDO
     10966                ENDDO
     10967             ELSE
     10968                DO  i = nxl, nxr
     10969                   DO  j = nys, nyn
     10970                      DO  k = nzb_do, nzt_do
     10971                         local_pf(i,j,k) = MERGE( nufp_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10972                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10973                      ENDDO
     10974                   ENDDO
     10975                ENDDO
     10976             ENDIF
     10977
     10978          CASE ( 'Ntot' )
     10979             IF ( av == 0 )  THEN
     10980                DO  i = nxl, nxr
     10981                   DO  j = nys, nyn
     10982                      DO  k = nzb_do, nzt_do
     10983                         temp_bin = 0.0_wp
     10984                         DO  ib = 1, nbins_aerosol
     10985                            temp_bin = temp_bin + aerosol_number(ib)%conc(k,j,i)
     10986                         ENDDO
     10987                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     10988                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10989                      ENDDO
     10990                   ENDDO
     10991                ENDDO
     10992             ELSE
     10993                DO  i = nxl, nxr
     10994                   DO  j = nys, nyn
     10995                      DO  k = nzb_do, nzt_do
     10996                         local_pf(i,j,k) = MERGE( ntot_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     10997                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     10998                      ENDDO
     10999                   ENDDO
     11000                ENDDO
     11001             ENDIF
     11002
     11003          CASE ( 'PM0.1' )
     11004             IF ( av == 0 )  THEN
     11005                DO  i = nxl, nxr
     11006                   DO  j = nys, nyn
     11007                      DO  k = nzb_do, nzt_do
     11008                         temp_bin = 0.0_wp
     11009                         DO  ib = 1, nbins_aerosol
     11010                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 0.1E-6_wp )  THEN
     11011                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     11012                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     11013                               ENDDO
     11014                            ENDIF
     11015                         ENDDO
     11016                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     11017                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     11018                      ENDDO
     11019                   ENDDO
     11020                ENDDO
     11021             ELSE
     11022                DO  i = nxl, nxr
     11023                   DO  j = nys, nyn
     11024                      DO  k = nzb_do, nzt_do
     11025                         local_pf(i,j,k) = MERGE( pm01_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     11026                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     11027                      ENDDO
     11028                   ENDDO
     11029                ENDDO
     11030             ENDIF
     11031
     11032          CASE ( 'PM2.5' )
     11033             IF ( av == 0 )  THEN
     11034                DO  i = nxl, nxr
     11035                   DO  j = nys, nyn
     11036                      DO  k = nzb_do, nzt_do
     11037                         temp_bin = 0.0_wp
     11038                         DO  ib = 1, nbins_aerosol
     11039                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 2.5E-6_wp )  THEN
     11040                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     11041                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     11042                               ENDDO
     11043                            ENDIF
     11044                         ENDDO
     11045                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     11046                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     11047                      ENDDO
     11048                   ENDDO
     11049                ENDDO
     11050             ELSE
     11051                DO  i = nxl, nxr
     11052                   DO  j = nys, nyn
     11053                      DO  k = nzb_do, nzt_do
     11054                         local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     11055                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     11056                      ENDDO
     11057                   ENDDO
     11058                ENDDO
     11059             ENDIF
     11060
     11061          CASE ( 'PM10' )
     11062             IF ( av == 0 )  THEN
     11063                DO  i = nxl, nxr
     11064                   DO  j = nys, nyn
     11065                      DO  k = nzb_do, nzt_do
     11066                         temp_bin = 0.0_wp
     11067                         DO  ib = 1, nbins_aerosol
     11068                            IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
     11069                               DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
     11070                                  temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
     11071                               ENDDO
     11072                            ENDIF
     11073                         ENDDO
     11074                         local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),         &
     11075                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     11076                      ENDDO
     11077                   ENDDO
     11078                ENDDO
     11079             ELSE
     11080                DO  i = nxl, nxr
     11081                   DO  j = nys, nyn
     11082                      DO  k = nzb_do, nzt_do
     11083                         local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ),   &
     11084                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
     11085                      ENDDO
     11086                   ENDDO
     11087                ENDDO
     11088             ENDIF
     11089
     11090          CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
     11091             IF ( is_used( prtcl, TRIM( variable(9:) ) ) )  THEN  ! 9: remove salsa_s_
     11092                found_index = get_index( prtcl, TRIM( variable(9:) ) )
     11093                IF ( av == 0 )  THEN
     11094                   DO  i = nxl, nxr
     11095                      DO  j = nys, nyn
     11096                         DO  k = nzb_do, nzt_do
     11097                            temp_bin = 0.0_wp
     11098                            DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
    1069811099                               temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    1069911100                            ENDDO
    10700                          ENDIF
     11101                            local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),      &
     11102                                                     BTEST( wall_flags_0(k,j,i), 0 ) )
     11103                         ENDDO
    1070111104                      ENDDO
    10702                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10703                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    1070411105                   ENDDO
    10705                 ENDDO
    10706              ENDDO
    10707           ELSE
    10708              DO  i = nxl, nxr
    10709                 DO  j = nys, nyn
    10710                    DO  k = nzb_do, nzt_do
    10711                       local_pf(i,j,k) = MERGE( pm25_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10712                                                BTEST( wall_flags_0(k,j,i), 0 ) )
     11106                ELSE
     11107!
     11108!--                9: remove salsa_s_ from the beginning
     11109                   IF ( TRIM( variable(9:) ) == 'BC' )   to_be_resorted => s_bc_av
     11110                   IF ( TRIM( variable(9:) ) == 'DU' )   to_be_resorted => s_du_av
     11111                   IF ( TRIM( variable(9:) ) == 'NH' )   to_be_resorted => s_nh_av
     11112                   IF ( TRIM( variable(9:) ) == 'NO' )   to_be_resorted => s_no_av
     11113                   IF ( TRIM( variable(9:) ) == 'OC' )   to_be_resorted => s_oc_av
     11114                   IF ( TRIM( variable(9:) ) == 'SO4' )  to_be_resorted => s_so4_av
     11115                   IF ( TRIM( variable(9:) ) == 'SS' )   to_be_resorted => s_ss_av
     11116                   DO  i = nxl, nxr
     11117                      DO  j = nys, nyn
     11118                         DO  k = nzb_do, nzt_do
     11119                            local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,      &
     11120                                                     KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
     11121                         ENDDO
     11122                      ENDDO
    1071311123                   ENDDO
    10714                 ENDDO
    10715              ENDDO
    10716           ENDIF
    10717 
    10718        CASE ( 'PM10' )
    10719           IF ( av == 0 )  THEN
    10720              DO  i = nxl, nxr
    10721                 DO  j = nys, nyn
    10722                    DO  k = nzb_do, nzt_do
    10723                       temp_bin = 0.0_wp
    10724                       DO  ib = 1, nbins_aerosol
    10725                          IF ( 2.0_wp * ra_dry(k,j,i,ib) <= 10.0E-6_wp )  THEN
    10726                             DO  ic = ib, nbins_aerosol * ncc, nbins_aerosol
    10727                                temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10728                             ENDDO
    10729                          ENDIF
    10730                       ENDDO
    10731                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10732                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10733                    ENDDO
    10734                 ENDDO
    10735              ENDDO
    10736           ELSE
    10737              DO  i = nxl, nxr
    10738                 DO  j = nys, nyn
    10739                    DO  k = nzb_do, nzt_do
    10740                       local_pf(i,j,k) = MERGE( pm10_av(k,j,i), REAL( fill_value, KIND = wp ),      &
    10741                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10742                    ENDDO
    10743                 ENDDO
    10744              ENDDO
    10745           ENDIF
    10746 
    10747        CASE ( 's_BC', 's_DU', 's_NH', 's_NO', 's_OC', 's_SO4', 's_SS' )
    10748           IF ( is_used( prtcl, TRIM( variable(3:) ) ) )  THEN
    10749              found_index = get_index( prtcl, TRIM( variable(3:) ) )
     11124                ENDIF
     11125             ENDIF
     11126
     11127          CASE ( 's_H2O' )
     11128             found_index = get_index( prtcl, 'H2O' )
    1075011129             IF ( av == 0 )  THEN
    1075111130                DO  i = nxl, nxr
     
    1076211141                ENDDO
    1076311142             ELSE
    10764                 IF ( TRIM( variable(3:) ) == 'BC' )   to_be_resorted => s_bc_av
    10765                 IF ( TRIM( variable(3:) ) == 'DU' )   to_be_resorted => s_du_av
    10766                 IF ( TRIM( variable(3:) ) == 'NH' )   to_be_resorted => s_nh_av
    10767                 IF ( TRIM( variable(3:) ) == 'NO' )   to_be_resorted => s_no_av
    10768                 IF ( TRIM( variable(3:) ) == 'OC' )   to_be_resorted => s_oc_av
    10769                 IF ( TRIM( variable(3:) ) == 'SO4' )  to_be_resorted => s_so4_av
    10770                 IF ( TRIM( variable(3:) ) == 'SS' )   to_be_resorted => s_ss_av
     11143                to_be_resorted => s_h2o_av
    1077111144                DO  i = nxl, nxr
    1077211145                   DO  j = nys, nyn
     
    1077811151                ENDDO
    1077911152             ENDIF
    10780           ENDIF
    10781 
    10782        CASE ( 's_H2O' )
    10783           found_index = get_index( prtcl, 'H2O' )
    10784           IF ( av == 0 )  THEN
    10785              DO  i = nxl, nxr
    10786                 DO  j = nys, nyn
    10787                    DO  k = nzb_do, nzt_do
    10788                       temp_bin = 0.0_wp
    10789                       DO  ic = ( found_index-1 ) * nbins_aerosol + 1, found_index * nbins_aerosol
    10790                          temp_bin = temp_bin + aerosol_mass(ic)%conc(k,j,i)
    10791                       ENDDO
    10792                       local_pf(i,j,k) = MERGE( temp_bin, REAL( fill_value, KIND = wp ),            &
    10793                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    10794                    ENDDO
    10795                 ENDDO
    10796              ENDDO
    10797           ELSE
    10798              to_be_resorted => s_h2o_av
    10799              DO  i = nxl, nxr
    10800                 DO  j = nys, nyn
    10801                    DO  k = nzb_do, nzt_do
    10802                       local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value,            &
    10803                                                KIND = wp ), BTEST( wall_flags_0(k,j,i), 0 ) )
    10804                    ENDDO
    10805                 ENDDO
    10806              ENDDO
    10807           ENDIF
    10808 
    10809        CASE DEFAULT
    10810           found = .FALSE.
    10811 
    10812     END SELECT
     11153
     11154          CASE DEFAULT
     11155             found = .FALSE.
     11156
     11157       END SELECT
     11158    ENDIF
    1081311159
    1081411160 END SUBROUTINE salsa_data_output_3d
     
    1083711183    CHARACTER(LEN=7) ::  vari      !< trimmed format of variable
    1083811184
    10839     INTEGER(iwp) ::  av              !<
    10840     INTEGER(iwp) ::  found_index     !< index of a chemical compound
    10841     INTEGER(iwp) ::  ib              !< loop index for aerosol size number bins
    10842     INTEGER(iwp) ::  ic              !< loop index for chemical components
    10843     INTEGER(iwp) ::  i               !< loop index in x-direction
    10844     INTEGER(iwp) ::  j               !< loop index in y-direction
    10845     INTEGER(iwp) ::  k               !< loop index in z-direction
     11185    INTEGER(iwp) ::  av             !<
     11186    INTEGER(iwp) ::  char_to_int    !< for converting character to integer
     11187    INTEGER(iwp) ::  found_index    !< index of a chemical compound
     11188    INTEGER(iwp) ::  ib             !< loop index for aerosol size number bins
     11189    INTEGER(iwp) ::  ic             !< loop index for chemical components
     11190    INTEGER(iwp) ::  i              !< loop index in x-direction
     11191    INTEGER(iwp) ::  j              !< loop index in y-direction
     11192    INTEGER(iwp) ::  k              !< loop index in z-direction
    1084611193    INTEGER(iwp) ::  mid            !< masked output running index
    10847     INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
     11194    INTEGER(iwp) ::  topo_top_ind   !< k index of highest horizontal surface
    1084811195
    1084911196    LOGICAL ::  found      !<
    1085011197    LOGICAL ::  resorted   !<
    1085111198
    10852     REAL(wp) ::  df        !< For calculating LDSA: fraction of particles
    10853                            !< depositing in the alveolar (or tracheobronchial)
    10854                            !< region of the lung. Depends on the particle size
     11199    REAL(wp) ::  df        !< For calculating LDSA: fraction of particles depositing in the alveolar
     11200                           !< (or tracheobronchial) region of the lung. Depends on the particle size
    1085511201    REAL(wp) ::  mean_d    !< Particle diameter in micrometres
    10856     REAL(wp) ::  nc        !< Particle number concentration in units 1/cm**3
    1085711202    REAL(wp) ::  temp_bin  !< temporary array for calculating output variables
    1085811203
    1085911204    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !<
    1086011205
     11206    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), TARGET ::  temp_array  !< temporary array
     11207
    1086111208    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< pointer
    1086211209
    10863     found     = .TRUE.
    10864     resorted  = .FALSE.
    10865     grid      = 's'
    10866     temp_bin  = 0.0_wp
    10867 
    10868     SELECT CASE ( TRIM( variable ) )
    10869 
    10870        CASE ( 'g_H2SO4', 'g_HNO3', 'g_NH3', 'g_OCNV',  'g_OCSV' )
    10871           vari = TRIM( variable )
    10872           IF ( av == 0 )  THEN
    10873              IF ( vari == 'g_H2SO4')  to_be_resorted => salsa_gas(1)%conc
    10874              IF ( vari == 'g_HNO3')   to_be_resorted => salsa_gas(2)%conc
    10875              IF ( vari == 'g_NH3')    to_be_resorted => salsa_gas(3)%conc
    10876              IF ( vari == 'g_OCNV')   to_be_resorted => salsa_gas(4)%conc
    10877              IF ( vari == 'g_OCSV')   to_be_resorted => salsa_gas(5)%conc
    10878           ELSE
    10879              IF ( vari == 'g_H2SO4') to_be_resorted => g_h2so4_av
    10880              IF ( vari == 'g_HNO3')  to_be_resorted => g_hno3_av
    10881              IF ( vari == 'g_NH3')   to_be_resorted => g_nh3_av
    10882              IF ( vari == 'g_OCNV')  to_be_resorted => g_ocnv_av
    10883              IF ( vari == 'g_OCSV')  to_be_resorted => g_ocsv_av
    10884           ENDIF
    10885 
    10886        CASE ( 'LDSA' )
    10887           IF ( av == 0 )  THEN
    10888