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

Last change on this file since 4403 was 4403, checked in by banzhafs, 20 months ago

chemistry model: implemented on-demand emission read mode for LOD 2

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