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

Last change on this file since 4577 was 4577, checked in by raasch, 7 months ago

further re-formatting to follow the PALM coding standard

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