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

Last change on this file since 4881 was 4881, checked in by forkel, 3 years ago

removed unused parameters and write statements in chemistry

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