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

Last change on this file since 4029 was 4029, checked in by raasch, 2 years ago

bugfix: decycling of chemistry species after nesting data transfer, exchange of ghost points and boundary conditions separated for chemical species and SALSA module, nest_chemistry option removed, netcdf variable NF90_NOFILL is used as argument instead of 1 in calls to NF90_DEF_VAR_FILL

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