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

Last change on this file since 3887 was 3887, checked in by schwenkel, 2 years ago

bugfix for chemistry_model_mod via introducing module_interface_exchange_horiz

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