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

Last change on this file since 4441 was 4441, checked in by suehring, 19 months ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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