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

Last change on this file since 4272 was 4272, checked in by schwenkel, 20 months ago

further modularization of boundary conditions: moving boundary conditions to their modules

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