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

Last change on this file since 4511 was 4511, checked in by raasch, 14 months ago

chemistry decycling replaced by explicit setting of lateral boundary conditions

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