Changeset 3864 for palm/trunk/SOURCE/salsa_mod.f90
 Timestamp:
 Apr 5, 2019 9:01:56 AM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/salsa_mod.f90
r3833 r3864 26 26 !  27 27 ! $Id$ 28 ! Major changes in formatting, performance and data input structure (see branch 29 ! the history for details) 30 !  Timedependent emissions enabled: lod=1 for yearly PM emissions that are 31 ! normalised depending on the time, and lod=2 for preprocessed emissions 32 ! (similar to the chemistry module). 33 !  Additionally, 'uniform' emissions allowed. This emission is set constant on 34 ! all horisontal upward facing surfaces and it is created based on parameters 35 ! surface_aerosol_flux, aerosol_flux_dpg/sigmag/mass_fracs_a/mass_fracs_b. 36 !  All emissions are now implemented as surface fluxes! No 3D sources anymore. 37 !  Update the emission information by calling salsa_emission_update if 38 ! skip_time_do_salsa >= time_since_reference_point and 39 ! next_aero_emission_update <= time_since_reference_point 40 !  Aerosol background concentrations read from PIDS_DYNAMIC. The vertical grid 41 ! must match the one applied in the model. 42 !  Gas emissions and background concentrations can be also read in in salsa_mod 43 ! if the chemistry module is not applied. 44 !  In deposition, information on the land use type can be now imported from 45 ! the land use model 46 !  Use SI units in PARIN, i.e. n_lognorm given in #/m3 and dpg in metres. 47 !  Apply 100 character line limit 48 !  Change all variable names from capital to lowercase letter 49 !  Change real exponents to integer if possible. If not, precalculate the value 50 ! value of exponent 51 !  Rename in1a to start_subrange_1a, fn2a to end_subrange_1a etc. 52 !  Rename nbins > nbins_aerosol, ncc_tot > ncomponents_mass and ngast > 53 ! ngases_salsa 54 !  Rename ibc to index_bc, idu to index_du etc. 55 !  Renamed loop indices b, c and sg to ib, ic and ig 56 !  run_salsa subroutine removed 57 !  Corrected a bud in salsa_driver: falsely applied ino instead of inh 58 !  Call salsa_tendency within salsa_prognostic_equations which is called in 59 ! module_interface_mod instead of prognostic_equations_mod 60 !  Removed tailing white spaces and unused variables 61 !  Change error message to start by PA instead of SA 62 ! 63 ! 3833 20190328 15:04:04Z forkel 28 64 ! added USE chem_gasphase_mod for nvar, nspec and spc_names 29 65 ! … … 71 107 !  72 108 ! @author Mona Kurppa (University of Helsinki) 73 ! 109 ! 74 110 ! 75 111 ! Description: 76 112 !  77 !> Sectional aerosol module for large scale applications SALSA 113 !> Sectional aerosol module for large scale applications SALSA 78 114 !> (Kokkola et al., 2008, ACP 8, 24692483). Solves the aerosol number and mass 79 !> concentration as well as chemical composition. Includes aerosol dynamic 115 !> concentration as well as chemical composition. Includes aerosol dynamic 80 116 !> processes: nucleation, condensation/evaporation of vapours, coagulation and 81 !> deposition on tree leaves, ground and roofs. 117 !> deposition on tree leaves, ground and roofs. 82 118 !> Implementation is based on formulations implemented in UCLALESSALSA except 83 !> for deposition which is based on parametrisations by Zhang et al. (2001, 84 !> Atmos. Environ. 35, 549560) or Petroff&Zhang (2010, Geosci. Model Dev. 3, 119 !> for deposition which is based on parametrisations by Zhang et al. (2001, 120 !> Atmos. Environ. 35, 549560) or Petroff&Zhang (2010, Geosci. Model Dev. 3, 85 121 !> 753769) 86 122 !> 87 !> @todo Implement turbulent inflow of aerosols in inflow_turbulence. 88 !> @todo Deposition on subgrid scale vegetation 89 !> @todo Deposition on vegetation calculated by default for deciduous broadleaf 90 !> trees 91 !> @todo Revise masked data output. There is a potential bug in case of 92 !> terrainfollowing masked output, according to data_output_mask. 93 !> @todo There are now improved interfaces for NetCDF data input which can be 94 !> used instead of get variable etc. 123 !> @todo Apply information from emission_stack_height to lift emission sources 124 !> @todo emission mode "parameterized", i.e. based on street type 95 125 !! 96 126 MODULE salsa_mod … … 98 128 USE basic_constants_and_equations_mod, & 99 129 ONLY: c_p, g, p_0, pi, r_d 100 130 101 131 USE chem_gasphase_mod, & 102 132 ONLY: nspec, nvar, spc_names … … 113 143 ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, & 114 144 nzb_s_inner, nz, nzt, wall_flags_0 115 145 116 146 USE kinds 117 147 118 148 USE pegrid 119 149 120 150 USE salsa_util_mod 121 151 … … 125 155 ! 126 156 ! Local constants: 127 INTEGER(iwp), PARAMETER :: ngast = 5 !< total number of gaseous tracers: 128 !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 129 !< 4 = OCNV (nonvolatile OC), 130 !< 5 = OCSV (semivolatile) 131 INTEGER(iwp), PARAMETER :: nmod = 7 !< number of modes for initialising 132 !< the aerosol size distribution 133 INTEGER(iwp), PARAMETER :: nreg = 2 !< Number of main size subranges 134 INTEGER(iwp), PARAMETER :: maxspec = 7 !< Max. number of aerosol species 135 ! 157 INTEGER(iwp), PARAMETER :: luc_urban = 8 !< default landuse type for urban: use desert! 158 INTEGER(iwp), PARAMETER :: ngases_salsa = 5 !< total number of gaseous tracers: 159 !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV 160 !< (nonvolatile OC), 5 = OCSV (semivolatile) 161 INTEGER(iwp), PARAMETER :: nmod = 7 !< number of modes for initialising the aerosol size 162 !< distribution 163 INTEGER(iwp), PARAMETER :: nreg = 2 !< Number of main size subranges 164 INTEGER(iwp), PARAMETER :: maxspec = 7 !< Max. number of aerosol species 165 INTEGER(iwp), PARAMETER :: season = 1 !< For dry depostion by Zhang et al.: 1 = summer, 166 !< 2 = autumn (no harvest yet), 3 = late autumn 167 !< (already frost), 4 = winter, 5 = transitional spring 168 ! 136 169 ! Universal constants 137 REAL(wp), PARAMETER :: abo = 1.380662E23_wp !< Boltzmann constant (J/K) 138 REAL(wp), PARAMETER :: alv = 2.260E+6_wp !< latent heat for H2O 139 !< vaporisation (J/kg) 140 REAL(wp), PARAMETER :: alv_d_rv = 4896.96865_wp !< alv / rv 141 REAL(wp), PARAMETER :: am_airmol = 4.8096E26_wp !< Average mass of one air 142 !< molecule (Jacobson, 143 !< 2005, Eq. 2.3) 144 REAL(wp), PARAMETER :: api6 = 0.5235988_wp !< pi / 6 145 REAL(wp), PARAMETER :: argas = 8.314409_wp !< Gas constant (J/(mol K)) 146 REAL(wp), PARAMETER :: argas_d_cpd = 8.281283865E3_wp !< argas per cpd 147 REAL(wp), PARAMETER :: avo = 6.02214E+23_wp !< Avogadro constant (1/mol) 148 REAL(wp), PARAMETER :: d_sa = 5.539376964394570E10_wp !< diameter of 149 !< condensing sulphuric 150 !< acid molecule (m) 151 REAL(wp), PARAMETER :: for_ppm_to_nconc = 7.243016311E+16_wp !< 152 !< ppm * avo / R (K/(Pa*m3)) 170 REAL(wp), PARAMETER :: abo = 1.380662E23_wp !< Boltzmann constant (J/K) 171 REAL(wp), PARAMETER :: alv = 2.260E+6_wp !< latent heat for H2O 172 !< vaporisation (J/kg) 173 REAL(wp), PARAMETER :: alv_d_rv = 4896.96865_wp !< alv / rv 174 REAL(wp), PARAMETER :: am_airmol = 4.8096E26_wp !< Average mass of one air 175 !< molecule (Jacobson, 176 !< 2005, Eq. 2.3) 177 REAL(wp), PARAMETER :: api6 = 0.5235988_wp !< pi / 6 178 REAL(wp), PARAMETER :: argas = 8.314409_wp !< Gas constant (J/(mol K)) 179 REAL(wp), PARAMETER :: argas_d_cpd = 8.281283865E3_wp !< argas per cpd 180 REAL(wp), PARAMETER :: avo = 6.02214E+23_wp !< Avogadro constant (1/mol) 181 REAL(wp), PARAMETER :: d_sa = 5.539376964394570E10_wp !< diameter of condensing sulphuric 182 !< acid molecule (m) 183 REAL(wp), PARAMETER :: for_ppm_to_nconc = 7.243016311E+16_wp !< ppm * avo / R (K/(Pa*m3)) 153 184 REAL(wp), PARAMETER :: epsoc = 0.15_wp !< water uptake of organic 154 !< material 155 REAL(wp), PARAMETER :: mclim = 1.0E23_wp !< mass concentration min 156 !< limit for aerosols (kg/m3) 157 REAL(wp), PARAMETER :: n3 = 158.79_wp !< Number of H2SO4 molecules in 158 !< 3 nm cluster if d_sa=5.54e10m 159 REAL(wp), PARAMETER :: nclim = 1.0_wp !< number concentration min limit 160 !< for aerosols and gases (#/m3) 161 REAL(wp), PARAMETER :: surfw0 = 0.073_wp !< surface tension of pure water 162 !< at ~ 293 K (J/m2) 163 REAL(wp), PARAMETER :: vclim = 1.0E24_wp !< volume concentration min 164 !< limit for aerosols (m3/m3) 185 !< material 186 REAL(wp), PARAMETER :: mclim = 1.0E23_wp !< mass concentration min limit (kg/m3) 187 REAL(wp), PARAMETER :: n3 = 158.79_wp !< Number of H2SO4 molecules in 3 nm cluster 188 !< if d_sa=5.54e10m 189 REAL(wp), PARAMETER :: nclim = 1.0_wp !< number concentration min limit (#/m3) 190 REAL(wp), PARAMETER :: surfw0 = 0.073_wp !< surface tension of water at 293 K (J/m2) 191 ! 165 192 ! Molar masses in kg/mol 166 193 REAL(wp), PARAMETER :: ambc = 12.0E3_wp !< black carbon (BC) 167 194 REAL(wp), PARAMETER :: amdair = 28.970E3_wp !< dry air 168 REAL(wp), PARAMETER :: amdu = 100.E3_wp !< mineral dust 195 REAL(wp), PARAMETER :: amdu = 100.E3_wp !< mineral dust 169 196 REAL(wp), PARAMETER :: amh2o = 18.0154E3_wp !< H2O 170 197 REAL(wp), PARAMETER :: amh2so4 = 98.06E3_wp !< H2SO4 … … 176 203 REAL(wp), PARAMETER :: amoc = 150.E3_wp !< organic carbon (OC) 177 204 REAL(wp), PARAMETER :: amss = 58.44E3_wp !< sea salt (NaCl) 205 ! 178 206 ! Densities in kg/m3 179 REAL(wp), PARAMETER :: arhobc = 2000.0_wp !< black carbon 180 REAL(wp), PARAMETER :: arhodu = 2650.0_wp !< mineral dust 181 REAL(wp), PARAMETER :: arhoh2o = 1000.0_wp !< H2O 182 REAL(wp), PARAMETER :: arhoh2so4 = 1830.0_wp !< SO4 183 REAL(wp), PARAMETER :: arhohno3 = 1479.0_wp !< HNO3 184 REAL(wp), PARAMETER :: arhonh3 = 1530.0_wp !< NH3 185 REAL(wp), PARAMETER :: arhooc = 2000.0_wp !< organic carbon 186 REAL(wp), PARAMETER :: arhoss = 2165.0_wp !< sea salt (NaCl) 207 REAL(wp), PARAMETER :: arhobc = 2000.0_wp !< black carbon 208 REAL(wp), PARAMETER :: arhodu = 2650.0_wp !< mineral dust 209 REAL(wp), PARAMETER :: arhoh2o = 1000.0_wp !< H2O 210 REAL(wp), PARAMETER :: arhoh2so4 = 1830.0_wp !< SO4 211 REAL(wp), PARAMETER :: arhohno3 = 1479.0_wp !< HNO3 212 REAL(wp), PARAMETER :: arhonh3 = 1530.0_wp !< NH3 213 REAL(wp), PARAMETER :: arhooc = 2000.0_wp !< organic carbon 214 REAL(wp), PARAMETER :: arhoss = 2165.0_wp !< sea salt (NaCl) 215 ! 187 216 ! Volume of molecule in m3/# 188 217 REAL(wp), PARAMETER :: amvh2o = amh2o /avo / arhoh2o !< H2O 189 218 REAL(wp), PARAMETER :: amvh2so4 = amh2so4 / avo / arhoh2so4 !< SO4 190 REAL(wp), PARAMETER :: amvhno3 = amhno3 / avo / arhohno3 !< HNO3 191 REAL(wp), PARAMETER :: amvnh3 = amnh3 / avo / arhonh3 !< NH3 219 REAL(wp), PARAMETER :: amvhno3 = amhno3 / avo / arhohno3 !< HNO3 220 REAL(wp), PARAMETER :: amvnh3 = amnh3 / avo / arhonh3 !< NH3 192 221 REAL(wp), PARAMETER :: amvoc = amoc / avo / arhooc !< OC 193 222 REAL(wp), PARAMETER :: amvss = amss / avo / arhoss !< sea salt 194 223 ! 224 ! Constants for the dry deposition model by Petroff and Zhang (2010): 225 ! obstacle characteristic dimension "L" (cm) (plane obstacle by default) and empirical constants 226 ! C_B, C_IN, C_IM, beta_IM and C_IT for each land use category (15, as in Zhang et al. (2001)) 227 REAL(wp), DIMENSION(1:15), PARAMETER :: l_p10 = & 228 (/0.15, 4.0, 0.15, 3.0, 3.0, 0.5, 3.0, 99., 0.5, 2.0, 1.0, 99., 99., 99., 3.0/) 229 REAL(wp), DIMENSION(1:15), PARAMETER :: c_b_p10 = & 230 (/0.887, 1.262, 0.887, 1.262, 1.262, 0.996, 0.996, 99., 0.7, 0.93, 0.996, 99., 99., 99., 1.262/) 231 REAL(wp), DIMENSION(1:15), PARAMETER :: c_in_p10 = & 232 (/0.81, 0.216, 0.81, 0.216, 0.216, 0.191, 0.162, 99., 0.7, 0.14, 0.162, 99., 99., 99., 0.216/) 233 REAL(wp), DIMENSION(1:15), PARAMETER :: c_im_p10 = & 234 (/0.162, 0.13, 0.162, 0.13, 0.13, 0.191, 0.081, 99., 0.191, 0.086, 0.081, 99., 99., 99., 0.13/) 235 REAL(wp), DIMENSION(1:15), PARAMETER :: beta_im_p10 = & 236 (/0.6, 0.47, 0.6, 0.47, 0.47, 0.47, 0.47, 99., 0.6, 0.47, 0.47, 99., 99., 99., 0.47/) 237 REAL(wp), DIMENSION(1:15), PARAMETER :: c_it_p10 = & 238 (/0.0, 0.056, 0.0, 0.056, 0.056, 0.042, 0.056, 99., 0.042, 0.014, 0.056, 99., 99., 99., 0.056/) 239 ! 240 ! Constants for the dry deposition model by Zhang et al. (2001): 241 ! empirical constants "alpha" and "gamma" and characteristic radius "A" for 242 ! each land use category (15) and season (5) 243 REAL(wp), DIMENSION(1:15), PARAMETER :: alpha_z01 = & 244 (/1.0, 0.6, 1.1, 0.8, 0.8, 1.2, 1.2, 50.0, 50.0, 1.3, 2.0, 50.0, 100.0, 100.0, 1.5/) 245 REAL(wp), DIMENSION(1:15), PARAMETER :: gamma_z01 = & 246 (/0.56, 0.58, 0.56, 0.56, 0.56, 0.54, 0.54, 0.54, 0.54, 0.54, 0.54, 0.54, 0.50, 0.50, 0.56/) 247 REAL(wp), DIMENSION(1:15,1:5), PARAMETER :: A_z01 = RESHAPE( (/& 248 2.0, 5.0, 2.0, 5.0, 5.0, 2.0, 2.0, 99., 99., 10.0, 10.0, 99., 99., 99., 10.0,& ! SC1 249 2.0, 5.0, 2.0, 5.0, 5.0, 2.0, 2.0, 99., 99., 10.0, 10.0, 99., 99., 99., 10.0,& ! SC2 250 2.0, 5.0, 5.0, 10.0, 5.0, 5.0, 5.0, 99., 99., 10.0, 10.0, 99., 99., 99., 10.0,& ! SC3 251 2.0, 5.0, 5.0, 10.0, 5.0, 5.0, 5.0, 99., 99., 10.0, 10.0, 99., 99., 99., 10.0,& ! SC4 252 2.0, 5.0, 2.0, 5.0, 5.0, 2.0, 2.0, 99., 99., 10.0, 10.0, 99., 99., 99., 10.0 & ! SC5 253 /), (/ 15, 5 /) ) 254 ! Land use categories (based on Z01 but the same applies here also for P10): 255 ! 1 = evergreen needleleaf trees, 256 ! 2 = evergreen broadleaf trees, 257 ! 3 = deciduous needleleaf trees, 258 ! 4 = deciduous broadleaf trees, 259 ! 5 = mixed broadleaf and needleleaf trees (deciduous broadleaf trees for P10), 260 ! 6 = grass (short grass for P10), 261 ! 7 = crops, mixed farming, 262 ! 8 = desert, 263 ! 9 = tundra, 264 ! 10 = shrubs and interrupted woodlands (thorn shrubs for P10), 265 ! 11 = wetland with plants (long grass for P10) 266 ! 12 = ice cap and glacier, 267 ! 13 = inland water (inland lake for P10) 268 ! 14 = ocean (water for P10), 269 ! 15 = urban 270 ! 271 ! SALSA variables: 272 CHARACTER(LEN=20) :: bc_salsa_b = 'neumann' !< bottom boundary condition 273 CHARACTER(LEN=20) :: bc_salsa_t = 'neumann' !< top boundary condition 274 CHARACTER(LEN=20) :: depo_pcm_par = 'zhang2001' !< or 'petroff2010' 275 CHARACTER(LEN=20) :: depo_pcm_type = 'deciduous_broadleaf' !< leaf type 276 CHARACTER(LEN=20) :: depo_surf_par = 'zhang2001' !< or 'petroff2010' 277 CHARACTER(LEN=100) :: input_file_dynamic = 'PIDS_DYNAMIC' !< file name for dynamic input 278 CHARACTER(LEN=100) :: input_file_salsa = 'PIDS_SALSA' !< file name for emission data 279 CHARACTER(LEN=20) :: salsa_emission_mode = 'no_emission' !< 'no_emission', 'uniform', 280 !< 'parameterized', 'read_from_file' 281 282 CHARACTER(LEN=20), DIMENSION(4) :: decycle_method = & 283 (/'dirichlet','dirichlet','dirichlet','dirichlet'/) 284 !< Decycling method at horizontal boundaries 285 !< 1=left, 2=right, 3=south, 4=north 286 !< dirichlet = initial profiles for the ghost and first 3 layers 287 !< neumann = zero gradient 288 289 CHARACTER(LEN=3), DIMENSION(maxspec) :: listspec = & !< Active aerosols 290 (/'SO4',' ',' ',' ',' ',' ',' '/) 291 292 INTEGER(iwp) :: depo_pcm_type_num = 0 !< index for the dry deposition type on the plant canopy 293 INTEGER(iwp) :: dots_salsa = 0 !< starting index for salsatimeseries 294 INTEGER(iwp) :: end_subrange_1a = 1 !< last index for bin subrange 1a 295 INTEGER(iwp) :: end_subrange_2a = 1 !< last index for bin subrange 2a 296 INTEGER(iwp) :: end_subrange_2b = 1 !< last index for bin subrange 2b 297 INTEGER(iwp) :: ibc_salsa_b !< index for the bottom boundary condition 298 INTEGER(iwp) :: ibc_salsa_t !< index for the top boundary condition 299 INTEGER(iwp) :: index_bc = 1 !< index for black carbon (BC) 300 INTEGER(iwp) :: index_du = 1 !< index for dust 301 INTEGER(iwp) :: igctyp = 0 !< Initial gas concentration type 302 !< 0 = uniform (read from PARIN) 303 !< 1 = read vertical profile from an input file 304 INTEGER(iwp) :: index_nh = 1 !< index for NH3 305 INTEGER(iwp) :: index_no = 1 !< index for HNO3 306 INTEGER(iwp) :: index_oc = 1 !< index for organic carbon (OC) 307 INTEGER(iwp) :: isdtyp = 0 !< Initial size distribution type 308 !< 0 = uniform (read from PARIN) 309 !< 1 = read vertical profile of the mode number 310 !< concentration from an input file 311 INTEGER(iwp) :: index_so4 = 1 !< index for SO4 or H2SO4 312 INTEGER(iwp) :: index_ss = 1 !< index for sea salt 313 INTEGER(iwp) :: lod_gas_emissions = 0 !< level of detail of the gaseous emission data 314 INTEGER(iwp) :: nbins_aerosol = 1 !< total number of size bins 315 INTEGER(iwp) :: ncc = 1 !< number of chemical components used 316 INTEGER(iwp) :: ncomponents_mass = 1 !< total number of chemical compounds (ncc+1) 317 !< if particle water is advected) 318 INTEGER(iwp) :: nj3 = 1 !< J3 parametrization (nucleation) 319 !< 1 = condensational sink (Kerminen&Kulmala, 2002) 320 !< 2 = coagulational sink (Lehtinen et al. 2007) 321 !< 3 = coagS+selfcoagulation (Anttila et al. 2010) 322 INTEGER(iwp) :: nsnucl = 0 !< Choice of the nucleation scheme: 323 !< 0 = off 324 !< 1 = binary nucleation 325 !< 2 = activation type nucleation 326 !< 3 = kinetic nucleation 327 !< 4 = ternary nucleation 328 !< 5 = nucleation with ORGANICs 329 !< 6 = activation type of nucleation with H2SO4+ORG 330 !< 7 = heteromolecular nucleation with H2SO4*ORG 331 !< 8 = homomolecular nucleation of H2SO4 332 !< + heteromolecular nucleation with H2SO4*ORG 333 !< 9 = homomolecular nucleation of H2SO4 and ORG 334 !< + heteromolecular nucleation with H2SO4*ORG 335 INTEGER(iwp) :: start_subrange_1a = 1 !< start index for bin subranges: subrange 1a 336 INTEGER(iwp) :: start_subrange_2a = 1 !< subrange 2a 337 INTEGER(iwp) :: start_subrange_2b = 1 !< subrange 2b 338 339 INTEGER(iwp), DIMENSION(nreg) :: nbin = (/ 3, 7/) !< Number of size bins per subrange: 1 & 2 340 341 INTEGER(iwp), DIMENSION(ngases_salsa) :: gas_index_chem = & 342 (/ 1, 1, 1, 1, 1/) !< gas indices in chemistry_model_mod 343 !< 1 = H2SO4, 2 = HNO3, 3 = NH3, 4 = OCNV, 5 = OCSV 344 INTEGER(iwp), DIMENSION(ngases_salsa) :: emission_index_chem !< gas indices in the gas emission file 345 346 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: k_topo_top !< vertical index of the topography top 195 347 ! 196 348 ! SALSA switches: 197 INTEGER(iwp) :: nj3 = 1 !< J3 parametrization (nucleation) 198 !< 1 = condensational sink (Kerminen&Kulmala, 2002) 199 !< 2 = coagulational sink (Lehtinen et al. 2007) 200 !< 3 = coagS+selfcoagulation (Anttila et al. 2010) 201 INTEGER(iwp) :: nsnucl = 0 !< Choice of the nucleation scheme: 202 !< 0 = off 203 !< 1 = binary nucleation 204 !< 2 = activation type nucleation 205 !< 3 = kinetic nucleation 206 !< 4 = ternary nucleation 207 !< 5 = nucleation with ORGANICs 208 !< 6 = activation type of nucleation with 209 !< H2SO4+ORG 210 !< 7 = heteromolecular nucleation with H2SO4*ORG 211 !< 8 = homomolecular nucleation of H2SO4 + 212 !< heteromolecular nucleation with H2SO4*ORG 213 !< 9 = homomolecular nucleation of H2SO4 and ORG 214 !< +heteromolecular nucleation with H2SO4*ORG 215 LOGICAL :: advect_particle_water = .TRUE. !< advect water concentration of 216 !< particles 217 LOGICAL :: decycle_lr = .FALSE. !< Undo cyclic boundary 218 !< conditions: left and right 219 LOGICAL :: decycle_ns = .FALSE. !< north and south boundaries 220 LOGICAL :: feedback_to_palm = .FALSE. !< allow feedback due to 221 !< hydration and/or condensation 222 !< of H20 223 LOGICAL :: no_insoluble = .FALSE. !< Switch to exclude insoluble 224 !< chemical components 225 LOGICAL :: read_restart_data_salsa = .FALSE. !< read restart data for salsa 226 LOGICAL :: salsa_gases_from_chem = .FALSE. !< Transfer the gaseous 227 !< components to SALSA from 228 !< from chemistry model 229 LOGICAL :: van_der_waals_coagc = .FALSE. !< Enhancement of coagulation 230 !< kernel by van der Waals and 231 !< viscous forces 232 LOGICAL :: write_binary_salsa = .FALSE. !< read binary for salsa 349 LOGICAL :: advect_particle_water = .TRUE. !< advect water concentration of particles 350 LOGICAL :: decycle_lr = .FALSE. !< Undo cyclic boundary conditions: left and right 351 LOGICAL :: decycle_ns = .FALSE. !< north and south boundaries 352 LOGICAL :: feedback_to_palm = .FALSE. !< allow feedback due to condensation of H2O 353 LOGICAL :: nest_salsa = .FALSE. !< apply nesting for salsa 354 LOGICAL :: no_insoluble = .FALSE. !< Switch to exclude insoluble chemical components 355 LOGICAL :: read_restart_data_salsa = .FALSE. !< read restart data for salsa 356 LOGICAL :: salsa_gases_from_chem = .FALSE. !< Transfer the gaseous components to SALSA from 357 !< from chemistry model 358 LOGICAL :: van_der_waals_coagc = .FALSE. !< Enhancement of coagulation kernel by van der 359 !< Waals and viscous forces 360 LOGICAL :: write_binary_salsa = .FALSE. !< read binary for salsa 361 ! 233 362 ! Process switches: nl* is read from the NAMELIST and is NOT changed. 234 ! ls* is the switch used and will get the value of nl* 363 ! ls* is the switch used and will get the value of nl* 235 364 ! except for special circumstances (spinup period etc.) 236 LOGICAL :: nlcoag = .FALSE. !< Coagulation master switch 237 LOGICAL :: lscoag = .FALSE. !< 238 LOGICAL :: nlcnd = .FALSE. !< Condensation master switch 239 LOGICAL :: lscnd = .FALSE. !< 240 LOGICAL :: nlcndgas = .FALSE. !< Condensation of precursor gases 241 LOGICAL :: lscndgas = .FALSE. !< 242 LOGICAL :: nlcndh2oae = .FALSE. !< Condensation of H2O on aerosol 243 LOGICAL :: lscndh2oae = .FALSE. !< particles (FALSE > equilibrium calc.) 244 LOGICAL :: nldepo = .FALSE. !< Deposition master switch 245 LOGICAL :: lsdepo = .FALSE. !< 246 LOGICAL :: nldepo_topo = .FALSE. !< Deposition on vegetation master switch 247 LOGICAL :: lsdepo_topo = .FALSE. !< 248 LOGICAL :: nldepo_vege = .FALSE. !< Deposition on walls master switch 249 LOGICAL :: lsdepo_vege = .FALSE. !< 250 LOGICAL :: nldistupdate = .TRUE. !< Size distribution update master switch 251 LOGICAL :: lsdistupdate = .FALSE. !< 252 ! 253 ! SALSA variables: 254 CHARACTER (LEN=20) :: bc_salsa_b = 'neumann' !< bottom boundary condition 255 CHARACTER (LEN=20) :: bc_salsa_t = 'neumann' !< top boundary condition 256 CHARACTER (LEN=20) :: depo_vege_type = 'zhang2001' !< or 'petroff2010' 257 CHARACTER (LEN=20) :: depo_topo_type = 'zhang2001' !< or 'petroff2010' 258 CHARACTER (LEN=20), DIMENSION(4) :: decycle_method = & 259 (/'dirichlet','dirichlet','dirichlet','dirichlet'/) 260 !< Decycling method at horizontal boundaries, 261 !< 1=left, 2=right, 3=south, 4=north 262 !< dirichlet = initial size distribution and 263 !< chemical composition set for the ghost and 264 !< first three layers 265 !< neumann = zero gradient 266 CHARACTER (LEN=3), DIMENSION(maxspec) :: listspec = & !< Active aerosols 267 (/'SO4',' ',' ',' ',' ',' ',' '/) 268 CHARACTER (LEN=20) :: salsa_source_mode = 'no_source' 269 !< 'read_from_file', 270 !< 'constant' or 'no_source' 271 INTEGER(iwp) :: dots_salsa = 0 !< starting index for salsatimeseries 272 INTEGER(iwp) :: fn1a = 1 !< last index for bin subranges: subrange 1a 273 INTEGER(iwp) :: fn2a = 1 !< subrange 2a 274 INTEGER(iwp) :: fn2b = 1 !< subrange 2b 275 INTEGER(iwp), DIMENSION(ngast) :: gas_index_chem = (/ 1, 1, 1, 1, 1/) !< 276 !< Index of gaseous compounds in the chemistry 277 !< model. In SALSA, 1 = H2SO4, 2 = HNO3, 278 !< 3 = NH3, 4 = OCNV, 5 = OCSV 279 INTEGER(iwp) :: ibc_salsa_b !< 280 INTEGER(iwp) :: ibc_salsa_t !< 281 INTEGER(iwp) :: igctyp = 0 !< Initial gas concentration type 282 !< 0 = uniform (use H2SO4_init, HNO3_init, 283 !< NH3_init, OCNV_init and OCSV_init) 284 !< 1 = read vertical profile from an input file 285 INTEGER(iwp) :: in1a = 1 !< start index for bin subranges: subrange 1a 286 INTEGER(iwp) :: in2a = 1 !< subrange 2a 287 INTEGER(iwp) :: in2b = 1 !< subrange 2b 288 INTEGER(iwp) :: isdtyp = 0 !< Initial size distribution type 289 !< 0 = uniform 290 !< 1 = read vertical profile of the mode number 291 !< concentration from an input file 292 INTEGER(iwp) :: ibc = 1 !< Indice for: black carbon (BC) 293 INTEGER(iwp) :: idu = 1 !< dust 294 INTEGER(iwp) :: inh = 1 !< NH3 295 INTEGER(iwp) :: ino = 1 !< HNO3 296 INTEGER(iwp) :: ioc = 1 !< organic carbon (OC) 297 INTEGER(iwp) :: iso4 = 1 !< SO4 or H2SO4 298 INTEGER(iwp) :: iss = 1 !< sea salt 299 INTEGER(iwp) :: lod_aero = 0 !< level of detail for aerosol emissions 300 INTEGER(iwp) :: lod_gases = 0 !< level of detail for gaseous emissions 301 INTEGER(iwp), DIMENSION(nreg) :: nbin = (/ 3, 7/) !< Number of size bins 302 !< for each aerosol size subrange 303 INTEGER(iwp) :: nbins = 1 !< total number of size bins 304 INTEGER(iwp) :: ncc = 1 !< number of chemical components used 305 INTEGER(iwp) :: ncc_tot = 1!< total number of chemical compounds (ncc+1 306 !< if particle water is advected) 307 REAL(wp) :: act_coeff = 1.0E7_wp !< Activation coefficient 308 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: emission_mass_fracs !< array for 309 !< aerosol composition per emission category 310 !< 1:SO4 2:OC 3:BC 4:DU 5:SS 6:NO 7:NH 311 REAL(wp) :: dt_salsa = 0.00001_wp !< Time step of SALSA 312 REAL(wp) :: H2SO4_init = nclim !< Init value for sulphuric acid gas 313 REAL(wp) :: HNO3_init = nclim !< Init value for nitric acid gas 314 REAL(wp) :: last_salsa_time = 0.0_wp !< time of the previous salsa 315 !< timestep 316 REAL(wp) :: nf2a = 1.0_wp !< Number fraction allocated to a 317 !< bins in subrange 2 318 !< (bbins will get 1nf2a) 319 REAL(wp) :: NH3_init = nclim !< Init value for ammonia gas 320 REAL(wp) :: OCNV_init = nclim !< Init value for nonvolatile 321 !< organic gases 322 REAL(wp) :: OCSV_init = nclim !< Init value for semivolatile 323 !< organic gases 365 LOGICAL :: nlcoag = .FALSE. !< Coagulation master switch 366 LOGICAL :: lscoag = .FALSE. !< 367 LOGICAL :: nlcnd = .FALSE. !< Condensation master switch 368 LOGICAL :: lscnd = .FALSE. !< 369 LOGICAL :: nlcndgas = .FALSE. !< Condensation of precursor gases 370 LOGICAL :: lscndgas = .FALSE. !< 371 LOGICAL :: nlcndh2oae = .FALSE. !< Condensation of H2O on aerosol 372 LOGICAL :: lscndh2oae = .FALSE. !< particles (FALSE > equilibrium calc.) 373 LOGICAL :: nldepo = .FALSE. !< Deposition master switch 374 LOGICAL :: lsdepo = .FALSE. !< 375 LOGICAL :: nldepo_surf = .FALSE. !< Deposition on vegetation master switch 376 LOGICAL :: lsdepo_surf = .FALSE. !< 377 LOGICAL :: nldepo_pcm = .FALSE. !< Deposition on walls master switch 378 LOGICAL :: lsdepo_pcm = .FALSE. !< 379 LOGICAL :: nldistupdate = .TRUE. !< Size distribution update master switch 380 LOGICAL :: lsdistupdate = .FALSE. !< 381 LOGICAL :: lspartition = .FALSE. !< Partition of HNO3 and NH3 382 383 REAL(wp) :: act_coeff = 1.0E7_wp !< Activation coefficient 384 REAL(wp) :: dt_salsa = 0.00001_wp !< Time step of SALSA 385 REAL(wp) :: h2so4_init = nclim !< Init value for sulphuric acid gas 386 REAL(wp) :: hno3_init = nclim !< Init value for nitric acid gas 387 REAL(wp) :: last_salsa_time = 0.0_wp !< previous salsa call 388 REAL(wp) :: next_aero_emission_update = 0.0_wp !< previous emission update 389 REAL(wp) :: next_gas_emission_update = 0.0_wp !< previous emission update 390 REAL(wp) :: nf2a = 1.0_wp !< Number fraction allocated to 2abins 391 REAL(wp) :: nh3_init = nclim !< Init value for ammonia gas 392 REAL(wp) :: ocnv_init = nclim !< Init value for nonvolatile organic gases 393 REAL(wp) :: ocsv_init = nclim !< Init value for semivolatile organic gases 394 REAL(wp) :: rhlim = 1.20_wp !< RH limit in %/100. Prevents unrealistical RH 395 REAL(wp) :: skip_time_do_salsa = 0.0_wp !< Starting time of SALSA (s) 396 ! 397 ! Initial lognormal size distribution: mode diameter (dpg, metres), 398 ! standard deviation (sigmag) and concentration (n_lognorm, #/m3) 399 REAL(wp), DIMENSION(nmod) :: dpg = & 400 (/0.013_wp, 0.054_wp, 0.86_wp, 0.2_wp, 0.2_wp, 0.2_wp, 0.2_wp/) 401 REAL(wp), DIMENSION(nmod) :: sigmag = & 402 (/1.8_wp, 2.16_wp, 2.21_wp, 2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/) 403 REAL(wp), DIMENSION(nmod) :: n_lognorm = & 404 (/1.04e+11_wp, 3.23E+10_wp, 5.4E+6_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) 405 ! 406 ! Initial mass fractions / chemical composition of the size distribution 407 REAL(wp), DIMENSION(maxspec) :: mass_fracs_a = & !< mass fractions between 408 (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for A bins 409 REAL(wp), DIMENSION(maxspec) :: mass_fracs_b = & !< mass fractions between 410 (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for B bins 324 411 REAL(wp), DIMENSION(nreg+1) :: reglim = & !< Min&max diameters of size subranges 325 412 (/ 3.0E9_wp, 5.0E8_wp, 1.0E5_wp/) 326 REAL(wp) :: rhlim = 1.20_wp !< RH limit in %/100. Prevents 327 !< unrealistically high RH in condensation 328 REAL(wp) :: skip_time_do_salsa = 0.0_wp !< Starting time of SALSA (s) 329 ! Initial lognormal size distribution: mode diameter (dpg, micrometres), 330 ! standard deviation (sigmag) and concentration (n_lognorm, #/cm3) 331 REAL(wp), DIMENSION(nmod) :: dpg = (/0.013_wp, 0.054_wp, 0.86_wp, & 332 0.2_wp, 0.2_wp, 0.2_wp, 0.2_wp/) 333 REAL(wp), DIMENSION(nmod) :: sigmag = (/1.8_wp, 2.16_wp, 2.21_wp, & 334 2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/) 335 REAL(wp), DIMENSION(nmod) :: n_lognorm = (/1.04e+5_wp, 3.23E+4_wp, 5.4_wp,& 336 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) 337 ! Initial mass fractions / chemical composition of the size distribution 338 REAL(wp), DIMENSION(maxspec) :: mass_fracs_a = & !< mass fractions between 339 (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for A bins 340 REAL(wp), DIMENSION(maxspec) :: mass_fracs_b = & !< mass fractions between 341 (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) !< aerosol species for B bins 342 343 REAL(wp), ALLOCATABLE, DIMENSION(:) :: bin_low_limits !< to deliver 344 !< information about 345 !< the lower 346 !< diameters per bin 347 REAL(wp), ALLOCATABLE, DIMENSION(:) :: nsect !< Background number 348 !< concentration per bin 349 REAL(wp), ALLOCATABLE, DIMENSION(:) :: massacc !< Mass accomodation 350 !< coefficients per bin 351 ! 352 ! SALSA derived datatypes: 353 ! 354 ! Prognostic variable: Aerosol size bin information (number (#/m3) and 355 ! mass (kg/m3) concentration) and the concentration of gaseous tracers (#/m3). 356 ! Gas tracers are contained sequentially in dimension 4 as: 357 ! 1. H2SO4, 2. HNO3, 3. NH3, 4. OCNV (nonvolatile organics), 358 ! 5. OCSV (semivolatile) 413 ! 414 ! Initial lognormal size distribution: mode diameter (dpg, metres), standard deviation (sigmag) 415 ! concentration (n_lognorm, #/m3) and mass fractions of all chemical components (listed in 416 ! listspec) for both a (soluble) and b (insoluble) bins. 417 REAL(wp), DIMENSION(nmod) :: aerosol_flux_dpg = & 418 (/0.013_wp, 0.054_wp, 0.86_wp, 0.2_wp, 0.2_wp, 0.2_wp, 0.2_wp/) 419 REAL(wp), DIMENSION(nmod) :: aerosol_flux_sigmag = & 420 (/1.8_wp, 2.16_wp, 2.21_wp, 2.0_wp, 2.0_wp, 2.0_wp, 2.0_wp/) 421 REAL(wp), DIMENSION(nmod) :: surface_aerosol_flux = & 422 (/1.04e+11_wp, 3.23E+10_wp, 5.4E+6_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) 423 REAL(wp), DIMENSION(maxspec) :: aerosol_flux_mass_fracs_a = & 424 (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) 425 REAL(wp), DIMENSION(maxspec) :: aerosol_flux_mass_fracs_b = & 426 (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/) 427 428 REAL(wp), DIMENSION(:), ALLOCATABLE :: bin_low_limits !< to deliver information about 429 !< the lower diameters per bin 430 REAL(wp), DIMENSION(:), ALLOCATABLE :: bc_am_t_val !< vertical gradient of: aerosol mass 431 REAL(wp), DIMENSION(:), ALLOCATABLE :: bc_an_t_val !< of: aerosol number 432 REAL(wp), DIMENSION(:), ALLOCATABLE :: bc_gt_t_val !< salsa gases near domain top 433 REAL(wp), DIMENSION(:), ALLOCATABLE :: gas_emission_time !< Time array in gas emission data (s) 434 REAL(wp), DIMENSION(:), ALLOCATABLE :: nsect !< Background number concentrations 435 REAL(wp), DIMENSION(:), ALLOCATABLE :: massacc !< Mass accomodation coefficients 436 ! 437 ! SALSA derived datatypes: 438 ! 439 ! For matching LSM and the deposition module surface types 440 TYPE match_lsm_depo 441 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: match 442 END TYPE match_lsm_depo 443 ! 444 ! Aerosol emission data attributes 445 TYPE salsa_emission_attribute_type 446 447 CHARACTER(LEN=25) :: units 448 449 CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cat_name !< 450 CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cc_name !< 451 CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: unit_time !< 452 CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names !< 453 454 INTEGER(iwp) :: lod = 0 !< level of detail 455 INTEGER(iwp) :: nbins = 10 !< number of aerosol size bins 456 INTEGER(iwp) :: ncat = 0 !< number of emission categories 457 INTEGER(iwp) :: ncc = 7 !< number of aerosol chemical components 458 INTEGER(iwp) :: nhoursyear = 0 !< number of hours: HOURLY mode 459 INTEGER(iwp) :: nmonthdayhour = 0 !< number of month days and hours: MDH mode 460 INTEGER(iwp) :: num_vars !< number of variables 461 INTEGER(iwp) :: nt = 0 !< number of time steps 462 INTEGER(iwp) :: nz = 0 !< number of vertical levels 463 INTEGER(iwp) :: tind !< time index for reference time in salsa emission data 464 465 INTEGER(iwp), DIMENSION(maxspec) :: cc_input_to_model !< 466 467 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: cat_index !< Index of emission categories 468 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: cc_index !< Index of chemical components 469 470 REAL(wp) :: conversion_factor !< unit conversion factor for aerosol emissions 471 472 REAL(wp), DIMENSION(:), ALLOCATABLE :: dmid !< mean diameters of size bins (m) 473 REAL(wp), DIMENSION(:), ALLOCATABLE :: rho !< average density (kg/m3) 474 REAL(wp), DIMENSION(:), ALLOCATABLE :: time !< time (s) 475 REAL(wp), DIMENSION(:), ALLOCATABLE :: time_factor !< emission time factor 476 REAL(wp), DIMENSION(:), ALLOCATABLE :: z !< height (m) 477 478 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: etf !< emission time factor 479 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: stack_height 480 481 END TYPE salsa_emission_attribute_type 482 ! 483 ! The default size distribution and mass composition per emission category: 484 ! 1 = traffic, 2 = road dust, 3 = wood combustion, 4 = other 485 ! Mass fractions: H2SO4, OC, BC, DU, SS, HNO3, NH3 486 TYPE salsa_emission_mode_type 487 488 INTEGER(iwp) :: ndm = 3 !< number of default modes 489 INTEGER(iwp) :: ndc = 4 !< number of default categories 490 491 CHARACTER(LEN=25), DIMENSION(1:4) :: cat_name_table = (/'traffic exhaust', & 492 'road dust ', & 493 'wood combustion', & 494 'other '/) 495 496 INTEGER(iwp), DIMENSION(1:4) :: cat_input_to_model !< 497 498 REAL(wp), DIMENSION(1:3) :: dpg_table = (/ 13.5E9_wp, 1.4E6_wp, 5.4E8_wp/) !< 499 REAL(wp), DIMENSION(1:3) :: ntot_table !< 500 REAL(wp), DIMENSION(1:3) :: sigmag_table = (/ 1.6_wp, 1.4_wp, 1.7_wp /) !< 501 502 REAL(wp), DIMENSION(1:maxspec,1:4) :: mass_frac_table = & !< 503 RESHAPE( (/ 0.04_wp, 0.48_wp, 0.48_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 504 0.0_wp, 0.05_wp, 0.0_wp, 0.95_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 505 0.0_wp, 0.5_wp, 0.5_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 506 0.0_wp, 0.5_wp, 0.5_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp & 507 /), (/maxspec,4/) ) 508 509 REAL(wp), DIMENSION(1:3,1:4) :: pm_frac_table = & !< rel. mass 510 RESHAPE( (/ 0.016_wp, 0.000_wp, 0.984_wp, & 511 0.000_wp, 1.000_wp, 0.000_wp, & 512 0.000_wp, 0.000_wp, 1.000_wp, & 513 1.000_wp, 0.000_wp, 1.000_wp & 514 /), (/3,4/) ) 515 516 END TYPE salsa_emission_mode_type 517 ! 518 ! Aerosol emission data values 519 TYPE salsa_emission_value_type 520 521 REAL(wp) :: fill !< fill value 522 523 REAL(wp), DIMENSION(:), ALLOCATABLE :: preproc_mass_fracs !< mass fractions 524 525 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: def_mass_fracs !< mass fractions per emis. category 526 527 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: def_data !< surface emission values in PM 528 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: preproc_data !< surface emission values per bin 529 530 END TYPE salsa_emission_value_type 531 ! 532 ! Prognostic variable: Aerosol size bin information (number (#/m3) and mass (kg/m3) concentration) 533 ! and the concentration of gaseous tracers (#/m3). Gas tracers are contained sequentially in 534 ! dimension 4 as: 535 ! 1. H2SO4, 2. HNO3, 3. NH3, 4. OCNV (nonvolatile organics), 5. OCSV (semivolatile) 359 536 TYPE salsa_variable 360 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: conc 361 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: conc_p 362 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: tconc_m 363 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: flux_s, diss_s 364 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: flux_l, diss_l 365 REAL(wp), ALLOCATABLE, DIMENSION(:) :: init 366 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: source 367 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: sums_ws_l 537 538 REAL(wp), ALLOCATABLE, DIMENSION(:) :: init !< 539 540 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diss_s !< 541 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: flux_s !< 542 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: source !< 543 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: sums_ws_l !< 544 545 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: diss_l !< 546 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: flux_l !< 547 548 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: conc !< 549 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: conc_p !< 550 REAL(wp), POINTER, DIMENSION(:,:,:), CONTIGUOUS :: tconc_m !< 551 368 552 END TYPE salsa_variable 369 370 ! Map bin indices between parallel size distributions 371 TYPE t_parallelbin 372 INTEGER(iwp) :: cur ! Index for current distribution 373 INTEGER(iwp) :: par ! Index for corresponding parallel distribution 374 END TYPE t_parallelbin 375 376 ! Datatype used to store information about the binned size distributions of 377 ! aerosols 553 ! 554 ! Datatype used to store information about the binned size distributions of aerosols 378 555 TYPE t_section 556 557 REAL(wp) :: dmid !< bin middle diameter (m) 379 558 REAL(wp) :: vhilim !< bin volume at the high limit 380 559 REAL(wp) :: vlolim !< bin volume at the low limit 381 560 REAL(wp) :: vratiohi !< volume ratio between the center and high limit 382 561 REAL(wp) :: vratiolo !< volume ratio between the center and low limit 383 REAL(wp) :: dmid !< bin middle diameter (m)384 562 !****************************************************** 385 563 ! ^ Do NOT change the stuff above after initialization ! 386 564 !****************************************************** 565 REAL(wp) :: core !< Volume of dry particle 387 566 REAL(wp) :: dwet !< Wet diameter or mean droplet diameter (m) 388 REAL(wp), DIMENSION(maxspec+1) :: volc !< Volume concentrations 389 !< (m^3/m^3) of aerosols + water. Since most of 390 !< the stuff in SALSA is hard coded, these *have to 391 !< be* in the order 392 !< 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O 567 REAL(wp) :: numc !< Number concentration of particles/droplets (#/m3) 393 568 REAL(wp) :: veqh2o !< Equilibrium H2O concentration for each particle 394 REAL(wp) :: numc !< Number concentration of particles/droplets (#/m3) 395 REAL(wp) :: core !< Volume of dry particle 396 END TYPE t_section 397 ! 398 ! Local aerosol properties in SALSA 399 TYPE(t_section), ALLOCATABLE :: aero(:) 400 ! 401 ! SALSA tracers: 402 ! Tracers as x = x(k,j,i,bin). The 4th dimension contains all the size bins 403 ! sequentially for each aerosol species + water. 404 ! 405 ! Prognostic tracers: 406 ! 407 ! Number concentration (#/m3) 408 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: aerosol_number 409 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_1 410 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_2 411 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_3 412 ! 413 ! Mass concentration (kg/m3) 414 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: aerosol_mass 415 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_1 416 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_2 417 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_3 418 ! 419 ! Gaseous tracers (#/m3) 420 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: salsa_gas 421 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_1 422 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_2 423 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_3 569 570 REAL(wp), DIMENSION(maxspec+1) :: volc !< Volume concentrations (m^3/m^3) of aerosols + 571 !< water. Since most of the stuff in SALSA is hard 572 !< coded, these *have to be* in the order 573 !< 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O 574 END TYPE t_section 575 576 TYPE(salsa_emission_attribute_type) :: aero_emission_att !< emission attributes 577 TYPE(salsa_emission_value_type) :: aero_emission !< emission values 578 TYPE(salsa_emission_mode_type) :: def_modes !< default emission modes 579 580 TYPE(t_section), DIMENSION(:), ALLOCATABLE :: aero !< local aerosol properties 581 582 TYPE(match_lsm_depo) :: lsm_to_depo_h 583 584 TYPE(match_lsm_depo), DIMENSION(0:3) :: lsm_to_depo_v 585 ! 586 ! SALSA variables: as x = x(k,j,i,bin). 587 ! The 4th dimension contains all the size bins sequentially for each aerosol species + water. 588 ! 589 ! Prognostic variables: 590 ! 591 ! Number concentration (#/m3) 592 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: aerosol_number !< 593 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_1 !< 594 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_2 !< 595 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nconc_3 !< 596 ! 597 ! Mass concentration (kg/m3) 598 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: aerosol_mass !< 599 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_1 !< 600 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_2 !< 601 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mconc_3 !< 602 ! 603 ! Gaseous concentrations (#/m3) 604 TYPE(salsa_variable), ALLOCATABLE, DIMENSION(:), TARGET :: salsa_gas !< 605 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_1 !< 606 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_2 !< 607 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: gconc_3 !< 424 608 ! 425 609 ! Diagnostic tracers 426 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: sedim_vd !< sedimentation 427 !< velocity per size 428 !< bin (m/s) 429 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: Ra_dry !< dry radius (m) 430 610 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: sedim_vd !< sedimentation velocity per bin (m/s) 611 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ra_dry !< aerosol dry radius (m) 612 431 613 ! Particle component index tables 432 TYPE(component_index) :: prtcl !< Contains "getIndex" which gives the index 433 !< for a given aerosol component name, i.e. 434 !< 1:SO4, 2:OC, 3:BC, 4:DU, 435 !< 5:SS, 6:NO, 7:NH, 8:H2O 436 ! 614 TYPE(component_index) :: prtcl !< Contains "getIndex" which gives the index for a given aerosol 615 !< component name: 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O 616 ! 437 617 ! Data output arrays: 618 ! 438 619 ! Gases: 439 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_H2SO4_av !< H2SO4 440 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_HNO3_av !< HNO3 441 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_NH3_av !< NH3 442 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_OCNV_av !< nonvola 443 !< tile OC 444 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_OCSV_av !< semivol. 445 !< OC 446 ! Integrated: 447 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: LDSA_av !< lung 448 !< deposited 449 !< surface area 450 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: Ntot_av !< total number 451 !< conc. 452 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: PM25_av !< PM2.5 453 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: PM10_av !< PM10 454 ! In the particle phase: 455 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_BC_av !< black carbon 456 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_DU_av !< dust 457 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_H2O_av !< liquid water 458 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_NH_av !< ammonia 459 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_NO_av !< nitrates 460 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_OC_av !< org. carbon 461 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_SO4_av !< sulphates 462 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_SS_av !< sea salt 463 ! Bins: 464 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mbins_av !< bin mass 465 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: Nbins_av !< bin number 466 467 620 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_h2so4_av !< H2SO4 621 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_hno3_av !< HNO3 622 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_nh3_av !< NH3 623 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_ocnv_av !< nonvolatile OC 624 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: g_ocsv_av !< semivolatile OC 625 ! 626 ! Integrated: 627 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: ldsa_av !< lungdeposited surface area 628 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: ntot_av !< total number concentration 629 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: pm25_av !< PM2.5 630 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: pm10_av !< PM10 631 ! 632 ! In the particle phase: 633 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_bc_av !< black carbon 634 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_du_av !< dust 635 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_h2o_av !< liquid water 636 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_nh_av !< ammonia 637 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_no_av !< nitrates 638 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_oc_av !< org. carbon 639 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_so4_av !< sulphates 640 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: s_ss_av !< sea salt 641 ! 642 ! Bin specific mass and number concentrations: 643 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: mbins_av !< bin mas 644 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: nbins_av !< bin number 645 468 646 ! 469 647 ! PALM interfaces: … … 474 652 MODULE PROCEDURE salsa_boundary_conds_decycle 475 653 END INTERFACE salsa_boundary_conds 476 ! 654 ! 477 655 ! Data output checks for 2D/3D data to be done in check_parameters 478 656 INTERFACE salsa_check_data_output 479 657 MODULE PROCEDURE salsa_check_data_output 480 658 END INTERFACE salsa_check_data_output 481 482 659 ! 483 660 ! Input parameter checks to be done in check_parameters … … 485 662 MODULE PROCEDURE salsa_check_parameters 486 663 END INTERFACE salsa_check_parameters 487 488 664 ! 489 665 ! Averaging of 3D data for output … … 491 667 MODULE PROCEDURE salsa_3d_data_averaging 492 668 END INTERFACE salsa_3d_data_averaging 493 494 669 ! 495 670 ! Data output of 2D quantities … … 497 672 MODULE PROCEDURE salsa_data_output_2d 498 673 END INTERFACE salsa_data_output_2d 499 500 674 ! 501 675 ! Data output of 3D data … … 503 677 MODULE PROCEDURE salsa_data_output_3d 504 678 END INTERFACE salsa_data_output_3d 505 506 679 ! 507 680 ! Data output of 3D data … … 509 682 MODULE PROCEDURE salsa_data_output_mask 510 683 END INTERFACE salsa_data_output_mask 511 512 684 ! 513 685 ! Definition of data output quantities … … 515 687 MODULE PROCEDURE salsa_define_netcdf_grid 516 688 END INTERFACE salsa_define_netcdf_grid 517 518 689 ! 519 690 ! Output of information to the header file … … 521 692 MODULE PROCEDURE salsa_header 522 693 END INTERFACE salsa_header 523 524 ! 525 ! Initialization actions 694 ! 695 ! Initialization actions 526 696 INTERFACE salsa_init 527 697 MODULE PROCEDURE salsa_init 528 698 END INTERFACE salsa_init 529 530 699 ! 531 700 ! Initialization of arrays … … 533 702 MODULE PROCEDURE salsa_init_arrays 534 703 END INTERFACE salsa_init_arrays 535 536 704 ! 537 705 ! Writing of binary output for restart runs !!! renaming?! … … 539 707 MODULE PROCEDURE salsa_wrd_local 540 708 END INTERFACE salsa_wrd_local 541 542 709 ! 543 710 ! Reading of NAMELIST parameters … … 545 712 MODULE PROCEDURE salsa_parin 546 713 END INTERFACE salsa_parin 547 548 714 ! 549 715 ! Reading of parameters for restart runs … … 551 717 MODULE PROCEDURE salsa_rrd_local 552 718 END INTERFACE salsa_rrd_local 553 554 719 ! 555 720 ! Swapping of time levels (required for prognostic variables) … … 557 722 MODULE PROCEDURE salsa_swap_timelevel 558 723 END INTERFACE salsa_swap_timelevel 559 724 ! 725 ! Interface between PALM and salsa 560 726 INTERFACE salsa_driver 561 727 MODULE PROCEDURE salsa_driver 562 728 END INTERFACE salsa_driver 563 729 ! 730 ! Prognostics equations for salsa variables 731 INTERFACE salsa_prognostic_equations 732 MODULE PROCEDURE salsa_prognostic_equations 733 MODULE PROCEDURE salsa_prognostic_equations_ij 734 END INTERFACE salsa_prognostic_equations 735 ! 736 ! Tendency salsa variables 564 737 INTERFACE salsa_tendency 565 738 MODULE PROCEDURE salsa_tendency 566 739 MODULE PROCEDURE salsa_tendency_ij 567 740 END INTERFACE salsa_tendency 568 569 570 741 742 571 743 SAVE 572 744 573 745 PRIVATE 574 746 ! 575 ! Public functions: 576 PUBLIC salsa_boundary_conds, salsa_check_data_output, & 577 salsa_check_parameters, salsa_3d_data_averaging, & 578 salsa_data_output_2d, salsa_data_output_3d, salsa_data_output_mask, & 579 salsa_define_netcdf_grid, salsa_diagnostics, salsa_driver, & 580 salsa_header, salsa_init, salsa_init_arrays, salsa_parin, & 581 salsa_rrd_local, salsa_swap_timelevel, salsa_tendency, & 582 salsa_wrd_local 747 ! Public functions: 748 PUBLIC salsa_boundary_conds, salsa_check_data_output, salsa_check_parameters, & 749 salsa_3d_data_averaging, salsa_data_output_2d, salsa_data_output_3d, & 750 salsa_data_output_mask, salsa_define_netcdf_grid, salsa_diagnostics, salsa_driver, & 751 salsa_emission_update, salsa_header, salsa_init, salsa_init_arrays, salsa_parin, & 752 salsa_rrd_local, salsa_swap_timelevel, salsa_prognostic_equations, salsa_wrd_local 583 753 ! 584 754 ! Public parameters, constants and initial values 585 PUBLIC dots_salsa, dt_salsa, last_salsa_time, lsdepo, salsa, & 586 salsa_gases_from_chem, skip_time_do_salsa 755 PUBLIC bc_am_t_val, bc_an_t_val, bc_gt_t_val, dots_salsa, dt_salsa, & 756 ibc_salsa_b, last_salsa_time, lsdepo, nest_salsa, salsa, salsa_gases_from_chem, & 757 skip_time_do_salsa 587 758 ! 588 759 ! Public prognostic variables 589 PUBLIC aerosol_mass, aerosol_number, fn2a, fn2b, gconc_2, in1a, in2b, & 590 mconc_2, nbins, ncc, ncc_tot, nclim, nconc_2, ngast, prtcl, Ra_dry, & 591 salsa_gas, sedim_vd 592 760 PUBLIC aerosol_mass, aerosol_number, gconc_2, mconc_2, nbins_aerosol, ncc, ncomponents_mass, & 761 nclim, nconc_2, ngases_salsa, prtcl, ra_dry, salsa_gas, sedim_vd 762 593 763 594 764 CONTAINS … … 603 773 IMPLICIT NONE 604 774 605 CHARACTER (LEN=80) :: line !< dummy string that contains the current line 606 !< of the parameter file 607 608 NAMELIST /salsa_parameters/ & 609 advect_particle_water, & ! Switch for advecting 610 ! particle water. If .FALSE., 611 ! equilibration is called at 612 ! each time step. 613 bc_salsa_b, & ! bottom boundary condition 614 bc_salsa_t, & ! top boundary condition 615 decycle_lr, & ! decycle SALSA components 616 decycle_method, & ! decycle method applied: 617 ! 1=left 2=right 3=south 4=north 618 decycle_ns, & ! decycle SALSA components 619 depo_vege_type, & ! Parametrisation type 620 depo_topo_type, & ! Parametrisation type 621 dpg, & ! Mean diameter for the initial 622 ! lognormal modes 623 dt_salsa, & ! SALSA timestep in seconds 624 feedback_to_palm, & ! allow feedback due to 625 ! hydration / condensation 626 H2SO4_init, & ! Init value for sulphuric acid 627 HNO3_init, & ! Init value for nitric acid 628 igctyp, & ! Initial gas concentration type 629 isdtyp, & ! Initial size distribution type 630 listspec, & ! List of actived aerosols 631 ! (string list) 632 mass_fracs_a, & ! Initial relative contribution 633 ! of each species to particle 634 ! volume in abins, 0 for unused 635 mass_fracs_b, & ! Initial relative contribution 636 ! of each species to particle 637 ! volume in bbins, 0 for unused 638 n_lognorm, & ! Number concentration for the 639 ! lognormal modes 640 nbin, & ! Number of size bins for 641 ! aerosol size subranges 1 & 2 642 nf2a, & ! Number fraction of particles 643 ! allocated to abins in 644 ! subrange 2 bbins will get 645 ! 1nf2a 646 NH3_init, & ! Init value for ammonia 647 nj3, & ! J3 parametrization 648 ! 1 = condensational sink 649 ! (Kerminen&Kulmala, 2002) 650 ! 2 = coagulational sink 651 ! (Lehtinen et al. 2007) 652 ! 3 = coagS+selfcoagulation 653 ! (Anttila et al. 2010) 654 nlcnd, & ! Condensation master switch 655 nlcndgas, & ! Condensation of gases 656 nlcndh2oae, & ! Condensation of H2O 657 nlcoag, & ! Coagulation master switch 658 nldepo, & ! Deposition master switch 659 nldepo_vege, & ! Deposition on vegetation 660 ! master switch 661 nldepo_topo, & ! Deposition on topo master 662 ! switch 663 nldistupdate, & ! Size distribution update 664 ! master switch 665 nsnucl, & ! Nucleation scheme: 666 ! 0 = off, 667 ! 1 = binary nucleation 668 ! 2 = activation type nucleation 669 ! 3 = kinetic nucleation 670 ! 4 = ternary nucleation 671 ! 5 = nucleation with organics 672 ! 6 = activation type of 673 ! nucleation with H2SO4+ORG 674 ! 7 = heteromolecular nucleation 675 ! with H2SO4*ORG 676 ! 8 = homomolecular nucleation 677 ! of H2SO4 + heteromolecular 678 ! nucleation with H2SO4*ORG 679 ! 9 = homomolecular nucleation 680 ! of H2SO4 and ORG + hetero 681 ! molecular nucleation with 682 ! H2SO4*ORG 683 OCNV_init, & ! Init value for nonvolatile 684 ! organic gases 685 OCSV_init, & ! Init value for semivolatile 686 ! organic gases 687 read_restart_data_salsa, & ! read restart data for 688 ! salsa 689 reglim, & ! Min&max diameter limits of 690 ! size subranges 691 salsa, & ! Master switch for SALSA 692 salsa_source_mode,& ! 'read_from_file' or 'constant' 693 ! or 'no_source' 694 sigmag, & ! stdev for the initial log 695 ! normal modes 696 skip_time_do_salsa, & ! Starting time of SALSA (s) 697 van_der_waals_coagc,& ! include van der Waals forces 698 write_binary_salsa ! Write binary for salsa 699 700 775 CHARACTER(LEN=80) :: line !< dummy string that contains the current line 776 !< of the parameter file 777 778 NAMELIST /salsa_parameters/ aerosol_flux_dpg, aerosol_flux_mass_fracs_a, & 779 aerosol_flux_mass_fracs_b, aerosol_flux_sigmag, & 780 advect_particle_water, bc_salsa_b, bc_salsa_t, decycle_lr, & 781 decycle_method, decycle_ns, depo_pcm_par, depo_pcm_type, & 782 depo_surf_par, dpg, dt_salsa, feedback_to_palm, h2so4_init, & 783 hno3_init, igctyp, isdtyp, listspec, mass_fracs_a, & 784 mass_fracs_b, n_lognorm, nbin, nest_salsa, nf2a, nh3_init, & 785 nj3, nlcnd, nlcndgas, nlcndh2oae, nlcoag, nldepo, nldepo_pcm, & 786 nldepo_surf, nldistupdate, nsnucl, ocnv_init, ocsv_init, & 787 read_restart_data_salsa, reglim, salsa, salsa_emission_mode, & 788 sigmag, skip_time_do_salsa, surface_aerosol_flux, & 789 van_der_waals_coagc, write_binary_salsa 790 701 791 line = ' ' 702 703 792 ! 704 793 ! Try to find salsa package … … 709 798 ENDDO 710 799 BACKSPACE ( 11 ) 711 712 800 ! 713 801 ! Read userdefined namelist 714 802 READ ( 11, salsa_parameters ) 715 716 803 ! 717 804 ! Enable salsa (salsa switch in modules.f90) … … 719 806 720 807 10 CONTINUE 721 808 722 809 END SUBROUTINE salsa_parin 723 810 724 725 811 !! 726 812 ! Description: … … 730 816 SUBROUTINE salsa_check_parameters 731 817 732 USE control_parameters, &818 USE control_parameters, & 733 819 ONLY: message_string 734 820 735 821 IMPLICIT NONE 736 822 737 823 ! 738 824 ! Checks go here (cf. check_parameters.f90). 739 825 IF ( salsa .AND. .NOT. humidity ) THEN 740 WRITE( message_string, * ) 'salsa = ', salsa, ' is ', & 741 'not allowed with humidity = ', humidity 742 CALL message( 'check_parameters', 'SA0009', 1, 2, 0, 6, 0 ) 826 WRITE( message_string, * ) 'salsa = ', salsa, ' is not allowed with humidity = ', humidity 827 CALL message( 'salsa_check_parameters', 'PA0594', 1, 2, 0, 6, 0 ) 743 828 ENDIF 744 829 745 830 IF ( bc_salsa_b == 'dirichlet' ) THEN 746 831 ibc_salsa_b = 0 … … 748 833 ibc_salsa_b = 1 749 834 ELSE 750 message_string = 'unknown boundary condition: bc_salsa_b = "' & 751 // TRIM( bc_salsa_t ) // '"' 752 CALL message( 'check_parameters', 'SA0011', 1, 2, 0, 6, 0 ) 835 message_string = 'unknown boundary condition: bc_salsa_b = "' // TRIM( bc_salsa_t ) // '"' 836 CALL message( 'salsa_check_parameters', 'PA0595', 1, 2, 0, 6, 0 ) 753 837 ENDIF 754 838 755 839 IF ( bc_salsa_t == 'dirichlet' ) THEN 756 840 ibc_salsa_t = 0 757 841 ELSEIF ( bc_salsa_t == 'neumann' ) THEN 758 842 ibc_salsa_t = 1 843 ELSEIF ( bc_salsa_t == 'nested' ) THEN 844 ibc_salsa_t = 2 759 845 ELSE 760 message_string = 'unknown boundary condition: bc_salsa_t = "' & 761 // TRIM( bc_salsa_t ) // '"' 762 CALL message( 'check_parameters', 'SA0012', 1, 2, 0, 6, 0 ) 846 message_string = 'unknown boundary condition: bc_salsa_t = "' // TRIM( bc_salsa_t ) // '"' 847 CALL message( 'salsa_check_parameters', 'PA0596', 1, 2, 0, 6, 0 ) 763 848 ENDIF 764 849 765 850 IF ( nj3 < 1 .OR. nj3 > 3 ) THEN 766 851 message_string = 'unknown nj3 (must be 13)' 767 CALL message( ' check_parameters', 'SA0044', 1, 2, 0, 6, 0 )852 CALL message( 'salsa_check_parameters', 'PA0597', 1, 2, 0, 6, 0 ) 768 853 ENDIF 769 854 855 IF ( salsa_emission_mode == 'read_from_file' .AND. ibc_salsa_b == 0 ) THEN 856 message_string = 'salsa_emission_mode == read_from_file requires bc_salsa_b = "Neumann"' 857 CALL message( 'salsa_check_parameters','PA0598', 1, 2, 0, 6, 0 ) 858 ENDIF 859 770 860 END SUBROUTINE salsa_check_parameters 771 861 … … 775 865 !  776 866 !> Subroutine defining appropriate grid for netcdf variables. 777 !> It is called out from subroutine netcdf. 867 !> It is called out from subroutine netcdf. 778 868 !> Same grid as for other scalars (see netcdf_interface_mod.f90) 779 869 !! 780 870 SUBROUTINE salsa_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) 781 871 782 872 IMPLICIT NONE 783 873 784 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<785 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<786 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<787 CHARACTER 788 874 CHARACTER(LEN=*), INTENT(OUT) :: grid_x !< 875 CHARACTER(LEN=*), INTENT(OUT) :: grid_y !< 876 CHARACTER(LEN=*), INTENT(OUT) :: grid_z !< 877 CHARACTER(LEN=*), INTENT(IN) :: var !< 878 789 879 LOGICAL, INTENT(OUT) :: found !< 790 880 791 881 found = .TRUE. 792 882 ! … … 794 884 795 885 IF ( var(1:2) == 'g_' ) THEN 796 grid_x = 'x' 797 grid_y = 'y' 798 grid_z = 'zu' 886 grid_x = 'x' 887 grid_y = 'y' 888 grid_z = 'zu' 799 889 ELSEIF ( var(1:4) == 'LDSA' ) THEN 800 grid_x = 'x' 801 grid_y = 'y' 890 grid_x = 'x' 891 grid_y = 'y' 802 892 grid_z = 'zu' 803 893 ELSEIF ( var(1:5) == 'm_bin' ) THEN 804 grid_x = 'x' 805 grid_y = 'y' 894 grid_x = 'x' 895 grid_y = 'y' 806 896 grid_z = 'zu' 807 897 ELSEIF ( var(1:5) == 'N_bin' ) THEN 808 grid_x = 'x' 809 grid_y = 'y' 898 grid_x = 'x' 899 grid_y = 'y' 810 900 grid_z = 'zu' 811 901 ELSEIF ( var(1:4) == 'Ntot' ) THEN 812 grid_x = 'x' 813 grid_y = 'y' 902 grid_x = 'x' 903 grid_y = 'y' 814 904 grid_z = 'zu' 815 905 ELSEIF ( var(1:2) == 'PM' ) THEN 816 grid_x = 'x' 817 grid_y = 'y' 906 grid_x = 'x' 907 grid_y = 'y' 818 908 grid_z = 'zu' 819 909 ELSEIF ( var(1:2) == 's_' ) THEN 820 grid_x = 'x' 821 grid_y = 'y' 910 grid_x = 'x' 911 grid_y = 'y' 822 912 grid_z = 'zu' 823 913 ELSE … … 830 920 END SUBROUTINE salsa_define_netcdf_grid 831 921 832 833 922 !! 834 923 ! Description: … … 846 935 WRITE( io, 2 ) skip_time_do_salsa 847 936 WRITE( io, 3 ) dt_salsa 848 WRITE( io, 12 ) SHAPE( aerosol_number(1)%conc ), nbins937 WRITE( io, 4 ) SHAPE( aerosol_number(1)%conc ), nbins_aerosol 849 938 IF ( advect_particle_water ) THEN 850 WRITE( io, 16 ) SHAPE( aerosol_mass(1)%conc ), ncc_tot*nbins,&939 WRITE( io, 5 ) SHAPE( aerosol_mass(1)%conc ), ncomponents_mass*nbins_aerosol, & 851 940 advect_particle_water 852 941 ELSE 853 WRITE( io, 16 ) SHAPE( aerosol_mass(1)%conc ), ncc*nbins, & 854 advect_particle_water 942 WRITE( io, 5 ) SHAPE( aerosol_mass(1)%conc ), ncc*nbins_aerosol, advect_particle_water 855 943 ENDIF 856 944 IF ( .NOT. salsa_gases_from_chem ) THEN 857 WRITE( io, 17 ) SHAPE( aerosol_mass(1)%conc ), ngast, & 858 salsa_gases_from_chem 945 WRITE( io, 6 ) SHAPE( aerosol_mass(1)%conc ), ngases_salsa, salsa_gases_from_chem 859 946 ENDIF 860 WRITE( io, 4 )947 WRITE( io, 7 ) 861 948 IF ( nsnucl > 0 ) THEN 862 WRITE( io, 5) nsnucl, nj3949 WRITE( io, 8 ) nsnucl, nj3 863 950 ENDIF 864 951 IF ( nlcoag ) THEN 865 WRITE( io, 6 )952 WRITE( io, 9 ) 866 953 ENDIF 867 954 IF ( nlcnd ) THEN 868 WRITE( io, 7 ) nlcndgas, nlcndh2oae 955 WRITE( io, 10 ) nlcndgas, nlcndh2oae 956 ENDIF 957 IF ( lspartition ) THEN 958 WRITE( io, 11 ) 869 959 ENDIF 870 960 IF ( nldepo ) THEN 871 WRITE( io, 1 4 ) nldepo_vege, nldepo_topo961 WRITE( io, 12 ) nldepo_pcm, nldepo_surf 872 962 ENDIF 873 WRITE( io, 8) reglim, nbin, bin_low_limits874 WRITE( io, 15) nsect875 WRITE( io, 1 3) ncc, listspec, mass_fracs_a, mass_fracs_b963 WRITE( io, 13 ) reglim, nbin, bin_low_limits 964 IF ( isdtyp == 0 ) WRITE( io, 14 ) nsect 965 WRITE( io, 15 ) ncc, listspec, mass_fracs_a, mass_fracs_b 876 966 IF ( .NOT. salsa_gases_from_chem ) THEN 877 WRITE( io, 18 ) ngast, H2SO4_init, HNO3_init, NH3_init, OCNV_init, & 878 OCSV_init 967 WRITE( io, 16 ) ngases_salsa, h2so4_init, hno3_init, nh3_init, ocnv_init, ocsv_init 879 968 ENDIF 880 WRITE( io, 9) isdtyp, igctyp969 WRITE( io, 17 ) isdtyp, igctyp 881 970 IF ( isdtyp == 0 ) THEN 882 WRITE( io, 1 0) dpg, sigmag, n_lognorm971 WRITE( io, 18 ) dpg, sigmag, n_lognorm 883 972 ELSE 884 WRITE( io, 1 1)973 WRITE( io, 19 ) 885 974 ENDIF 886 887 888 1 FORMAT (//' SALSA information:'/ & 975 IF ( nest_salsa ) WRITE( io, 20 ) nest_salsa 976 WRITE( io, 21 ) salsa_emission_mode 977 978 979 1 FORMAT (//' SALSA information:'/ & 889 980 ' '/) 890 981 2 FORMAT (' Starts at: skip_time_do_salsa = ', F10.2, ' s') 891 982 3 FORMAT (/' Timestep: dt_salsa = ', F6.2, ' s') 892 12 FORMAT (/' Array shape (z,y,x,bins):'/&983 4 FORMAT (/' Array shape (z,y,x,bins):'/ & 893 984 ' aerosol_number: ', 4(I3)) 894 16 FORMAT (/' aerosol_mass: ', 4(I3),/&985 5 FORMAT (/' aerosol_mass: ', 4(I3),/ & 895 986 ' (advect_particle_water = ', L1, ')') 896 17 FORMAT (' salsa_gas: ', 4(I3),/&987 6 FORMAT (' salsa_gas: ', 4(I3),/ & 897 988 ' (salsa_gases_from_chem = ', L1, ')') 898 4 FORMAT (/' Aerosol dynamic processes included: ') 899 5 FORMAT (/' nucleation (scheme = ', I1, ' and J3 parametrization = ',& 900 I1, ')') 901 6 FORMAT (/' coagulation') 902 7 FORMAT (/' condensation (of precursor gases = ', L1, & 903 ' and water vapour = ', L1, ')' ) 904 14 FORMAT (/' dry deposition (on vegetation = ', L1, & 905 ' and on topography = ', L1, ')') 906 8 FORMAT (/' Aerosol bin subrange limits (in metres): ', 3(ES10.2E3), / & 907 ' Number of size bins for each aerosol subrange: ', 2I3,/ & 989 7 FORMAT (/' Aerosol dynamic processes included: ') 990 8 FORMAT (/' nucleation (scheme = ', I1, ' and J3 parametrization = ', I1, ')') 991 9 FORMAT (/' coagulation') 992 10 FORMAT (/' condensation (of precursor gases = ', L1, ' and water vapour = ', L1, ')' ) 993 11 FORMAT (/' dissolutional growth by HNO3 and NH3') 994 12 FORMAT (/' dry deposition (on vegetation = ', L1, ' and on topography = ', L1, ')') 995 13 FORMAT (/' Aerosol bin subrange limits (in metres): ', 3(ES10.2E3), / & 996 ' Number of size bins for each aerosol subrange: ', 2I3,/ & 908 997 ' Aerosol bin limits (in metres): ', 9(ES10.2E3)) 909 15 FORMAT (' Initial number concentration in bins at the lowest level', & 910 ' (#/m**3):', 9(ES10.2E3)) 911 13 FORMAT (/' Number of chemical components used: ', I1,/ & 912 ' Species: ',7(A6),/ & 913 ' Initial relative contribution of each species to particle', & 914 ' volume in:',/ & 915 ' abins: ', 7(F6.3),/ & 998 14 FORMAT (' Initial number concentration in bins at the lowest level (#/m**3):', 9(ES10.2E3)) 999 15 FORMAT (/' Number of chemical components used: ', I1,/ & 1000 ' Species: ',7(A6),/ & 1001 ' Initial relative contribution of each species to particle volume in:',/ & 1002 ' abins: ', 7(F6.3),/ & 916 1003 ' bbins: ', 7(F6.3)) 917 1 8 FORMAT (/' Number of gaseous tracers used: ', I1,/&918 ' Initial gas concentrations:',/ &919 ' H2SO4: ',ES12.4E3, ' #/m**3',/ &920 ' HNO3: ',ES12.4E3, ' #/m**3',/ &921 ' NH3: ',ES12.4E3, ' #/m**3',/ &922 ' OCNV: ',ES12.4E3, ' #/m**3',/ &1004 16 FORMAT (/' Number of gaseous tracers used: ', I1,/ & 1005 ' Initial gas concentrations:',/ & 1006 ' H2SO4: ',ES12.4E3, ' #/m**3',/ & 1007 ' HNO3: ',ES12.4E3, ' #/m**3',/ & 1008 ' NH3: ',ES12.4E3, ' #/m**3',/ & 1009 ' OCNV: ',ES12.4E3, ' #/m**3',/ & 923 1010 ' OCSV: ',ES12.4E3, ' #/m**3') 924 9 FORMAT (/' Initialising concentrations: ', /&925 ' Aerosol size distribution: isdtyp = ', I1,/ &1011 17 FORMAT (/' Initialising concentrations: ', / & 1012 ' Aerosol size distribution: isdtyp = ', I1,/ & 926 1013 ' Gas concentrations: igctyp = ', I1 ) 927 10 FORMAT ( ' Mode diametres: dpg(nmod) = ', 7(F7.3),/ & 928 ' Standard deviation: sigmag(nmod) = ', 7(F7.2),/ & 929 ' Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3) ) 930 11 FORMAT (/' Size distribution read from a file.') 1014 18 FORMAT ( ' Mode diametres: dpg(nmod) = ', 7(F7.3), ' (m)', / & 1015 ' Standard deviation: sigmag(nmod) = ', 7(F7.2),/ & 1016 ' Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3), ' (#/m3)' ) 1017 19 FORMAT (/' Size distribution read from a file.') 1018 20 FORMAT (/' Nesting for salsa variables: ', L1 ) 1019 21 FORMAT (/' Emissions: salsa_emission_mode = ', A ) 931 1020 932 1021 END SUBROUTINE salsa_header … … 938 1027 !! 939 1028 SUBROUTINE salsa_init_arrays 940 941 USE surface_mod, & 942 ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 943 surf_usm_v 1029 1030 USE chem_gasphase_mod, & 1031 ONLY: nvar 1032 1033 USE surface_mod, & 1034 ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v 944 1035 945 1036 IMPLICIT NONE 946 947 INTEGER(iwp) :: gases_available !< Number of available gas components in 948 !< the chemistry model 949 INTEGER(iwp) :: i !< loop index for allocating 950 INTEGER(iwp) :: l !< loop index for allocating: surfaces 951 INTEGER(iwp) :: lsp !< loop index for chem species in the chemistry model 952 1037 1038 INTEGER(iwp) :: gases_available !< Number of available gas components in the chemistry model 1039 INTEGER(iwp) :: i !< loop index for allocating 1040 INTEGER(iwp) :: l !< loop index for allocating: surfaces 1041 INTEGER(iwp) :: lsp !< loop index for chem species in the chemistry model 1042 953 1043 gases_available = 0 954 955 1044 ! 956 1045 ! Allocate prognostic variables (see salsa_swap_timelevel) 957 958 1046 ! 959 1047 ! Set derived indices: 960 ! (This does the same as the subroutine salsa_initialize in SALSA/ 961 ! UCLALESSALSA) 962 in1a = 1 ! 1st index of subrange 1a 963 in2a = in1a + nbin(1) ! 1st index of subrange 2a 964 fn1a = in2a  1 ! last index of subrange 1a 965 fn2a = fn1a + nbin(2) ! last index of subrange 2a 966 967 ! 968 ! If the fraction of insoluble aerosols in subrange 2 is zero: do not allocate 969 ! arrays for them 1048 ! (This does the same as the subroutine salsa_initialize in SALSA/UCLALESSALSA) 1049 start_subrange_1a = 1 ! 1st index of subrange 1a 1050 start_subrange_2a = start_subrange_1a + nbin(1) ! 1st index of subrange 2a 1051 end_subrange_1a = start_subrange_2a  1 ! last index of subrange 1a 1052 end_subrange_2a = end_subrange_1a + nbin(2) ! last index of subrange 2a 1053 1054 ! 1055 ! If the fraction of insoluble aerosols in subrange 2 is zero: do not allocate arrays for them 970 1056 IF ( nf2a > 0.999999_wp .AND. SUM( mass_fracs_b ) < 0.00001_wp ) THEN 971 1057 no_insoluble = .TRUE. 972 in2b = fn2a+1! 1st index of subrange 2b973 fn2b = fn2a! last index of subrange 2b1058 start_subrange_2b = end_subrange_2a+1 ! 1st index of subrange 2b 1059 end_subrange_2b = end_subrange_2a ! last index of subrange 2b 974 1060 ELSE 975 in2b = in2a + nbin(2)! 1st index of subrange 2b976 fn2b = fn2a + nbin(2)! last index of subrange 2b1061 start_subrange_2b = start_subrange_2a + nbin(2) ! 1st index of subrange 2b 1062 end_subrange_2b = end_subrange_2a + nbin(2) ! last index of subrange 2b 977 1063 ENDIF 978 979 980 nbins = fn2b ! total number of aerosol size bins 981 ! 1064 1065 nbins_aerosol = end_subrange_2b ! total number of aerosol size bins 1066 ! 982 1067 ! Create index tables for different aerosol components 983 1068 CALL component_index_constructor( prtcl, ncc, maxspec, listspec ) 984 985 nc c_tot= ncc986 IF ( advect_particle_water ) nc c_tot= ncc + 1 ! Add water987 1069 1070 ncomponents_mass = ncc 1071 IF ( advect_particle_water ) ncomponents_mass = ncc + 1 ! Add water 1072 988 1073 ! 989 1074 ! Allocate: 990 ALLOCATE( aero(nbins), bin_low_limits(nbins), nsect(nbins), massacc(nbins) ) 991 IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins) ) 992 ALLOCATE( Ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins) ) 993 994 ! 1075 ALLOCATE( aero(nbins_aerosol), bc_am_t_val(nbins_aerosol*ncomponents_mass), & 1076 bc_an_t_val(ngases_salsa), bc_gt_t_val(nbins_aerosol), bin_low_limits(nbins_aerosol),& 1077 nsect(nbins_aerosol), massacc(nbins_aerosol) ) 1078 ALLOCATE( k_topo_top(nysg:nyng,nxlg:nxrg) ) 1079 IF ( nldepo ) ALLOCATE( sedim_vd(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1080 ALLOCATE( ra_dry(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1081 1082 ! 995 1083 ! Aerosol number concentration 996 ALLOCATE( aerosol_number(nbins ) )997 ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins ),&998 nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins ),&999 nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins ) )1084 ALLOCATE( aerosol_number(nbins_aerosol) ) 1085 ALLOCATE( nconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol), & 1086 nconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol), & 1087 nconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nbins_aerosol) ) 1000 1088 nconc_1 = 0.0_wp 1001 1089 nconc_2 = 0.0_wp 1002 1090 nconc_3 = 0.0_wp 1003 1004 DO i = 1, nbins 1091 1092 DO i = 1, nbins_aerosol 1005 1093 aerosol_number(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,i) 1006 1094 aerosol_number(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,i) … … 1012 1100 aerosol_number(i)%init(nzb:nzt+1), & 1013 1101 aerosol_number(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task1) ) 1014 ENDDO 1015 1016 ! 1017 ! Aerosol mass concentration 1018 ALLOCATE( aerosol_mass(nc c_tot*nbins) )1019 ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nc c_tot*nbins),&1020 mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nc c_tot*nbins),&1021 mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nc c_tot*nbins) )1102 ENDDO 1103 1104 ! 1105 ! Aerosol mass concentration 1106 ALLOCATE( aerosol_mass(ncomponents_mass*nbins_aerosol) ) 1107 ALLOCATE( mconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol), & 1108 mconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol), & 1109 mconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ncomponents_mass*nbins_aerosol) ) 1022 1110 mconc_1 = 0.0_wp 1023 1111 mconc_2 = 0.0_wp 1024 1112 mconc_3 = 0.0_wp 1025 1026 DO i = 1, nc c_tot*nbins1113 1114 DO i = 1, ncomponents_mass*nbins_aerosol 1027 1115 aerosol_mass(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,i) 1028 1116 aerosol_mass(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,i) 1029 aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i) 1030 ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task1), &1031 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task1), &1032 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task1), &1033 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task1), &1034 aerosol_mass(i)%init(nzb:nzt+1), &1117 aerosol_mass(i)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_3(:,:,:,i) 1118 ALLOCATE( aerosol_mass(i)%flux_s(nzb+1:nzt,0:threads_per_task1), & 1119 aerosol_mass(i)%diss_s(nzb+1:nzt,0:threads_per_task1), & 1120 aerosol_mass(i)%flux_l(nzb+1:nzt,nys:nyn,0:threads_per_task1), & 1121 aerosol_mass(i)%diss_l(nzb+1:nzt,nys:nyn,0:threads_per_task1), & 1122 aerosol_mass(i)%init(nzb:nzt+1), & 1035 1123 aerosol_mass(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task1) ) 1036 1124 ENDDO 1037 1038 ! 1039 ! Surface fluxes: answs = aerosol number, amsws = aerosol mass 1125 1126 ! 1127 ! Surface fluxes: answs = aerosol number, amsws = aerosol mass 1040 1128 ! 1041 1129 ! Horizontal surfaces: default type 1042 1130 DO l = 0, 2 ! upward (l=0), downward (l=1) and model top (l=2) 1043 ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins ) )1044 ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins *ncc_tot) )1131 ALLOCATE( surf_def_h(l)%answs( 1:surf_def_h(l)%ns, nbins_aerosol ) ) 1132 ALLOCATE( surf_def_h(l)%amsws( 1:surf_def_h(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1045 1133 surf_def_h(l)%answs = 0.0_wp 1046 1134 surf_def_h(l)%amsws = 0.0_wp 1047 1135 ENDDO 1048 ! Horizontal surfaces: natural type 1049 IF ( land_surface ) THEN 1050 ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins ) ) 1051 ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins*ncc_tot ) ) 1052 surf_lsm_h%answs = 0.0_wp 1053 surf_lsm_h%amsws = 0.0_wp 1054 ENDIF 1055 ! Horizontal surfaces: urban type 1056 IF ( urban_surface ) THEN 1057 ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins ) ) 1058 ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins*ncc_tot ) ) 1059 surf_usm_h%answs = 0.0_wp 1060 surf_usm_h%amsws = 0.0_wp 1061 ENDIF 1062 ! 1063 ! Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and 1064 ! westward (l=3) facing 1065 DO l = 0, 3 1066 ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins ) ) 1136 ! 1137 ! Horizontal surfaces: natural type 1138 ALLOCATE( surf_lsm_h%answs( 1:surf_lsm_h%ns, nbins_aerosol ) ) 1139 ALLOCATE( surf_lsm_h%amsws( 1:surf_lsm_h%ns, nbins_aerosol*ncomponents_mass ) ) 1140 surf_lsm_h%answs = 0.0_wp 1141 surf_lsm_h%amsws = 0.0_wp 1142 ! 1143 ! Horizontal surfaces: urban type 1144 ALLOCATE( surf_usm_h%answs( 1:surf_usm_h%ns, nbins_aerosol ) ) 1145 ALLOCATE( surf_usm_h%amsws( 1:surf_usm_h%ns, nbins_aerosol*ncomponents_mass ) ) 1146 surf_usm_h%answs = 0.0_wp 1147 surf_usm_h%amsws = 0.0_wp 1148 1149 ! 1150 ! Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and westward (l=3) facing 1151 DO l = 0, 3 1152 ALLOCATE( surf_def_v(l)%answs( 1:surf_def_v(l)%ns, nbins_aerosol ) ) 1067 1153 surf_def_v(l)%answs = 0.0_wp 1068 ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins *ncc_tot) )1154 ALLOCATE( surf_def_v(l)%amsws( 1:surf_def_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1069 1155 surf_def_v(l)%amsws = 0.0_wp 1070 1071 IF ( land_surface) THEN 1072 ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins ) ) 1073 surf_lsm_v(l)%answs = 0.0_wp 1074 ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins*ncc_tot ) ) 1075 surf_lsm_v(l)%amsws = 0.0_wp 1076 ENDIF 1077 1078 IF ( urban_surface ) THEN 1079 ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins ) ) 1080 surf_usm_v(l)%answs = 0.0_wp 1081 ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins*ncc_tot ) ) 1082 surf_usm_v(l)%amsws = 0.0_wp 1083 ENDIF 1084 ENDDO 1085 1156 1157 ALLOCATE( surf_lsm_v(l)%answs( 1:surf_lsm_v(l)%ns, nbins_aerosol ) ) 1158 surf_lsm_v(l)%answs = 0.0_wp 1159 ALLOCATE( surf_lsm_v(l)%amsws( 1:surf_lsm_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1160 surf_lsm_v(l)%amsws = 0.0_wp 1161 1162 ALLOCATE( surf_usm_v(l)%answs( 1:surf_usm_v(l)%ns, nbins_aerosol ) ) 1163 surf_usm_v(l)%answs = 0.0_wp 1164 ALLOCATE( surf_usm_v(l)%amsws( 1:surf_usm_v(l)%ns, nbins_aerosol*ncomponents_mass ) ) 1165 surf_usm_v(l)%amsws = 0.0_wp 1166 1167 ENDDO 1168 1086 1169 ! 1087 1170 ! Concentration of gaseous tracers (1. SO4, 2. HNO3, 3. NH3, 4. OCNV, 5. OCSV) … … 1091 1174 ! allocate salsa_gas array. 1092 1175 1093 IF ( air_chemistry ) THEN 1176 IF ( air_chemistry ) THEN 1094 1177 DO lsp = 1, nvar 1095 IF ( TRIM( chem_species(lsp)%name ) == 'H2SO4' ) THEN 1096 gases_available = gases_available + 1 1097 gas_index_chem(1) = lsp 1098 ELSEIF ( TRIM( chem_species(lsp)%name ) == 'HNO3' ) THEN 1099 gases_available = gases_available + 1 1100 gas_index_chem(2) = lsp 1101 ELSEIF ( TRIM( chem_species(lsp)%name ) == 'NH3' ) THEN 1102 gases_available = gases_available + 1 1103 gas_index_chem(3) = lsp 1104 ELSEIF ( TRIM( chem_species(lsp)%name ) == 'OCNV' ) THEN 1105 gases_available = gases_available + 1 1106 gas_index_chem(4) = lsp 1107 ELSEIF ( TRIM( chem_species(lsp)%name ) == 'OCSV' ) THEN 1108 gases_available = gases_available + 1 1109 gas_index_chem(5) = lsp 1110 ENDIF 1178 SELECT CASE ( TRIM( chem_species(lsp)%name ) ) 1179 CASE ( 'H2SO4', 'h2so4' ) 1180 gases_available = gases_available + 1 1181 gas_index_chem(1) = lsp 1182 CASE ( 'HNO3', 'hno3' ) 1183 gases_available = gases_available + 1 1184 gas_index_chem(2) = lsp 1185 CASE ( 'NH3', 'nh3' ) 1186 gases_available = gases_available + 1 1187 gas_index_chem(3) = lsp 1188 CASE ( 'OCNV', 'ocnv' ) 1189 gases_available = gases_available + 1 1190 gas_index_chem(4) = lsp 1191 CASE ( 'OCSV', 'ocsv' ) 1192 gases_available = gases_available + 1 1193 gas_index_chem(5) = lsp 1194 END SELECT 1111 1195 ENDDO 1112 1196 1113 IF ( gases_available == ngas t) THEN1197 IF ( gases_available == ngases_salsa ) THEN 1114 1198 salsa_gases_from_chem = .TRUE. 1115 1199 ELSE 1116 WRITE( message_string, * ) 'SALSA is run together with chemistry '// & 1117 'but not all gaseous components are '// & 1118 'provided by kpp (H2SO4, HNO3, NH3, '// & 1119 'OCNV, OCSC)' 1120 CALL message( 'check_parameters', 'SA0024', 1, 2, 0, 6, 0 ) 1200 WRITE( message_string, * ) 'SALSA is run together with chemistry but not all gaseous '// & 1201 'components are provided by kpp (H2SO4, HNO3, NH3, OCNV, OCSV)' 1202 CALL message( 'check_parameters', 'PA0599', 1, 2, 0, 6, 0 ) 1121 1203 ENDIF 1122 1204 1123 1205 ELSE 1124 1206 1125 ALLOCATE( salsa_gas(ngas t) )1126 ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngas t), &1127 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngas t), &1128 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngas t) )1207 ALLOCATE( salsa_gas(ngases_salsa) ) 1208 ALLOCATE( gconc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa), & 1209 gconc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa), & 1210 gconc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg,ngases_salsa) ) 1129 1211 gconc_1 = 0.0_wp 1130 1212 gconc_2 = 0.0_wp 1131 1213 gconc_3 = 0.0_wp 1132 1133 DO i = 1, ngas t1214 1215 DO i = 1, ngases_salsa 1134 1216 salsa_gas(i)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,i) 1135 1217 salsa_gas(i)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,i) … … 1141 1223 salsa_gas(i)%init(nzb:nzt+1), & 1142 1224 salsa_gas(i)%sums_ws_l(nzb:nzt+1,0:threads_per_task1) ) 1143 ENDDO 1225 ENDDO 1144 1226 ! 1145 1227 ! Surface fluxes: gtsws = gaseous tracer flux … … 1147 1229 ! Horizontal surfaces: default type 1148 1230 DO l = 0, 2 ! upward (l=0), downward (l=1) and model top (l=2) 1149 ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngas t) )1231 ALLOCATE( surf_def_h(l)%gtsws( 1:surf_def_h(l)%ns, ngases_salsa ) ) 1150 1232 surf_def_h(l)%gtsws = 0.0_wp 1151 1233 ENDDO 1152 ! Horizontal surfaces: natural type 1153 IF ( land_surface ) THEN 1154 ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngast ) ) 1155 surf_lsm_h%gtsws = 0.0_wp 1156 ENDIF 1157 ! Horizontal surfaces: urban type 1158 IF ( urban_surface ) THEN 1159 ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngast ) ) 1160 surf_usm_h%gtsws = 0.0_wp 1161 ENDIF 1234 ! Horizontal surfaces: natural type 1235 ALLOCATE( surf_lsm_h%gtsws( 1:surf_lsm_h%ns, ngases_salsa ) ) 1236 surf_lsm_h%gtsws = 0.0_wp 1237 ! Horizontal surfaces: urban type 1238 ALLOCATE( surf_usm_h%gtsws( 1:surf_usm_h%ns, ngases_salsa ) ) 1239 surf_usm_h%gtsws = 0.0_wp 1162 1240 ! 1163 1241 ! Vertical surfaces: northward (l=0), southward (l=1), eastward (l=2) and 1164 1242 ! westward (l=3) facing 1165 DO l = 0, 3 1166 ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngas t) )1243 DO l = 0, 3 1244 ALLOCATE( surf_def_v(l)%gtsws( 1:surf_def_v(l)%ns, ngases_salsa ) ) 1167 1245 surf_def_v(l)%gtsws = 0.0_wp 1168 IF ( land_surface ) THEN 1169 ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngast ) ) 1170 surf_lsm_v(l)%gtsws = 0.0_wp 1171 ENDIF 1172 IF ( urban_surface ) THEN 1173 ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngast ) ) 1174 surf_usm_v(l)%gtsws = 0.0_wp 1175 ENDIF 1246 ALLOCATE( surf_lsm_v(l)%gtsws( 1:surf_lsm_v(l)%ns, ngases_salsa ) ) 1247 surf_lsm_v(l)%gtsws = 0.0_wp 1248 ALLOCATE( surf_usm_v(l)%gtsws( 1:surf_usm_v(l)%ns, ngases_salsa ) ) 1249 surf_usm_v(l)%gtsws = 0.0_wp 1176 1250 ENDDO 1177 1251 ENDIF 1178 1252 1179 1253 END SUBROUTINE salsa_init_arrays 1180 1254 … … 1182 1256 ! Description: 1183 1257 !  1184 !> Initialization of SALSA. Based on salsa_initialize in UCLALESSALSA. 1258 !> Initialization of SALSA. Based on salsa_initialize in UCLALESSALSA. 1185 1259 !> Subroutines salsa_initialize, SALSAinit and DiagInitAero in UCLALESSALSA are 1186 1260 !> also merged here. … … 1189 1263 1190 1264 IMPLICIT NONE 1191 1192 INTEGER(iwp) :: b 1193 INTEGER(iwp) :: c 1194 INTEGER(iwp) :: g 1195 INTEGER(iwp) :: i 1196 INTEGER(iwp) :: j 1197 1198 CALL location_message( 'initializing SALSA model', .TRUE. ) 1199 1265 1266 INTEGER(iwp) :: i !< 1267 INTEGER(iwp) :: ib !< loop index for aerosol number bins 1268 INTEGER(iwp) :: ic !< loop index for aerosol mass bins 1269 INTEGER(iwp) :: ig !< loop index for gases 1270 INTEGER(iwp) :: ii !< index for indexing 1271 INTEGER(iwp) :: j !< 1272 1273 CALL location_message( 'initializing salsa (sectional aerosol module )', .TRUE. ) 1274 1200 1275 bin_low_limits = 0.0_wp 1276 k_topo_top = 0 1201 1277 nsect = 0.0_wp 1202 massacc = 1.0_wp 1203 1278 massacc = 1.0_wp 1279 1204 1280 ! 1205 1281 ! Indices for chemical components used (1 = not used) 1206 i = 01282 ii = 0 1207 1283 IF ( is_used( prtcl, 'SO4' ) ) THEN 1208 i so4 = get_index( prtcl,'SO4' )1209 i =i + 11284 index_so4 = get_index( prtcl,'SO4' ) 1285 ii = ii + 1 1210 1286 ENDIF 1211 1287 IF ( is_used( prtcl,'OC' ) ) THEN 1212 i oc = get_index(prtcl, 'OC')1213 i =i + 11288 index_oc = get_index(prtcl, 'OC') 1289 ii = ii + 1 1214 1290 ENDIF 1215 1291 IF ( is_used( prtcl, 'BC' ) ) THEN 1216 i bc = get_index( prtcl, 'BC' )1217 i =i + 11292 index_bc = get_index( prtcl, 'BC' ) 1293 ii = ii + 1 1218 1294 ENDIF 1219 1295 IF ( is_used( prtcl, 'DU' ) ) THEN 1220 i du = get_index( prtcl, 'DU' )1221 i =i + 11296 index_du = get_index( prtcl, 'DU' ) 1297 ii = ii + 1 1222 1298 ENDIF 1223 1299 IF ( is_used( prtcl, 'SS' ) ) THEN 1224 i ss = get_index( prtcl, 'SS' )1225 i =i + 11300 index_ss = get_index( prtcl, 'SS' ) 1301 ii = ii + 1 1226 1302 ENDIF 1227 1303 IF ( is_used( prtcl, 'NO' ) ) THEN 1228 in o = get_index( prtcl, 'NO' )1229 i =i + 11304 index_no = get_index( prtcl, 'NO' ) 1305 ii = ii + 1 1230 1306 ENDIF 1231 1307 IF ( is_used( prtcl, 'NH' ) ) THEN 1232 in h = get_index( prtcl, 'NH' )1233 i =i + 11308 index_nh = get_index( prtcl, 'NH' ) 1309 ii = ii + 1 1234 1310 ENDIF 1235 ! 1311 ! 1236 1312 ! All species must be known 1237 IF ( i /= ncc ) THEN 1238 message_string = 'Unknown aerosol species/component(s) given in the' // & 1239 ' initialization' 1240 CALL message( 'salsa_mod: salsa_init', 'SA0020', 1, 2, 0, 6, 0 ) 1313 IF ( ii /= ncc ) THEN 1314 message_string = 'Unknown aerosol species/component(s) given in the initialization' 1315 CALL message( 'salsa_mod: salsa_init', 'PA0600', 1, 2, 0, 6, 0 ) 1241 1316 ENDIF 1242 1317 ! 1318 ! Partition and dissolutional growth by gaseous HNO3 and NH3 1319 IF ( index_no > 0 .AND. index_nh > 0 .AND. index_so4 > 0 ) lspartition = .TRUE. 1243 1320 ! 1244 1321 ! Initialise … … 1249 1326 aero(:)%numc = nclim 1250 1327 aero(:)%core = 1.0E10_wp 1251 DO c = 1, maxspec+1 ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O1252 aero(:)%volc( c) = 0.0_wp1328 DO ic = 1, maxspec+1 ! 1:SO4, 2:OC, 3:BC, 4:DU, 5:SS, 6:NO, 7:NH, 8:H2O 1329 aero(:)%volc(ic) = 0.0_wp 1253 1330 ENDDO 1254 1255 IF ( nldepo ) sedim_vd = 0.0_wp 1256 1257 DO b = 1, nbins1258 IF ( .NOT. read_restart_data_salsa ) aerosol_number( b)%conc = nclim1259 aerosol_number( b)%conc_p = 0.0_wp1260 aerosol_number( b)%tconc_m = 0.0_wp1261 aerosol_number( b)%flux_s = 0.0_wp1262 aerosol_number( b)%diss_s = 0.0_wp1263 aerosol_number( b)%flux_l = 0.0_wp1264 aerosol_number( b)%diss_l = 0.0_wp1265 aerosol_number( b)%init = nclim1266 aerosol_number( b)%sums_ws_l = 0.0_wp1331 1332 IF ( nldepo ) sedim_vd = 0.0_wp 1333 1334 DO ib = 1, nbins_aerosol 1335 IF ( .NOT. read_restart_data_salsa ) aerosol_number(ib)%conc = nclim 1336 aerosol_number(ib)%conc_p = 0.0_wp 1337 aerosol_number(ib)%tconc_m = 0.0_wp 1338 aerosol_number(ib)%flux_s = 0.0_wp 1339 aerosol_number(ib)%diss_s = 0.0_wp 1340 aerosol_number(ib)%flux_l = 0.0_wp 1341 aerosol_number(ib)%diss_l = 0.0_wp 1342 aerosol_number(ib)%init = nclim 1343 aerosol_number(ib)%sums_ws_l = 0.0_wp 1267 1344 ENDDO 1268 DO c = 1, ncc_tot*nbins1269 IF ( .NOT. read_restart_data_salsa ) aerosol_mass( c)%conc = mclim1270 aerosol_mass( c)%conc_p = 0.0_wp1271 aerosol_mass( c)%tconc_m = 0.0_wp1272 aerosol_mass( c)%flux_s = 0.0_wp1273 aerosol_mass( c)%diss_s = 0.0_wp1274 aerosol_mass( c)%flux_l = 0.0_wp1275 aerosol_mass( c)%diss_l = 0.0_wp1276 aerosol_mass( c)%init = mclim1277 aerosol_mass( c)%sums_ws_l = 0.0_wp1345 DO ic = 1, ncomponents_mass*nbins_aerosol 1346 IF ( .NOT. read_restart_data_salsa ) aerosol_mass(ic)%conc = mclim 1347 aerosol_mass(ic)%conc_p = 0.0_wp 1348 aerosol_mass(ic)%tconc_m = 0.0_wp 1349 aerosol_mass(ic)%flux_s = 0.0_wp 1350 aerosol_mass(ic)%diss_s = 0.0_wp 1351 aerosol_mass(ic)%flux_l = 0.0_wp 1352 aerosol_mass(ic)%diss_l = 0.0_wp 1353 aerosol_mass(ic)%init = mclim 1354 aerosol_mass(ic)%sums_ws_l = 0.0_wp 1278 1355 ENDDO 1279 1356 1280 1357 IF ( .NOT. salsa_gases_from_chem ) THEN 1281 DO g = 1, ngast1282 salsa_gas( g)%conc_p = 0.0_wp1283 salsa_gas( g)%tconc_m = 0.0_wp1284 salsa_gas( g)%flux_s = 0.0_wp1285 salsa_gas( g)%diss_s = 0.0_wp1286 salsa_gas( g)%flux_l = 0.0_wp1287 salsa_gas( g)%diss_l = 0.0_wp1288 salsa_gas( g)%sums_ws_l = 0.0_wp1358 DO ig = 1, ngases_salsa 1359 salsa_gas(ig)%conc_p = 0.0_wp 1360 salsa_gas(ig)%tconc_m = 0.0_wp 1361 salsa_gas(ig)%flux_s = 0.0_wp 1362 salsa_gas(ig)%diss_s = 0.0_wp 1363 salsa_gas(ig)%flux_l = 0.0_wp 1364 salsa_gas(ig)%diss_l = 0.0_wp 1365 salsa_gas(ig)%sums_ws_l = 0.0_wp 1289 1366 ENDDO 1290 1367 IF ( .NOT. read_restart_data_salsa ) THEN 1291 salsa_gas(1)%conc = H2SO4_init1292 salsa_gas(2)%conc = HNO3_init1293 salsa_gas(3)%conc = NH3_init1294 salsa_gas(4)%conc = OCNV_init1295 salsa_gas(5)%conc = OCSV_init1368 salsa_gas(1)%conc = h2so4_init 1369 salsa_gas(2)%conc = hno3_init 1370 salsa_gas(3)%conc = nh3_init 1371 salsa_gas(4)%conc = ocnv_init 1372 salsa_gas(5)%conc = ocsv_init 1296 1373 ENDIF 1297 1374 ! 1298 1375 ! Set initial value for gas compound tracers and initial values 1299 salsa_gas(1)%init = H2SO4_init1300 salsa_gas(2)%init = HNO3_init1301 salsa_gas(3)%init = NH3_init1302 salsa_gas(4)%init = OCNV_init1303 salsa_gas(5)%init = OCSV_init1376 salsa_gas(1)%init = h2so4_init 1377 salsa_gas(2)%init = hno3_init 1378 salsa_gas(3)%init = nh3_init 1379 salsa_gas(4)%init = ocnv_init 1380 salsa_gas(5)%init = ocsv_init 1304 1381 ENDIF 1305 1382 ! 1306 1383 ! Aerosol radius in each bin: dry and wet (m) 1307 Ra_dry = 1.0E10_wp1308 ! 1309 ! Initialise aerosol tracers 1384 ra_dry = 1.0E10_wp 1385 ! 1386 ! Initialise aerosol tracers 1310 1387 aero(:)%vhilim = 0.0_wp 1311 1388 aero(:)%vlolim = 0.0_wp … … 1315 1392 ! 1316 1393 ! Initialise the sectional particle size distribution 1317 CALL set_sizebins() 1318 ! 1319 ! Initialise locationdependent aerosol size distributions and 1320 ! chemical compositions: 1321 CALL aerosol_init 1322 ! 1323 ! Initalisation run of SALSA 1394 CALL set_sizebins 1395 ! 1396 ! Initialise locationdependent aerosol size distributions and chemical compositions: 1397 CALL aerosol_init 1398 ! 1399 ! Initalisation run of SALSA + calculate the vertical top index of the topography 1324 1400 DO i = nxl, nxr 1325 1401 DO j = nys, nyn 1402 1403 k_topo_top(j,i) = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,j,i), 12 ) ), DIM = 1 )  1 1404 1326 1405 CALL salsa_driver( i, j, 1 ) 1327 1406 CALL salsa_diagnostics( i, j ) 1328 1407 ENDDO 1329 ENDDO 1330 ! 1331 ! Set the aerosol and gas sources 1332 IF ( salsa_source_mode == 'read_from_file' ) THEN 1333 CALL salsa_set_source 1408 ENDDO 1409 ! 1410 ! Initialise the deposition scheme and surface types 1411 IF ( nldepo ) CALL init_deposition 1412 1413 IF ( salsa_emission_mode /= 'no_emission' ) THEN 1414 ! 1415 ! Read in and initialize emissions 1416 CALL salsa_emission_setup( .TRUE. ) 1417 IF ( .NOT. salsa_gases_from_chem .AND. salsa_emission_mode == 'read_from_file' ) THEN 1418 CALL salsa_gas_emission_setup( .TRUE. ) 1419 ENDIF 1334 1420 ENDIF 1335 1421 1336 1422 CALL location_message( 'finished', .TRUE. ) 1337 1423 1338 1424 END SUBROUTINE salsa_init 1339 1425 … … 1363 1449 !! 1364 1450 SUBROUTINE set_sizebins 1365 1451 1366 1452 IMPLICIT NONE 1367 ! 1368 ! Local variables 1369 INTEGER(iwp) :: cc1370 INTEGER(iwp) :: dd 1371 REAL(wp) :: ratio_d !< ratio of the upper and lower diameter of subranges1372 ! 1373 ! vlolim&vhilim: min & max *dry* volumes [fxm] 1453 1454 INTEGER(iwp) :: cc !< running index 1455 INTEGER(iwp) :: dd !< running index 1456 1457 REAL(wp) :: ratio_d !< ratio of the upper and lower diameter of subranges 1458 ! 1459 ! vlolim&vhilim: min & max *dry* volumes [fxm] 1374 1460 ! dmid: bin mid *dry* diameter (m) 1375 1461 ! vratiolo&vratiohi: volume ratio between the center and low/high limit … … 1377 1463 ! 1) Size subrange 1: 1378 1464 ratio_d = reglim(2) / reglim(1) ! section spacing (m) 1379 DO cc = in1a,fn1a 1380 aero(cc)%vlolim = api6 * ( reglim(1) * ratio_d ** & 1381 ( REAL( cc1 ) / nbin(1) ) ) ** 3.0_wp 1382 aero(cc)%vhilim = api6 * ( reglim(1) * ratio_d ** & 1383 ( REAL( cc ) / nbin(1) ) ) ** 3.0_wp 1384 aero(cc)%dmid = SQRT( ( aero(cc)%vhilim / api6 ) ** ( 1.0_wp / 3.0_wp ) & 1385 * ( aero(cc)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp ) ) 1386 aero(cc)%vratiohi = aero(cc)%vhilim / ( api6 * aero(cc)%dmid ** 3.0_wp ) 1387 aero(cc)%vratiolo = aero(cc)%vlolim / ( api6 * aero(cc)%dmid ** 3.0_wp ) 1465 DO cc = start_subrange_1a, end_subrange_1a 1466 aero(cc)%vlolim = api6 * ( reglim(1) * ratio_d**( REAL( cc1 ) / nbin(1) ) )**3 1467 aero(cc)%vhilim = api6 * ( reglim(1) * ratio_d**( REAL( cc ) / nbin(1) ) )**3 1468 aero(cc)%dmid = SQRT( ( aero(cc)%vhilim / api6 )**0.33333333_wp * & 1469 ( aero(cc)%vlolim / api6 )**0.33333333_wp ) 1470 aero(cc)%vratiohi = aero(cc)%vhilim / ( api6 * aero(cc)%dmid**3 ) 1471 aero(cc)%vratiolo = aero(cc)%vlolim / ( api6 * aero(cc)%dmid**3 ) 1388 1472 ENDDO 1389 1473 ! 1390 ! 2) Size subrange 2: 1474 ! 2) Size subrange 2: 1391 1475 ! 2.1) Subsubrange 2a: high hygroscopicity 1392 1476 ratio_d = reglim(3) / reglim(2) ! section spacing 1393 DO dd = in2a, fn2a 1394 cc = dd  in2a 1395 aero(dd)%vlolim = api6 * ( reglim(2) * ratio_d ** & 1396 ( REAL( cc ) / nbin(2) ) ) ** 3.0_wp 1397 aero(dd)%vhilim = api6 * ( reglim(2) * ratio_d ** & 1398 ( REAL( cc+1 ) / nbin(2) ) ) ** 3.0_wp 1399 aero(dd)%dmid = SQRT( ( aero(dd)%vhilim / api6 ) ** ( 1.0_wp / 3.0_wp ) & 1400 * ( aero(dd)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp ) ) 1401 aero(dd)%vratiohi = aero(dd)%vhilim / ( api6 * aero(dd)%dmid ** 3.0_wp ) 1402 aero(dd)%vratiolo = aero(dd)%vlolim / ( api6 * aero(dd)%dmid ** 3.0_wp ) 1477 DO dd = start_subrange_2a, end_subrange_2a 1478 cc = dd  start_subrange_2a 1479 aero(dd)%vlolim = api6 * ( reglim(2) * ratio_d**( REAL( cc ) / nbin(2) ) )**3 1480 aero(dd)%vhilim = api6 * ( reglim(2) * ratio_d**( REAL( cc+1 ) / nbin(2) ) )**3 1481 aero(dd)%dmid = SQRT( ( aero(dd)%vhilim / api6 )**0.33333333_wp * & 1482 ( aero(dd)%vlolim / api6 )**0.33333333_wp ) 1483 aero(dd)%vratiohi = aero(dd)%vhilim / ( api6 * aero(dd)%dmid**3 ) 1484 aero(dd)%vratiolo = aero(dd)%vlolim / ( api6 * aero(dd)%dmid**3 ) 1403 1485 ENDDO 1404 ! 1486 ! 1405 1487 ! 2.2) Subsubrange 2b: low hygroscopicity 1406 1488 IF ( .NOT. no_insoluble ) THEN 1407 aero( in2b:fn2b)%vlolim = aero(in2a:fn2a)%vlolim1408 aero( in2b:fn2b)%vhilim = aero(in2a:fn2a)%vhilim1409 aero( in2b:fn2b)%dmid = aero(in2a:fn2a)%dmid1410 aero( in2b:fn2b)%vratiohi = aero(in2a:fn2a)%vratiohi1411 aero( in2b:fn2b)%vratiolo = aero(in2a:fn2a)%vratiolo1489 aero(start_subrange_2b:end_subrange_2b)%vlolim = aero(start_subrange_2a:end_subrange_2a)%vlolim 1490 aero(start_subrange_2b:end_subrange_2b)%vhilim = aero(start_subrange_2a:end_subrange_2a)%vhilim 1491 aero(start_subrange_2b:end_subrange_2b)%dmid = aero(start_subrange_2a:end_subrange_2a)%dmid 1492 aero(start_subrange_2b:end_subrange_2b)%vratiohi = aero(start_subrange_2a:end_subrange_2a)%vratiohi 1493 aero(start_subrange_2b:end_subrange_2b)%vratiolo = aero(start_subrange_2a:end_subrange_2a)%vratiolo 1412 1494 ENDIF 1413 ! 1414 ! Initialize the wet diameter with the bin dry diameter to avoid numerical 1415 ! problems later 1495 ! 1496 ! Initialize the wet diameter with the bin dry diameter to avoid numerical problems later 1416 1497 aero(:)%dwet = aero(:)%dmid 1417 1498 ! 1418 ! Save bin limits (lower diameter) to be delivered to the host modelif needed1419 DO cc = 1, nbins 1420 bin_low_limits(cc) = ( aero(cc)%vlolim / api6 )** ( 1.0_wp / 3.0_wp )1421 ENDDO 1422 1499 ! Save bin limits (lower diameter) to be delivered to PALM if needed 1500 DO cc = 1, nbins_aerosol 1501 bin_low_limits(cc) = ( aero(cc)%vlolim / api6 )**0.33333333_wp 1502 ENDDO 1503 1423 1504 END SUBROUTINE set_sizebins 1424 1505 1425 1506 !! 1426 1507 ! Description: … … 1434 1515 !! 1435 1516 SUBROUTINE aerosol_init 1436 1437 USE arrays_3d, & 1438 ONLY: zu 1439 1440 ! USE NETCDF 1441 1442 USE netcdf_data_input_mod, & 1443 ONLY: get_attribute, get_variable, & 1444 netcdf_data_input_get_dimension_length, open_read_file 1445 1517 1518 USE netcdf_data_input_mod, & 1519 ONLY: get_attribute, get_variable, netcdf_data_input_get_dimension_length, open_read_file 1520 1446 1521 IMPLICIT NONE 1447 1448 INTEGER(iwp) :: b !< loop index: size bins 1449 INTEGER(iwp) :: c !< loop index: chemical components 1450 INTEGER(iwp) :: ee !< index: end 1451 INTEGER(iwp) :: g !< loop index: gases 1452 INTEGER(iwp) :: i !< loop index: xdirection 1453 INTEGER(iwp) :: id_faero !< NetCDF id of PIDS_SALSA 1454 INTEGER(iwp) :: id_fchem !< NetCDF id of PIDS_CHEM 1455 INTEGER(iwp) :: j !< loop index: ydirection 1456 INTEGER(iwp) :: k !< loop index: zdirection 1457 INTEGER(iwp) :: kk !< loop index: zdirection 1458 INTEGER(iwp) :: nz_file !< Number of gridpoints in file (heights) 1459 INTEGER(iwp) :: prunmode 1460 INTEGER(iwp) :: ss !< index: start 1461 LOGICAL :: netcdf_extend = .FALSE. !< Flag indicating wether netcdf 1462 !< topography input file or not 1463 REAL(wp), DIMENSION(nbins) :: core !< size of the bin mid aerosol particle, 1464 REAL(wp) :: flag !< flag to mask topography grid points 1465 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pr_gas !< gas profiles 1466 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pr_mass_fracs_a !< mass fraction 1467 !< profiles: a 1468 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pr_mass_fracs_b !< and b 1469 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pr_nsect !< sectional size 1470 !< distribution profile 1471 REAL(wp), DIMENSION(nbins) :: nsect !< size distribution (#/m3) 1472 REAL(wp), DIMENSION(0:nz+1,nbins) :: pndist !< size dist as a function 1473 !< of height (#/m3) 1474 REAL(wp), DIMENSION(0:nz+1) :: pnf2a !< number fraction: bins 2a 1475 REAL(wp), DIMENSION(0:nz+1,maxspec) :: pvf2a !< mass distributions of 1476 !< aerosol species for a 1477 REAL(wp), DIMENSION(0:nz+1,maxspec) :: pvf2b !< and bbins 1478 REAL(wp), DIMENSION(0:nz+1) :: pvfOC1a !< mass fraction between 1479 !< SO4 and OC in 1a 1480 REAL(wp), DIMENSION(:), ALLOCATABLE :: pr_z 1481 1522 1523 CHARACTER(LEN=25), DIMENSION(:), ALLOCATABLE :: cc_name !< chemical component name 1524 1525 INTEGER(iwp) :: ee !< index: end 1526 INTEGER(iwp) :: i !< loop index: xdirection 1527 INTEGER(iwp) :: ib !< loop index: size bins 1528 INTEGER(iwp) :: ic !< loop index: chemical components 1529 INTEGER(iwp) :: id_dyn !< NetCDF id of PIDS_DYNAMIC_SALSA 1530 INTEGER(iwp) :: ig !< loop index: gases 1531 INTEGER(iwp) :: j !< loop index: ydirection 1532 INTEGER(iwp) :: k !< loop index: zdirection 1533 INTEGER(iwp) :: lod_aero !< level of detail of inital aerosol concentrations 1534 INTEGER(iwp) :: pr_nbins !< Number of aerosol size bins in file 1535 INTEGER(iwp) :: pr_ncc !< Number of aerosol chemical components in file 1536 INTEGER(iwp) :: pr_nz !< Number of vertical gridpoints in file 1537 INTEGER(iwp) :: prunmode !< running mode of SALSA 1538 INTEGER(iwp) :: ss !< index: start 1539 1540 INTEGER(iwp), DIMENSION(maxspec) :: cc_input_to_model 1541 1542 LOGICAL :: netcdf_extend = .FALSE. !< Flag: netcdf file exists 1543 1544 REAL(wp) :: flag !< flag to mask topography grid points 1545 1546 REAL(wp), DIMENSION(nbins_aerosol) :: core !< size of the bin mid aerosol particle 1547 REAL(wp), DIMENSION(nbins_aerosol) :: nsect !< size distribution (#/m3) 1548 1549 REAL(wp), DIMENSION(0:nz+1) :: pnf2a !< number fraction in 2a 1550 REAL(wp), DIMENSION(0:nz+1) :: pmfoc1a !< mass fraction of OC in 1a 1551 1552 REAL(wp), DIMENSION(0:nz+1,nbins_aerosol) :: pndist !< size dist as a function of height (#/m3) 1553 REAL(wp), DIMENSION(0:nz+1,maxspec) :: pmf2a !< mass distributions in subrange 2a 1554 REAL(wp), DIMENSION(0:nz+1,maxspec) :: pmf2b !< mass distributions in subrange 2b 1555 1556 REAL(wp), DIMENSION(:), ALLOCATABLE :: pr_dmid !< vertical profile of aerosol bin diameters 1557 REAL(wp), DIMENSION(:), ALLOCATABLE :: pr_z !< z levels of profiles 1558 1559 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pr_mass_fracs_a !< mass fraction: a 1560 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pr_mass_fracs_b !< and b 1561 1562 cc_input_to_model = 0 1482 1563 prunmode = 1 1483 1564 ! 1484 1565 ! Bin mean aerosol particle volume (m3) 1485 1566 core(:) = 0.0_wp 1486 core(1:nbins ) = api6 * aero(1:nbins)%dmid ** 3.0_wp1487 ! 1488 ! Set concentrations to zero 1567 core(1:nbins_aerosol) = api6 * aero(1:nbins_aerosol)%dmid**3 1568 ! 1569 ! Set concentrations to zero 1489 1570 nsect(:) = 0.0_wp 1490 1571 pndist(:,:) = 0.0_wp 1491 pnf2a(:) = nf2a 1492 p vf2a(:,:) = 0.0_wp1493 p vf2b(:,:) = 0.0_wp1494 p vfOC1a(:) = 0.0_wp1572 pnf2a(:) = nf2a 1573 pmf2a(:,:) = 0.0_wp 1574 pmf2b(:,:) = 0.0_wp 1575 pmfoc1a(:) = 0.0_wp 1495 1576 1496 1577 IF ( isdtyp == 1 ) THEN 1497 1578 ! 1498 ! Read input profiles from PIDS_ SALSA1579 ! Read input profiles from PIDS_DYNAMIC_SALSA 1499 1580 #if defined( __netcdf ) 1500 ! 1501 ! Locationdependent size distributions and compositions. 1502 INQUIRE( FILE ='PIDS_SALSA'// TRIM( coupling_char ), EXIST=netcdf_extend )1581 ! 1582 ! Locationdependent size distributions and compositions. 1583 INQUIRE( FILE = TRIM( input_file_dynamic ) // TRIM( coupling_char ), EXIST = netcdf_extend ) 1503 1584 IF ( netcdf_extend ) THEN 1504 1585 ! 1505 ! Open file in readonly mode 1506 CALL open_read_file( 'PIDS_SALSA' // TRIM( coupling_char ), id_faero ) 1507 ! 1508 ! Input heights 1509 CALL netcdf_data_input_get_dimension_length( id_faero, nz_file, & 1510 "profile_z" ) 1511 1512 ALLOCATE( pr_z(nz_file), pr_mass_fracs_a(maxspec,nz_file), & 1513 pr_mass_fracs_b(maxspec,nz_file), pr_nsect(nbins,nz_file) ) 1514 CALL get_variable( id_faero, 'profile_z', pr_z ) 1515 ! 1516 ! Mass fracs profile: 1: H2SO4 (sulphuric acid), 2: OC (organic carbon), 1517 ! 3: BC (black carbon), 4: DU (dust), 1518 ! 5: SS (sea salt), 6: HNO3 (nitric acid), 1519 ! 7: NH3 (ammonia) 1520 CALL get_variable( id_faero, "profile_mass_fracs_a", pr_mass_fracs_a,& 1521 0, nz_file1, 0, maxspec1 ) 1522 CALL get_variable( id_faero, "profile_mass_fracs_b", pr_mass_fracs_b,& 1523 0, nz_file1, 0, maxspec1 ) 1524 CALL get_variable( id_faero, "profile_nsect", pr_nsect, 0, nz_file1,& 1525 0, nbins1 ) 1526 1527 kk = 1 1528 DO k = nzb, nz+1 1529 IF ( kk < nz_file ) THEN 1530 DO WHILE ( pr_z(kk+1) <= zu(k) ) 1531 kk = kk + 1 1532 IF ( kk == nz_file ) EXIT 1533 ENDDO 1586 ! Open file in readonly mode 1587 CALL open_read_file( input_file_dynamic // TRIM( coupling_char ), id_dyn ) 1588 ! 1589 ! Inquire dimensions: 1590 CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' ) 1591 IF ( pr_nz /= nz ) THEN 1592 WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//& 1593 'the number of numeric grid points.' 1594 CALL message( 'aerosol_init', 'PA0601', 1, 2, 0, 6, 0 ) 1595 ENDIF 1596 CALL netcdf_data_input_get_dimension_length( id_dyn, pr_ncc, 'composition_index' ) 1597 ! 1598 ! Allocate memory 1599 ALLOCATE( pr_z(1:pr_nz), pr_mass_fracs_a(nzb:nzt+1,pr_ncc), & 1600 pr_mass_fracs_b(nzb:nzt+1,pr_ncc) ) 1601 pr_mass_fracs_a = 0.0_wp 1602 pr_mass_fracs_b = 0.0_wp 1603 ! 1604 ! Read vertical levels 1605 CALL get_variable( id_dyn, 'z', pr_z ) 1606 ! 1607 ! Read name and index of chemical components 1608 CALL get_variable( id_dyn, 'composition_name', cc_name, pr_ncc ) 1609 DO ic = 1, pr_ncc 1610 SELECT CASE ( TRIM( cc_name(ic) ) ) 1611 CASE ( 'H2SO4', 'SO4', 'h2so4', 'so4' ) 1612 cc_input_to_model(1) = ic 1613 CASE ( 'OC', 'oc' ) 1614 cc_input_to_model(2) = ic 1615 CASE ( 'BC', 'bc' ) 1616 cc_input_to_model(3) = ic 1617 CASE ( 'DU', 'du' ) 1618 cc_input_to_model(4) = ic 1619 CASE ( 'SS', 'ss' ) 1620 cc_input_to_model(5) = ic 1621 CASE ( 'HNO3', 'hno3', 'NO', 'no' ) 1622 cc_input_to_model(6) = ic 1623 CASE ( 'NH3', 'nh3', 'NH', 'nh' ) 1624 cc_input_to_model(7) = ic 1625 END SELECT 1626 ENDDO 1627 1628 IF ( SUM( cc_input_to_model ) == 0 ) THEN 1629 message_string = 'None of the aerosol chemical components in ' // TRIM( & 1630 input_file_dynamic ) // ' correspond to ones applied in SALSA.' 1631 CALL message( 'salsa_mod: aerosol_init', 'PA0602', 2, 2, 0, 6, 0 ) 1632 ENDIF 1633 ! 1634 ! Vertical profiles of mass fractions of different chemical components: 1635 CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_a', pr_mass_fracs_a, & 1636 0, pr_ncc1, 0, pr_nz1 ) 1637 CALL get_variable( id_dyn, 'init_atmosphere_mass_fracs_b', pr_mass_fracs_b, & 1638 0, pr_ncc1, 0, pr_nz1 ) 1639 ! 1640 ! Match the input data with the chemical composition applied in the model 1641 DO ic = 1, maxspec 1642 ss = cc_input_to_model(ic) 1643 IF ( ss == 0 ) CYCLE 1644 pmf2a(nzb+1:nzt+1,ic) = pr_mass_fracs_a(nzb:nzt,ss) 1645 pmf2b(nzb+1:nzt+1,ic) = pr_mass_fracs_b(nzb:nzt,ss) 1646 ENDDO 1647 ! 1648 ! Aerosol concentrations: lod=1 (total PM) or lod=2 (sectional number size distribution) 1649 CALL get_attribute( id_dyn, 'lod', lod_aero, .FALSE., 'init_atmosphere_aerosol' ) 1650 IF ( lod_aero /= 2 ) THEN 1651 message_string = 'Currently only lod=2 accepted for init_atmosphere_aerosol' 1652 CALL message( 'salsa_mod: aerosol_init', 'PA0603', 2, 2, 0, 6, 0 ) 1653 ELSE 1654 ! 1655 ! Bin mean diameters in the input file 1656 CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nbins, 'Dmid') 1657 IF ( pr_nbins /= nbins_aerosol ) THEN 1658 message_string = 'Number of size bins in init_atmosphere_aerosol does not match ' & 1659 // 'with that applied in the model' 1660 CALL message( 'salsa_mod: aerosol_init', 'PA0604', 2, 2, 0, 6, 0 ) 1534 1661 ENDIF 1535 IF ( kk < nz_file ) THEN 1536 ! 1537 ! Set initial value for gas compound tracers and initial values 1538 pvf2a(k,:) = pr_mass_fracs_a(:,kk) + ( zu(k)  pr_z(kk) ) / ( & 1539 pr_z(kk+1)  pr_z(kk) ) * ( pr_mass_fracs_a(:,kk+1)& 1540  pr_mass_fracs_a(:,kk) ) 1541 pvf2b(k,:) = pr_mass_fracs_b(:,kk) + ( zu(k)  pr_z(kk) ) / ( & 1542 pr_z(kk+1)  pr_z(kk) ) * ( pr_mass_fracs_b(:,kk+1)& 1543  pr_mass_fracs_b(:,kk) ) 1544 pndist(k,:) = pr_nsect(:,kk) + ( zu(k)  pr_z(kk) ) / ( & 1545 pr_z(kk+1)  pr_z(kk) ) * ( pr_nsect(:,kk+1)  & 1546 pr_nsect(:,kk) ) 1547 ELSE 1548 pvf2a(k,:) = pr_mass_fracs_a(:,kk) 1549 pvf2b(k,:) = pr_mass_fracs_b(:,kk) 1550 pndist(k,:) = pr_nsect(:,kk) 1662 1663 ALLOCATE( pr_dmid(pr_nbins) ) 1664 pr_dmid = 0.0_wp 1665 1666 CALL get_variable( id_dyn, 'Dmid', pr_dmid ) 1667 ! 1668 ! Check whether the sectional representation conform to the one 1669 ! applied in the model 1670 IF ( ANY( ABS( ( aero(1:nbins_aerosol)%dmid  pr_dmid ) / & 1671 aero(1:nbins_aerosol)%dmid ) > 0.1_wp ) ) THEN 1672 message_string = 'Mean diameters of the aerosol size bins ' // TRIM( & 1673 input_file_dynamic ) // ' in do not conform to the sectional '// & 1674 'representation of the model.' 1675 CALL message( 'salsa_mod: aerosol_init', 'PA0605', 2, 2, 0, 6, 0 ) 1551 1676 ENDIF 1552 IF ( iso4 < 0 ) THEN 1553 pvf2a(k,1) = 0.0_wp 1554 pvf2b(k,1) = 0.0_wp 1555 ENDIF 1556 IF ( ioc < 0 ) THEN 1557 pvf2a(k,2) = 0.0_wp 1558 pvf2b(k,2) = 0.0_wp 1559 ENDIF 1560 IF ( ibc < 0 ) THEN 1561 pvf2a(k,3) = 0.0_wp 1562 pvf2b(k,3) = 0.0_wp 1563 ENDIF 1564 IF ( idu < 0 ) THEN 1565 pvf2a(k,4) = 0.0_wp 1566 pvf2b(k,4) = 0.0_wp 1567 ENDIF 1568 IF ( iss < 0 ) THEN 1569 pvf2a(k,5) = 0.0_wp 1570 pvf2b(k,5) = 0.0_wp 1571 ENDIF 1572 IF ( ino < 0 ) THEN 1573 pvf2a(k,6) = 0.0_wp 1574 pvf2b(k,6) = 0.0_wp 1575 ENDIF 1576 IF ( inh < 0 ) THEN 1577 pvf2a(k,7) = 0.0_wp 1578 pvf2b(k,7) = 0.0_wp 1579 ENDIF 1580 ! 1581 ! Then normalise the mass fraction so that SUM = 1 1582 pvf2a(k,:) = pvf2a(k,:) / SUM( pvf2a(k,:) ) 1583 IF ( SUM( pvf2b(k,:) ) > 0.0_wp ) pvf2b(k,:) = pvf2b(k,:) / & 1584 SUM( pvf2b(k,:) ) 1585 ENDDO 1586 DEALLOCATE( pr_z, pr_mass_fracs_a, pr_mass_fracs_b, pr_nsect ) 1677 ! 1678 ! Inital aerosol concentrations 1679 CALL get_variable( id_dyn, 'init_atmosphere_aerosol', pndist(nzb+1:nzt,:), & 1680 0, pr_nbins1, 0, pr_nz1 ) 1681 ENDIF 1682 ! 1683 ! Set bottom and top boundary condition (Neumann) 1684 pmf2a(nzb,:) = pmf2a(nzb+1,:) 1685 pmf2a(nzt+1,:) = pmf2a(nzt,:) 1686 pmf2b(nzb,:) = pmf2b(nzb+1,:) 1687 pmf2b(nzt+1,:) = pmf2b(nzt,:) 1688 pndist(nzb,:) = pndist(nzb+1,:) 1689 pndist(nzt+1,:) = pndist(nzt,:) 1690 1691 IF ( index_so4 < 0 ) THEN 1692 pmf2a(:,1) = 0.0_wp 1693 pmf2b(:,1) = 0.0_wp 1694 ENDIF 1695 IF ( index_oc < 0 ) THEN 1696 pmf2a(:,2) = 0.0_wp 1697 pmf2b(:,2) = 0.0_wp 1698 ENDIF 1699 IF ( index_bc < 0 ) THEN 1700 pmf2a(:,3) = 0.0_wp 1701 pmf2b(:,3) = 0.0_wp 1702 ENDIF 1703 IF ( index_du < 0 ) THEN 1704 pmf2a(:,4) = 0.0_wp 1705 pmf2b(:,4) = 0.0_wp 1706 ENDIF 1707 IF ( index_ss < 0 ) THEN 1708 pmf2a(:,5) = 0.0_wp 1709 pmf2b(:,5) = 0.0_wp 1710 ENDIF 1711 IF ( index_no < 0 ) THEN 1712 pmf2a(:,6) = 0.0_wp 1713 pmf2b(:,6) = 0.0_wp 1714 ENDIF 1715 IF ( index_nh < 0 ) THEN 1716 pmf2a(:,7) = 0.0_wp 1717 pmf2b(:,7) = 0.0_wp 1718 ENDIF 1719 1720 IF ( SUM( pmf2a ) < 0.00001_wp .AND. SUM( pmf2b ) < 0.00001_wp ) THEN 1721 message_string = 'Error in initialising mass fractions of chemical components. ' // & 1722 'Check that all chemical components are included in parameter file!' 1723 CALL message( 'salsa_mod: aerosol_init', 'PA0606', 2, 2, 0, 6, 0 ) 1724 ENDIF 1725 ! 1726 ! Then normalise the mass fraction so that SUM = 1 1727 DO k = nzb, nzt+1 1728 pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) ) 1729 IF ( SUM( pmf2b(k,:) ) > 0.0_wp ) pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) ) 1730 ENDDO 1731 1732 DEALLOCATE( pr_z, pr_mass_fracs_a, pr_mass_fracs_b ) 1733 1587 1734 ELSE 1588 message_string = 'Input file '// TRIM( 'PIDS_SALSA' ) // & 1589 TRIM( coupling_char ) // ' for SALSA missing!' 1590 CALL message( 'salsa_mod: aerosol_init', 'SA0032', 1, 2, 0, 6, 0 ) 1591 ENDIF ! netcdf_extend 1735 message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) // & 1736 ' for SALSA missing!' 1737 CALL message( 'salsa_mod: aerosol_init', 'PA0607', 1, 2, 0, 6, 0 ) 1738 1739 ENDIF ! netcdf_extend 1740 1741 #else 1742 message_string = 'isdtyp = 1 but preprocessor directive __netcdf is not used in compiling!' 1743 CALL message( 'salsa_mod: aerosol_init', 'PA0608', 1, 2, 0, 6, 0 ) 1744 1592 1745 #endif 1593 1746 1594 1747 ELSEIF ( isdtyp == 0 ) THEN 1595 1748 ! 1596 1749 ! Mass fractions for species in a and bbins 1597 IF ( iso4 > 0 ) THEN 1598 pvf2a(:,1) = mass_fracs_a(iso4) 1599 pvf2b(:,1) = mass_fracs_b(iso4) 1600 ENDIF 1601 IF ( ioc > 0 ) THEN 1602 pvf2a(:,2) = mass_fracs_a(ioc) 1603 pvf2b(:,2) = mass_fracs_b(ioc) 1604 ENDIF 1605 IF ( ibc > 0 ) THEN 1606 pvf2a(:,3) = mass_fracs_a(ibc) 1607 pvf2b(:,3) = mass_fracs_b(ibc) 1608 ENDIF 1609 IF ( idu > 0 ) THEN 1610 pvf2a(:,4) = mass_fracs_a(idu) 1611 pvf2b(:,4) = mass_fracs_b(idu) 1612 ENDIF 1613 IF ( iss > 0 ) THEN 1614 pvf2a(:,5) = mass_fracs_a(iss) 1615 pvf2b(:,5) = mass_fracs_b(iss) 1616 ENDIF 1617 IF ( ino > 0 ) THEN 1618 pvf2a(:,6) = mass_fracs_a(ino) 1619 pvf2b(:,6) = mass_fracs_b(ino) 1620 ENDIF 1621 IF ( inh > 0 ) THEN 1622 pvf2a(:,7) = mass_fracs_a(inh) 1623 pvf2b(:,7) = mass_fracs_b(inh) 1624 ENDIF 1625 DO k = nzb, nz+1 1626 pvf2a(k,:) = pvf2a(k,:) / SUM( pvf2a(k,:) ) 1627 IF ( SUM( pvf2b(k,:) ) > 0.0_wp ) pvf2b(k,:) = pvf2b(k,:) / & 1628 SUM( pvf2b(k,:) ) 1750 IF ( index_so4 > 0 ) THEN 1751 pmf2a(:,1) = mass_fracs_a(index_so4) 1752 pmf2b(:,1) = mass_fracs_b(index_so4) 1753 ENDIF 1754 IF ( index_oc > 0 ) THEN 1755 pmf2a(:,2) = mass_fracs_a(index_oc) 1756 pmf2b(:,2) = mass_fracs_b(index_oc) 1757 ENDIF 1758 IF ( index_bc > 0 ) THEN 1759 pmf2a(:,3) = mass_fracs_a(index_bc) 1760 pmf2b(:,3) = mass_fracs_b(index_bc) 1761 ENDIF 1762 IF ( index_du > 0 ) THEN 1763 pmf2a(:,4) = mass_fracs_a(index_du) 1764 pmf2b(:,4) = mass_fracs_b(index_du) 1765 ENDIF 1766 IF ( index_ss > 0 ) THEN 1767 pmf2a(:,5) = mass_fracs_a(index_ss) 1768 pmf2b(:,5) = mass_fracs_b(index_ss) 1769 ENDIF 1770 IF ( index_no > 0 ) THEN 1771 pmf2a(:,6) = mass_fracs_a(index_no) 1772 pmf2b(:,6) = mass_fracs_b(index_no) 1773 ENDIF 1774 IF ( index_nh > 0 ) THEN 1775 pmf2a(:,7) = mass_fracs_a(index_nh) 1776 pmf2b(:,7) = mass_fracs_b(index_nh) 1777 ENDIF 1778 DO k = nzb, nzt+1 1779 pmf2a(k,:) = pmf2a(k,:) / SUM( pmf2a(k,:) ) 1780 IF ( SUM( pmf2b(k,:) ) > 0.0_wp ) pmf2b(k,:) = pmf2b(k,:) / SUM( pmf2b(k,:) ) 1629 1781 ENDDO 1630 1782 1631 1783 CALL size_distribution( n_lognorm, dpg, sigmag, nsect ) 1632 1784 ! 1633 1785 ! Normalize by the given total number concentration 1634 nsect = nsect * SUM( n_lognorm ) * 1.0E+6_wp / SUM( nsect )1635 DO b = in1a, fn2b1636 pndist(:, b) = nsect(b)1786 nsect = nsect * SUM( n_lognorm ) / SUM( nsect ) 1787 DO ib = start_subrange_1a, end_subrange_2b 1788 pndist(:,ib) = nsect(ib) 1637 1789 ENDDO 1638 1790 ENDIF 1639 1791 1640 1792 IF ( igctyp == 1 ) THEN 1641 1793 ! 1642 ! Read input profiles from PIDS_CHEM 1794 ! Read input profiles from PIDS_CHEM 1643 1795 #if defined( __netcdf ) 1644 ! 1645 ! Locationdependent size distributions and compositions. 1646 INQUIRE( FILE ='PIDS_CHEM' // TRIM( coupling_char ), EXIST=netcdf_extend )1796 ! 1797 ! Locationdependent size distributions and compositions. 1798 INQUIRE( FILE = TRIM( input_file_dynamic ) // TRIM( coupling_char ), EXIST = netcdf_extend ) 1647 1799 IF ( netcdf_extend .AND. .NOT. salsa_gases_from_chem ) THEN 1648 1800 ! 1649 ! Open file in readonly mode 1650 CALL open_read_file( 'PIDS_CHEM' // TRIM( coupling_char ), id_fchem ) 1651 ! 1652 ! Input heights 1653 CALL netcdf_data_input_get_dimension_length( id_fchem, nz_file, & 1654 "profile_z" ) 1655 ALLOCATE( pr_z(nz_file), pr_gas(ngast,nz_file) ) 1656 CALL get_variable( id_fchem, 'profile_z', pr_z ) 1657 ! 1658 ! Gases: 1659 CALL get_variable( id_fchem, "profile_H2SO4", pr_gas(1,:) ) 1660 CALL get_variable( id_fchem, "profile_HNO3", pr_gas(2,:) ) 1661 CALL get_variable( id_fchem, "profile_NH3", pr_gas(3,:) ) 1662 CALL get_variable( id_fchem, "profile_OCNV", pr_gas(4,:) ) 1663 CALL get_variable( id_fchem, "profile_OCSV", pr_gas(5,:) ) 1664 1665 kk = 1 1666 DO k = nzb, nz+1 1667 IF ( kk < nz_file ) THEN 1668 DO WHILE ( pr_z(kk+1) <= zu(k) ) 1669 kk = kk + 1 1670 IF ( kk == nz_file ) EXIT 1671 ENDDO 1672 ENDIF 1673 IF ( kk < nz_file ) THEN 1674 ! 1675 ! Set initial value for gas compound tracers and initial values 1676 DO g = 1, ngast 1677 salsa_gas(g)%init(k) = pr_gas(g,kk) + ( zu(k)  pr_z(kk) ) & 1678 / ( pr_z(kk+1)  pr_z(kk) ) * & 1679 ( pr_gas(g,kk+1)  pr_gas(g,kk) ) 1680 salsa_gas(g)%conc(k,:,:) = salsa_gas(g)%init(k) 1681 ENDDO 1682 ELSE 1683 DO g = 1, ngast 1684 salsa_gas(g)%init(k) = pr_gas(g,kk) 1685 salsa_gas(g)%conc(k,:,:) = salsa_gas(g)%init(k) 1686 ENDDO 1687 ENDIF 1801 ! Open file in readonly mode 1802 CALL open_read_file( input_file_dynamic // TRIM( coupling_char ), id_dyn ) 1803 ! 1804 ! Inquire dimensions: 1805 CALL netcdf_data_input_get_dimension_length( id_dyn, pr_nz, 'z' ) 1806 IF ( pr_nz /= nz ) THEN 1807 WRITE( message_string, * ) 'Number of inifor horizontal grid points does not match '//& 1808 'the number of numeric grid points.' 1809 CALL message( 'aerosol_init', 'PA0609', 1, 2, 0, 6, 0 ) 1810 ENDIF 1811 ! 1812 ! Read vertical profiles of gases: 1813 CALL get_variable( id_dyn, 'init_atmosphere_h2so4', salsa_gas(1)%init(nzb+1:nzt) ) 1814 CALL get_variable( id_dyn, 'init_atmosphere_hno3', salsa_gas(2)%init(nzb+1:nzt) ) 1815 CALL get_variable( id_dyn, 'init_atmosphere_nh3', salsa_gas(3)%init(nzb+1:nzt) ) 1816 CALL get_variable( id_dyn, 'init_atmosphere_ocnv', salsa_gas(4)%init(nzb+1:nzt) ) 1817 CALL get_variable( id_dyn, 'init_atmosphere_ocsv', salsa_gas(5)%init(nzb+1:nzt) ) 1818 ! 1819 ! Set Neumann top and surface boundary condition for initial + initialise concentrations 1820 DO ig = 1, ngases_salsa 1821 salsa_gas(ig)%init(nzb) = salsa_gas(ig)%init(nzb+1) 1822 salsa_gas(ig)%init(nzt+1) = salsa_gas(ig)%init(nzt) 1823 DO k = nzb, nzt+1 1824 salsa_gas(ig)%conc(k,:,:) = salsa_gas(ig)%init(k) 1825 ENDDO 1688 1826 ENDDO 1689 1690 DEALLOCATE( pr_z, pr_gas ) 1827 1691 1828 ELSEIF ( .NOT. netcdf_extend .AND. .NOT. salsa_gases_from_chem ) THEN 1692 message_string = 'Input file '// TRIM( 'PIDS_CHEM' ) // & 1693 TRIM( coupling_char ) // ' for SALSA missing!' 1694 CALL message( 'salsa_mod: aerosol_init', 'SA0033', 1, 2, 0, 6, 0 ) 1695 ENDIF ! netcdf_extend 1829 message_string = 'Input file '// TRIM( input_file_dynamic ) // TRIM( coupling_char ) // & 1830 ' for SALSA missing!' 1831 CALL message( 'salsa_mod: aerosol_init', 'PA0610', 1, 2, 0, 6, 0 ) 1832 ENDIF ! netcdf_extend 1833 #else 1834 message_string = 'igctyp = 1 but preprocessor directive __netcdf is not used in compiling!' 1835 CALL message( 'salsa_mod: aerosol_init', 'PA0611', 1, 2, 0, 6, 0 ) 1836 1696 1837 #endif 1697 1838 1698 1839 ENDIF 1699 1700 IF ( ioc > 0 .AND. iso4 > 0 ) THEN 1701 ! Both are there, so use the given "massDistrA" 1702 pvfOC1a(:) = pvf2a(:,2) / ( pvf2a(:,2) + pvf2a(:,1) ) ! Normalize 1703 ELSEIF ( ioc > 0 ) THEN 1704 ! Pure organic carbon 1705 pvfOC1a(:) = 1.0_wp 1706 ELSEIF ( iso4 > 0 ) THEN 1707 ! Pure SO4 1708 pvfOC1a(:) = 0.0_wp 1840 ! 1841 ! Both SO4 and OC are included, so use the given mass fractions 1842 IF ( index_oc > 0 .AND. index_so4 > 0 ) THEN 1843 pmfoc1a(:) = pmf2a(:,2) / ( pmf2a(:,2) + pmf2a(:,1) ) ! Normalize 1844 ! 1845 ! Pure organic carbon 1846 ELSEIF ( index_oc > 0 ) THEN 1847 pmfoc1a(:) = 1.0_wp 1848 ! 1849 ! Pure SO4 1850 ELSEIF ( index_so4 > 0 ) THEN 1851 pmfoc1a(:) = 0.0_wp 1852 1709 1853 ELSE 1710 1854 message_string = 'Either OC or SO4 must be active for aerosol region 1a!' 1711 CALL message( 'salsa_mod: aerosol_init', ' SA0021', 1, 2, 0, 6, 0 )1712 ENDIF 1713 1855 CALL message( 'salsa_mod: aerosol_init', 'PA0612', 1, 2, 0, 6, 0 ) 1856 ENDIF 1857 1714 1858 ! 1715 1859 ! Initialize concentrations 1716 DO i = nxlg, nxrg 1860 DO i = nxlg, nxrg 1717 1861 DO j = nysg, nyng 1718 1862 DO k = nzb, nzt+1 1719 1863 ! 1720 ! Predetermine flag to mask topography 1864 ! Predetermine flag to mask topography 1721 1865 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1722 ! 1866 ! 1723 1867 ! a) Number concentrations 1724 ! 1725 DO b = in1a, fn1a1726 aerosol_number( b)%conc(k,j,i) = pndist(k,b) * flag1868 ! Region 1: 1869 DO ib = start_subrange_1a, end_subrange_1a 1870 aerosol_number(ib)%conc(k,j,i) = pndist(k,ib) * flag 1727 1871 IF ( prunmode == 1 ) THEN 1728 aerosol_number( b)%init = pndist(:,b)1872 aerosol_number(ib)%init = pndist(:,ib) 1729 1873 ENDIF 1730 1874 ENDDO 1731 ! 1732 ! 1875 ! 1876 ! Region 2: 1733 1877 IF ( nreg > 1 ) THEN 1734 DO b = in2a, fn2a 1735 aerosol_number(b)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * & 1736 pndist(k,b) * flag 1878 DO ib = start_subrange_2a, end_subrange_2a 1879 aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, pnf2a(k) ) * pndist(k,ib) * flag 1737 1880 IF ( prunmode == 1 ) THEN 1738 aerosol_number( b)%init = MAX( 0.0_wp, nf2a ) * pndist(:,b)1881 aerosol_number(ib)%init = MAX( 0.0_wp, nf2a ) * pndist(:,ib) 1739 1882 ENDIF 1740 1883 ENDDO 1741 1884 IF ( .NOT. no_insoluble ) THEN 1742 DO b = in2b, fn2b1743 IF ( pnf2a(k) < 1.0_wp ) THEN 1744 aerosol_number( b)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp&1745  pnf2a(k) ) * pndist(k,b) * flag1885 DO ib = start_subrange_2b, end_subrange_2b 1886 IF ( pnf2a(k) < 1.0_wp ) THEN 1887 aerosol_number(ib)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp  pnf2a(k) ) * & 1888 pndist(k,ib) * flag 1746 1889 IF ( prunmode == 1 ) THEN 1747 aerosol_number(b)%init = MAX( 0.0_wp, 1.0_wp  & 1748 nf2a ) * pndist(:,b) 1890 aerosol_number(ib)%init = MAX( 0.0_wp, 1.0_wp  nf2a ) * pndist(:,ib) 1749 1891 ENDIF 1750 1892 ENDIF … … 1755 1897 ! b) Aerosol mass concentrations 1756 1898 ! bin subrange 1: done here separately due to the SO4/OC convention 1899 ! 1757 1900 ! SO4: 1758 IF ( iso4 > 0 ) THEN 1759 ss = ( iso4  1 ) * nbins + in1a !< start 1760 ee = ( iso4  1 ) * nbins + fn1a !< end 1761 b = in1a 1762 DO c = ss, ee 1763 aerosol_mass(c)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp  & 1764 pvfOC1a(k) ) * pndist(k,b) * & 1765 core(b) * arhoh2so4 * flag 1901 IF ( index_so4 > 0 ) THEN 1902 ss = ( index_so4  1 ) * nbins_aerosol + start_subrange_1a !< start 1903 ee = ( index_so4  1 ) * nbins_aerosol + end_subrange_1a !< end 1904 ib = start_subrange_1a 1905 DO ic = ss, ee 1906 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, 1.0_wp  pmfoc1a(k) ) * pndist(k,ib)& 1907 * core(ib) * arhoh2so4 * flag 1766 1908 IF ( prunmode == 1 ) THEN 1767 aerosol_mass(c)%init = MAX( 0.0_wp, 1.0_wp  MAXVAL( & 1768 pvfOC1a ) ) * pndist(:,b) * & 1769 core(b) * arhoh2so4 1909 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, 1.0_wp  pmfoc1a(k) ) * pndist(k,ib) & 1910 * core(ib) * arhoh2so4 1770 1911 ENDIF 1771 b =b+11912 ib = ib+1 1772 1913 ENDDO 1773 1914 ENDIF 1915 ! 1774 1916 ! OC: 1775 IF ( i oc > 0 ) THEN1776 ss = ( i oc  1 ) * nbins + in1a !< start1777 ee = ( i oc  1 ) * nbins + fn1a !< end1778 b = in1a1779 DO c = ss, ee1780 aerosol_mass( c)%conc(k,j,i) = MAX( 0.0_wp, pvfOC1a(k) ) *&1781 pndist(k,b) * core(b) * arhooc * flag1917 IF ( index_oc > 0 ) THEN 1918 ss = ( index_oc  1 ) * nbins_aerosol + start_subrange_1a !< start 1919 ee = ( index_oc  1 ) * nbins_aerosol + end_subrange_1a !< end 1920 ib = start_subrange_1a 1921 DO ic = ss, ee 1922 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) * & 1923 core(ib) * arhooc * flag 1782 1924 IF ( prunmode == 1 ) THEN 1783 aerosol_mass( c)%init = MAX( 0.0_wp, MAXVAL( pvfOC1a ) )&1784 * pndist(:,b) * core(b) * arhooc1925 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmfoc1a(k) ) * pndist(k,ib) * & 1926 core(ib) * arhooc 1785 1927 ENDIF 1786 b =b+11928 ib = ib+1 1787 1929 ENDDO 1788 1930 ENDIF 1789 1790 prunmode = 3 ! Init only once1791 1792 1931 ENDDO !< k 1932 1933 prunmode = 3 ! Init only once 1934 1793 1935 ENDDO !< j 1794 1936 ENDDO !< i 1795 1937 1796 1938 ! 1797 1939 ! c) Aerosol mass concentrations 1798 1940 ! bin subrange 2: 1799 1941 IF ( nreg > 1 ) THEN 1800 1801 IF ( iso4 > 0 ) THEN 1802 CALL set_aero_mass( iso4, pvf2a(:,1), pvf2b(:,1), pnf2a, pndist, & 1803 core, arhoh2so4 ) 1804 ENDIF 1805 IF ( ioc > 0 ) THEN 1806 CALL set_aero_mass( ioc, pvf2a(:,2), pvf2b(:,2), pnf2a, pndist, core,& 1807 arhooc ) 1808 ENDIF 1809 IF ( ibc > 0 ) THEN 1810 CALL set_aero_mass( ibc, pvf2a(:,3), pvf2b(:,3), pnf2a, pndist, core,& 1811 arhobc ) 1812 ENDIF 1813 IF ( idu > 0 ) THEN 1814 CALL set_aero_mass( idu, pvf2a(:,4), pvf2b(:,4), pnf2a, pndist, core,& 1815 arhodu ) 1816 ENDIF 1817 IF ( iss > 0 ) THEN 1818 CALL set_aero_mass( iss, pvf2a(:,5), pvf2b(:,5), pnf2a, pndist, core,& 1819 arhoss ) 1820 ENDIF 1821 IF ( ino > 0 ) THEN 1822 CALL set_aero_mass( ino, pvf2a(:,6), pvf2b(:,6), pnf2a, pndist, core,& 1823 arhohno3 ) 1824 ENDIF 1825 IF ( inh > 0 ) THEN 1826 CALL set_aero_mass( inh, pvf2a(:,7), pvf2b(:,7), pnf2a, pndist, core,& 1827 arhonh3 ) 1942 1943 IF ( index_so4 > 0 ) THEN 1944 CALL set_aero_mass( index_so4, pmf2a(:,1), pmf2b(:,1), pnf2a, pndist, core, arhoh2so4 ) 1945 ENDIF 1946 IF ( index_oc > 0 ) THEN 1947 CALL set_aero_mass( index_oc, pmf2a(:,2), pmf2b(:,2), pnf2a, pndist, core, arhooc ) 1948 ENDIF 1949 IF ( index_bc > 0 ) THEN 1950 CALL set_aero_mass( index_bc, pmf2a(:,3), pmf2b(:,3), pnf2a, pndist, core, arhobc ) 1951 ENDIF 1952 IF ( index_du > 0 ) THEN 1953 CALL set_aero_mass( index_du, pmf2a(:,4), pmf2b(:,4), pnf2a, pndist, core, arhodu ) 1954 ENDIF 1955 IF ( index_ss > 0 ) THEN 1956 CALL set_aero_mass( index_ss, pmf2a(:,5), pmf2b(:,5), pnf2a, pndist, core, arhoss ) 1957 ENDIF 1958 IF ( index_no > 0 ) THEN 1959 CALL set_aero_mass( index_no, pmf2a(:,6), pmf2b(:,6), pnf2a, pndist, core, arhohno3 ) 1960 ENDIF 1961 IF ( index_nh > 0 ) THEN 1962 CALL set_aero_mass( index_nh, pmf2a(:,7), pmf2b(:,7), pnf2a, pndist, core, arhonh3 ) 1828 1963 ENDIF 1829 1964 1830 1965 ENDIF 1831 1966 1832 1967 END SUBROUTINE aerosol_init 1833 1968 1834 1969 !! 1835 1970 ! Description: … … 1839 1974 !! 1840 1975 SUBROUTINE size_distribution( in_ntot, in_dpg, in_sigma, psd_sect ) 1841 1976 1842 1977 IMPLICIT NONE 1843 1844 ! Lognormal size distribution: modes 1845 REAL(wp), DIMENSION(:), INTENT(in) :: in_dpg !< geometric mean diameter 1846 !< (micrometres) 1847 REAL(wp), DIMENSION(:), INTENT(in) :: in_ntot !< number conc. (#/cm3) 1978 1979 INTEGER(iwp) :: ib !< running index: bin 1980 INTEGER(iwp) :: iteration !< running index: iteration 1981 1982 REAL(wp) :: d1 !< particle diameter (m, dummy) 1983 REAL(wp) :: d2 !< particle diameter (m, dummy) 1984 REAL(wp) :: delta_d !< (d2d1)/10 1985 REAL(wp) :: deltadp !< bin width 1986 REAL(wp) :: dmidi !< ( d1 + d2 ) / 2 1987 1988 REAL(wp), DIMENSION(:), INTENT(in) :: in_dpg !< geometric mean diameter (m) 1989 REAL(wp), DIMENSION(:), INTENT(in) :: in_ntot !< number conc. (#/m3) 1848 1990 REAL(wp), DIMENSION(:), INTENT(in) :: in_sigma !< standard deviation 1849 REAL(wp), DIMENSION(:), INTENT(inout) :: psd_sect !< sectional size 1850 !< distribution 1851 INTEGER(iwp) :: b !< running index: bin 1852 INTEGER(iwp) :: ib !< running index: iteration 1853 REAL(wp) :: d1 !< particle diameter (m, dummy) 1854 REAL(wp) :: d2 !< particle diameter (m, dummy) 1855 REAL(wp) :: delta_d !< (d2d1)/10 1856 REAL(wp) :: deltadp !< bin width 1857 REAL(wp) :: dmidi !< ( d1 + d2 ) / 2 1858 1859 DO b = in1a, fn2b !< aerosol size bins 1860 psd_sect(b) = 0.0_wp 1991 1992 REAL(wp), DIMENSION(:), INTENT(inout) :: psd_sect !< sectional size distribution 1993 1994 DO ib = start_subrange_1a, end_subrange_2b 1995 psd_sect(ib) = 0.0_wp 1996 ! 1861 1997 ! Particle diameter at the low limit (largest in the bin) (m) 1862 d1 = ( aero(b)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp ) 1998 d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp 1999 ! 1863 2000 ! Particle diameter at the high limit (smallest in the bin) (m) 1864 d2 = ( aero(b)%vhilim / api6 ) ** ( 1.0_wp / 3.0_wp ) 2001 d2 = ( aero(ib)%vhilim / api6 )**0.33333333_wp 2002 ! 1865 2003 ! Span of particle diameter in a bin (m) 1866 delta_d = ( d2  d1 ) / 10.0_wp1867 !  Iterate:1868 DO ib = 1, 10 1869 d1 = ( aero(b)%vlolim / api6 ) ** ( 1.0_wp / 3.0_wp ) + ( ib  1) &1870 2004 delta_d = 0.1_wp * ( d2  d1 ) 2005 ! 2006 ! Iterate: 2007 DO iteration = 1, 10 2008 d1 = ( aero(ib)%vlolim / api6 )**0.33333333_wp + ( ib  1) * delta_d 1871 2009 d2 = d1 + delta_d 1872 dmidi = ( d1 + d2 ) / 2.0_wp2010 dmidi = 0.5_wp * ( d1 + d2 ) 1873 2011 deltadp = LOG10( d2 / d1 ) 1874 2012 ! 1875 2013 ! Size distribution 1876 2014 ! in_ntot = total number, total area, or total volume concentration 1877 2015 ! in_dpg = geometricmean number, area, or volume diameter 1878 2016 ! n(k) = number, area, or volume concentration in a bin 1879 ! n_lognorm and dpg converted to units of #/m3 and m 1880 psd_sect(b) = psd_sect(b) + SUM( in_ntot * 1.0E+6_wp * deltadp / & 1881 ( SQRT( 2.0_wp * pi ) * LOG10( in_sigma ) ) * & 1882 EXP( LOG10( dmidi / ( 1.0E6_wp * in_dpg ) )**2.0_wp / & 1883 ( 2.0_wp * LOG10( in_sigma ) ** 2.0_wp ) ) ) 1884 2017 psd_sect(ib) = psd_sect(ib) + SUM( in_ntot * deltadp / ( SQRT( 2.0_wp * pi ) * & 2018 LOG10( in_sigma ) ) * EXP( LOG10( dmidi / in_dpg )**2.0_wp / & 2019 ( 2.0_wp * LOG10( in_sigma ) ** 2.0_wp ) ) ) 2020 1885 2021 ENDDO 1886 2022 ENDDO 1887 2023 1888 2024 END SUBROUTINE size_distribution 1889 2025 … … 1895 2031 !> Tomi Raatikainen, FMI, 29.2.2016 1896 2032 !! 1897 SUBROUTINE set_aero_mass( ispec, p pvf2a, ppvf2b, ppnf2a, ppndist, pcore, prho )1898 2033 SUBROUTINE set_aero_mass( ispec, pmf2a, pmf2b, pnf2a, pndist, pcore, prho ) 2034 1899 2035 IMPLICIT NONE 1900 2036 2037 INTEGER(iwp) :: ee !< index: end 2038 INTEGER(iwp) :: i !< loop index 2039 INTEGER(iwp) :: ib !< loop index 2040 INTEGER(iwp) :: ic !< loop index 2041 INTEGER(iwp) :: j !< loop index 2042 INTEGER(iwp) :: k !< loop index 2043 INTEGER(iwp) :: prunmode !< 1 = initialise 2044 INTEGER(iwp) :: ss !< index: start 2045 1901 2046 INTEGER(iwp), INTENT(in) :: ispec !< Aerosol species index 1902 REAL(wp), INTENT(in) :: pcore(nbins) !< Aerosol bin mid core volume 1903 REAL(wp), INTENT(in) :: ppndist(0:nz+1,nbins) !< Aerosol size distribution 1904 REAL(wp), INTENT(in) :: ppnf2a(0:nz+1) !< Number fraction for 2a 1905 REAL(wp), INTENT(in) :: ppvf2a(0:nz+1) !< Mass distributions for a 1906 REAL(wp), INTENT(in) :: ppvf2b(0:nz+1) !< and b bins 2047 2048 REAL(wp) :: flag !< flag to mask topography grid points 2049 1907 2050 REAL(wp), INTENT(in) :: prho !< Aerosol density 1908 INTEGER(iwp) :: b !< loop index 1909 INTEGER(iwp) :: c !< loop index 1910 INTEGER(iwp) :: ee !< index: end 1911 INTEGER(iwp) :: i !< loop index 1912 INTEGER(iwp) :: j !< loop index 1913 INTEGER(iwp) :: k !< loop index 1914 INTEGER(iwp) :: prunmode !< 1 = initialise 1915 INTEGER(iwp) :: ss !< index: start 1916 REAL(wp) :: flag !< flag to mask topography grid points 1917 2051 2052 REAL(wp), DIMENSION(nbins_aerosol), INTENT(in) :: pcore !< Aerosol bin mid core volume 2053 REAL(wp), DIMENSION(0:nz+1), INTENT(in) :: pnf2a !< Number fraction for 2a 2054 REAL(wp), DIMENSION(0:nz+1), INTENT(in) :: pmf2a !< Mass distributions for a 2055 REAL(wp), DIMENSION(0:nz+1), INTENT(in) :: pmf2b !< and b bins 2056 2057 REAL(wp), DIMENSION(0:nz+1,nbins_aerosol), INTENT(in) :: pndist !< Aerosol size distribution 2058 1918 2059 prunmode = 1 1919 1920 DO i = nxlg, nxrg 2060 2061 DO i = nxlg, nxrg 1921 2062 DO j = nysg, nyng 1922 DO k = nzb, nzt+1 2063 DO k = nzb, nzt+1 1923 2064 ! 1924 2065 ! Predetermine flag to mask topography 1925 2066 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1926 ! 2067 ! 1927 2068 ! Regime 2a: 1928 ss = ( ispec  1 ) * nbins + in2a1929 ee = ( ispec  1 ) * nbins + fn2a1930 b = in2a1931 DO c = ss, ee1932 aerosol_mass( c)%conc(k,j,i) = MAX( 0.0_wp, ppvf2a(k) ) *&1933 ppnf2a(k) * ppndist(k,b) * pcore(b) * prho * flag2069 ss = ( ispec  1 ) * nbins_aerosol + start_subrange_2a 2070 ee = ( ispec  1 ) * nbins_aerosol + end_subrange_2a 2071 ib = start_subrange_2a 2072 DO ic = ss, ee 2073 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) * & 2074 pcore(ib) * prho * flag 1934 2075 IF ( prunmode == 1 ) THEN 1935 aerosol_mass(c)%init = MAX( 0.0_wp, MAXVAL( ppvf2a(:) ) ) * & 1936 MAXVAL( ppnf2a ) * pcore(b) * prho * & 1937 MAXVAL( ppndist(:,b) ) 2076 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2a(k) ) * pnf2a(k) * pndist(k,ib) * & 2077 pcore(ib) * prho 1938 2078 ENDIF 1939 b = b+12079 ib = ib + 1 1940 2080 ENDDO 2081 ! 1941 2082 ! Regime 2b: 1942 2083 IF ( .NOT. no_insoluble ) THEN 1943 ss = ( ispec  1 ) * nbins + in2b 1944 ee = ( ispec  1 ) * nbins + fn2b 1945 b = in2a 1946 DO c = ss, ee 1947 aerosol_mass(c)%conc(k,j,i) = MAX( 0.0_wp, ppvf2b(k) ) * ( & 1948 1.0_wp  ppnf2a(k) ) * ppndist(k,b) * & 1949 pcore(b) * prho * flag 2084 ss = ( ispec  1 ) * nbins_aerosol + start_subrange_2b 2085 ee = ( ispec  1 ) * nbins_aerosol + end_subrange_2b 2086 ib = start_subrange_2a 2087 DO ic = ss, ee 2088 aerosol_mass(ic)%conc(k,j,i) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp  pnf2a(k) ) *& 2089 pndist(k,ib) * pcore(ib) * prho * flag 1950 2090 IF ( prunmode == 1 ) THEN 1951 aerosol_mass(c)%init = MAX( 0.0_wp, MAXVAL( ppvf2b(:) ) )& 1952 * ( 1.0_wp  MAXVAL( ppnf2a ) ) * & 1953 MAXVAL( ppndist(:,b) ) * pcore(b) * prho 2091 aerosol_mass(ic)%init(k) = MAX( 0.0_wp, pmf2b(k) ) * ( 1.0_wp  pnf2a(k) ) * & 2092 pndist(k,ib) * pcore(ib) * prho 1954 2093 ENDIF 1955 b = b+1 1956 ENDDO 2094 ib = ib + 1 2095 ENDDO ! c 2096 1957 2097 ENDIF 1958 prunmode = 3 ! Init only once 1959 ENDDO 2098 ENDDO ! k 2099 2100 prunmode = 3 ! Init only once 2101 2102 ENDDO ! j 2103 ENDDO ! i 2104 2105 END SUBROUTINE set_aero_mass 2106 2107 !! 2108 ! Description: 2109 !  2110 !> Initialise the matching between surface types in LSM and deposition models. 2111 !> Do the matching based on Zhang et al. (2001). Atmos. Environ. 35, 549560 2112 !> (here referred as Z01). 2113 !! 2114 SUBROUTINE init_deposition 2115 2116 USE surface_mod, & 2117 ONLY: surf_lsm_h, surf_lsm_v 2118 2119 IMPLICIT NONE 2120 2121 INTEGER(iwp) :: l !< loop index for vertical surfaces 2122 2123 IF ( nldepo_surf .AND. land_surface ) THEN 2124 2125 ALLOCATE( lsm_to_depo_h%match(1:surf_lsm_h%ns) ) 2126 lsm_to_depo_h%match = 0 2127 CALL match_lsm_zhang( surf_lsm_h, lsm_to_depo_h%match ) 2128 2129 DO l = 0, 3 2130 ALLOCATE( lsm_to_depo_v(l)%match(1:surf_lsm_v(l)%ns) ) 2131 lsm_to_depo_v(l)%match = 0 2132 CALL match_lsm_zhang( surf_lsm_v(l), lsm_to_depo_v(l)%match ) 1960 2133 ENDDO 2134 ENDIF 2135 2136 IF ( nldepo_pcm ) THEN 2137 SELECT CASE ( depo_pcm_type ) 2138 CASE ( 'evergreen_needleleaf' ) 2139 depo_pcm_type_num = 1 2140 CASE ( 'evergreen_broadleaf' ) 2141 depo_pcm_type_num = 2 2142 CASE ( 'deciduous_needleleaf' ) 2143 depo_pcm_type_num = 3 2144 CASE ( 'deciduous_broadleaf' ) 2145 depo_pcm_type_num = 4 2146 CASE DEFAULT 2147 message_string = 'depo_pcm_type not set correctly.' 2148 CALL message( 'salsa_mod: init_deposition', 'PA0613', 1, 2, 0, 6, 0 ) 2149 END SELECT 2150 ENDIF 2151 2152 END SUBROUTINE init_deposition 2153 2154 !! 2155 ! Description: 2156 !  2157 !> Match the surface types in PALM and Zhang et al. 2001 deposition module 2158 !! 2159 SUBROUTINE match_lsm_zhang( surf, match_array ) 2160 2161 USE surface_mod, & 2162 ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_type 2163 2164 IMPLICIT NONE 2165 2166 INTEGER(iwp) :: m !< index for surface elements 2167 INTEGER(iwp) :: pav_type_palm !< pavement type in PALM 2168 INTEGER(iwp) :: vege_type_palm !< vegetation type in PALM 2169 INTEGER(iwp) :: water_type_palm !< water type in PALM 2170 2171 INTEGER(iwp), DIMENSION(:), INTENT(inout) :: match_array !< array matching 2172 !< the surface types 2173 TYPE(surf_type), INTENT(in) :: surf !< respective surface type 2174 2175 DO m = 1, surf%ns 2176 2177 IF ( surf%frac(ind_veg_wall,m) > 0 ) THEN 2178 vege_type_palm = surf%vegetation_type(m) 2179 SELECT CASE ( vege_type_palm ) 2180 CASE ( 0 ) 2181 message_string = 'No vegetation type defined.' 2182 CALL message( 'salsa_mod: init_depo_surfaces', 'PA0614', 1, 2, 0, 6, 0 ) 2183 CASE ( 1 ) ! bare soil 2184 match_array(m) = 6 ! grass in Z01 2185 CASE ( 2 ) ! crops, mixed farming 2186 match_array(m) = 7 ! crops, mixed farming Z01 2187 CASE ( 3 ) ! short grass 2188 match_array(m) = 6 ! grass in Z01 2189 CASE ( 4 ) ! evergreen needleleaf trees 2190 match_array(m) = 1 ! evergreen needleleaf trees in Z01 2191 CASE ( 5 ) ! deciduous needleleaf trees 2192 match_array(m) = 3 ! deciduous needleleaf trees in Z01 2193 CASE ( 6 ) ! evergreen broadleaf trees 2194 match_array(m) = 2 ! evergreen broadleaf trees in Z01 2195 CASE ( 7 ) ! deciduous broadleaf trees 2196 match_array(m) = 4 ! deciduous broadleaf trees in Z01 2197 CASE ( 8 ) ! tall grass 2198 match_array(m) = 6 ! grass in Z01 2199 CASE ( 9 ) ! desert 2200 match_array(m) = 8 ! desert in Z01 2201 CASE ( 10 ) ! tundra 2202 match_array(m) = 9 ! tundra in Z01 2203 CASE ( 11 ) ! irrigated crops 2204 match_array(m) = 7 ! crops, mixed farming Z01 2205 CASE ( 12 ) ! semidesert 2206 match_array(m) = 8 ! desert in Z01 2207 CASE ( 13 ) ! ice caps and glaciers 2208 match_array(m) = 12 ! ice cap and glacier in Z01 2209 CASE ( 14 ) ! bogs and marshes 2210 match_array(m) = 11 ! wetland with plants in Z01 2211 CASE ( 15 ) ! evergreen shrubs 2212 match_array(m) = 10 ! shrubs and interrupted woodlands in Z01 2213 CASE ( 16 ) ! deciduous shrubs 2214 match_array(m) = 10 ! shrubs and interrupted woodlands in Z01 2215 CASE ( 17 ) ! mixed forest/woodland 2216 match_array(m) = 5 ! mixed broadleaf and needleleaf trees in Z01 2217 CASE ( 18 ) ! interrupted forest 2218 match_array(m) = 10 ! shrubs and interrupted woodlands in Z01 2219 END SELECT 2220 ENDIF 2221 2222 IF ( surf%frac(ind_pav_green,m) > 0 ) THEN 2223 pav_type_palm = surf%pavement_type(m) 2224 IF ( pav_type_palm == 0 ) THEN ! error 2225 message_string = 'No pavement type defined.' 2226 CALL message( 'salsa_mod: match_lsm_zhang', 'PA0615', 1, 2, 0, 6, 0 ) 2227 ELSEIF ( pav_type_palm > 0 .AND. pav_type_palm <= 15 ) THEN 2228 match_array(m) = 15 ! urban in Z01 2229 ENDIF 2230 ENDIF 2231 2232 IF ( surf%frac(ind_wat_win,m) > 0 ) THEN 2233 water_type_palm = surf%water_type(m) 2234 IF ( water_type_palm == 0 ) THEN ! error 2235 message_string = 'No water type defined.' 2236 CALL message( 'salsa_mod: match_lsm_zhang', 'PA0616', 1, 2, 0, 6, 0 ) 2237 ELSEIF ( water_type_palm == 3 ) THEN 2238 match_array(m) = 14 ! ocean in Z01 2239 ELSEIF ( water_type_palm == 1 .OR. water_type_palm == 2 .OR. water_type_palm == 4 & 2240 .OR. water_type_palm == 5 ) THEN 2241 match_array(m) = 13 ! inland water in Z01 2242 ENDIF 2243 ENDIF 2244 1961 2245 ENDDO 1962 END SUBROUTINE set_aero_mass 2246 2247 END SUBROUTINE match_lsm_zhang 1963 2248 1964 2249 !! … … 1971 2256 IMPLICIT NONE 1972 2257 2258 INTEGER(iwp) :: ib !< 2259 INTEGER(iwp) :: ic !< 2260 INTEGER(iwp) :: icc !< 2261 INTEGER(iwp) :: ig !< 2262 1973 2263 INTEGER(iwp), INTENT(IN) :: mod_count !< 1974 INTEGER(iwp) :: b !< 1975 INTEGER(iwp) :: c !< 1976 INTEGER(iwp) :: cc !< 1977 INTEGER(iwp) :: g !< 1978 1979 IF ( simulated_time >= time_since_reference_point ) THEN 1980 1981 SELECT CASE ( mod_count ) 1982 1983 CASE ( 0 ) 1984 1985 DO b = 1, nbins 1986 aerosol_number(b)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 1987 nconc_1(:,:,:,b) 1988 aerosol_number(b)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 1989 nconc_2(:,:,:,b) 1990 DO c = 1, ncc_tot 1991 cc = ( c1 ) * nbins + b ! required due to possible Intel18 bug 1992 aerosol_mass(cc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 1993 mconc_1(:,:,:,cc) 1994 aerosol_mass(cc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 1995 mconc_2(:,:,:,cc) 2264 2265 IF ( time_since_reference_point >= skip_time_do_salsa ) THEN 2266 2267 SELECT CASE ( mod_count ) 2268 2269 CASE ( 0 ) 2270 2271 DO ib = 1, nbins_aerosol 2272 aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,ib) 2273 aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,ib) 2274 2275 DO ic = 1, ncomponents_mass 2276 icc = ( ic1 ) * nbins_aerosol + ib 2277 aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,icc) 2278 aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,icc) 2279 ENDDO 1996 2280 ENDDO 1997 ENDDO 1998 1999 IF ( .NOT. salsa_gases_from_chem ) THEN 2000 DO g = 1, ngast 2001 salsa_gas(g)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2002 gconc_1(:,:,:,g) 2003 salsa_gas(g)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2004 gconc_2(:,:,:,g) 2281 2282 IF ( .NOT. salsa_gases_from_chem ) THEN 2283 DO ig = 1, ngases_salsa 2284 salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,ig) 2285 salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,ig) 2286 ENDDO 2287 ENDIF 2288 2289 CASE ( 1 ) 2290 2291 DO ib = 1, nbins_aerosol 2292 aerosol_number(ib)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_2(:,:,:,ib) 2293 aerosol_number(ib)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => nconc_1(:,:,:,ib) 2294 DO ic = 1, ncomponents_mass 2295 icc = ( ic1 ) * nbins_aerosol + ib 2296 aerosol_mass(icc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_2(:,:,:,icc) 2297 aerosol_mass(icc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => mconc_1(:,:,:,icc) 2298 ENDDO 2005 2299 ENDDO 2006 ENDIF 2007 2008 CASE ( 1 ) 2009 2010 DO b = 1, nbins 2011 aerosol_number(b)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2012 nconc_2(:,:,:,b) 2013 aerosol_number(b)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2014 nconc_1(:,:,:,b) 2015 DO c = 1, ncc_tot 2016 cc = ( c1 ) * nbins + b ! required due to possible Intel18 bug 2017 aerosol_mass(cc)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2018 mconc_2(:,:,:,cc) 2019 aerosol_mass(cc)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2020 mconc_1(:,:,:,cc) 2021 ENDDO 2022 ENDDO 2023 2024 IF ( .NOT. salsa_gases_from_chem ) THEN 2025 DO g = 1, ngast 2026 salsa_gas(g)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2027 gconc_2(:,:,:,g) 2028 salsa_gas(g)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => & 2029 gconc_1(:,:,:,g) 2030 ENDDO 2031 ENDIF 2032 2033 END SELECT 2300 2301 IF ( .NOT. salsa_gases_from_chem ) THEN 2302 DO ig = 1, ngases_salsa 2303 salsa_gas(ig)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_2(:,:,:,ig) 2304 salsa_gas(ig)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => gconc_1(:,:,:,ig) 2305 ENDDO 2306 ENDIF 2307 2308 END SELECT 2034 2309 2035 2310 ENDIF … … 2043 2318 !> This routine reads the respective restart data. 2044 2319 !! 2045 SUBROUTINE salsa_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & 2046 nxr_on_file, nynf, nync, nyn_on_file, nysf, & 2047 nysc, nys_on_file, tmp_3d, found ) 2048 2049 2320 SUBROUTINE salsa_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, & 2321 nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found ) 2322 2050 2323 IMPLICIT NONE 2051 2052 INTEGER(iwp) :: b !<2053 INTEGER(iwp) :: c!<2054 INTEGER(iwp) :: g!<2055 INTEGER(iwp) :: k !<2324 2325 INTEGER(iwp) :: ib !< 2326 INTEGER(iwp) :: ic !< 2327 INTEGER(iwp) :: ig !< 2328 INTEGER(iwp) :: k !< 2056 2329 INTEGER(iwp) :: nxlc !< 2057 2330 INTEGER(iwp) :: nxlf !< … … 2067 2340 INTEGER(iwp) :: nys_on_file !< 2068 2341 2069 LOGICAL, INTENT(OUT) :: found 2342 LOGICAL, INTENT(OUT) :: found !< 2070 2343 2071 2344 REAL(wp), & 2072 2345 DIMENSION(nzb:nzt+1,nys_on_filenbgp:nyn_on_file+nbgp,nxl_on_filenbgp:nxr_on_file+nbgp) :: tmp_3d !< 2073 2346 2074 2347 found = .FALSE. 2075 2348 2076 2349 IF ( read_restart_data_salsa ) THEN 2077 2350 2078 2351 SELECT CASE ( restart_string(1:length) ) 2079 2352 2080 2353 CASE ( 'aerosol_number' ) 2081 DO b = 1, nbins2354 DO ib = 1, nbins_aerosol 2082 2355 IF ( k == 1 ) READ ( 13 ) tmp_3d 2083 aerosol_number( b)%conc(:,nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = &2084 tmp_3d(:,nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp)2356 aerosol_number(ib)%conc(:,nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 2357 tmp_3d(:,nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 2085 2358 found = .TRUE. 2086 2359 ENDDO 2087 2360 2088 2361 CASE ( 'aerosol_mass' ) 2089 DO c = 1, ncc_tot * nbins2362 DO ic = 1, ncomponents_mass * nbins_aerosol 2090 2363 IF ( k == 1 ) READ ( 13 ) tmp_3d 2091 aerosol_mass( c)%conc(:,nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = &2092 tmp_3d(:,nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp)2364 aerosol_mass(ic)%conc(:,nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 2365 tmp_3d(:,nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 2093 2366 found = .TRUE. 2094 2367 ENDDO 2095 2368 2096 2369 CASE ( 'salsa_gas' ) 2097 DO g = 1, ngast2370 DO ig = 1, ngases_salsa 2098 2371 IF ( k == 1 ) READ ( 13 ) tmp_3d 2099 salsa_gas( g)%conc(:,nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = &2100 tmp_3d(:,nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp)2372 salsa_gas(ig)%conc(:,nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 2373 tmp_3d(:,nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 2101 2374 found = .TRUE. 2102 2375 ENDDO 2103 2376 2104 2377 CASE DEFAULT 2105 2378 found = .FALSE. 2106 2379 2107 2380 END SELECT 2108 2381 ENDIF 2109 2382 2110 2383 END SUBROUTINE salsa_rrd_local 2111 2112 2384 2113 2385 !! … … 2122 2394 2123 2395 IMPLICIT NONE 2124 2125 INTEGER(iwp) :: b !<2126 INTEGER(iwp) :: c!<2127 INTEGER(iwp) :: g !<2128 2396 2397 INTEGER(iwp) :: ib !< 2398 INTEGER(iwp) :: ic !< 2399 INTEGER(iwp) :: ig !< 2400 2129 2401 IF ( write_binary .AND. write_binary_salsa ) THEN 2130 2402 2131 2403 CALL wrd_write_string( 'aerosol_number' ) 2132 DO b = 1, nbins2133 WRITE ( 14 ) aerosol_number( b)%conc2404 DO ib = 1, nbins_aerosol 2405 WRITE ( 14 ) aerosol_number(ib)%conc 2134 2406 ENDDO 2135 2407 2136 2408 CALL wrd_write_string( 'aerosol_mass' ) 2137 DO c = 1, nbins*ncc_tot2138 WRITE ( 14 ) aerosol_mass( c)%conc2409 DO ic = 1, nbins_aerosol * ncomponents_mass 2410 WRITE ( 14 ) aerosol_mass(ic)%conc 2139 2411 ENDDO 2140 2412 2141 2413 CALL wrd_write_string( 'salsa_gas' ) 2142 DO g = 1, ngast2143 WRITE ( 14 ) salsa_gas( g)%conc2414 DO ig = 1, ngases_salsa 2415 WRITE ( 14 ) salsa_gas(ig)%conc 2144 2416 ENDDO 2145 2417 2146 2418 ENDIF 2147 2148 END SUBROUTINE salsa_wrd_local 2149 2419 2420 END SUBROUTINE salsa_wrd_local 2150 2421 2151 2422 !! … … 2165 2436 SUBROUTINE salsa_driver( i, j, prunmode ) 2166 2437 2167 USE arrays_3d, &2168 ONLY: pt_p, q_p, rho_air_zw,u, v, w2169 2170 USE plant_canopy_model_mod, &2438 USE arrays_3d, & 2439 ONLY: pt_p, q_p, u, v, w 2440 2441 USE plant_canopy_model_mod, & 2171 2442 ONLY: lad_s 2172 2173 USE surface_mod, & 2174 ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 2175 surf_usm_v 2176 2443 2444 USE surface_mod, & 2445 ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v 2446 2177 2447 IMPLICIT NONE 2178 2179 INTEGER(iwp), INTENT(in) :: i !< loop index 2180 INTEGER(iwp), INTENT(in) :: j !< loop index 2181 INTEGER(iwp), INTENT(in) :: prunmode !< 1: Initialization call 2182 !< 2: Spinup period call 2183 !< 3: Regular runtime call 2184 ! Local variables 2185 TYPE(t_section), DIMENSION(fn2b) :: aero_old !< helper array 2186 INTEGER(iwp) :: bb !< loop index 2187 INTEGER(iwp) :: cc !< loop index 2188 INTEGER(iwp) :: endi !< end index 2189 INTEGER(iwp) :: k_wall !< vertical index of topography top 2190 INTEGER(iwp) :: k !< loop index 2191 INTEGER(iwp) :: l !< loop index 2192 INTEGER(iwp) :: nc_h2o !< index of H2O in the prtcl index table 2193 INTEGER(iwp) :: ss !< loop index 2194 INTEGER(iwp) :: str !< start index 2195 INTEGER(iwp) :: vc !< default index in prtcl 2196 REAL(wp) :: cw_old !< previous H2O mixing ratio 2197 REAL(wp) :: flag !< flag to mask topography grid points 2198 REAL(wp), DIMENSION(nzb:nzt+1) :: in_adn !< air density (kg/m3) 2199 REAL(wp), DIMENSION(nzb:nzt+1) :: in_cs !< H2O sat. vapour conc. 2200 REAL(wp), DIMENSION(nzb:nzt+1) :: in_cw !< H2O vapour concentration 2201 REAL(wp) :: in_lad !< leaf area density (m2/m3) 2202 REAL(wp), DIMENSION(nzb:nzt+1) :: in_p !< pressure (Pa) 2203 REAL(wp) :: in_rh !< relative humidity 2204 REAL(wp), DIMENSION(nzb:nzt+1) :: in_t !< temperature (K) 2205 REAL(wp), DIMENSION(nzb:nzt+1) :: in_u !< wind magnitude (m/s) 2206 REAL(wp), DIMENSION(nzb:nzt+1) :: kvis !< kinematic viscosity of air(m2/s) 2207 REAL(wp), DIMENSION(nzb:nzt+1,fn2b) :: Sc !< particle Schmidt number 2208 REAL(wp), DIMENSION(nzb:nzt+1,fn2b) :: vd !< particle fall seed (m/s, 2209 !< sedimentation velocity) 2210 REAL(wp), DIMENSION(nzb:nzt+1) :: ppm_to_nconc !< Conversion factor 2211 !< from ppm to #/m3 2212 REAL(wp) :: zgso4 !< SO4 2213 REAL(wp) :: zghno3 !< HNO3 2214 REAL(wp) :: zgnh3 !< NH3 2215 REAL(wp) :: zgocnv !< nonvolatile OC 2216 REAL(wp) :: zgocsv !< semivolatile OC 2217 2448 2449 INTEGER(iwp) :: endi !< end index 2450 INTEGER(iwp) :: ib !< loop index 2451 INTEGER(iwp) :: ic !< loop index 2452 INTEGER(iwp) :: ig !< loop index 2453 INTEGER(iwp) :: k_wall !< vertical index of topography top 2454 INTEGER(iwp) :: k !< loop index 2455 INTEGER(iwp) :: l !< loop index 2456 INTEGER(iwp) :: nc_h2o !< index of H2O in the prtcl index table 2457 INTEGER(iwp) :: ss !< loop index 2458 INTEGER(iwp) :: str !< start index 2459 INTEGER(iwp) :: vc !< default index in prtcl 2460 2461 INTEGER(iwp), INTENT(in) :: i !< loop index 2462 INTEGER(iwp), INTENT(in) :: j !< loop index 2463 INTEGER(iwp), INTENT(in) :: prunmode !< 1: Initialization, 2: Spinup, 3: Regular runtime 2464 2465 REAL(wp) :: cw_old !< previous H2O mixing ratio 2466 REAL(wp) :: flag !< flag to mask topography grid points 2467 REAL(wp) :: in_lad !< leaf area density (m2/m3) 2468 REAL(wp) :: in_rh !< relative humidity 2469 REAL(wp) :: zgso4 !< SO4 2470 REAL(wp) :: zghno3 !< HNO3 2471 REAL(wp) :: zgnh3 !< NH3 2472 REAL(wp) :: zgocnv !< nonvolatile OC 2473 REAL(wp) :: zgocsv !< semivolatile OC 2474 2475 REAL(wp), DIMENSION(nzb:nzt+1) :: in_adn !< air density (kg/m3) 2476 REAL(wp), DIMENSION(nzb:nzt+1) :: in_cs !< H2O sat. vapour conc. 2477 REAL(wp), DIMENSION(nzb:nzt+1) :: in_cw !< H2O vapour concentration 2478 REAL(wp), DIMENSION(nzb:nzt+1) :: in_p !< pressure (Pa) 2479 REAL(wp), DIMENSION(nzb:nzt+1) :: in_t !< temperature (K) 2480 REAL(wp), DIMENSION(nzb:nzt+1) :: in_u !< wind magnitude (m/s) 2481 REAL(wp), DIMENSION(nzb:nzt+1) :: kvis !< kinematic viscosity of air(m2/s) 2482 REAL(wp), DIMENSION(nzb:nzt+1) :: ppm_to_nconc !< Conversion factor from ppm to #/m3 2483 2484 REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) :: schmidt_num !< particle Schmidt number 2485 REAL(wp), DIMENSION(nzb:nzt+1,nbins_aerosol) :: vd !< particle fall seed (m/s) 2486 2487 TYPE(t_section), DIMENSION(nbins_aerosol) :: aero_old !< helper array 2488 2218 2489 aero_old(:)%numc = 0.0_wp 2219 in_adn = 0.0_wp2220 in_cs = 0.0_wp2221 in_cw = 0.0_wp2222 2490 in_lad = 0.0_wp 2223 in_rh = 0.0_wp2224 in_p = 0.0_wp2225 in_t = 0.0_wp2226 2491 in_u = 0.0_wp 2227 2492 kvis = 0.0_wp 2228 Sc= 0.0_wp2493 schmidt_num = 0.0_wp 2229 2494 vd = 0.0_wp 2230 ppm_to_nconc = 1.0_wp2231 2495 zgso4 = nclim 2232 2496 zghno3 = nclim … … 2234 2498 zgocnv = nclim 2235 2499 zgocsv = nclim 2236 2237 ! 2500 ! 2238 2501 ! Aerosol number is always set, but mass can be uninitialized 2239 DO cc = 1, nbins2240 aero( cc)%volc= 0.0_wp2241 aero_old( cc)%volc= 0.0_wp2242 ENDDO 2243 ! 2502 DO ib = 1, nbins_aerosol 2503 aero(ib)%volc(:) = 0.0_wp 2504 aero_old(ib)%volc(:) = 0.0_wp 2505 ENDDO 2506 ! 2244 2507 ! Set the salsa runtime config (How to make this more efficient?) 2245 2508 CALL set_salsa_runtime( prunmode ) 2246 ! 2509 ! 2247 2510 ! Calculate thermodynamic quantities needed in SALSA 2248 CALL salsa_thrm_ij( i, j, p_ij=in_p, temp_ij=in_t, cw_ij=in_cw, & 2249 cs_ij=in_cs, adn_ij=in_adn ) 2511 CALL salsa_thrm_ij( i, j, p_ij=in_p, temp_ij=in_t, cw_ij=in_cw, cs_ij=in_cs, adn_ij=in_adn ) 2250 2512 ! 2251 2513 ! Magnitude of wind: needed for deposition 2252 2514 IF ( lsdepo ) THEN 2253 in_u(nzb+1:nzt) = SQRT( & 2254 ( 0.5_wp * ( u(nzb+1:nzt,j,i) + u(nzb+1:nzt,j,i+1) ) )**2 + & 2255 ( 0.5_wp * ( v(nzb+1:nzt,j,i) + v(nzb+1:nzt,j+1,i) ) )**2 + & 2256 ( 0.5_wp * ( w(nzb:nzt1,j,i) + w(nzb+1:nzt,j, i) ) )**2 ) 2515 in_u(nzb+1:nzt) = SQRT( ( 0.5_wp * ( u(nzb+1:nzt,j,i) + u(nzb+1:nzt,j,i+1) ) )**2 + & 2516 ( 0.5_wp * ( v(nzb+1:nzt,j,i) + v(nzb+1:nzt,j+1,i) ) )**2 + & 2517 ( 0.5_wp * ( w(nzb:nzt1,j,i) + w(nzb+1:nzt,j, i) ) )**2 ) 2257 2518 ENDIF 2258 2519 ! 2259 2520 ! Calculate conversion factors for gas concentrations 2260 ppm_to_nconc = for_ppm_to_nconc * in_p / in_t2521 ppm_to_nconc(:) = for_ppm_to_nconc * in_p(:) / in_t(:) 2261 2522 ! 2262 2523 ! Determine topographytop index on scalar grid 2263 k_wall = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,j,i), 12 ) ), & 2264 DIM = 1 )  1 2265 2524 k_wall = k_topo_top(j,i) 2525 2266 2526 DO k = nzb+1, nzt 2267 2527 ! 2268 2528 ! Predetermine flag to mask topography 2269 2529 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2270 ! 2271 ! Do not run inside buildings 2272 IF ( flag == 0.0_wp ) CYCLE 2273 ! 2274 ! Wind velocity for dry depositon on vegetation 2275 IF ( lsdepo_vege .AND. plant_canopy ) THEN 2276 in_lad = lad_s(kk_wall,j,i) 2277 ENDIF 2530 ! 2531 ! Wind velocity for dry depositon on vegetation 2532 IF ( lsdepo_pcm .AND. plant_canopy ) THEN 2533 in_lad = lad_s( MAX( kk_wall,0 ),j,i) 2534 ENDIF 2278 2535 ! 2279 2536 ! For initialization and spinup, limit the RH with the parameter rhlim … … 2284 2541 ENDIF 2285 2542 cw_old = in_cw(k) !* in_adn(k) 2286 ! 2543 ! 2287 2544 ! Set volume concentrations: 2288 2545 ! Sulphate (SO4) or sulphuric acid H2SO4 2289 IF ( i so4 > 0 ) THEN2546 IF ( index_so4 > 0 ) THEN 2290 2547 vc = 1 2291 str = ( i so41 ) * nbins+ 1 ! start index2292 endi = i so4 * nbins! end index2293 cc = 12548 str = ( index_so41 ) * nbins_aerosol + 1 ! start index 2549 endi = index_so4 * nbins_aerosol ! end index 2550 ic = 1 2294 2551 DO ss = str, endi 2295 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2so42296 cc = cc+12552 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2so4 2553 ic = ic+1 2297 2554 ENDDO 2298 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2299 ENDIF 2300 2555 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2556 ENDIF 2557 ! 2301 2558 ! Organic carbon (OC) compounds 2302 IF ( i oc > 0 ) THEN2559 IF ( index_oc > 0 ) THEN 2303 2560 vc = 2 2304 str = ( i oc1 ) * nbins+ 12305 endi = i oc * nbins2306 cc = 12561 str = ( index_oc1 ) * nbins_aerosol + 1 2562 endi = index_oc * nbins_aerosol 2563 ic = 1 2307 2564 DO ss = str, endi 2308 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhooc2309 cc = cc+12565 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhooc 2566 ic = ic+1 2310 2567 ENDDO 2311 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2312 ENDIF 2313 2568 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2569 ENDIF 2570 ! 2314 2571 ! Black carbon (BC) 2315 IF ( i bc > 0 ) THEN2572 IF ( index_bc > 0 ) THEN 2316 2573 vc = 3 2317 str = ( i bc1 ) * nbins + 1 + fn1a2318 endi = i bc * nbins2319 cc = 1 + fn1a2574 str = ( index_bc1 ) * nbins_aerosol + 1 + end_subrange_1a 2575 endi = index_bc * nbins_aerosol 2576 ic = 1 + end_subrange_1a 2320 2577 DO ss = str, endi 2321 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhobc2322 cc = cc+12323 ENDDO 2324 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2325 ENDIF 2326 2578 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhobc 2579 ic = ic+1 2580 ENDDO 2581 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2582 ENDIF 2583 ! 2327 2584 ! Dust (DU) 2328 IF ( i du > 0 ) THEN2585 IF ( index_du > 0 ) THEN 2329 2586 vc = 4 2330 str = ( i du1 ) * nbins + 1 + fn1a2331 endi = i du * nbins2332 cc = 1 + fn1a2587 str = ( index_du1 ) * nbins_aerosol + 1 + end_subrange_1a 2588 endi = index_du * nbins_aerosol 2589 ic = 1 + end_subrange_1a 2333 2590 DO ss = str, endi 2334 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhodu2335 cc = cc+12591 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhodu 2592 ic = ic+1 2336 2593 ENDDO 2337 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2338 ENDIF 2339 2594 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2595 ENDIF 2596 ! 2340 2597 ! Sea salt (SS) 2341 IF ( i ss > 0 ) THEN2598 IF ( index_ss > 0 ) THEN 2342 2599 vc = 5 2343 str = ( i ss1 ) * nbins + 1 + fn1a2344 endi = i ss * nbins2345 cc = 1 + fn1a2600 str = ( index_ss1 ) * nbins_aerosol + 1 + end_subrange_1a 2601 endi = index_ss * nbins_aerosol 2602 ic = 1 + end_subrange_1a 2346 2603 DO ss = str, endi 2347 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoss2348 cc = cc+12604 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoss 2605 ic = ic+1 2349 2606 ENDDO 2350 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2351 ENDIF 2352 2607 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2608 ENDIF 2609 ! 2353 2610 ! Nitrate (NO(3)) or nitric acid HNO3 2354 IF ( in o > 0 ) THEN2611 IF ( index_no > 0 ) THEN 2355 2612 vc = 6 2356 str = ( in o1 ) * nbins+ 12357 endi = in o * nbins2358 cc = 12613 str = ( index_no1 ) * nbins_aerosol + 1 2614 endi = index_no * nbins_aerosol 2615 ic = 1 2359 2616 DO ss = str, endi 2360 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhohno32361 cc = cc+12617 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhohno3 2618 ic = ic+1 2362 2619 ENDDO 2363 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2364 ENDIF 2365 2620 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2621 ENDIF 2622 ! 2366 2623 ! Ammonium (NH(4+)) or ammonia NH3 2367 IF ( in h > 0 ) THEN2624 IF ( index_nh > 0 ) THEN 2368 2625 vc = 7 2369 str = ( in h1 ) * nbins+ 12370 endi = in h * nbins2371 cc = 12626 str = ( index_nh1 ) * nbins_aerosol + 1 2627 endi = index_nh * nbins_aerosol 2628 ic = 1 2372 2629 DO ss = str, endi 2373 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhonh32374 cc = cc+12630 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhonh3 2631 ic = ic+1 2375 2632 ENDDO 2376 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2377 ENDIF 2378 2633 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2634 ENDIF 2635 ! 2379 2636 ! Water (always used) 2380 2637 nc_h2o = get_index( prtcl,'H2O' ) 2381 2638 vc = 8 2382 str = ( nc_h2o1 ) * nbins + 12383 endi = nc_h2o * nbins 2384 cc = 12639 str = ( nc_h2o1 ) * nbins_aerosol + 1 2640 endi = nc_h2o * nbins_aerosol 2641 ic = 1 2385 2642 IF ( advect_particle_water ) THEN 2386 2643 DO ss = str, endi 2387 aero( cc)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2o2388 cc = cc+12644 aero(ic)%volc(vc) = aerosol_mass(ss)%conc(k,j,i) / arhoh2o 2645 ic = ic+1 2389 2646 ENDDO 2390 2647 ELSE 2391 aero(1:nbins )%volc(vc) = mclim2392 ENDIF 2393 aero_old(1:nbins )%volc(vc) = aero(1:nbins)%volc(vc)2648 aero(1:nbins_aerosol)%volc(vc) = mclim 2649 ENDIF 2650 aero_old(1:nbins_aerosol)%volc(vc) = aero(1:nbins_aerosol)%volc(vc) 2394 2651 ! 2395 2652 ! Number concentrations (numc) and particle sizes 2396 2653 ! (dwet = wet diameter, core = dry volume) 2397 DO bb = 1, nbins 2398 aero(bb)%numc = aerosol_number(bb)%conc(k,j,i) 2399 aero_old(bb)%numc = aero(bb)%numc 2400 IF ( aero(bb)%numc > nclim ) THEN 2401 aero(bb)%dwet = ( SUM( aero(bb)%volc(:) ) / aero(bb)%numc / api6 )& 2402 **( 1.0_wp / 3.0_wp ) 2403 aero(bb)%core = SUM( aero(bb)%volc(1:7) ) / aero(bb)%numc 2654 DO ib = 1, nbins_aerosol 2655 aero(ib)%numc = aerosol_number(ib)%conc(k,j,i) 2656 aero_old(ib)%numc = aero(ib)%numc 2657 IF ( aero(ib)%numc > nclim ) THEN 2658 aero(ib)%dwet = ( SUM( aero(ib)%volc(:) ) / aero(ib)%numc / api6 )**0.33333333_wp 2659 aero(ib)%core = SUM( aero(ib)%volc(1:7) ) / aero(ib)%numc 2404 2660 ELSE 2405 aero( bb)%dwet = aero(bb)%dmid2406 aero( bb)%core = api6 * ( aero(bb)%dwet ) ** 3.0_wp2661 aero(ib)%dwet = aero(ib)%dmid 2662 aero(ib)%core = api6 * ( aero(ib)%dwet )**3 2407 2663 ENDIF 2408 2664 ENDDO 2409 ! 2665 ! 2410 2666 ! On EACH call of salsa_driver, calculate the ambient sizes of 2411 2667 ! particles by equilibrating soluble fraction of particles with water … … 2417 2673 ! 2418 2674 ! Gaseous tracer concentrations in #/m3 2419 IF ( salsa_gases_from_chem ) THEN 2420 ! 2675 IF ( salsa_gases_from_chem ) THEN 2676 ! 2421 2677 ! Convert concentrations in ppm to #/m3 2422 2678 zgso4 = chem_species(gas_index_chem(1))%conc(k,j,i) * ppm_to_nconc(k) 2423 2679 zghno3 = chem_species(gas_index_chem(2))%conc(k,j,i) * ppm_to_nconc(k) 2424 2680 zgnh3 = chem_species(gas_index_chem(3))%conc(k,j,i) * ppm_to_nconc(k) 2425 zgocnv = chem_species(gas_index_chem(4))%conc(k,j,i) * ppm_to_nconc(k) 2426 zgocsv = chem_species(gas_index_chem(5))%conc(k,j,i) * ppm_to_nconc(k) 2681 zgocnv = chem_species(gas_index_chem(4))%conc(k,j,i) * ppm_to_nconc(k) 2682 zgocsv = chem_species(gas_index_chem(5))%conc(k,j,i) * ppm_to_nconc(k) 2427 2683 ELSE 2428 zgso4 = salsa_gas(1)%conc(k,j,i) 2429 zghno3 = salsa_gas(2)%conc(k,j,i) 2430 zgnh3 = salsa_gas(3)%conc(k,j,i) 2431 zgocnv = salsa_gas(4)%conc(k,j,i) 2684 zgso4 = salsa_gas(1)%conc(k,j,i) 2685 zghno3 = salsa_gas(2)%conc(k,j,i) 2686 zgnh3 = salsa_gas(3)%conc(k,j,i) 2687 zgocnv = salsa_gas(4)%conc(k,j,i) 2432 2688 zgocsv = salsa_gas(5)%conc(k,j,i) 2433 ENDIF 2434 ! 2435 ! ***************************************! 2436 ! Run SALSA ! 2437 ! ***************************************! 2438 CALL run_salsa( in_p(k), in_cw(k), in_cs(k), in_t(k), in_u(k), & 2439 in_adn(k), in_lad, zgso4, zgocnv, zgocsv, zghno3, zgnh3,& 2440 aero, prtcl, kvis(k), Sc(k,:), vd(k,:), dt_salsa ) 2441 ! ***************************************! 2689 ENDIF 2690 ! 2691 ! Calculate aerosol processes: 2692 ! ********************************************************************************************* 2693 ! 2694 ! Coagulation 2695 IF ( lscoag ) THEN 2696 CALL coagulation( aero, dt_salsa, in_t(k), in_p(k) ) 2697 ENDIF 2698 ! 2699 ! Condensation 2700 IF ( lscnd ) THEN 2701 CALL condensation( aero, zgso4, zgocnv, zgocsv, zghno3, zgnh3, in_cw(k), in_cs(k), & 2702 in_t(k), in_p(k), dt_salsa, prtcl ) 2703 ENDIF 2704 ! 2705 ! Deposition 2706 IF ( lsdepo ) THEN 2707 CALL deposition( aero, in_t(k), in_adn(k), in_u(k), in_lad, kvis(k), schmidt_num(k,:), & 2708 vd(k,:) ) 2709 ENDIF 2710 ! 2711 ! Size distribution bin update 2712 IF ( lsdistupdate ) THEN 2713 CALL distr_update( aero ) 2714 ENDIF 2715 ! ********************************************************************************************* 2716 2442 2717 IF ( lsdepo ) sedim_vd(k,j,i,:) = vd(k,:) 2443 ! 2718 ! 2444 2719 ! Calculate changes in concentrations 2445 DO bb = 1, nbins2446 aerosol_number( bb)%conc(k,j,i) = aerosol_number(bb)%conc(k,j,i)&2447 + ( aero(bb)%numc  aero_old(bb)%numc ) * flag2720 DO ib = 1, nbins_aerosol 2721 aerosol_number(ib)%conc(k,j,i) = aerosol_number(ib)%conc(k,j,i) + ( aero(ib)%numc  & 2722 aero_old(ib)%numc ) * flag 2448 2723 ENDDO 2449 2450 IF ( i so4 > 0 ) THEN2724 2725 IF ( index_so4 > 0 ) THEN 2451 2726 vc = 1 2452 str = ( i so41 ) * nbins+ 12453 endi = i so4 * nbins2454 cc = 12727 str = ( index_so41 ) * nbins_aerosol + 1 2728 endi = index_so4 * nbins_aerosol 2729 ic = 1 2455 2730 DO ss = str, endi 2456 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2457 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2458 * arhoh2so4 * flag 2459 cc = cc+1 2731 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2732 aero_old(ic)%volc(vc) ) * arhoh2so4 * flag 2733 ic = ic+1 2460 2734 ENDDO 2461 2735 ENDIF 2462 2463 IF ( i oc > 0 ) THEN2736 2737 IF ( index_oc > 0 ) THEN 2464 2738 vc = 2 2465 str = ( i oc1 ) * nbins+ 12466 endi = i oc * nbins2467 cc = 12739 str = ( index_oc1 ) * nbins_aerosol + 1 2740 endi = index_oc * nbins_aerosol 2741 ic = 1 2468 2742 DO ss = str, endi 2469 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2470 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2471 * arhooc * flag 2472 cc = cc+1 2743 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2744 aero_old(ic)%volc(vc) ) * arhooc * flag 2745 ic = ic+1 2473 2746 ENDDO 2474 2747 ENDIF 2475 2476 IF ( i bc > 0 ) THEN2748 2749 IF ( index_bc > 0 ) THEN 2477 2750 vc = 3 2478 str = ( i bc1 ) * nbins + 1 + fn1a2479 endi = i bc * nbins2480 cc = 1 + fn1a2751 str = ( index_bc1 ) * nbins_aerosol + 1 + end_subrange_1a 2752 endi = index_bc * nbins_aerosol 2753 ic = 1 + end_subrange_1a 2481 2754 DO ss = str, endi 2482 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2483 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2484 * arhobc * flag 2485 cc = cc+1 2755 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2756 aero_old(ic)%volc(vc) ) * arhobc * flag 2757 ic = ic+1 2486 2758 ENDDO 2487 2759 ENDIF 2488 2489 IF ( i du > 0 ) THEN2760 2761 IF ( index_du > 0 ) THEN 2490 2762 vc = 4 2491 str = ( i du1 ) * nbins + 1 + fn1a2492 endi = i du * nbins2493 cc = 1 + fn1a2763 str = ( index_du1 ) * nbins_aerosol + 1 + end_subrange_1a 2764 endi = index_du * nbins_aerosol 2765 ic = 1 + end_subrange_1a 2494 2766 DO ss = str, endi 2495 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2496 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2497 * arhodu * flag 2498 cc = cc+1 2767 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2768 aero_old(ic)%volc(vc) ) * arhodu * flag 2769 ic = ic+1 2499 2770 ENDDO 2500 2771 ENDIF 2501 2502 IF ( i ss > 0 ) THEN2772 2773 IF ( index_ss > 0 ) THEN 2503 2774 vc = 5 2504 str = ( i ss1 ) * nbins + 1 + fn1a2505 endi = i ss * nbins2506 cc = 1 + fn1a2775 str = ( index_ss1 ) * nbins_aerosol + 1 + end_subrange_1a 2776 endi = index_ss * nbins_aerosol 2777 ic = 1 + end_subrange_1a 2507 2778 DO ss = str, endi 2508 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2509 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2510 * arhoss * flag 2511 cc = cc+1 2779 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2780 aero_old(ic)%volc(vc) ) * arhoss * flag 2781 ic = ic+1 2512 2782 ENDDO 2513 2783 ENDIF 2514 2515 IF ( in o > 0 ) THEN2784 2785 IF ( index_no > 0 ) THEN 2516 2786 vc = 6 2517 str = ( in o1 ) * nbins+ 12518 endi = in o * nbins2519 cc = 12787 str = ( index_no1 ) * nbins_aerosol + 1 2788 endi = index_no * nbins_aerosol 2789 ic = 1 2520 2790 DO ss = str, endi 2521 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2522 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2523 * arhohno3 * flag 2524 cc = cc+1 2791 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2792 aero_old(ic)%volc(vc) ) * arhohno3 * flag 2793 ic = ic+1 2525 2794 ENDDO 2526 2795 ENDIF 2527 2528 IF ( in h > 0 ) THEN2796 2797 IF ( index_nh > 0 ) THEN 2529 2798 vc = 7 2530 str = ( in o1 ) * nbins+ 12531 endi = in o * nbins2532 cc = 12799 str = ( index_nh1 ) * nbins_aerosol + 1 2800 endi = index_nh * nbins_aerosol 2801 ic = 1 2533 2802 DO ss = str, endi 2534 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2535 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2536 * arhonh3 * flag 2537 cc = cc+1 2803 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2804 aero_old(ic)%volc(vc) ) * arhonh3 * flag 2805 ic = ic+1 2538 2806 ENDDO 2539 2807 ENDIF 2540 2808 2541 2809 IF ( advect_particle_water ) THEN 2542 2810 nc_h2o = get_index( prtcl,'H2O' ) 2543 2811 vc = 8 2544 str = ( nc_h2o1 ) * nbins + 12545 endi = nc_h2o * nbins 2546 cc = 12812 str = ( nc_h2o1 ) * nbins_aerosol + 1 2813 endi = nc_h2o * nbins_aerosol 2814 ic = 1 2547 2815 DO ss = str, endi 2548 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) & 2549 + ( aero(cc)%volc(vc)  aero_old(cc)%volc(vc) ) & 2550 * arhoh2o * flag 2816 aerosol_mass(ss)%conc(k,j,i) = aerosol_mass(ss)%conc(k,j,i) + ( aero(ic)%volc(vc)  & 2817 aero_old(ic)%volc(vc) ) * arhoh2o * flag 2551 2818 IF ( prunmode == 1 ) THEN 2552 aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), & 2553 aerosol_mass(ss)%conc(k,j,i) ) 2819 aerosol_mass(ss)%init(k) = MAX( aerosol_mass(ss)%init(k), & 2820 aerosol_mass(ss)%conc(k,j,i) ) 2821 IF ( k == nzb+1 ) THEN 2822 aerosol_mass(ss)%init(k1) = 0.0_wp 2823 ELSEIF ( k == nzt ) THEN 2824 aerosol_mass(ss)%init(k+1) = aerosol_mass(ss)%init(k) 2825 ENDIF 2554 2826 ENDIF 2555 cc = cc+12827 ic = ic+1 2556 2828 ENDDO 2557 2829 ENDIF 2558 2830 ! 2559 2831 ! Condensation of precursor gases 2560 2832 IF ( lscndgas ) THEN 2561 IF ( salsa_gases_from_chem ) THEN 2562 ! 2833 IF ( salsa_gases_from_chem ) THEN 2834 ! 2563 2835 ! SO4 (or H2SO4) 2564 chem_species( gas_index_chem(1) )%conc(k,j,i) = & 2565 chem_species( gas_index_chem(1) )%conc(k,j,i) + & 2566 ( zgso4 / ppm_to_nconc(k)  & 2567 chem_species( gas_index_chem(1) )%conc(k,j,i) ) * flag 2568 ! 2836 ig = gas_index_chem(1) 2837 chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgso4 / & 2838 ppm_to_nconc(k)  chem_species(ig)%conc(k,j,i) ) * flag 2839 ! 2569 2840 ! HNO3 2570 chem_species( gas_index_chem(2) )%conc(k,j,i) = & 2571 chem_species( gas_index_chem(2) )%conc(k,j,i) + & 2572 ( zghno3 / ppm_to_nconc(k)  & 2573 chem_species( gas_index_chem(2) )%conc(k,j,i) ) * flag 2574 ! 2841 ig = gas_index_chem(2) 2842 chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zghno3 / & 2843 ppm_to_nconc(k)  chem_species(ig)%conc(k,j,i) ) * flag 2844 ! 2575 2845 ! NH3 2576 chem_species( gas_index_chem(3) )%conc(k,j,i) = & 2577 chem_species( gas_index_chem(3) )%conc(k,j,i) + & 2578 ( zgnh3 / ppm_to_nconc(k)  & 2579 chem_species( gas_index_chem(3) )%conc(k,j,i) ) * flag 2580 ! 2846 ig = gas_index_chem(3) 2847 chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgnh3 / & 2848 ppm_to_nconc(k)  chem_species(ig)%conc(k,j,i) ) * flag 2849 ! 2581 2850 ! nonvolatile OC 2582 chem_species( gas_index_chem(4) )%conc(k,j,i) = & 2583 chem_species( gas_index_chem(4) )%conc(k,j,i) + & 2584 ( zgocnv / ppm_to_nconc(k)  & 2585 chem_species( gas_index_chem(4) )%conc(k,j,i) ) * flag 2586 ! 2851 ig = gas_index_chem(4) 2852 chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgocnv / & 2853 ppm_to_nconc(k)  chem_species(ig)%conc(k,j,i) ) * flag 2854 ! 2587 2855 ! semivolatile OC 2588 chem_species( gas_index_chem(5) )%conc(k,j,i) = & 2589 chem_species( gas_index_chem(5) )%conc(k,j,i) + & 2590 ( zgocsv / ppm_to_nconc(k)  & 2591 chem_species( gas_index_chem(5) )%conc(k,j,i) ) * flag 2592 2856 ig = gas_index_chem(5) 2857 chem_species(ig)%conc(k,j,i) = chem_species(ig)%conc(k,j,i) + ( zgocsv / & 2858 ppm_to_nconc(k)  chem_species(ig)%conc(k,j,i) ) * flag 2859 2593 2860 ELSE 2594 ! 2861 ! 2595 2862 ! SO4 (or H2SO4) 2596 salsa_gas(1)%conc(k,j,i) = salsa_gas(1)%conc(k,j,i) + ( zgso4  &2597 2598 ! 2863 salsa_gas(1)%conc(k,j,i) = salsa_gas(1)%conc(k,j,i) + ( zgso4  & 2864 salsa_gas(1)%conc(k,j,i) ) * flag 2865 ! 2599 2866 ! HNO3 2600 salsa_gas(2)%conc(k,j,i) = salsa_gas(2)%conc(k,j,i) + ( zghno3  &2601 2602 ! 2867 salsa_gas(2)%conc(k,j,i) = salsa_gas(2)%conc(k,j,i) + ( zghno3  & 2868 salsa_gas(2)%conc(k,j,i) ) * flag 2869 ! 2603 2870 ! NH3 2604 salsa_gas(3)%conc(k,j,i) = salsa_gas(3)%conc(k,j,i) + ( zgnh3  &2605 2606 ! 2871 salsa_gas(3)%conc(k,j,i) = salsa_gas(3)%conc(k,j,i) + ( zgnh3  & 2872 salsa_gas(3)%conc(k,j,i) ) * flag 2873 ! 2607 2874 ! nonvolatile OC 2608 salsa_gas(4)%conc(k,j,i) = salsa_gas(4)%conc(k,j,i) + ( zgocnv  &2609 2610 ! 2875 salsa_gas(4)%conc(k,j,i) = salsa_gas(4)%conc(k,j,i) + ( zgocnv  & 2876 salsa_gas(4)%conc(k,j,i) ) * flag 2877 ! 2611 2878 ! semivolatile OC 2612 salsa_gas(5)%conc(k,j,i) = salsa_gas(5)%conc(k,j,i) + ( zgocsv  &2613 2879 salsa_gas(5)%conc(k,j,i) = salsa_gas(5)%conc(k,j,i) + ( zgocsv  & 2880 salsa_gas(5)%conc(k,j,i) ) * flag 2614 2881 ENDIF 2615 2882 ENDIF 2616 ! 2883 ! 2617 2884 ! Tendency of water vapour mixing ratio is obtained from the 2618 2885 ! change in RH during SALSA run. This releases heat and changes pt. … … 2621 2888 ! 2622 2889 IF ( feedback_to_palm ) THEN 2623 q_p(k,j,i) = q_p(k,j,i) + 1.0_wp / ( in_cw(k) * in_adn(k) + 1.0_wp ) & 2624 ** 2.0_wp * ( in_cw(k)  cw_old ) * in_adn(k) 2625 pt_p(k,j,i) = pt_p(k,j,i) + alv / c_p * ( in_cw(k)  cw_old ) * & 2626 in_adn(k) / ( in_cw(k) / in_adn(k) + 1.0_wp ) ** 2.0_wp& 2627 * pt_p(k,j,i) / in_t(k) 2628 ENDIF 2629 2890 q_p(k,j,i) = q_p(k,j,i) + 1.0_wp / ( in_cw(k) * in_adn(k) + 1.0_wp )**2 * & 2891 ( in_cw(k)  cw_old ) * in_adn(k) * flag 2892 pt_p(k,j,i) = pt_p(k,j,i) + alv / c_p * ( in_cw(k)  cw_old ) * in_adn(k) / ( in_cw(k) / & 2893 in_adn(k) + 1.0_wp )**2 * pt_p(k,j,i) / in_t(k) * flag 2894 ENDIF 2895 2630 2896 ENDDO ! k 2631 ! 2632 ! Set surfaces and wall fluxes due to deposition 2633 IF ( lsdepo _topo.AND. prunmode == 3 ) THEN2897 ! 2898 ! Set surfaces and wall fluxes due to deposition 2899 IF ( lsdepo .AND. lsdepo_surf .AND. prunmode == 3 ) THEN 2634 2900 IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN 2635 CALL depo_ topo( i, j, surf_def_h(0), vd, Sc, kvis, in_u, rho_air_zw)2901 CALL depo_surf( i, j, surf_def_h(0), vd, schmidt_num, kvis, in_u, .TRUE. ) 2636 2902 DO l = 0, 3 2637 CALL depo_topo( i, j, surf_def_v(l), vd, Sc, kvis, in_u, & 2638 rho_air_zw**0.0_wp ) 2903 CALL depo_surf( i, j, surf_def_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l ) 2639 2904 ENDDO 2640 2905 ELSE 2641 CALL depo_ topo( i, j, surf_usm_h, vd, Sc, kvis, in_u, rho_air_zw)2906 CALL depo_surf( i, j, surf_usm_h, vd, schmidt_num, kvis, in_u, .TRUE. ) 2642 2907 DO l = 0, 3 2643 CALL depo_topo( i, j, surf_usm_v(l), vd, Sc, kvis, in_u, & 2644 rho_air_zw**0.0_wp ) 2908 CALL depo_surf( i, j, surf_usm_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l ) 2645 2909 ENDDO 2646 CALL depo_ topo( i, j, surf_lsm_h, vd, Sc, kvis, in_u, rho_air_zw)2910 CALL depo_surf( i, j, surf_lsm_h, vd, schmidt_num, kvis, in_u, .TRUE. ) 2647 2911 DO l = 0, 3 2648 CALL depo_topo( i, j, surf_lsm_v(l), vd, Sc, kvis, in_u, & 2649 rho_air_zw**0.0_wp ) 2912 CALL depo_surf( i, j, surf_lsm_v(l), vd, schmidt_num, kvis, in_u, .FALSE., l ) 2650 2913 ENDDO 2651 2914 ENDIF 2652 2915 ENDIF 2653 2916 2654 2917 END SUBROUTINE salsa_driver 2655 2918 2656 !!2657 ! Description:2658 ! 2659 !> The SALSA subroutine2660 !> Modified for the new aerosol datatype,2661 !> Juha Tonttila, FMI, 2014.2662 !> Only aerosol processes included, Mona Kurppa, UHel, 20172663 !!2664 SUBROUTINE run_salsa( ppres, pcw, pcs, ptemp, mag_u, adn, lad, pc_h2so4, &2665 pc_ocnv, pc_ocsv, pc_hno3, pc_nh3, paero, prtcl, kvis, &2666 Sc, vc, ptstep )2667 2668 IMPLICIT NONE2669 !2670 ! Input parameters and variables2671 REAL(wp), INTENT(in) :: adn !< air density (kg/m3)2672 REAL(wp), INTENT(in) :: lad !< leaf area density (m2/m3)2673 REAL(wp), INTENT(in) :: mag_u !< magnitude of wind (m/s)2674 REAL(wp), INTENT(in) :: ppres !< atmospheric pressure at each grid2675 !< point (Pa)2676 REAL(wp), INTENT(in) :: ptemp !< temperature at each grid point (K)2677 REAL(wp), INTENT(in) :: ptstep !< time step of salsa processes (s)2678 TYPE(component_index), INTENT(in) :: prtcl !< part. component index table2679 !2680 ! Input variables that are changed within:2681 REAL(wp), INTENT(inout) :: kvis !< kinematic viscosity of air (m2/s)2682 REAL(wp), INTENT(inout) :: Sc(:) !< particle Schmidt number2683 REAL(wp), INTENT(inout) :: vc(:) !< particle fall speed (m/s,2684 !< sedimentation velocity)2685 ! Gas phase concentrations at each grid point (#/m3)2686 REAL(wp), INTENT(inout) :: pc_h2so4 !< sulphuric acid2687 REAL(wp), INTENT(inout) :: pc_hno3 !< nitric acid2688 REAL(wp), INTENT(inout) :: pc_nh3 !< ammonia2689 REAL(wp), INTENT(inout) :: pc_ocnv !< nonvolatile OC2690 REAL(wp), INTENT(inout) :: pc_ocsv !< semivolatile OC2691 REAL(wp), INTENT(inout) :: pcs !< Saturation concentration of water2692 !< vapour (kg/m3)2693 REAL(wp), INTENT(inout) :: pcw !< Water vapour concentration (kg/m3)2694 TYPE(t_section), INTENT(inout) :: paero(fn2b)2695 !2696 ! Coagulation2697 IF ( lscoag ) THEN2698 CALL coagulation( paero, ptstep, ptemp, ppres )2699 ENDIF2700 !2701 ! Condensation2702 IF ( lscnd ) THEN2703 CALL condensation( paero, pc_h2so4, pc_ocnv, pc_ocsv, pc_hno3, pc_nh3, &2704 pcw, pcs, ptemp, ppres, ptstep, prtcl )2705 ENDIF2706 !2707 ! Deposition2708 IF ( lsdepo ) THEN2709 CALL deposition( paero, ptemp, adn, mag_u, lad, kvis, Sc, vc )2710 ENDIF2711 !2712 ! Size distribution bin update2713 ! Mona: why done 3 times in SALSAstandalone?2714 IF ( lsdistupdate ) THEN2715 CALL distr_update( paero )2716 ENDIF2717 2718 END SUBROUTINE run_salsa2719 2720 2919 !! 2721 2920 ! Description: … … 2727 2926 !! 2728 2927 SUBROUTINE set_salsa_runtime( prunmode ) 2729 2928 2730 2929 IMPLICIT NONE 2731 2930 2732 2931 INTEGER(iwp), INTENT(in) :: prunmode 2733 2932 2734 2933 SELECT CASE(prunmode) 2735 2934 … … 2740 2939 lscndh2oae = .FALSE. 2741 2940 lsdepo = .FALSE. 2742 lsdepo_ vege= .FALSE.2743 lsdepo_ topo= .FALSE.2941 lsdepo_pcm = .FALSE. 2942 lsdepo_surf = .FALSE. 2744 2943 lsdistupdate = .TRUE. 2944 lspartition = .FALSE. 2745 2945 2746 2946 CASE(2) !< Spinup period … … 2756 2956 lscndh2oae = nlcndh2oae 2757 2957 lsdepo = nldepo 2758 lsdepo_ vege = nldepo_vege2759 lsdepo_ topo = nldepo_topo2958 lsdepo_pcm = nldepo_pcm 2959 lsdepo_surf = nldepo_surf 2760 2960 lsdistupdate = nldistupdate 2761 2762 2961 END SELECT 2763 2962 2764 2963 2765 END SUBROUTINE set_salsa_runtime 2964 END SUBROUTINE set_salsa_runtime 2766 2965 2767 2966 !! … … 2770 2969 !> Calculates the absolute temperature (using hydrostatic pressure), saturation 2771 2970 !> vapour pressure and mixing ratio over water, relative humidity and air 2772 !> density needed in the SALSA model. 2971 !> density needed in the SALSA model. 2773 2972 !> NOTE, no saturation adjustment takes place > the resulting water vapour 2774 2973 !> mixing ratio can be supersaturated, allowing the microphysical calculations … … 2779 2978 !! 2780 2979 SUBROUTINE salsa_thrm_ij( i, j, p_ij, temp_ij, cw_ij, cs_ij, adn_ij ) 2781 2782 USE arrays_3d, &2783 ONLY: p , pt, q, zu2784 2785 USE basic_constants_and_equations_mod, &2980 2981 USE arrays_3d, & 2982 ONLY: pt, q, zu 2983 2984 USE basic_constants_and_equations_mod, & 2786 2985 ONLY: barometric_formula, exner_function, ideal_gas_law_rho, magnus 2787 2788 USE control_parameters, &2986 2987 USE control_parameters, & 2789 2988 ONLY: pt_surface, surface_pressure 2790 2989 2791 2990 IMPLICIT NONE 2792 2793 INTEGER(iwp), INTENT(in) :: i 2794 INTEGER(iwp), INTENT(in) :: j 2795 REAL(wp), DIMENSION(:), INTENT(inout) :: adn_ij 2796 REAL(wp), DIMENSION(:), INTENT(inout) :: p_ij 2797 REAL(wp), DIMENSION(:), INTENT(inout) :: temp_ij 2798 REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: cw_ij 2799 REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: cs_ij 2800 REAL(wp), DIMENSION(nzb:nzt+1) :: e_s !< saturation vapour pressure 2801 !< over water (Pa) 2802 REAL(wp) :: t_surface !< absolute surface temperature (K) 2803 ! 2804 ! Pressure p_ijk (Pa) = hydrostatic pressure + perturbation pressure (p) 2991 2992 INTEGER(iwp), INTENT(in) :: i !< 2993 INTEGER(iwp), INTENT(in) :: j !< 2994 2995 REAL(wp) :: t_surface !< absolute surface temperature (K) 2996 2997 REAL(wp), DIMENSION(nzb:nzt+1) :: e_s !< saturation vapour pressure over water (Pa) 2998 2999 REAL(wp), DIMENSION(:), INTENT(inout) :: adn_ij !< air density (kg/m3) 3000 REAL(wp), DIMENSION(:), INTENT(inout) :: p_ij !< air pressure (Pa) 3001 REAL(wp), DIMENSION(:), INTENT(inout) :: temp_ij !< air temperature (K) 3002 3003 REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: cw_ij !< water vapour concentration (kg/m3) 3004 REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: cs_ij !< saturation water vap. conc.(kg/m3) 3005 ! 3006 ! Pressure p_ijk (Pa) = hydrostatic pressure 2805 3007 t_surface = pt_surface * exner_function( surface_pressure * 100.0_wp ) 2806 p_ij(:) = barometric_formula( zu, t_surface, surface_pressure * 100.0_wp ) & 2807 + p(:,j,i) 2808 ! 3008 p_ij(:) = barometric_formula( zu, t_surface, surface_pressure * 100.0_wp ) 3009 ! 2809 3010 ! Absolute ambient temperature (K) 2810 temp_ij(:) = pt(:,j,i) * exner_function( p_ij(:) ) 3011 temp_ij(:) = pt(:,j,i) * exner_function( p_ij(:) ) 2811 3012 ! 2812 3013 ! Air density … … 2815 3016 ! Water vapour concentration r_v (kg/m3) 2816 3017 IF ( PRESENT( cw_ij ) ) THEN 2817 cw_ij(:) = ( q(:,j,i) / ( 1.0_wp  q(:,j,i) ) ) * adn_ij(:) 3018 cw_ij(:) = ( q(:,j,i) / ( 1.0_wp  q(:,j,i) ) ) * adn_ij(:) 2818 3019 ENDIF 2819 3020 ! … … 2822 3023 e_s(:) = 611.0_wp * EXP( alv_d_rv * ( 3.6609E3_wp  1.0_wp / & 2823 3024 temp_ij(:) ) )! magnus( temp_ij(:) ) 2824 cs_ij(:) = ( 0.622_wp * e_s / ( p_ij(:)  e_s(:) ) ) * adn_ij(:) 3025 cs_ij(:) = ( 0.622_wp * e_s / ( p_ij(:)  e_s(:) ) ) * adn_ij(:) 2825 3026 ENDIF 2826 3027 2827 END SUBROUTINE salsa_thrm_ij 3028 END SUBROUTINE salsa_thrm_ij 2828 3029 2829 3030 !! … … 2862 3063 !> fxm: crashes if no sulphate or sea salt 2863 3064 !> fxm: do we really need to consider Kelvin effect for subrange 2 2864 !! 3065 !! 2865 3066 SUBROUTINE equilibration( prh, ptemp, paero, init ) 2866 3067 2867 3068 IMPLICIT NONE 2868 ! 2869 ! Input variables 2870 LOGICAL, INTENT(in) :: init !< TRUE: Initialization call 2871 !< FALSE: Normal runtime: update water 2872 !< content only for 1a 2873 REAL(wp), INTENT(in) :: prh !< relative humidity [01] 2874 REAL(wp), INTENT(in) :: ptemp !< temperature (K) 2875 ! 2876 ! Output variables 2877 TYPE(t_section), INTENT(inout) :: paero(fn2b) 2878 ! 2879 ! Local 2880 INTEGER(iwp) :: b !< loop index 3069 3070 INTEGER(iwp) :: ib !< loop index 2881 3071 INTEGER(iwp) :: counti !< loop index 2882 REAL(wp) :: zaw !< water activity [01] 2883 REAL(wp) :: zbinmol(7) !< binary molality of each components (mol/kg) 2884 REAL(wp) :: zcore !< Volume of dry particle 2885 REAL(wp) :: zdold !< Old diameter 2886 REAL(wp) :: zdwet !< Wet diameter or mean droplet diameter 2887 REAL(wp) :: zke !< Kelvin term in the KÃ¶hler equation 2888 REAL(wp) :: zlwc !< liquid water content [kg/m3air] 2889 REAL(wp) :: zrh !< Relative humidity 2890 REAL(wp) :: zvpart(7) !< volume of chem. compounds in one particle 2891 3072 3073 LOGICAL, INTENT(in) :: init !< TRUE: Initialization, FALSE: Normal runtime: update water 3074 !< content only for 1a 3075 3076 REAL(wp) :: zaw !< water activity [01] 3077 REAL(wp) :: zcore !< Volume of dry particle 3078 REAL(wp) :: zdold !< Old diameter 3079 REAL(wp) :: zdwet !< Wet diameter or mean droplet diameter 3080 REAL(wp) :: zke !< Kelvin term in the KÃ¶hler equation 3081 REAL(wp) :: zlwc !< liquid water content [kg/m3air] 3082 REAL(wp) :: zrh !< Relative humidity 3083 3084 REAL(wp), DIMENSION(maxspec) :: zbinmol !< binary molality of each components (mol/kg) 3085 REAL(wp), DIMENSION(maxspec) :: zvpart !< volume of chem. compounds in one particle 3086 3087 REAL(wp), INTENT(in) :: prh !< relative humidity [01] 3088 REAL(wp), INTENT(in) :: ptemp !< temperature (K) 3089 3090 TYPE(t_section), DIMENSION(nbins_aerosol), INTENT(inout) :: paero !< aerosol properties 3091 2892 3092 zaw = 0.0_wp 2893 zbinmol = 0.0_wp2894 zcore = 0.0_wp2895 zdold = 0.0_wp2896 zdwet = 0.0_wp2897 3093 zlwc = 0.0_wp 2898 zrh = 0.0_wp 2899 2900 ! 3094 ! 2901 3095 ! Relative humidity: 2902 3096 zrh = prh 2903 3097 zrh = MAX( zrh, 0.05_wp ) 2904 zrh = MIN( zrh, 0.98_wp) 3098 zrh = MIN( zrh, 0.98_wp) 2905 3099 ! 2906 3100 ! 1) Regime 1: sulphate and partly watersoluble OC. Done for every CALL 2907 DO b = in1a, fn1a ! size bin2908 3101 DO ib = start_subrange_1a, end_subrange_1a ! size bin 3102 2909 3103 zbinmol = 0.0_wp 2910 zdold = 1.0_wp 3104 zdold = 1.0_wp 2911 3105 zke = 1.02_wp 2912 2913 IF ( paero( b)%numc > nclim ) THEN3106 3107 IF ( paero(ib)%numc > nclim ) THEN 2914 3108 ! 2915 3109 ! Volume in one particle 2916 3110 zvpart = 0.0_wp 2917 zvpart(1:2) = paero( b)%volc(1:2) / paero(b)%numc2918 zvpart(6:7) = paero( b)%volc(6:7) / paero(b)%numc2919 ! 3111 zvpart(1:2) = paero(ib)%volc(1:2) / paero(ib)%numc 3112 zvpart(6:7) = paero(ib)%volc(6:7) / paero(ib)%numc 3113 ! 2920 3114 ! Total volume and wet diameter of one dry particle 2921 3115 zcore = SUM( zvpart(1:2) ) 2922 zdwet = paero( b)%dwet2923 3116 zdwet = paero(ib)%dwet 3117 2924 3118 counti = 0 2925 DO WHILE ( ABS( zdwet / zdold  1.0_wp ) > 1.0E2_wp ) 2926 3119 DO WHILE ( ABS( zdwet / zdold  1.0_wp ) > 1.0E2_wp ) 3120 2927 3121 zdold = MAX( zdwet, 1.0E20_wp ) 2928 3122 zaw = MAX( 1.0E3_wp, zrh / zke ) ! To avoid underflow 2929 ! 3123 ! 2930 3124 ! Binary molalities (mol/kg): 2931 3125 ! Sulphate 2932 zbinmol(1) = 1.1065495E+2_wp  3.6759197E+2_wp * zaw & 2933 + 5.0462934E+2_wp * zaw**2.0_wp & 2934  3.1543839E+2_wp * zaw**3.0_wp & 2935 + 6.770824E+1_wp * zaw**4.0_wp 2936 ! Organic carbon 2937 zbinmol(2) = 1.0_wp / ( zaw * amh2o )  1.0_wp / amh2o 2938 ! Nitric acid 2939 zbinmol(6) = 2.306844303E+1_wp  3.563608869E+1_wp * zaw & 2940  6.210577919E+1_wp * zaw**2.0_wp & 2941 + 5.510176187E+2_wp * zaw**3.0_wp & 2942  1.460055286E+3_wp * zaw**4.0_wp & 2943 + 1.894467542E+3_wp * zaw**5.0_wp & 2944  1.220611402E+3_wp * zaw**6.0_wp & 2945 + 3.098597737E+2_wp * zaw**7.0_wp 2946 ! 2947 ! Calculate the liquid water content (kg/m3air) using ZSR (see e.g. 2948 ! Eq. 10.98 in Seinfeld and Pandis (2006)) 2949 zlwc = ( paero(b)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / & 2950 zbinmol(1) + epsoc * paero(b)%volc(2) * ( arhooc / amoc ) & 2951 / zbinmol(2) + ( paero(b)%volc(6) * ( arhohno3/amhno3 ) ) & 2952 / zbinmol(6) 2953 ! 2954 ! Particle wet diameter (m) 2955 zdwet = ( zlwc / paero(b)%numc / arhoh2o / api6 + & 2956 ( SUM( zvpart(6:7) ) / api6 ) + & 2957 zcore / api6 )**( 1.0_wp / 3.0_wp ) 2958 ! 3126 zbinmol(1) = 1.1065495E+2_wp  3.6759197E+2_wp * zaw + 5.0462934E+2_wp * zaw**2  & 3127 3.1543839E+2_wp * zaw**3 + 6.770824E+1_wp * zaw**4 3128 ! Organic carbon 3129 zbinmol(2) = 1.0_wp / ( zaw * amh2o )  1.0_wp / amh2o 3130 ! Nitric acid 3131 zbinmol(6) = 2.306844303E+1_wp  3.563608869E+1_wp * zaw  6.210577919E+1_wp * zaw**2 & 3132 + 5.510176187E+2_wp * zaw**3  1.460055286E+3_wp * zaw**4 & 3133 + 1.894467542E+3_wp * zaw**5  1.220611402E+3_wp * zaw**6 & 3134 + 3.098597737E+2_wp * zaw**7 3135 ! 3136 ! Calculate the liquid water content (kg/m3air) using ZSR (see e.g. Eq. 10.98 in 3137 ! Seinfeld and Pandis (2006)) 3138 zlwc = ( paero(ib)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / zbinmol(1) + & 3139 epsoc * paero(ib)%volc(2) * ( arhooc / amoc ) / zbinmol(2) + & 3140 ( paero(ib)%volc(6) * ( arhohno3/amhno3 ) ) / zbinmol(6) 3141 ! 3142 ! Particle wet diameter (m) 3143 zdwet = ( zlwc / paero(ib)%numc / arhoh2o / api6 + ( SUM( zvpart(6:7) ) / api6 ) + & 3144 zcore / api6 )**0.33333333_wp 3145 ! 2959 3146 ! Kelvin effect (Eq. 10.85 in in Seinfeld and Pandis (2006)). Avoid 2960 3147 ! overflow. 2961 zke = EXP( MIN( 50.0_wp, & 2962 4.0_wp * surfw0 * amvh2so4 / ( abo * ptemp * zdwet ) ) ) 2963 3148 zke = EXP( MIN( 50.0_wp, 4.0_wp * surfw0 * amvh2so4 / ( abo * ptemp * zdwet ) ) ) 3149 2964 3150 counti = counti + 1 2965 3151 IF ( counti > 1000 ) THEN 2966 3152 message_string = 'Subrange 1: no convergence!' 2967 CALL message( 'salsa_mod: equilibration', 'SA0042', & 2968 1, 2, 0, 6, 0 ) 3153 CALL message( 'salsa_mod: equilibration', 'PA0617', 1, 2, 0, 6, 0 ) 2969 3154 ENDIF 2970 3155 ENDDO 2971 ! 3156 ! 2972 3157 ! Instead of lwc, use the volume concentration of water from now on 2973 3158 ! (easy to convert...) 2974 paero( b)%volc(8) = zlwc / arhoh2o2975 ! 3159 paero(ib)%volc(8) = zlwc / arhoh2o 3160 ! 2976 3161 ! If this is initialization, update the core and wet diameter 2977 3162 IF ( init ) THEN 2978 paero( b)%dwet = zdwet2979 paero( b)%core = zcore3163 paero(ib)%dwet = zdwet 3164 paero(ib)%core = zcore 2980 3165 ENDIF 2981 3166 2982 3167 ELSE 2983 3168 ! If initialization 2984 ! 1.2) empty bins given bin average values 3169 ! 1.2) empty bins given bin average values 2985 3170 IF ( init ) THEN 2986 paero( b)%dwet = paero(b)%dmid2987 paero( b)%core = api6 * paero(b)%dmid ** 3.0_wp3171 paero(ib)%dwet = paero(ib)%dmid 3172 paero(ib)%core = api6 * paero(ib)%dmid**3 2988 3173 ENDIF 2989 2990 ENDIF 2991 2992 ENDDO !<b3174 3175 ENDIF 3176 3177 ENDDO ! ib 2993 3178 ! 2994 3179 ! 2) Regime 2a: sulphate, OC, BC and sea salt … … 2996 3181 ! are computed via condensation 2997 3182 IF ( init ) THEN 2998 DO b = in2a, fn2b2999 3183 DO ib = start_subrange_2a, end_subrange_2b 3184 ! 3000 3185 ! Initialize 3001 3186 zke = 1.02_wp 3002 3187 zbinmol = 0.0_wp 3003 3188 zdold = 1.0_wp 3004 ! 3189 ! 3005 3190 ! 1) Particle properties calculated for nonempty bins 3006 IF ( paero( b)%numc > nclim ) THEN3007 ! 3191 IF ( paero(ib)%numc > nclim ) THEN 3192 ! 3008 3193 ! Volume in one particle [fxm] 3009 3194 zvpart = 0.0_wp 3010 zvpart(1:7) = paero( b)%volc(1:7) / paero(b)%numc3011 ! 3012 ! Total volume and wet diameter of one dry particle [fxm] 3195 zvpart(1:7) = paero(ib)%volc(1:7) / paero(ib)%numc 3196 ! 3197 ! Total volume and wet diameter of one dry particle [fxm] 3013 3198 zcore = SUM( zvpart(1:5) ) 3014 zdwet = paero( b)%dwet3199 zdwet = paero(ib)%dwet 3015 3200 3016 3201 counti = 0 3017 3202 DO WHILE ( ABS( zdwet / zdold  1.0_wp ) > 1.0E12_wp ) 3018 3203 3019 3204 zdold = MAX( zdwet, 1.0E20_wp ) 3020 3205 zaw = zrh / zke 3021 ! 3206 ! 3022 3207 ! Binary molalities (mol/kg): 3023 3208 ! Sulphate 3024 zbinmol(1) = 1.1065495E+2_wp  3.6759197E+2_wp * zaw & 3025 + 5.0462934E+2_wp * zaw**2  3.1543839E+2_wp * zaw**3 & 3026 + 6.770824E+1_wp * zaw**4 3027 ! Organic carbon 3028 zbinmol(2) = 1.0_wp / ( zaw * amh2o )  1.0_wp / amh2o 3209 zbinmol(1) = 1.1065495E+2_wp  3.6759197E+2_wp * zaw + 5.0462934E+2_wp * zaw**2  & 3210 3.1543839E+2_wp * zaw**3 + 6.770824E+1_wp * zaw**4 3211 ! Organic carbon 3212 zbinmol(2) = 1.0_wp / ( zaw * amh2o )  1.0_wp / amh2o 3029 3213 ! Nitric acid 3030 zbinmol(6) = 2.306844303E+1_wp  3.563608869E+1_wp * zaw & 3031  6.210577919E+1_wp * zaw**2 + 5.510176187E+2_wp * zaw**3 & 3032  1.460055286E+3_wp * zaw**4 + 1.894467542E+3_wp * zaw**5 & 3033  1.220611402E+3_wp * zaw**6 + 3.098597737E+2_wp * zaw**7 3034 ! Sea salt (natrium chloride) 3035 zbinmol(5) = 5.875248E+1_wp  1.8781997E+2_wp * zaw & 3036 + 2.7211377E+2_wp * zaw**2  1.8458287E+2_wp * zaw**3 & 3037 + 4.153689E+1_wp * zaw**4 3038 ! 3214 zbinmol(6) = 2.306844303E+1_wp  3.563608869E+1_wp * zaw  & 3215 6.210577919E+1_wp * zaw**2 + 5.510176187E+2_wp * zaw**3  & 3216 1.460055286E+3_wp * zaw**4 + 1.894467542E+3_wp * zaw**5  & 3217 1.220611402E+3_wp * zaw**6 + 3.098597737E+2_wp * zaw**7 3218 ! Sea salt (natrium chloride) 3219 zbinmol(5) = 5.875248E+1_wp  1.8781997E+2_wp * zaw + 2.7211377E+2_wp * zaw**2  & 3220 1.8458287E+2_wp * zaw**3 + 4.153689E+1_wp * zaw**4 3221 ! 3039 3222 ! Calculate the liquid water content (kg/m3air) 3040 zlwc = ( paero(b)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / & 3041 zbinmol(1) + epsoc * ( paero(b)%volc(2) * ( arhooc / & 3042 amoc ) ) / zbinmol(2) + ( paero(b)%volc(6) * ( arhohno3 & 3043 / amhno3 ) ) / zbinmol(6) + ( paero(b)%volc(5) * & 3044 ( arhoss / amss ) ) / zbinmol(5) 3045 3046 ! Particle wet radius (m) 3047 zdwet = ( zlwc / paero(b)%numc / arhoh2o / api6 + & 3048 ( SUM( zvpart(6:7) ) / api6 ) + & 3049 zcore / api6 ) ** ( 1.0_wp / 3.0_wp ) 3050 ! 3223 zlwc = ( paero(ib)%volc(1) * ( arhoh2so4 / amh2so4 ) ) / zbinmol(1) + & 3224 epsoc * ( paero(ib)%volc(2) * ( arhooc / amoc ) ) / zbinmol(2) + & 3225 ( paero(ib)%volc(6) * ( arhohno3 / amhno3 ) ) / zbinmol(6)