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

Last change on this file since 4899 was 4899, checked in by raasch, 11 months ago

small adjustments regarding r4897, bugfix for r4898 (missing preprocessor directive added)

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