source: palm/trunk/SOURCE/chemistry_model_mod.f90

Last change on this file was 4901, checked in by banzhafs, 4 months ago

removed unnecessary informative message and unused code from chemistry_model_mod

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