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

Last change on this file since 4544 was 4544, checked in by raasch, 12 months ago

conc_av changed from pointer to allocatable array, array spec_conc_av removed

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