source: palm/trunk/SOURCE/chemistry_model_mod.f90 @ 4852

Last change on this file since 4852 was 4852, checked in by raasch, 3 years ago

re-numbering if message IDs

  • Property svn:keywords set to Id
File size: 269.6 KB
Line 
1!> @file chemistry_model_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 2017-2021 Leibniz Universitaet Hannover
17! Copyright 2017-2021 Karlsruhe Institute of Technology
18! Copyright 2017-2021 Freie Universitaet Berlin
19!--------------------------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: chemistry_model_mod.f90 4852 2021-01-25 07:18:30Z raasch $
28! re-numbering of message IDs
29!
30! 4847 2021-01-19 13:44:45Z raasch
31! error message IDs changed from CM to the default PA
32!
33! 4843 2021-01-15 15:22:11Z raasch
34! local namelist parameter added to switch off the module although the respective module namelist
35! appears in the namelist file
36!
37! 4833 2021-01-07 08:57:54Z raasch
38! openmp bugfix for time measurements of non advective processes
39!
40! 4828 2021-01-05 11:21:41Z Giersch
41! reading of namelist file and actions in case of namelist errors revised so that statement labels
42! and goto statements are not required any more
43!
44! 4768 2020-11-02 19:11:23Z suehring
45! Enable 3D data output also with 64-bit precision
46!
47! 4731 2020-10-07 13:25:11Z schwenkel
48! Move exchange_horiz from time_integration to modules
49!
50! 4671 2020-09-09 20:27:58Z pavelkrc
51! Implementation of downward facing USM and LSM surfaces
52!
53! 4637 2020-08-07 07:49:33Z suehring
54! Avoid usage of omp_lib, instead declare omp_get_thread_num explicitly
55!
56! 4636 2020-08-06 14:10:12Z resler
57! Fix bugs in OpenMP directives
58!
59! 4636 2020-08-06 14:10:12Z suehring
60! - Bugfix in variable name separation in profile-output initialization
61! - Bugfix in counting the number of chemistry profiles
62!
63! 4581 2020-06-29 08:49:58Z suehring
64! Enable output of resolved-scale vertical fluxes of chemical species.
65!
66! 4577 2020-06-25 09:53:58Z raasch
67! further re-formatting to follow the PALM coding standard
68!
69! 4559 2020-06-11 08:51:48Z raasch
70! file re-formatted to follow the PALM coding standard
71!
72! 4550 2020-05-29 15:22:13Z raasch
73! bugfix for reading local restart data
74!
75! 4544 2020-05-21 14:43:05Z raasch
76! conc_av changed from pointer to allocatable array, array spec_conc_av removed
77!
78! 4542 2020-05-19 15:45:12Z raasch
79! redundant if statement removed
80!
81! 4535 2020-05-15 12:07:23Z raasch
82! bugfix for restart data format query
83!
84! 4517 2020-05-03 14:29:30Z raasch
85! added restart with MPI-IO
86!
87! 4511 2020-04-30 12:20:40Z raasch
88! decycling replaced by explicit setting of lateral boundary conditions
89!
90! 4487 2020-04-03 09:38:20Z raasch
91! bugfix for subroutine calls that contain the decycle_chem switches as arguments
92!
93! 4481 2020-03-31 18:55:54Z maronga
94! use statement for exchange horiz added,
95! bugfix for call of exchange horiz 2d
96!
97! 4442 2020-03-04 19:21:13Z suehring
98! Change order of dimension in surface array %frac to allow for better
99! vectorization.
100!
101! 4441 2020-03-04 19:20:35Z suehring
102! in subroutine chem_init (ECC)
103! - allows different init paths emission data for legacy mode emission and on-demand mode in
104!   subroutine chem_init_internal (ECC)
105! - reads netcdf file only when legacy mode is activated (i.e., emiss_read_legacy_mode = .TRUE.)
106!   otherwise file is read once at the beginning to obtain header information, and emission data
107!   are extracted on an on-demand basis
108!
109! 4372 2020-01-14 10:20:35Z banzhafs
110! chem_parin : added handler for new namelist item emiss_legacy_read_mode (ECC) added messages
111! CM0465 - legacy read mode selection
112! CM0466 - legacy read mode force selection
113! CM0467 - new read mode selection
114!
115! 4370 2020-01-10 14:00:44Z raasch
116! vector directives added to force vectorization on Intel19 compiler
117!
118! 4346 2019-12-18 11:55:56Z motisi
119! Introduction of wall_flags_total_0, which currently sets bits based on static topography
120! information used in wall_flags_static_0
121!
122! 4329 2019-12-10 15:46:36Z motisi
123! Renamed wall_flags_0 to wall_flags_static_0
124!
125! 4306 2019-11-25 12:04:48Z banzhafs
126! Corretion for r4304 commit
127!
128! 4304 2019-11-25 10:43:03Z banzhafs
129! Precision clean-up in drydepo_aero_zhang_vd subroutine
130!
131! 4292 2019-11-11 13:04:50Z banzhafs
132! Bugfix for r4290
133!
134! 4290 2019-11-11 12:06:14Z banzhafs
135! Bugfix in sedimentation resistance calculation in drydepo_aero_zhang_vd subroutine
136!
137! 4273 2019-10-24 13:40:54Z monakurppa
138! Add logical switches nesting_chem and nesting_offline_chem (both .TRUE. by default)
139!
140! 4272 2019-10-23 15:18:57Z schwenkel
141! Further modularization of boundary conditions: moved boundary conditions to respective modules
142!
143! 4268 2019-10-17 11:29:38Z schwenkel
144! Moving module specific boundary conditions from time_integration to module
145!
146! 4230 2019-09-11 13:58:14Z suehring
147! Bugfix, initialize mean profiles also in restart runs. Also initialize array used for Runge-Kutta
148! tendecies in restart runs.
149!
150! 4227 2019-09-10 18:04:34Z gronemeier
151! implement new palm_date_time_mod
152!
153! 4182 2019-08-22 15:20:23Z scharf
154! Corrected "Former revisions" section
155!
156! 4167 2019-08-16 11:01:48Z suehring
157! Changed behaviour of masked output over surface to follow terrain and ignore buildings
158! (J.Resler, T.Gronemeier)
159!
160! 4166 2019-08-16 07:54:21Z resler
161! Bugfix in decycling
162!
163! 4115 2019-07-24 12:50:49Z suehring
164! Fix faulty IF statement in decycling initialization
165!
166! 4110 2019-07-22 17:05:21Z suehring
167! - Decycling boundary conditions are only set at the ghost points not on the prognostic grid points
168! - Allocation and initialization of special advection flags cs_advc_flags_s used for chemistry.
169!   These are exclusively used for chemical species in order to distinguish from the usually-used
170!   flags which might be different when decycling is applied in combination with cyclic boundary
171!   conditions.
172!   Moreover, cs_advc_flags_s considers extended zones around buildings where first-order upwind
173!   scheme is applied for the horizontal advection terms, in order to overcome high concentration
174!   peaks due to stationary numerical oscillations caused by horizontal advection discretization.
175!
176! 4109 2019-07-22 17:00:34Z suehring
177! Slightly revise setting of boundary conditions at horizontal walls, use data-structure offset
178! index instead of pre-calculate it for each facing
179!
180! 4080 2019-07-09 18:17:37Z suehring
181! Restore accidantly removed limitation to positive values
182!
183! 4079 2019-07-09 18:04:41Z suehring
184! Application of monotonic flux limiter for the vertical scalar advection up to the topography top
185! (only for the cache-optimized version at the moment).
186!
187! 4069 2019-07-01 14:05:51Z Giersch
188! Masked output running index mid has been introduced as a local variable to avoid runtime error
189! (Loop variable has been modified) in time_integration
190!
191! 4029 2019-06-14 14:04:35Z raasch
192! nest_chemistry option removed
193!
194! 4004 2019-05-24 11:32:38Z suehring
195! in subroutine chem_parin check emiss_lod / mod_emis only when emissions_anthropogenic is activated
196! in namelist (E.C. Chan)
197!
198! 3968 2019-05-13 11:04:01Z suehring
199! - added "emiss_lod" which serves the same function as "mode_emis" both will be synchronized with
200!   emiss_lod having pirority over mode_emis (see informational messages)
201! - modified existing error message and introduced new informational messages
202!    - CM0436 - now also applies to invalid LOD settings
203!    - CM0463 - emiss_lod take precedence in case of conflict with mod_emis
204!    - CM0464 - emiss_lod valued assigned based on mode_emis if undefined
205!
206! 3930 2019-04-24 14:57:18Z forkel
207! Changed chem_non_transport_physics to chem_non_advective_processes
208!
209! 3929 2019-04-24 12:52:08Z banzhafs
210! Correct/complete module_interface introduction for chemistry model
211! Add subroutine chem_exchange_horiz_bounds
212! Bug fix deposition
213!
214! 3784 2019-03-05 14:16:20Z banzhafs
215! 2D output of emission fluxes
216!
217! 3784 2019-03-05 14:16:20Z banzhafs
218! Bugfix, uncomment erroneous commented variable used for dry deposition.
219! Bugfix in 3D emission output.
220!
221! 3784 2019-03-05 14:16:20Z banzhafs
222! Changes related to global restructuring of location messages and introduction
223! of additional debug messages
224!
225! 3784 2019-03-05 14:16:20Z banzhafs
226! some formatting of the deposition code
227!
228! 3784 2019-03-05 14:16:20Z banzhafs
229! some formatting
230!
231! 3784 2019-03-05 14:16:20Z banzhafs
232! added cs_mech to USE chem_gasphase_mod
233!
234! 3784 2019-03-05 14:16:20Z banzhafs
235! renamed get_mechanismname to get_mechanism_name
236! renamed do_emiss to emissions_anthropogenic and do_depo to deposition_dry (ecc)
237!
238! 3784 2019-03-05 14:16:20Z banzhafs
239! Unused variables removed/taken care of
240!
241! 3784 2019-03-05 14:16:20Z forkel
242! Replaced READ from unit 10 by CALL get_mechanismname also in chem_header
243!
244! 3783 2019-03-05 13:23:50Z forkel
245! Removed forgotte write statements an some unused variables (did not touch the parts related to
246! deposition)
247!
248! 3780 2019-03-05 11:19:45Z forkel
249! Removed READ from unit 10, added CALL get_mechanismname
250!
251! 3767 2019-02-27 08:18:02Z raasch
252! unused variable for file index removed from rrd-subroutines parameter list
253!
254! 3738 2019-02-12 17:00:45Z suehring
255! Clean-up debug prints
256!
257! 3737 2019-02-12 16:57:06Z suehring
258! Enable mesoscale offline nesting for chemistry variables as well as initialization of chemistry
259! via dynamic input file.
260!
261! 3719 2019-02-06 13:10:18Z kanani
262! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call to prognostic_equations
263! for better overview
264!
265! 3700 2019-01-26 17:03:42Z knoop
266! Some interface calls moved to module_interface + cleanup
267!
268! 3664 2019-01-09 14:00:37Z forkel
269! Replaced misplaced location message by @todo
270!
271! 3654 2019-01-07 16:31:57Z suehring
272! Disable misplaced location message
273!
274! 3652 2019-01-07 15:29:59Z forkel
275! Checks added for chemistry mechanism, parameter chem_mechanism added (basit)
276!
277! 2718 2018-01-02 08:49:38Z maronga
278! Initial revision
279!
280!
281!
282!
283! Authors:
284! --------
285! @author Renate Forkel
286! @author Farah Kanani-Suehring
287! @author Klaus Ketelsen
288! @author Basit Khan
289! @author Sabine Banzhaf
290!
291!
292!--------------------------------------------------------------------------------------------------!
293! Description:
294! ------------
295!> Chemistry model for PALM-4U
296!> @todo Extend chem_species type by nspec and nvar as addititional elements (RF)
297!> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF)
298!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not allowed to use the
299!>       chemistry model in a precursor run and additionally not using it in a main run
300!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
301!> @todo Separate boundary conditions for each chem spcs to be implemented
302!> @todo Currently only total concentration are calculated. Resolved, parameterized and chemistry
303!>       fluxes although partially and some completely coded but are not operational/activated in
304!>       this version. bK.
305!> @todo slight differences in passive scalar and chem spcs when chem reactions turned off. Need to
306!>       be fixed. bK
307!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
308!> @todo chemistry error messages
309!
310!--------------------------------------------------------------------------------------------------!
311
312 MODULE chemistry_model_mod
313
314    USE advec_s_pw_mod,                                                                            &
315         ONLY:  advec_s_pw
316
317    USE advec_s_up_mod,                                                                            &
318         ONLY:  advec_s_up
319
320    USE advec_ws,                                                                                  &
321         ONLY:  advec_s_ws, ws_init_flags_scalar
322
323    USE diffusion_s_mod,                                                                           &
324         ONLY:  diffusion_s
325
326    USE kinds,                                                                                     &
327         ONLY:  iwp, wp
328
329    USE indices,                                                                                   &
330         ONLY:  advc_flags_s,                                                                      &
331                nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, nzb, nzt,            &
332                wall_flags_total_0
333
334    USE pegrid,                                                                                    &
335         ONLY: myid, threads_per_task
336
337    USE bulk_cloud_model_mod,                                                                      &
338         ONLY:  bulk_cloud_model
339
340    USE control_parameters,                                                                        &
341         ONLY:  air_chemistry,                                                                     &
342                bc_dirichlet_l,                                                                    &
343                bc_dirichlet_n,                                                                    &
344                bc_dirichlet_r,                                                                    &
345                bc_dirichlet_s,                                                                    &
346                bc_radiation_l,                                                                    &
347                bc_lr_cyc,                                                                         &
348                bc_ns_cyc,                                                                         &
349                bc_radiation_n,                                                                    &
350                bc_radiation_r,                                                                    &
351                bc_radiation_s,                                                                    &
352                debug_output,                                                                      &
353                dt_3d,                                                                             &
354                humidity,                                                                          &
355                initializing_actions,                                                              &
356                intermediate_timestep_count,                                                       &
357                intermediate_timestep_count_max,                                                   &
358                max_pr_user,                                                                       &
359                message_string,                                                                    &
360                monotonic_limiter_z,                                                               &
361                omega,                                                                             &
362                restart_data_format_output,                                                        &
363                scalar_advec,                                                                      &
364                timestep_scheme,                                                                   &
365                tsc,                                                                               &
366                use_prescribed_profile_data,                                                       &
367                ws_scheme_sca
368
369    USE arrays_3d,                                                                                 &
370         ONLY:  exner, hyp, pt, q, ql, rdf_sc, tend, zu
371
372    USE chem_gasphase_mod,                                                                         &
373         ONLY:  atol, chem_gasphase_integrate, cs_mech, get_mechanism_name, nkppctrl,              &
374         nmaxfixsteps, nphot, nreact, nspec, nvar, phot_names, rtol, spc_names, t_steps, vl_dim
375
376    USE chem_modules
377
378    USE chem_photolysis_mod,                                                                       &
379        ONLY:  photolysis_control
380
381    USE cpulog,                                                                                    &
382        ONLY:  cpu_log, log_point_s
383
384    USE restart_data_mpi_io_mod,                                                                   &
385        ONLY:  rrd_mpi_io, rd_mpi_io_check_array, wrd_mpi_io
386
387    USE statistics
388
389    USE surface_mod,                                                                               &
390         ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
391
392    IMPLICIT NONE
393
394    PRIVATE
395    SAVE
396
397    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
398
399    REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
400
401    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
402                                                                            !< (e.g starting internal timestep of solver)
403
404    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1  !< pointer for swapping of timelevels for conc
405    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2  !< pointer for swapping of timelevels for conc
406    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3  !< pointer for swapping of timelevels for conc
407    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies
408                                                                            !< (only 1 timelevel required) (e.g. solver type)
409
410!
411!-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
412    INTEGER(iwp), PARAMETER ::  nlu_dep = 15            !< Number of DEPAC landuse classes (lu's)
413    INTEGER(iwp), PARAMETER ::  ncmp = 10               !< Number of DEPAC gas components
414    INTEGER(iwp), PARAMETER ::  nposp = 69              !< Number of possible species for deposition
415!
416!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1
417    INTEGER(iwp) ::  ilu_grass              = 1
418    INTEGER(iwp) ::  ilu_arable             = 2
419    INTEGER(iwp) ::  ilu_permanent_crops    = 3
420    INTEGER(iwp) ::  ilu_coniferous_forest  = 4
421    INTEGER(iwp) ::  ilu_deciduous_forest   = 5
422    INTEGER(iwp) ::  ilu_water_sea          = 6
423    INTEGER(iwp) ::  ilu_urban              = 7
424    INTEGER(iwp) ::  ilu_other              = 8
425    INTEGER(iwp) ::  ilu_desert             = 9
426    INTEGER(iwp) ::  ilu_ice                = 10
427    INTEGER(iwp) ::  ilu_savanna            = 11
428    INTEGER(iwp) ::  ilu_tropical_forest    = 12
429    INTEGER(iwp) ::  ilu_water_inland       = 13
430    INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14
431    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15
432
433!
434!-- NH3/SO2 ratio regimes:
435    INTEGER(iwp), PARAMETER ::  iratns_low      = 1       !< low ratio NH3/SO2
436    INTEGER(iwp), PARAMETER ::  iratns_high     = 2       !< high ratio NH3/SO2
437    INTEGER(iwp), PARAMETER ::  iratns_very_low = 3       !< very low ratio NH3/SO2
438!
439!-- Default:
440    INTEGER, PARAMETER ::  iratns_default = iratns_low
441!
442!-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
443    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   = (/ 0.009, 0.009, 0.009, 0.006, 0.006,    &
444                    -999.0, -999., 0.009, -999.0, -999.0, 0.009, 0.006, -999.0, 0.009, 0.008 /)*4.57
445!
446!-- Set temperatures per land use for f_temp
447    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0,  12.0,  12.0,  0.0,  0.0, -999.0,  &
448                                      -999.0, 12.0, -999.0, -999.0, 12.0,  0.0, -999.0, 12.0,  8.0/)
449    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999.0,   &
450                                     -999.0, 26.0, -999.0, -999.0, 26.0, 20.0, -999.0, 26.0, 24.0 /)
451    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999.0,   &
452                                    -999.0, 40.0, -999.0, -999.0,  40.0, 35.0, -999.0, 40.0, 39.0 /)
453!
454!-- Set f_min:
455    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999.0,     &
456                                       -999.0, 0.01, -999.0, -999.0, 0.01, 0.1, -999.0, 0.01, 0.04/)
457
458!
459!-- Set maximal conductance (m/s)
460!-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
461    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270.0, 300.0, 300.0, 140.0, 150.0,      &
462              -999.0, -999.0, 270.0, -999.0, -999.0, 270.0, 150.0, -999.0, 300.0, 422.0 /) / 41000.0
463!
464!-- Set max, min for vapour pressure deficit vpd
465    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/ 1.3, 0.9, 0.9, 0.5, 1.0, -999.0,      &
466                                           -999.0, 1.3, -999.0, -999.0, 1.3, 1.0, -999.0, 0.9, 2.8/)
467    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_min = (/ 3.0, 2.8, 2.8, 3.0, 3.25, -999.0,     &
468                                         -999.0, 3.0, -999.0, -999.0, 3.0,  3.25, -999.0, 2.8, 4.5/)
469
470    PUBLIC nreact
471    PUBLIC nspec               !< number of gas phase chemical species including constant compound (e.g. N2)
472    PUBLIC nvar                !< number of variable gas phase chemical species (nvar <= nspec)
473    PUBLIC spc_names           !< names of gas phase chemical species (come from KPP) (come from KPP)
474    PUBLIC spec_conc_2
475!
476!-- Interface section
477    INTERFACE chem_actions
478       MODULE PROCEDURE chem_actions
479       MODULE PROCEDURE chem_actions_ij
480    END INTERFACE chem_actions
481
482    INTERFACE chem_3d_data_averaging
483       MODULE PROCEDURE chem_3d_data_averaging
484    END INTERFACE chem_3d_data_averaging
485
486    INTERFACE chem_boundary_conditions
487       MODULE PROCEDURE chem_boundary_conditions
488    END INTERFACE chem_boundary_conditions
489
490    INTERFACE chem_check_data_output
491       MODULE PROCEDURE chem_check_data_output
492    END INTERFACE chem_check_data_output
493
494    INTERFACE chem_data_output_2d
495       MODULE PROCEDURE chem_data_output_2d
496    END INTERFACE chem_data_output_2d
497
498    INTERFACE chem_data_output_3d
499       MODULE PROCEDURE chem_data_output_3d
500    END INTERFACE chem_data_output_3d
501
502    INTERFACE chem_data_output_mask
503       MODULE PROCEDURE chem_data_output_mask
504    END INTERFACE chem_data_output_mask
505
506    INTERFACE chem_check_data_output_pr
507       MODULE PROCEDURE chem_check_data_output_pr
508    END INTERFACE chem_check_data_output_pr
509
510    INTERFACE chem_check_parameters
511       MODULE PROCEDURE chem_check_parameters
512    END INTERFACE chem_check_parameters
513
514    INTERFACE chem_define_netcdf_grid
515       MODULE PROCEDURE chem_define_netcdf_grid
516    END INTERFACE chem_define_netcdf_grid
517
518    INTERFACE chem_header
519       MODULE PROCEDURE chem_header
520    END INTERFACE chem_header
521
522    INTERFACE chem_init_arrays
523       MODULE PROCEDURE chem_init_arrays
524    END INTERFACE chem_init_arrays
525
526    INTERFACE chem_init
527       MODULE PROCEDURE chem_init
528    END INTERFACE chem_init
529
530    INTERFACE chem_init_profiles
531       MODULE PROCEDURE chem_init_profiles
532    END INTERFACE chem_init_profiles
533
534    INTERFACE chem_integrate
535       MODULE PROCEDURE chem_integrate_ij
536    END INTERFACE chem_integrate
537
538    INTERFACE chem_parin
539       MODULE PROCEDURE chem_parin
540    END INTERFACE chem_parin
541
542    INTERFACE chem_non_advective_processes
543       MODULE PROCEDURE chem_non_advective_processes
544       MODULE PROCEDURE chem_non_advective_processes_ij
545    END INTERFACE chem_non_advective_processes
546
547    INTERFACE chem_exchange_horiz_bounds
548       MODULE PROCEDURE chem_exchange_horiz_bounds
549    END INTERFACE chem_exchange_horiz_bounds
550
551    INTERFACE chem_prognostic_equations
552       MODULE PROCEDURE chem_prognostic_equations
553       MODULE PROCEDURE chem_prognostic_equations_ij
554    END INTERFACE chem_prognostic_equations
555
556    INTERFACE chem_rrd_local
557       MODULE PROCEDURE chem_rrd_local_ftn
558       MODULE PROCEDURE chem_rrd_local_mpi
559    END INTERFACE chem_rrd_local
560
561    INTERFACE chem_statistics
562       MODULE PROCEDURE chem_statistics
563    END INTERFACE chem_statistics
564
565    INTERFACE chem_swap_timelevel
566       MODULE PROCEDURE chem_swap_timelevel
567    END INTERFACE chem_swap_timelevel
568
569    INTERFACE chem_wrd_local
570       MODULE PROCEDURE chem_wrd_local
571    END INTERFACE chem_wrd_local
572
573    INTERFACE chem_depo
574       MODULE PROCEDURE chem_depo
575    END INTERFACE chem_depo
576
577    INTERFACE drydepos_gas_depac
578       MODULE PROCEDURE drydepos_gas_depac
579    END INTERFACE drydepos_gas_depac
580
581    INTERFACE rc_special
582       MODULE PROCEDURE rc_special
583    END INTERFACE rc_special
584
585    INTERFACE  rc_gw
586       MODULE PROCEDURE rc_gw
587    END INTERFACE rc_gw
588
589    INTERFACE rw_so2
590       MODULE PROCEDURE rw_so2
591    END INTERFACE rw_so2
592
593    INTERFACE rw_nh3_sutton
594       MODULE PROCEDURE rw_nh3_sutton
595    END INTERFACE rw_nh3_sutton
596
597    INTERFACE rw_constant
598       MODULE PROCEDURE rw_constant
599    END INTERFACE rw_constant
600
601    INTERFACE rc_gstom
602       MODULE PROCEDURE rc_gstom
603    END INTERFACE rc_gstom
604
605    INTERFACE rc_gstom_emb
606       MODULE PROCEDURE rc_gstom_emb
607    END INTERFACE rc_gstom_emb
608
609    INTERFACE par_dir_diff
610       MODULE PROCEDURE par_dir_diff
611    END INTERFACE par_dir_diff
612
613    INTERFACE rc_get_vpd
614       MODULE PROCEDURE rc_get_vpd
615    END INTERFACE rc_get_vpd
616
617    INTERFACE rc_gsoil_eff
618       MODULE PROCEDURE rc_gsoil_eff
619    END INTERFACE rc_gsoil_eff
620
621    INTERFACE rc_rinc
622       MODULE PROCEDURE rc_rinc
623    END INTERFACE rc_rinc
624
625    INTERFACE rc_rctot
626       MODULE PROCEDURE rc_rctot
627    END INTERFACE rc_rctot
628
629    INTERFACE drydepo_aero_zhang_vd
630       MODULE PROCEDURE drydepo_aero_zhang_vd
631    END INTERFACE drydepo_aero_zhang_vd
632
633    INTERFACE get_rb_cell
634       MODULE PROCEDURE  get_rb_cell
635    END INTERFACE get_rb_cell
636
637
638
639    PUBLIC chem_3d_data_averaging, chem_boundary_conditions, chem_check_data_output,               &
640           chem_check_data_output_pr, chem_check_parameters, chem_data_output_2d,                  &
641           chem_data_output_3d, chem_data_output_mask, chem_define_netcdf_grid, chem_header,       &
642           chem_init, chem_init_arrays, chem_init_profiles, chem_integrate, chem_parin,            &
643           chem_actions, chem_prognostic_equations, chem_rrd_local, chem_statistics,               &
644           chem_swap_timelevel, chem_wrd_local, chem_depo, chem_non_advective_processes,           &
645           chem_exchange_horiz_bounds
646
647 CONTAINS
648
649
650!--------------------------------------------------------------------------------------------------!
651! Description:
652! ------------
653!> Subroutine for averaging 3D data of chemical species. Due to the fact that the averaged chem
654!> arrays are allocated in chem_init, no if-query concerning the allocation is required (in any
655!> mode). Attention: If you just specify an averaged output quantity in the _p3dr file during
656!> restarts the first output includes the time between the beginning of the restart run and the
657!> first output time (not necessarily the whole averaging_interval you have specified in your
658!> _p3d/_p3dr file ).
659!--------------------------------------------------------------------------------------------------!
660 SUBROUTINE chem_3d_data_averaging( mode, variable )
661
662    USE control_parameters
663
664    USE exchange_horiz_mod,                                                                        &
665        ONLY:  exchange_horiz_2d
666
667
668    CHARACTER (LEN=*) ::  mode     !<
669    CHARACTER (LEN=*) ::  variable !<
670
671    LOGICAL ::  match_def  !< flag indicating default-type surface
672    LOGICAL ::  match_lsm  !< flag indicating natural-type surface
673    LOGICAL ::  match_usm  !< flag indicating urban-type surface
674
675    INTEGER(iwp) ::  i                  !< grid index x direction
676    INTEGER(iwp) ::  j                  !< grid index y direction
677    INTEGER(iwp) ::  k                  !< grid index z direction
678    INTEGER(iwp) ::  lsp                !< running index for chem spcs
679    INTEGER(iwp) ::  m                  !< running index surface type
680
681    IF ( ( variable(1:3) == 'kc_'  .OR.  variable(1:3) == 'em_' )  )  THEN
682
683       IF ( mode == 'allocate' )  THEN
684
685          DO  lsp = 1, nspec
686             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
687                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
688                IF ( .NOT. ALLOCATED( chem_species(lsp)%conc_av ) )  THEN
689                   ALLOCATE( chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
690                   chem_species(lsp)%conc_av = 0.0_wp
691                ENDIF
692             ENDIF
693          ENDDO
694
695       ELSEIF ( mode == 'sum' )  THEN
696
697          DO  lsp = 1, nspec
698             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
699                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
700                DO  i = nxlg, nxrg
701                   DO  j = nysg, nyng
702                      DO  k = nzb, nzt+1
703                         chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) +     &
704                                                            chem_species(lsp)%conc(k,j,i)
705                      ENDDO
706                   ENDDO
707                ENDDO
708             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
709                DO  i = nxl, nxr
710                   DO  j = nys, nyn
711                      match_def = surf_def_h(0)%start_index(j,i) <= surf_def_h(0)%end_index(j,i)
712                      match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
713                      match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
714
715                      IF ( match_def )  THEN
716                         m = surf_def_h(0)%end_index(j,i)
717                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
718                                                           surf_def_h(0)%cssws(lsp,m)
719                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
720                         m = surf_lsm_h(0)%end_index(j,i)
721                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
722                                                           surf_lsm_h(0)%cssws(lsp,m)
723                      ELSEIF ( match_usm )  THEN
724                         m = surf_usm_h(0)%end_index(j,i)
725                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) +       &
726                                                           surf_usm_h(0)%cssws(lsp,m)
727                      ENDIF
728                   ENDDO
729                ENDDO
730             ENDIF
731          ENDDO
732
733       ELSEIF ( mode == 'average' )  THEN
734
735          DO  lsp = 1, nspec
736             IF ( TRIM( variable(1:3) ) == 'kc_'  .AND.                                            &
737                  TRIM( variable(4:) )  == TRIM( chem_species(lsp)%name ) )  THEN
738                DO  i = nxlg, nxrg
739                   DO  j = nysg, nyng
740                      DO  k = nzb, nzt+1
741                         chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) /     &
742                                                            REAL( average_count_3d, KIND = wp )
743                      ENDDO
744                   ENDDO
745                ENDDO
746
747             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
748                DO  i = nxlg, nxrg
749                   DO  j = nysg, nyng
750                      chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) /          &
751                                                        REAL( average_count_3d, KIND = wp )
752                   ENDDO
753                ENDDO
754                CALL exchange_horiz_2d( chem_species(lsp)%cssws_av )
755             ENDIF
756          ENDDO
757       ENDIF
758
759    ENDIF
760
761 END SUBROUTINE chem_3d_data_averaging
762
763
764!--------------------------------------------------------------------------------------------------!
765! Description:
766! ------------
767!> Subroutine to initialize and set all boundary conditions for chemical species
768!--------------------------------------------------------------------------------------------------!
769 SUBROUTINE chem_boundary_conditions( horizontal_conditions_only )
770
771    USE arrays_3d,                                                                                 &
772        ONLY:  dzu
773
774    USE surface_mod,                                                                               &
775        ONLY:  bc_h
776
777    INTEGER(iwp) ::  i                            !< grid index x direction.
778    INTEGER(iwp) ::  j                            !< grid index y direction.
779    INTEGER(iwp) ::  k                            !< grid index z direction.
780    INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
781    INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
782    INTEGER(iwp) ::  m                            !< running index surface elements.
783
784    LOGICAL, OPTIONAL ::  horizontal_conditions_only  !< switch to set horizontal bc only
785
786
787    IF ( .NOT. PRESENT( horizontal_conditions_only ) )  THEN
788!
789!--    Boundary condtions for chemical species at horizontal walls
790       DO  lsp = 1, nspec
791
792          IF ( ibc_cs_b == 0 )  THEN
793             DO  l = 0, 1
794                !$OMP PARALLEL DO PRIVATE( i, j, k )
795                DO  m = 1, bc_h(l)%ns
796                    i = bc_h(l)%i(m)
797                    j = bc_h(l)%j(m)
798                    k = bc_h(l)%k(m)
799                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =                                  &
800                                                          chem_species(lsp)%conc(k+bc_h(l)%koff,j,i)
801                ENDDO
802             ENDDO
803
804          ELSEIF ( ibc_cs_b == 1 )  THEN
805!
806!--          In boundary_conds there is som extra loop over m here for passive tracer
807!>           TODO: clarify the meaning of the above comment. Explain in more detail or remove it. (Siggi)
808             DO  l = 0, 1
809                !$OMP PARALLEL DO PRIVATE( i, j, k )
810                DO m = 1, bc_h(l)%ns
811                   i = bc_h(l)%i(m)
812                   j = bc_h(l)%j(m)
813                   k = bc_h(l)%k(m)
814                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = chem_species(lsp)%conc_p(k,j,i)
815                ENDDO
816             ENDDO
817          ENDIF
818
819       ENDDO    ! end lsp loop
820
821!
822!--    Top boundary conditions for chemical species - Should this not be done for all species?
823!>     TODO: This question also needs to be clarified. I guess it can be removed because the loop
824!>           already runs over all species? (Siggi)
825       DO  lsp = 1, nspec
826          IF ( ibc_cs_t == 0 )  THEN
827             chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
828          ELSEIF ( ibc_cs_t == 1 )  THEN
829             chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
830          ELSEIF ( ibc_cs_t == 2 )  THEN
831             chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) +             &
832                                                   bc_cs_t_val(lsp) * dzu(nzt+1)
833          ENDIF
834       ENDDO
835
836!
837!--    Lateral boundary conditions.
838!--    Dirichlet conditions have been already set when chem_species concentration is initialized.
839!--    The initially set value is not touched during time integration, hence, this boundary value
840!--    remains at a constant value.
841!--    Neumann conditions:
842       IF ( bc_radiation_cs_s )  THEN
843          DO  lsp = 1, nspec
844             chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
845          ENDDO
846       ENDIF
847       IF ( bc_radiation_cs_n )  THEN
848          DO  lsp = 1, nspec
849             chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
850          ENDDO
851       ENDIF
852       IF ( bc_radiation_cs_l )  THEN
853          DO  lsp = 1, nspec
854             chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
855          ENDDO
856       ENDIF
857       IF ( bc_radiation_cs_r )  THEN
858          DO  lsp = 1, nspec
859             chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
860          ENDDO
861       ENDIF
862
863!--    For testing: set Dirichlet conditions for all boundaries
864!      IF ( bc_dirichlet_cs_s )  THEN
865!         DO  lsp = 1, nspec
866!            DO  i = nxlg, nxrg
867!               DO  j = nysg, nys-1
868!                  DO  k = nzb, nzt
869!                     IF ( k /= nzb )  THEN
870!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
871!                     ELSE
872!                        flag = 1.0
873!                     ENDIF
874!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
875!                  ENDDO
876!               ENDDO
877!            ENDDO
878!         ENDDO
879!      ENDIF
880!      IF ( bc_dirichlet_cs_n )  THEN
881!         DO  lsp = 1, nspec
882!            DO  i = nxlg, nxrg
883!               DO  j = nyn+1, nyng
884!                  DO  k = nzb, nzt
885!                     IF ( k /= nzb )  THEN
886!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
887!                     ELSE
888!                        flag = 1.0
889!                     ENDIF
890!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
891!                  ENDDO
892!               ENDDO
893!            ENDDO
894!         ENDDO
895!      ENDIF
896!      IF ( bc_dirichlet_cs_l )  THEN
897!         DO  lsp = 1, nspec
898!            DO  i = nxlg, nxl-1
899!               DO  j = nysg, nyng
900!                  DO  k = nzb, nzt
901!                     IF ( k /= nzb )  THEN
902!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
903!                     ELSE
904!                        flag = 1.0
905!                     ENDIF
906!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
907!                  ENDDO
908!               ENDDO
909!            ENDDO
910!         ENDDO
911!      ENDIF
912!      IF ( bc_dirichlet_cs_r )  THEN
913!         DO  lsp = 1, nspec
914!            DO  i = nxr+1, nxrg
915!               DO  j = nysg, nyng
916!                  DO  k = nzb, nzt
917!                     IF ( k /= nzb )  THEN
918!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
919!                     ELSE
920!                        flag = 1.0
921!                     ENDIF
922!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
923!                  ENDDO
924!               ENDDO
925!            ENDDO
926!         ENDDO
927!      ENDIF
928
929    ELSE
930!
931!--    Lateral Neumann booundary conditions for timelevel t.
932!--    This branch is executed when routine is called after the non-advective processes / before the
933!--    prognostic equations.
934       IF ( bc_radiation_cs_s )  THEN
935          DO  lsp = 1, nspec
936             chem_species(lsp)%conc(:,nys-1,:) = chem_species(lsp)%conc(:,nys,:)
937          ENDDO
938       ENDIF
939       IF ( bc_radiation_cs_n )  THEN
940          DO  lsp = 1, nspec
941             chem_species(lsp)%conc(:,nyn+1,:) = chem_species(lsp)%conc(:,nyn,:)
942          ENDDO
943       ENDIF
944       IF ( bc_radiation_cs_l )  THEN
945          DO  lsp = 1, nspec
946             chem_species(lsp)%conc(:,:,nxl-1) = chem_species(lsp)%conc(:,:,nxl)
947          ENDDO
948       ENDIF
949       IF ( bc_radiation_cs_r )  THEN
950          DO  lsp = 1, nspec
951             chem_species(lsp)%conc(:,:,nxr+1) = chem_species(lsp)%conc(:,:,nxr)
952          ENDDO
953       ENDIF
954
955!--    For testing: set Dirichlet conditions for all boundaries
956!      IF ( bc_dirichlet_cs_s )  THEN
957!         DO  lsp = 1, nspec
958!            DO  i = nxlg, nxrg
959!               DO  j = nysg, nys-1
960!                  DO  k = nzb, nzt
961!                     IF ( k /= nzb )  THEN
962!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
963!                     ELSE
964!                        flag = 1.0
965!                     ENDIF
966!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
967!                  ENDDO
968!               ENDDO
969!            ENDDO
970!         ENDDO
971!      ENDIF
972!      IF ( bc_dirichlet_cs_n )  THEN
973!         DO  lsp = 1, nspec
974!            DO  i = nxlg, nxrg
975!               DO  j = nyn+1, nyng
976!                  DO  k = nzb, nzt
977!                     IF ( k /= nzb )  THEN
978!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
979!                     ELSE
980!                        flag = 1.0
981!                     ENDIF
982!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
983!                  ENDDO
984!               ENDDO
985!            ENDDO
986!         ENDDO
987!      ENDIF
988!      IF ( bc_dirichlet_cs_l )  THEN
989!         DO  lsp = 1, nspec
990!            DO  i = nxlg, nxl-1
991!               DO  j = nysg, nyng
992!                  DO  k = nzb, nzt
993!                     IF ( k /= nzb )  THEN
994!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
995!                     ELSE
996!                        flag = 1.0
997!                     ENDIF
998!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
999!                  ENDDO
1000!               ENDDO
1001!            ENDDO
1002!         ENDDO
1003!      ENDIF
1004!      IF ( bc_dirichlet_cs_r )  THEN
1005!         DO  lsp = 1, nspec
1006!            DO  i = nxr+1, nxrg
1007!               DO  j = nysg, nyng
1008!                  DO  k = nzb, nzt
1009!                     IF ( k /= nzb )  THEN
1010!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
1011!                     ELSE
1012!                        flag = 1.0
1013!                     ENDIF
1014!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
1015!                  ENDDO
1016!               ENDDO
1017!            ENDDO
1018!         ENDDO
1019!      ENDIF
1020
1021    ENDIF
1022
1023 END SUBROUTINE chem_boundary_conditions
1024
1025
1026!--------------------------------------------------------------------------------------------------!
1027! Description:
1028! ------------
1029!> Subroutine for checking data output for chemical species
1030!--------------------------------------------------------------------------------------------------!
1031 SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
1032
1033
1034    CHARACTER (LEN=*) ::  unit     !<
1035    CHARACTER (LEN=*) ::  var      !<
1036
1037    INTEGER(iwp) ::  i
1038    INTEGER(iwp) ::  ilen
1039    INTEGER(iwp) ::  lsp
1040    INTEGER(iwp) ::  k
1041
1042    CHARACTER(LEN=16)    ::  spec_name
1043
1044!
1045!-- Next statement is to avoid compiler warnings about unused variables
1046    IF ( ( i + ilen + k ) > 0  .OR.  var(1:1) == ' ' )  CONTINUE
1047
1048    unit = 'illegal'
1049
1050    spec_name = TRIM( var(4:) )             !< var 1:3 is 'kc_' or 'em_'.
1051
1052    IF ( TRIM( var(1:3) ) == 'em_' )  THEN
1053       DO  lsp=1,nspec
1054          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1055             unit = 'mol m-2 s-1'
1056          ENDIF
1057!
1058!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code as passive species
1059!--       (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): set unit to micrograms per m**3 for PM10
1060!--       and PM25 (PM2.5)
1061          IF (spec_name(1:2) == 'PM')  THEN
1062             unit = 'kg m-2 s-1'
1063          ENDIF
1064       ENDDO
1065
1066    ELSE
1067
1068       DO  lsp=1,nspec
1069          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1070             unit = 'ppm'
1071          ENDIF
1072!
1073!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code as passive species
1074!--       (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): set unit to kilograms per m**3 for PM10
1075!--       and PM25 (PM2.5)
1076          IF (spec_name(1:2) == 'PM')  THEN
1077            unit = 'kg m-3'
1078          ENDIF
1079       ENDDO
1080
1081       DO  lsp=1,nphot
1082          IF ( TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
1083             unit = 'sec-1'
1084          ENDIF
1085       ENDDO
1086    ENDIF
1087
1088
1089    RETURN
1090 END SUBROUTINE chem_check_data_output
1091
1092
1093!--------------------------------------------------------------------------------------------------!
1094! Description:
1095! ------------
1096!> Subroutine for checking data output of profiles for chemistry model
1097!--------------------------------------------------------------------------------------------------!
1098 SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
1099
1100    USE arrays_3d
1101
1102    USE control_parameters,                                                                        &
1103        ONLY:  data_output_pr, message_string
1104
1105    USE profil_parameter
1106
1107    USE statistics
1108
1109
1110    CHARACTER (LEN=*)  ::  unit      !< unit
1111    CHARACTER (LEN=*)  ::  variable  !< variable name
1112    CHARACTER (LEN=*)  ::  dopr_unit !< unit
1113    CHARACTER (LEN=16) ::  spec_name !< species name extracted from output string
1114
1115    INTEGER(iwp) ::  index_start     !< start index of the species name in data-output string
1116    INTEGER(iwp) ::  index_end       !< end index of the species name in data-output string
1117    INTEGER(iwp) ::  var_count       !< number of data-output quantity
1118    INTEGER(iwp) ::  lsp             !< running index over species
1119
1120    SELECT CASE ( TRIM( variable(1:3 ) ) )
1121
1122       CASE ( 'kc_' )
1123
1124          IF ( .NOT. air_chemistry )  THEN
1125             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //          &
1126                              ' is not implemented for air_chemistry = .FALSE.'
1127             CALL message( 'chem_check_data_output_pr', 'PA0293', 1, 2, 0, 6, 0 )
1128
1129          ENDIF
1130!
1131!--       Output of total fluxes is not allowed to date.
1132          IF ( TRIM( variable(1:4) ) == 'kc_w' )  THEN
1133             IF ( TRIM( variable(1:5) ) /= 'kc_w*'  .AND.  TRIM( variable(1:5) ) /= 'kc_w"' )  THEN
1134                message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
1135                              ' is currently not implemented. Please output resolved- and '//      &
1136                              'subgrid-scale fluxes individually to obtain the total flux.'
1137                CALL message( 'chem_check_data_output_pr', 'PA0487', 1, 2, 0, 6, 0 )
1138             ENDIF
1139          ENDIF
1140!
1141!--       Check for profile output of first-order moments, i.e. variable(4:) equals a species name.
1142          spec_name = TRIM( variable(4:) )
1143          DO  lsp = 1, nspec
1144             IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1145                cs_pr_count_sp                 = cs_pr_count_sp + 1
1146                cs_pr_index_sp(cs_pr_count_sp) = lsp
1147                dopr_index(var_count)          = pr_palm + cs_pr_count_sp +                        &
1148                                                 cs_pr_count_fl_sgs + cs_pr_count_fl_res
1149                dopr_unit                      = 'ppm'
1150                IF ( spec_name(1:2) == 'PM')  THEN
1151                   dopr_unit = 'kg m-3'
1152                ENDIF
1153                hom(:,2,dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
1154                unit                              = dopr_unit
1155
1156                hom_index_spec(cs_pr_count_sp)    = dopr_index(var_count)
1157             ENDIF
1158          ENDDO
1159!
1160!--       Check for profile output of fluxes. variable(index_start:index_end) equals a species name.
1161!--       Start with SGS components.
1162          IF ( TRIM( variable(1:5) ) == 'kc_w"' )  THEN
1163             DO  lsp = 1, nspec
1164                index_end   = LEN( TRIM( variable ) ) - 1
1165                index_start = 6
1166                spec_name = TRIM( variable(index_start:index_end) )
1167                IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1168                   cs_pr_count_fl_sgs                     = cs_pr_count_fl_sgs + 1
1169                   cs_pr_index_fl_sgs(cs_pr_count_fl_sgs) = lsp
1170                   dopr_index(var_count)                  = pr_palm + cs_pr_count_sp +          &
1171                                                            cs_pr_count_fl_sgs +                &
1172                                                            cs_pr_count_fl_res
1173                   dopr_unit                              = 'm ppm s-1'
1174                   IF ( spec_name(1:2) == 'PM')  THEN
1175                      dopr_unit = 'kg m-2 s-1'
1176                   ENDIF
1177                   hom(:,2,dopr_index(var_count),:)     = SPREAD( zu, 2, statistic_regions+1 )
1178                   unit                                 = dopr_unit
1179
1180                   hom_index_fl_sgs(cs_pr_count_fl_sgs) = dopr_index(var_count)
1181                ENDIF
1182             ENDDO
1183          ENDIF
1184!
1185!--       Proceed with resolved-scale fluxes.
1186          IF ( TRIM( variable(1:5) ) == 'kc_w*' )  THEN
1187             spec_name = TRIM( variable(6:) )
1188             DO  lsp = 1, nspec
1189                index_start = 6
1190                index_end   = LEN( TRIM( variable ) ) - 1
1191                spec_name = TRIM( variable(index_start:index_end) )
1192                IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1193                   cs_pr_count_fl_res                     = cs_pr_count_fl_res + 1
1194                   cs_pr_index_fl_res(cs_pr_count_fl_res) = lsp
1195                   dopr_index(var_count)                  = pr_palm + cs_pr_count_sp +             &
1196                                                            cs_pr_count_fl_sgs +                   &
1197                                                            cs_pr_count_fl_res
1198                   dopr_unit                              = 'm ppm s-1'
1199                   IF ( spec_name(1:2) == 'PM')  THEN
1200                      dopr_unit = 'kg m-2 s-1'
1201                   ENDIF
1202                   hom(:,2, dopr_index(var_count),:)    = SPREAD( zu, 2, statistic_regions+1 )
1203                   unit                                 = dopr_unit
1204
1205                   hom_index_fl_res(cs_pr_count_fl_res) = dopr_index(var_count)
1206                ENDIF
1207             ENDDO
1208          ENDIF
1209       CASE DEFAULT
1210          unit = 'illegal'
1211
1212    END SELECT
1213
1214 END SUBROUTINE chem_check_data_output_pr
1215
1216
1217!--------------------------------------------------------------------------------------------------!
1218! Description:
1219! ------------
1220!> Check parameters routine for chemistry_model_mod
1221!--------------------------------------------------------------------------------------------------!
1222 SUBROUTINE chem_check_parameters
1223
1224    USE control_parameters,                                                                        &
1225        ONLY:  bc_lr, bc_ns
1226
1227    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
1228    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
1229
1230    LOGICAL  ::  found
1231
1232!
1233!-- Check for chemical reactions status
1234    IF ( chem_gasphase_on )  THEN
1235       message_string = 'Chemical reactions: ON'
1236       CALL message( 'chem_check_parameters', 'PA0517', 0, 0, 0, 6, 0 )
1237    ELSEIF ( .NOT. ( chem_gasphase_on ) )  THEN
1238       message_string = 'Chemical reactions: OFF'
1239       CALL message( 'chem_check_parameters', 'PA0517', 0, 0, 0, 6, 0 )
1240    ENDIF
1241!
1242!-- Check for chemistry time-step
1243    IF ( call_chem_at_all_substeps )  THEN
1244       message_string = 'Chemistry is calculated at all meteorology time-step'
1245       CALL message( 'chem_check_parameters', 'PA0522', 0, 0, 0, 6, 0 )
1246    ELSEIF ( .NOT. ( call_chem_at_all_substeps ) )  THEN
1247       message_string = 'Sub-time-steps are skipped for chemistry time-steps'
1248       CALL message( 'chem_check_parameters', 'PA0523', 0, 0, 0, 6, 0 )
1249    ENDIF
1250!
1251!-- Check for photolysis scheme
1252    IF ( ( photolysis_scheme /= 'simple' ) .AND. ( photolysis_scheme /= 'constant' )  )  THEN
1253       message_string = 'Incorrect photolysis scheme selected, please check spelling'
1254       CALL message( 'chem_check_parameters', 'PA0593', 1, 2, 0, 6, 0 )
1255    ENDIF
1256!
1257!-- Check for chemical mechanism used
1258    CALL get_mechanism_name
1259    IF ( chem_mechanism /= TRIM( cs_mech ) )  THEN
1260       message_string =                                                                            &
1261       'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
1262       CALL message( 'chem_check_parameters', 'PA0675', 1, 2, 0, 6, 0 )
1263    ENDIF
1264
1265!
1266!-- Check bottom boundary condition and set internal steering parameter
1267    IF ( bc_cs_b == 'dirichlet' )  THEN
1268       ibc_cs_b = 0
1269    ELSEIF ( bc_cs_b == 'neumann' )  THEN
1270       ibc_cs_b = 1
1271    ELSE
1272       message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'
1273       CALL message( 'chem_check_parameters', 'PA0730', 1, 2, 0, 6, 0 )
1274    ENDIF
1275!
1276!-- Check top boundary condition and set internal steering parameter
1277    IF ( bc_cs_t == 'dirichlet' )  THEN
1278       ibc_cs_t = 0
1279    ELSEIF ( bc_cs_t == 'neumann' )  THEN
1280       ibc_cs_t = 1
1281    ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
1282       ibc_cs_t = 2
1283    ELSEIF ( bc_cs_t == 'nested' )  THEN
1284       ibc_cs_t = 3
1285    ELSE
1286       message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'
1287       CALL message( 'chem_check_parameters', 'PA0731', 1, 2, 0, 6, 0 )
1288    ENDIF
1289
1290!
1291!-- If nesting_chem = .F., set top boundary condition to its default value
1292    IF ( .NOT. nesting_chem  .AND.  ibc_cs_t == 3  )  THEN
1293       ibc_cs_t = 2
1294       bc_cs_t = 'initial_gradient'
1295    ENDIF
1296
1297!
1298!-- Check left and right boundary conditions. First set default value if not set by user.
1299    IF ( bc_cs_l == 'undefined' )  THEN
1300       IF ( bc_lr == 'cyclic' )  THEN
1301          bc_cs_l = 'cyclic'
1302       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
1303          bc_cs_l = 'dirichlet'
1304       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1305          bc_cs_l = 'neumann'
1306       ENDIF
1307    ENDIF
1308    IF ( bc_cs_r == 'undefined' )  THEN
1309       IF ( bc_lr == 'cyclic' )  THEN
1310          bc_cs_r = 'cyclic'
1311       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
1312          bc_cs_r = 'neumann'
1313       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1314          bc_cs_r = 'dirichlet'
1315       ENDIF
1316    ENDIF
1317    IF ( bc_cs_l /= 'dirichlet'  .AND.  bc_cs_l /= 'neumann'  .AND.  bc_cs_l /= 'cyclic' )  THEN
1318       message_string = 'unknown boundary condition: bc_cs_l = "' // TRIM( bc_cs_l ) // '"'
1319       CALL message( 'chem_check_parameters','PA0505', 1, 2, 0, 6, 0 )
1320    ENDIF
1321    IF ( bc_cs_r /= 'dirichlet'  .AND.  bc_cs_r /= 'neumann'  .AND.  bc_cs_r /= 'cyclic' )  THEN
1322       message_string = 'unknown boundary condition: bc_cs_r = "' // TRIM( bc_cs_r ) // '"'
1323       CALL message( 'chem_check_parameters','PA0505', 1, 2, 0, 6, 0 )
1324    ENDIF
1325!
1326!-- Check north and south boundary conditions. First set default value if not set by user.
1327    IF ( bc_cs_n == 'undefined' )  THEN
1328       IF ( bc_ns == 'cyclic' )  THEN
1329          bc_cs_n = 'cyclic'
1330       ELSEIF ( bc_ns == 'dirichlet/radiation' )  THEN
1331          bc_cs_n = 'dirichlet'
1332       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1333          bc_cs_n = 'neumann'
1334       ENDIF
1335    ENDIF
1336    IF ( bc_cs_s == 'undefined' )  THEN
1337       IF ( bc_ns == 'cyclic' )  THEN
1338          bc_cs_s = 'cyclic'
1339       ELSEIF ( bc_ns == 'dirichlet/radiation' )  THEN
1340          bc_cs_s = 'neumann'
1341       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1342          bc_cs_s = 'dirichlet'
1343       ENDIF
1344    ENDIF
1345    IF ( bc_cs_n /= 'dirichlet'  .AND.  bc_cs_n /= 'neumann'  .AND.  bc_cs_n /= 'cyclic' )  THEN
1346       message_string = 'unknown boundary condition: bc_cs_n = "' // TRIM( bc_cs_n ) // '"'
1347       CALL message( 'chem_check_parameters','PA0505', 1, 2, 0, 6, 0 )
1348    ENDIF
1349    IF ( bc_cs_s /= 'dirichlet'  .AND.  bc_cs_s /= 'neumann'  .AND.  bc_cs_s /= 'cyclic' )  THEN
1350       message_string = 'unknown boundary condition: bc_cs_s = "' // TRIM( bc_cs_s ) // '"'
1351       CALL message( 'chem_check_parameters','PA0505', 1, 2, 0, 6, 0 )
1352    ENDIF
1353!
1354!-- Cyclic conditions must be set identically at opposing boundaries
1355    IF ( ( bc_cs_l == 'cyclic' .AND. bc_cs_r /= 'cyclic' )  .OR.                                   &
1356         ( bc_cs_r == 'cyclic' .AND. bc_cs_l /= 'cyclic' ) )  THEN
1357       message_string = 'boundary conditions bc_cs_l and bc_cs_r must both be cyclic or non-cyclic'
1358       CALL message( 'chem_check_parameters','PA0716', 1, 2, 0, 6, 0 )
1359    ENDIF
1360    IF ( ( bc_cs_n == 'cyclic' .AND. bc_cs_s /= 'cyclic' )  .OR.                                   &
1361         ( bc_cs_s == 'cyclic' .AND. bc_cs_n /= 'cyclic' ) )  THEN
1362       message_string = 'boundary conditions bc_cs_n and bc_cs_s must both be cyclic or non-cyclic'
1363       CALL message( 'chem_check_parameters','PA0716', 1, 2, 0, 6, 0 )
1364    ENDIF
1365!
1366!-- Set the switches that control application of horizontal boundary conditions at the boundaries
1367!-- of the total domain
1368    IF ( bc_cs_n == 'dirichlet'  .AND.  nyn == ny )  bc_dirichlet_cs_n = .TRUE.
1369    IF ( bc_cs_n == 'neumann'    .AND.  nyn == ny )  bc_radiation_cs_n = .TRUE.
1370    IF ( bc_cs_s == 'dirichlet'  .AND.  nys ==  0 )  bc_dirichlet_cs_s = .TRUE.
1371    IF ( bc_cs_s == 'neumann'    .AND.  nys ==  0 )  bc_radiation_cs_s = .TRUE.
1372    IF ( bc_cs_l == 'dirichlet'  .AND.  nxl ==  0 )  bc_dirichlet_cs_l = .TRUE.
1373    IF ( bc_cs_l == 'neumann'    .AND.  nxl ==  0 )  bc_radiation_cs_l = .TRUE.
1374    IF ( bc_cs_r == 'dirichlet'  .AND.  nxr == nx )  bc_dirichlet_cs_r = .TRUE.
1375    IF ( bc_cs_r == 'neumann'    .AND.  nxr == nx )  bc_radiation_cs_r = .TRUE.
1376
1377!
1378!-- Set the communicator to be used for ghost layer data exchange
1379!-- 1: cyclic, 2: cyclic along x, 3: cyclic along y, 4: non-cyclic
1380    IF ( bc_cs_l == 'cyclic' )  THEN
1381       IF ( bc_cs_s == 'cyclic' )  THEN
1382          communicator_chem = 1
1383       ELSE
1384          communicator_chem = 2
1385       ENDIF
1386    ELSE
1387       IF ( bc_cs_s == 'cyclic' )  THEN
1388          communicator_chem = 3
1389       ELSE
1390          communicator_chem = 4
1391       ENDIF
1392    ENDIF
1393
1394!
1395!-- chem_check_parameters is called before the array chem_species is allocated!
1396!-- temporary switch of this part of the check
1397!    RETURN                !bK commented
1398!> TODO: this workaround definitely needs to be removed from here!!!
1399    CALL chem_init_internal
1400!
1401!-- Check for initial chem species input
1402    lsp_usr = 1
1403    lsp     = 1
1404    DO WHILE ( cs_name (lsp_usr) /= 'novalue')
1405       found = .FALSE.
1406       DO  lsp = 1, nvar
1407          IF ( TRIM( cs_name (lsp_usr) ) == TRIM( chem_species(lsp)%name) )  THEN
1408             found = .TRUE.
1409             EXIT
1410          ENDIF
1411       ENDDO
1412       IF ( .NOT.  found )  THEN
1413          message_string = 'Unused/incorrect input for initial surface value: ' //     &
1414                            TRIM( cs_name(lsp_usr) )
1415          CALL message( 'chem_check_parameters', 'PA0732', 1, 2, 0, 6, 0 )
1416       ENDIF
1417       lsp_usr = lsp_usr + 1
1418    ENDDO
1419!
1420!-- Check for surface  emission flux chem species
1421    lsp_usr = 1
1422    lsp     = 1
1423    DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
1424       found = .FALSE.
1425       DO  lsp = 1, nvar
1426          IF ( TRIM( surface_csflux_name (lsp_usr) ) == TRIM( chem_species(lsp)%name ) )  THEN
1427             found = .TRUE.
1428             EXIT
1429          ENDIF
1430       ENDDO
1431       IF ( .NOT.  found )  THEN
1432          message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: '  &
1433                            // TRIM( surface_csflux_name(lsp_usr) )
1434          CALL message( 'chem_check_parameters', 'PA0733', 1, 2, 0, 6, 0 )
1435       ENDIF
1436       lsp_usr = lsp_usr + 1
1437    ENDDO
1438
1439 END SUBROUTINE chem_check_parameters
1440
1441
1442!--------------------------------------------------------------------------------------------------!
1443! Description:
1444! ------------
1445!> Subroutine defining 2D output variables for chemical species
1446!> @todo: Remove "mode" from argument list, not used.
1447!--------------------------------------------------------------------------------------------------!
1448 SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do, &
1449                                 fill_value )
1450
1451
1452    CHARACTER (LEN=*) ::  grid       !<
1453    CHARACTER (LEN=*) ::  mode       !<
1454    CHARACTER (LEN=*) ::  variable   !<
1455
1456    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or
1457                                     !< time-averaged data
1458    INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
1459    INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
1460
1461    LOGICAL      ::  found           !<
1462    LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross
1463                                     !< sections)
1464
1465    REAL(wp)     ::  fill_value
1466
1467    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf
1468
1469!
1470!-- local variables.
1471    CHARACTER(LEN=16)    ::  spec_name
1472    INTEGER(iwp) ::  lsp
1473    INTEGER(iwp) ::  i               !< grid index along x-direction
1474    INTEGER(iwp) ::  j               !< grid index along y-direction
1475    INTEGER(iwp) ::  k               !< grid index along z-direction
1476    INTEGER(iwp) ::  m               !< running indices for surfaces
1477    INTEGER(iwp) ::  char_len        !< length of a character string
1478!
1479!-- Next statement is to avoid compiler warnings about unused variables
1480    IF ( mode(1:1) == ' '  .OR.  two_d )  CONTINUE
1481
1482    found = .FALSE.
1483    char_len  = LEN_TRIM( variable )
1484
1485    spec_name = TRIM( variable(4:char_len-3) )
1486!
1487!-- Output of emission values, i.e. surface fluxes cssws.
1488    IF ( variable(1:3) == 'em_' )  THEN
1489
1490       local_pf = 0.0_wp
1491
1492       DO  lsp = 1, nvar
1493          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1494!
1495!--          No average output for now.
1496!--          !!! IT NEEDS TO RETHINK - with fully 3D structure, only lower (upper)
1497!--          !!! upward facing horizontal surfaces should be taken into account here
1498             DO  m = 1, surf_lsm_h(0)%ns
1499                local_pf(surf_lsm_h(0)%i(m),surf_lsm_h(0)%j(m),nzb+1) =                                  &
1500                                                   local_pf(surf_lsm_h(0)%i(m),surf_lsm_h(0)%j(m),nzb+1) &
1501                                                   + surf_lsm_h(0)%cssws(lsp,m)
1502             ENDDO
1503             DO  m = 1, surf_usm_h(0)%ns
1504                   local_pf(surf_usm_h(0)%i(m),surf_usm_h(0)%j(m),nzb+1) =                               &
1505                                                   local_pf(surf_usm_h(0)%i(m),surf_usm_h(0)%j(m),nzb+1) &
1506                                                   + surf_usm_h(0)%cssws(lsp,m)
1507             ENDDO
1508             grid = 'zu'
1509             found = .TRUE.
1510          ENDIF
1511       ENDDO
1512
1513    ELSE
1514
1515       DO  lsp=1,nspec
1516          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                          &
1517                ( (variable(char_len-2:) == '_xy')  .OR.                                           &
1518                  (variable(char_len-2:) == '_xz')  .OR.                                           &
1519                  (variable(char_len-2:) == '_yz') ) )  THEN
1520!
1521!--   todo: remove or replace by "CALL message" mechanism (kanani)
1522!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM( variable )  //       &
1523!                                                             TRIM( chem_species(lsp)%name )
1524             IF (av == 0)  THEN
1525                DO  i = nxl, nxr
1526                   DO  j = nys, nyn
1527                      DO  k = nzb_do, nzt_do
1528                           local_pf(i,j,k) = MERGE(                                                &
1529                                              chem_species(lsp)%conc(k,j,i),                       &
1530                                              REAL( fill_value, KIND = wp ),                       &
1531                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
1532                      ENDDO
1533                   ENDDO
1534                ENDDO
1535
1536             ELSE
1537                DO  i = nxl, nxr
1538                   DO  j = nys, nyn
1539                      DO  k = nzb_do, nzt_do
1540                           local_pf(i,j,k) = MERGE(                                                &
1541                                              chem_species(lsp)%conc_av(k,j,i),                    &
1542                                              REAL( fill_value, KIND = wp ),                       &
1543                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
1544                      ENDDO
1545                   ENDDO
1546                ENDDO
1547             ENDIF
1548             grid = 'zu'
1549             found = .TRUE.
1550          ENDIF
1551       ENDDO
1552    ENDIF
1553
1554    RETURN
1555
1556 END SUBROUTINE chem_data_output_2d
1557
1558
1559!--------------------------------------------------------------------------------------------------!
1560! Description:
1561! ------------
1562!> Subroutine defining 3D output variables for chemical species
1563!--------------------------------------------------------------------------------------------------!
1564 SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1565
1566
1567    USE surface_mod
1568
1569    CHARACTER (LEN=*)    ::  variable     !<
1570
1571    INTEGER(iwp)         ::  av         !<
1572    INTEGER(iwp) ::  nzb_do             !< lower limit of the data output (usually 0)
1573    INTEGER(iwp) ::  nzt_do             !< vertical upper limit of the data output (usually nz_do3d)
1574
1575    LOGICAL      ::  found                !<
1576
1577    REAL(wp)             ::  fill_value   !<
1578
1579    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf
1580!
1581!-- Local variables
1582    CHARACTER(LEN=16)    ::  spec_name
1583
1584    INTEGER(iwp)         ::  i
1585    INTEGER(iwp)         ::  j
1586    INTEGER(iwp)         ::  k
1587    INTEGER(iwp)         ::  m       !< running indices for surfaces
1588    INTEGER(iwp)         ::  l
1589    INTEGER(iwp)         ::  lsp     !< running index for chem spcs
1590
1591
1592    found = .FALSE.
1593    IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) )  THEN
1594       RETURN
1595    ENDIF
1596
1597    spec_name = TRIM( variable(4:) )
1598
1599    IF ( variable(1:3) == 'em_' )  THEN
1600
1601       DO  lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
1602          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1603
1604             local_pf = 0.0_wp
1605!
1606!--          no average for now
1607             DO l = 0, 1
1608                DO  m = 1, surf_usm_h(l)%ns
1609                   local_pf(surf_usm_h(l)%i(m),surf_usm_h(l)%j(m),surf_usm_h(l)%k(m)) =            &
1610                                local_pf(surf_usm_h(l)%i(m),surf_usm_h(l)%j(m),surf_usm_h(l)%k(m)) &
1611                                + surf_usm_h(l)%cssws(lsp,m)
1612                ENDDO
1613                DO  m = 1, surf_lsm_h(l)%ns
1614                   local_pf(surf_lsm_h(l)%i(m),surf_lsm_h(l)%j(m),surf_lsm_h(l)%k(m)) =            &
1615                                local_pf(surf_lsm_h(l)%i(m),surf_lsm_h(l)%j(m),surf_lsm_h(l)%k(m)) &
1616                                + surf_lsm_h(l)%cssws(lsp,m)
1617                ENDDO
1618             ENDDO
1619             DO  l = 0, 3
1620                DO  m = 1, surf_usm_v(l)%ns
1621                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) =            &
1622                                local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) &
1623                                + surf_usm_v(l)%cssws(lsp,m)
1624                ENDDO
1625                DO  m = 1, surf_lsm_v(l)%ns
1626                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) =            &
1627                                local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) &
1628                                + surf_lsm_v(l)%cssws(lsp,m)
1629                ENDDO
1630             ENDDO
1631             found = .TRUE.
1632          ENDIF
1633       ENDDO
1634    ELSE
1635      DO  lsp = 1, nspec
1636         IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1637!
1638!--   todo: remove or replace by "CALL message" mechanism (kanani)
1639!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
1640!                                                           TRIM( chem_species(lsp)%name )
1641            IF (av == 0)  THEN
1642               DO  i = nxl, nxr
1643                  DO  j = nys, nyn
1644                     DO  k = nzb_do, nzt_do
1645                         local_pf(i,j,k) = MERGE(                                                  &
1646                                             chem_species(lsp)%conc(k,j,i),                        &
1647                                             REAL( fill_value, KIND = wp ),                        &
1648                                             BTEST( wall_flags_total_0(k,j,i), 0 ) )
1649                     ENDDO
1650                  ENDDO
1651               ENDDO
1652
1653            ELSE
1654
1655               DO  i = nxl, nxr
1656                  DO  j = nys, nyn
1657                     DO  k = nzb_do, nzt_do
1658                         local_pf(i,j,k) = MERGE(                                                  &
1659                                             chem_species(lsp)%conc_av(k,j,i),                     &
1660                                             REAL( fill_value, KIND = wp ),                        &
1661                                             BTEST( wall_flags_total_0(k,j,i), 0 ) )
1662                     ENDDO
1663                  ENDDO
1664               ENDDO
1665            ENDIF
1666            found = .TRUE.
1667         ENDIF
1668      ENDDO
1669    ENDIF
1670
1671    RETURN
1672
1673 END SUBROUTINE chem_data_output_3d
1674
1675
1676!--------------------------------------------------------------------------------------------------!
1677! Description:
1678! ------------
1679!> Subroutine defining mask output variables for chemical species
1680!--------------------------------------------------------------------------------------------------!
1681 SUBROUTINE chem_data_output_mask( av, variable, found, local_pf, mid )
1682
1683
1684    USE control_parameters
1685
1686    REAL(wp), PARAMETER ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
1687
1688    CHARACTER(LEN=16) ::  spec_name
1689    CHARACTER(LEN=*)  ::  variable    !<
1690
1691    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or
1692                                     !< time-averaged data
1693    INTEGER(iwp) ::  i               !< grid index along x-direction
1694    INTEGER(iwp) ::  im              !< loop index for masked variables
1695    INTEGER(iwp) ::  j               !< grid index along y-direction
1696    INTEGER(iwp) ::  jm              !< loop index for masked variables
1697    INTEGER(iwp) ::  k               !< grid index along z-direction
1698    INTEGER(iwp) ::  kk              !< masked output index along z-direction
1699    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
1700    INTEGER(iwp) ::  lsp
1701    INTEGER(iwp) ::  mid             !< masked output running index
1702
1703    LOGICAL ::  found
1704
1705    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !<
1706
1707
1708!
1709!-- Local variables.
1710
1711    spec_name = TRIM( variable(4:) )
1712    found = .FALSE.
1713
1714    DO  lsp=1,nspec
1715       IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1716!
1717!-- todo: remove or replace by "CALL message" mechanism (kanani)
1718!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
1719!                                                        TRIM( chem_species(lsp)%name )
1720          IF (av == 0)  THEN
1721             IF ( .NOT. mask_surface(mid) )  THEN
1722
1723                DO  i = 1, mask_size_l(mid,1)
1724                   DO  j = 1, mask_size_l(mid,2)
1725                      DO  k = 1, mask_size(mid,3)
1726                          local_pf(i,j,k) = chem_species(lsp)%conc( mask_k(mid,k),                 &
1727                                                                    mask_j(mid,j),                 &
1728                                                                    mask_i(mid,i)      )
1729                      ENDDO
1730                   ENDDO
1731                ENDDO
1732
1733             ELSE
1734!
1735!--             Terrain-following masked output
1736                DO  i = 1, mask_size_l(mid,1)
1737                   DO  j = 1, mask_size_l(mid,2)
1738!
1739!--                   Get k index of the highest terraing surface
1740                      im = mask_i(mid,i)
1741                      jm = mask_j(mid,j)
1742                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ),        &
1743                                    DIM = 1 ) - 1
1744                      DO  k = 1, mask_size_l(mid,3)
1745                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
1746!
1747!--                      Set value if not in building
1748                         IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
1749                            local_pf(i,j,k) = fill_value
1750                         ELSE
1751                            local_pf(i,j,k) = chem_species(lsp)%conc(kk,jm,im)
1752                         ENDIF
1753                      ENDDO
1754                   ENDDO
1755                ENDDO
1756
1757             ENDIF
1758          ELSE
1759             IF ( .NOT. mask_surface(mid) )  THEN
1760
1761                DO  i = 1, mask_size_l(mid,1)
1762                   DO  j = 1, mask_size_l(mid,2)
1763                      DO  k =  1, mask_size_l(mid,3)
1764                          local_pf(i,j,k) = chem_species(lsp)%conc_av(  mask_k(mid,k),             &
1765                                                                        mask_j(mid,j),             &
1766                                                                        mask_i(mid,i)         )
1767                      ENDDO
1768                   ENDDO
1769                ENDDO
1770
1771             ELSE
1772!
1773!--             Terrain-following masked output
1774                DO  i = 1, mask_size_l(mid,1)
1775                   DO  j = 1, mask_size_l(mid,2)
1776!
1777!--                   Get k index of the highest terraing surface
1778                      im = mask_i(mid,i)
1779                      jm = mask_j(mid,j)
1780                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )),         &
1781                                    DIM = 1 ) - 1
1782                      DO  k = 1, mask_size_l(mid,3)
1783                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
1784!
1785!--                      Set value if not in building
1786                         IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
1787                            local_pf(i,j,k) = fill_value
1788                         ELSE
1789                            local_pf(i,j,k) = chem_species(lsp)%conc_av(kk,jm,im)
1790                         ENDIF
1791                      ENDDO
1792                   ENDDO
1793                ENDDO
1794
1795             ENDIF
1796
1797          ENDIF
1798          found = .TRUE.
1799          EXIT
1800       ENDIF
1801    ENDDO
1802
1803    RETURN
1804
1805 END SUBROUTINE chem_data_output_mask
1806
1807
1808!--------------------------------------------------------------------------------------------------!
1809! Description:
1810! ------------
1811!> Subroutine defining appropriate grid for netcdf variables.
1812!> It is called out from subroutine netcdf.
1813!--------------------------------------------------------------------------------------------------!
1814 SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1815
1816
1817    CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
1818
1819    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
1820    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
1821    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
1822
1823    LOGICAL, INTENT(OUT)           ::  found        !<
1824
1825    found  = .TRUE.
1826
1827    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                !< always the same grid for
1828                                                                         !< chemistry variables
1829       grid_x = 'x'
1830       grid_y = 'y'
1831       grid_z = 'zu'
1832    ELSE
1833       found  = .FALSE.
1834       grid_x = 'none'
1835       grid_y = 'none'
1836       grid_z = 'none'
1837    ENDIF
1838
1839
1840 END SUBROUTINE chem_define_netcdf_grid
1841
1842
1843!--------------------------------------------------------------------------------------------------!
1844! Description:
1845! ------------
1846!> Subroutine defining header output for chemistry model
1847!--------------------------------------------------------------------------------------------------!
1848 SUBROUTINE chem_header( io )
1849
1850    CHARACTER (LEN=80)  :: docsflux_chr
1851    CHARACTER (LEN=80)  :: docsinit_chr
1852
1853    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1854
1855    INTEGER(iwp)  :: cs_fixed
1856    INTEGER(iwp)  :: lsp                       !< running index for chem spcs
1857
1858!
1859!-- Get name of chemical mechanism from chem_gasphase_mod
1860    CALL get_mechanism_name
1861!
1862!-- Write chemistry model  header
1863    WRITE( io, 1 )
1864!
1865!-- Gasphase reaction status
1866    IF ( chem_gasphase_on )  THEN
1867       WRITE( io, 2 )
1868    ELSE
1869       WRITE( io, 3 )
1870    ENDIF
1871!
1872!-- Chemistry time-step
1873    WRITE ( io, 4 ) cs_time_step
1874!
1875!-- Emission mode info
1876!-- At the moment the evaluation is done with both emiss_lod and mode_emis but once salsa has been
1877!-- migrated to emiss_lod the .OR. mode_emis conditions can be removed (ecc 20190513)
1878    IF     ( ( emiss_lod == 1 )  .OR.  ( mode_emis == 'DEFAULT' ) )        THEN
1879        WRITE ( io, 5 )
1880    ELSEIF ( ( emiss_lod == 0 )  .OR.  ( mode_emis == 'PARAMETERIZED' ) )  THEN
1881        WRITE ( io, 6 )
1882    ELSEIF ( ( emiss_lod == 2 )  .OR.  ( mode_emis == 'PRE-PROCESSED' ) )  THEN
1883        WRITE ( io, 7 )
1884    ENDIF
1885!
1886!-- Photolysis scheme info
1887    IF ( photolysis_scheme == "simple" )  THEN
1888       WRITE( io, 8 )
1889    ELSEIF (photolysis_scheme == "constant" )  THEN
1890       WRITE( io, 9 )
1891    ENDIF
1892!
1893!-- Emission flux info
1894    lsp = 1
1895    docsflux_chr ='Chemical species for surface emission flux: '
1896    DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
1897       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','
1898       IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
1899          WRITE ( io, 10 )  docsflux_chr
1900          docsflux_chr = '       '
1901       ENDIF
1902       lsp = lsp + 1
1903    ENDDO
1904
1905    IF ( docsflux_chr /= '' )  THEN
1906       WRITE ( io, 10 )  docsflux_chr
1907    ENDIF
1908!
1909!-- Initialization of Surface and profile chemical species
1910    lsp = 1
1911    docsinit_chr ='Chemical species for initial surface and profile emissions: '
1912    DO WHILE ( cs_name(lsp) /= 'novalue' )
1913       docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ','
1914       IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
1915          WRITE ( io, 11 )  docsinit_chr
1916          docsinit_chr = '       '
1917       ENDIF
1918       lsp = lsp + 1
1919    ENDDO
1920
1921    IF ( docsinit_chr /= '' )  THEN
1922       WRITE ( io, 11 )  docsinit_chr
1923    ENDIF
1924
1925    IF ( nesting_chem )  WRITE( io, 12 )  nesting_chem
1926    IF ( nesting_offline_chem )  WRITE( io, 13 )  nesting_offline_chem
1927
1928    WRITE( io, 14 )  TRIM( bc_cs_b ), TRIM( bc_cs_t ), TRIM( bc_cs_s ), TRIM( bc_cs_n ),           &
1929                     TRIM( bc_cs_l ), TRIM( bc_cs_r )
1930
1931!
1932!-- Number of variable and fix chemical species and number of reactions
1933    cs_fixed = nspec - nvar
1934    WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech
1935    WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
1936    WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
1937    WRITE ( io, * ) '   --> Total number of reactions : ', nreact
1938
1939
19401   FORMAT (//' Chemistry model information:'/' ----------------------------'/)
19412   FORMAT ('    --> Chemical reactions are turned on')
19423   FORMAT ('    --> Chemical reactions are turned off')
19434   FORMAT ('    --> Time-step for chemical species: ',F6.2, ' s')
19445   FORMAT ('    --> Emission mode = DEFAULT ')
19456   FORMAT ('    --> Emission mode = PARAMETERIZED ')
19467   FORMAT ('    --> Emission mode = PRE-PROCESSED ')
19478   FORMAT ('    --> Photolysis scheme used =  simple ')
19489   FORMAT ('    --> Photolysis scheme used =  constant ')
194910  FORMAT (/'    ',A)
195011  FORMAT (/'    ',A)
195112  FORMAT (/'   Nesting for chemistry variables: ', L1 )
195213  FORMAT (/'   Offline nesting for chemistry variables: ', L1 )
195314  FORMAT (/'   Boundary conditions for chemical species:', /                                     &
1954             '      bottom/top:   ',A10,' / ',A10, /                                               &
1955             '      north/south:  ',A10,' / ',A10, /                                               &
1956             '      left/right:   ',A10,' / ',A10)
1957
1958 END SUBROUTINE chem_header
1959
1960
1961!--------------------------------------------------------------------------------------------------!
1962! Description:
1963! ------------
1964!> Subroutine initializating chemistry_model_mod specific arrays
1965!--------------------------------------------------------------------------------------------------!
1966 SUBROUTINE chem_init_arrays
1967!
1968!-- Please use this place to allocate required arrays
1969
1970 END SUBROUTINE chem_init_arrays
1971
1972
1973!--------------------------------------------------------------------------------------------------!
1974! Description:
1975! ------------
1976!> Subroutine initializating chemistry_model_mod
1977!--------------------------------------------------------------------------------------------------!
1978 SUBROUTINE chem_init
1979
1980!
1981!-- 20200203 (ECC)
1982!-- introduced additional interfaces for on-demand emission update
1983
1984!    USE chem_emissions_mod,                                                                        &
1985!        ONLY:  chem_emissions_init
1986
1987    USE chem_emissions_mod,                                                                        &
1988        ONLY:  chem_emissions_header_init, chem_emissions_init
1989
1990    USE netcdf_data_input_mod,                                                                     &
1991        ONLY:  init_3d
1992
1993
1994    INTEGER(iwp) ::  i !< running index x dimension
1995    INTEGER(iwp) ::  j !< running index y dimension
1996    INTEGER(iwp) ::  n !< running index for chemical species
1997
1998
1999    IF ( debug_output )  CALL debug_message( 'chem_init', 'start' )
2000!
2001!-- Next statement is to avoid compiler warning about unused variables
2002    IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub +     &
2003           ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest +        &
2004           ilu_urban ) == 0 )  CONTINUE
2005
2006!
2007!-- 20200203 (ECC)
2008!-- Calls specific emisisons initialization subroutines for legacy mode and on-demand mode
2009
2010!    IF ( emissions_anthropogenic )  CALL chem_emissions_init
2011
2012    IF  ( emissions_anthropogenic )  THEN
2013
2014       IF  ( emiss_read_legacy_mode )  THEN
2015          CALL chem_emissions_init
2016       ELSE
2017          CALL chem_emissions_header_init
2018       ENDIF
2019
2020    ENDIF
2021
2022
2023!
2024!-- Chemistry variables will be initialized if availabe from dynamic input file. Note, it is
2025!-- possible to initialize only part of the chemistry variables from dynamic input.
2026    IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
2027       DO  n = 1, nspec
2028          IF ( init_3d%from_file_chem(n) )  THEN
2029             DO  i = nxlg, nxrg
2030                DO  j = nysg, nyng
2031                   chem_species(n)%conc(:,j,i) = init_3d%chem_init(:,n)
2032                ENDDO
2033             ENDDO
2034          ENDIF
2035       ENDDO
2036    ENDIF
2037
2038    IF ( debug_output )  CALL debug_message( 'chem_init', 'end' )
2039
2040 END SUBROUTINE chem_init
2041
2042
2043!--------------------------------------------------------------------------------------------------!
2044! Description:
2045! ------------
2046!> Subroutine initializating chemistry_model_mod
2047!> internal workaround for chem_species dependency in chem_check_parameters
2048!--------------------------------------------------------------------------------------------------!
2049 SUBROUTINE chem_init_internal
2050
2051    USE pegrid
2052
2053    USE netcdf_data_input_mod,                                                                     &
2054        ONLY:  chem_emis, chem_emis_att, input_pids_dynamic, init_3d,                              &
2055               netcdf_data_input_chemistry_data
2056
2057!
2058!-- Local variables
2059    INTEGER(iwp) ::  i                 !< running index for for horiz numerical grid points
2060    INTEGER(iwp) ::  j                 !< running index for for horiz numerical grid points
2061    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
2062    INTEGER(iwp) ::  lsp               !< running index for chem spcs
2063
2064    REAL(wp)     ::  flag              !< flag for masking topography/building grid points
2065!
2066!-- 20200203 ECC
2067!-- reads netcdf data only under legacy mode
2068
2069!    IF ( emissions_anthropogenic )  THEN
2070!       CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis )
2071!    ENDIF
2072
2073    IF ( emissions_anthropogenic )  THEN
2074       IF ( emiss_read_legacy_mode )  THEN
2075          CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis )
2076       ENDIF
2077    ENDIF
2078
2079!
2080!-- Allocate memory for chemical species
2081    ALLOCATE( chem_species(nspec) )
2082    ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
2083    ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
2084    ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
2085    ALLOCATE( phot_frequen(nphot) )
2086    ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
2087    ALLOCATE( bc_cs_t_val(nspec) )
2088!
2089!-- Initialize arrays
2090    spec_conc_1 (:,:,:,:) = 0.0_wp
2091    spec_conc_2 (:,:,:,:) = 0.0_wp
2092    spec_conc_3 (:,:,:,:) = 0.0_wp
2093
2094!
2095!-- Allocate array to store locally summed-up resolved-scale vertical fluxes.
2096    IF ( scalar_advec == 'ws-scheme' )  THEN
2097       ALLOCATE( sums_ws_l(nzb:nzt+1,0:threads_per_task-1,nspec) )
2098       sums_ws_l = 0.0_wp
2099    ENDIF
2100
2101    DO  lsp = 1, nspec
2102       chem_species(lsp)%name    = spc_names(lsp)
2103
2104       chem_species(lsp)%conc   (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_1 (:,:,:,lsp)
2105       chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_2 (:,:,:,lsp)
2106       chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
2107
2108       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))
2109       chem_species(lsp)%cssws_av    = 0.0_wp
2110!
2111!--    The following block can be useful when emission module is not applied. &
2112!--    If emission module is applied the following block will be overwritten.
2113       ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))
2114       ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))
2115       ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
2116       ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
2117       chem_species(lsp)%flux_s_cs = 0.0_wp
2118       chem_species(lsp)%flux_l_cs = 0.0_wp
2119       chem_species(lsp)%diss_s_cs = 0.0_wp
2120       chem_species(lsp)%diss_l_cs = 0.0_wp
2121!
2122!--   Allocate memory for initial concentration profiles (concentration values come from namelist)
2123!--   (@todo (FK): Because of this, chem_init is called in palm before check_parameters, since
2124!--                conc_pr_init is used there.
2125!--                We have to find another solution since chem_init should eventually be called from
2126!--                init_3d_model!!)
2127       ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
2128       chem_species(lsp)%conc_pr_init(:) = 0.0_wp
2129
2130    ENDDO
2131
2132!
2133!-- For some passive scalars decycling may be enabled. This case, the lateral boundary conditions
2134!-- are non-cyclic for these scalars (chemical species and aerosols), while the other scalars may
2135!-- have cyclic boundary conditions. However, large gradients near the boundaries may produce
2136!-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
2137!-- applied near these boundaries.
2138!-- To get rid-off this, set-up additional flags that control the order of the scalar advection
2139!-- scheme near the lateral boundaries for passive scalars with decycling.
2140    IF ( scalar_advec == 'ws-scheme' )  THEN
2141       ALLOCATE( cs_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2142!
2143!--    In case of decyling, set Neumann boundary conditions for wall_flags_total_0 bit 31 instead of
2144!--    cyclic boundary conditions.
2145!--    Bit 31 is used to identify extended degradation zones (please see following comment).
2146!--    Note, since several also other modules like Salsa or other future one may access this bit but
2147!--    may have other boundary conditions, the original value of wall_flags_total_0 bit 31 must not
2148!--    be modified. Hence, store the boundary conditions directly on cs_advc_flags_s.
2149!--    cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used to
2150!--    control the numerical order.
2151!--    Initialize with flag 31 only.
2152       cs_advc_flags_s = 0
2153       cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0, BTEST( wall_flags_total_0, 31 ) )
2154
2155       IF ( bc_dirichlet_cs_n .OR. bc_dirichlet_cs_s )  THEN
2156          IF ( nys == 0  )  THEN
2157             DO  i = 1, nbgp
2158                cs_advc_flags_s(:,nys-i,:) = MERGE(                                                &
2159                                                    IBSET( cs_advc_flags_s(:,nys,:), 31 ),         &
2160                                                    IBCLR( cs_advc_flags_s(:,nys,:), 31 ),         &
2161                                                    BTEST( cs_advc_flags_s(:,nys,:), 31 )          &
2162                                                  )
2163             ENDDO
2164          ENDIF
2165          IF ( nyn == ny )  THEN
2166             DO  i = 1, nbgp
2167                cs_advc_flags_s(:,nyn+i,:) = MERGE(                                                &
2168                                                    IBSET( cs_advc_flags_s(:,nyn,:), 31 ),         &
2169                                                    IBCLR( cs_advc_flags_s(:,nyn,:), 31 ),         &
2170                                                    BTEST( cs_advc_flags_s(:,nyn,:), 31 )          &
2171                                                  )
2172             ENDDO
2173          ENDIF
2174       ENDIF
2175       IF ( bc_dirichlet_cs_l .OR. bc_dirichlet_cs_r )  THEN
2176          IF ( nxl == 0  )  THEN
2177             DO  i = 1, nbgp
2178                cs_advc_flags_s(:,:,nxl-i) = MERGE(                                                &
2179                                                    IBSET( cs_advc_flags_s(:,:,nxl), 31 ),         &
2180                                                    IBCLR( cs_advc_flags_s(:,:,nxl), 31 ),         &
2181                                                    BTEST( cs_advc_flags_s(:,:,nxl), 31 )          &
2182                                                  )
2183             ENDDO
2184          ENDIF
2185          IF ( nxr == nx )  THEN
2186             DO  i = 1, nbgp
2187                cs_advc_flags_s(:,:,nxr+i) = MERGE(                                                &
2188                                                    IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
2189                                                    IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
2190                                                    BTEST( cs_advc_flags_s(:,:,nxr), 31 )  &
2191                                                  )
2192             ENDDO
2193          ENDIF
2194
2195       ENDIF
2196!
2197!--    To initialize advection flags appropriately, pass the boundary flags.
2198!--    The last argument indicates that a passive scalar is treated, where the horizontal advection
2199!--    terms are degraded already 2 grid points before the lateral boundary to avoid stationary
2200!--    oscillations at large-gradients.
2201!--    Also, extended degradation zones are applied, where horizontal advection of passive scalars
2202!--    is discretized by first-order scheme at all grid points that in the vicinity of buildings
2203!--    (<= 3 grid points). Even though no building is within the numerical stencil, first-order
2204!--    scheme is used.
2205!--    At fourth and fifth grid point the order of the horizontal advection scheme
2206!--    is successively upgraded.
2207!--    These extended degradation zones are used to avoid stationary numerical oscillations, which
2208!--    are responsible for high concentration maxima that may appear under shear-free stable
2209!--    conditions.
2210       CALL ws_init_flags_scalar( bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                      &
2211                                  bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                      &
2212                                  bc_dirichlet_cs_r  .OR.  bc_radiation_cs_r,                      &
2213                                  bc_dirichlet_cs_s  .OR.  bc_radiation_cs_s,                      &
2214                                  cs_advc_flags_s, .TRUE. )
2215    ENDIF
2216!
2217!-- Initial concentration of profiles is prescribed by parameters cs_profile and cs_heights in the
2218!-- namelist &chemistry_parameters.
2219    CALL chem_init_profiles
2220!
2221!-- In case there is dynamic input file, create a list of names for chemistry initial input files.
2222!-- Also, initialize array that indicates whether the respective variable is on file or not.
2223    IF ( input_pids_dynamic )  THEN
2224       ALLOCATE( init_3d%var_names_chem(1:nspec) )
2225       ALLOCATE( init_3d%from_file_chem(1:nspec) )
2226       init_3d%from_file_chem(:) = .FALSE.
2227
2228       DO  lsp = 1, nspec
2229          init_3d%var_names_chem(lsp) = init_3d%init_char // TRIM( chem_species(lsp)%name )
2230       ENDDO
2231    ENDIF
2232!
2233!-- Initialize model variables
2234    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
2235         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
2236!
2237!--    First model run of a possible job queue.
2238!--    Initial profiles of the variables must be computed.
2239       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2240!
2241!--       Transfer initial profiles to the arrays of the 3D model
2242!--       Concentrations within buildings are set to zero.
2243          DO  lsp = 1, nspec
2244             DO  i = nxlg, nxrg
2245                DO  j = nysg, nyng
2246                   DO lpr_lev = 1, nz + 1
2247                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(lpr_lev,j,i), 0 ) )
2248                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)&
2249                                                            * flag
2250                   ENDDO
2251                ENDDO
2252             ENDDO
2253          ENDDO
2254
2255       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles') /= 0 )  THEN
2256
2257          DO  lsp = 1, nspec
2258             DO  i = nxlg, nxrg
2259                DO  j = nysg, nyng
2260                   DO  lpr_lev = nzb, nz+1
2261                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(lpr_lev,j,i), 0 ) )
2262                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)&
2263                                                            * flag
2264                   ENDDO
2265                ENDDO
2266             ENDDO
2267          ENDDO
2268
2269       ENDIF
2270!
2271!--    If required, change the surface chem spcs at the start of the 3D run
2272       IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN
2273          DO  lsp = 1, nspec
2274             DO  i = nxlg, nxrg
2275                DO  j = nysg, nyng
2276                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb,j,i), 0 ) )
2277                   chem_species(lsp)%conc(nzb,j,i) = chem_species(lsp)%conc(nzb,j,i) +             &
2278                                                     cs_surface_initial_change(lsp) * flag
2279                ENDDO
2280             ENDDO
2281          ENDDO
2282       ENDIF
2283
2284    ENDIF
2285!
2286!-- Initial old and new time levels. Note, this has to be done also in restart runs
2287    DO  lsp = 1, nvar
2288       chem_species(lsp)%tconc_m = 0.0_wp
2289       chem_species(lsp)%conc_p  = chem_species(lsp)%conc
2290    ENDDO
2291
2292    DO  lsp = 1, nphot
2293       phot_frequen(lsp)%name = phot_names(lsp)
2294!
2295!-- todo: remove or replace by "CALL message" mechanism (kanani)
2296!--       IF( myid == 0 )  THEN
2297!--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
2298!--       ENDIF
2299       phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
2300    ENDDO
2301
2302!    CALL photolysis_init   ! probably also required for restart
2303
2304    RETURN
2305
2306 END SUBROUTINE chem_init_internal
2307
2308
2309!--------------------------------------------------------------------------------------------------!
2310! Description:
2311! ------------
2312!> Subroutine defining initial vertical profiles of chemical species (given by namelist parameters
2313!> chem_profiles and chem_heights)  --> which should work analogically to parameters u_profile,
2314!> v_profile and uv_heights)
2315!--------------------------------------------------------------------------------------------------!
2316 SUBROUTINE chem_init_profiles
2317
2318    USE chem_modules
2319
2320!
2321!-- Local variables
2322    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
2323    INTEGER ::  lsp        !< running index for number of species in derived data type species_def
2324    INTEGER ::  lsp_usr    !< running index for number of species (user defined)  in cs_names,
2325                           !< cs_profiles etc
2326    INTEGER ::  npr_lev    !< the next available profile lev
2327!
2328!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
2329!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
2330!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
2331!-- "cs_heights" are prescribed, their values will!override the constant profile given by
2332!-- "cs_surface".
2333!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2334    lsp_usr = 1
2335    DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
2336       DO  lsp = 1, nspec                                !
2337!
2338!--       Create initial profile (conc_pr_init) for each chemical species
2339          IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
2340             IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
2341!
2342!--            Set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of
2343!--            each species
2344                DO lpr_lev = 0, nzt+1
2345                   chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
2346                ENDDO
2347             ELSE
2348                IF ( cs_heights(1,1) /= 0.0_wp )  THEN
2349                   message_string = 'The surface value of cs_heights must be 0.0'
2350                   CALL message( 'chem_check_parameters', 'PA0734', 1, 2, 0, 6, 0 )
2351                ENDIF
2352
2353                use_prescribed_profile_data = .TRUE.
2354
2355                npr_lev = 1
2356!                chem_species(lsp)%conc_pr_init(0) = 0.0_wp
2357                DO  lpr_lev = 1, nz+1
2358                   IF ( npr_lev < 100 )  THEN
2359                      DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
2360                         npr_lev = npr_lev + 1
2361                         IF ( npr_lev == 100 )  THEN
2362                            message_string = 'number of chem spcs exceeding the limit'
2363                            CALL message( 'chem_check_parameters', 'PA0735', 1, 2, 0, 6, 0 )
2364                            EXIT
2365                         ENDIF
2366                      ENDDO
2367                   ENDIF
2368                   IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
2369                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +     &
2370                           ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                        &
2371                           ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *&
2372                           ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
2373                   ELSE
2374                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
2375                   ENDIF
2376                ENDDO
2377             ENDIF
2378!
2379!--       If a profile is prescribed explicity using cs_profiles and cs_heights, then
2380!--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based on the
2381!--       cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
2382          ENDIF
2383
2384       ENDDO
2385
2386       lsp_usr = lsp_usr + 1
2387    ENDDO
2388!     ENDIF
2389
2390 END SUBROUTINE chem_init_profiles
2391
2392
2393!--------------------------------------------------------------------------------------------------!
2394! Description:
2395! ------------
2396!> Subroutine to integrate chemical species in the given chemical mechanism
2397!--------------------------------------------------------------------------------------------------!
2398 SUBROUTINE chem_integrate_ij( i, j )
2399
2400    USE statistics,                                                                                &
2401        ONLY:  weight_pres
2402
2403    USE control_parameters,                                                                        &
2404        ONLY:  dt_3d, intermediate_timestep_count, time_since_reference_point
2405
2406    REAL(wp), PARAMETER              ::  fr2ppm  = 1.0e6_wp              !< Conversion factor fraction to ppm
2407!    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp              !< Mole mass of dry air
2408!    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp           !< Mole mass of water vapor
2409    REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp           !< standard pressure (Pa)
2410    REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp             !< Conversion factor ppm to fraction
2411    REAL(wp), PARAMETER              ::  t_std   = 273.15_wp             !< standard pressure (Pa)
2412    REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp           !< Mole volume (22.414 l) in cm^3
2413    REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp           !< Avogadro number (molecules/mol)
2414
2415    INTEGER,INTENT(IN)       :: i
2416    INTEGER,INTENT(IN)       :: j
2417!
2418!-- Local variables
2419    INTEGER(iwp) ::  lph                                 !< running index for photolysis frequencies
2420    INTEGER(iwp) ::  lsp                                 !< running index for chem spcs.
2421
2422    INTEGER, DIMENSION(20)        :: istatus
2423
2424    INTEGER,DIMENSION(nzb+1:nzt)  :: nacc          !< Number of accepted steps
2425    INTEGER,DIMENSION(nzb+1:nzt)  :: nrej          !< Number of rejected steps
2426
2427    REAL(wp)                         ::  conv                                !< conversion factor
2428    REAL(kind=wp)                    ::  dt_chem
2429
2430    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
2431
2432    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
2433    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
2434                                                                              !< molecules cm^{-3} and ppm
2435    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
2436    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
2437
2438    REAL(kind=wp), DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
2439    REAL(kind=wp), DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
2440
2441!
2442!-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
2443    IF (chem_gasphase_on)  THEN
2444       nacc = 0
2445       nrej = 0
2446
2447       tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
2448!
2449!--    Convert ppm to molecules/cm**3
2450!--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp *
2451!--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )
2452       conv = ppm2fr * xna / vmolcm
2453       tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
2454       tmp_fact_i = 1.0_wp/tmp_fact
2455
2456       IF ( humidity )  THEN
2457          IF ( bulk_cloud_model )  THEN
2458             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *                              &
2459                             xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2460          ELSE
2461             tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2462          ENDIF
2463       ELSE
2464          tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:)   !< Constant value for q if
2465                                                                      !< water vapor is not computed
2466       ENDIF
2467
2468       DO  lsp = 1,nspec
2469          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
2470       ENDDO
2471
2472       DO lph = 1,nphot
2473          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)
2474       ENDDO
2475!
2476!--    Compute length of time step
2477       IF ( call_chem_at_all_substeps )  THEN
2478          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
2479       ELSE
2480          dt_chem = dt_3d
2481       ENDIF
2482
2483       cs_time_step = dt_chem
2484
2485       IF ( MAXVAL( rcntrl ) > 0.0 )  THEN    ! Only if rcntrl is set
2486          IF( time_since_reference_point <= 2*dt_3d)  THEN
2487             rcntrl_local = 0
2488          ELSE
2489             rcntrl_local = rcntrl
2490          ENDIF
2491       ELSE
2492          rcntrl_local = 0
2493       END IF
2494
2495       CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot,   &
2496            icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
2497
2498       DO  lsp = 1,nspec
2499          chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
2500       ENDDO
2501
2502
2503    ENDIF
2504
2505    RETURN
2506 END SUBROUTINE chem_integrate_ij
2507
2508
2509!--------------------------------------------------------------------------------------------------!
2510! Description:
2511! ------------
2512!> Subroutine defining parin for &chemistry_parameters for chemistry model
2513!--------------------------------------------------------------------------------------------------!
2514 SUBROUTINE chem_parin
2515
2516    USE chem_modules
2517    USE control_parameters
2518    USE pegrid
2519    USE statistics
2520
2521    CHARACTER(LEN=100) ::  line                        !< dummy string that contains the current
2522                                                       !< line of the parameter file
2523    CHARACTER(LEN=8)   ::  solver_type
2524
2525    INTEGER(iwp) ::  i                                 !<
2526    INTEGER(iwp) ::  io_status                         !< Status after reading the namelist file
2527    INTEGER(iwp) ::  max_pr_cs_tmp                     !<
2528
2529    LOGICAL ::  switch_off_module = .FALSE.  !< local namelist parameter to switch off the module
2530                                             !< although the respective module namelist appears in
2531                                             !< the namelist file
2532
2533    REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0
2534                                                       !< automatic stepping
2535
2536
2537    NAMELIST /chemistry_parameters/                                                                &
2538         bc_cs_b,                                                                                  &
2539         bc_cs_l,                                                                                  &
2540         bc_cs_n,                                                                                  &
2541         bc_cs_r,                                                                                  &
2542         bc_cs_s,                                                                                  &
2543         bc_cs_t,                                                                                  &
2544         call_chem_at_all_substeps,                                                                &
2545         chem_debug0,                                                                              &
2546         chem_debug1,                                                                              &
2547         chem_debug2,                                                                              &
2548         chem_gasphase_on,                                                                         &
2549         chem_mechanism,                                                                           &
2550         cs_heights,                                                                               &
2551         cs_name,                                                                                  &
2552         cs_profile,                                                                               &
2553         cs_surface,                                                                               &
2554         cs_surface_initial_change,                                                                &
2555         cs_vertical_gradient_level,                                                               &
2556         daytype_mdh,                                                                              &
2557         deposition_dry,                                                                           &
2558         emissions_anthropogenic,                                                                  &
2559         emiss_lod,                                                                                &
2560         emiss_factor_main,                                                                        &
2561         emiss_factor_side,                                                                        &
2562         emiss_read_legacy_mode,                                                                   &
2563         icntrl,                                                                                   &
2564         main_street_id,                                                                           &
2565         max_street_id,                                                                            &
2566         mode_emis,                                                                                &
2567         my_steps,                                                                                 &
2568         nesting_chem,                                                                             &
2569         nesting_offline_chem,                                                                     &
2570         rcntrl,                                                                                   &
2571         side_street_id,                                                                           &
2572         photolysis_scheme,                                                                        &
2573         wall_csflux,                                                                              &
2574         cs_vertical_gradient,                                                                     &
2575         top_csflux,                                                                               &
2576         surface_csflux,                                                                           &
2577         surface_csflux_name,                                                                      &
2578         switch_off_module,                                                                        &
2579         time_fac_type
2580!
2581!-- Analogically to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
2582!-- so this way we could prescribe a specific flux value for each species
2583    !>  chemistry_parameters for initial profiles
2584    !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
2585    !>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
2586    !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
2587    !>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
2588    !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
2589    !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
2590
2591!
2592!-- Read chem namelist
2593    icntrl    = 0
2594    rcntrl    = 0.0_wp
2595    my_steps  = 0.0_wp
2596    photolysis_scheme = 'simple'
2597    atol = 1.0_wp
2598    rtol = 0.01_wp
2599!
2600!-- Move to the beginning of the namelist file and try to find and read the namelist named
2601!-- chemistry_parameters.
2602    REWIND( 11 )
2603    READ( 11, chemistry_parameters, IOSTAT=io_status )
2604!
2605!-- Action depending on the READ status
2606    IF ( io_status == 0 )  THEN
2607!
2608!      chemistry_parameters namelist was found and read correctly. Switch on chemistry model.
2609       IF ( .NOT. switch_off_module )  air_chemistry = .TRUE.
2610
2611    ELSEIF ( io_status > 0 )  THEN
2612!
2613!--    chemistry_parameters namelist was found, but contained errors. Print an error message
2614!--    including the line that caused the problem.
2615       BACKSPACE( 11 )
2616       READ( 11 , '(A)') line
2617       CALL parin_fail_message( 'chemistry_parameters', line )
2618
2619    ENDIF
2620
2621!
2622!-- Synchronize emiss_lod and mod_emis only if emissions_anthropogenic is activated in the namelist.
2623!-- Otherwise their values are "don't care"
2624    IF ( emissions_anthropogenic )  THEN
2625
2626!
2627!--    Check for emission mode for chem species
2628       IF ( emiss_lod < 0 )  THEN   !- if LOD not defined in namelist
2629          IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.                                          &
2630               ( mode_emis /= 'DEFAULT'        )    .AND.                                          &
2631               ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
2632             message_string = 'Incorrect mode_emiss  option select. Please check spelling'
2633             CALL message( 'chem_parin', 'PA0736', 1, 2, 0, 6, 0 )
2634          ENDIF
2635       ELSE
2636          IF ( ( emiss_lod /= 0 )    .AND.                                                         &
2637               ( emiss_lod /= 1 )    .AND.                                                         &
2638               ( emiss_lod /= 2 ) )  THEN
2639             message_string = 'Invalid value for emiss_lod (0, 1, or 2)'
2640             CALL message( 'chem_parin', 'PA0737', 1, 2, 0, 6, 0 )
2641          ENDIF
2642       ENDIF
2643
2644!
2645! For reference (ecc)
2646!    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
2647!       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
2648!       CALL message( 'chem_parin', 'PA0736', 1, 2, 0, 6, 0 )
2649!    ENDIF
2650
2651!
2652!-- Conflict resolution for emiss_lod and mode_emis
2653!-- 1) if emiss_lod is defined, have mode_emis assume same setting as emiss_lod
2654!-- 2) if emiss_lod it not defined, have emiss_lod assuem same setting as mode_emis
2655!-- this check is in place to retain backward compatibility with salsa until the code is migrated
2656!-- completed to emiss_lod
2657!-- note that
2658       IF  ( emiss_lod >= 0 ) THEN
2659
2660          SELECT CASE  ( emiss_lod )
2661             CASE (0)  !- parameterized mode
2662                mode_emis = 'PARAMETERIZED'
2663             CASE (1)  !- default mode
2664                mode_emis = 'DEFAULT'
2665             CASE (2)  !- preprocessed mode
2666                mode_emis = 'PRE-PROCESSED'
2667          END SELECT
2668
2669          message_string = 'Synchronizing mode_emis to defined emiss_lod&'              //         &
2670                           'NOTE - mode_emis will be depreciated in future releases&'   //         &
2671                           'please use emiss_lod to define emission mode'
2672          CALL message( 'chem_parin', 'PA0738', 0, 0, 0, 6, 0 )
2673
2674       ELSE ! if emiss_lod is not set
2675
2676          SELECT CASE ( mode_emis )
2677             CASE ('PARAMETERIZED')
2678                emiss_lod = 0
2679             CASE ('DEFAULT')
2680                emiss_lod = 1
2681             CASE ('PRE-PROCESSED')
2682                emiss_lod = 2
2683          END SELECT
2684
2685          message_string = 'emiss_lod undefined.  Using existing mod_emiss setting&'    //         &
2686                           'NOTE - mode_emis will be depreciated in future releases&'   //         &
2687                           'please use emiss_lod to define emission mode'
2688          CALL message( 'chem_parin', 'PA0739', 0, 0, 0, 6, 0 )
2689       ENDIF
2690
2691!
2692!-- (ECC) input check for emission read mode.
2693!--    legacy : business as usual (everything read / set up at start of run)
2694!--    new    : emission based on timestamp, and for lod2 data is loaded on an hourly basis
2695
2696!
2697!-- (ECC) handler for emiss_read_legacy_mode
2698!-- * emiss_read_legacy_mode is defaulted to TRUE
2699!-- * if emiss_read_legacy_mode is TRUE and LOD is 0 or 1,
2700!--       force emission_read_legacy_mode to TRUE (not yet implemented)
2701
2702       IF ( emiss_read_legacy_mode )  THEN       !< notify legacy read mode
2703
2704          message_string = 'Legacy emission read mode activated&'           //                     &
2705                           'All emissions data will be loaded '             //                     &
2706                           'prior to start of simulation'
2707          CALL message( 'chem_parin', 'PA0740', 0, 0, 0, 6, 0 )
2708
2709       ELSE                                     !< if new read mode selected
2710
2711          IF ( emiss_lod < 2 )  THEN            !< check LOD compatibility
2712
2713             message_string = 'New emission read mode '                    //                      &
2714                              'currently unavailable for LODs 0 and 1.&'   //                      &
2715                              'Reverting to legacy emission read mode'
2716             CALL message( 'chem_parin', 'PA0741', 0, 0, 0, 6, 0 )
2717
2718             emiss_read_legacy_mode = .TRUE.
2719
2720          ELSE                                  !< notify new read mode
2721
2722             message_string = 'New emission read mode activated&'          //                      &
2723                              'LOD 2 emissions will be updated on-demand ' //                      &
2724                              'according to indicated timestamps'
2725             CALL message( 'chem_parin', 'PA0742', 0, 0, 0, 6, 0 )
2726
2727          ENDIF
2728
2729       ENDIF ! if emiss_read_legacy_mode
2730
2731
2732    ENDIF  ! if emissions_anthropengic
2733
2734
2735    t_steps = my_steps
2736!
2737!-- Determine the number of user-defined profiles and append them to the standard data output
2738!-- (data_output_pr)
2739    max_pr_cs_tmp = 0
2740    i = 1
2741
2742    DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= SIZE( data_output_pr ) )
2743       IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' )  THEN
2744          max_pr_cs_tmp = max_pr_cs_tmp+1
2745       ENDIF
2746       i = i +1
2747    ENDDO
2748
2749    IF ( max_pr_cs_tmp > 0 )  THEN
2750       cs_pr_namelist_found = .TRUE.
2751       max_pr_cs = max_pr_cs_tmp
2752    ENDIF
2753!
2754!-- Set Solver Type
2755    IF ( icntrl(3) == 0 )  THEN
2756       solver_type = 'rodas3'           !Default
2757    ELSEIF ( icntrl(3) == 1 )  THEN
2758       solver_type = 'ros2'
2759    ELSEIF ( icntrl(3) == 2 )  THEN
2760       solver_type = 'ros3'
2761    ELSEIF ( icntrl(3) == 3 )  THEN
2762       solver_type = 'ro4'
2763    ELSEIF ( icntrl(3) == 4 )  THEN
2764       solver_type = 'rodas3'
2765    ELSEIF ( icntrl(3) == 5 )  THEN
2766       solver_type = 'rodas4'
2767    ELSEIF ( icntrl(3) == 6 )  THEN
2768       solver_type = 'Rang3'
2769    ELSE
2770       message_string = 'illegal Rosenbrock-solver type'
2771       CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
2772    END IF
2773
2774!
2775!--   todo: remove or replace by "CALL message" mechanism (kanani)
2776!       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM( solver_type )
2777!kk    Has to be changed to right calling sequence
2778!        IF(myid == 0)  THEN
2779!           write(9,*) ' '
2780!           write(9,*) 'kpp setup '
2781!           write(9,*) ' '
2782!           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM( solver_type )
2783!           write(9,*) ' '
2784!           write(9,*) '    Hstart  = ',rcntrl(3)
2785!           write(9,*) '    FacMin  = ',rcntrl(4)
2786!           write(9,*) '    FacMax  = ',rcntrl(5)
2787!           write(9,*) ' '
2788!           IF(vl_dim > 1)  THEN
2789!              write(9,*) '    Vector mode                   vektor length = ',vl_dim
2790!           ELSE
2791!              write(9,*) '    Scalar mode'
2792!           ENDIF
2793!           write(9,*) ' '
2794!        END IF
2795
2796    RETURN
2797
2798 END SUBROUTINE chem_parin
2799
2800
2801!--------------------------------------------------------------------------------------------------!
2802! Description:
2803! ------------
2804!> Call for all grid points
2805!--------------------------------------------------------------------------------------------------!
2806    SUBROUTINE chem_actions( location )
2807
2808
2809    CHARACTER (LEN=*), INTENT(IN) ::  location !< call location string
2810
2811    SELECT CASE ( location )
2812
2813       CASE ( 'before_prognostic_equations' )
2814!
2815!--       Chemical reactions and deposition
2816          IF ( chem_gasphase_on )  THEN
2817!
2818!--          If required, calculate photolysis frequencies -
2819!--          UNFINISHED: Why not before the intermediate timestep loop?
2820             IF ( intermediate_timestep_count ==  1 )  THEN
2821                CALL photolysis_control
2822             ENDIF
2823
2824          ENDIF
2825
2826       CASE ( 'before_timestep' )
2827!
2828!--       Set array used to sum-up resolved scale fluxes to zero.
2829          IF ( ws_scheme_sca )  THEN
2830             sums_ws_l = 0.0_wp
2831          ENDIF
2832
2833       CASE DEFAULT
2834          CONTINUE
2835
2836    END SELECT
2837
2838    END SUBROUTINE chem_actions
2839
2840
2841!--------------------------------------------------------------------------------------------------!
2842! Description:
2843! ------------
2844!> Call for grid points i,j
2845!--------------------------------------------------------------------------------------------------!
2846
2847    SUBROUTINE chem_actions_ij( i, j, location )
2848
2849    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
2850
2851    INTEGER(iwp)  ::  dummy                     !< call location string
2852
2853    INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
2854    INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
2855
2856    IF ( air_chemistry    )   dummy = i + j
2857
2858    SELECT CASE ( location )
2859
2860       CASE DEFAULT
2861          CONTINUE
2862
2863    END SELECT
2864
2865
2866    END SUBROUTINE chem_actions_ij
2867
2868
2869!--------------------------------------------------------------------------------------------------!
2870! Description:
2871! ------------
2872!> Call for all grid points
2873!--------------------------------------------------------------------------------------------------!
2874    SUBROUTINE chem_non_advective_processes()
2875
2876
2877      INTEGER(iwp) ::  i  !<
2878      INTEGER(iwp) ::  j  !<
2879
2880!
2881!--   Calculation of chemical reactions and deposition.
2882      IF ( intermediate_timestep_count == 1  .OR.  call_chem_at_all_substeps )  THEN
2883
2884         IF ( chem_gasphase_on )  THEN
2885            CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2886            !$OMP PARALLEL PRIVATE (i,j)
2887            !$OMP DO schedule(static,1)
2888            DO  i = nxl, nxr
2889               DO  j = nys, nyn
2890                  CALL chem_integrate( i, j )
2891               ENDDO
2892            ENDDO
2893            !$OMP END PARALLEL
2894            CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2895         ENDIF
2896
2897         IF ( deposition_dry )  THEN
2898            CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2899            DO  i = nxl, nxr
2900               DO  j = nys, nyn
2901                  CALL chem_depo( i, j )
2902               ENDDO
2903            ENDDO
2904            CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2905         ENDIF
2906
2907      ENDIF
2908
2909    END SUBROUTINE chem_non_advective_processes
2910
2911
2912!--------------------------------------------------------------------------------------------------!
2913! Description:
2914! ------------
2915!> Call for grid points i,j
2916!--------------------------------------------------------------------------------------------------!
2917 SUBROUTINE chem_non_advective_processes_ij( i, j )
2918
2919
2920    INTEGER(iwp), INTENT(IN) ::  i  !< grid index in x-direction
2921    INTEGER(iwp), INTENT(IN) ::  j  !< grid index in y-direction
2922
2923!
2924!-- Calculation of chemical reactions and deposition.
2925!-- !> TODO
2926!-- WARNING: time measurements within i,j loops degrade performance because cpu-log routine is
2927!--          called extremely often. Furthermore, because of that, the counter for this measurement
2928!--          will get an extremely huge value, too.
2929!--          Should be removed, because in other _ij routines no measurements are done because of
2930!--          the reasons given above.
2931    IF ( intermediate_timestep_count == 1  .OR.  call_chem_at_all_substeps )  THEN
2932
2933       IF ( chem_gasphase_on )  THEN
2934          !$OMP MASTER
2935          CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2936          !$OMP END MASTER
2937          CALL chem_integrate( i, j )
2938          !$OMP MASTER
2939          CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2940          !$OMP END MASTER
2941       ENDIF
2942
2943       IF ( deposition_dry )  THEN
2944          !$OMP MASTER
2945          CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2946          !$OMP END MASTER
2947          CALL chem_depo( i, j )
2948          !$OMP MASTER
2949          CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2950          !$OMP END MASTER
2951       ENDIF
2952
2953    ENDIF
2954
2955 END SUBROUTINE chem_non_advective_processes_ij
2956
2957
2958!--------------------------------------------------------------------------------------------------!
2959! Description:
2960! ------------
2961!> routine for exchange horiz of chemical quantities
2962!--------------------------------------------------------------------------------------------------!
2963 SUBROUTINE chem_exchange_horiz_bounds( location )
2964
2965    USE exchange_horiz_mod,                                                                        &
2966        ONLY:  exchange_horiz
2967
2968   INTEGER(iwp) ::  lsp       !<
2969   INTEGER(iwp) ::  n
2970
2971   CHARACTER (LEN=*), INTENT(IN) ::  location !< call location string
2972
2973   SELECT CASE ( location )
2974
2975       CASE ( 'before_prognostic_equation' )
2976!
2977!--       Loop over chemical species
2978          CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
2979          DO  lsp = 1, nvar
2980             CALL exchange_horiz( chem_species(lsp)%conc, nbgp,                                    &
2981                                  alternative_communicator = communicator_chem )
2982          ENDDO
2983
2984          CALL chem_boundary_conditions( horizontal_conditions_only = .TRUE. )
2985
2986          CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
2987
2988       CASE ( 'after_prognostic_equation' )
2989
2990          IF ( air_chemistry )  THEN
2991             DO  n = 1, nvar
2992                CALL exchange_horiz( chem_species(n)%conc_p, nbgp,                                 &
2993                                     alternative_communicator = communicator_chem )
2994             ENDDO
2995          ENDIF
2996
2997       CASE ( 'after_anterpolation' )
2998
2999          IF ( air_chemistry )  THEN
3000             DO  n = 1, nvar
3001                CALL exchange_horiz( chem_species(n)%conc, nbgp,                                   &
3002                                     alternative_communicator = communicator_chem )
3003             ENDDO
3004          ENDIF
3005
3006    END SELECT
3007
3008 END SUBROUTINE chem_exchange_horiz_bounds
3009
3010
3011!--------------------------------------------------------------------------------------------------!
3012! Description:
3013! ------------
3014!> Subroutine calculating prognostic equations for chemical species (vector-optimized).
3015!> Routine is called separately for each chemical species over a loop from prognostic_equations.
3016!--------------------------------------------------------------------------------------------------!
3017 SUBROUTINE chem_prognostic_equations()
3018
3019
3020    INTEGER ::  i   !< running index
3021    INTEGER ::  j   !< running index
3022    INTEGER ::  k   !< running index
3023
3024    INTEGER(iwp) ::  ilsp   !<
3025
3026
3027    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
3028
3029    DO  ilsp = 1, nvar
3030!
3031!--    Tendency terms for chemical species
3032       tend = 0.0_wp
3033!
3034!--    Advection terms
3035       IF ( timestep_scheme(1:5) == 'runge' )  THEN
3036          IF ( ws_scheme_sca )  THEN
3037             sums_wschs_ws_l(nzb:,0:) => sums_ws_l(:,:,ilsp)
3038             CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc',                      &
3039                              bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                          &
3040                              bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                          &
3041                              bc_dirichlet_cs_r  .OR.  bc_radiation_cs_r,                          &
3042                              bc_dirichlet_cs_s  .OR.  bc_radiation_cs_s )
3043          ELSE
3044             CALL advec_s_pw( chem_species(ilsp)%conc )
3045          ENDIF
3046       ELSE
3047          CALL advec_s_up( chem_species(ilsp)%conc )
3048       ENDIF
3049!
3050!--    Diffusion terms  (the last three arguments are zero)
3051       CALL diffusion_s( chem_species(ilsp)%conc,                                                  &
3052            surf_def_h(0)%cssws(ilsp,:),                                                           &
3053            surf_def_h(1)%cssws(ilsp,:),                                                           &
3054            surf_def_h(2)%cssws(ilsp,:),                                                           &
3055            surf_lsm_h(0)%cssws(ilsp,:),                                                           &
3056            surf_lsm_h(1)%cssws(ilsp,:),                                                           &
3057            surf_usm_h(0)%cssws(ilsp,:),                                                           &
3058            surf_usm_h(1)%cssws(ilsp,:),                                                           &
3059            surf_def_v(0)%cssws(ilsp,:),                                                           &
3060            surf_def_v(1)%cssws(ilsp,:),                                                           &
3061            surf_def_v(2)%cssws(ilsp,:),                                                           &
3062            surf_def_v(3)%cssws(ilsp,:),                                                           &
3063            surf_lsm_v(0)%cssws(ilsp,:),                                                           &
3064            surf_lsm_v(1)%cssws(ilsp,:),                                                           &
3065            surf_lsm_v(2)%cssws(ilsp,:),                                                           &
3066            surf_lsm_v(3)%cssws(ilsp,:),                                                           &
3067            surf_usm_v(0)%cssws(ilsp,:),                                                           &
3068            surf_usm_v(1)%cssws(ilsp,:),                                                           &
3069            surf_usm_v(2)%cssws(ilsp,:),                                                           &
3070            surf_usm_v(3)%cssws(ilsp,:) )
3071!
3072!--    Prognostic equation for chemical species
3073       DO  i = nxl, nxr
3074          DO  j = nys, nyn
3075             !following directive is required to vectorize on Intel19
3076             !DIR$ IVDEP
3077             DO  k = nzb+1, nzt
3078                chem_species(ilsp)%conc_p(k,j,i) =   chem_species(ilsp)%conc(k,j,i)                &
3079                     + ( dt_3d  *                                                                  &
3080                     (   tsc(2) * tend(k,j,i)                                                      &
3081                     + tsc(3) * chem_species(ilsp)%tconc_m(k,j,i)                                  &
3082                     )                                                                             &
3083                     - tsc(5) * rdf_sc(k)                                                          &
3084                     * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )     &
3085                     )                                                                             &
3086                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
3087
3088                IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
3089                   chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)
3090                ENDIF
3091             ENDDO
3092          ENDDO
3093       ENDDO
3094!
3095!--    Calculate tendencies for the next Runge-Kutta step
3096       IF ( timestep_scheme(1:5) == 'runge' )  THEN
3097          IF ( intermediate_timestep_count == 1 )  THEN
3098             DO  i = nxl, nxr
3099                DO  j = nys, nyn
3100                   DO  k = nzb+1, nzt
3101                      chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
3102                   ENDDO
3103                ENDDO
3104             ENDDO
3105          ELSEIF ( intermediate_timestep_count < &
3106               intermediate_timestep_count_max )  THEN
3107             DO  i = nxl, nxr
3108                DO  j = nys, nyn
3109                   DO  k = nzb+1, nzt
3110                      chem_species(ilsp)%tconc_m(k,j,i) = - 9.5625_wp * tend(k,j,i)                &
3111                                                     + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
3112                   ENDDO
3113                ENDDO
3114             ENDDO
3115          ENDIF
3116       ENDIF
3117
3118    ENDDO
3119
3120    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )
3121
3122 END SUBROUTINE chem_prognostic_equations
3123
3124
3125!--------------------------------------------------------------------------------------------------!
3126! Description:
3127! ------------
3128!> Subroutine calculating prognostic equations for chemical species (cache-optimized).
3129!> Routine is called separately for each chemical species over a loop from prognostic_equations.
3130!--------------------------------------------------------------------------------------------------!
3131 SUBROUTINE chem_prognostic_equations_ij( i, j, i_omp_start, tn )
3132
3133
3134    INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn
3135
3136    INTEGER(iwp) :: ilsp
3137!
3138!-- local variables
3139
3140    INTEGER :: k
3141
3142    DO  ilsp = 1, nvar
3143!
3144!--    Tendency-terms for chem spcs.
3145       tend(:,j,i) = 0.0_wp
3146!
3147!--    Advection terms
3148       IF ( timestep_scheme(1:5) == 'runge' )  THEN
3149          IF ( ws_scheme_sca )  THEN
3150             sums_wschs_ws_l(nzb:,0:) => sums_ws_l(nzb:nzt+1,0:threads_per_task-1,ilsp)
3151             CALL advec_s_ws( cs_advc_flags_s,                                                     &
3152                              i,                                                                   &
3153                              j,                                                                   &
3154                              chem_species(ilsp)%conc,                                             &
3155                              'kc',                                                                &
3156                              chem_species(ilsp)%flux_s_cs,                                        &
3157                              chem_species(ilsp)%diss_s_cs,                                        &
3158                              chem_species(ilsp)%flux_l_cs,                                        &
3159                              chem_species(ilsp)%diss_l_cs,                                        &
3160                              i_omp_start,                                                         &
3161                              tn,                                                                  &
3162                              bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                          &
3163                              bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                          &
3164                              bc_dirichlet_cs_r  .OR.  bc_radiation_cs_r,                          &
3165                              bc_dirichlet_cs_s  .OR.  bc_radiation_cs_s,                          &
3166                              monotonic_limiter_z )
3167          ELSE
3168             CALL advec_s_pw( i, j, chem_species(ilsp)%conc )
3169          ENDIF
3170       ELSE
3171          CALL advec_s_up( i, j, chem_species(ilsp)%conc )
3172       ENDIF
3173!
3174!--    Diffusion terms (the last three arguments are zero)
3175       CALL diffusion_s( i, j, chem_species(ilsp)%conc,                                            &
3176            surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),                              &
3177            surf_def_h(2)%cssws(ilsp,:),                                                           &
3178            surf_lsm_h(0)%cssws(ilsp,:), surf_lsm_h(1)%cssws(ilsp,:),                              &
3179            surf_usm_h(0)%cssws(ilsp,:), surf_usm_h(1)%cssws(ilsp,:),                              &
3180            surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),                              &
3181            surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),                              &
3182            surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),                              &
3183            surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),                              &
3184            surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),                              &
3185            surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
3186!
3187!--    Prognostic equation for chem spcs
3188       DO  k = nzb+1, nzt
3189          chem_species(ilsp)%conc_p(k,j,i) = chem_species(ilsp)%conc(k,j,i) + ( dt_3d  *           &
3190               ( tsc(2) * tend(k,j,i) +                                                            &
3191               tsc(3) * chem_species(ilsp)%tconc_m(k,j,i) )                                        &
3192               - tsc(5) * rdf_sc(k)                                                                &
3193               * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )           &
3194               )                                                                                   &
3195               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )                      &
3196               )
3197
3198          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
3199             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6
3200          ENDIF
3201       ENDDO
3202!
3203!--    Calculate tendencies for the next Runge-Kutta step
3204       IF ( timestep_scheme(1:5) == 'runge' )  THEN
3205          IF ( intermediate_timestep_count == 1 )  THEN
3206             DO  k = nzb+1, nzt
3207                chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
3208             ENDDO
3209          ELSEIF ( intermediate_timestep_count < &
3210               intermediate_timestep_count_max )  THEN
3211             DO  k = nzb+1, nzt
3212                chem_species(ilsp)%tconc_m(k,j,i) = -9.5625_wp * tend(k,j,i) +                     &
3213                                                     5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
3214             ENDDO
3215          ENDIF
3216       ENDIF
3217
3218    ENDDO
3219
3220 END SUBROUTINE chem_prognostic_equations_ij
3221
3222
3223!--------------------------------------------------------------------------------------------------!
3224! Description:
3225! ------------
3226!> Read module-specific local restart data arrays (Fortran binary format).
3227!--------------------------------------------------------------------------------------------------!
3228 SUBROUTINE chem_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,   &
3229                                nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
3230
3231    USE control_parameters
3232
3233    INTEGER(iwp) ::  k               !<
3234    INTEGER(iwp) ::  lsp             !<
3235    INTEGER(iwp) ::  nxlc            !<
3236    INTEGER(iwp) ::  nxlf            !<
3237    INTEGER(iwp) ::  nxl_on_file     !<
3238    INTEGER(iwp) ::  nxrc            !<
3239    INTEGER(iwp) ::  nxrf            !<
3240    INTEGER(iwp) ::  nxr_on_file     !<
3241    INTEGER(iwp) ::  nync            !<
3242    INTEGER(iwp) ::  nynf            !<
3243    INTEGER(iwp) ::  nyn_on_file     !<
3244    INTEGER(iwp) ::  nysc            !<
3245    INTEGER(iwp) ::  nysf            !<
3246    INTEGER(iwp) ::  nys_on_file     !<
3247
3248    LOGICAL, INTENT(OUT) :: found
3249
3250    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) &
3251                 :: tmp_3d   !< 3D array to temp store data
3252
3253
3254    found = .FALSE.
3255
3256
3257    IF ( ALLOCATED( chem_species ) )  THEN
3258
3259       DO  lsp = 1, nspec
3260
3261          IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) )  THEN
3262
3263             IF ( k == 1 )  READ ( 13 )  tmp_3d
3264             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                   &
3265                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3266             found = .TRUE.
3267
3268          ELSEIF (restart_string(1:length) == TRIM( chem_species(lsp)%name ) // '_av' )  THEN
3269
3270             IF ( .NOT. ALLOCATED( chem_species(lsp)%conc_av ) )  THEN
3271                ALLOCATE( chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg ) )
3272             ENDIF
3273             IF ( k == 1 )  READ ( 13 )  tmp_3d
3274             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                &
3275                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3276             found = .TRUE.
3277
3278          ENDIF
3279
3280       ENDDO
3281
3282    ENDIF
3283
3284 END SUBROUTINE chem_rrd_local_ftn
3285
3286
3287!--------------------------------------------------------------------------------------------------!
3288! Description:
3289! ------------
3290!> Read module-specific local restart data arrays (Fortran binary format).
3291!--------------------------------------------------------------------------------------------------!
3292 SUBROUTINE chem_rrd_local_mpi
3293
3294    IMPLICIT NONE
3295
3296    INTEGER(iwp) ::  lsp  !<
3297
3298    LOGICAL      ::  array_found  !<
3299
3300
3301    DO  lsp = 1, nspec
3302
3303       CALL rrd_mpi_io( TRIM( chem_species(lsp)%name ), chem_species(lsp)%conc )
3304
3305       CALL rd_mpi_io_check_array( TRIM( chem_species(lsp)%name )//'_av' , found = array_found )
3306       IF ( array_found )  THEN
3307          IF ( .NOT. ALLOCATED( chem_species(lsp)%conc_av ) )  THEN
3308             ALLOCATE( chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3309          ENDIF
3310          CALL rrd_mpi_io( TRIM( chem_species(lsp)%name )//'_av', chem_species(lsp)%conc_av )
3311       ENDIF
3312
3313    ENDDO
3314
3315 END SUBROUTINE chem_rrd_local_mpi
3316
3317
3318!--------------------------------------------------------------------------------------------------!
3319!> Description:
3320!> Calculation of horizontally averaged profiles
3321!> This routine is called for every statistic region (sr) defined by the user,
3322!> but at least for the region "total domain" (sr=0).
3323!> quantities.
3324!--------------------------------------------------------------------------------------------------!
3325 SUBROUTINE chem_statistics( mode, sr, tn )
3326
3327
3328    USE arrays_3d
3329
3330    USE statistics
3331
3332    CHARACTER (LEN=*) ::  mode   !<
3333
3334    INTEGER(iwp) ::  i                   !< running index on x-axis
3335    INTEGER(iwp) ::  j                   !< running index on y-axis
3336    INTEGER(iwp) ::  k                   !< vertical index counter
3337    INTEGER(iwp) ::  l                   !< direction index counter
3338    INTEGER(iwp) ::  lpr                 !< running index chem spcs
3339    INTEGER(iwp) ::  m                   !< running index for surface elements
3340!$  INTEGER(iwp) ::  omp_get_thread_num  !< intrinsic OMP function
3341    INTEGER(iwp) ::  sr                  !< statistical region
3342    INTEGER(iwp) ::  surf_e              !< end surface index
3343    INTEGER(iwp) ::  surf_s              !< start surface index
3344    INTEGER(iwp) ::  tn                  !< thread number
3345
3346    REAL(wp)                                            ::  flag     !< topography masking flag
3347    REAL(wp), DIMENSION(nzb:nzt+1,0:threads_per_task-1) ::  sums_tmp !< temporary array used to sum-up profiles
3348
3349    IF ( mode == 'profiles' )  THEN
3350!
3351!--    Sum-up profiles for the species
3352       tn = 0
3353       !$OMP PARALLEL PRIVATE( i, j, k, tn, lpr, sums_tmp )
3354       !$ tn = omp_get_thread_num()
3355       !$OMP DO
3356       DO  lpr = 1, cs_pr_count_sp
3357          sums_tmp(:,tn) = 0.0_wp
3358          DO  i = nxl, nxr
3359             DO  j = nys, nyn
3360                DO  k = nzb, nzt+1
3361                   sums_tmp(k,tn) = sums_tmp(k,tn) +                                               &
3362                        chem_species(cs_pr_index_sp(lpr))%conc(k,j,i) *                            &
3363                        rmask(j,i,sr)  *                                                           &
3364                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) )
3365                ENDDO
3366             ENDDO
3367          ENDDO
3368          sums_l(nzb:nzt+1,hom_index_spec(lpr),tn) = sums_tmp(nzb:nzt+1,tn)
3369       ENDDO
3370       !$OMP END PARALLEL
3371!
3372!--    Sum-up profiles for vertical fluxes of the the species. Note, in case of WS5 scheme the
3373!--    profiles of resolved-scale fluxes have been already summed-up, while resolved-scale fluxes
3374!--    need to be calculated in case of PW scheme.
3375!--    For summation employ a temporary array.
3376       !$OMP PARALLEL PRIVATE( i, j, k, tn, lpr, sums_tmp, flag )
3377       !$ tn = omp_get_thread_num()
3378       !$OMP DO
3379       DO  lpr = 1, cs_pr_count_fl_sgs
3380          sums_tmp(:,tn) = 0.0_wp
3381          DO  i = nxl, nxr
3382             DO  j = nys, nyn
3383                DO  k = nzb, nzt
3384                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 23 ) ) *        &
3385                          MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 9  ) )
3386                   sums_tmp(k,tn) = sums_tmp(k,tn) -                                               &
3387                        0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                                       &
3388                               * ( chem_species(cs_pr_index_fl_sgs(lpr))%conc(k+1,j,i) -           &
3389                                   chem_species(cs_pr_index_fl_sgs(lpr))%conc(k,j,i) ) *           &
3390                        ddzu(k+1) * rmask(j,i,sr)  * flag
3391                ENDDO
3392!
3393!--             Add surface fluxes (?Is the order mandatory or could it be done in one cycle?)
3394                DO l = 0, 1
3395                   surf_s = surf_def_h(0)%start_index(j,i)
3396                   surf_e = surf_def_h(0)%end_index(j,i)
3397                   DO  m = surf_s, surf_e
3398                      k   = surf_def_h(0)%k(m) + surf_def_h(0)%koff
3399                      sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
3400                   ENDDO
3401                ENDDO
3402                DO l = 0, 1
3403                   surf_s = surf_lsm_h(0)%start_index(j,i)
3404                   surf_e = surf_lsm_h(0)%end_index(j,i)
3405                   DO  m = surf_s, surf_e
3406                      k   = surf_lsm_h(0)%k(m) + surf_lsm_h(0)%koff
3407                      sums_tmp(k,tn) = sums_tmp(k,tn) + surf_lsm_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
3408                   ENDDO
3409                ENDDO
3410                DO l = 0, 1
3411                   surf_s = surf_usm_h(0)%start_index(j,i)
3412                   surf_e = surf_usm_h(0)%end_index(j,i)
3413                   DO  m = surf_s, surf_e
3414                      k   = surf_usm_h(0)%k(m) + surf_usm_h(0)%koff
3415                      sums_tmp(k,tn) = sums_tmp(k,tn) + surf_usm_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
3416                   ENDDO
3417                ENDDO
3418             ENDDO
3419          ENDDO
3420          sums_l(nzb:nzt+1,hom_index_fl_sgs(lpr),tn) = sums_tmp(nzb:nzt+1,tn)
3421       ENDDO
3422       !$OMP END PARALLEL
3423!
3424!--    Resolved-scale fluxes from the WS5 scheme
3425       IF ( ws_scheme_sca )  THEN
3426          !$OMP PARALLEL PRIVATE( tn, lpr )
3427          !$ tn = omp_get_thread_num()
3428          !$OMP DO
3429          DO  lpr = 1, cs_pr_count_fl_res
3430             sums_l(nzb:nzt+1,hom_index_fl_res(lpr),tn) =                                           &
3431                                                      sums_ws_l(nzb:nzt+1,tn,cs_pr_index_fl_res(lpr))
3432          ENDDO
3433          !$OMP END PARALLEL
3434       ENDIF
3435
3436    ELSEIF ( mode == 'time_series' )  THEN
3437!      @todo
3438    ENDIF
3439
3440 END SUBROUTINE chem_statistics
3441
3442
3443!--------------------------------------------------------------------------------------------------!
3444! Description:
3445! ------------
3446!> Subroutine for swapping of timelevels for chemical species called out from subroutine
3447!> swap_timelevel
3448!--------------------------------------------------------------------------------------------------!
3449
3450
3451 SUBROUTINE chem_swap_timelevel( level )
3452
3453
3454    INTEGER(iwp), INTENT(IN) ::  level
3455!
3456!-- Local variables
3457    INTEGER(iwp)             ::  lsp
3458
3459
3460    IF ( level == 0 )  THEN
3461       DO  lsp=1, nvar
3462          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
3463          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
3464       ENDDO
3465    ELSE
3466       DO  lsp=1, nvar
3467          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
3468          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
3469       ENDDO
3470    ENDIF
3471
3472    RETURN
3473 END SUBROUTINE chem_swap_timelevel
3474
3475
3476!--------------------------------------------------------------------------------------------------!
3477! Description:
3478! ------------
3479!> Subroutine to write restart data for chemistry model
3480!--------------------------------------------------------------------------------------------------!
3481 SUBROUTINE chem_wrd_local
3482
3483
3484    INTEGER(iwp) ::  lsp  !< running index for chem spcs.
3485
3486    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
3487
3488       DO  lsp = 1, nspec
3489          CALL wrd_write_string( TRIM( chem_species(lsp)%name ) )
3490          WRITE ( 14 )  chem_species(lsp)%conc
3491          IF ( ALLOCATED( chem_species(lsp)%conc_av ) )  THEN
3492             CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
3493             WRITE ( 14 )  chem_species(lsp)%conc_av
3494          ENDIF
3495       ENDDO
3496
3497    ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
3498
3499       DO  lsp = 1, nspec
3500          CALL wrd_mpi_io( TRIM( chem_species(lsp)%name ), chem_species(lsp)%conc )
3501          IF ( ALLOCATED( chem_species(lsp)%conc_av ) )  THEN
3502             CALL wrd_mpi_io( TRIM( chem_species(lsp)%name ) // '_av', chem_species(lsp)%conc_av )
3503          ENDIF
3504       ENDDO
3505
3506    ENDIF
3507
3508 END SUBROUTINE chem_wrd_local
3509
3510
3511!--------------------------------------------------------------------------------------------------!
3512! Description:
3513! ------------
3514!> Subroutine to calculate the deposition of gases and PMs. For now deposition only takes place on
3515!> lsm and usm horizontal upward faceing surfaces. Default surfaces are NOT considered.
3516!> The deposition of particlesis derived following Zhang et al., 2001, gases are deposited using
3517!>  the DEPAC module (van Zanten et al., 2010).
3518!>
3519!> @TODO: Consider deposition on vertical surfaces
3520!> @TODO: Consider overlaying horizontal surfaces
3521!> @TODO: Consider resolved vegetation
3522!> @TODO: Check error messages
3523!--------------------------------------------------------------------------------------------------!
3524 SUBROUTINE chem_depo( i, j )
3525
3526    USE control_parameters,                                                                        &
3527         ONLY:  dt_3d, intermediate_timestep_count, latitude, time_since_reference_point
3528
3529    USE arrays_3d,                                                                                 &
3530         ONLY:  dzw, rho_air_zw
3531
3532    USE palm_date_time_mod,                                                                        &
3533         ONLY:  get_date_time
3534
3535    USE surface_mod,                                                                               &
3536         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h, surf_type, surf_usm_h
3537
3538    USE radiation_model_mod,                                                                       &
3539         ONLY:  cos_zenith
3540
3541
3542    INTEGER(iwp) ::  day_of_year                   !< current day of the year
3543    INTEGER(iwp), INTENT(IN) ::  i
3544    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
3545    INTEGER(iwp), INTENT(IN) ::  j
3546    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
3547    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
3548    INTEGER(iwp) ::  luv_palm                      !< index of PALM LSM vegetation_type at current
3549                                                   !< surface element
3550    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current
3551                                                   !< surface element
3552    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current
3553                                                   !< surface element
3554    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current
3555                                                   !< surface element
3556    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current
3557                                                   !< surface element
3558    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface
3559                                                   !< element
3560    INTEGER(iwp) ::  luv_dep                       !< matching DEPAC LU to luv_palm
3561    INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
3562    INTEGER(iwp) ::  luw_dep                       !< matching DEPAC LU to luw_palm
3563    INTEGER(iwp) ::  luu_dep                       !< matching DEPAC LU to luu_palm
3564    INTEGER(iwp) ::  lug_dep                       !< matching DEPAC LU to lug_palm
3565    INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
3566    INTEGER(iwp) ::  m                             !< index for horizontal surfaces
3567    INTEGER(iwp) ::  pspec                         !< running index
3568!
3569!-- Vegetation                                               !<  Assign PALM classes to DEPAC land
3570                                                             !<  use classes
3571    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM
3572    INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
3573    INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
3574    INTEGER(iwp) ::  ind_luv_s_grass = 3                     !<  assigned to ilu_grass
3575    INTEGER(iwp) ::  ind_luv_ev_needle_trees = 4             !<  assigned to ilu_coniferous_forest
3576    INTEGER(iwp) ::  ind_luv_de_needle_trees = 5             !<  assigned to ilu_coniferous_forest
3577    INTEGER(iwp) ::  ind_luv_ev_broad_trees = 6              !<  assigned to ilu_tropical_forest
3578    INTEGER(iwp) ::  ind_luv_de_broad_trees = 7              !<  assigned to ilu_deciduous_forest
3579    INTEGER(iwp) ::  ind_luv_t_grass = 8                     !<  assigned to ilu_grass
3580    INTEGER(iwp) ::  ind_luv_desert = 9                      !<  assigned to ilu_desert
3581    INTEGER(iwp) ::  ind_luv_tundra = 10                     !<  assigned to ilu_other
3582    INTEGER(iwp) ::  ind_luv_irr_crops = 11                  !<  assigned to ilu_arable
3583    INTEGER(iwp) ::  ind_luv_semidesert = 12                 !<  assigned to ilu_other
3584    INTEGER(iwp) ::  ind_luv_ice = 13                        !<  assigned to ilu_ice
3585    INTEGER(iwp) ::  ind_luv_marsh = 14                      !<  assigned to ilu_other
3586    INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
3587    INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
3588    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest(ave(decid+conif))
3589    INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to ilu_other (ave(other+decid))
3590!
3591!-- Water
3592    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM
3593    INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
3594    INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
3595    INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
3596    INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
3597    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland
3598!
3599!-- Pavement
3600    INTEGER(iwp) ::  ind_lup_user = 0                        !<  ERROR as no class given in PALM
3601    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !<  assigned to ilu_desert
3602    INTEGER(iwp) ::  ind_lup_asph = 2                        !<  assigned to ilu_desert
3603    INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
3604    INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
3605    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert
3606    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
3607    INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
3608    INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
3609    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert
3610    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
3611    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
3612    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !<  assigned to ilu_desert
3613    INTEGER(iwp) ::  ind_lup_tartan = 13                     !<  assigned to ilu_desert
3614    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !<  assigned to ilu_desert
3615    INTEGER(iwp) ::  ind_lup_clay = 15                       !<  assigned to ilu_desert
3616!
3617!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
3618    INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
3619    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
3620    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
3621    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
3622    INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
3623
3624
3625    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
3626    REAL(wp) ::  dens              !< density at layer k at i,j
3627    REAL(wp) ::  dh                !< vertical grid size
3628    REAL(wp) ::  diffusivity       !< diffusivity
3629    REAL(wp) ::  dt_chem           !< length of chem time step
3630    REAL(wp) ::  dt_dh             !< dt_chem/dh
3631    REAL(wp) ::  inv_dh            !< inverse of vertical grid size
3632    REAL(wp) ::  lai               !< leaf area index at current surface element
3633    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
3634    REAL(wp) ::  qv_tmp            !< surface mixing ratio at current surface element
3635    REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
3636    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
3637    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
3638    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
3639    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
3640    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be
3641                                   !< lai + 1
3642    REAL(wp) ::  slinnfac
3643    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface
3644                                   !< element
3645    REAL(wp) ::  temp_tmp          !< temperatur at i,j,k
3646    REAL(wp) ::  ts                !< surface temperatur in degrees celsius
3647    REAL(wp) ::  ustar_surf        !< ustar at current surface element
3648    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
3649    REAL(wp) ::  visc              !< Viscosity
3650    REAL(wp) ::  vs                !< Sedimentation velocity
3651    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
3652
3653    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
3654    REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
3655    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface
3656                                                !< element
3657    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface
3658                                                !< element
3659    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface
3660                                                !< element
3661    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface
3662                                                !< element
3663    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface
3664                                                !< element
3665    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to
3666                                                !< zero for all species!
3667    REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
3668
3669
3670!
3671!-- Particle parameters (PM10 (1), PM25 (2)) partsize (diameter in m), rhopart (density in kg/m3),
3672!-- slipcor (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
3673    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/                        &
3674                                                         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &   !<  1
3675                                                         0.7e-6_wp, 1.14e3_wp, 1.082_wp  &   !<  2
3676                                                         /), (/ 3, 2 /) )
3677
3678    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
3679    LOGICAL ::  match_usm     !< flag indicating urban-type surface
3680
3681!
3682!-- List of names of possible tracers
3683    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/                                         &
3684         'NO2           ', &    !< NO2
3685         'NO            ', &    !< NO
3686         'O3            ', &    !< O3
3687         'CO            ', &    !< CO
3688         'form          ', &    !< FORM
3689         'ald           ', &    !< ALD
3690         'pan           ', &    !< PAN
3691         'mgly          ', &    !< MGLY
3692         'par           ', &    !< PAR
3693         'ole           ', &    !< OLE
3694         'eth           ', &    !< ETH
3695         'tol           ', &    !< TOL
3696         'cres          ', &    !< CRES
3697         'xyl           ', &    !< XYL
3698         'SO4a_f        ', &    !< SO4a_f
3699         'SO2           ', &    !< SO2
3700         'HNO2          ', &    !< HNO2
3701         'CH4           ', &    !< CH4
3702         'NH3           ', &    !< NH3
3703         'NO3           ', &    !< NO3
3704         'OH            ', &    !< OH
3705         'HO2           ', &    !< HO2
3706         'N2O5          ', &    !< N2O5
3707         'SO4a_c        ', &    !< SO4a_c
3708         'NH4a_f        ', &    !< NH4a_f
3709         'NO3a_f        ', &    !< NO3a_f
3710         'NO3a_c        ', &    !< NO3a_c
3711         'C2O3          ', &    !< C2O3
3712         'XO2           ', &    !< XO2
3713         'XO2N          ', &    !< XO2N
3714         'cro           ', &    !< CRO
3715         'HNO3          ', &    !< HNO3
3716         'H2O2          ', &    !< H2O2
3717         'iso           ', &    !< ISO
3718         'ispd          ', &    !< ISPD
3719         'to2           ', &    !< TO2
3720         'open          ', &    !< OPEN
3721         'terp          ', &    !< TERP
3722         'ec_f          ', &    !< EC_f
3723         'ec_c          ', &    !< EC_c
3724         'pom_f         ', &    !< POM_f
3725         'pom_c         ', &    !< POM_c
3726         'ppm_f         ', &    !< PPM_f
3727         'ppm_c         ', &    !< PPM_c
3728         'na_ff         ', &    !< Na_ff
3729         'na_f          ', &    !< Na_f
3730         'na_c          ', &    !< Na_c
3731         'na_cc         ', &    !< Na_cc
3732         'na_ccc        ', &    !< Na_ccc
3733         'dust_ff       ', &    !< dust_ff
3734         'dust_f        ', &    !< dust_f
3735         'dust_c        ', &    !< dust_c
3736         'dust_cc       ', &    !< dust_cc
3737         'dust_ccc      ', &    !< dust_ccc
3738         'tpm10         ', &    !< tpm10
3739         'tpm25         ', &    !< tpm25
3740         'tss           ', &    !< tss
3741         'tdust         ', &    !< tdust
3742         'tc            ', &    !< tc
3743         'tcg           ', &    !< tcg
3744         'tsoa          ', &    !< tsoa
3745         'tnmvoc        ', &    !< tnmvoc
3746         'SOxa          ', &    !< SOxa
3747         'NOya          ', &    !< NOya
3748         'NHxa          ', &    !< NHxa
3749         'NO2_obs       ', &    !< NO2_obs
3750         'tpm10_biascorr', &    !< tpm10_biascorr
3751         'tpm25_biascorr', &    !< tpm25_biascorr
3752         'O3_biascorr   ' /)    !< o3_biascorr
3753!
3754!-- Tracer mole mass:
3755    REAL(wp), PARAMETER ::  specmolm(nposp) = (/                                                   &
3756         xm_O * 2 + xm_N, &                         !< NO2
3757         xm_O + xm_N, &                             !< NO
3758         xm_O * 3, &                                !< O3
3759         xm_C + xm_O, &                             !< CO
3760         xm_H * 2 + xm_C + xm_O, &                  !< FORM
3761         xm_H * 3 + xm_C * 2 + xm_O, &              !< ALD
3762         xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   !< PAN
3763         xm_H * 4 + xm_C * 3 + xm_O * 2, &          !< MGLY
3764         xm_H * 3 + xm_C, &                         !< PAR
3765         xm_H * 3 + xm_C * 2, &                     !< OLE
3766         xm_H * 4 + xm_C * 2, &                     !< ETH
3767         xm_H * 8 + xm_C * 7, &                     !< TOL
3768         xm_H * 8 + xm_C * 7 + xm_O, &              !< CRES
3769         xm_H * 10 + xm_C * 8, &                    !< XYL
3770         xm_S + xm_O * 4, &                         !< SO4a_f
3771         xm_S + xm_O * 2, &                         !< SO2
3772         xm_H + xm_O * 2 + xm_N, &                  !< HNO2
3773         xm_H * 4 + xm_C, &                         !< CH4
3774         xm_H * 3 + xm_N, &                         !< NH3
3775         xm_O * 3 + xm_N, &                         !< NO3
3776         xm_H + xm_O, &                             !< OH
3777         xm_H + xm_O * 2, &                         !< HO2
3778         xm_O * 5 + xm_N * 2, &                     !< N2O5
3779         xm_S + xm_O * 4, &                         !< SO4a_c
3780         xm_H * 4 + xm_N, &                         !< NH4a_f
3781         xm_O * 3 + xm_N, &                         !< NO3a_f
3782         xm_O * 3 + xm_N, &                         !< NO3a_c
3783         xm_C * 2 + xm_O * 3, &                     !< C2O3
3784         xm_dummy, &                                !< XO2
3785         xm_dummy, &                                !< XO2N
3786         xm_dummy, &                                !< CRO
3787         xm_H + xm_O * 3 + xm_N, &                  !< HNO3
3788         xm_H * 2 + xm_O * 2, &                     !< H2O2
3789         xm_H * 8 + xm_C * 5, &                     !< ISO
3790         xm_dummy, &                                !< ISPD
3791         xm_dummy, &                                !< TO2
3792         xm_dummy, &                                !< OPEN
3793         xm_H * 16 + xm_C * 10, &                   !< TERP
3794         xm_dummy, &                                !< EC_f
3795         xm_dummy, &                                !< EC_c
3796         xm_dummy, &                                !< POM_f
3797         xm_dummy, &                                !< POM_c
3798         xm_dummy, &                                !< PPM_f
3799         xm_dummy, &                                !< PPM_c
3800         xm_Na, &                                   !< Na_ff
3801         xm_Na, &                                   !< Na_f
3802         xm_Na, &                                   !< Na_c
3803         xm_Na, &                                   !< Na_cc
3804         xm_Na, &                                   !< Na_ccc
3805         xm_dummy, &                                !< dust_ff
3806         xm_dummy, &                                !< dust_f
3807         xm_dummy, &                                !< dust_c
3808         xm_dummy, &                                !< dust_cc
3809         xm_dummy, &                                !< dust_ccc
3810         xm_dummy, &                                !< tpm10
3811         xm_dummy, &                                !< tpm25
3812         xm_dummy, &                                !< tss
3813         xm_dummy, &                                !< tdust
3814         xm_dummy, &                                !< tc
3815         xm_dummy, &                                !< tcg
3816         xm_dummy, &                                !< tsoa
3817         xm_dummy, &                                !< tnmvoc
3818         xm_dummy, &                                !< SOxa
3819         xm_dummy, &                                !< NOya
3820         xm_dummy, &                                !< NHxa
3821         xm_O * 2 + xm_N, &                         !< NO2_obs
3822         xm_dummy, &                                !< tpm10_biascorr
3823         xm_dummy, &                                !< tpm25_biascorr
3824         xm_O * 3 /)                                !< o3_biascorr
3825!
3826!-- Get current day of the year
3827    CALL get_date_time( time_since_reference_point, day_of_year = day_of_year )
3828!
3829!-- Initialize surface element m
3830    m = 0
3831    k = 0
3832!
3833!-- LSM or USM horizintal upward facing surface present at i,j:
3834!-- Default surfaces are NOT considered for deposition
3835    match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
3836    match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
3837!
3838!--For LSM surfaces
3839    IF ( match_lsm )  THEN
3840!
3841!--    Get surface element information at i,j:
3842       m = surf_lsm_h(0)%start_index(j,i)
3843       k = surf_lsm_h(0)%k(m)
3844!
3845!--    Get needed variables for surface element m
3846       ustar_surf  = surf_lsm_h(0)%us(m)
3847       z0h_surf    = surf_lsm_h(0)%z0h(m)
3848       r_aero_surf = surf_lsm_h(0)%r_a(m)
3849       solar_rad   = surf_lsm_h(0)%rad_sw_dir(m) + surf_lsm_h(0)%rad_sw_dif(m)
3850
3851       lai = surf_lsm_h(0)%lai(m)
3852       sai = lai + 1
3853!
3854!--    For small grid spacing neglect R_a
3855       IF ( dzw(k) <= 1.0 )  THEN
3856          r_aero_surf = 0.0_wp
3857       ENDIF
3858!
3859!--    Initialize lu's
3860       luv_palm = 0
3861       luv_dep = 0
3862       lup_palm = 0
3863       lup_dep = 0
3864       luw_palm = 0
3865       luw_dep = 0
3866!
3867!--    Initialize budgets
3868       bud_luv = 0.0_wp
3869       bud_lup = 0.0_wp
3870       bud_luw = 0.0_wp
3871!
3872!--    Get land use for i,j and assign to DEPAC lu
3873       IF ( surf_lsm_h(0)%frac(m,ind_veg_wall) > 0 )  THEN
3874          luv_palm = surf_lsm_h(0)%vegetation_type(m)
3875          IF ( luv_palm == ind_luv_user )  THEN
3876             message_string = 'No lsm-vegetation type defined. Please define vegetation type to' //&
3877                              ' enable deposition calculation'
3878             CALL message( 'chem_depo', 'PA0743', 1, 2, 0, 6, 0 )
3879          ELSEIF ( luv_palm == ind_luv_b_soil )  THEN
3880             luv_dep = 9
3881          ELSEIF ( luv_palm == ind_luv_mixed_crops )  THEN
3882             luv_dep = 2
3883          ELSEIF ( luv_palm == ind_luv_s_grass )  THEN
3884             luv_dep = 1
3885          ELSEIF ( luv_palm == ind_luv_ev_needle_trees )  THEN
3886             luv_dep = 4
3887          ELSEIF ( luv_palm == ind_luv_de_needle_trees )  THEN
3888             luv_dep = 4
3889          ELSEIF ( luv_palm == ind_luv_ev_broad_trees )  THEN
3890             luv_dep = 12
3891          ELSEIF ( luv_palm == ind_luv_de_broad_trees )  THEN
3892             luv_dep = 5
3893          ELSEIF ( luv_palm == ind_luv_t_grass )  THEN
3894             luv_dep = 1
3895          ELSEIF ( luv_palm == ind_luv_desert )  THEN
3896             luv_dep = 9
3897          ELSEIF ( luv_palm == ind_luv_tundra )  THEN
3898             luv_dep = 8
3899          ELSEIF ( luv_palm == ind_luv_irr_crops )  THEN
3900             luv_dep = 2
3901          ELSEIF ( luv_palm == ind_luv_semidesert )  THEN
3902             luv_dep = 8
3903          ELSEIF ( luv_palm == ind_luv_ice )  THEN
3904             luv_dep = 10
3905          ELSEIF ( luv_palm == ind_luv_marsh )  THEN
3906             luv_dep = 8
3907          ELSEIF ( luv_palm == ind_luv_ev_shrubs )  THEN
3908             luv_dep = 14
3909          ELSEIF ( luv_palm == ind_luv_de_shrubs )  THEN
3910             luv_dep = 14
3911          ELSEIF ( luv_palm == ind_luv_mixed_forest )  THEN
3912             luv_dep = 4
3913          ELSEIF ( luv_palm == ind_luv_intrup_forest )  THEN
3914             luv_dep = 8
3915          ENDIF
3916       ENDIF
3917
3918       IF ( surf_lsm_h(0)%frac(m,ind_pav_green) > 0 )  THEN
3919          lup_palm = surf_lsm_h(0)%pavement_type(m)
3920          IF ( lup_palm == ind_lup_user )  THEN
3921             message_string = 'No lsm-pavement type defined. Please define pavement type to ' //   &
3922                              'enable deposition calculation'
3923             CALL message( 'chem_depo', 'PA0743', 1, 2, 0, 6, 0 )
3924          ELSEIF ( lup_palm == ind_lup_asph_conc )  THEN
3925             lup_dep = 9
3926          ELSEIF ( lup_palm == ind_lup_asph )  THEN
3927             lup_dep = 9
3928          ELSEIF ( lup_palm == ind_lup_conc )  THEN
3929             lup_dep = 9
3930          ELSEIF ( lup_palm == ind_lup_sett )  THEN
3931             lup_dep = 9
3932          ELSEIF ( lup_palm == ind_lup_pav_stones )  THEN
3933             lup_dep = 9
3934          ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
3935             lup_dep = 9
3936          ELSEIF ( lup_palm == ind_lup_metal )  THEN
3937             lup_dep = 9
3938          ELSEIF ( lup_palm == ind_lup_wood )  THEN
3939             lup_dep = 9
3940          ELSEIF ( lup_palm == ind_lup_gravel )  THEN
3941             lup_dep = 9
3942          ELSEIF ( lup_palm == ind_lup_f_gravel )  THEN
3943             lup_dep = 9
3944          ELSEIF ( lup_palm == ind_lup_pebblest )  THEN
3945             lup_dep = 9
3946          ELSEIF ( lup_palm == ind_lup_woodchips )  THEN
3947             lup_dep = 9
3948          ELSEIF ( lup_palm == ind_lup_tartan )  THEN
3949             lup_dep = 9
3950          ELSEIF ( lup_palm == ind_lup_art_turf )  THEN
3951             lup_dep = 9
3952          ELSEIF ( lup_palm == ind_lup_clay )  THEN
3953             lup_dep = 9
3954          ENDIF
3955       ENDIF
3956
3957       IF ( surf_lsm_h(0)%frac(m,ind_wat_win) > 0 )  THEN
3958          luw_palm = surf_lsm_h(0)%water_type(m)
3959          IF ( luw_palm == ind_luw_user )  THEN
3960             message_string = 'No lsm-water type defined. Please define water type to enable' //   &
3961                              ' deposition calculation'
3962             CALL message( 'chem_depo', 'PA0743', 1, 2, 0, 6, 0 )
3963          ELSEIF ( luw_palm ==  ind_luw_lake )  THEN
3964             luw_dep = 13
3965          ELSEIF ( luw_palm == ind_luw_river )  THEN
3966             luw_dep = 13
3967          ELSEIF ( luw_palm == ind_luw_ocean )  THEN
3968             luw_dep = 6
3969          ELSEIF ( luw_palm == ind_luw_pond )  THEN
3970             luw_dep = 13
3971          ELSEIF ( luw_palm == ind_luw_fountain )  THEN
3972             luw_dep = 13
3973          ENDIF
3974       ENDIF
3975!
3976!--    Set wetness indicator to dry or wet for lsm vegetation or pavement
3977       IF ( surf_lsm_h(0)%c_liq(m) > 0 )  THEN
3978          nwet = 1
3979       ELSE
3980          nwet = 0
3981       ENDIF
3982!
3983!--    Compute length of time step
3984       IF ( call_chem_at_all_substeps )  THEN
3985          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3986       ELSE
3987          dt_chem = dt_3d
3988       ENDIF
3989
3990       dh = dzw(k)
3991       inv_dh = 1.0_wp / dh
3992       dt_dh = dt_chem/dh
3993!
3994!--    Concentration at i,j,k
3995       DO  lsp = 1, nspec
3996          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
3997       ENDDO
3998
3999!--    Temperature at i,j,k
4000       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
4001
4002       ts       = temp_tmp - 273.15  !< in degrees celcius
4003!
4004!--    Viscosity of air
4005       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
4006!
4007!--    Air density at k
4008       dens = rho_air_zw(k)
4009!
4010!--    Calculate relative humidity from specific humidity for DEPAC
4011       qv_tmp = MAX( q(k,j,i), 0.0_wp)
4012       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
4013!
4014!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
4015!-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
4016!
4017!--    Vegetation
4018       IF ( surf_lsm_h(0)%frac(m,ind_veg_wall) > 0 )  THEN
4019
4020!
4021!--       No vegetation on bare soil, desert or ice:
4022          IF ( ( luv_palm == ind_luv_b_soil )  .OR.  ( luv_palm == ind_luv_desert )  .OR.          &
4023               ( luv_palm == ind_luv_ice ) ) THEN
4024
4025             lai = 0.0_wp
4026             sai = 0.0_wp
4027
4028          ENDIF
4029
4030          slinnfac = 1.0_wp
4031!
4032!--       Get deposition velocity vd
4033          DO  lsp = 1, nvar
4034!
4035!--          Initialize
4036             vs = 0.0_wp
4037             vd_lu = 0.0_wp
4038             rs = 0.0_wp
4039             rb = 0.0_wp
4040             rc_tot = 0.0_wp
4041             IF ( spc_names(lsp) == 'PM10' )  THEN
4042                part_type = 1
4043!
4044!--             Sedimentation velocity
4045                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
4046                                                        particle_pars(ind_p_size, part_type),      &
4047                                                        particle_pars(ind_p_slip, part_type),      &
4048                                                        visc)
4049
4050                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
4051                                            vs,                                                    &
4052                                            particle_pars(ind_p_size, part_type),                  &
4053                                            particle_pars(ind_p_slip, part_type),                  &
4054                                            nwet, temp_tmp, dens, visc,                            &
4055                                            luv_dep,                                               &
4056                                            r_aero_surf, ustar_surf )
4057
4058                bud_luv(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
4059
4060
4061             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4062                part_type = 2
4063!
4064!--             Sedimentation velocity
4065                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
4066                                                        particle_pars(ind_p_size, part_type),      &
4067                                                        particle_pars(ind_p_slip, part_type),      &
4068                                                        visc)
4069
4070                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
4071                                            vs,                                                    &
4072                                            particle_pars(ind_p_size, part_type),                  &
4073                                            particle_pars(ind_p_slip, part_type),                  &
4074                                            nwet, temp_tmp, dens, visc,                            &
4075                                            luv_dep ,                                              &
4076                                            r_aero_surf, ustar_surf )
4077
4078                bud_luv(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
4079
4080             ELSE !< GASES
4081!
4082!--             Read spc_name of current species for gas parameter
4083                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4084                   i_pspec = 0
4085                   DO  pspec = 1, nposp
4086                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4087                         i_pspec = pspec
4088                      END IF
4089                   ENDDO
4090
4091                ELSE
4092!
4093!--             For now species not deposited
4094                   CYCLE
4095                ENDIF
4096!
4097!--             Factor used for conversion from ppb to ug/m3 :
4098!--             ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4099!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4100!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
4101!--             thus:
4102!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4103!--             Use density at k:
4104                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
4105!
4106!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4107                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4108                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4109!
4110!--             Diffusivity for DEPAC relevant gases
4111!--             Use default value
4112                diffusivity            = 0.11e-4
4113!
4114!--             Overwrite with known coefficients of diffusivity from Massman (1998)
4115                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
4116                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4117                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4118                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4119                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4120                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4121                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4122!
4123!--             Get quasi-laminar boundary layer resistance rb:
4124                CALL get_rb_cell( ( luv_dep == ilu_water_sea )  .OR.                               &
4125                                  ( luv_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
4126                                  diffusivity, rb )
4127!
4128!--             Get rc_tot
4129                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
4130                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luv_dep,  &
4131                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
4132                                         diffusivity, r_aero_surf , rb )
4133!
4134!--             Calculate budget
4135                IF ( rc_tot <= 0.0 )  THEN
4136
4137                   bud_luv(lsp) = 0.0_wp
4138
4139                ELSE
4140
4141                   vd_lu = 1.0_wp / ( r_aero_surf + rb + rc_tot )
4142
4143                   bud_luv(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
4144                                    ( 1.0_wp - EXP( -vd_lu * dt_dh ) ) * dh
4145                ENDIF
4146
4147             ENDIF
4148          ENDDO
4149       ENDIF
4150!
4151!--    Pavement
4152       IF ( surf_lsm_h(0)%frac(m,ind_pav_green) > 0 )  THEN
4153!
4154!--       No vegetation on pavements:
4155          lai = 0.0_wp
4156          sai = 0.0_wp
4157
4158          slinnfac = 1.0_wp
4159!
4160!--       Get vd
4161          DO  lsp = 1, nvar
4162!
4163!--       Initialize
4164             vs = 0.0_wp
4165             vd_lu = 0.0_wp
4166             rs = 0.0_wp
4167             rb = 0.0_wp
4168             rc_tot = 0.0_wp
4169             IF ( spc_names(lsp) == 'PM10' )  THEN
4170                part_type = 1
4171!
4172!--             Sedimentation velocity:
4173                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
4174                                                        particle_pars(ind_p_size, part_type),      &
4175                                                        particle_pars(ind_p_slip, part_type),      &
4176                                                        visc)
4177
4178                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
4179                                            vs,                                                    &
4180                                            particle_pars(ind_p_size, part_type),                  &
4181                                            particle_pars(ind_p_slip, part_type),                  &
4182                                            nwet, temp_tmp, dens, visc,                            &
4183                                            lup_dep,                                               &
4184                                            r_aero_surf, ustar_surf )
4185
4186                bud_lup(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
4187
4188
4189             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4190                part_type = 2
4191!
4192!--             Sedimentation velocity:
4193                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
4194                                                        particle_pars(ind_p_size, part_type),      &
4195                                                        particle_pars(ind_p_slip, part_type),      &
4196                                                        visc)
4197
4198                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
4199                                            vs,                                                    &
4200                                            particle_pars(ind_p_size, part_type),                  &
4201                                            particle_pars(ind_p_slip, part_type),                  &
4202                                            nwet, temp_tmp, dens, visc,                            &
4203                                            lup_dep,                                               &
4204                                            r_aero_surf, ustar_surf )
4205
4206                bud_lup(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
4207
4208             ELSE  !<GASES
4209!
4210!--             Read spc_name of current species for gas parameter
4211                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4212                   i_pspec = 0
4213                   DO  pspec = 1, nposp
4214                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4215                         i_pspec = pspec
4216                      END IF
4217                   ENDDO
4218
4219                ELSE
4220!
4221!--                For now species not deposited
4222                   CYCLE
4223                ENDIF
4224!
4225!--             Factor used for conversion from ppb to ug/m3 :
4226!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4227!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4228!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4229!--             thus:
4230!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4231!--             Use density at lowest layer:
4232                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
4233!
4234!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4235                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4236                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4237!
4238!--             Diffusivity for DEPAC relevant gases
4239!--             Use default value
4240                diffusivity            = 0.11e-4
4241!
4242!--             Overwrite with known coefficients of diffusivity from Massman (1998)
4243                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
4244                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4245                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4246                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4247                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4248                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4249                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4250!
4251!--             Get quasi-laminar boundary layer resistance rb:
4252                CALL get_rb_cell( ( lup_dep == ilu_water_sea )  .OR.                               &
4253                                  ( lup_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
4254                                    diffusivity, rb )
4255!
4256!--             Get rc_tot
4257                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
4258                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep,  &
4259                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
4260                                         diffusivity, r_aero_surf , rb )
4261!
4262!--             Calculate budget
4263                IF ( rc_tot <= 0.0 )  THEN
4264                   bud_lup(lsp) = 0.0_wp
4265                ELSE
4266                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4267                   bud_lup(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
4268                                  ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
4269                ENDIF
4270
4271             ENDIF
4272          ENDDO
4273       ENDIF
4274!
4275!--    Water
4276       IF ( surf_lsm_h(0)%frac(m,ind_wat_win) > 0 )  THEN
4277!
4278!--       No vegetation on water:
4279          lai = 0.0_wp
4280          sai = 0.0_wp
4281          slinnfac = 1.0_wp
4282!
4283!--       Get vd
4284          DO  lsp = 1, nvar
4285!
4286!--          Initialize
4287             vs = 0.0_wp
4288             vd_lu = 0.0_wp
4289             rs = 0.0_wp
4290             rb = 0.0_wp
4291             rc_tot = 0.0_wp
4292             IF ( spc_names(lsp) == 'PM10' )  THEN
4293                part_type = 1
4294!
4295!--             Sedimentation velocity:
4296                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
4297                                                        particle_pars(ind_p_size, part_type),      &
4298                                                        particle_pars(ind_p_slip, part_type),      &
4299                                                        visc)
4300
4301                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
4302                                            vs,                                                    &
4303                                            particle_pars(ind_p_size, part_type),                  &
4304                                            particle_pars(ind_p_slip, part_type),                  &
4305                                            nwet, temp_tmp, dens, visc,                            &
4306                                            luw_dep,                                               &
4307                                            r_aero_surf, ustar_surf )
4308
4309                bud_luw(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
4310
4311             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4312                part_type = 2
4313!
4314!--             Sedimentation velocity:
4315                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type),      &
4316                                                        particle_pars(ind_p_size, part_type), &
4317                                                        particle_pars(ind_p_slip, part_type), &
4318                                                        visc)
4319
4320                CALL drydepo_aero_zhang_vd( vd_lu, rs,                                             &
4321                                            vs,                                                    &
4322                                            particle_pars(ind_p_size, part_type),                  &
4323                                            particle_pars(ind_p_slip, part_type),                  &
4324                                            nwet, temp_tmp, dens, visc,                            &
4325                                            luw_dep,                                               &
4326                                            r_aero_surf, ustar_surf )
4327
4328                bud_luw(lsp) = - conc_ijk(lsp) * ( 1.0_wp - EXP( - vd_lu * dt_dh ) ) * dh
4329
4330             ELSE  !<GASES
4331!
4332!--             Read spc_name of current species for gas PARAMETER
4333                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
4334                   i_pspec = 0
4335                   DO  pspec = 1, nposp
4336                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4337                         i_pspec = pspec
4338                      END IF
4339                   ENDDO
4340
4341                ELSE
4342!
4343!--                For now species not deposited
4344                   CYCLE
4345                ENDIF
4346!
4347!--             Factor used for conversion from ppb to ug/m3 :
4348!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4349!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4350!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4351!--             thus:
4352!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4353!--             Use density at lowest layer:
4354                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
4355!
4356!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4357!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4358                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
4359!
4360!--             Diffusivity for DEPAC relevant gases
4361!--             Use default value
4362                diffusivity            = 0.11e-4
4363!
4364!--             Overwrite with known coefficients of diffusivity from Massman (1998)
4365                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4
4366                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4367                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4368                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4369                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4370                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4371                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4372!
4373!--             Get quasi-laminar boundary layer resistance rb:
4374                CALL get_rb_cell( ( luw_dep == ilu_water_sea )  .OR.                               &
4375                                  ( luw_dep == ilu_water_inland ), z0h_surf, ustar_surf,           &
4376                                    diffusivity, rb )
4377
4378!--             Get rc_tot
4379                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,    &
4380                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep,  &
4381                                         2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3,       &
4382                                         diffusivity, r_aero_surf , rb )
4383!
4384!--             Calculate budget
4385                IF ( rc_tot <= 0.0 )  THEN
4386
4387                   bud_luw(lsp) = 0.0_wp
4388
4389                ELSE
4390
4391                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4392
4393                   bud_luw(lsp) = - ( conc_ijk(lsp) - ccomp_tot(lsp) ) *                           &
4394                                  ( 1.0_wp - EXP(-vd_lu * dt_dh ) ) * dh
4395                ENDIF
4396
4397             ENDIF
4398          ENDDO
4399       ENDIF
4400
4401
4402       bud = 0.0_wp
4403!
4404!--    Calculate overall budget for surface m and adapt concentration
4405       DO  lsp = 1, nspec
4406
4407          bud(lsp) = surf_lsm_h(0)%frac(m,ind_veg_wall)  * bud_luv(lsp) +                             &
4408                     surf_lsm_h(0)%frac(m,ind_pav_green) * bud_lup(lsp) +                             &
4409                     surf_lsm_h(0)%frac(m,ind_wat_win)   * bud_luw(lsp)
4410!
4411!--       Compute new concentration:
4412          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
4413
4414          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
4415
4416       ENDDO
4417
4418    ENDIF
4419!
4420!-- For USM surfaces
4421    IF ( match_usm )  THEN
4422!
4423!--    Get surface element information at i,j:
4424       m = surf_usm_h(0)%start_index(j,i)
4425       k = surf_usm_h(0)%k(m)
4426!
4427!--    Get needed variables for surface element m
4428       ustar_surf  = surf_usm_h(0)%us(m)
4429       z0h_surf    = surf_usm_h(0)%z0h(m)
4430       r_aero_surf = surf_usm_h(0)%r_a(m)
4431       solar_rad   = surf_usm_h(0)%rad_sw_dir(m) + surf_usm_h(0)%rad_sw_dif(m)
4432       lai = surf_usm_h(0)%lai(m)
4433       sai = lai + 1
4434!
4435!--    For small grid spacing neglect R_a
4436       IF ( dzw(k) <= 1.0 )  THEN
4437          r_aero_surf = 0.0_wp
4438       ENDIF
4439!
4440!--    Initialize lu's
4441       luu_palm = 0
4442       luu_dep = 0
4443       lug_palm = 0
4444       lug_dep = 0
4445       lud_palm = 0
4446       lud_dep = 0
4447!
4448!--    Initialize budgets
4449       bud_luu  = 0.0_wp
4450       bud_lug = 0.0_wp
4451       bud_lud = 0.0_wp
4452!
4453!--    Get land use for i,j and assign to DEPAC lu
4454       IF ( surf_usm_h(0)%frac(m,ind_pav_green) > 0 )  THEN
4455!
4456!--       For green urban surfaces (e.g. green roofs assume LU short grass
4457          lug_palm = ind_luv_s_grass
4458          IF ( lug_palm == ind_luv_user )  THEN
4459             message_string = 'No usm-vegetation type defined. Please define vegetation type to' //&
4460                              ' enable deposition calculation'
4461             CALL message( 'chem_depo', 'PA0743', 1, 2, 0, 6, 0 )
4462          ELSEIF ( lug_palm == ind_luv_b_soil )  THEN
4463             lug_dep = 9
4464          ELSEIF ( lug_palm == ind_luv_mixed_crops )  THEN
4465             lug_dep = 2
4466          ELSEIF ( lug_palm == ind_luv_s_grass )  THEN
4467             lug_dep = 1
4468          ELSEIF ( lug_palm == ind_luv_ev_needle_trees )  THEN
4469             lug_dep = 4
4470          ELSEIF ( lug_palm == ind_luv_de_needle_trees )  THEN
4471             lug_dep = 4
4472          ELSEIF ( lug_palm == ind_luv_ev_broad_trees )  THEN
4473             lug_dep = 12
4474          ELSEIF ( lug_palm == ind_luv_de_broad_trees )  THEN
4475             lug_dep = 5
4476          ELSEIF ( lug_palm == ind_luv_t_grass )  THEN
4477             lug_dep = 1
4478          ELSEIF ( lug_palm == ind_luv_desert )  THEN
4479             lug_dep = 9
4480          ELSEIF ( lug_palm == ind_luv_tundra  )  THEN
4481             lug_dep = 8
4482          ELSEIF ( lug_palm == ind_luv_irr_crops )  THEN
4483             lug_dep = 2
4484          ELSEIF ( lug_palm == ind_luv_semidesert )  THEN
4485             lug_dep = 8
4486          ELSEIF ( lug_palm == ind_luv_ice )  THEN
4487             lug_dep = 10
4488          ELSEIF ( lug_palm == ind_luv_marsh )  THEN
4489             lug_dep = 8
4490          ELSEIF ( lug_palm == ind_luv_ev_shrubs )  THEN
4491             lug_dep = 14
4492          ELSEIF ( lug_palm == ind_luv_de_shrubs  )  THEN
4493             lug_dep = 14
4494          ELSEIF ( lug_palm == ind_luv_mixed_forest )  THEN
4495             lug_dep = 4
4496          ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
4497             lug_dep = 8
4498          ENDIF
4499       ENDIF
4500
4501       IF ( surf_usm_h(0)%frac(m,ind_veg_wall) > 0 )  THEN
4502!
4503!--       For walls in USM assume concrete walls/roofs,
4504!--       assumed LU class desert as also assumed for pavements in LSM
4505          luu_palm = ind_lup_conc
4506          IF ( luu_palm == ind_lup_user )  THEN
4507             message_string = 'No pavement type defined. Please define pavement type to enable' // &
4508                              ' deposition calculation'
4509             CALL message( 'chem_depo', 'PA0745', 1, 2, 0, 6, 0 )
4510          ELSEIF ( luu_palm == ind_lup_asph_conc )  THEN
4511             luu_dep = 9
4512          ELSEIF ( luu_palm == ind_lup_asph )  THEN
4513             luu_dep = 9
4514          ELSEIF ( luu_palm ==  ind_lup_conc )  THEN
4515             luu_dep = 9
4516          ELSEIF ( luu_palm ==  ind_lup_sett )  THEN
4517             luu_dep = 9
4518          ELSEIF ( luu_palm == ind_lup_pav_stones )  THEN
4519             luu_dep = 9
4520          ELSEIF ( luu_palm == ind_lup_cobblest )  THEN
4521             luu_dep = 9
4522          ELSEIF ( luu_palm == ind_lup_metal )  THEN
4523             luu_dep = 9
4524          ELSEIF ( luu_palm == ind_lup_wood )  THEN
4525             luu_dep = 9
4526          ELSEIF ( luu_palm == ind_lup_gravel )  THEN
4527             luu_dep = 9
4528          ELSEIF ( luu_palm == ind_lup_f_gravel )  THEN
4529             luu_dep = 9
4530          ELSEIF ( luu_palm == ind_lup_pebblest )  THEN
4531             luu_dep = 9
4532          ELSEIF ( luu_palm == ind_lup_woodchips )  THEN
4533             luu_dep = 9
4534          ELSEIF ( luu_palm == ind_lup_tartan )  THEN
4535             luu_dep = 9
4536          ELSEIF ( luu_palm == ind_lup_art_turf )  THEN
4537             luu_dep = 9
4538          ELSEIF ( luu_palm == ind_lup_clay )  THEN
4539             luu_dep = 9
4540          ENDIF
4541       ENDIF
4542
4543       IF ( surf_usm_h(0)%frac(m,ind_wat_win) > 0 )  THEN
4544!
4545!--       For windows in USM assume metal as this is as close as we get, assumed LU class desert as
4546!--       also assumed for pavements in LSM.
4547          lud_palm = ind_lup_metal
4548          IF ( lud_palm == ind_lup_user )  THEN
4549             message_string = 'No usm-pavement type defined. Please define pavement type to ' //   &
4550                              'enable deposition calculation'
4551             CALL message( 'chem_depo', 'PA0743', 1, 2, 0, 6, 0 )
4552          ELSEIF ( lud_palm == ind_lup_asph_conc )  THEN
4553             lud_dep = 9
4554          ELSEIF ( lud_palm == ind_lup_asph )  THEN
4555             lud_dep = 9
4556          ELSEIF ( lud_palm ==  ind_lup_conc )  THEN
4557             lud_dep = 9
4558          ELSEIF ( lud_palm ==  ind_lup_sett )  THEN
4559             lud_dep = 9
4560          ELSEIF ( lud_palm == ind_lup_pav_stones )  THEN
4561             lud_dep = 9
4562          ELSEIF ( lud_palm == ind_lup_cobblest )  THEN
4563             lud_dep = 9
4564          ELSEIF ( lud_palm == ind_lup_metal )  THEN
4565             lud_dep = 9
4566          ELSEIF ( lud_palm == ind_lup_wood )  THEN
4567             lud_dep = 9
4568          ELSEIF ( lud_palm == ind_lup_gravel )  THEN
4569             lud_dep = 9
4570          ELSEIF ( lud_palm == ind_lup_f_gravel )  THEN
4571             lud_dep = 9
4572          ELSEIF ( lud_palm == ind_lup_pebblest )  THEN
4573             lud_dep = 9
4574          ELSEIF ( lud_palm == ind_lup_woodchips )  THEN
4575             lud_dep = 9
4576          ELSEIF ( lud_palm == ind_lup_tartan )  THEN
4577             lud_dep = 9
4578          ELSEIF ( lud_palm == ind_lup_art_turf )  THEN
4579             lud_dep = 9
4580          ELSEIF ( lud_palm == ind_lup_clay )  THEN
4581             lud_dep = 9
4582          ENDIF
4583       ENDIF
4584!
4585!--    @TODO: Activate these lines as soon as new ebsolver branch is merged:
4586!--    Set wetness indicator to dry or wet for usm vegetation or pavement
4587       !IF ( surf_usm_h(0)%c_liq(m) > 0 )  THEN
4588       !   nwet = 1
4589       !ELSE
4590       nwet = 0
4591       !ENDIF
4592!
4593!--    Compute length of time step
4594       IF ( call_chem_at_all_substeps )  THEN
4595          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
4596       ELSE
4597          dt_chem = dt_3d
4598       ENDIF
4599
4600       dh = dzw(k)
4601       inv_dh = 1.0_wp / dh
4602       dt_dh = dt_chem / dh
4603!
4604!--    Concentration at i,j,k
4605