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

Last change on this file since 4304 was 4304, checked in by banzhafs, 17 months ago

Precision clean-up in drydepo_aero_zhang_vd subroutine in chemistry_model_mod

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