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

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

further re-numbering of message IDs

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