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

Last change on this file since 4581 was 4581, checked in by suehring, 11 months ago

mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill

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