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

Last change on this file since 4887 was 4887, checked in by banzhafs, 3 years ago

Unnecessary comments removed from chemistry_model_mod and chem_emissions_mod

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